text:
using SimJulia

function and_callback(ev::AbstractEvent)
  println("Both events are triggered: $(value(ev))")
end

function or_callback(ev::AbstractEvent)
  println("One of both events is triggered: $(value(ev))")
end

function or_callback_succeed(ev::AbstractEvent, ev2::AbstractEvent)
  println("One of both events is triggered: $(value(ev))")
  succeed(ev2)
end

function or_callback_fail(ev::AbstractEvent, ev2::AbstractEvent)
  println("One of both events is triggered: $(value(ev))")
  fail(ev2, TestException())
end

sim = Simulation()
ev1 = timeout(sim, 1)
ev2 = Event(sim)
@callback and_callback(ev1 & ev2)
@callback or_callback_succeed(ev1 | ev2, ev2)
run(sim)

sim = Simulation()
ev1 = timeout(sim, 1, value=TestException())
ev2 = Event(sim)
@callback or_callback_fail(ev1 | ev2, ev2)
@callback and_callback(ev1 & ev2)
run(sim)

sim = Simulation()
ev1 = timeout(sim)
ev2 = timeout(sim, value=TestException())
@callback or_callback(ev1 | ev2)
run(sim)

sim = Simulation()
ev1 = timeout(sim)
ev2 = timeout(sim)
@callback or_callback(ev1 | ev2)
run(sim)

sim = Simulation()
ev1 = timeout(sim, 1)
ev2 = Event(sim)
ev3 = timeout(sim, 2)
@callback and_callback(AllOf(ev1, ev2, ev3))
@callback or_callback_succeed(AnyOf(ev1, ev2, ev3), ev2)
run(sim)

sim = Simulation()
ev1 = timeout(sim, 1, value=TestException())
ev2 = Event(sim)
ev3 = Event(sim)
@callback or_callback_fail(AnyOf(ev1, ev2, ev3), ev2)
@callback and_callback(AllOf(ev1, ev2, ev3))
run(sim)

sim = Simulation()
ev1 = timeout(sim)
ev2 = timeout(sim, value=TestException())
ev3 = timeout(sim)
@callback or_callback(AnyOf(ev1, ev2, ev3))
run(sim)

sim = Simulation()
ev1 = timeout(sim)
ev2 = timeout(sim)
ev3 = timeout(sim)
@callback or_callback(AnyOf(ev1, ev2, ev3))
run(sim)


text:
function ungetindex!(dx::AbstractArray, x::AbstractArray, ds, i...)
    dx[i...] = ds
    return dx
end

function ungetindex!(dx::AbstractArray, x::AbstractArray, ds, i::AbstractArray{Int})
    dx[i] = ds
    return dx
end

function ungetindex(x::AbstractArray, ds, i...)
    dx = zero(x)
    return ungetindex!(dx, x, ds, i...)
end

function ungetindex(x::Tuple, dy, I...)
    dx = map(1:length(x)) do i
        i in I ? dy : zero(x[i])
    end
    return dx
end

function sum_grad(x::AbstractArray, ds)
    dx = similar(x)
    dx .= ds
    return dx
end

# function mean_grad(x::AbstractArray, ds)
#     dx = similar(x)
#     dx .= ds ./ length(x)
#     return dx
# end

function ∇mean(x::AbstractArray, dy, dims=1:ndims(x))
    dx = similar(x)
    dx .= dy ./ prod(size(x, d) for d in dims)
    return dx
end

function sum_dropdims(x::AbstractArray, dims)
    return dropdims(sum(x; dims=dims); dims=dims)
end

# unbroadcast from Flux
# in in-place version we can consider sum!(similar(x), ds),
# but we need to carefully measure performance in each case

# reshape Δ to be consistent with x
trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

function unbroadcast(x::AbstractArray, Δ)
    if size(x) == size(Δ)
        return Δ
    elseif length(x) == length(Δ)
        return trim(x, Δ)
    else
        sum_dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))
        return trim(x, sum(Δ, dims=sum_dims))
    end
end

unbroadcast(x::Number, Δ) = sum(Δ)

function unbroadcast_prod_x(x::AbstractArray, y::AbstractArray, Δ)
    if size(x) == size(Δ)
        return Δ .* y
    elseif length(x) == length(Δ)
        return trim(x, Δ .* y)
    else
        sum_dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))
        return trim(x, sum(Δ.* y, dims=sum_dims))
    end
end
unbroadcast_prod_y(x::AbstractArray, y::AbstractArray, Δ) = unbroadcast_prod_x(y, x, Δ)

device_like(example, a) = (device = guess_device([example]); device(a))
unbroadcast_prod_x(x::Number, y::AbstractArray, Δ) = unbroadcast_prod_x(device_like(y, [x]), y, Δ)[1]
unbroadcast_prod_x(x::AbstractArray, y::Number, Δ) = unbroadcast_prod_x(x, device_like(x, [y]), Δ)
unbroadcast_prod_y(x::AbstractArray, y::Number, Δ) = unbroadcast_prod_y(x, device_like(x, [y]), Δ)[1]
unbroadcast_prod_y(x::Number, y::AbstractArray, Δ) = unbroadcast_prod_y(device_like(y, [x]), y, Δ)

untranspose_vec(ds::Transpose{T, <:AbstractVector{T}}) where T = transpose(ds)
untranspose_vec(ds::Adjoint{T, <:AbstractVector{T}}) where T = adjoint(ds)
untranspose_vec(ds::AbstractMatrix) = dropdims(transpose(ds); dims=2)

function unvcat(dy::AbstractArray, n::Int, arrs::AbstractArray...)
    a = arrs[n]
    from = n == 1 ? 1 : sum(size(arr, 1) for arr in arrs[1:n-1]) + 1
    to = from + size(a, 1) - 1
    return dy[from:to, [(:) for i=1:length(size(dy)) - 1]...]
end

function unhcat(dy::AbstractArray, n::Int, arrs::AbstractArray...)
    a = arrs[n]
    from = n == 1 ? 1 : sum(size(arr, 2) for arr in arrs[1:n-1]) + 1
    to = from + size(a, 2) - 1
    return dy[:, from:to, [(:) for i=1:length(size(dy)) - 2]...]
end

function uncat(dy::AbstractArray, n::Int, arrs::AbstractArray...; dims)
    @assert(dims isa Integer, "Can only undo cat() over a single dimension, " *
            "but dimensions $dims were provided")
    dim = dims
    a = arrs[n]
    from = n == 1 ? 1 : sum(size(arr, dim) for arr in arrs[1:n-1]) + 1
    to = from + size(a, dim) - 1
    return dy[[(:) for i=1:dim - 1]..., from:to, [(:) for i=1:length(size(dy)) - dim]...]
end

namedtuple(names, values) = NamedTuple{names}(values)

function rev_perm(perm::NTuple{N, Int}) where N
    rperm = Vector{Int}(undef, length(perm))
    for (i, j) in enumerate(perm)
        rperm[j] = i
    end
    return tuple(rperm...)
end

function ∇permutedims(dy, perm)
    rperm = rev_perm(perm)
    return permutedims(dy, rperm)
end

ext: jl


text:
using Blpapi
using CSV

# Verbindungsaufbau
con = createConn()

t = floor((now() - Dates.Month(0)), Dates.Day(1))
startDate = floor((t - Dates.Day(1)), Dates.Minute(30))

    # Instrumen deklaration
inst = instrument("RXU8 Comdty", "TRADE", startDate, t)

    #request data
BarFrame(con, inst, 30)

z = bundDaten[[:Dates, :Volume]]

sort(z, cols = [:Volume])[(end-19):end,:]

i
TS(con, inst)

bdh(connection.session, "IBM US Equity", "TRADE", startDate, t)

bdh(connection.session, ["IBM US Equity"], ["PX_LAST"], "20180101", "20180331")
bdh(connection.session, [inst.ticker], "PX_Last", "20180101", "20180331")

barsize = 30

# Pulling all from BBG
Response = bar(connection.session, inst.ticker, inst.fields, string(inst.startDate), string(inst.endDate); interval=barsize)

#sorting into DataFrame

df = DataFrame()
df[:Dates] = Response.datetimes
df[:Open] = zeros(length(Response.datetimes))
df[:High] = zeros(length(Response.datetimes))
df[:Low] = zeros(length(Response.datetimes))
df[:Close] = zeros(length(Response.datetimes))
df[:Volume] = zeros(length(Response.datetimes))

for i in 1:length(Response.datetimes)
df[i,:Open] = Response[Response.datetimes[i]].openVar
df[i,:High] = Response[Response.datetimes[i]].highVar
df[i,:Low] = Response[Response.datetimes[i]].lowVar
df[i,:Close] = Response[Response.datetimes[i]].closeVar
df[i,:Volume] = Response[Response.datetimes[i]].volumeVar
end

using TimeZones

df

d = df[:,1]

zd = ZonedDateTime.(d, tz"UTC")

finalDT = astimezone.(zd, tz"Europe/Zurich")

df[:,1] = finalDT

df

Date(finalDT[1], dateformat"y-m-d")

con = createConn()
"PX_High", "PX_Low", "PX_Last"
# For historical end of day, mind the date format
inst = instrument("RXU8 Comdty", ["PX_Open", "PX_High", "PX_Low", "PX_Last"], "20180501", "20180625")

#function TS(con = connection, inst = ibm)

#endData = Array()

for i in (1:length(inst.fields))

    response = bdh(con.session, [inst.ticker], [inst.fields[i]], inst.startDate, inst.endDate; periodicitySelection = "DAILY")
    data = response[inst.ticker, inst.fields[i]]

    dformat = Dates.DateFormat("y-m-d")
    dates = similar(data[:,1], Date)
    prices = similar(data[:,2], Float64)

    for j in (1:length(data[:,1]))

        dates[j] = Date(data[j,1], dformat)
        prices[j] = data[j,2]

    end

    k = sortperm(dates)
    dates = dates[k]
    prices = prices[k]

    #if i > 1
    #    data2 = similar(data)

    data = TimeArray(dates, prices, [inst.fields[i]])

    # if i == 1
    #     dataFinal = similar(1:length(data), Float64)
    #     dataFinal = data
    # end
    #
    # if i > 1
    #     dataFinal = merge(dataFinal, data)
    # end
    #

end

connection = createConn()
bund = instrument("RXU8 Comdty", ["PX_Open", "PX_High", "PX_Low", "PX_Last"], "20180601", "20180627")

function TS(con = connection, inst = ibm, field = "PX_Last")

  response = bdh(con.session, [inst.ticker], [field], inst.startDate, inst.endDate; periodicitySelection = "DAILY")
  data = response[inst.ticker, field]

  dformat = Dates.DateFormat("y-m-d")
  dates = similar(data[:,1], Date)
  prices = similar(data[:,2], Float64)

  for i in (1:length(data[:,1]))
    dates[i] = Date(data[i,1], dformat)
    prices[i] = data[i,2]
  end

  k = sortperm(dates)
  dates = dates[k]
  prices = prices[k]

  data = TimeArray(dates, prices, [field])

end

TS(connection, bund, bund.fields[1])

data = TS(connection, bund)

colnames(data)

Dates.firstdayofweek(Date(now()-Dates.Day(5))) + Dates.Day(1)
Dates.firstdayofweek(Date(now()-Dates.Day(5)))

data[Dates.firstdayofweek(Date(now()-Dates.Day(5)))]
data[Dates.firstdayofweek(Date(now()-Dates.Day(5))) + Dates.Day(1)]

data.timestamp

#end
for

    #refresh end time to now

    #call BarFrame

    #write result to csv./ or call this function from excel

function barUpdate(contract, con = con)

    t = now()
    startDate = t - Dates.Month(2)

    inst = instrument(contract.ticker, contract.fieldnames, startDate, t)

    BarFrame(con, inst, 30)

end

quantile(barUpdate()[(end-300):end,:Volume],0.95)
log(bund2[:High]./bund2[:Low]) ./ log(bund2[:Close])

(bund2[end,:High] - bund2[end,:Low]) / maximum(bund2[:High] - bund2[:Low])

for i in keys(instContracts)
       contract = instContracts[i]
       barUpdate(contract)
end

function volIndicatorXL(commodity, days_back)

    inst = bbgTS.contracts[commodity]
    inst.startDate = now() - Dates.Month(2)
    inst.endDate = now()
    df = BarFrame(bbgTS.conInternal, inst, bbgTS.n)
    df = Array{Any}(df[Date(df[:,:Dates]) .== Date(now())-Dates.Day(days_back), [:Dates, :Volume, :RotFact]])
    df[:,1] = toXLDate.(df[:,1])
    return df

end


text:
"""
    RPModel{S<:AbstractString, T<:AbstractFloat, A<:AbstractMatrix{T}, H<:Integer}

Random projection model. It constructs from a document term matrix (DTM)
a model that can be used to embed documents in a random sub-space. The model requires
that the document term matrix be a `DocumentTermMatrix{T<:AbstractFloat}` because
the elements of the matrices resulted projection operation are floating point
numbers and these have to match or be convertible to type `T`. The approach is
based on the effects of the
[Johnson-Lindenstrauss lemma](https://en.wikipedia.org/wiki/Johnson%E2%80%93Lindenstrauss_lemma).

# Fields
  * `vocab::Vector{S}` a vector with all the words in the corpus
  * `vocab_hash::OrderedDict{S,H}` a word to index in the random projection maatrix mapping
  * `R::A` the random projection matrix
  * `stats::Symbol` the statistical measure to use for word importances in documents. Available values are: `:count` (term count), `:tf` (term frequency), `:tfidf` (default, term frequency-inverse document frequency) and `:bm25` (Okapi BM25)
  * `idf::Vector{T}` inverse document frequencies for the words in the vocabulary
  * `nwords::T` averge number of words in a document
  * `ngram_complexity::Int` ngram complexity
  * `κ::Int` the `κ` parameter of the BM25 statistic
  * `β::Float64` the `β` parameter of the BM25 statistic
  * `project::Bool` specifies whether the model actually performs the projection or not; it is false if the number of dimensions provided is zero or negative

# References:
  * [Kaski 1998](http://www.academia.edu/371863/Dimensionality_Reduction_by_Random_Mapping_Fast_Similarity_Computation_for_Clustering)
  * [Achlioptas 2001](https://users.soe.ucsc.edu/~optas/papers/jl.pdf)
  * [Li et al. 2006](http://web.stanford.edu/~hastie/Papers/Ping/KDD06_rp.pdf)
"""
struct RPModel{S<:AbstractString, T<:AbstractFloat, A<:AbstractMatrix{T}, H<:Integer}
    vocab::Vector{S}        # vocabulary
    vocab_hash::OrderedDict{S,H}   # term to column index in V
    R::A                    # projection matrix
    stats::Symbol           # term/document importance
    idf::Vector{T}          # inverse document frequencies
    nwords::T               # average words/document in corpus
    ngram_complexity::Int   # ngram complexity
    κ::Int                  # κ parameter for Okapi BM25 (used if stats==:bm25)
    β::Float64              # β parameter for Okapi BM25 (used if stats==:bm25)
    project::Bool
