text:
# doc,jl - api

"""
    create(::AppliSales.Order, invoice_id::String)::UnpaidInvoice

    create(::UnpaidInvoice, ::AppliGeneralLedger.BankStatement)::PaidInvoice

- Create an UnpaidInvoice from an AppliSales.Order.
- Create a PaidInvoice from an UnpaidInvoice and a BankStatement.

@see also [`conv2entry`](@ref)

# Example - create an UnpaidInvoice
```jldoctest
julia> using AppliInvoicing

julia> using AppliSales

julia> orders = AppliSales.process()

julia> invnbr = 1000

julia> invoices = [create(order, "A" * string(global invnbr += 1)) for order in orders]
```

# Example - create a PaidInvoice
```jldoctest
julia> using Dates

julia> using AppliSales

julia> using AppliGeneralLedger

julia> using AppliInvoicing
[ Info: Precompiling AppliInvoicing [3941c6da-33b5-11ea-2884-afa98fed5e3b]

julia> const PATH_CSV = "./bank.csv"
"./bank.csv"

julia> orders = AppliSales.process();

julia> journal_entries_1 = AppliInvoicing.process(orders);

julia> unpaid_invoices = retrieve_unpaid_invoices();

julia> stms = AppliInvoicing.read_bank_statements(PATH_CSV);

julia> journal_entries_2 = AppliInvoicing.process(unpaid_invoices, stms);

julia> paid_invoices = retrieve_paid_invoices();

julia> paid_invoices[1]
AppliInvoicing.PaidInvoice("A1002", AppliInvoicing.MetaInvoice("17583855613870738754", "LS", 2020-03-31T10:53:17.153, "€", 1.0), AppliInvoicing.Header("A1002", "Duck City Chronicals", "1185 Seven Seas Dr", "FL 32830", "Lake Buena Vista", "USA", "DD-001", "Mickey Mouse", "mickey@duckcity.com"), AppliInvoicing.OpentrainingItem("Learn Smiling", 2019-08-30T00:00:00, 1000.0, ["Mini Mouse", "Goofy"], 0.21), AppliInvoicing.BankStatement(2020-01-15, "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0))

julia> cmd = `rm test_invoicing.sqlite`;

julia> run(cmd);
```
"""
function create end

ext: jl


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

using MXNet

export _build_conv_autoencoder
function _build_conv_autoencoder(properties, input_size) #input size is 256, the length of the trace
  batch_size = properties["batch_size"]
  filter_counts = properties["conv_filters"]
  filter_lengths = properties["conv_lengths"]
  pool_sizes = properties["pool_size"]
  assert(length(filter_counts) == length(filter_lengths) == length(pool_sizes))
  pool_type = properties["pool_type"]
  act_type = properties["activation"]
  conv_dropout = properties["conv_dropout"]
  fc_layers = properties["fc"]
  dropout = properties["dropout"]

  # last_conv_data_length = Int(input_size / prod(pool_sizes))

  X = mx.Variable(:data)
  Y = mx.Variable(:label)
  # Convolutional layers
  X = mx.reshape(X, (1, input_size, 1, batch_size))
  for i in 1:length(filter_counts)
    X = conv_layer("conv_$i", X, filter_counts[i], filter_lengths[i], act_type, pool_sizes[i], pool_type, conv_dropout)
  end
  X = mx.Flatten(X)
  # Fully connected layers
  for i in 1:(length(fc_layers)-1)
    X = fc_layer("fc_$i", X, fc_layers[i], act_type, dropout)
  end
  X = fc_layer("latent", X, fc_layers[end], act_type, dropout)
  return _build_conv_decoder(X, Y, properties, input_size)
end

function _build_conv_decoder(full_size, properties, input_size)
  X = mx.Variable(:data)
  Y = mx.Variable(:label)
  return _build_conv_decoder(X, Y, properties, full_size)
end

export _build_conv_decoder
function _build_conv_decoder(X::mx.SymbolicNode, Y::mx.SymbolicNode, properties, full_size)
  batch_size = properties["batch_size"]
  filter_counts = properties["conv_filters"]
  filter_lengths = properties["conv_lengths"]
  assert(length(filter_counts) == length(filter_lengths))
  pool_sizes = convert(Vector{Int64}, properties["pool_size"])
  pool_type = properties["pool_type"]
  act_type = properties["activation"]
  conv_dropout = properties["conv_dropout"]
  fc_layers = properties["fc"]
  dropout = properties["dropout"]
  last_conv_data_length = Int(full_size / prod(pool_sizes))

  for i in (length(fc_layers)-1):-1:1
    X = fc_layer("dec_fc_$i", X, fc_layers[i], act_type, dropout)
  end
  if length(filter_counts) > 0
    X = fc_layer("blow_up", X, filter_counts[end] * last_conv_data_length, act_type, conv_dropout)
        X = mx.reshape(X, (1, last_conv_data_length, Int(filter_counts[end]), batch_size))
  else
    X = fc_layer("blow_up", X, full_size, act_type, conv_dropout)
  end

  # Deconvolutions
  for i in length(filter_counts):-1:2 # goes till 2 inclusive
    X = deconv_layer("deconv_$i", X, filter_counts[i-1], filter_lengths[i], act_type, pool_sizes[i], conv_dropout)
  end
  if length(filter_counts) > 0
    X = deconv_layer("deconv_1", X, 1, filter_lengths[1], nothing, pool_sizes[1], 0)
    # output should have shape (features, batch_size)
    X = mx.Flatten(X, name=:out) # (batch_size, width, height=1)
  end
    # X = mx.FullyConnected(X, num_hidden=input_size, name=:out)
  loss = mx.LinearRegressionOutput(X, Y, name=:loss)
  return loss, X
end

export autoencoder
"""
Parses training and validation data from *data* then calls build() method.
"""
function autoencoder(env::DLEnv, data; id="autoencoder", action::Symbol=:auto, train_key="train", xval_key="xval")
  if action == :auto
    action = decide_best_action(network(env,id))
    info(env, 2, "$id: auto-selected action is $action")
  end

  n = network(env, id)

  if action != :load
    training_waveforms = waveforms(data[:set=>train_key])
    xval_data = data[:set=>xval_key]
    #training_waveforms = waveforms(augment>0 ? augment_data(data[:set=>train_key], copy_n_times=augment, risetime_threshold=risetime_threshold, add_noise=true) : data[:set=>train_key])
    #xval_data = augment>0 ? augment_data(data[:set=>xval_key], copy_n_times=augment, risetime_threshold=risetime_threshold, add_noise=true) : data[:set=>xval_key]
    if eventcount(xval_data) < n["batch_size"]
      n["batch_size"] = eventcount(xval_data)
      info("Cross validation set only has $(eventcount(xval_data)) data points. Adjusting bach size accordingly.")
    end
    xval_waveforms = waveforms(xval_data)
    train_provider = mx.ArrayDataProvider(:data => training_waveforms,
        :label => training_waveforms, batch_size=n["batch_size"])
    eval_provider = mx.ArrayDataProvider(:data => xval_waveforms,
        :label => xval_waveforms, batch_size=n["batch_size"])
  else
    train_provider = nothing
    eval_provider = nothing
  end
  build(n, action, train_provider, eval_provider, _build_conv_autoencoder; verbosity=get_verbosity(env), weight_init=false)
  return n
end

export decoder
function decoder(env::DLEnv, latent_data::EventCollection,
  target_data::EventCollection; id="decoder", action::Symbol=:auto, train_key=:train, xval_key=:xval)

  if action == :auto
    action = decide_best_action(network(env,id))
    info(env, 2, "$id: auto-selected action is $action")
  end

  n = network(env, id)

  train_provider = mx.ArrayDataProvider(:data => waveforms(latent_data[:set=>train_key]),
      :label => waveforms(target_data[:set=>train_key]), batch_size=n["batch_size"])
  eval_provider = mx.ArrayDataProvider(:data => waveforms(latent_data[:set=>xval_key]),
      :label => waveforms(target_data[:set=>xval_key]), batch_size=n["batch_size"])

  full_size = sample_size(target_data)
  build(n, action, train_provider, eval_provider, (p, s) -> _build_conv_decoder(full_size, p, s);
  verbosity=get_verbosity(env), weight_init=false)
  return n
end

export encode
function encode(events::EventLibrary, n::NetworkInfo; log=false)
  log && info("$(n.name): encoding '$(events[:name])'...")
  model = n.model
  model = subnetwork(model.arch, model.arg_params, model.aux_params, "latent", true, n.context)
  batch_size = n["batch_size"]

  result = copy(events)

  if eventcount(events) > 0
        all_waveforms = waveforms(events)
        all_encoded = []
        pred_batch_size = batch_size*10
        # process in batches, else kernel crashes
        batch = zeros(size(all_waveforms,1), pred_batch_size)

        for i in 1:pred_batch_size:size(all_waveforms,2)
            batch_range = i:min(size(all_waveforms,2),(i+pred_batch_size-1))
            batch[:, 1:length(batch_range)] = all_waveforms[:,batch_range]
            provider = mx.ArrayDataProvider(:data => batch, batch_size=batch_size)
            push!(all_encoded, mx.predict(model, provider, verbosity=0)[:,1:length(batch_range)])
        end
        result.waveforms = hcat(all_encoded...)
  else
    result.waveforms = zeros(Float32, 0, 0)
  end

  setname!(result, result[:name]*"_encoded")
  push_classifier!(result, "Autoencoder")
  return result
end

function encode(data::DLData, n::NetworkInfo)
  mapvalues(data, encode, n)
end

export decode
function decode(compact::EventLibrary, n::NetworkInfo, pulse_size; log=false)
  log && info("$(n.name): decoding '$(compact[:name])'...")
  X = mx.Variable(:data)
  Y = mx.Variable(:label) # not needed because no training
  loss, X = _build_conv_decoder(X, Y, n.config, pulse_size)

  model = subnetwork(n.model, mx.FeedForward(loss, context=n.context))

  batch_size=n["batch_size"]
  if eventcount(compact) < batch_size
      waveforms = hcat(compact.waveforms, fill(0, size(compact.waveforms, 1), batch_size-size(compact.waveforms, 2)))
    else
      waveforms = compact.waveforms
  end

  provider = mx.ArrayDataProvider(:data => waveforms, batch_size=batch_size)
  transformed = mx.predict(model, provider, verbosity=0)

  result = copy(compact)
  result.waveforms = transformed
  setname!(result, result[:name]*"_decoded")
  push_classifier!(result, "Autoencoder")
  return result
end

function decode(data::DLData, n::NetworkInfo, pulse_size)
  mapvalues(data, decode, n, pulse_size)
end

export mse
function mse(events::EventLibrary, n::NetworkInfo)
  provider = mx.ArrayDataProvider(:data => events.waveforms,
      :label => events.waveforms, batch_size=n["batch_size"])
  mse_result = eval(n.model, provider, mx.MSE())
  return mse_result[1][2]
end

function mse(data::DLData, n::NetworkInfo)
  mapvalues(data, mse, n)
end

export encode_decode
function encode_decode(events::EventLibrary, n::NetworkInfo)
  batch_size=n["batch_size"]
  provider = padded_array_provider(:data, waveforms(events), batch_size)
  reconst = mx.predict(n.model, provider, verbosity=0)
    if eventcount(events) < batch_size
        reconst = reconst[:,1:eventcount(events)]
    end

  result = copy(events)
  result.waveforms = reconst
  setname!(result, "$(result[:name])_reconst")
  push_classifier!(result, "Autoencoder")
  return result
end

function encode_decode(data::DLData, n::NetworkInfo)
  mapvalues(data, encode_decode, n)
end


text:
# Determine the input scitype of an array of detectors
to_input_scitype(detectors) = MLJ.glb(MLJ.input_scitype.(detectors)...)

# Determine the supported composite detector type
ex_to_eltype(type_symbol) = type_symbol == :Unsupervised ? UnsupervisedDetector :
                                Union{SupervisedDetector, UnsupervisedDetector}

# transform an expression specifying the type (unsupervised or supervised) to the corresponding composite types
# e.g. type_ex = Unsupervised, target_type = :"" would lead to a struct called :UnsupervisedCompositeDetector
# that is a subtype of :UnsupervisedDetectorComposite, if you would add target_type = :Probabilistic, the resulting
# struct would be called :ProbabilisticUnsupervisedCompositeDetector
function ex_to_types(type_ex, target_type::Symbol)
    type_symbol = Symbol(type_ex)
    detector_type = ex_to_eltype(type_symbol)
    composite_type = Symbol(target_type, type_symbol, :Detector, :Composite)
    type_composite = Symbol(target_type, type_symbol, :Composite, :Detector)
    detector_type, composite_type, type_composite
end