end

function RPModel(dtm::DocumentTermMatrix{T};
                 k::Int=size(dtm.dtm, 1),
				 density::Float64=1/sqrt(k),
                 stats::Symbol=:tfidf,
                 ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                 κ::Int=BM25_KAPPA,
                 β::Float64=BM25_BETA
                ) where T<:AbstractFloat
    m, n = size(dtm.dtm)
    zeroval = zero(T)
    # Checks
    length(dtm.terms) != m &&
        throw(DimensionMismatch("Dimensions inside dtm are inconsistent."))
    if !(stats in [:count, :tf, :tfidf, :bm25])
        @warn "stats has to be either :tf, :tfidf or :bm25; defaulting to :tfidf"
        stats = :tfidf
    end
    # Calculate inverse document frequency, mean document size
    documents_containing_term = vec(sum(dtm.dtm .> 0, dims=2)) .+ one(T)
    idf = log.(n ./ documents_containing_term) .+ one(T)
    nwords = mean(sum(dtm.dtm, dims=1))
    R = random_projection_matrix(k, m, T, density)
    project = ifelse(k > 0, true, false)
    # Return the model
    return RPModel(dtm.terms, dtm.row_indices, R, stats, idf,
                   nwords, ngram_complexity, κ, β, project)
end

function RPModel(dtm::DocumentTermMatrix{T}; kwargs...) where T<:Integer
    throw(ErrorException(
        """A random projection model requires a that the document term matrix
        be a DocumentTermMatrix{<:AbstractFloat}!"""))
end

"""
    random_projection_matrix(k::Int, m::Int, eltype::Type{T<:AbstractFloat}, density::Float64)

Builds a `k`×`m` sparse random projection matrix with elements of type `T` and
a non-zero element frequency of `density`. `k` and `m` are the output and input
dimensionalities.

# Matrix Probabilities
If we note `s = 1 / density`, the components of the random matrix are drawn from:
- `-sqrt(s) / sqrt(k)` with probability `1/2s`
- `0` with probability `1 - 1/s`
- `+sqrt(s) / sqrt(k)`   with probability `1/2s`

# No projection hack
If `k<=0` no projection is performed and the function returns an identity matrix
sized `m`×`m` with elements of type `T`. This is useful if one does not want to
embed documents but rather calculate term frequencies, BM25 and other statistical
indicators (similar to `dtv`).
"""
function random_projection_matrix(k::Int, m::Int, eltype::Type{T}, density::Float64
                                 ) where T<:AbstractFloat
    if k <= 0
        # No projection, return the identity matrix
        R = spdiagm(0 => ones(T, m))
        return R
    else
        R = zeros(T, k, m)
        s = 1/density
        is_pos = 0.0
        is_neg = 0.0
        v = sqrt(s/k)
        pmin = 1/(2*s)
        @inbounds for j in 1:m
            for i in 1:k
                p = rand()
                if p < pmin
                    R[i,j] = v
                elseif p > 1-pmin
                    R[i,j] = -v
                end
            end
        end
        return sparse(R)
    end
end

function Base.show(io::IO, rpm::RPModel{S,T,A,H}) where {S,T,A,H}
    len_vecs, num_terms = size(rpm.R)
    str_proj = ifelse(rpm.project, "Random Projection model",
                      "Identity Projection")
    print(io, "$str_proj ($(rpm.stats)), " *
          "$(num_terms) terms, dimensionality $(len_vecs), $(T) vectors")
end

"""
    rp(X [;k=m, density=1/sqrt(k), stats=:tfidf, ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, κ=2, β=0.75])

Constructs a random projection model. The input `X` can be a `Corpus` or a `DocumentTermMatrix`
with `m` words in the lexicon. The model does not store the corpus or DTM document embeddings,
just the projection matrix. Use `?RPModel` for more details.
"""
function rp(dtm::DocumentTermMatrix{T};
            k::Int=size(dtm.dtm, 1),
            density::Float64=1/sqrt(k),
            stats::Symbol=:tfidf,
            ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
            κ::Int=BM25_KAPPA,
            β::Float64=BM25_BETA
           ) where T<:AbstractFloat
    RPModel(dtm, k=k, density=density, stats=stats,
            ngram_complexity=ngram_complexity, κ=κ, β=β)
end

function rp(crps::Corpus,
            ::Type{T} = DEFAULT_FLOAT_TYPE;
            k::Int=length(crps.lexicon),
            density::Float64=1/sqrt(k),
            stats::Symbol=:tfidf,
            ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
            κ::Int=BM25_KAPPA,
            β::Float64=BM25_BETA
           ) where T<:AbstractFloat
    lex = create_lexicon(crps, ngram_complexity)
    rp(DocumentTermMatrix{T}(crps, lex, ngram_complexity=ngram_complexity),
       k=k, density=density, stats=stats, ngram_complexity=ngram_complexity, κ=κ, β=β)
end

"""
    vocabulary(rpm)

Return the vocabulary as a vector of words of the random projection model `rpm`.
"""
vocabulary(rpm::RPModel) = rpm.vocab

"""
    in_vocabulary(rpm, word)

Return `true` if `word` is part of the vocabulary of the random projection
model `rpm` and `false` otherwise.
"""
in_vocabulary(rpm::RPModel, word::AbstractString) = word in rpm.vocab

"""
    size(rpm)

Return a tuple containing the input data and projection sub-space
dimensionalities of the random projection model `rpm`.
"""
size(rpm::RPModel) = size(rpm.R, 2), size(rpm.R, 1)

"""
    index(rpm, word)

Return the index of `word` from the random projection model `rpm`.
"""
index(rpm::RPModel, word) = rpm.vocab_hash[word]

"""
    get_vector(rpm, word)

Returns the random projection vector corresponding to `word` in the
random projection model `rpm`.
"""
function get_vector(rpm::RPModel{S,T,A,H}, word) where {S,T,A,H}
    default = zeros(T, size(rpm.R, 1))
    idx = get(rpm.vocab_hash, word, 0)
    if idx == 0
        return default
    else
        return rpm.R[:, idx]
    end
end

"""
    embed_document(rpm, doc)

Return the vector representation of `doc`, obtained using the
random projection model `rpm`. `doc` can be an `AbstractDocument`,
`Corpus` or DTV or DTM.
"""
embed_document(rpm::RPModel{S,T,A,H}, doc::AbstractDocument) where {S,T,A,H} =
    # Hijack vocabulary hash to use as lexicon (only the keys needed)
    embed_document(rpm, dtv(doc, rpm.vocab_hash, T,
                            ngram_complexity=rpm.ngram_complexity,
                            lex_is_row_indices=true))

embed_document(rpm::RPModel{S,T,A,H}, doc::AbstractString) where {S,T,A,H} =
    embed_document(rpm, NGramDocument{S}(doc, DocumentMetadata(), rpm.ngram_complexity))

embed_document(rpm::RPModel{S,T,A,H}, doc::AbstractVector{S2}) where {S,T,A,H,S2<:AbstractString} =
    embed_document(rpm, TokenDocument{S}(doc))

# Actual embedding function: takes as input the random projection model `rpm` and a document
# term vector `dtv`. Returns the representation of `dtv` in the embedding space.
function embed_document(rpm::RPModel{S,T,A,H}, dtv::AbstractVector{T}) where {S,T,A,H}
    words_in_document = sum(dtv)
    # Calculate document vector
    if rpm.stats == :count
        v = dtv
    elseif rpm.stats == :tf
        v = sqrt.(dtv ./ max(words_in_document, one(T)))
    elseif rpm.stats == :tfidf
        v = sqrt.(dtv ./ max(words_in_document, one(T))) .* rpm.idf
    elseif rpm.stats == :bm25
        k = T(rpm.κ)
        b = T(rpm.β)
        tf = sqrt.(dtv ./ max(words_in_document, one(T)))
        v = rpm.idf .* (k + 1) .* tf ./
                       (k * (one(T) - b + b * words_in_document/rpm.nwords) .+ tf)
    end
    # Embed
    local d̂
    if rpm.project
        d̂ = rpm.R * v  # embed
    else
        d̂ = v
    end
    d̂ = d̂ ./ (norm(d̂,2) .+ eps(T))  # normalize
    return d̂
end

function embed_document(rpm::RPModel{S,T,A,H}, dtm::DocumentTermMatrix{T}) where {S,T,A,H}
    n = size(dtm.dtm,1)
    k = size(rpm.R, 1)
    if rpm.stats == :count
        X = dtm.dtm
    elseif rpm.stats == :tf
        X = tf(dtm)
    elseif rpm.stats == :tfidf
        X = tf_idf(dtm)
    elseif rpm.stats == :bm25
        X = bm_25(dtm, κ=rpm.κ, β=rpm.β)
    end
    local U
    if rpm.project
        U = rpm.R * X
    else
        U = X
    end
    U ./= (sqrt.(sum(U.^2, dims=1)) .+ eps(T))
    return U
end

function embed_document(rpm::RPModel{S,T,A,H}, crps::Corpus) where {S,T,A,H}
    lex = create_lexicon(crps, rpm.ngram_complexity)
    embed_document(rpm, DocumentTermMatrix{T}(crps, lex, ngram_complexity=rpm.ngram_complexity))
end

"""
    save_rp_model(rpm, filename)

Saves an random projection model `rpm` to disc in file `filename`.
"""
function save_rp_model(rpm::RPModel{S,T,A,H}, filename::AbstractString) where {S,T,A,H}
    k, nwords = size(rpm.R)
    open(filename, "w") do fid
        println(fid, "Random Projection Model saved at $(Dates.now())")
        println(fid, "$nwords $k")  # number of words, k
        println(fid, rpm.project)
        println(fid, rpm.stats)
        writedlm(fid, rpm.idf', " ")
        println(fid, rpm.nwords)
        println(fid, rpm.ngram_complexity)
        println(fid, rpm.κ)
        println(fid, rpm.β)
        # Vocabulary
        mv = Matrix{String}(undef, 1, nwords)
        mv[1,:] = rpm.vocab
        writedlm(fid, mv, " ")
        # Random projection matrix
        writedlm(fid, rpm.R, " ")
    end
end

"""
    load_rp_model(filename, eltype; [sparse=true])

Loads an random projection model from `filename` into an random projection model object.
The projection matrix element type is specified by `eltype` (default `DEFAULT_FLOAT_TYPE`)
while the keyword argument `sparse` specifies whether the matrix should be sparse or not.
"""
function load_rp_model(filename::AbstractString, ::Type{T}=DEFAULT_FLOAT_TYPE;
                       sparse::Bool=true) where T<: AbstractFloat
    # Matrix type for random projection model
    A = ifelse(sparse, SparseMatrixCSC{T, Int}, Matrix{T})
    # Define parsed variables local to outer scope of do statement
    local vocab, vocab_hash, R, stats, idf, nwords, ngram_complexity, κ, β, project
    open(filename, "r") do fid
        readline(fid)  # first line, header
        line = readline(fid)
        vocab_size, k = map(x -> parse(Int, x), split(line, ' '))
        # Preallocate
        if sparse
            R = SparseMatrixCSC{T, Int}(UniformScaling(0), k, vocab_size)
        else
            R = zeros(T, k, vocab_size)
        end
        # Start parsing the rest of the file
        project = parse(Bool, strip(readline(fid)))
        stats = Symbol(strip(readline(fid)))
        idf = map(x->parse(T, x), split(readline(fid), ' '))
        nwords = parse(T, readline(fid))
        ngram_complexity = parse(Int, readline(fid))
        κ = parse(Int, readline(fid))
        β = parse(Float64, readline(fid))
        vocab = Vector{String}(split(readline(fid),' '))
        vocab_hash = OrderedDict((v,i) for (i,v) in enumerate(vocab))
        for i in 1:k
            R[i,:] = map(x->parse(T,x), split(readline(fid), ' '))
        end
    end
    return RPModel{String,T,A,Int}(vocab, vocab_hash, R, stats, idf,
                                      nwords, ngram_complexity, κ, β, project)
end


text:
export Actor
using Flux

struct Actor
   state_size::Int8
   action_size::Int8
   layer_size = Int8[5,6,4]
   learning_rate::Float64=0.001
   batch_size::Int8=6
   tau::Float64= 0.01
   """ Connstructor for the Actor Network
     state_size:An integer denoting the dimensionality of the states
            in the current problem
     action_size: An integer denoting the dimensionality of the
            actions in the current problem
      hidden_units: An iterable defining the number of hidden units in
            each layer
      learning_rate: A fload denoting the speed at which the network
           will learn. default
      batch_size: An integer denoting the batch size
      tau: A flot denoting the rate at which the target model will
            track the main model. Formally, the tracking function is defined as:
              target_weights = tau * main_weights + (1 - tau) * target_weights """

   function generate_model(s::Int8=state_size,I::Array{Float64})
       """Generates the model based on the hyperparameters defined in the constructor"""
       input_layer = I
       weights1 = rand(5,6)
       weights2 = rand(6,4)
       model = Flux.chain(Dense(param(weights1),relu),Dense(param(weights2),sigmoid));
       return model(I) ,weights1, I
    end
    function train(states::Array{Float64},action_gradients::Array{Float64})
         """Updates the weights of the main network
            states: The states of the input to the network
            action_gradients: The gradients of the actions to update the
            network"""
         model1,model1_weight,model1_input = generate_model()
         for i = 0:100
            model1_weight .-= learning_rate .*action_gradients
            function train_target_model(m::Array({Float64},5,6))
            """ Updates the weights of the target network to slowly track the main
        network.
        The speed at which the target network tracks the main network is
        defined by tau, given in the constructor to this class. Formally,
        the tracking function is defined as:
            target_weights = tau * main_weights + (1 - tau) * target_weights"""
                main_weights = m
                target_model,target_model_weight, target_state = generate_model()
                target_weights = target_model_weight
                target_weights = tau.*main_weight .+ (1 - tau).*target_weight
                target_model_weight = target_weights
             end
         end
    end


text:
using TemplateMatching
using Documenter

DocMeta.setdocmeta!(TemplateMatching, :DocTestSetup, :(using TemplateMatching); recursive=true)

makedocs(;
    modules=[TemplateMatching],
    authors="Stefano Campanella <15182642+stefanocampanella@users.noreply.github.com> and contributors",
    repo="https://github.com/stefanocampanella/TemplateMatching.jl/blob/{commit}{path}#{line}",
    sitename="TemplateMatching.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://stefanocampanella.github.io/TemplateMatching.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/stefanocampanella/TemplateMatching.jl",
    devbranch="master",
)


text:
using Test, Limbdark, LinearAlgebra, Statistics, ForwardDiff, DiffResults, Random

include("loglinspace.jl")

include("test_cel_vec.jl")
include("test_s2_gradient.jl")
include("test_transit_poly.jl")
include("test_Mn.jl")

# Skip plots for CI tests
skip_plots = true
include("test_transit_poly_gradient2.jl")
include("test_integrate_transit_simpson_vec.jl")


text:
#------------------------------------------------------------------------------
"""
    printmat([fh::IO],x,width=10,prec=3,NoPrinting=false,htmlQ=false)

Print all elements of matrix with predefined formatting.

# Input
- `fh::IO`:           (optional) file handle. If not supplied, prints to screen
- `x::Array`:         string, date or array to print
- `width::Int`:       (optional) scalar, width of printed cells. [10]
- `prec::Int`:        (optional) scalar, precision of printed cells. []
- `NoPrinting::Bool`: (optional) bool, true: no printing, just return formatted string
- `hmtlQ::Bool`:      (optional) bool, true: format as htmlQ cells

# Output
- str         (if NoPrinting) string, (otherwise nothing)

# Examples. Try the printing the following arrays:
- x = [11 12;21 22]
- x = Any[1 "ab"; Date(2018,10,7) 3.14]

# Uses
- fmtNumPs

# To do
- use Dict() for the options, width etc?
- sort out dispatch and keyword arguments

Paul.Soderlind@unisg.ch

"""
function printmat(fh::IO,x,width=10,prec=3,NoPrinting=false,htmlQ=false)

  if isa(x,Union{String,Date,DateTime,Missing})  #these types need special treatment
    str = string(lpad(x,width),"\n")
    if NoPrinting
      return str
    else
      print(fh,str,"\n")
      return nothing
    end
  elseif isa(x,Nothing)
    return nothing
  end

  if ndims(x) > 2
    @warn("more than 2 dimensions")
    return nothing
  end

  (m,n) = (size(x,1),size(x,2))

  iob = IOBuffer()
  for i = 1:m                #loop over lines
    for j = 1:n                #loop over columns
      if isa(x[i,j],AbstractFloat)        #Float
        write(iob,fmtNumPs(x[i,j],width,prec,"right",htmlQ))
      elseif isa(x[i,j],Nothing)           #Nothing
        htmlQ ? write(iob,"",lpad("",width),"	") : write(iob,lpad("",width))
      else                                #other types (Integer,Missing,String,Date,...)
        htmlQ ? write(iob,"",lpad(x[i,j],width),"	") : write(iob,lpad(x[i,j],width))
      end
    end
    write(iob,"\n")            #newline
  end
  str = String(take!(iob))

  if NoPrinting                              #no printing, just return str
    return str
  else                                       #print, return nothing
    print(fh,str,"\n")
    return nothing
  end

end
                  #when fh is not supplied: printing to screen
printmat(x,width=10,prec=3,NoPrinting=false,htmlQ=false) = printmat(stdout::IO,
         x,width,prec,NoPrinting,htmlQ)
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    fmtNumPs(z,width=10,prec=2,justify="right",htmlQ=false)

Formats a scalar and creates a string of it.

The Formatting.jl package provides more elegant solutions:
fmt  = FormatSpec(string(">",width,".",prec,"f"))   #right justified, else "<"
fmt = FormatSpec(string(">",wid,"d"))               #for Int
str  = Formatting.fmt(fmt1,z))

"""
function fmtNumPs(z,width=10,prec=2,justify="right",htmlQ=false)
  if prec > 0                        #if decimal number
    z   = round(z,digits=prec)       #101.23
    str = split(string(z),'.')
    if length(str) > 1
      strR  = string(".",rpad(str[2],prec,"0"))   #.23
      strLR = string(str[1],strR)                 #"101" * ".23"
    else                                          #eg. NaN, missing
      strLR = string(z)
    end
  else
    if typeof(z) <: AbstractFloat                  #Floats
      z = round(Int,z)
    end
    strLR = string(z)
  end
  if justify == "left"
    strLR = rpad(strLR,width)
  else
    strLR = lpad(strLR,width)
  end
  htmlQ && (strLR = string("",strLR,"	"))
  return strLR
end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
global printlnPsOpt = Dict("width"=>10,"prec"=>3)    #options for printlnPs()
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    set_printlnPsOpt(opt)

Set global options (width and prec) to printlnPs.

# Notice
This function is needed if printlnPs is part of a module.
Otherwise we could set printlnPsOpt directly.

"""
function set_printlnPsOpt(opt)
   global printlnPsOpt
   printlnPsOpt = opt
   println("printlnPsOpt is now $printlnPsOpt")
end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    printlnPs([fh::IO],z...)

Subsitute for println, with predefined formatting.

# Input
- `fh::IO`:    (optional) file handle. If not supplied, prints to screen
- `z::String`: string, numbers and arrays to print

# Refers to
- printlnPsOpt, a global Dict("width"=>10,"prec"=>3)

The formatting can be set globally by defining a dictionary in the calling scope.

Paul.Soderlind@unisg.ch

"""
function printlnPs(fh::IO,z...)

  global printlnPsOpt

  width = try                    #getting defaults from dictionary
    get(printlnPsOpt,"width",10)
  catch                          #if Dict() printlnPsOpt is not defined
    10
  end
  prec = try
    get(printlnPsOpt,"prec",3)
  catch
    3
  end

  for x in z                              #loop over inputs in z...
    if isa(x,Union{String,Date,DateTime,Missing})
      print(fh,lpad(x,width))
    elseif isa(x,Nothing)
      print(fh,"")
    else                                         #other types
      iob = IOBuffer()
      for i = 1:length(x)
        if isa(x[i],AbstractFloat)               #Float
          write(iob,fmtNumPs(x[i],width,prec,"right"))
        elseif isa(x[i],Nothing)                 #Nothing
          write(iob,lpad("",width))
        else                                     #Integer, etc
          write(iob,lpad(x[i],width))
        end
      end
      print(fh,String(take!(iob)))
    end
  end

  print(fh,"\n")

end
                      #when fh is not supplied: printing to screen
printlnPs(z...) = printlnPs(stdout::IO,z...)
#------------------------------------------------------------------------------

max_stars_repo_name: AlTaken/EmpiricalFinancePhD

max_issues_repo_name: AlTaken/EmpiricalFinancePhD

max_forks_repo_name: AlTaken/EmpiricalFinancePhD


text:
# Compute statistics on entity data

using JLD, PyPlot, SparseArrays

# Load the data file
file_dir = joinpath(Base.source_dir(),"./Entity.jld")
#E = loadassoc(file_dir)
E = load(file_dir)["E"]

# Calculate number of entities in each category, then count the number of times each entity occurs.
print(sum(logical(col2type(E,"/")),1))
En = sum(logical(E),1)

# Plot the log-log plot of location frequencies. Notice the power-law distribution.
row,entity,count = find(En)
An = Assoc(count,entity,1)
figure()
loglog(sum(getadj(An[:,StartsWith("LOCATION/,")]),dims=2) ,"o")


text:
mutable struct TrackedReal{T<:Real} <: Real
  data::T
  tracker::Tracked{T}
end

TrackedReal(x::Real) = TrackedReal(x, Tracked{typeof(x)}(Call(), zero(x)))

data(x::TrackedReal) = x.data
tracker(x::TrackedReal) = x.tracker

track(f::Call, x::Real) = TrackedReal(x, Tracked{typeof(x)}(f, zero(x)))

function back!(x::TrackedReal; once = true)
    isinf(x) && error("Loss is Inf")
    isnan(x) && error("Loss is NaN")
    return back!(x, 1, once = once)
end

function update!(x::TrackedReal, Δ)
  x.data += data(Δ)
  tracker(x).grad = 0
  return x
end

function Base.show(io::IO, x::TrackedReal)
  T = get(io, :typeinfo, Any)
  show(io, data(x))
  T <: TrackedReal || print(io, " (tracked)")
end

Base.decompose(x::TrackedReal) = Base.decompose(data(x))

Base.copy(x::TrackedReal) = x

Base.convert(::Type{TrackedReal{T}}, x::TrackedReal{T}) where T = x

Base.convert(::Type{TrackedReal{T}}, x::Real) where T = TrackedReal(convert(T, x))

Base.convert(::Type{TrackedReal{T}}, x::TrackedReal{S}) where {T,S} =
  error("Not implemented: convert tracked $S to tracked $T")

for op in [:(==), :≈, :<, :(<=)]
  @eval Base.$op(x::TrackedReal, y::Real) = Base.$op(data(x), y)
  @eval Base.$op(x::Real, y::TrackedReal) = Base.$op(x, data(y))
  @eval Base.$op(x::TrackedReal, y::TrackedReal) = Base.$op(data(x), data(y))
end

Base.eps(x::TrackedReal) = eps(data(x))
Base.eps(::Type{TrackedReal{T}}) where T = eps(T)

for f in :[isinf, isnan, isfinite].args
  @eval Base.$f(x::TrackedReal) = Base.$f(data(x))
end

Base.Printf.fix_dec(x::TrackedReal, n::Int, a...) = Base.Printf.fix_dec(data(x), n, a...)

Base.float(x::TrackedReal) = x

Base.promote_rule(::Type{TrackedReal{S}},::Type{T}) where {S,T} =
  TrackedReal{promote_type(S,T)}

using Random

for f in :[rand, randn, randexp].args
  @eval Random.$f(rng::AbstractRNG,::Type{TrackedReal{T}}) where {T} = param(rand(rng,T))
end

using DiffRules, SpecialFunctions, NaNMath

for (M, f, arity) in DiffRules.diffrules()
  arity == 1 || continue
  @eval begin
    @grad $M.$f(a::Real) =
      $M.$f(data(a)), Δ -> (Δ * $(DiffRules.diffrule(M, f, :a)),)
    $M.$f(a::TrackedReal) = track($M.$f, a)
  end
end

# Work around zero(π) not working, for some reason
_zero(::Irrational) = nothing
_zero(x) = zero(x)

for (M, f, arity) in DiffRules.diffrules()
  arity == 2 || continue
  da, db = DiffRules.diffrule(M, f, :a, :b)
  f = :($M.$f)
  @eval begin
    @grad $f(a::TrackedReal, b::TrackedReal) = $f(data(a), data(b)), Δ -> (Δ * $da, Δ * $db)
    @grad $f(a::TrackedReal, b::Real) = $f(data(a), b), Δ -> (Δ * $da, _zero(b))
    @grad $f(a::Real, b::TrackedReal) = $f(a, data(b)), Δ -> (_zero(a), Δ * $db)
    $f(a::TrackedReal, b::TrackedReal)  = track($f, a, b)
    $f(a::TrackedReal, b::Real) = track($f, a, b)
    $f(a::Real, b::TrackedReal) = track($f, a, b)
  end
end

# Eliminating ambiguity
import Base:^

^(a::TrackedReal, b::Integer) = track(^, a, b)

# Hack for conversions

(T::Type{<:Real})(x::Dual) = Dual(T(x.value), map(T, x.partials.values))

# Tuples

struct TrackedTuple{T<:Tuple}
  data::T
  tracker::Tracked{T}
end

data(xs::TrackedTuple) = xs.data
tracker(xs::TrackedTuple) = xs.tracker

accum!(x::Tuple, Δ::Tuple) = accum!.(x, Δ)
init_grad(x::Tuple) = init_grad.(x)
zero_grad!(x::Tuple) = zero_grad!.(x)

track(f::Call, xs::Tuple) = TrackedTuple(xs, Tracked{typeof(xs)}(f, zero.(xs)))

function Base.show(io::IO, xs::TrackedTuple)
  show(io, data(xs))
  print(io, " (tracked)")
end

Base.length(x::TrackedTuple) = length(data(x))

Base.getindex(xs::TrackedTuple, i::Integer) = track(getindex, xs, i)

@grad function getindex(xs::TrackedTuple, i)
  data(xs)[i], Δ -> (ntuple(j -> i == j ? Δ : 0, length(xs)), nothing)
end

# Array collection

function collect(xs)
  xs = Base.collect(xs)
  track(Call(collect, (tracker.(xs),)), data.(xs))
end

function scan(c::Call{typeof(collect)})
  foreach(scan, c.args[1])
end

function back_(c::Call{typeof(collect)}, Δ, once)
  foreach((x, d) -> back(x, d, once), c.args[1], data(Δ))
end

function back_(g::Grads, c::Call{typeof(collect)}, Δ)
  foreach((x, Δ) -> back(g, x, Δ), c.args[1], Δ)
end


text: Dict{String, Any}("after" => Dict{String, Any}("value" => "3141.5927", "type" => "float"), "exponent" => Dict{String, Any}("value" => "3.0e14", "type" => "float"), "before" => Dict{String, Any}("value" => "3141.5927", "type" => "float"))


text:
export NTget_g, NTget_w, NTget_M
export NesterovToddstep, NTpredcorr, NTgetcentralpathpoint

get_nc(z::PointPrimalDual) = get_nc(z.x)
get_nc(x::PointE) = sum(x.dims) + length(x.vec)

mu(z::PointPrimalDual) = dot(z.x, z.s) / (get_nc(z))

function delta(pb::SDCOContext, z::PointPrimalDual)
    g = NTget_g(pb, z.x, z.s, check=true)
    mu_z = mu(z)

    return 0.5 * norm(hadamard(g, sqrt(mu_z)*inv(z.x) - 1/sqrt(mu_z) * z.s, g, transposefirst=true))
end

function phi_mu(z::PointPrimalDual, mu::Float64)
    for i in eachindex(z.x.mats)
        if (eigmin(z.x.mats[i]) <= 0) || (eigmin(z.s.mats[i]) <= 0)
            return Inf
        end
    end

    for i in eachindex(z.x.vec)
        if (z.x.vec[i] <= 0) || (z.s.vec[i] <= 0)
            return Inf
        end
    end

    return dot(z.x, z.s) + mu * psi(hadamard(z.x, z.s))
end

function psi(x::PointE)
    res = 0

    for mat in x.mats
        res += -logdet(mat)
    end

    for xi in x.vec
        res += -log(xi)
    end
    return res
end

function NTgetcentralpathpoint(pb::SDCOContext, z0::PointPrimalDual{T, Dense{T}}; outlev=1) where T<:Number
    theta = 1/3
    tau = 0.99
    omega = 1e-4

    z = z0

    (outlev>0) && printstyled("\n---delta(pb, z) = $(delta(pb, z))\n", color = :light_cyan)

    it = 0
    while (delta(pb, z) > theta) && (it < 20)
        (outlev>0) && printstyled("\n--------- Iteration i = $it\n", color = :light_cyan)
        (outlev>0) && printstyled("--- delta(pb, z) = $(delta(pb, z))\n", color = :light_cyan)

        mu_z = mu(z)
        dc = NesterovToddstep(pb, z, mu_z)

        (outlev>1) && println("-- ", delta(pb, z), " < ", sqrt( (2*tau^2) / (1+2*tau^2) ))
        if delta(pb, z) < sqrt( (2*tau^2) / (1+2*tau^2) )
            alpha = 1
        else
            i = 0
            alpha = 2^(-i)

            while phi_mu(z+alpha*dc, mu_z) > (phi_mu(z, mu_z) - 4*omega*mu_z*alpha*delta(pb, z)^2)

                i += 1

                (i > 100) && error("i = $i...")
            end
        end
        (outlev > 1) && @show alpha

        z = z+alpha * dc

        it += 1

    end

    (outlev>0) && printstyled("\n--- delta(pb, z) = $(delta(pb, z))\n", color = :green)

    return z
end

function NTpredcorr(pb::SDCOContext, z::PointPrimalDual{T, Dense{T}}; outlev = 1) where T<:Number
    z_curr = z
    epsilon = 1e-15
    maxit = 20
    i = 0

    nc = get_nc(z.x)
    K = ceil((1 + sqrt(1+ 13*nc/2)) * log(mu(z) / epsilon))
    @show K

    while (i < maxit) && (mu(z) > epsilon)
        if outlev > 0
            printstyled("\n--------- Iteration i = $i\n", color = :light_cyan)
            printstyled("Starting point - central point:             ", mu(z), "\n", color=:light_cyan)
            printstyled("Starting point - distance to central path:  ", delta(pb, z), "\n", color=:light_cyan)
        end

        # centering step
        mu_z = mu(z)
        dc = NesterovToddstep(pb, z, mu_z)

        zc = z + dc

        if outlev > 0
            printstyled("\nCentered step - central point:              ", mu(zc), "\n", color=:light_cyan)
            printstyled("Centered step - distance to central path:   ", delta(pb, zc), "\n", color=:light_cyan)
        end

        # affine step
        da = NesterovToddstep(pb, zc, 0.)

        alpha = NT_getalpha(pb, zc, da)

        (outlev > 0) && println("alpha = ", alpha)
        (outlev > 1) && printstyled("mu(z) - mu(zc) : ", abs(mu(z) - mu(zc)), "\n", color=:red)

        z = zc + alpha*da
        (outlev > 1) && printstyled("mu(z+) / mu(zc) - (1 - alpha): ", abs(mu(z) / mu(zc) - (1-alpha)), "\n", color=:red)

        i+=1
    end

    @show mu(z_curr), K
    (outlev > 0) && printstyled("\n\nDone, mu = $(mu(z))\n\n", color=:green)
    (outlev > 1) && print_pointsummary(pb, z)

    return z
end

function NesterovToddstep(pb::SDCOContext{T}, z::PointPrimalDual{T, Dense{T}}, mu::Float64; outlev=0) where T<:Number
    io = stdout

    (outlev >= 1) && println(io, "------ NesterovToddstep ------")

    # Step 1. Compute necessary quatities
    g = NTget_g(pb, z.x, z.s)
    w = NTget_w(pb, g)
    sinv = inv(z.s)
    muxs = z.x - product(sinv, mu)

    # Step 2. Solve M.dy = A(x - \mu inv(s))
    M = NTget_M(pb, g)
    rhs = evaluate(pb.A, muxs)

    # Choleski fact, solve
    dy = cholesky(M) \ rhs

    # Step 3.
    ds = - evaluate(pb.A, dy)

    dx = - muxs - hadamard(w, ds, w)

    # Step 4.
    xstep = z.x+dx
    ystep = z.y+dy
    sstep = z.s+ds

    if outlev >= 1
        @printf("||A*(dy) + ds||                    = %.3e\n", norm(evaluate(pb.A, dy) + ds))
        @printf("||A(dx)||                          = %.3e\n", norm(evaluate(pb.A, dx)))
        @printf("||dx + w.ds.w - mu inv(s) + x||    = %.3e\n", norm(dx + hadamard(w, ds, w) - mu * inv(z.s) + z.x))
        @printf("System actually solved:\n")
        @printf("||A(w . A*(dy) . w) - A(x-muinv(s))|| = %.3e\n", norm(evaluate(pb.A, hadamard(w, evaluate(pb.A, dy), w))
                                                                    - evaluate(pb.A, z.x - mu * inv(z.s))))
        @printf("||A*(dy) + ds||                       = %.3e\n", norm(evaluate(pb.A, dy) + ds))
        @printf("||dx + w.ds.w - mu inv(s) + x||       = %.3e\n", norm(dx + hadamard(w, ds, w) - mu * inv(z.s) + z.x))
    end

    # Check that z = (x, y, s) is strict feasable
    # Check that z+dz is strict feasable too.

    return PointPrimalDual(dx, dy, ds)
end

function NesterovToddstep!(pb::SDCOContext{T}, z::PointPrimalDual{T, Dense{T}}, mu::Float64, d::PointPrimalDual) where T<:Number
    # Step 1. Compute necessary quatities
    g = NTget_g(pb, z.x, z.s)
    w = NTget_w(pb, g)

    muxs = z.x - product(inv(z.s), mu)

    # Step 2. Solve M.dy = A(x - \mu inv(s))
    M = NTget_M(pb, g)
    rhs = evaluate(pb.A, muxs)

    # Choleski fact, solve
    d.y = cholesky(M) \ rhs

    # Step 3.
    d.s = - evaluate(pb.A, d.y)

    d.x = - muxs - hadamard(w, d.s, w)

    nothing
end

function NTget_g(pb::SDCOContext, x::PointE{T, Dense{T}}, s::PointE{T, Dense{T}}; check=true) where T<:Number
    mats = Vector{Dense{T}}([zeros(T, n, n) for n in x.dims])

    ## Deal with SDP part
    for j in eachindex(x.mats)
        xj_fact = LinearAlgebra.cholesky(Symmetric(x.mats[j]), check=check)
        sj_fact = LinearAlgebra.cholesky(Symmetric(s.mats[j]), check=check)

        svd_xs = svd(sj_fact.U * xj_fact.L)

        mats[j] = xj_fact.L * transpose(svd_xs.Vt) * Diagonal( 1 ./ sqrt.(svd_xs.S))
    end

    ## Vectorial part
    vec = sqrt.( sqrt.(x.vec) ./ sqrt.(s.vec) )

    return PointE(mats, vec)
end

function NTget_w(pb::SDCOContext, g::PointE{T, Dense{T}}) where T<:Number
    mats = Vector{Dense{T}}([zeros(T, n, n) for n in g.dims])

    for (i, mat) in enumerate(g.mats)
        mats[i] = mat * transpose(mat)
    end

    vec = g.vec .* g.vec

    return PointE(mats, vec)
end

function NTget_M(pb::SDCOContext, g::PointE{T}) where T<:Number
    m = pb.m
    M = zeros(m, m)

    vectors = [PointE(g.dims, length(g.vec), T, Dense{T}) for i=1:m]
    for i=1:m
        hadamard!(g, pb.A[i], g, vectors[i], transposefirst = true)
    end

    for i=1:m, j=1:m
        M[i, j] = dot(vectors[i], vectors[j])
    end
    return M
end

function NT_getalpha(pb, zc, da)
    g = NTget_g(pb, zc.x, zc.s)
    gtr = transpose(g)
    dasdax = hadamard(da.s, da.x)
    Enorm = norm(hadamard(inv(g), dasdax, g) + hadamard(gtr, dasdax, inv(gtr)))
    mu_zc = mu(zc)

    return 2 / (1 + sqrt(1+13*Enorm/(2*mu_zc)))
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, DocumenterTools, JuliaExamplePkg

makedocs(sitename="JuliaExamplePkg.jl", format = Documenter.HTML(prettyurls = false))

deploydocs(
    repo = "github.com/migarstka/JuliaExamplePkg.jl.git",
)


text:
# Linear Prediction Coefficients (LPC)

function estimate(def::LinearPredictionCoef, x::AbstractArray;
                  use_mgcep::Bool=false, kargs...)
    order = param_order(def)
    data = if use_mgcep
        _mgcep(x, order, 0.0, -1.0; otype=5, kargs...)
    else
        SPTK.lpc(x, order; kargs...)
    end
    SpectralParamState(def, data, false, true)
end

function lpc(x::AbstractArray, order=25; kargs...)
    def = LinearPredictionCoef(order)
    estimate(def, x; kargs...)
end


text:
function save_mesh(::Val{:TECPLOT}, io::IO, m::AirfoilMesh)
    NI, NJ = size(m.grid)
    println(io, """T="Exported from AirfoilMeshers.jl"
    VARIABLES = "X" "Y"
    ZONE T="FLOW" I=$NI J=$NJ DATAPACKING=POINT""")
    for p in m.grid
        @printf io "%24.16E %24.16E\n" p[1] p[2]
    end
    return nothing
end


text:
using DataFrames

# import Base.isequal, Base.show, Base.isless, Base.parse

# Implements code related to atomic shells and subshells
# The SubShell and AtomicSubShell structures

"""
    Shell

Represents the K, L, M, N etc atomic shell
"""
struct Shell
    n::Int
    function Shell(n::Int)
        @assert (n >= 1) && (n <= 10)
        new(n)
    end
end

name(sh::Shell) = "$('J'+sh.n)"
Base.show(io::IO, sh::Shell) = print(io, name(sh))
n(sh::Shell) = sh.n
Base.hash(sh::Shell, h::UInt)::UInt = hash(sh.n, h)
Base.isequal(sh1::Shell, sh2::Shell) = sh1.n == sh2.n
Base.isless(sh1::Shell, sh2::Shell) = sh1.n > sh2.n # In binding energy order

const KShell = Shell(1)
const LShell = Shell(2)
const MShell = Shell(3)
const NShell = Shell(4)

"""
    SubShell

Represents one of the various subshells in an atom.  (See `AtomicSubShell` to combine
the element with a `SubShell`.)
Member data items are index::Int where 1=>K, 2=>L1, ..., 36=>P11.
Construct using SubShell(name::AbstractString) where name = "K", "L1"...,"P11"
"""
struct SubShell
    index::Int
    SubShell(idx::Int) =
        ((idx >= 1) && (idx <= length(subshellnames))) ? new(idx) :
        error("Unknown sub-shell: Index = $(idx)")
end

Base.show(io::IO, ss::SubShell) = print(io, subshellnames[ss.index])

Base.hash(ss::SubShell, h::UInt)::UInt = hash(ss.index, h)
Base.isequal(sh1::SubShell, sh2::SubShell) = sh1.index == sh2.index
Base.isless(sh1::SubShell, sh2::SubShell) = sh1.index > sh2.index # Outer to inner, approx E order

const subshellnames = (
    "K",
    ( "L$i" for i in 1:3 )...,
    ( "M$i" for i in 1:5 )...,
    ( "N$i" for i in 1:7 )...,
    ( "O$i" for i in 1:9 )...,
    ( "P$i" for i in 1:11 )...,
    ( "Q$i" for i in 1:3 )...
)

"""
    allsubshells

A tuple containing all K, L, M, N and O sub-shells
"""
const allsubshells = SubShell.(eachindex(subshellnames))

const K = allsubshells[1]
const L1 = allsubshells[2]
const L2 = allsubshells[3]
const L3 = allsubshells[4]
const M1 = allsubshells[5]
const M2 = allsubshells[6]
const M3 = allsubshells[7]
const M4 = allsubshells[8]
const M5 = allsubshells[9]
const N1 = allsubshells[10]
const N2 = allsubshells[11]
const N3 = allsubshells[12]
const N4 = allsubshells[13]
const N5 = allsubshells[14]
const N6 = allsubshells[15]
const N7 = allsubshells[16]

subshell(idx::Int) = allsubshells[idx]

"""
    capacity(ss::SubShell)

Electron capacity for the specified sub-shell.  This is the total number of electrons
that the sub-shell can hold, not the actual number a specific ground-state element may
have in that sub-shell.
"""
capacity(ss::SubShell) = convert(Int, 2 * j(ss)) + 1

n(ss::SubShell) = (
    1,
    2,
    2,
    2,
    3,
    3,
    3,
    3,
    3,
    4,
    4,
    4,
    4,
    4,
    4,
    4,
    5,
    5,
    5,
    5,
    5,
    5,
    5,
    5,
    5,
    6,
    6,
    6,
    6,
    6,
    6,
    6,
    6,
    6,
    6,
    6,
    7,
    7,
    7,
)[ss.index]

"""
    l(ss::SubShell)

Orbital angular momentum quantum number
"""
l(ss::SubShell) = (
    0,
    0,
    1,
    1,
    0,
    1,
    1,
    2,
    2,
    0,
    1,
    1,
    2,
    2,
    3,
    3,
    0,
    1,
    1,
    2,
    2,
    3,
    3,
    4,
    4,
    0,
    1,
    1,
    2,
    2,
    3,
    3,
    4,
    4,
    5,
    5,
    0,
    1,
    1,
)[ss.index]

"""
     j(ss::SubShell)

Total angular momentum quantum number
"""
j(ss::SubShell) = (
    1 // 2,
    1 // 2,
    1 // 2,
    3 // 2,
    1 // 2,
    1 // 2,
    3 // 2,
    3 // 2,
    5 // 2,
    1 // 2,
    1 // 2,
    3 // 2,
    3 // 2,
    5 // 2,
    5 // 2,
    7 // 2,
    1 // 2,
    1 // 2,
    3 // 2,
    3 // 2,
    5 // 2,
    5 // 2,
    7 // 2,
    7 // 2,
    9 // 2,
    1 // 2,
    1 // 2,
    3 // 2,
    3 // 2,
    5 // 2,
    5 // 2,
    7 // 2,
    7 // 2,
    9 // 2,
    9 // 2,
    11 // 2,
    1 // 2,
    1 // 2,
    3 // 2,
)[ss.index]

function NeXLUncertainties.asa(::Type{DataFrame}, vss::AbstractVector{SubShell})
    css = sort(vss, rev=true)
    return DataFrame(
        SubShell = css,
        Shell = shell.(css),
        n = n.(css),
        𝓁 = l.(css),
        j = j.(css),
        Capacity = capacity.(css)
    )
end

"""
    shell(sh::SubShell)

Returns the appropriate Shell object.

Example:

    julia> shell(n"M5")
    Shell[M]
"""
shell(sh::SubShell)::Shell = Shell(n(sh))

firstsubshell(sh::Shell) = allsubshells[findfirst(ss -> n(ss) == n(sh), allsubshells)]
lastsubshell(sh::Shell) = allsubshells[findlast(ss -> n(ss) == n(sh), allsubshells)]

"""
    ksubshells

All K sub-shells ( K )
"""
const ksubshells = tuple(filter(sh -> n(sh) == 1, collect(allsubshells))...)

"""
    lsubshells

All L sub-shells ( L1, L2, L3 )
"""
const lsubshells = tuple(filter(sh -> n(sh) == 2, collect(allsubshells))...)

"""
    msubshells

All M sub-shells ( M1, M2,.., M5 )
"""
const msubshells = tuple(filter(sh -> n(sh) == 3, collect(allsubshells))...)

"""
    nsubshells

All N sub-shells ( N1, N2,.., N7 ) ]
"""
const nsubshells = tuple(filter(sh -> n(sh) == 4, collect(allsubshells))...)

"""
    osubshells

All O sub-shells  ( O1, O2,.., O9 )
"""
const osubshells = tuple(filter(sh -> n(sh) == 5, collect(allsubshells))...)

const shelltosubshells = Dict{Shell,Tuple{Vararg{SubShell}}}(
    Shell(1) => ksubshells,
    Shell(2) => lsubshells,
    Shell(3) => msubshells,
    Shell(4) => nsubshells,
    Shell(5) => osubshells,
)

# AtomicSubShell functions

"""
    AtomicSubShell

Represents a specific ground-state occupied sub-shell in a specific element.
"""
struct AtomicSubShell
    z::Int
    subshell::SubShell
    function AtomicSubShell(z::Int, ss::SubShell)
        (!hasedge(z, ss.index)) && error("The $(symbol(elements[z])) $(ss) sub-shell is not occupied in the ground state.")
        return new(z, ss)
    end
    AtomicSubShell(elm::Element, ss::SubShell) = AtomicSubShell(z(elm), ss)
end

jumpratio(ass::AtomicSubShell) = jumpratio(ass.z, ass.subshell.index, FFASTDB)

"""
     element(ass::AtomicSubShell)

The Element associated with the specified sub-shell.
"""
element(ass::AtomicSubShell) = elements[ass.z]

Base.isequal(ass1::AtomicSubShell, ass2::AtomicSubShell) =
    (ass1.z == ass2.z) && isequal(ass1.subshell, ass2.subshell)

Base.isless(ass1::AtomicSubShell, ass2::AtomicSubShell) =
    return ass1.z == ass2.z ? isless(ass1.subshell, ass2.subshell) : isless(ass1.z, ass2.z)

Base.show(io::IO, ass::AtomicSubShell) =
    print(io, "$(elements[ass.z].symbol) $(ass.subshell)")

"""
    shell(ass::AtomicSubShell)

Example:

    julia> shell(n"Fe L3")
    Shell(L)
"""
shell(ass::AtomicSubShell) = shell(ass.subshell)

subshell(ass::AtomicSubShell) = ass.subshell
subshellindex(ass::AtomicSubShell) = ass.subshell.index

"""
    atomicsubshell(elm::Element, ss::SubShell)::AtomicSubShell

Construct an AtomicSubShell from from an Element and a SubShell.
"""
atomicsubshell(elm::Element, ss::SubShell) = AtomicSubShell(z(elm), ss)

"""
    has(elm::Element, s::SubShell) =

Is the specified sub-shell occupied by one or more electrons in a ground-state
atom of the specified element?
"""
has(elm::Element, ss::SubShell) = hasedge(z(elm), ss.index)

"""
     atomicsubshells(elm::Element, maxE=1.0e6)::Vector{AtomicSubShell}

 Returns a Vector containing all AtomicSubShell structs associated with the
 specified element with less than the specified energy (in eV).

Example:

    julia> atomicsubshells(n"Fe",1.0e3)
    8-element Array{AtomicSubShell,1}:
     Fe M3
     Fe L3
     Fe M5
     Fe L1
     Fe L2
     Fe M1
     Fe M4
     Fe M2
 """
atomicsubshells(
    elm::Element,
    maxE = 1.0e6,
    eety::Type{<:NeXLAlgorithm} = FFASTDB,
) = AtomicSubShell[ atomicsubshell(elm, SubShell(sh)) #
        for sh in filter(sh->NeXLCore.edgeenergy(z(elm), sh, eety) < maxE, subshellindices(z(elm))) ]

atomicsubshells(
    ss::SubShell,
) = AtomicSubShell[atomicsubshell(elm, ss) for elm in filter(e -> has(e, ss), elements[1:92])]

z(ass::AtomicSubShell) = ass.z
n(ass::AtomicSubShell) = n(ass.subshell)

function NeXLUncertainties.asa(::Type{DataFrame}, vass::AbstractVector{AtomicSubShell})
    cva = sort(vass)
    return DataFrame(
        AtomicSubShell = cva,
        SubShell = subshell.(cva),
        Energy = energy.(cva),
        ICX_U2 = map(a->ionizationcrosssection(a, 2.0*energy(a)), cva),
        JumpRatio = jumpratio.(cva),
        FluorYield = fluorescenceyield.(cva)
    )
end

capacity(ass::AtomicSubShell) = capacity(ass.subshell)

struct Pouchou1991 <: NeXLAlgorithm end

"""
    relativeionizationcrosssection(z::Int, ss::Int, ev::AbstractFloat)
    relativeionizationcrosssection(ass::AtomicSubShell, ev::AbstractFloat, ::Type{Pouchou1991})

An approximate expression based of Pouchou and Pouchoir's 1991 (Green Book) expression
for the ionization crosssection plus an additional factor for sub-shell capacity.

Example:

    > (/)(map(e->NeXLCore.relativeionizationcrosssection(n"Fe K",e),[10.0e3,20.0e3])...)
    0.5982578301818324
"""
function relativeionizationcrosssection(
    ass::AtomicSubShell,
    ev::AbstractFloat,
    ::Type{Pouchou1991},
)
    u = ev / energy(ass)
    ss = ass.subshell.index
    m = ss == 1 ? 0.86 + 0.12 * exp(-(0.2 * ass.z)^2) : (ss <= 4 ? 0.82 : 0.78)
    return capacity(ass.subshell) * log(u) / ((energy(ass)^2) * (u^m))
end
relativeionizationcrosssection(ass::AtomicSubShell, ev::AbstractFloat) =
    relativeionizationcrosssection(ass, ev, Pouchou1991)

struct Bambynek1972 <: NeXLAlgorithm end

"""
    meanfluorescenceyield(elm::Element, sh::Shell)
    meanfluorescenceyield(elm::Element, sh::Shell, ::Type{Bambynek1972})

Mean fluorescence yields from Bambynek in Reviews of Modern Physics 44(4), 1972
"""
function meanfluorescenceyield(elm::Element, sh::Shell, ::Type{Bambynek1972})
    fluorYields = (
        (0, 0, 0),
        (0, 0, 0),
        (0, 0, 0),
        (0.00045, 0, 0),
        (0.00101, 0, 0),
        (0.00198, 0, 0),
        (0.00351, 0, 0),
        (0.00579, 0, 0),
        (0.00901, 0, 0),
        (0.0134, 0, 0),
        (0.0192, 0, 0),
        (0.0265, 0, 0),
        (0.0357, 0, 0),
        (0.0469, 0, 0),
        (0.0603, 0, 0),
        (0.076, 0, 0),
        (0.0941, 0, 0),
        (0.115, 0, 0),
        (0.138, 0, 0),
        (0.163, 0.00067, 0),
        (0.19, 0.00092, 0),
        (0.219, 0.00124, 0),
        (0.249, 0.00163, 0),
        (0.281, 0.0021, 0),
        (0.314, 0.00267, 0),
        (0.347, 0.00335, 0),
        (0.381, 0.00415, 0),
        (0.414, 0.00507, 0),
        (0.446, 0.00614, 0),
        (0.479, 0.00736, 0),
        (0.51, 0.00875, 0),
        (0.54, 0.0103, 0),
        (0.568, 0.0121, 0),
        (0.596, 0.014, 0),
        (0.622, 0.0162, 0),
        (0.646, 0.0186, 0),
        (0.669, 0.0213, 0),
        (0.691, 0.0242, 0),
        (0.711, 0.0274, 0),
        (0.73, 0.0308, 0),
        (0.747, 0.0345, 0),
        (0.764, 0.0385, 0),
        (0.779, 0.0428, 0),
        (0.793, 0.0474, 0),
        (0.806, 0.0523, 0),
        (0.818, 0.0575, 0),
        (0.83, 0.063, 0),
        (0.84, 0.0689, 0),
        (0.85, 0.075, 0),
        (0.859, 0.0815, 0),
        (0.867, 0.0882, 0),
        (0.875, 0.0953, 0),
        (0.882, 0.103, 0),
        (0.888, 0.11, 0),
        (0.895, 0.118, 0),
        (0.9, 0.126, 0),
        (0.906, 0.135, 0.00111),
        (0.911, 0.144, 0.00115),
        (0.915, 0.153, 0.0012),
        (0.92, 0.162, 0.00124),
        (0.924, 0.171, 0.00129),
        (0.927, 0.181, 0.00133),
        (0.931, 0.19, 0.00137),
        (0.934, 0.2, 0.00141),
        (0.937, 0.21, 0.00145),
        (0.94, 0.22, 0.00149),
        (0.943, 0.231, 0.00153),
        (0.947, 0.251, 0.0016),
        (0.945, 0.241, 0.00156),
        (0.95, 0.262, 0.00163),
        (0.952, 0.272, 0.00165),
        (0.954, 0.283, 0.00168),
        (0.956, 0.294, 0.0017),
        (0.957, 0.304, 0.00172),
        (0.959, 0.315, 0.00174),
        (0.96, 0.325, 0.00175),
        (0.962, 0.335, 0.00177),
        (0.963, 0.346, 0.00177),
        (0.964, 0.356, 0.00178),
        (0.966, 0.366, 0.00178),
        (0.967, 0.376, 0.00178),
        (0.968, 0.386, 0.00177),
        (0.969, 0.396, 0.001777),
        (0.97, 0.406, 0.00175),
        (0.971, 0.415, 0.00174),
        (0.972, 0.425, 0.00172),
        (0.972, 0.434, 0.0017),
        (0.973, 0.443, 0.00167),
        (0.974, 0.452, 0.00164),
        (0.975, 0.461, 0.00161),
        (0.975, 0.469, 0.00158),
        (0.976, 0.478, 0.00155),
        (0.9757, 0.3958764, 0.00135),
        (0.9767, 0.40004352, 0.00132),
        (0.9777, 0.40421064, 0.00129),
        (0.9787, 0.40837776, 0.00126),
    )
    return fluorYields[z(elm)][n(sh)]
end
meanfluorescenceyield(elm::Element, sh::Shell) =
    meanfluorescenceyield(elm, sh, Bambynek1972)


text:
using Documenter, Matte

makedocs(
    modules = [Matte],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/favicon.ico"]
    ),
    authors = "Angus Moore",
    sitename = "Matte.jl",
    pages = Any["index.md",
        "Getting started" => "g01-getting-started.md",
        "Guides" => [
            "Introduction to Matte" => "g02-intro.md",
            "Building UIs" => "g03-ui.md",
            "Server logic" => "g04-server.md",
            "Displaying data" => [
                "DataFrames" => "g05a-dataframes.md",
                "Plots" => "g05b-plots.md"
            ],
            "State and side-effects" => [
                "Server-side state" => "g06a-server-side-state.md",
                "Side effects and manually updating the UI" => "g06b-side-effects.md"
            ],
            "Dealing with long-running computations" => "g07-long-running-computations.md"
        ],
        "Reference" => [
            "UI Layouts" => "r01-layouts.md",
            "Input elements" => "r02-input-elements.md",
            "Output elements" => "r03-output-elements.md",
            "Controlling style" => "r04-ui-style.md",
            "Running and creating Matte apps" => "r05-running-creating.md"
        ]
    ],
    checkdocs = :exports
)

deploydocs(
    repo = "github.com/angusmoore/Matte.jl.git",
)


text:
"""
	struct ExtractString{T} <: AbstractExtractor
		n::Int
		b::Int
		m::Int
	end

Represents `String` as `n-`grams (`NGramMatrix` from `Mill.jl`) with base `b` and modulo `m`.
"""
struct ExtractString <: AbstractExtractor
	n::Int
	b::Int
	m::Int
end

ExtractString() = ExtractString(3, 256, 2053)
(s::ExtractString)(v::String) = ArrayNode(Mill.NGramMatrix([v], s.n, s.b, s.m))
(s::ExtractString)(v::Vector{String}) = ArrayNode(Mill.NGramMatrix(v, s.n, s.b, s.m))
(s::ExtractString)(v::AbstractString) = s(String(v))
(s::ExtractString)(v::S) where {S<:Union{Missing, Nothing}} = ArrayNode(Mill.NGramMatrix(missing, s.n, s.b, s.m))
(s::ExtractString)(v::ExtractEmpty) = ArrayNode(Mill.NGramMatrix(Vector{String}(), s.n, s.b, s.m))
(s::ExtractString)(v) = s(missing)
(s::ExtractString)(v::Symbol) = s(String(v))
# todo: add jldoctest here so it's better
"""
	extractscalar(Type{String}, n = 3, b = 256, m = 2053)

represents strings as ngrams with
- `n` (the degree of ngram),
- `b` base of string,
- `m` modulo on index of the token to reduce dimension

	extractscalar(Type{Number}, m = 0, s = 1)

extracts number subtracting `m` and multiplying by `s`
"""
extractscalar(::Type{T}, n = 3, b = 256, m = 2053) where {T<:AbstractString} = ExtractString(n, b, m)

Base.hash(e::ExtractString, h::UInt) = hash((e.n, e.b, e.m), h)
Base.:(==)(e1::ExtractString, e2::ExtractString) = e1.n === e2.n && e1.b === e2.b && e1.m === e2.m

# to be compatible with the number version
extractscalar(::Type{T}, e::Entry) where {T<:AbstractString} = ExtractString()


text:
if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

file = string("reference/smc_cloud_fix=true_version=", ver, ".jld2")
cloud = load(file, "cloud")
split_cloud(file, 2)
rejoined_cloud = join_cloud(file, 2)

@testset "Test split and join clouds" begin
    @test cloud.particles           == rejoined_cloud.particles
    @test SMC.get_vals(cloud)       == SMC.get_vals(rejoined_cloud)
    @test SMC.get_loglh(cloud)      == SMC.get_loglh(rejoined_cloud)
    @test SMC.get_old_loglh(cloud)   == SMC.get_old_loglh(rejoined_cloud)
    @test SMC.get_logpost(cloud)    == SMC.get_logpost(rejoined_cloud)
    @test cloud.ESS                 == rejoined_cloud.ESS
    @test cloud.c                   == rejoined_cloud.c
    @test cloud.stage_index         == rejoined_cloud.stage_index
    @test cloud.total_sampling_time == rejoined_cloud.total_sampling_time
    @test cloud.accept              == rejoined_cloud.accept
    @test cloud.n_Φ                 == rejoined_cloud.n_Φ
    @test cloud.resamples           == rejoined_cloud.resamples
    @test cloud.tempering_schedule  == rejoined_cloud.tempering_schedule
end

rm("reference/smc_cloud_fix=true_version=$(ver)_part1.jld2")
rm("reference/smc_cloud_fix=true_version=$(ver)_part2.jld2")


text:
struct DIIS
    trial_vector::Array{Float64,3}
    residual_vector::Array{Float64,3}
    L::Matrix{Float64}
    depth::Int
    cycle::MVector{1,Int}
    function DIIS(n::Integer, depth::Integer=10)
        trial_vector = zeros(n, n, depth)
        residual_vector =  zeros(n, n, depth)
        L = zeros(depth+1, depth+1)
        L[1, 2:end] .= 1
        L[2:end, 1] .= 1
        new(trial_vector, residual_vector, L, depth, [0])
    end
end

function update!(
    diis::DIIS,
    trial::AbstractMatrix{Float64},
    residual::AbstractMatrix{Float64},
    )
    ptr = mod(diis.cycle[], diis.depth) + 1
    diis.trial_vector[:, :, ptr] = trial
    diis.residual_vector[:, :, ptr] = residual

    n = diis.cycle[] < diis.depth ? ptr : diis.depth
    for i = 1:n
        diis.L[i+1, ptr+1] = diis.residual_vector[:, :, i] ⋅ residual
        diis.L[ptr+1, i+1] = diis.L[i+1, ptr+1]
    end

    coeff = inv(diis.L[1:n+1, 1:n+1])[2:n+1, 1]

    m = LinearAlgebra.checksquare(trial)
    for i = 1:m, j = i:m
        trial[i, j] = 0.0
        for k = 1:n
            trial[i, j] += coeff[k] * diis.trial_vector[i, j, k]
        end
        trial[j, i] = trial[i, j]
    end
    diis.cycle[] += 1
end


text:
module ConstrainedSystems

using LinearMaps
using KrylovKit
using RecursiveArrayTools
#using IterativeSolvers
#using UnPack
using Reexport
@reexport using OrdinaryDiffEq

import OrdinaryDiffEq: OrdinaryDiffEqAlgorithm, alg_order, alg_cache,
                    OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                    initialize!, perform_step!, @muladd, @unpack, constvalue,
                    full_cache, @..

import OrdinaryDiffEq.DiffEqBase: AbstractDiffEqLinearOperator,
                                  DEFAULT_UPDATE_FUNC, has_exp,
                                  AbstractODEFunction, isinplace, numargs

import RecursiveArrayTools: recursivecopy, recursivecopy!, recursive_mean

using LinearAlgebra
import LinearAlgebra: ldiv!, mul!, *, \, I

import Base: size, eltype, *, /, +, -

export SaddleSystem, SaddleVector, state, constraint, aux_state, linear_map
export solvector, mainvector
export SchurSolverType, Direct, CG, GMRES, Iterative

include("vectors.jl")
include("saddlepoint/saddlesystems.jl")
include("saddlepoint/linearmaps.jl")
include("saddlepoint/arithmetic.jl")
include("saddlepoint/testproblems.jl")

include("timemarching/types.jl")
include("timemarching/misc_utils.jl")
include("timemarching/timesaddlesystems.jl")
include("timemarching/algorithms.jl")
include("timemarching/testproblems.jl")

end


text:
"""
```
    AdaptiveEqualization <: AbstractHistogramAdjustmentAlgorithm
    AdaptiveEqualization(; nbins = 256, minval = 0, maxval = 1, rblocks = 8, cblocks = 8, clip = 0.1)

    adjust_histogram([T,] img, f::AdaptiveEqualization)
    adjust_histogram!([out,] img, f::AdaptiveEqualization)
```

Performs Contrast Limited Adaptive Histogram Equalisation (CLAHE) on the input
image. It differs from ordinary histogram equalization in the respect that the
adaptive method computes several histograms, each corresponding to a distinct
section of the image, and uses them to redistribute the lightness values of the
image. It is therefore suitable for improving the local contrast and enhancing
the definitions of edges in each region of an image.

# Details

Histogram equalisation was initially conceived to  improve the contrast in a
single-channel grayscale image. The method transforms the distribution of the
intensities in an image so that they are as uniform as possible [1]. The natural
justification for uniformity is that the image has better contrast  if the
intensity levels of an image span a wide range on the intensity scale. As it
turns out, the necessary transformation is a mapping based on the cumulative
histogram---see [Equalization](@ref) for more details.

A natural extension of histogram equalisation is to apply the contrast
enhancement locally rather than globally [2]. Conceptually, one can imagine that
the process involves partitioning the image into a grid of rectangular regions
and applying histogram equalisation based on the local CDF of each contextual
region. However, to smooth the transition of the pixels from one contextual
region to another,  the mapping of a pixel is not necessarily done soley based
on the local CDF of its contextual region. Rather, the mapping of a pixel may be
interpolated based on the CDF of its contextual region, and the CDFs of the
immediate neighbouring regions.

In adaptive histogram equalisation the image ``\\mathbf{G}`` is partitioned into
``P \\times Q`` equisized submatrices,
```math
\\mathbf{G} =  \\begin{bmatrix}
\\mathbf{G}_{11} & \\mathbf{G}_{12} & \\ldots & \\mathbf{G}_{1C} \\\\
\\mathbf{G}_{21} & \\mathbf{G}_{22} & \\ldots & \\mathbf{G}_{2C} \\\\
\\vdots & \\vdots & \\ldots & \\vdots \\\\
\\mathbf{G}_{R1} & \\mathbf{G}_{R2} & \\ldots & \\mathbf{G}_{RC} \\\\
\\end{bmatrix}.
```

For each submatrix ``\\mathbf{G}_{rc}``, one computes a concomitant CDF, which we
shall denote by ``T_{rc}(G_{i,j})``. In the most general case, we will require
four CDFs
```math
\\begin{aligned}
T_1(v)  & \\triangleq  T_{rc}(G_{i,j}) \\\\
T_2(v)  & \\triangleq  T_{(r+1)c}(G_{i,j}) \\\\
T_3(v)  & \\triangleq  T_{(r+1)(c+1)}(G_{i,j}) \\\\
T_4(v)  & \\triangleq  T_{r(c+1)}(G_{i,j}).
\\end{aligned}
```

In order to determine which particular CDFs will be
used in the interpolation step, it is useful to (i) introduce the function

```math
\\Phi(\\mathbf{G}_{rc}) = \\left(  \\phi_{rc},  \\phi'_{rc}\\right) \\triangleq \\left(rP - \\frac{P}{2}, cQ - \\frac{Q}{2} \\right),
```

(ii) form the sequences  ``\\left(\\phi_{11}, \\phi_{21}, \\ldots, \\phi_{R1} \\right)``
and ``\\left(\\phi'_{11}, \\phi'_{12}, \\ldots, \\phi'_{1C} \\right)``, and (iii) define

```math
\\begin{aligned}
t  & \\triangleq  \\frac{i - \\phi_{r1}}{\\phi_{(r+1)1} - \\phi_{r1} } \\\\
u  & \\triangleq  \\frac{j - \\phi'_{1c}}{\\phi'_{1(c+1)} - \\phi'_{1c} }.
\\end{aligned}
```

#### Case I (Interior)
For a  pixel ``G_{i,j}`` in the range

```math
P - \\frac{P}{2} \\le i \\le RP - \\frac{P}{2}  \\quad \\text{and}  \\quad  Q - \\frac{Q}{2} \\le j \\le CQ - \\frac{Q}{2}.
```
values of ``r`` and ``c`` are implicitly defined by the solution to the inequalities
```math
\\phi_{r1} \\le i < \\phi_{(r+1)1}  \\quad \\text{and}  \\quad  \\phi'_{1c} \\le j < \\phi'_{1(c+1)}.
```

The *bilinearly interpolated* transformation that maps an intensity ``v`` at location ``(i,j)`` in the image
to an intensity ``v'`` is given by [3]

```math
v' \\triangleq \\bar{T}(v)  = (1-t) (1-u)T_1(G_{i,j}) + t(1-u)T_2(G_{i,j}) + tuT_3(G_{i,j}) +(1-t)uT_4(G_{i,j}).
```

#### Case II (Vertical Border)

For a  pixel ``G_{i,j}`` in the range

```math
P - \\frac{P}{2} \\le i \\le RP - \\frac{P}{2}  \\quad \\text{and}  \\quad  1 \\le j < Q - \\frac{Q}{2}  \\; \\cup \\;   CQ - \\frac{Q}{2} < j \\le CQ,
```

``r`` is implicitly defined by the solution to the inequality ``\\phi_{r1} \\le i < \\phi_{(r+1)1}``, while

```math
c = \\begin{cases}
   1, & \\text{if }  \\quad  1 \\le j < Q - \\frac{Q}{2}  \\\\
   C, & \\text{if } \\quad   CQ - \\frac{Q}{2} < j \\le CQ.
\\end{cases}
```

The *linearly interpolated* transformation that maps an intensity ``v`` at location ``(i,j)`` in the image
to an intensity ``v'`` is given by

```math
v' \\triangleq \\bar{T}(v)  = (1-t)T_1(G_{i,j}) + tT_2(G_{i,j}).
```

#### Case III (Horizontal Border)

For a  pixel ``G_{i,j}`` in the range

```math
1 \\le i < P - \\frac{P}{2}  \\;\\cup \\;   RP - \\frac{P}{2} < i \\le RP    \\quad \\text{and}  \\quad  Q - \\frac{Q}{2} \\le j \\le CQ - \\frac{Q}{2},
```

``c`` is implicitly defined by the solution to the inequality ``\\phi'_{1c} \\le j < \\phi'_{1(c+1)}``, while
```math
r = \\begin{cases}
   1, & \\text{if }  \\quad  1 \\le i < P - \\frac{P}{2}  \\\\
   R, & \\text{if } \\quad   RP - \\frac{P}{2} < i \\le RP .
\\end{cases}
```

The *linearly interpolated* transformation that maps an intensity ``v`` at location ``(i,j)`` in the image
to an intensity ``v'`` is given by

```math
v' \\triangleq \\bar{T}(v)  = (1-u)T_1(G_{i,j}) + uT_4(G_{i,j}).
```

#### Case IV (Corners)
For a  pixel ``G_{i,j}`` in the range

```math
1 \\le i < \\frac{P}{2}  \\;\\cup \\; RP - \\frac{P}{2} < i \\le RP   \\quad \\text{and}  \\quad  1 \\le j < CQ -  \\frac{Q}{2} \\; \\cup \\;   CQ - \\frac{Q}{2} < j \\le CQ ,
```
we have
```math
r = \\begin{cases}
   1, & \\text{if }  \\quad  1 \\le i < P - \\frac{P}{2}  \\\\
   R, & \\text{if } \\quad   RP - \\frac{P}{2} < i \\le RP
\\end{cases}
 \\quad \\text{and}  \\quad
c = \\begin{cases}
   1, & \\text{if }  \\quad  1 \\le j < Q - \\frac{Q}{2}  \\\\
   C, & \\text{if } \\quad   CQ - \\frac{Q}{2} < j \\le CQ.
\\end{cases}
```
The transformation that maps an intensity ``v`` at location ``(i,j)`` in the image
to an intensity ``v'`` is given by

```math
v' \\triangleq \\bar{T}(v)  = T_1(G_{i,j}).
```

## Limiting Contrast
An unfortunate side-effect of contrast enhancement is that it has a tendency to
amplify the level of noise in an image, especially when the magnitude of the
contrast enhancement is very high. The magnitude of contrast enhancement is
associated with the gradient of ``T(\\cdot)``, because the  gradient determines the
extent to which consecutive input intensities are stretched across the
grey-level spectrum. One can diminish the level of noise amplification by
limiting the magnitude of the contrast enhancement, that is, by limiting the
magnitude of the gradient.

Since the derivative of ``T(\\cdot)`` is the empirical density ``\\hat{f}_{G}``,
the slope of the mapping function at any input intensity is proportional to the
height of the histogram  ``\\hat{f}_{G}`` at that intensity.  Therefore,
limiting the slope of the local mapping function is equivalent to clipping the
height of the histogram. A detailed description of the  implementation  details
of the clipping process can be found in [2].

# Options

Various options for the parameters of this function are described in more detail
below.

## Choices for `img`

The function can handle a variety of input types. The returned image
depends on the input type.

For coloured images, the input is converted to
[YIQ](https://en.wikipedia.org/wiki/YIQ) type and the Y channel is equalised.
This is the combined with the I and Q channels and the resulting image converted
to the same type as the input.

## Choices for `nbins` in `AdaptiveEqualization`

You can specify the total number of bins in the histogram of each local region.

## Choices for `rblocks` and `cblocks` in `AdaptiveEqualization`

The `rblocks` and `cblocks` specify the number of blocks to divide the input
image into in each direction. By default both values are set to eight.

## Choices for `clip` in `AdaptiveEqualization`

The `clip` parameter must be a value between 0 and 1. It defines an implicit
threshold at which a histogram is clipped. Counts that exceed the threshold are
redistributed as equally as possible so that no bin exceeds the threshold limit.
A value of zero means no clipping, whereas a value of one sets the threshold at
the smallest feasible bin limit. A bin limit is feasible if all bin counts can be
redistributed such that no bin count exceeds the limit. In practice, a `clip` value
of zero corresponds to maximal contrast enhancement, whereas a `clip` value of
one corredponds to minimal contrast enhancement. The default value is `0.1`.

## Choices for `minval` and `maxval` in `AdaptiveEqualization`

If `minval` and `maxval` are specified then intensities are equalized to the range
[`minval`, `maxval`]. The default values are 0 and 1.

# Example

```julia

using TestImages, FileIO, ImageView

img =  testimage("mandril_gray")
imgeq = adjust_histogram(img, AdativeEqualization(nbins = 256, rblocks = 4, cblocks = 4, clip = 0.2))

imshow(img)
imshow(imgeq)
```

# References
1. R. C. Gonzalez and R. E. Woods. *Digital Image Processing (3rd Edition)*.  Upper Saddle River, NJ, USA: Prentice-Hall,  2006.
2. S. M. Pizer, E. P. Amburn, J. D. Austin, R. Cromartie, A. Geselowitz, T. Greer, B. ter Haar Romeny, J. B. Zimmerman and K. Zuiderveld “Adaptive histogram equalization and its variations,” *Computer Vision, Graphics, and Image Processing*, vol. 38, no. 1, p. 99, Apr. 1987. [10.1016/S0734-189X(87)80186-X](https://doi.org/10.1016/s0734-189x(87)80156-1)
3. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.  *Numerical Recipes: The Art of Scientific Computing (3rd Edition)*. New York, NY, USA: Cambridge University Press, 2007.
"""
@with_kw struct AdaptiveEqualization{T₁ <: Union{Real,AbstractGray},
                                        T₂ <: Union{Real,AbstractGray},
                                        T₃ <: Real} <: AbstractHistogramAdjustmentAlgorithm
    nbins::Int = 256
    minval::T₁ = 0.0
    maxval::T₂ = 1.0
    rblocks::Int = 8
    cblocks::Int = 8
    clip::T₃ = 0.1
end

function (f::AdaptiveEqualization)(out::GenericGrayImage, img::GenericGrayImage)
    validate_parameters(f)
    height, width = length.(axes(img))
    # If necessary, resize the image so that the requested number of blocks fit exactly.
    resized_height = ceil(Int, height / (2 * f.rblocks)) * 2 * f.rblocks
    resized_width = ceil(Int, width / (2 * f.cblocks)) * 2 * f.cblocks
    must_resize = (resized_height != height) || (resized_width != width) ? true : false
    if must_resize
        img_tmp = imresize(img, (resized_height, resized_width))
        out_tmp = copy(img_tmp)
    else
        img_tmp = img
        out_tmp = out
    end

    # Determine the relevant pixel coordinates for each block.
    block_height = resized_height ÷  f.rblocks
    block_width =  resized_width ÷  f.cblocks
    r_intervals = [StepRange((r - 1) * block_height + 1, 1, r * block_height) for r = 1:f.rblocks]
    c_intervals = [StepRange((c - 1) * block_width + 1, 1, c * block_width) for c = 1:f.cblocks]

    function construct_hist(roi::AbstractArray)
        edges, histogram = build_histogram(roi, f.nbins, minval = f.minval, maxval = f.maxval)
    end

    function construct_cdf(edges::AbstractArray, histogram::AbstractArray)
        lb = first(axes(histogram, 1))
        ub = last(axes(histogram, 1))
        csum = cumsum(histogram[lb:ub])
        cdf = csum / csum[end]
        return edges, cdf
    end

    # Construct a histogram for each block in the image.
    block_hist = [construct_hist(view(img_tmp, r_intervals[r], c_intervals[c], :)) for r = 1:f.rblocks, c = 1:f.cblocks]

    # Redistribute histogram counts in accordance with the clip weight.
    map!(x-> (x[1], clip_histogram!(x[2], f.clip)), block_hist, block_hist)

    # Construct a CDF for each block in the image.
    block_cdf = [construct_cdf(block_hist[r,c]...) for r = 1:f.rblocks, c = 1:f.cblocks]

    block_centroid_r = [r * block_height - block_height ÷ 2 for r = 1:f.rblocks]
    block_centroid_c = [c * block_width - block_width ÷ 2 for c = 1:f.cblocks]
    intensity_range = (f.minval, f.maxval)
    # Transform pixels using linear and bilinear interpolation of the block CDFs.
    transform_image!(out_tmp, img_tmp, block_centroid_r, block_centroid_c,
                     block_width, block_height, intensity_range,
                     f.cblocks, f.rblocks, block_cdf)
    out .= must_resize ?  imresize(out_tmp, (height, width)) : out_tmp
    return out
end

function validate_parameters(f::AdaptiveEqualization)
    !(0 <= f.clip <= 1) && throw(ArgumentError("The parameter `clip` must be in the range [0..1]."))
    !(1 <= f.rblocks && 1 <= f.cblocks) && throw(ArgumentError("The parameters `rblocks` and `cblocks` must be greater than 0."))
end

function (f::AdaptiveEqualization)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    T = eltype(img)
    yiq = convert.(YIQ, img)
    yiq_view = channelview(yiq)
    #=
       TODO: Understand the cause and solution of this error.
       When I pass a view I run into this error on Julia 1.1.
       ERROR: ArgumentError: an array of type `Base.ReinterpretArray` shares memory with another argument and must
       make a preventative copy of itself in order to maintain consistent semantics,
       but `copy(A)` returns a new array of type `Array{Float64,3}`. To fix, implement:
       `Base.unaliascopy(A::Base.ReinterpretArray)::typeof(A)`
    =#
    #adjust_histogram!(view(yiq_view,1,:,:), f)
    y = comp1.(yiq)
    adjust_histogram!(y, f)
    yiq_view[1, :, :] .= y
    out .= convert.(T, yiq)
end

(f::AdaptiveEqualization)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))

function clip_histogram!(histogram::AbstractArray, clip_weight::Number)
    limit, unfilled_bin_count = determine_threshold(histogram, clip_weight)
    initial_excess = determine_excess!(histogram, limit)
    remaining_excess = perform_initial_redistribution!(histogram, limit, initial_excess, unfilled_bin_count)
    perform_iterative_redistribution!(histogram, limit, remaining_excess)
    return histogram
end

function determine_excess!(histogram::AbstractArray, limit::Number)
    excess = zero(limit)
    for n in eachindex(histogram)
        val = histogram[n]
        if val > limit
            excess = excess + (val - limit)
            histogram[n] = limit
        end
    end
    return excess
end

function determine_threshold(histogram::AbstractArray, clip_weight::Number)
    sorted_hist = sort(collect(histogram), rev = true)
    required_capacity = zeros(axes(sorted_hist))
    available_capacity = zeros(axes(sorted_hist))
    # If we pick an empirical bin count as the actual limit, determine how many
    # counts we will need to move (required_capacity) and how much space do
    # we have in the bins (available_capacity).
    for i = 1:length(sorted_hist)
        required_capacity[i] = sum(sorted_hist[1:i] .- sorted_hist[i])
        available_capacity[i] = sum(sorted_hist[i] .- sorted_hist[i:end])
    end
    # Find the smallest feasible bin limit.
    largest_feasible_index = 1
    for i in eachindex(required_capacity)
        if available_capacity[i] < required_capacity[i]
            largest_feasible_index = i - 1
            break
        end
    end
    # The target bin limit is a convex combination of the highest bin count
    # and the smallest feasible bin count.
    smallest_feasible_limit = sorted_hist[largest_feasible_index]
    target_limit = (1 - clip_weight) * first(sorted_hist) + (clip_weight) * smallest_feasible_limit
    # Set the limit to be the empirical count that is close to the target limit.
    chosen_index = searchsortedlast(sorted_hist, round(Int,target_limit); rev = true)
    threshold = sorted_hist[chosen_index]
    # Knowing how many bins have spare capacity will help us distribute the
    # excess bin counts more efficiently in the perform_initial_redistribution!
    # function.
    unfilled_bin_count  = length(sorted_hist) - chosen_index
    return threshold, unfilled_bin_count
end

function perform_initial_redistribution!(histogram::AbstractArray, limit::Number, excess::Number, N::Integer)
    m = excess ÷ N
    for n in eachindex(histogram)
        val = histogram[n]
        if excess > 0
            if val < limit - m
                histogram[n] = histogram[n] + m
                excess = excess - m
            elseif val < limit
                excess = excess - (limit - val)
                histogram[n] = limit
            end
        end
    end
    return excess
end

function perform_iterative_redistribution!(histogram::AbstractArray, limit::Number, excess::Number)
    while excess > 0
        for n in eachindex(histogram)
            val = histogram[n]
            if excess > 0
                if val < limit
                    excess = excess - 1
                    histogram[n] = histogram[n]  + 1
                end
            end
        end
    end
end

function apply_cdf_transform(val::Union{Real,AbstractGray}, minval::Union{Real,AbstractGray}, maxval::Union{Real,AbstractGray}, edges::AbstractArray, cdf::AbstractArray)
    first_edge = first(edges)
    inv_step_size = 1 / step(edges)
    scale = (maxval - minval) / (cdf[end] - first(cdf))
    if val >= edges[end]
        newval = cdf[end]
    elseif val < first_edge
        newval = first(cdf)
    else
        index = floor(Int, (val - first_edge) * inv_step_size) + 1
        newval = cdf[index]
    end
    # Scale the new intensity value to so that it lies in the range [minval, maxval].
    newval = minval + (newval - first(cdf)) * scale
end

function transform_image!(out, img, block_centroid_r, block_centroid_c, block_width, block_height, intensity_range, cblocks, rblocks, block_cdf)
    height, width = length.(axes(out))
    r₀ = first(block_centroid_r) + 1
    r₁ = last(block_centroid_r) - 1
    c₀ = first(block_centroid_c) + 1
    c₁ = last(block_centroid_c) - 1
    block_dimensions = (block_width, block_height)

    # Interior
    bounds = (r₀:r₁, c₀:c₁)
    block_centroids = (block_centroid_r, block_centroid_c)
    transform_interior!(out, img, bounds, block_centroids, block_dimensions,
                        intensity_range, block_cdf)
    # West
    bounds = (r₀:r₁, 1:(c₀-1))
    block_c_idx = 1
    transform_vertical_strip!(out, img, bounds, block_centroid_r, block_c_idx,
                              block_height, intensity_range, block_cdf)
    # East
    bounds = (r₀:r₁, (c₁+1):width)
    block_c_idx = cblocks
    transform_vertical_strip!(out, img, bounds, block_centroid_r, block_c_idx,
                              block_height, intensity_range, block_cdf)
    # North
    bounds = (1:r₀-1, c₀:c₁)
    block_r_idx = 1
    transform_horizontal_strip!(out, img, bounds, block_centroid_c, block_r_idx,
                                block_width, intensity_range, block_cdf)
    # South
    bounds = (r₁+1:height, c₀:c₁)
    block_r_idx = rblocks
    transform_horizontal_strip!(out, img, bounds, block_centroid_c, block_r_idx,
                                block_width, intensity_range, block_cdf)
    # North-West
    bounds = (1:(r₀ - 1), 1:(c₀ - 1))
    block_r_idx = 1
    block_c_idx = 1
    transform_corner!(out, img, bounds, block_r_idx, block_c_idx,
                      intensity_range, block_cdf)
    # North-East
    bounds = (1:(r₀ - 1), (c₁ + 1):width)
    block_r_idx = 1
    block_c_idx = cblocks
    transform_corner!(out, img, bounds, block_r_idx, block_c_idx,
                      intensity_range, block_cdf)
    # South-West
    bounds = ((r₁ + 1):height, 1:(c₀ - 1))
    block_r_idx = rblocks
    block_c_idx = 1
    transform_corner!(out, img, bounds, block_r_idx , block_c_idx,
                      intensity_range, block_cdf)
    # South-East
    bounds = ((r₁ + 1):height, (c₁ + 1):width)
    block_r_idx = rblocks
    block_c_idx = cblocks
    transform_corner!(out, img, bounds, block_r_idx, block_c_idx,
                      intensity_range, block_cdf)
end

function transform_interior!(out, img, bounds, block_centroids, block_dimensions, intensity_range, block_cdf)
    rows, cols = bounds
    block_centroid_r, block_centroid_c = block_centroids
    block_width, block_height = block_dimensions
    minval, maxval = intensity_range
    inv_block_height = 1 / block_height
    inv_block_width = 1 / block_width
    for r in rows
        for c in cols
            rᵢ = round(Int, r * inv_block_height)
            cᵢ = round(Int, c * inv_block_width)
            t = (r - block_centroid_r[rᵢ]) / block_height
            u = (c - block_centroid_c[cᵢ]) / block_width
            T₁ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
            T₂ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ + 1, cᵢ]...)
            T₃ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ + 1, cᵢ + 1]...)
            T₄ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ + 1]...)
            interpolated_val = (1 - t)*(1 - u)*T₁ + t*(1 - u)*T₂  + t*u*T₃ + (1 - t)*u*T₄
            out[r,c] = eltype(img) <: Integer ? ceil(interpolated_val) : interpolated_val
        end
    end
end

function transform_vertical_strip!(out, img, bounds, block_centroid_r, cᵢ, block_height, intensity_range, block_cdf)
    rows, cols = bounds
    minval, maxval = intensity_range
    inv_block_height = 1 / block_height
    for r in rows
        for c in cols
            rᵢ = round(Int, r * inv_block_height)
            t = (r - block_centroid_r[rᵢ]) / block_height
            T₁ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
            T₂ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ + 1, cᵢ]...)
            interpolated_val = (1 - t)*T₁ + t*T₂
            out[r,c] = eltype(img) <: Integer ? ceil(interpolated_val) : interpolated_val
        end
    end
end

function transform_horizontal_strip!(out, img, bounds, block_centroid_c, rᵢ, block_width, intensity_range, block_cdf)
    rows, cols = bounds
    minval, maxval = intensity_range
    inv_block_width = 1 / block_width
    for r in rows
        for c in cols
            cᵢ = round(Int, c * inv_block_width)
            u = (c - block_centroid_c[cᵢ]) / block_width
            T₁ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
            T₂ = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ + 1]...)
            interpolated_val = (1 - u)*T₁ + u*T₂
            out[r,c] = eltype(img) <: Integer ? ceil(interpolated_val) : interpolated_val
        end
    end
end

function transform_corner!(out, img, bounds, rᵢ, cᵢ, intensity_range, block_cdf)
    rows, cols = bounds
    minval, maxval = intensity_range
    for r in rows
        for c in cols
            val = apply_cdf_transform(img[r,c], minval, maxval, block_cdf[rᵢ, cᵢ]...)
            out[r,c] = eltype(img) <: Integer ? ceil(val) : val
        end
    end
end


text:
struct IOs{T,I}
    path::T
    folder::T
    filename::T
    iplot::I
end

function setup_output(path, folder)
    filename = "file"
    iplot = Int32(0)
    OUT = IOs(path, folder, filename, iplot)
    mkpath(joinpath(path, folder))
    return OUT, iplot
end

function savedata(OUT, Upolar, Ucartesian, T, η, 𝓒, ρ, ε, F, FSE,
    nθ, nr, particle_fields, particle_info, time2save::Float64, ::Val{Isotropic})

    # unpack
    path, folder, filename, iplot =
        OUT.path, OUT.folder, OUT.filename, OUT.iplot

    # output path
    filename = string(filename,"_",iplot)
    final_path = string(joinpath(pwd(),path,folder,filename),".h5")

    # Prepare F
    Fxx = [F[i,j][1,1] for i in axes(F,1), j in axes(F,2) ]
    Fzz = [F[i,j][2,2] for i in axes(F,1), j in axes(F,2) ]
    Fxz = [F[i,j][1,2] for i in axes(F,1), j in axes(F,2) ]
    Fzx = [F[i,j][2,1] for i in axes(F,1), j in axes(F,2) ]

    # Prepare FSE
    a1 = [FSE[i,j].a1 for i in axes(FSE,1), j in axes(FSE,2)]
    a2 = [FSE[i,j].a2 for i in axes(FSE,1), j in axes(FSE,2)]
    x1 = [FSE[i,j].x1 for i in axes(FSE,1), j in axes(FSE,2)]
    x2 = [FSE[i,j].x2 for i in axes(FSE,1), j in axes(FSE,2)]
    y1 = [FSE[i,j].y1 for i in axes(FSE,1), j in axes(FSE,2)]
    y2 = [FSE[i,j].y2 for i in axes(FSE,1), j in axes(FSE,2)]

    # Prepare output variables
    Uθ,Ur = getvelocity(Upolar)
    Ux,Uz = getvelocity(Ucartesian)

    # Unpack particle info
    np = length(particle_info)
    xp = [particle_info[i].CCart.x for i in 1:np]
    zp = [particle_info[i].CCart.z for i in 1:np]
    t = [particle_info[i].t for i in 1:np]
    Tp = particle_fields.T

    # Save output
    h5open(final_path, "w") do file

        # create groups
        MESH = create_group(file, "MESH") # create a group
        VAR = create_group(file, "VAR") # create a group
        PART = create_group(file, "PART") # create a group
        Time = create_group(file, "Time") # create a group

        # time
        Time["t"] = time2save

        # mesh variables
        MESH["nθ"] = nθ # create a scalar dataset inside the group
        MESH["nr"] = nr # create a scalar dataset inside the group

        # physical variables
        save_viscosity!(VAR, η)
        VAR["strain_xx"] = ε.xx
        VAR["strain_zz"] = ε.zz
        VAR["strain_xz"] = ε.xz
        VAR["Ux"] = Ux
        VAR["Uz"] = Uz
        VAR["Utheta"] = Uθ
        VAR["Ur"] = Ur
        VAR["T"] = T
        # VAR["nu"] = η.val
        VAR["density"] = ρ
        VAR["Fxx"] = Fxx
        VAR["Fzz"] = Fzz
        VAR["Fxz"] = Fxz
        VAR["Fzx"] = Fzx
        VAR["a1"] = a1
        VAR["a2"] = a2
        VAR["x1"] = x1
        VAR["x2"] = x2
        VAR["y1"] = y1
        VAR["y2"] = y2

#         # particle variables
#         PART["x"] = xp
#         PART["z"] = zp
#         PART["t"] = t
#         PART["T"] = Tp

    end
end

function save_viscosity!(VAR, η::TemperatureDependantPlastic)
    VAR["nu_node"] = η.node
    VAR["nu_ip"] = η.ip
end

function save_viscosity!(VAR, η)
    VAR["nu_node"] = η.val
end

function savedata(OUT, Upolar, Ucartesian, T, η, 𝓒, ρ, ε, F, FSE, nθ, nr,
    particle_fields, particle_info, time2save, ::Val{Anisotropic})

    # unpack
    path, folder, filename, iplot =
        OUT.path, OUT.folder, OUT.filename, OUT.iplot

    # output path
    filename = string(filename,"_",iplot)
    final_path = string(joinpath(pwd(),path,folder,filename),".h5")

    # Prepare output variables
    Uθ,Ur = getvelocity(Upolar)
    Ux,Uz = getvelocity(Ucartesian)

    # Prepare F
    Fxx = [F[i,j][1,1] for i in axes(F,1), j in axes(F,2) ]
    Fzz = [F[i,j][2,2] for i in axes(F,1), j in axes(F,2) ]
    Fxz = [F[i,j][1,2] for i in axes(F,1), j in axes(F,2) ]
    Fzx = [F[i,j][2,1] for i in axes(F,1), j in axes(F,2) ]

    # Prepare FSE
    a1 = [FSE[i,j].a1 for i in axes(FSE,1), j in axes(FSE,2)]
    a2 = [FSE[i,j].a2 for i in axes(FSE,1), j in axes(FSE,2)]
    x1 = [FSE[i,j].x1 for i in axes(FSE,1), j in axes(FSE,2)]
    x2 = [FSE[i,j].x2 for i in axes(FSE,1), j in axes(FSE,2)]
    y1 = [FSE[i,j].y1 for i in axes(FSE,1), j in axes(FSE,2)]
    y2 = [FSE[i,j].y2 for i in axes(FSE,1), j in axes(FSE,2)]

    # Unpack particle info
    np = length(particle_info)
    xp = [particle_info[i].CCart.x for i in 1:np]
    zp = [particle_info[i].CCart.z for i in 1:np]
    t = [particle_info[i].t for i in 1:np]
    Tp = particle_fields.T

    # Save output
    h5open(final_path, "w") do file

        # create groups
        MESH = create_group(file, "MESH") # create a group
        VAR = create_group(file, "VAR") # create a group
        PART = create_group(file, "PART") # create a group
        Time = create_group(file, "Time") # create a group

        # time
        Time["t"] = time2save

        # mesh variables
        MESH["nθ"] = nθ                    # create a scalar dataset inside the group
        MESH["nr"] = nr                     # create a scalar dataset inside the group

        # physical variables
        save_viscosity!(VAR, η)
        VAR["strain_xx"] = ε.xx
        VAR["strain_zz"] = ε.zz
        VAR["strain_xz"] = ε.xz
        VAR["Ux"] = Ux
        VAR["Uz"] = Uz
        VAR["Utheta"] = Uθ
        VAR["Ur"] = Ur
        VAR["T"] = T
        VAR["nu"] = η.val
        VAR["density"] = ρ
        VAR["nu11"] = 𝓒.η11
        VAR["nu33"] = 𝓒.η33
        VAR["nu55"] = 𝓒.η55
        VAR["nu13"] = 𝓒.η13
        VAR["nu15"] = 𝓒.η15
        VAR["nu35"] = 𝓒.η35

        VAR["Fxx"] = Fxx
        VAR["Fzz"] = Fzz
        VAR["Fxz"] = Fxz
        VAR["Fzx"] = Fzx

        VAR["a1"] = a1
        VAR["a2"] = a2
        VAR["x1"] = x1
        VAR["x2"] = x2
        VAR["y1"] = y1
        VAR["y2"] = y2

#         # particle variables
#         PART["x"] = xp
#         PART["z"] = zp
#         PART["t"] = t
#         PART["T"] = Tp

    end
end

function reloader(fname)
    fid = h5open(fname,"r")

    # -- Load fields
    V = fid["VAR"]
    T = read(V,"T")

    Fxx = read(V["Fxx"])
    Fzz = read(V["Fzz"])
    Fxz = read(V["Fxz"])
    Fzx = read(V["Fzx"])

    F = [ [Fxx[u] Fxz[u]; Fzx[u] Fzz[u] ] for u in CartesianIndices(Fxx)]
    return T, F
end

function load_particles(fname)
    fid = h5open(fname,"r")

    # -- Load fields
    V = fid["PART"]
    px = read(V,"x")
    pz = read(V,"z")

    return px, pz
end

function setup_metrics(gr, path, folder; load = false)
    ScratchNu = ScratchNusselt(gr)
    stats_file = joinpath(pwd(), path, folder, "statistics.txt")
    if !load && isfile(stats_file)
        rm(stats_file)
    end
    return ScratchNu, stats_file
end

function write_stats(U, T, np, gr, Time, ScratchNu, stats_file)
    # V = ones(size(S))
    # A = volume_integral_cartesian(V, gr) # total area
    r_in, r_out = 1.22, 2.22
    Area = π*(r_out^2-r_in^2)
    S = @. @views U[1:2:end]^2+U[2:2:end]^2
    Urms = sqrt(volume_integral_cartesian(S, gr)./ Area) # velocity root mean square
    Tav = volume_integral_cartesian(T, gr)./ Area # volume averaged temperature

    # -- Nusselt number
    Ttop_minus, Tbot_minus = T[ScratchNu.itop_minus], T[ScratchNu.ibot_minus]
    ∇T_top = @. (0.0 - Ttop_minus)/ScratchNu.Δr
    ∇T_bot = @. (Tbot_minus - 1.0)/ScratchNu.Δr
    Nu_top = ScratchNu.A * ScratchNu.Δθ * sum(∇T_top) * r_out
    Nu_bot = ScratchNu.A * ScratchNu.Δθ * sum(∇T_bot) * r_in

    open(stats_file, "a") do io
        println(io, Float32(Time), " ", Float32(Urms), " ", Float32(Tav), " ", Float32(Nu_top), " ", Float32(Nu_bot), " ", np)
    end

    println("Urms = ", Urms, " ⟨T⟩ = ", Tav, " Nusselt top =", Nu_top, " Nusselt bot =", Nu_bot, " #particles =", np)
end

struct Nusselt{T,N}
    A::T
    f::T
    Δr::T
    Δθ::T
    itop_minus::Vector{N}
    ibot_minus::Vector{N}
end

function ScratchNusselt(gr)
    θ, r = gr.θ, gr.r
    r_in, r_out = 1.22, 2.22

    runique = unique(r); sort!(runique)
    itop_minus, ibot_minus = findall(x->x==runique[end-1], r), findall(x->x==runique[2], r)
    θtop_minus = view(θ, itop_minus)

    Δr = runique[2] - runique[1]
    Δθ = θtop_minus[1] - θtop_minus[2]
    f = r_in/r_out
    A = log(f)/(2π*r_out*(1-f))

    Nusselt(A,
            f,
            Δr,
            Δθ,
            itop_minus,
            ibot_minus)

end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
using Random
using MLJ
using SigmaRidgeRegression

id_design = BlockCovarianceDesign([IdentityCovarianceDesign(), IdentityCovarianceDesign()])
ps = [100; 100]
grp = GroupedFeatures(ps)
design = set_groups(id_design, grp)
αs = [1.0; 1.0]

ntest = 20_000
ntrain = 400

ridge_sim = GroupRidgeSimulationSettings(;
    groups = grp,
    ntrain = ntrain,
    ntest = ntest,
    Σ = design,
    response_model = RandomLinearResponseModel(; αs = αs, grp = grp),
)

Random.seed!(1)
sim_res = simulate(ridge_sim)

@test length.(sim_res.resampling_idx[1]) == (ntrain, ntest)

X_train = sim_res.X[sim_res.resampling_idx[1][1], :]
Y_train = sim_res.Y[sim_res.resampling_idx[1][1]]

X_test = sim_res.X[sim_res.resampling_idx[1][2], :]
Y_test = sim_res.Y[sim_res.resampling_idx[1][2]]

λs = [1.0; 2.0]
multiridge =
    MultiGroupRidgeRegressor(; groups = grp, λs = λs, center = false, scale = false)

mse = mse_ridge(
    StatsBase.fit(MultiGroupRidgeRegressor(; groups = grp, λs = λs), X_train, Y_train, grp),
    X_test,
    Y_test,
)

_mach = machine(multiridge, MLJ.table(sim_res.X), sim_res.Y)
_eval = evaluate!(_mach, resampling = sim_res.resampling_idx, measure = l2)

@test _eval.measurement[1] == mse