# create the struct expression for a scoring composite model from the input symbols
function score_composite(detector_type, composite_type, type_composite)
    return quote
        mutable struct $type_composite{detector_names, input_scitype} <: MLJ.$composite_type
            detectors::Vector{<:$detector_type}
            normalize::Function
            combine::Function

            function $type_composite(detectors, detector_names, normalize, combine)
                input_scitypes = to_input_scitype(detectors)
                new{detector_names, input_scitypes}(detectors, normalize, combine)
            end
        end
    end
end

# create the struct expression for a classifying composite model from the input symbols
function class_composite(detector_type, composite_type, type_composite)
    return quote
        mutable struct $type_composite{detector_names, input_scitype} <: MLJ.$composite_type
            detectors::Vector{<:$detector_type}
            normalize::Function
            combine::Function
            classify::Function

            function $type_composite(detectors, detector_names, normalize, combine, classify)
                input_scitypes = to_input_scitype(detectors)
                new{detector_names, input_scitypes}(detectors, normalize, combine, classify)
            end
        end
    end
end

macro CompositeDetector(type_ex)
    detector_type, composite_type, type_composite = ex_to_types(type_ex, Symbol())
    score_composite(detector_type, composite_type, type_composite)
end

macro ProbabilisticDetector(type_ex)
    detector_type, composite_type, type_composite = ex_to_types(type_ex, :Probabilistic)
    score_composite(detector_type, composite_type, type_composite)
end

macro DeterministicDetector(type_ex)
    detector_type, composite_type, type_composite = ex_to_types(type_ex, :Deterministic)
    class_composite(detector_type, composite_type, type_composite)
end

@CompositeDetector Supervised
@CompositeDetector Unsupervised
@ProbabilisticDetector Supervised
@ProbabilisticDetector Unsupervised
@DeterministicDetector Supervised
@DeterministicDetector Unsupervised

# TODO: enable programmatic construction of multiple detectors without names
const ERR_SPECIFY_NAMES = ArgumentError(
    "When specifying more than one detector, also specify names, as in "*
    "`ProbabilisiticDetector(detector1=KNNDetector(k=1), detector2=KNNDetector(k=5))`. ")
warn_ignore_detector_names(detector, named_detectors) =
    "Wrapping the single detector `$detector`. Ignoring $named_detectors. "

# check and transform the arguments given to composite detectors
function extract_detector_args(args, named_args)
    nt = NamedTuple(named_args)
    length(args) < 2 || throw(ERR_SPECIFY_NAMES)
    if length(args) == 1
        detector_names = (:detector,)
        detectors = [only(args),]
        isempty(nt) || @warn warn_ignore_detector_names(only(args), nt)
    else
        detector_names = keys(nt)
        detectors = collect(nt)
    end
    (detector_names, detectors)
end

# fails if it cannot convert to a vector of detectors
function return_composite(target_symbol::Symbol, detectors, args)
    detector_type = eltype(detectors)
    unsupervised = eval(Symbol(target_symbol, :Unsupervised, :Composite, :Detector))
    supervised = eval(Symbol(target_symbol, :Supervised, :Composite, :Detector))
    return detector_type <: UnsupervisedDetector ? unsupervised(detectors, args...) :
           detector_type <: SupervisedDetector ? supervised(detectors, args...) :
           supervised(convert(Vector{OD.Detector}, detectors), args...) # mixed case
end

"""
    ProbabilisticDetector(unnamed_detectors...;
                          normalize,
                          combine,
                          named_detectors...)

Transform one or more raw detectors into a single probabilistic detector (that returns outlier probabilities).
"""
function ProbabilisticDetector(args...; normalize=scale_minmax, combine=combine_mean, named_detectors...)
    detector_names, detectors = extract_detector_args(args, named_detectors)
    args = (detector_names, normalize, combine)
    return_composite(:Probabilistic, detectors, args)
end

"""
    DeterministicDetector(unnamed_detectors...;
                          normalize,
                          combine,
                          named_detectors...)

Transform one or more raw detectors into a single deterministic detector (that returns inlier and outlier classes).
"""
function DeterministicDetector(args...; normalize=scale_minmax, combine=combine_mean,
                               classify=classify_quantile(DEFAULT_THRESHOLD), named_detectors...)
    detector_names, detectors = extract_detector_args(args, named_detectors)
    args = (detector_names, normalize, combine, classify)
    return_composite(:Deterministic, detectors, args)
end

"""
    CompositeDetector(unnamed_detectors...;
                      normalize,
                      combine,
                      named_detectors...)

Transform one or more raw detectors into a single composite detector (that returns raw outlier scores).
"""
function CompositeDetector(args...; normalize=scale_minmax, combine=combine_mean, named_detectors...)
    detector_names, detectors = extract_detector_args(args, named_detectors)
    args = (detector_names, normalize, combine)
    return_composite(Symbol(), detectors, args)
end

# extend the return values with the training scores
function return_with_scores!(network_mach, model, verbosity, scores_train, X)
    fitresult, cache, report = MLJ.return!(network_mach, model, verbosity)
    report = merge(report, (scores = scores_train(X),))
    return fitresult, cache, report
end

# augment the test scores with the training scores from the fit result (report)
augment(Xs) = mach -> MLJ.node((mach, Xs) -> (mach.report.scores, MLJ.transform(mach, Xs)), mach, Xs)
augment_scores(model, Xs) = augment(Xs).(map(d -> MLJ.machine(d, Xs), getfield(model, :detectors)))
augment_scores(model, Xs, ys) = augment(Xs).(map(d -> MLJ.machine(d, Xs, ys), getfield(model, :detectors)))
apply_to_scores(normalize, combine, augmented_scores) = begin
    _to_scores(scores...) = to_scores(normalize, combine, scores...)
    scores = MLJ.node(_to_scores, augmented_scores...)
end

function unsupervised_transformer(model, Xs)
    augmented_scores = augment_scores(model, Xs)
    scores = apply_to_scores(model.normalize, model.combine, augmented_scores)
    return scores
end

function supervised_transformer(model, Xs, ys)
    augmented_scores = augment_scores(model, Xs, ys)
    scores = apply_to_scores(model.normalize, model.combine, augmented_scores)
    return scores
end

function MLJ.fit(model::ProbabilisticUnsupervisedCompositeDetector, verbosity::Int64, X)
    Xs = MLJ.source(X)
    scores = unsupervised_transformer(model, Xs)
    scores_train, scores_test = first(scores), last(scores)
    network_mach = MLJ.machine(ProbabilisticUnsupervisedDetector(), Xs,
                               predict=to_univariate_finite(scores_test),
                               transform=scores_test)
    return_with_scores!(network_mach, model, verbosity, scores_train, X)
end

function MLJ.fit(model::DeterministicUnsupervisedCompositeDetector, verbosity::Int64, X)
    Xs = MLJ.source(X)
    scores = unsupervised_transformer(model, Xs)
    scores_train, scores_test = first(scores), last(scores)
    classes = MLJ.node(model.classify, scores) |> last
    network_mach = MLJ.machine(DeterministicUnsupervisedDetector(), Xs,
                               predict=to_categorical(classes),
                               transform=scores_test)
    return_with_scores!(network_mach, model, verbosity, scores_train, X)
end

function MLJ.fit(model::ProbabilisticSupervisedCompositeDetector, verbosity::Int64, X, y)
    Xs, ys = MLJ.source(X), MLJ.source(y)
    scores = supervised_transformer(model, Xs, ys)
    scores_train, scores_test = first(scores), last(scores)
    network_mach = MLJ.machine(ProbabilisticSupervisedDetector(), Xs, ys,
                               predict=to_univariate_finite(scores_test),
                               transform=scores_test)
    return_with_scores!(network_mach, model, verbosity, scores_train, X)
end

function MLJ.fit(model::DeterministicSupervisedCompositeDetector, verbosity::Int64, X, y)
    Xs, ys = MLJ.source(X), MLJ.source(y)
    scores = supervised_transformer(model, Xs, ys)
    scores_train, scores_test = first(scores), last(scores)
    classes = MLJ.node(model.classify, scores) |> last
    network_mach = MLJ.machine(DeterministicSupervisedDetector(), Xs, ys,
                               predict=to_categorical(classes),
                               transform=scores_test)
    return_with_scores!(network_mach, model, verbosity, scores_train, X)
end

function MLJ.fit(model::UnsupervisedCompositeDetector, verbosity::Integer, X)
    Xs = MLJ.source(X)
    scores = unsupervised_transformer(model, Xs)
    scores_train, scores_test = first(scores), last(scores)
    network_mach = MLJ.machine(UnsupervisedDetector(), Xs, transform=scores_test)
    return_with_scores!(network_mach, model, verbosity, scores_train, X)
end

function MLJ.fit(model::SupervisedCompositeDetector, verbosity::Integer, X, y)
    Xs, ys = MLJ.source(X), MLJ.source(y)
    scores = supervised_transformer(model, Xs, ys)
    scores_train, scores_test = first(scores), last(scores)
    network_mach = MLJ.machine(SupervisedDetector(), Xs, ys, transform=scores_test)
    return_with_scores!(network_mach, model, verbosity, scores_train, X)
end

ProbabilisticDetectorUnion{detector_names} = Union{
    ProbabilisticUnsupervisedCompositeDetector{detector_names},
    ProbabilisticSupervisedCompositeDetector{detector_names}}

DeterministicDetectorUnion{detector_names} = Union{
    DeterministicUnsupervisedCompositeDetector{detector_names},
    DeterministicSupervisedCompositeDetector{detector_names}}

CompositeDetectorUnion{detector_names} = Union{
    UnsupervisedCompositeDetector{detector_names},
    SupervisedCompositeDetector{detector_names}}

MLJ.input_scitype(::Type{<:Union{
    ProbabilisticUnsupervisedCompositeDetector{N,I},
    DeterministicUnsupervisedCompositeDetector{N,I},
    ProbabilisticSupervisedCompositeDetector{N,I},
    DeterministicSupervisedCompositeDetector{N,I},
    UnsupervisedCompositeDetector{N,I},
    SupervisedCompositeDetector{N,I}
}}) where {N,I} = I

Base.propertynames(::Union{ProbabilisticDetectorUnion{detector_names},
                           CompositeDetectorUnion{detector_names}}) where detector_names =
    tuple(:normalize, :combine, detector_names...)

Base.propertynames(::DeterministicDetectorUnion{detector_names}) where detector_names =
    tuple(:normalize, :combine, :classify, detector_names...)

ERR_NO_PROPERTY(name) = error("CompositeDetector has no property $name")

function Base.getproperty(model::Union{ProbabilisticDetectorUnion{detector_names},
                                       CompositeDetectorUnion{detector_names}}, name::Symbol) where detector_names
    name === :normalize && return getfield(model, :normalize)
    name === :combine && return getfield(model, :combine)
    models = getfield(model, :detectors)
    for j in eachindex(detector_names)
        name === detector_names[j] && return models[j]
    end
    ERR_NO_PROPERTY(name)
end

function Base.getproperty(model::DeterministicDetectorUnion{detector_names}, name::Symbol) where detector_names
    name === :normalize && return getfield(model, :normalize)
    name === :combine && return getfield(model, :combine)
    name === :classify && return getfield(model, :classify)
    models = getfield(model, :detectors)
    for j in eachindex(detector_names)
        name === detector_names[j] && return models[j]
    end
    ERR_NO_PROPERTY(name)
end

function Base.setproperty!(model::Union{ProbabilisticDetectorUnion{detector_names},
                                        CompositeDetectorUnion{detector_names}}, name::Symbol, val) where detector_names
    name === :normalize && return setfield!(model, :normalize, val)
    name === :combine && return setfield!(model, :combine, val)
    idx = findfirst(==(name), detector_names)
    idx isa Nothing || return getfield(model, :detectors)[idx] = val
    ERR_NO_PROPERTY(name)
end

function Base.setproperty!(model::DeterministicDetectorUnion{detector_names}, name::Symbol, val) where detector_names
    name === :normalize && return setfield!(model, :normalize, val)
    name === :combine && return setfield!(model, :combine, val)
    name === :classify && return setfield!(model, :classify, val)
    idx = findfirst(==(name), detector_names)
    idx isa Nothing || return getfield(model, :detectors)[idx] = val
    ERR_NO_PROPERTY(name)
end


text:
module BaseFilters

using Random
using Dates
using DataFrames
using Statistics

using ..Utils
using ..AbsTypes

import ..AbsTypes: fit!, transform!

export fit!,transform!
export OneHotEncoder, Imputer, Wrapper

"""
    OneHotEncoder(Dict(
       # Nominal columns
       :nominal_columns => Int[],

       # Nominal column values map. Key is column index, value is list of
       # possible values for that column.
       :nominal_column_values_map => Dict{Int,Any}()
    ))

Transforms myinstances with nominal features into one-hot form
and coerces the instance matrix to be of element type Float64.

Implements `fit!` and `transform`.
"""
mutable struct OneHotEncoder <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function OneHotEncoder(args::Dict=Dict{Symbol,Any}())
      default_args = Dict{Symbol,Any}(
         :name => "ohe",
         :nominal_columns => Int[],
         # Nominal column values map. Key is column index, value is list of
         # possible values for that column.
         :nominal_column_values_map => Dict{Int,Any}()
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

function OneHotEncoder(name::String; opt...)
   OneHotEncoder(Dict(:name=>name,Dict(pairs(opt))...))
end

function fit!(ohe::OneHotEncoder, myinstances::DataFrame, labels::Vector=[])
  # Obtain nominal columns
  nominal_columns = ohe.model[:nominal_columns]
  if nominal_columns == Int[]
    nominal_columns,_ = find_catnum_columns(myinstances)
  end

  # Obtain unique values for each nominal column
  nominal_column_values_map = ohe.model[:nominal_column_values_map]
  if nominal_column_values_map == Dict{Int,Any}()
    for column in nominal_columns
      nominal_column_values_map[column] = unique(myinstances[:, column])
    end
  end

  # Create model
  ohe.model[:nominal_columns] = nominal_columns
  ohe.model[:nominal_column_values_map] = nominal_column_values_map
end

function transform!(ohe::OneHotEncoder, pinstances::DataFrame)
  myinstances = deepcopy(pinstances)
  nominal_columns = ohe.model[:nominal_columns]
  nominal_column_values_map = ohe.model[:nominal_column_values_map]

  # Create new transformed instance matrix of type Float64
  num_rows = size(myinstances, 1)
  num_columns = (size(myinstances, 2) - length(nominal_columns))
  if !isempty(nominal_column_values_map)
    num_columns += sum(map(x -> length(x), values(nominal_column_values_map)))
  end
  transformed_instances = zeros(Float64, num_rows, num_columns)

  # Fill transformed instance matrix
  col_start_index = 1
  for column in 1:size(myinstances, 2)
    if !in(column, nominal_columns)
      transformed_instances[:, col_start_index] = myinstances[:, column]
      col_start_index += 1
    else
      col_values = nominal_column_values_map[column]
      for row in 1:size(myinstances, 1)
        entry_value = myinstances[row, column]
        entry_value_index = findfirst(isequal(entry_value),col_values)
        if entry_value_index == 0 || entry_value_index == nothing
          @warn "Unseen value found in OneHotEncoder,
                for entry ($row, $column) = $(entry_value).
                Patching value to $(col_values[1])."
          entry_value_index = 1
        end
        entry_column = (col_start_index - 1) + entry_value_index
        transformed_instances[row, entry_column] = 1
      end
      col_start_index += length(nominal_column_values_map[column])
    end
  end

  return transformed_instances |> DataFrame
end

"""
    Imputer(
       Dict(
          # Imputation strategy.
          # Statistic that takes a vector such as mean or median.
          :strategy => mean
       )
    )

Imputes NaN values from Float64 features.

Implements `fit!` and `transform`.
"""
mutable struct Imputer <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function Imputer(args=Dict())
      default_args = Dict{Symbol,Any}(
         :name => "imputer",
         # Imputation strategy.
         # Statistic that takes a vector such as mean or median.
         :strategy => mean
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

function Imputer(name::String;opt...)
   Imputer(Dict(:name=>name,Dict(pairs(opt))...))
end

function fit!(imp::Imputer, myinstances::DataFrame, labels::Vector=[])
   nothing
end

function transform!(imp::Imputer, myinstances::DataFrame)
  new_instances = deepcopy(myinstances)
  strategy = imp.model[:strategy]

  for column in 1:size(myinstances, 2)
    column_values = myinstances[:, column]
    col_eltype = infer_eltype(column_values)

    if <:(col_eltype, Real)
      na_rows = map(x -> isnan(x), column_values)
      if any(na_rows)
        fill_value = strategy(column_values[.!na_rows])
        new_instances[na_rows, column] .= fill_value
      end
    end
  end
  return new_instances |> DataFrame
end

"""
    Wrapper(
       default_args = Dict(
          :name => "ohe-wrapper",
          # Transformer to call.
          :transformer => OneHotEncoder(),
          # Transformer args.
          :transformer_args => Dict()
       )
    )

Wraps around a transformer.

Implements `fit!` and `transform`.
"""
mutable struct Wrapper <: Transformer
  name::String
  model::Dict

  function Wrapper(args=Dict())
    default_args = Dict(
      :name => "ohe-wrapper",
      # Transformer to call.
      :transformer => OneHotEncoder(),
      # Transformer args.
      :transformer_args => Dict()
    )
    cargs=nested_dict_merge(default_args,args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],cargs)
  end
end

function Wrapper(transformer::Machine;opt...)
   Wrapper(Dict(:transfomer => transformer,
                :transformer_args => Dict(pairs(opt)))
          )
end

function fit!(wrapper::Wrapper, myinstances::DataFrame, labels::Vector)
  transformer_args = wrapper.model[:transformer_args]
  transformer = createtransformer(
    wrapper.model[:transformer],transformer_args
  )

  if transformer_args != Dict()
    transformer_args = mergedict(transformer.model, transformer_args)
  end
  fit!(transformer, myinstances, labels)

  wrapper.model[:transformer] = transformer
  wrapper.model[:transformer_args] = transformer_args
end

function transform!(wrapper::Wrapper, myinstances::DataFrame)
  transformer = wrapper.model[:transformer]
  return transform!(transformer, myinstances)
end

"""
    createtransformer(prototype::Transformer, args=Dict())

Create transformer

- `prototype`: prototype transformer to base new transformer on
- `options`: additional options to override prototype's options

Returns: new transformer.
"""
function createtransformer(prototype::Transformer, args=Dict())
  new_args = copy(prototype.model)
  if args != Dict()
    new_args = mergedict(new_args, args)
  end

  prototype_type = typeof(prototype)
  return prototype_type(new_args)
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    iscompatible(model, task)

Check whether or not `model` can be used for
learning `task`.
"""
iscompatible(model::MI.Model, task::LearningTask) = false
iscompatible(model::MI.Model, task::RegressionTask) =
  issupervised(model) && (MI.target_scitype(model) == AbstractVector{Continuous})
iscompatible(model::MI.Model, task::ClassificationTask) =
  issupervised(model) && (MI.target_scitype(model) == AbstractVector{<:Finite})
iscompatible(model::MI.Model, task::ClusteringTask) = !issupervised(model)

"""
    isprobabilistic(model)

Check whether or not `model` is probabilistic.
"""
isprobabilistic(model::MI.Model) = false
isprobabilistic(model::MI.Probabilistic) = true

"""
    issupervised(model)

Check whether or not `model` is supervised.
"""
issupervised(model::MI.Model) = false
issupervised(model::MI.Supervised) = true

"""
    defaultloss(val)
    defaultloss(scitype)

Default loss for value `val` or its scientific type `scitype`.
"""
defaultloss(val) = defaultloss(scitype(val))
defaultloss(::Type{<:Infinite}) = L2DistLoss()
defaultloss(::Type{<:Finite}) = MisclassLoss()


text:
#installs packages
using Pkg
#stats
Pkg.add("StatsFuns")
Pkg.add("StaticArrays")
Pkg.add("StatsBase")

#computation
Pkg.add("Polynomials")
Pkg.add("SpecialFunctions")
Pkg.add("PDMats")
Pkg.add("Distributions")
Pkg.add("Random")
Pkg.add("Sobol")
Pkg.add("Distances")
Pkg.add("LinearAlgebra")
Pkg.add("FastGaussQuadrature")
Pkg.add("SparseArrays")
Pkg.add("JuMP")
Pkg.add("Ipopt")

#data analysis
Pkg.add("Plots")
Pkg.add("PyPlot")
Pkg.add("MATLAB")
Pkg.add("DataFrames")
Pkg.add("CSV")

#timing
Pkg.add("TimerOutputs")

#optimization
Pkg.add("Cubature")
Pkg.add("Optim")
Pkg.add("NLSolversBase")
Pkg.add("Roots")
Pkg.add("GaussianProcesses")

#tools
Pkg.add("ArgParse")

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# test_all_types.jl

@testset "Test All Types" begin
    db = DBInterface.connect(DuckDB.DB)
    con = DBInterface.connect(db)

    df = DataFrame(DBInterface.execute(con, "SELECT * FROM test_all_types()"))

    @test isequal(df.bool, [false, true, missing])
    @test isequal(df.tinyint, [-128, 127, missing])
    @test isequal(df.smallint, [-32768, 32767, missing])
    @test isequal(df.int, [-2147483648, 2147483647, missing])
    @test isequal(df.bigint, [-9223372036854775808, 9223372036854775807, missing])
    @test isequal(
        df.hugeint,
        [-170141183460469231731687303715884105727, 170141183460469231731687303715884105727, missing]
    )
    @test isequal(df.utinyint, [0, 255, missing])
    @test isequal(df.usmallint, [0, 65535, missing])
    @test isequal(df.uint, [0, 4294967295, missing])
    @test isequal(df.ubigint, [0, 18446744073709551615, missing])
    @test isequal(df.float, [-3.4028235f38, 3.4028235f38, missing])
    @test isequal(df.double, [-1.7976931348623157e308, 1.7976931348623157e308, missing])
    @test isequal(df.dec_4_1, [-999.9, 999.9, missing])
    @test isequal(df.dec_9_4, [-99999.9999, 99999.9999, missing])
    @test isequal(df.dec_18_6, [-999999999999.999999, 999999999999.999999, missing])
    @test isequal(
        df.dec38_10,
        [-9999999999999999999999999999.9999999999, 9999999999999999999999999999.9999999999, missing]
    )
    @test isequal(
        df.dec38_10,
        [-9999999999999999999999999999.9999999999, 9999999999999999999999999999.9999999999, missing]
    )
    @test isequal(
        df.dec38_10,
        [-9999999999999999999999999999.9999999999, 9999999999999999999999999999.9999999999, missing]
    )
    @test isequal(df.json, ["🦆🦆🦆🦆🦆🦆", "goose", missing])
    @test isequal(df.small_enum, ["DUCK_DUCK_ENUM", "GOOSE", missing])
    @test isequal(df.medium_enum, ["enum_0", "enum_299", missing])
    @test isequal(df.large_enum, ["enum_0", "enum_69999", missing])
    @test isequal(df.date, [Dates.Date(-5877641, 6, 23), Dates.Date(5881580, 7, 11), missing])
    @test isequal(df.time, [Dates.Time(0, 0, 0), Dates.Time(23, 59, 59, 999, 999), missing])
    @test isequal(df.time_tz, [Dates.Time(0, 0, 0), Dates.Time(23, 59, 59, 999, 999), missing])
    @test isequal(
        df.timestamp,
        [Dates.DateTime(-290308, 12, 22, 0, 0, 0), Dates.DateTime(294247, 1, 10, 4, 0, 54, 775), missing]
    )
    @test isequal(
        df.timestamp_tz,
        [Dates.DateTime(-290308, 12, 22, 0, 0, 0), Dates.DateTime(294247, 1, 10, 4, 0, 54, 775), missing]
    )
    @test isequal(
        df.timestamp_s,
        [Dates.DateTime(-290308, 12, 22, 0, 0, 0), Dates.DateTime(294247, 1, 10, 4, 0, 54, 0), missing]
    )
    @test isequal(
        df.timestamp_ms,
        [Dates.DateTime(-290308, 12, 22, 0, 0, 0), Dates.DateTime(294247, 1, 10, 4, 0, 54, 775), missing]
    )
    @test isequal(
        df.timestamp_ns,
        [Dates.DateTime(1677, 9, 21, 0, 12, 43, 146), Dates.DateTime(2262, 4, 11, 23, 47, 16, 854), missing]
    )
    @test isequal(
        df.interval,
        [
            Dates.CompoundPeriod(Dates.Month(0), Dates.Day(0), Dates.Microsecond(0)),
            Dates.CompoundPeriod(Dates.Month(999), Dates.Day(999), Dates.Microsecond(999999999)),
            missing
        ]
    )
    @test isequal(df.varchar, ["🦆🦆🦆🦆🦆🦆", "goose", missing])
    @test isequal(
        df.blob,
        [
            UInt8[
                0x74,
                0x68,
                0x69,
                0x73,
                0x69,
                0x73,
                0x61,
                0x6c,
                0x6f,
                0x6e,
                0x67,
                0x62,
                0x6c,
                0x6f,
                0x62,
                0x00,
                0x77,
                0x69,
                0x74,
                0x68,
                0x6e,
                0x75,
                0x6c,
                0x6c,
                0x62,
                0x79,
                0x74,
                0x65,
                0x73
            ],
            UInt8[0x5c, 0x78, 0x30, 0x30, 0x5c, 0x78, 0x30, 0x30, 0x5c, 0x78, 0x30, 0x30, 0x61],
            missing
        ]
    )
    @test isequal(df.uuid, [UUID(1), UUID(UInt128(340282366920938463463374607431768211455)), missing])
    @test isequal(df.int_array, [[], [42, 999, missing, missing, -42], missing])
    @test isequal(df.varchar_array, [[], ["🦆🦆🦆🦆🦆🦆", "goose", missing, ""], missing])
    @test isequal(
        df.nested_int_array,
        [[], [[], [42, 999, missing, missing, -42], missing, [], [42, 999, missing, missing, -42]], missing]
    )
    @test isequal(df.struct, [(a = missing, b = missing), (a = 42, b = "🦆🦆🦆🦆🦆🦆"), missing])
    @test isequal(
        df.struct_of_arrays,
        [
            (a = missing, b = missing),
            (a = [42, 999, missing, missing, -42], b = ["🦆🦆🦆🦆🦆🦆", "goose", missing, ""]),
            missing
        ]
    )
    @test isequal(df.array_of_structs, [[], [(a = missing, b = missing), (a = 42, b = "🦆🦆🦆🦆🦆🦆"), missing], missing])
    @test isequal(df.map, [Dict(), Dict("key1" => "🦆🦆🦆🦆🦆🦆", "key2" => "goose"), missing])
end


text:
module ImageProcessing

using Images

export createImage, loadImage, setImageProperties
export restrictResolutionToLevel
export getSize, getSpatialDomain, getPixelSpacing
export regImage

struct regImage
           data::ImageMeta
           voxelsize::Array{Float64, 1}
           shift::Array{Float64, 1}
end

# load image and convert it into gray image
function loadImage(pathToFile;
                   voxelsize = [1.0,1.0], shift = [0.0, 0.0])

    # load image
    image = Images.load(pathToFile)
    #image = permutedims(image, (2,1))

    # convert image to gray image
    image = Images.Gray.(image)

    # convert image data to float (easier access)
    image = convert(Array{Float64,2},image)

    # set image properties

    return regImage(image, voxelsize, shift)

end

function createImage(imageData;
                     voxelsize = [1.0,1.0], shift = [0.0, 0.0])

    # load image
    #image = grayim(copy(imageData))

    return regImage(imageData, voxelsize, shift)

end

function restrictResolutionToLevel(image::regImage,level)

    # copy restrict image, there seems to be a bug in the image copy method, take care of the deepcopy here
    im = deepcopy(image.data)
    maxlevel = 0;
    for l=1:level
        if( (width(im)>2) && (height(im)>2) )
            im = restrict(im)
        end
    end

    voxelsizeNew = zeros(2)
    for xy = 1:2
        voxelsizeNew[xy] = (image.voxelsize[xy] * size(image.data,xy))/size(im,xy)
    end

    return regImage(im,voxelsizeNew,image.shift)

end

function getSize(I::regImage)
  return [height(I.data), width(I.data)]
end

end


text:
"""
    GridWalks

This module simulates random walks on grids. We support simulation of a collection of
agents walking over the same grid graph, simultaneously.
"""
module GridWalks
using Distributions: Sampleable
using LightGraphs: AbstractGraph, grid
using StochasticProcesses
using StochasticProcesses.RandomWalks

export

# methods
gridwalk

# include scripts
include("gridwalks.jl")
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
module Functors

using MacroTools

export @functor, fmap

include("functor.jl")

end # module


text:
import Base: iterate

closedform_systems(rt::RecurrenceTemplate, iter) =
    (ClosedFormTemplate(rt, part) for part in (isnothing(iter) ? partitions(size(rt)) : iter))

# TODO: How to deal with templates where a permutation is useless?
recurrence_systems(vars::Vector{Symbol}, shape::MatrixShape, iter=permutations(vars); kwargs...) =
    (RecurrenceTemplate(copy(perm), shape; kwargs...) for perm in (isnothing(iter) ? permutations(vars) : iter))

synthesis_problems(inv::Invariant, recsys_iter, part_iter; kwargs...) =
    (SynthesisProblem(inv, r, c) for r in recsys_iter for c in closedform_systems(r, part_iter))

synthesis_problems(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape, perm_iter, part_iter; kwargs...) =
    synthesis_problems(inv, recurrence_systems(vars, shape, perm_iter; kwargs...), part_iter)

# ------------------------------------------------------------------------------

function strategy_all(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape; kwargs...)
    @assert shape == UpperTriangular || shape == UnitUpperTriangular || shape == UserSpecific
    synthesis_problems(inv, vars, shape, nothing, nothing; kwargs...)
end

function strategy_permutations(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape, roots::Vector{Int}; kwargs...)
    @assert shape == UpperTriangular || shape == UnitUpperTriangular || shape == UserSpecific
    synthesis_problems(inv, vars, shape, nothing, [roots]; kwargs...)
end

function strategy_partitions(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape; kwargs...)
    rec = RecurrenceTemplate(copy(vars), shape; kwargs...)
    synthesis_problems(inv, [rec], nothing; kwargs...)
end

function strategy_fixed(inv::Invariant, vars::Vector{Symbol}, shape::MatrixShape, roots::Vector{Int}; kwargs...)
    rec = RecurrenceTemplate(vars, shape; kwargs...)
    synthesis_problems(inv, [rec], [roots]; kwargs...)
end

strategy_mixed(inv::Invariant, vars::Vector{Symbol}; kwargs...) =
    Iterators.flatten([
        strategy_all(inv, vars, UnitUpperTriangular; kwargs...),
        strategy_all(inv, vars, UpperTriangular; kwargs...),
        strategy_partitions(inv, vars, FullSymbolic; kwargs...)
    ])

# ------------------------------------------------------------------------------

struct Solutions
    solver::NLSolver
    problem::SynthesisProblem
    maxsol::Number
    timeout::Int

    function Solutions(p::SynthesisProblem; maxsol::Number=1, timeout::Int=2, solver::Type{<:NLSolver}=Z3)
        @assert maxsol >= 1
        new(create_solver(p, solver), p, maxsol, timeout)
    end
end

iterate(ms::Solutions) = iterate(ms, 0)

function iterate(ms::Solutions, state)
    (state === nothing || state >= ms.maxsol) && return
    model, result = solve(ms.problem, ms.solver, timeout=ms.timeout)
    model === nothing && return result, nothing
    next_constraints!(ms.solver, ms.problem, model)
    result, state+1
end

function next_constraints!(s::NLSolver, sp::SynthesisProblem, m::NLModel)
    vs = [Symbol(string(v)) for v in body(sp) if !isconstant(v)]
    cs = Clause([Constraint{NEQ}(:($v - $(m[v]))) for v in vs])
    constraints!(s, ClauseSet(cs))
end

# ------------------------------------------------------------------------------

solutions(strategy; kwargs...) =
    Iterators.flatten(Solutions(problem; kwargs...) for problem in strategy)

models(strategy; kwargs...) =
    Iterators.filter(Absynth.issat, solutions(strategy; kwargs...))

function synth(inv::Invariant; timeout=2, solver=Z3, kwargs...)
    progress = ProgressUnknown("I'm trying hard:")
    strategy = strategy_mixed(inv, program_variables(inv); kwargs...)
    for s in solutions(strategy; maxsol=1, timeout=timeout, solver=solver)
        if issat(s)
            finish!(progress)
            return s.recsystem
        end
        next!(progress)
    end
    finish!(progress)
end


text:
using Distributions
using Turing
using Test
using Random

Random.seed!(125)

x = [1.5 2.0]

@model pmmhtest(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  for i in 1:length(x)
    x[i] ~ Normal(m, sqrt(s))
  end
  s, m
end

# PMMH with Gaussian proposal
GKernel(var) = (x) -> Normal(x, sqrt.(var))
alg = PMMH(1000, SMC(20, :m), MH(1,(:s, GKernel(1))))
chain = sample(pmmhtest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.2
@test mean(chain[:m]) ≈ 7/6 atol=0.1

# PMMH with prior as proposal
alg = PMMH(1000, SMC(20, :m), MH(1,:s))
chain = sample(pmmhtest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.1
@test mean(chain[:m]) ≈ 7/6 atol=0.1

# PIMH
alg = PIMH(1000, SMC(20))
chain = sample(pmmhtest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.15
@test mean(chain[:m]) ≈ 7/6 atol=0.1


text:
# # Dihedral symmetry of the Robinson form

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Symmetry/dihedral_symmetry_of_the_robinson_form.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Symmetry/dihedral_symmetry_of_the_robinson_form.ipynb)
# **Adapted from**: Example 5.4 of [GP04]
#
# [GP04] Gatermann, Karin and Parrilo, Pablo A.
# *Symmetry groups, semidefinite programs, and sums of squares*.
# Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128.

using Test #src

# We start by defining the Dihedral group of order 8.
# This group is isomorphic to the following permutation group:

using PermutationGroups
d = perm"(1, 2, 3, 4)"
c = perm"(1, 3)"
G = PermGroup([c, d])

# We could rely on this isomorphism to define this group.
# However, in order to illustrate how to do symmetry reduction with a custom group,
# we show in this example what should be implemented to define a new group.

import GroupsCore

struct DihedralGroup <: GroupsCore.Group
    n::Int
end

struct DihedralElement <: GroupsCore.GroupElement
    n::Int
    reflection::Bool
    id::Int
end

# Implementing GroupsCore API:

Base.one(G::DihedralGroup) = DihedralElement(G.n, false, 0)

Base.eltype(::Type{DihedralGroup}) = DihedralElement
Base.IteratorSize(::Type{DihedralGroup}) = Base.HasLength()

function Base.iterate(G::DihedralGroup, prev::DihedralElement=DihedralElement(G.n, false, -1))
    if prev.id + 1 >= G.n
        if prev.reflection
            return nothing
        else
            next = DihedralElement(G.n, true, 0)
        end
    else
        next = DihedralElement(G.n, prev.reflection, prev.id + 1)
    end
    return next, next
end

GroupsCore.order(::Type{T}, G::DihedralGroup) where {T} = convert(T, 2G.n)
GroupsCore.gens(G::DihedralGroup) = [DihedralElement(G.n, false, 1), DihedralElement(G.n, true, 0)]

# Base.rand not needed for our purposes here

Base.parent(g::DihedralElement) = DihedralGroup(g.n)
function Base.:(==)(g::DihedralElement, h::DihedralElement)
    return g.n == h.n && g.reflection == h.reflection && g.id == h.id
end

function Base.inv(el::DihedralElement)
    if el.reflection || iszero(el.id)
        return el
    else
        return DihedralElement(el.n, false, el.n - el.id)
    end
end
function Base.:*(a::DihedralElement, b::DihedralElement)
    a.n == b.n || error("Cannot multiply elements from different Dihedral groups")
    id = mod(a.reflection ? a.id - b.id : a.id + b.id, a.n)
    return DihedralElement(a.n, a.reflection != b.reflection, id)
end

Base.copy(a::DihedralElement) = DihedralElement(a.n, a.reflection, a.id)

# optional functions:
function GroupsCore.order(el::DihedralElement)
    if el.reflection
        return 2
    else
        if iszero(el.id)
            return 1
        else
            return div(el.n, gcd(el.n, el.id))
        end
    end
end

# The Robinson form is invariant under the following action of the Dihedral group on monomials:
# The action of each element of the groups is to map the variables `x, y` to:
#
# | id | rotation | reflection |
# |----|----------|------------|
# | 0  | x, y     | y, x       |
# | 1  | -y, x    | -x, y      |
# | 2  | -x, -y   | -y, -x     |
# | 3  | y, -x    | x, -y      |

using SumOfSquares
using DynamicPolynomials
@polyvar x y
struct DihedralAction <: Symmetry.OnMonomials end
import SymbolicWedderburn
SymbolicWedderburn.coeff_type(::DihedralAction) = Float64
function SymbolicWedderburn.action(::DihedralAction, el::DihedralElement, mono::AbstractMonomial)
    if iseven(el.reflection + el.id)
        var_x, var_y = x, y
    else
        var_x, var_y = y, x
    end
    sign_x = 1 <= el.id <= 2 ? -1 : 1
    sign_y = 2 <= el.id ? -1 : 1
    return mono([x, y] => [sign_x * var_x, sign_y * var_y])
end

poly = x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1

# We can verify that `poly` is indeed invariant under the action of each element of the group as follows.

G = DihedralGroup(4)
for g in G
    @show SymbolicWedderburn.action(DihedralAction(), g, poly)
    @test SymbolicWedderburn.action(DihedralAction(), g, poly) == poly #src
end

# We can exploit this symmetry for reducing the problem using the `SymmetricIdeal` certificate as follows:

import CSDP
function solve(G)
    solver = CSDP.Optimizer
    model = Model(solver)
    @variable(model, t)
    @objective(model, Max, t)
    pattern = Symmetry.Pattern(G, DihedralAction())
    con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
    optimize!(model)
    @test value(t) ≈ -3825/4096 rtol=1e-2 #src
    @show value(t)

    g = gram_matrix(con_ref).sub_gram_matrices #src
    @test length(g) == 5                       #src
    @test g[1].basis.polynomials == [x*y^2, x, x^3] #src
    @test g[2].basis.polynomials == [x^2*y, y, y^3] #src
    for i in 1:2                        #src
        I = 3:-1:1                      #src
        Q = g[i].Q[I, I]                #src
        @test size(Q) == (3, 3)         #src
        @test Q[1, 1] ≈  1    rtol=1e-2 #src
        @test Q[1, 2] ≈ -5/8  rtol=1e-2 #src
        @test Q[1, 3] ≈ -1    rtol=1e-2 #src
        @test Q[2, 2] ≈ 25/64 rtol=1e-2 #src
        @test Q[2, 3] ≈  5/8  rtol=1e-2 #src
        @test Q[3, 3] ≈  1    rtol=1e-2 #src
    end #src
    @test length(g[3].basis.polynomials) == 2 #src
    @test g[3].basis.polynomials[1] == 1.0 #src
    @test g[3].basis.polynomials[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
    @test size(g[3].Q) == (2, 2)             #src
    @test g[3].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src
    @test g[3].Q[1, 2] ≈ 0.983 rtol=1e-2  #src
    @test g[3].Q[2, 2] ≈ 1/2 rtol=1e-2 #src
    @test g[4].basis.polynomials == [x * y] #src
    @test size(g[4].Q) == (1, 1)       #src
    @test g[4].Q[1, 1] ≈ 0   atol=1e-2 #src
    @test length(g[5].basis.polynomials) == 1 #src
    @test g[5].basis.polynomials[1] ≈ -(√2/2)x^2 + (√2/2)y^2 #src
    @test size(g[5].Q) == (1, 1)       #src
    @test g[5].Q[1, 1] ≈ 0   atol=1e-2 #src
    for g in gram_matrix(con_ref).sub_gram_matrices
        println(g.basis.polynomials)
    end
end
solve(G)

# We notice that we indeed find `-3825/4096` and that symmetry was exploited.

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
const DATE_PATTERN = r"^([0-3][0-9])([0-1][0-9])([0-9]{2})$"
const TIME_PATTERN = r"^([0-2][0-9])([0-5][0-9])([0-5][0-9])$"
# const DATETIME_PATTERN = r"^([0-3][0-9])([0-1][0-9])([0-9]{2}) ([0-2][0-9])([0-5][0-9])([0-5][0-9])$"

# const THREE_LETTER_CODE_PATTERN = r"^[A-Z0-9]{3}$"
# const MANUFACTURER_CODE_PATTERN = THREE_LETTER_CODE_PATTERN
# const LOGGER_ID_PATTERN = THREE_LETTER_CODE_PATTERN
# const EXTENSION_CODE_PATTERN = THREE_LETTER_CODE_PATTERN


text:
using Test
using DataFlowTasks
using DataFlowTasks.TiledFactorization
using LinearAlgebra
using BenchmarkTools
using VectorizationBase

nc = min(Int(VectorizationBase.num_cores()), Threads.nthreads())

m         = 4000
# create an SPD matrix
A = rand(m,m)
A = (A + adjoint(A))/2
A = A + m*I
nt = Threads.nthreads()

@info "OpenBLAS cholesky"
@info BLAS.get_config()
BLAS.set_num_threads(nc)
@info "NBLAS = $(BLAS.get_num_threads())"
tnat = @btime cholesky!(B) setup=(B=copy(A)) evals=1

@info "pseudotiled cholesky"
BLAS.set_num_threads(1)
@info "NBLAS = $(BLAS.get_num_threads())"
tdf = @btime TiledFactorization.cholesky!(B) setup=(B=copy(A)) evals=1
tfj = @btime TiledFactorization._cholesky_forkjoin!(B) setup=(B=PseudoTiledMatrix(copy(A),150))

DataFlowTasks.TASKCOUNTER[] = 0

# F = TiledFactorization._cholesky_forkjoin!(PseudoTiledMatrix(copy(A),256))
F = TiledFactorization.cholesky(A)

BLAS.set_num_threads(8)
er = norm(F.L*F.U-A,Inf)

F = cholesky!(copy(A))
er_blas = norm(F.L*F.U-A,Inf)

using MKL
@info "MKL cholesky"
@info BLAS.get_config()
BLAS.set_num_threads(nc)
@info "NBLAS = $(BLAS.get_num_threads())"
tnat = @btime cholesky!(B) setup=(B=copy(A)) evals=1

# for m×m tiled matrix, there should be 1/6*(m^3-m) + 1/2*(m^2-m) + m tasks created
@info "Number of tasks = $(DataFlowTasks.TASKCOUNTER[])"
@info "Number of threads = $(nt)"
@info "Number of cores = $(nc)"
@info er,er_blas


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using FactCheck
using Latte

facts("Testing Math Layers") do
    context("Sum") do
        net = Net(8)
        data1, data_value1 = MemoryDataLayer(net, :data1, (4, 4, 2))
        data2, data_value2 = MemoryDataLayer(net, :data2, (4, 4, 2))
        data_value1[:]    = rand(Float32, size(data_value1)...) * 256
        data_value2[:]    = rand(Float32, size(data_value2)...) * 256
        add = AddLayer(:add, net, data1, data2)
        init(net)
        forward(net)
        actual = get_buffer(net, :addvalue)
        expected = data_value1 + data_value2
        @fact actual --> roughly(expected)
    end
end

FactCheck.exitstatus()


text:
export  coeffs,
        nw,
        dim,
        deg,
        multi2uni,
        getentry,
        issymmetric,
        integrate

dim(op::AbstractOrthoPoly) = op.deg + 1
dim(mop::MultiOrthoPoly) = mop.dim

deg(op::AbstractOrthoPoly) = op.deg
deg(mop::MultiOrthoPoly) = mop.deg

"""
```
nw(q::EmptyQuad)
nw(q::AbstractQuad)
nw(opq::AbstractOrthoPoly)
nw(opq::Vector{<:AbstractOrthoPoly})
nw(mop::MultiOrthoPoly)
```
returns nodes and weights in matrix form
"""
nw(quad::EmptyQuad) = Array{Float64}(undef,0,2)

function nw(quad::AbstractQuad)
    [quad.nodes quad.weights]
end

nw(op::AbstractOrthoPoly) = nw(op.quad)

function nw(quads::Vector{<:AbstractQuad})
    nodes = [ quad.nodes for quad in quads]
    weights = [ quad.weights for quad in quads]
    return nodes, weights
end

function nw(ops::Vector{<:AbstractOrthoPoly})
    quad = [ op.quad for op in ops ]
    nw(quad)
end

nw(mop::MultiOrthoPoly) = nw(mop.uni)

"""
```
coeffs(op::AbstractOrthoPoly)
coeffs(op::Vector{<:AbstractOrthoPoly})
coeffs(mop::MultiOrthoPoly)
```
returns recurrence coefficients of in matrix form
"""
function coeffs(op::AbstractOrthoPoly)
    [op.α op.β]
end

function coeffs(op::Vector{<:AbstractOrthoPoly})
    a = [ p.α for p in op]
    b = [ p.β for p in op]
    return a,b
end

coeffs(mop::MultiOrthoPoly) = coeffs(mop.uni)

"""
```
integrate(f::Function,nodes::Vector{<:Real},weights::Vector{<:Real})
integrate(f::Function,q::AbstractQuad)
integrate(f::Function,opq::AbstractOrthoPoly)
```
integrate function `f` using quadrature rule specified via `nodes`, `weights`.
For example ``\\int_0^1 6x^5 = 1`` can be solved as follows:

```@repl
julia> opq = Uniform01OrthoPoly(3) # a quadrature rule is added by default
julia> integrate(x -> 6x^5, opq)
0.9999999999999993
```

!!! note
- function ``f`` is assumed to return a scalar.
- interval of integration is "hidden" in `nodes`.
"""
function integrate(f::Function, nodes::Vector{<:Real}, weights::Vector{<:Real})
    dot(weights, f.(nodes))
end

function integrate(f::Function, quad::AbstractQuad)
    typeof(quad) == EmptyQuad && throw(DomainError(quad, "supplied an empty quadrature"))
    integrate(f, quad.nodes, quad.weights)
end

integrate(f::Function, op::AbstractOrthoPoly) = integrate(f, op.quad)

"""
```
issymmetric(m::AbstractMeasure)
issymmetric(op::AbstractOrthoPoly)
```
Is the measure symmetric (around any point in the domain)?
"""
issymmetric(m::AbstractMeasure) = m.symmetric
issymmetric(op::AbstractOrthoPoly) = issymmetric(op.measure)

function multi2uni(a::Vector{<:Int}, ind::Matrix{<:Int})
    minimum(a) < 0 && throw(DomainError(a, "no negative degrees allowed"))
    l, p = size(ind) # p-variate basis
    m = length(a) # dimension of scalar product
    l -= 1 # (l+1)-dimensional basis
    maximum(a) > l && throw(DomainError(a, "not enough elements in multi-index (requested: $(maximum(a)), max: $l)"))
    A = zeros(Int64,p,m)
    for (i, a_) in enumerate(a)
        A[:,i] = ind[a_+1,:]
    end
    return A
end

function getentry(a::Vector{<:Int}, T::SparseVector{<:Real,<:Int}, ind::Matrix{<:Int}, dim::Int)
    m = length(a)
    l = size(ind,1)-1
    minimum(a) < 0 && throw(DomainError(a, "no negative degrees allowed"))
    maximum(a) > l && throw(DomainError(a, "not enough elements in multi-index (requested: $(maximum(a)), max: $l)"))
    m != dim && throw(DomainError(m, "length $m of provided index $a is inconsistent with dimension $(dim) of multi-index"))
    # a .+= 1
    sort!(a)

    sp_ind = 1 + reduce(+,[idx*l^(m-i) for (i,idx) in enumerate(a)])
    return T[sp_ind]
end


text:
using Pkg
Pkg.activate(".")
using PlotlyJS, LsqFit, CSV, DataFrames, Dates, HTTP

# Timeframe in days to use for fitting the curve
timeframe = 60

# Download CSV File
println("Downloading new data...")
table = CSV.File(HTTP.get("https://impfdashboard.de/static/data/germany_vaccinations_timeseries_v2.tsv").body)
println("Done!")
df = table |> DataFrame # Save to DataFrame

# Extract dates and number of vaccinations
dates = df[!,"date"]
dates_numeric = Dates.value.(dates - dates[1])
vaccinations = df[!,"personen_erst_kumulativ"]

# Construct including future dates
dates_future = copy(dates)
days_in_future = 1:90;
append!(dates_future, [dates[end] + Day(i) for i in days_in_future])
dates_future_numeric = Dates.value.(dates_future - dates_future[1])

# Create exponential model function
# @. model(t, p) = p[1] * exp.(p[2] * t) .+ p[3]
# @. model(t, p) = p[1] * exp(p[2] * t) * t * t
# @. model(t, p) = p[1] * (t)^3 + p[3]
function model(t, p)
    p[1] .* t .+  p[2];
end

# Think of some random start value, TODO: improve to use more realistic values
# p0 = [2.0, 2.0, 10000.0]
p0 = [35e3, 2.0]

# Fit the model
fit = curve_fit(model, dates_numeric[end-timeframe:end], vaccinations[end-timeframe:end], p0)

# Get the model parameters
params_model = fit.param

# extrapolate
vax_extrapolated = model(dates_future_numeric, params_model)

# estimate vaccinations per day (1st derivative)
vax_per_day =
    params_model[1] .+ 0.0 .* dates_future_numeric;

# find index where model predicts >80e6 vaccinations
#markline_idx = findfirst(x -> x > 73e6, vax_extrapolated) # Vaccination rate too low at the moment, doesn't work (Sept 2021)

# Use plotly
# plotlyjs()

# Make background transparent for export
if isinteractive()
    bg_col = "#FFF"
else
    bg_col = "rgba(0,0,0,0)"
end

# plot the curve!

layout = Layout(
    title=attr(
        text="Vaccination Model Germany " *"
Last update: " * string(Dates.format(now(), "Y-m-d, HH:MM")),
    ),
    xaxis=attr(
        title=attr(text="Date"),
        showgrid=true,
        zeroline=true,
        rangeslider=attr(
            # visible=true,
            visible=false,
            yaxis=attr(rangemode="auto"),
        ),
        range=[dates[end-timeframe] - Day(10), dates_future[end] + Day(10)],
    ),
    yaxis=attr(
        title=attr(
            text="Number of People with min. 1 Dose",
            standoff=10),
        zeroline=false,
        automargin=true,
        rangemode="tozero",
        scaleanchor="y2",
        scaleratio=1 / 10,
        range=[0, maximum(vax_extrapolated * 1.1)]
    ),
    yaxis2=attr(
        title=attr(
            text="Vaccinations per Day",
            standoff=10),
        overlaying="y",
        side="right",
        scaleanchor="y",
        scaleratio=10,
        constraintoward="bottom",
        rangemode="tozero",

    ),
    paper_bgcolor=bg_col,
    plot_bgcolor=bg_col,
    legend=attr(
        orientation="h",
        x="0.5",
        xanchor="center",
    ),
    hovermode="x unified",
    hoverlabel=attr(namelength=-1)
);

p_actual = scatter(
    x=dates,
    y=vaccinations,
    mode="markers",
    name="Actual Vaccinations (Cumulative)",
    marker=attr(symbol="circle-open"),
    marker_color="blue"
)

p_model = scatter(
    x=dates_future,
    y=vax_extrapolated,
    linewidth=0.5,
    name="Model (a*x + b)",
    marker_color="dodgerblue"
)

p_vpd = scatter(
    x=dates,
    y=df.dosen_erst_differenz_zum_vortag,
    mode="markers",
    name="Actual Vaccinations per Day",
    marker=attr(symbol="diamond-open"),
    yaxis="y2",
    marker_color="red"
)

p_model_vpd = scatter(
    x=dates_future,
    y=vax_per_day,
    linewidth=0.5,
    name="Vaccinations per Day (Model)",
    yaxis="y2",
    marker_color="crimson"
)

# Plot horizontal and vertical line for >73 Mio
#p_lines = scatter(
#    x=vcat(dates_future[1:markline_idx], dates_future[markline_idx]),
#    y=vcat(vax_extrapolated[markline_idx] * ones(markline_idx), 0),
#    name=">73-Mio First-Dose",
#    marker_color="green"
#)

p = plot(
    [p_actual, p_model, p_vpd, p_model_vpd],#, p_lines],
    layout,
    options=Dict(:responsive => true),
    )

if isinteractive()
    display(p)
end

# Save the plot to a file
savehtml(p, "Plot_Vax.html", :remote)

# Create model params file
res_file = open("results.html", "w");
#println(res_file, "
>73 Mio First-Dose Vaccinations on ", dates_future[markline_idx], "
")
println(res_file, "Model Parameters:
")
println(res_file, "
a = ", params_model[1], "
")
println(res_file, "
b = ", params_model[2], "
")
#println(res_file, "
c = ", params_model[3], "
")
#println(res_file, "
d = ", params_model[4], "
")
println(res_file, "
")

close(res_file)

max_stars_repo_name: JuliDi/vax-mod-ger

max_issues_repo_name: JuliDi/vax-mod-ger

max_forks_repo_name: JuliDi/vax-mod-ger


text:
import Dierckx #for cubic spline interplation

function findLocalMaxima(signal::Vector; angle=0.0, includeEdge=:false)::Tuple{Array{Int64,1},Array{Float64,1}}
  Inds = Int64[]
  signal = real.(exp(-1im*angle)*signal)
  if includeEdge #include maxima on edges
    if length(signal)>1
      if signal[1]>signal[2]
        push!(Inds,1)
      end
      for i=2:length(signal)-1
        if signal[i-1]<signal[i]>signal[i+1]
          push!(Inds,i)
        end
      end
      if signal[end]>signal[end-1]
        push!(Inds,length(signal))
      end
    end
  else #ignore maxima on edges
    for i=2:length(signal)-1
      if signal[i-1]<signal[i]>signal[i+1]
        push!(Inds,i)
      end
    end
  end
  Vals = [ signal[idx[1]] for idx in Inds]
  return IndsVals = (Inds, Vals)
end

function SIFT(r::Vector{Float64}; siftStopThresh=1e-6, includeEdge=:true)
#References: Huang, Norden E., et al. "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis." Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences 454.1971 (1998): 903-995.
  L = length(r)
  eStart = round(Int,L/4); eStop = round(Int,3L/4); eL = length(eStart:eStop)
  while true
    maxIndsVals = findLocalMaxima(r, includeEdge=includeEdge)
    if length(maxIndsVals[1]) < 4; break; end
    maxInterpolator = Dierckx.Spline1D(maxIndsVals[1], maxIndsVals[2])
    u = maxInterpolator(1:L)
    minIndsVals = findLocalMaxima(-r, includeEdge=includeEdge)
    if length(minIndsVals[1]) < 4; break; end
    minInterpolator = Dierckx.Spline1D(minIndsVals[1], -minIndsVals[2])
    l = minInterpolator(1:L)
    e = (u+l)./2
    r -= e
    mav = sum(abs,e[eStart:eStop])/eL
    if mav < siftStopThresh; break; end
  return r
  end
end

function EMD(x::Vector{Float64}; emdStopThresh = 1e-2, siftStopThresh = 1e-6, includeEdge = :true)::Array{Vector{Float64},1}
#References: Huang, Norden E., et al. "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis." Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences 454.1971 (1998): 903-995.
  IMF = Vector{Float64}[]
  L = length(x); M = maximum(abs.(x))
  eStart = round(Int,L/4); eStop = round(Int,3L/4); eL = length(eStart:eStop)
  while true
    φ = SIFT(x, siftStopThresh=siftStopThresh, includeEdge=includeEdge)
    if isnothing(φ); break; end
    if maximum(abs.(φ)) > 3M; break; end
    x -= φ
    push!(IMF, φ)
    mav = sum(abs,x[eStart:eStop])/eL
    if mav < emdStopThresh; break; end
  end
  return push!(IMF, x)
end

function ℂSIFT(r::Vector{ComplexF64}; D = 4, siftStopThresh = 1e-3, includeEdge = :true)
#References: Rilling, Gabriel, et al. "Bivariate empirical mode decomposition." IEEE signal processing letters 14.12 (2007): 936-939.
  L = length(r)
  eStart = round(Int,L/4); eStop = round(Int,3L/4); eL = length(eStart:eStop)
  err = false
  while true
    e = zeros(ComplexF64,L)
    for d in 1:2D
      ϑ = 2π*(d-1)/2D
      maxIndsVals = findLocalMaxima(r, angle=ϑ, includeEdge=includeEdge)
      if length(maxIndsVals[1]) < 4; err = true; break; end
      maxInterpolator = Dierckx.Spline1D(maxIndsVals[1], maxIndsVals[2])
      e += exp(1im*ϑ)*maxInterpolator(1:L)
    end
    if err; break; end
    e /= 2D
    r -= e
    mav = sum(abs,e[eStart:eStop])/eL
    if mav < siftStopThresh; break; end
  end
  return r
end

function ℂEMD(z::Vector{ComplexF64}; emdStopThresh=1e-2, D=4, siftStopThresh=1e-3, includeEdge=:true)::Array{Vector{ComplexF64},1}
#References: Rilling, Gabriel, et al. "Bivariate empirical mode decomposition." IEEE signal processing letters 14.12 (2007): 936-939.
  IMF = Vector{ComplexF64}[]
  L = length(z); M = maximum(abs.(z))
  eStart = round(Int,L/4); eStop = round(Int,3L/4); eL = length(eStart:eStop)
  while true
    φ = ℂSIFT(z, D=D, siftStopThresh=siftStopThresh, includeEdge=includeEdge)
    if maximum(abs.(φ)) > 3M; break; end
    z -= φ
    push!(IMF, φ)
    mav = sum(abs,z[eStart:eStop])/eL
    if mav < emdStopThresh; break; end
  end
  return push!(IMF, z)
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
# this implements the 1D p-laplacian (p = 4)
#      F(u) = ∑_{i=1}^{N} h (W(u_i') - ∑_{i=1}^{N-1} h u_i
#  where u_i' = (u_i - u_{i-1})/h
#  plap: implements the functional without boundary condition
#  preconditioner is a discrete laplacian, which defines a metric
#     equivalent (in the limit h → 0) to that induced by the hessian, but
#     does not approximate the hessian explicitly.
@testset "Preconditioning" begin
    plap(U; n=length(U)) = (n-1) * sum((0.1 .+ diff(U).^2).^2) - sum(U) / (n-1)
    plap1(U; n=length(U), dU = diff(U), dW = 4 .* (0.1 .+ dU.^2) .* dU) =
                            (n - 1) .* ([0.0; dW] .- [dW; 0.0]) .- ones(n) / (n-1)
    precond(x::Vector) = precond(length(x))
    precond(n::Number) = Optim.InverseDiagonal(diag(spdiagm(-1 => -ones(n-1), 0 => 2*ones(n), 1 => -ones(n-1)) * (n+1)))
    f(X) = plap([0;X;0])
    g!(G, X) = copyto!(G, (plap1([0;X;0]))[2:end-1])

    GRTOL = 1e-6

    debug_printing && println("Test a basic preconditioning example")
    for N in (10, 50, 250)
        debug_printing && println("N = ", N)
        initial_x = zeros(N)
        Plap = precond(initial_x)
        ID = nothing
        for optimizer in (GradientDescent, ConjugateGradient, LBFGS)
            for (P, wwo) in zip((ID, Plap), (" WITHOUT", " WITH"))
                results = Optim.optimize(f, g!, copy(initial_x),
                                         optimizer(P = P),
                                         Optim.Options(g_tol = GRTOL, allow_f_increases = true, iterations=250000))
                debug_printing && println(optimizer, wwo,
                                          " preconditioning : g_calls = ", Optim.g_calls(results),
                                          ", f_calls = ", Optim.f_calls(results))
                if (optimizer == GradientDescent) && (N > 15) && (P == ID)
                    debug_printing && println("    (gradient descent is not expected to converge)")
                else
                    @test Optim.converged(results)
                end
            end
        end
    end
end


text:
using Test
using LinearAlgebra, SparseArrays, Random

using LuxurySparse, StaticArrays
import LuxurySparse: staticize
using StaticArrays: SVector, SMatrix

Random.seed!(2)

@testset "staticize" begin
    # permmatrix
    m = pmrand(ComplexF64, 4)
    sm = m |> staticize
    @test sm isa SPermMatrix{4, ComplexF64}
    @test sm.perm isa SVector
    @test sm.vals isa SVector
    @test sm.perm == m.perm
    @test sm.vals == m.vals
    dm = sm |> dynamicize
    @test dm isa PermMatrix{ComplexF64}
    @test dm.perm isa Vector
    @test dm.vals isa Vector
    @test dm.perm == m.perm
    @test dm.vals == m.vals

    # csc
    m = sprand(ComplexF64, 4,4, 0.5)
    sm = m |> staticize
    @test sm.colptr isa SVector
    @test sm.rowval isa SVector
    @test sm.nzval isa SVector
    @test sm.nzval == m.nzval
    @test sm.rowval == m.rowval
    @test sm.colptr == m.colptr

    dm = sm |> dynamicize
    @test dm.colptr isa Vector
    @test dm.rowval isa Vector
    @test dm.nzval isa Vector
    @test dm.nzval == m.nzval
    @test dm.rowval == m.rowval
    @test dm.colptr == m.colptr

    # diagonal
    m = Diagonal(randn(ComplexF64, 4))
    sm = m |> staticize
    @test sm.diag isa SVector
    @test sm.diag == m.diag
    dm = sm |> dynamicize
    @test dm.diag isa Vector
    @test dm.diag == m.diag

    # dense vector
    m = randn(ComplexF64, 4)
    sm = m |> staticize
    @test sm isa SVector
    @test sm == m
    dm = sm |> dynamicize
    @test dm isa Vector
    @test dm == m

    # dense matrix
    m = randn(ComplexF64, 4, 4)
    sm = m |> staticize
    @test sm isa SMatrix
    @test sm == m
    dm = sm |> dynamicize
    @test dm isa Matrix
    @test dm == m

    m = @SMatrix rand(2, 2)
    @test staticize(m) === m
    m = rand(2, 2)
    @test dynamicize(m) === m
end

@testset "SSparseMatrixCSC" begin
    m = sprand(ComplexF64, 4,4, 0.5)
    sm = m |> staticize
    @test nnz(sm) == nnz(m)
    @test nonzeros(sm) == nonzeros(m)
    @test isdense(sm) == false
    @test issparse(sm)
    @test dropzeros!(sm) == sm
    @test sm == m
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 88cfc46c-f46d-11ea-21df-59c027e4b86b
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add(["Images", "ImageIO", "ImageMagick", "ImageTransformations", "ImageFiltering"])
	using Images, LinearAlgebra, SparseArrays
	md"Libraries added"
end

# ╔═╡ 045fc73c-f5bb-11ea-210d-c3a88c7b5c08
md"# Reversing a convolution"

# ╔═╡ ce1dd8fc-f5bb-11ea-3cf4-4b20b1250bb1
md"## Functions"

# ╔═╡ f1cf67d8-f5bc-11ea-0015-192f0ebb1fdf
md"### Kernel"

# ╔═╡ f284c890-f5bb-11ea-30e8-abbc1f458cec
function Kernel_uniform(shape)
	value = 1/prod(shape)
	return fill(value, shape)
end

# ╔═╡ fcae7582-f5bb-11ea-2b79-43f8cd3178c3
function get_pad(kernel)
	return (size(kernel).-1).÷2
end

# ╔═╡ e6b903c2-f5bc-11ea-0c2c-03e099517d90
md"### Convolution and Reverse"

# ╔═╡ fb654aee-f5bc-11ea-35f4-3397194c8722
md"## Example"

# ╔═╡ 116fe702-f5bd-11ea-0ba9-c55263d026ba
md"### Importing image"

# ╔═╡ a3f814ac-f46c-11ea-33d7-81d4c435aba2
begin
	i = 2

	if i==1
		filename = "monkey"
		large_image = load("monkey.jpg")
		image = imresize(large_image; ratio=0.15)
	elseif i==2
		filename = "lena"
		large_image = load("lena.jpg")
		image = imresize(large_image; ratio=0.35)
	end

	#image = Gray.(image)
	#image = Gray.(rand(Float64, (5, 5)))
end

# ╔═╡ e731fd9c-f475-11ea-0123-b5df0503c9ba
size(image)

# ╔═╡ 5965dbc0-f4fa-11ea-3d29-536288bec4f6
typeof(image)

# ╔═╡ 1eb50ff0-f5bd-11ea-34a9-65ceeb74698a
md"### Initializing filter"

# ╔═╡ 820976ec-f476-11ea-0113-d909496fa398
begin
	#kernel = Kernel_uniform((3, 3))
	kernel = Kernel.gaussian((3, 3))

	pad = Pad(:replicate, get_pad(kernel))
	#pad = Fill(0, size(kernel))

	size(kernel)
end

# ╔═╡ 110e881e-f5bc-11ea-0a20-3dabe51098a3
function f2mat(img, f, pad=Pad(:replicate, get_pad(kernel)))
	input_shape = size(img)
	input_size = prod(input_shape)

	indices = reshape(1:input_size, input_shape)

	pad_size = get_pad(f)
	indices = padarray(indices, pad)

	indices = parent(indices)
	f = parent(f)

	output_shape = size(indices) .- size(f) .+ (1, 1)
	output_size = prod(output_shape)

	mat = sparse([], [], Float64[], input_size, output_size)

	for j=1:output_shape[2], i=1:output_shape[1]
		for n=1:size(f, 2), m=1:size(f, 1)
			input_idx = indices[i+m-1, j+n-1]
			if input_idx > 0
				output_idx = i + (j-1)*output_shape[1]
				mat[input_idx, output_idx] += f[m, n]
			end
		end
	end

	return mat'
end

# ╔═╡ 9966750a-f5bc-11ea-2402-f16ff519f56d
function convolution(image, kernel, pad=Pad(:replicate, get_pad(kernel)))
	matrix = f2mat(image, kernel, pad)
	mono = (length(size(channelview(image)))==2)

	if mono
		return reshape(matrix*vec(image), size(image))
	else
		channels = channelview(image)
		new_image = [reshape(matrix*vec(channels[i, :, :]), size(image)) for i=1:3]
		return colorview(RGB, new_image[1], new_image[2], new_image[3])
	end
end

# ╔═╡ a053dae2-f5bc-11ea-11fb-7115b2489efa
function get_original(image, kernel, pad=Pad(:replicate, get_pad(kernel)))
	matrix = f2mat(image, kernel, pad)

	mono = (length(size(channelview(image)))==2)

	if mono
		channel = channelview(image)
		estim = matrix\vec(channel)
		estim = [(x>1) ? 1 : (x<0) ? 0 : x for x in estim]	# image values in [0, 1]
		return Gray.(reshape(estim, size(image)))
	else
		channels = channelview(image)
		new_image = zeros(size(channels))
		for i=1:3
			estim = matrix\vec(channels[i, :, :])
			estim = [(x>1) ? 1 : (x<0) ? 0 : x for x in estim] # in [0, 1]
			new_image[i, :, :] = estim
		end
		return colorview(RGB, new_image)
	end
end

# ╔═╡ 3d477c82-f5bd-11ea-2571-f5cdb9fc1525
md"### Filtered image"

# ╔═╡ 9017eb36-f47d-11ea-0ae8-f3b3bb85bf53
filtered_pad = imfilter(image, kernel, pad, Algorithm.FIR())

# ╔═╡ d89f858a-f47d-11ea-20fc-ed418cc5653b
size(filtered_pad)

# ╔═╡ 3626ce80-f482-11ea-15f0-497b0ee9b3a4
#matrix = f2mat(image, kernel, "zero")

# ╔═╡ 30f93266-f4fa-11ea-35d5-b7bd56cb73e3
#convolution(image, kernel, "replicate")

# ╔═╡ 52ae79ae-f5bd-11ea-1d9b-a5963c9bcd16
md"### Reversing the convolution"

# ╔═╡ 4bd86730-f512-11ea-1c03-bfd6d008e37a
begin
	estimate = false

	if !estimate
		estimated_kernel = kernel
		estimated_pad = pad
	else
		estimated_kernel = Kernel.gaussian((3, 3))
		#estimated_kernel = Kernel_uniform((7, 7))
		#estimated_pad = Fill(0, get_pad(estimated_kernel))
		estimated_pad = Pad(:replicate, get_pad(estimated_kernel))
	end

	recovered = get_original(filtered_pad, estimated_kernel, estimated_pad)
end

# ╔═╡ 68b18282-f5bd-11ea-1801-933fb61400ab
md"#### Original | Filtered | Recovered"

# ╔═╡ 4eb55d58-f514-11ea-064c-2b8a08b3c1eb
[image filtered_pad recovered]

# ╔═╡ 18642c9c-f5bf-11ea-2d7f-0956bb770cc7
begin
	save(filename * "_original.jpg", image)
	save(filename * "_blurred.jpg", filtered_pad)
	save(filename * "_recovered.jpg", recovered)
	md"Files saved"
end

# ╔═╡ d6714f18-f60b-11ea-15b6-517d7006d9e4
md"old functions inside"
#=
begin
	function get_b(img, f, pad, mat, idx)
		pad_size = get_pad(f)

		if pad=="replicate"
			img_padded = padarray(img, Pad(:replicate,pad_size))
			img_padded[1:end-pad_size[1], 1:end-pad_size[2]] .= 0
		else
			return 0
		end

		mono = (length(size(channelview(img)))==2)

		if mono
			b = mat*vec(img)
			return reshape(b[idx], size(img))
		else
			channels = channelview(img)
			b = zeros(size(img))
			for i=1:3
				b_i = mat*vec(channels[i, :, :])
				b[i, :, :] = reshape(b_i[idx], size(img))
			end
			return colorview(RGB, b)
		end
	end

	function f2mat2(img, f, pad="zero")
		# https://stackoverflow.com/questions/16798888/2-d-convolution-as-a-matrix-matrix-multiplication

		pad_size = get_pad(f)
		mask = padarray(fill(true, size(img)), Fill(false, pad_size))

		input_shape = size(mask)
		output_shape = input_shape .- size(f) .+ (1, 1)
		input_size = prod(input_shape)
		output_size = prod(output_shape)

		mat = sparse([], [], Float64[], input_size, output_size)

		delta = input_shape[1]-size(f, 1)
		f_flat = sparse([], [], Float64[], input_shape[1], size(f, 2))

		f_flat[1:size(f, 1), :] = parent(f)
		f_flat = sparsevec(f_flat)[1:end-delta]

		delta = input_shape[1] - size(f, 1) + 1

		for i=1:size(mat, 2)÷delta, j=1:delta
			ii = (i-1)*input_shape[1]+j
			jj = (i-1)*delta+j
			mat[ii:ii+length(f_flat)-1, jj] = f_flat
		end

		mask_flat = vec(mask)
		#b = get_b(img, f, pad, mat, mask_flat)

		idx = [i for i=1:length(mask) if mask[i]]
		mat = mat[idx, :]

		return mat'
	end
end
=#

# ╔═╡ Cell order:
# ╟─045fc73c-f5bb-11ea-210d-c3a88c7b5c08
# ╟─ce1dd8fc-f5bb-11ea-3cf4-4b20b1250bb1
# ╟─88cfc46c-f46d-11ea-21df-59c027e4b86b
# ╟─f1cf67d8-f5bc-11ea-0015-192f0ebb1fdf
# ╟─f284c890-f5bb-11ea-30e8-abbc1f458cec
# ╟─fcae7582-f5bb-11ea-2b79-43f8cd3178c3
# ╟─e6b903c2-f5bc-11ea-0c2c-03e099517d90
# ╟─110e881e-f5bc-11ea-0a20-3dabe51098a3
# ╟─9966750a-f5bc-11ea-2402-f16ff519f56d
# ╟─a053dae2-f5bc-11ea-11fb-7115b2489efa
# ╟─fb654aee-f5bc-11ea-35f4-3397194c8722
# ╟─116fe702-f5bd-11ea-0ba9-c55263d026ba
# ╠═a3f814ac-f46c-11ea-33d7-81d4c435aba2
# ╟─e731fd9c-f475-11ea-0123-b5df0503c9ba
# ╟─5965dbc0-f4fa-11ea-3d29-536288bec4f6
# ╟─1eb50ff0-f5bd-11ea-34a9-65ceeb74698a
# ╠═820976ec-f476-11ea-0113-d909496fa398
# ╟─3d477c82-f5bd-11ea-2571-f5cdb9fc1525
# ╠═9017eb36-f47d-11ea-0ae8-f3b3bb85bf53
# ╟─d89f858a-f47d-11ea-20fc-ed418cc5653b
# ╠═3626ce80-f482-11ea-15f0-497b0ee9b3a4
# ╠═30f93266-f4fa-11ea-35d5-b7bd56cb73e3
# ╟─52ae79ae-f5bd-11ea-1d9b-a5963c9bcd16
# ╠═4bd86730-f512-11ea-1c03-bfd6d008e37a
# ╟─68b18282-f5bd-11ea-1801-933fb61400ab
# ╠═4eb55d58-f514-11ea-064c-2b8a08b3c1eb
# ╟─18642c9c-f5bf-11ea-2d7f-0956bb770cc7
# ╟─d6714f18-f60b-11ea-15b6-517d7006d9e4


text:
using DiscretisedFluidQueues
using LinearAlgebra, SparseArrays, StableRNGs
using Test
include("SFM_operators.jl")

T = [-2.5 2 0.5; 1 -2 1; 1 2 -3]
T_nz = T - [0.01 0 0;0 0 0;0 0 0]
T_warn = T - [0.000001 0 0;0 0 0;0 0 0]

C = [0.0; 2.0; -3.0]
m = -1 .+ 2*Int.(DiscretisedFluidQueues._strictly_pos.(C))
S = PhaseSet(C)

model = BoundedFluidQueue(T,S,12.0)

T_aug = [-2.5 0 2 0.5; 0 -2.5 2 0.5; 1 0 -2 1; 0 1 2 -3]
C_aug = [0.0;-0.0;C[2:3]]
m_aug = -1 .+ 2*Int.(DiscretisedFluidQueues._strictly_pos.(C_aug))
S_aug = PhaseSet(C_aug)
am = BoundedFluidQueue(T_aug,S_aug,12.0)

nds_vec = Array{Any,1}(undef,2)
nds_vec[1] = [0.0;3.0;4.0;12.0]
nds_vec[2] = 0.0:4.0:12.0
for nds in 1:2 # probably overkill to test everything twice but oh well!
    global nodes=nds_vec[nds]

    global nbases=3
    global dgmesh=DGMesh(nodes,nbases)

    global am=augment_model(model)

    global fv_order=3
    global fvmesh=FVMesh(nodes,fv_order)

    global order=3
    global frapmesh=FRAPMesh(nodes,order)

    @testset begin

        include("models.jl")

        include("mesh.jl")

        include("generators.jl")

    # test 5_SFM_operators

        include("me_tools.jl")
        include("polynomials.jl")

        include("distributions.jl")

        include("time_integration.jl")

        include("sims.jl")

        include("numerics.jl")

        include("limiters.jl")

    # more testing for lazy_generators now with new modularised code
    # testing for DiscretisedFluidQueue
    # etc...
    end
end


text:
using CUDA
using CUDA.CUSPARSE
using ExaPF
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using FiniteDiff

import ExaPF: ParsePSSE, PowerSystem, IndexSet

@testset "PenaltyEvaluators" begin
    @testset "Inactive constraints" begin
        datafile = joinpath(dirname(@__FILE__), "..", "data", "case9.m")
        nlp = ExaPF.ReducedSpaceEvaluator(datafile)
        u0 = ExaPF.initial(nlp)
        # Build penalty evaluator
        pen = ExaPF.PenaltyEvaluator(nlp, u0)

        u = u0
        # Update nlp to stay on manifold
        ExaPF.update!(pen, u)
        # Compute objective
        c = ExaPF.objective(pen, u)
        c_ref = ExaPF.objective(nlp, u)
        @test isa(c, Real)
        # For case9.m all constraints are inactive and penalty are equal to zero
        @test iszero(pen.infeasibility)
        @test c == c_ref

        # Compute gradient of objective
        g = similar(u)
        g_ref = similar(u)
        fill!(g, 0)
        ExaPF.gradient!(pen, g, u)
        fill!(g_ref, 0)
        ExaPF.gradient!(nlp, g_ref, u)
        @test isequal(g_ref, g)
        # Update penalty weigth
        ExaPF.update_penalty!(pen)

        # Test reset
        ExaPF.reset!(pen)
    end
    @testset "Active constraints" begin
        datafile = joinpath(dirname(@__FILE__), "..", "data", "case57.m")
        nlp = ExaPF.ReducedSpaceEvaluator(datafile)
        u0 = ExaPF.initial(nlp)

        # Build penalty evaluator
        for scaling in [true, false]
            pen = ExaPF.PenaltyEvaluator(nlp, u0; scale=scaling)
            u = nlp.u_min
            # Update nlp to stay on manifold
            ExaPF.update!(pen, u)
            # Compute objective
            c = ExaPF.objective(pen, u)
            c_ref = ExaPF.objective(nlp, u)
            @test isa(c, Real)
            # For case57.m some constraints are active, so penalty are >= 0
            @test c > c_ref

            # Update penalty weigth with a large factor to have
            # a meaningful derivative check
            ExaPF.update_penalty!(pen, η=100.)
            # Compute gradient of objective
            g = similar(u)
            fill!(g, 0)
            ExaPF.gradient!(pen, g, u)
            # Compare with finite differences
            function reduced_cost(u_)
                ExaPF.update!(pen, u_)
                return ExaPF.objective(pen, u_)
            end
            grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, u)
            @test isapprox(grad_fd, g, rtol=1e-6)
        end
    end
end


text:
# -------------------------------------------------------------------------------
# Primal problem
# -------------------------------------------------------------------------------
@with_kw_noshow struct Primal{I<:Integer, B<:BitArray{1}, V<:AbstractVector, A<:AbstractMatrix} <: AbstractData
    X::A
    y::B
    ind_pos::V = findall(y)
    ind_neg::V = findall(.~y)
    dim::I     = size(X,2)
    n::I       = length(y)
    npos::I    = length(ind_pos)
    nneg::I    = length(ind_neg)
end

Primal(X::AbstractMatrix, y::AbstractVector) = Primal(X = X, y = Bool.(y))

# -------------------------------------------------------------------------------
# Dual problem
# -------------------------------------------------------------------------------
abstract type DualType; end

@with_kw_noshow struct DTrain{I<:Integer, V<:AbstractVector, P<:AbstractVector} <: DualType
    ind_pos::V
    ind_neg::V
    inv_perm::P
    n::I    = length(inv_perm)
    npos::I = length(ind_pos)
    nneg::I = length(ind_neg)
end

DTrain(ipos, ineg, iperm) = DTrain(ind_pos = ipos, ind_neg = ineg, inv_perm = iperm)

@with_kw_noshow struct DValidation{I<:Integer, V<:AbstractVector, P<:AbstractVector} <: DualType
    ind_pos::V
    ind_neg::V
    inv_perm::P
    n::I    = length(inv_perm)
    npos::I = length(ind_pos)
    nneg::I = length(ind_neg)
end

DValidation(ipos, ineg, iperm) = DValidation(ind_pos = ipos, ind_neg = ineg, inv_perm = iperm)

struct DTest{I<:Integer} <: DualType
    n::I
end

# type Dual
@with_kw_noshow struct Dual{T<:DualType, I<:Integer, V<:AbstractVector, A<:AbstractMatrix} <: AbstractData
    type::T
    io::IO

    K::A
    nα::I
    nβ::I
    n::I

    ind_α::V = 1:nα
    ind_β::V = nα .+ (1:nβ)
end

Dual(type::DualType, io::IO, K::AbstractMatrix, nα::Int, nβ::Int, n::Int) =
    Dual(type = type, io = io, K = K, nα = nα, nβ = nβ, n = n)

# train data
function Dual(
    M::AbstractModel,
    Xtrain::AbstractMatrix,
    ytrain::BitArray{1};
    kwargs...
)

    K, n, nα, nβ, ind_pos, ind_neg, inv_perm = kernelmatrix(M, Xtrain, ytrain; kwargs...)

    type = DTrain(ind_pos, ind_neg, inv_perm)
    io   = IOBuffer()
    close(io)

    return Dual(type, io, K, nα, nβ, n)
end

# validation data
function Dual(
    M::AbstractModel,
    Xtrain::AbstractMatrix,
    ytrain::BitArray{1},
    Xvalid::AbstractMatrix,
    yvalid::BitArray{1};
    kwargs...
)

    K, n, nα, nβ, ind_pos, ind_neg, inv_perm = kernelmatrix(M, Xtrain, ytrain, Xvalid, yvalid; kwargs...)

    type = DValidation(ind_pos, ind_neg, inv_perm)
    io   = IOBuffer()
    close(io)

    return Dual(type, io, K, nα, nβ, n)
end

# test data
function Dual(
    M::AbstractModel,
    Xtrain::AbstractMatrix,
    ytrain::BitArray{1},
    Xtest::AbstractMatrix;
    kwargs...
)

    K, n, nα, nβ = kernelmatrix(M, Xtrain, ytrain, Xtest; kwargs...)

    type = DTest(n)
    io   = IOBuffer()
    close(io)

    return Dual(type, io, K, nα, nβ, n)
end

# Dual dataset load function
function Dual(file::AbstractString; kwargs...)

    type_int, io, out = load_kernelmatrix(file; kwargs...)

    if type_int == 0
        K, n, nα, nβ, ind_pos, ind_neg, inv_perm = out
        type = DTrain(ind_pos, ind_neg, inv_perm)

    elseif type_int == 1
        K, n, nα, nβ, ind_pos, ind_neg, inv_perm = out
        type = DValidation(ind_pos, ind_neg, inv_perm)

    elseif type_int == 2
        K, n, nα, nβ = out
        type = DTest(n)
    end

    return Dual(type, io, K, nα, nβ, n)
end


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

using SeparableOptimization
using PiecewiseQuadratics
using LinearAlgebra
using SparseArrays
using Documenter
using Test

@testset "Doctests" begin
    DocMeta.setdocmeta!(SeparableOptimization, :DocTestSetup,
                        :(using SeparableOptimization, PiecewiseQuadratics); recursive=true)
    doctest(SeparableOptimization)
end

@testset "ADMM Structs" begin
    include("test_admm_structs.jl")
end

@testset "ADMM" begin
    include("test_admm.jl")
end

@testset "Examples" begin
    include("test_examples.jl")
end


text:
using Test, PolynomialBases
import SymEngine

tol = 5.e-14

for p in 0:6
    basis_symbo = LobattoLegendre(p, SymEngine.Basic)
    basis_float = LobattoLegendre(p, Float64)
    @test maximum(abs.( float.(basis_symbo.nodes) - basis_float.nodes )) < tol
    @test maximum(abs.( float.(basis_symbo.weights) - basis_float.weights )) < tol
    @test maximum(abs.( float.(basis_symbo.baryweights) - basis_float.baryweights )) < tol
    @test maximum(abs.( float.(basis_symbo.D) - basis_float.D, )) < tol
end
@test_throws ArgumentError LobattoLegendre(7, SymEngine.Basic)

for p in 0:4
    basis_symbo = GaussLegendre(p, SymEngine.Basic)
    basis_float = GaussLegendre(p, Float64)
    @test maximum(abs.( float.(basis_symbo.nodes) - basis_float.nodes )) < tol
    @test maximum(abs.( float.(basis_symbo.weights) - basis_float.weights )) < tol
    @test maximum(abs.( float.(basis_symbo.baryweights) - basis_float.baryweights )) < tol
    @test maximum(abs.( float.(basis_symbo.D) - basis_float.D, )) < tol
end
@test_throws ArgumentError GaussLegendre(5, SymEngine.Basic)

interpolation_matrix([-1, 1], LobattoLegendre(4, SymEngine.Basic))
interpolation_matrix([-1, 1], GaussLegendre(4, SymEngine.Basic))


text:
module Cache

using ...SimpleTasks.Types

import Base

export Service, exists, put!, get, delete!

"""
    haskey(cache::CacheService, key::AbstractString)

Returns true if the key is found in the cache.
"""
function Base.haskey(cache::CacheService, key::AbstractString)
    error("exists is unimplemented for $cache")
end

"""
    put!(cache::CacheService, key::AbstractString)

Reads the value from IO to store into the cache under the key.
"""
function put!(cache::CacheService, key::AbstractString, value_io::IO)
    error("put is unimplemented for $cache")
end
function Base.setindex!(cache::CacheService, value_io::IO, key::AbstractString)
    put!(cache, key, value_io)
end

"""
    get(cache::CacheService, key::AbstractString)

Get an IO stream to the object that is cached with this key.
Throws KeyError if key doesn't exist.
"""
function get(cache::CacheService, key::AbstractString)
    error("get is unimplemented for $cache")
end
function Base.getindex(cache::CacheService, key::AbstractString)
    get(cache, key)
end

"""
    delete!(cache::CacheService, key::AbstractString)

Delete the cached value for the given key.
Returns cache and does nothing if key is not found (this is what Julia
`Dict` does)
"""
function Base.delete!(cache::CacheService, key::AbstractString)
    error("delete! is unimplemented for $cache")
end

"""
    clear!(cache::CacheService)

Delete all cached items.
"""
function clear!(cache::CacheService)
    error("clear! is unimplemented for $cache")
end

end #module Cache

max_forks_repo_name: seung-lab/SimpleTasks


text:
# commands
const CONNECT = 0x10
const CONNACK = 0x20
const PUBLISH = 0x30
const PUBACK = 0x40
const PUBREC = 0x50
const PUBREL = 0x60
const PUBCOMP = 0x70
const SUBSCRIBE = 0x80
const SUBACK = 0x90
const UNSUBSCRIBE = 0xA0
const UNSUBACK = 0xB0
const PINGREQ = 0xC0
const PINGRESP = 0xD0
const DISCONNECT = 0xE0
const QOS_1_BYTE = 0x02
const QOS_2_BYTE = 0x04

mutable struct TestFileHandler <: IO
    out_channel::Channel{UInt8}
    in_channel::Channel{UInt8}
    closenotify::Condition
    new_input::Condition
    closed::Bool
end

function TestFileHandler()
    th = TestFileHandler(Channel{UInt8}(256), Channel{UInt8}(256), Condition(), Condition(), false)
    put_from_file(th, "data/input/connack.dat")
    return th
end

function read(fh::TestFileHandler, t)
    return take!(fh.in_channel)
end

function read(fh::TestFileHandler, t::Type{UInt8})
    return take!(fh.in_channel)
end

function read(fh::TestFileHandler, length::Integer)
    data = Vector{UInt8}()
    for i = 1:length
        append!(data, take!(fh.in_channel))
    end
    return data
end

function close(fh::TestFileHandler)
    fh.closed = true
    notify(fh.closenotify)
end

function write(fh::TestFileHandler, data::UInt8)
    put!(fh.out_channel, data)
end

function put(fh::TestFileHandler, data::Array{UInt8})
  for i in data
      put!(fh.in_channel, i)
  end
end

function put(fh::TestFileHandler, data::SubArray)
  for i in data
      put!(fh.in_channel, i)
  end
end

function put(fh::TestFileHandler, data::UInt16)
    buffer = PipeBuffer()
    write(buffer, data)
    put(fh, reverse(read(buffer), 1)) #TODO use take! instead?
end

function put(fh::TestFileHandler, data::UInt8)
  put!(fh.in_channel, data)
end

function put_from_file(fh::TestFileHandler, filename)
  put(fh, read_all_to_arr(filename))
end

function get_mid_index(data::Array{UInt8})
  cmd = data[1] & 0xF0
  if cmd == CONNACK || cmd == PINGRESP || cmd == CONNECT || cmd == DISCONNECT || cmd == PINGREQ
    return -1
  elseif cmd == PUBLISH
    qos = data[1] & 0x06
    if qos == QOS_1_BYTE || qos == QOS_2_BYTE
      buffer = PipeBuffer()
      write(buffer, data[4])
      write(buffer, data[3])
      topic_len = read_len(buffer)
      return 5 + topic_len
    else
      # QOS_0 has no message id
      return -1
    end
  else
    # all other packets have m_id as their 3rd and 4th byte
    return 3
  end
end

function put_from_file(fh::TestFileHandler, filename, messageId::UInt16)
  data = read_all_to_arr(filename)
  mid_index = get_mid_index(data)
  if mid_index > 0
    put(fh, view(data, 1:mid_index - 1))
    put(fh, messageId)
    put(fh, view(data, mid_index + 2:length(data)))
  else
    put(fh, data)
  end
end

function read_all_to_arr(filename)
    file = open(filename, "r")
    data = Vector{UInt8}()
    while !eof(file)
        append!(data, read(file, UInt8))
    end
    close(file)
    return data
end

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT