text:
#####
##### `LearnLogger`
#####

"""
    LearnLogger

A struct that wraps a `TensorBoardLogger.TBLogger` in order to enforce the following:

- all values logged to Tensorboard should be accessible to the `post_epoch_callback`
  argument to [`learn!`](@ref)
- all values that are cached during [`learn!`](@ref) should be logged to Tensorboard

To access values logged to a `LearnLogger` instance, inspect the instance's `logged` field.
"""
struct LearnLogger
    path::String
    tensorboard_logger::TensorBoardLogger.TBLogger
    logged::Dict{String,Vector{Any}}
end

function LearnLogger(path, run_name; kwargs...)
    tensorboard_logger = TBLogger(joinpath(path, run_name); kwargs...)
    return LearnLogger(path, tensorboard_logger, Dict{String,Any}())
end

function log_event!(logger::LearnLogger, value)
    logged = string(now(), " | ", value)
    TensorBoardLogger.log_text(logger.tensorboard_logger, "events", logged)
    return logged
end

function log_plot!(logger::LearnLogger, field::AbstractString, plot, plot_data)
    values = get!(() -> Any[], logger.logged, field)
    push!(values, plot_data)
    TensorBoardLogger.log_image(logger.tensorboard_logger, field, plot; step=length(values))
    return plot
end

function log_value!(logger::LearnLogger, field::AbstractString, value)
    values = get!(() -> Any[], logger.logged, field)
    push!(values, value)
    TensorBoardLogger.log_value(logger.tensorboard_logger, field, value;
                                step=length(values))
    return value
end

function log_resource_info!(logger, section::AbstractString, info::ResourceInfo;
                            suffix::AbstractString="")
    log_value!(logger, section * "/time_in_seconds" * suffix, info.time_in_seconds)
    log_value!(logger, section * "/gc_time_in_seconds" * suffix, info.gc_time_in_seconds)
    log_value!(logger, section * "/allocations" * suffix, info.allocations)
    log_value!(logger, section * "/memory_in_mb" * suffix, info.memory_in_mb)
    return info
end

function log_resource_info!(f, logger, section::AbstractString; suffix::AbstractString="")
    result, resource_info = call_with_resource_info(f)
    log_resource_info!(logger, section, resource_info; suffix=suffix)
    return result
end

"""
    Base.flush(logger)

Persist possibly transient logger state.
"""
Base.flush(logger::LearnLogger) = nothing

"""
    forwarding_task = forward_logs(channel, logger::LearnLogger)

Forwards logs with values supported by `TensorBoardLogger` to `logger::LearnLogger`:
- string events of type `AbstractString`
- scalars of type `Union{Real,Complex}`
- plots that `TensorBoardLogger` can convert to raster images

returns the `forwarding_task:::Task` that does the forwarding.
To cleanly stop forwarding, `close(channel)` and `wait(forwarding_task)`.

outbox is a Channel or RemoteChannel of Pair{String, Any}
field names starting with "__plot__" forward to TensorBoardLogger.log_image
"""
function forward_logs(outbox, logger::LearnLogger)
    @async try
        while true
            (field, value) = take!(outbox)
            if typeof(value) <: AbstractString
                log_event!(logger, value)
            elseif startswith(field, "__plot__")
                original_field = field[9:end]
                values = get!(() -> Any[], logger.logged, original_field)
                TensorBoardLogger.log_image(logger.tensorboard_logger, original_field,
                                            value; step=length(values))
            elseif typeof(value) <: Union{Real,Complex}
                log_value!(logger, field, value)
            end
        end
    catch e
        if !(isa(e, InvalidStateException) && e.state == :closed)
            @error "error forwarding logs, STOPPING FORWARDING!" exception = (e,
                                                                              catch_backtrace())
        end
    end
end

#####
##### `predict!!`
#####

"""
    predict!(model::AbstractClassifier,
             predicted_soft_labels::AbstractMatrix,
             batches, logger::LearnLogger;
             logger_prefix::AbstractString)

Return `mean_loss` of all `batches` after using `model` to predict their soft labels
and storing those results in `predicted_soft_labels`.

The following quantities are logged to `logger`:
 - `<logger_prefix>/loss_per_batch`
 - `<logger_prefix>/mean_loss_per_epoch`
 - `<logger_prefix>/\$resource_per_batch`

Where...

- `model` is a model that outputs soft labels when called on a batch of `batches`,
  `model(batch)`.

- `predicted_soft_labels` is a matrix whose columns correspond to classes and
  whose rows correspond to samples in batches, and which is filled in with soft-label
  predictions.

- `batches` is an iterable of batches, where each element of
  the iterable takes the form `(batch, votes_locations)`. Internally, `batch` is
  passed to [`loss_and_prediction`](@ref) as `loss_and_prediction(model, batch...)`.

 """
function predict!(model::AbstractClassifier, predicted_soft_labels::AbstractMatrix, batches,
                  logger; logger_prefix::AbstractString)
    losses = Float32[]
    for (batch, votes_locations) in batches
        batch_loss = log_resource_info!(logger, logger_prefix; suffix="_per_batch") do
            batch_loss, soft_label_batch = loss_and_prediction(model, batch...)
            for (i, soft_label) in enumerate(eachcol(soft_label_batch))
                predicted_soft_labels[votes_locations[i], :] = soft_label
            end
            return batch_loss
        end
        log_value!(logger, logger_prefix * "/loss_per_batch", batch_loss)
        push!(losses, batch_loss)
    end
    mean_loss = mean(losses)
    log_value!(logger, logger_prefix * "/mean_loss_per_epoch", mean_loss)
    return mean_loss
end

#####
##### `evaluate!`
#####

"""
    evaluate!(predicted_hard_labels::AbstractVector,
              predicted_soft_labels::AbstractMatrix,
              elected_hard_labels::AbstractVector,
              classes, logger::LearnLogger;
              logger_prefix, logger_suffix,
              votes::Union{Nothing,AbstractMatrix}=nothing,
              thresholds=0.0:0.01:1.0,
              optimal_threshold_class::Union{Nothing,Integer}=nothing)

Return `nothing` after computing and logging a battery of classifier performance
metrics that each compare `predicted_soft_labels` and/or `predicted_hard_labels`
agaist `elected_hard_labels`.

The following quantities are logged to `logger`:
    - `<logger_prefix>/metrics<logger_suffix>`
    - `<logger_prefix>/\$resource<logger_suffix>`

Where...

- `predicted_soft_labels` is a matrix of soft labels whose columns correspond to
  classes and whose rows correspond to samples in the evaluation set.

- `predicted_hard_labels` is a vector of hard labels where the `i`th element
  is the hard label predicted by the model for sample `i` in the evaulation set.

- `elected_hard_labels` is a vector of hard labels where the `i`th element
  is the hard label elected as "ground truth" for sample `i` in the evaulation set.

- `thresholds` are the range of thresholds used by metrics (e.g. PR curves) that
  are calculated on the `predicted_soft_labels` for a range of thresholds.

- `votes` is a matrix of hard labels whose columns correspond to voters and whose
  rows correspond to the samples in the test set that have been voted on. If
  `votes[sample, voter]` is not a valid hard label for `model`, then `voter` will
  simply be considered to have not assigned a hard label to `sample`.

- `optimal_threshold_class` is the class index (`1` or `2`) for which to calculate
  an optimal threshold for converting the `predicted_soft_labels` to
  `predicted_hard_labels`. If present, the input `predicted_hard_labels` will be
  ignored and new `predicted_hard_labels` will be recalculated from the new threshold.
  This is only a valid parameter when `length(classes) == 2`
"""
function evaluate!(predicted_hard_labels::AbstractVector,
                   predicted_soft_labels::AbstractMatrix,
                   elected_hard_labels::AbstractVector, classes, logger;
                   logger_prefix::AbstractString, logger_suffix::AbstractString="",
                   votes::Union{Nothing,AbstractMatrix}=nothing, thresholds=0.0:0.01:1.0,
                   optimal_threshold_class::Union{Nothing,Integer}=nothing)
    _validate_threshold_class(optimal_threshold_class, classes)

    log_resource_info!(logger, logger_prefix; suffix=logger_suffix) do
        plot, plot_data = evaluation_metrics_plot(predicted_hard_labels,
                                                  predicted_soft_labels,
                                                  elected_hard_labels, classes, thresholds;
                                                  votes=votes,
                                                  optimal_threshold_class=optimal_threshold_class)
        log_plot!(logger, logger_prefix * "/metrics" * logger_suffix, plot, plot_data)
        if haskey(plot_data, "spearman_correlation")
            log_value!(logger, logger_prefix * "/spearman_correlation" * logger_suffix,
                       plot_data["spearman_correlation"].ρ)
        end
    end
    return nothing
end

function _calculate_stratified_ea_kappas(predicted_hard_labels, elected_hard_labels,
                                         class_count, strata)
    groups = reduce(∪, strata)
    kappas = Pair{String,Any}[]
    for group in groups
        index = group .∈ strata
        predicted = predicted_hard_labels[index]
        elected = elected_hard_labels[index]
        k = _calculate_ea_kappas(predicted, elected, class_count)
        push!(kappas,
              group => (per_class=k.per_class, multiclass=k.multiclass, n=sum(index)))
    end
    return sort(kappas; by=p -> last(p).multiclass)
end

"""
    _calculate_ea_kappas(predicted_hard_labels, elected_hard_labels, classes)

Return `NamedTuple` with keys `:per_class`, `:multiclass` containing the Cohen's
Kappa per-class and over all classes, respectively. The value of output key
`:per_class` is an `Array` such that item `i` is the Cohen's kappa calculated
for class `i`.

Where...

- `predicted_hard_labels` is a vector of hard labels where the `i`th element
  is the hard label predicted by the model for sample `i` in the evaulation set.

- `elected_hard_labels` is a vector of hard labels where the `i`th element
  is the hard label elected as "ground truth" for sample `i` in the evaulation set.

- `class_count` is the number of possible classes.

"""
function _calculate_ea_kappas(predicted_hard_labels, elected_hard_labels, class_count)
    multiclass = first(cohens_kappa(class_count,
                                    zip(predicted_hard_labels, elected_hard_labels)))

    CLASS_VS_ALL_CLASS_COUNT = 2
    per_class = map(1:class_count) do class_index
        predicted = ((label == class_index) + 1 for label in predicted_hard_labels)
        elected = ((label == class_index) + 1 for label in elected_hard_labels)
        return first(cohens_kappa(CLASS_VS_ALL_CLASS_COUNT, zip(predicted, elected)))
    end
    return (per_class=per_class, multiclass=multiclass)
end

"""
    _calculate_ira_kappas(votes, classes)

Return `NamedTuple` with keys `:per_class`, `:multiclass` containing the Cohen's
Kappa for inter-rater agreement (IRA) per-class and over all classes, respectively.
The value of output key `:per_class` is an `Array` such that item `i` is the
IRA kappa calculated for class `i`.

Where...

- `votes` is a matrix of hard labels whose columns correspond to voters and whose
  rows correspond to the samples in the test set that have been voted on. If
  `votes[sample, voter]` is not a valid hard label for `model`, then `voter` will
  simply be considered to have not assigned a hard label to `sample`.

- `classes` all possible classes voted on.

Returns `nothing` if `votes` has only a single voter (i.e., a single column) or if
no two voters rated the same sample. Note that vote entries of `0` are taken to
mean that the voter did not rate that sample.
"""
function _calculate_ira_kappas(votes, classes)
    (isnothing(votes) || size(votes, 2) < 2) && return nothing  # no votes given or only one expert

    all_hard_label_pairs = Array{Int}(undef, 0, 2)
    num_voters = size(votes, 2)
    for i_voter in 1:(num_voters - 1)
        for j_voter in (i_voter + 1):num_voters
            all_hard_label_pairs = vcat(all_hard_label_pairs, votes[:, [i_voter, j_voter]])
        end
    end
    hard_label_pairs = filter(row -> all(row .!= 0), collect(eachrow(all_hard_label_pairs)))
    length(hard_label_pairs) > 0 || return nothing  # No common observations voted on
    length(hard_label_pairs) < 10 &&
        @warn "...only $(length(hard_label_pairs)) in common, potentially questionable IRA results"

    multiclass_ira = first(cohens_kappa(length(classes), hard_label_pairs))

    CLASS_VS_ALL_CLASS_COUNT = 2
    per_class_ira = map(1:length(classes)) do class_index
        class_v_other_hard_label_pair = map(row -> 1 .+ (row .== class_index),
                                            hard_label_pairs)
        return first(cohens_kappa(CLASS_VS_ALL_CLASS_COUNT, class_v_other_hard_label_pair))
    end
    return (per_class=per_class_ira, multiclass=multiclass_ira)
end

function _spearman_corr(predicted_soft_labels, elected_soft_labels)
    n = length(predicted_soft_labels)
    ρ = StatsBase.corspearman(predicted_soft_labels, elected_soft_labels)
    if isnan(ρ)
        @warn "Uh oh, correlation is NaN! Probably because StatsBase.corspearman(...)
               returns NaN when a set of labels is all the same!"
        # Note: accounted for in https://github.com/JuliaStats/HypothesisTests.jl/pull/53/files;
        # probably not worth implementing here until we need it (at which point maybe
        # it will be ready in HypothesisTests!)
    end

    # 95% confidence interval calculated according to
    # https://stats.stackexchange.com/questions/18887/how-to-calculate-a-confidence-interval-for-spearmans-rank-correlation
    stderr = 1.0 / sqrt(n - 3)
    delta = 1.96 * stderr
    ci_lower = tanh(atanh(ρ) - delta)
    ci_upper = tanh(atanh(ρ) + delta)
    return (ρ=ρ, n=n, ci_lower=round(ci_lower; digits=3),
            ci_upper=round(ci_upper; digits=3))
end

"""
    _calculate_spearman_correlation(predicted_soft_labels, votes, classes)

Return `NamedTuple` with keys `:ρ`, `:n`, `:ci_lower`, and `ci_upper` that are
the Spearman correlation constant ρ and its 95% confidence interval bounds.
Only valid for binary classification problems (i.e., `length(classes) == 2`)

Where...

- `predicted_soft_labels` is a matrix of soft labels whose columns correspond to
  the two classes and whose rows correspond to the samples in the test set that have been
  classified. For a given sample, the two class column values must sum to 1 (i.e.,
  softmax has been applied to the classification output).

- `votes` is a matrix of hard labels whose columns correspond to voters and whose
  rows correspond to the samples in the test set that have been voted on. If
  `votes[sample, voter]` is not a valid hard label for `model`, then `voter` will
  simply be considered to have not assigned a hard label to `sample`. May contain
  a single voter (i.e., a single column).

- `classes` are the two classes voted on.
"""
function _calculate_spearman_correlation(predicted_soft_labels, votes, classes)
    length(classes) > 2 && throw(ArgumentError("Only valid for 2-class problems"))
    if !all(x -> x ≈ 1, sum(predicted_soft_labels; dims=2))
        @info predicted_soft_labels
        throw(ArgumentError("Input probabiliities fail softmax assumption"))
    end

    class_index = 1 # Note: Result will be the same whether class 1 or class 2
    elected_soft_labels = Vector{Float64}()
    for sample_votes in eachrow(votes)
        actual_sample_votes = filter(v -> v in Set([1, 2]), sample_votes)
        push!(elected_soft_labels, mean(actual_sample_votes .== class_index))
    end
    return _spearman_corr(predicted_soft_labels[:, class_index], elected_soft_labels)
end

function _calculate_optimal_threshold_from_discrimination_calibration(predicted_soft_labels,
                                                                      votes; thresholds,
                                                                      class_of_interest_index)
    elected_probabilities = _elected_probabilities(votes, class_of_interest_index)
    bin_count = min(size(votes, 2) + 1, 10)
    per_threshold_curves = map(thresholds) do thresh
        return calibration_curve(elected_probabilities,
                                 predicted_soft_labels[:, class_of_interest_index] .>=
                                 thresh; bin_count=bin_count)
    end
    i_min = argmin([c.mean_squared_error for c in per_threshold_curves])
    curve = per_threshold_curves[i_min]
    return (threshold=collect(thresholds)[i_min], mse=curve.mean_squared_error,
            plot_curve_data=(mean.(curve.bins), curve.fractions))
end

function _elected_probabilities(votes, class_of_interest_index)
    elected_probabilities = Vector{Float64}()
    for sample_votes in eachrow(votes)
        actual_sample_votes = filter(v -> v in Set([1, 2]), sample_votes)
        push!(elected_probabilities, mean(actual_sample_votes .== class_of_interest_index))
    end
    return elected_probabilities
end

function _calculate_voter_discrimination_calibration(votes; class_of_interest_index)
    elected_probabilities = _elected_probabilities(votes, class_of_interest_index)

    bin_count = min(size(votes, 2) + 1, 10)
    per_voter_calibration_curves = map(1:size(votes, 2)) do i_voter
        return calibration_curve(elected_probabilities,
                                 votes[:, i_voter] .== class_of_interest_index;
                                 bin_count=bin_count)
    end

    return (mse=map(curve -> curve.mean_squared_error, per_voter_calibration_curves),
            plot_curve_data=map(curve -> (mean.(curve.bins), curve.fractions),
                                per_voter_calibration_curves))
end

function _get_optimal_threshold_from_ROC(per_class_roc_curves; thresholds,
                                         class_of_interest_index)
    dist = (p1, p2) -> sqrt((p1[1] - p2[1])^2 + (p1[2] - p2[2])^2)
    min = Inf
    curr_counter = 1
    opt_point = nothing
    threshold_idx = 1
    for point in zip(per_class_roc_curves[class_of_interest_index][1],
            per_class_roc_curves[class_of_interest_index][2])
        d = dist((0, 1), point)
        if d < min
            min = d
            threshold_idx = curr_counter
            opt_point = point
        end
        curr_counter += 1
    end
    return collect(thresholds)[threshold_idx]
end

function _validate_threshold_class(optimal_threshold_class, classes)
    isnothing(optimal_threshold_class) && return nothing
    length(classes) == 2 ||
        throw(ArgumentError("Only valid for binary classification problems"))
    optimal_threshold_class in Set([1, 2]) ||
        throw(ArgumentError("Invalid threshold class"))
    return nothing
end

"""
    evaluation_metrics_plot(predicted_hard_labels::AbstractVector,
                            predicted_soft_labels::AbstractMatrix,
                            elected_hard_labels::AbstractVector,
                            classes,
                            thresholds=0.0:0.01:1.0;
                            votes::Union{Nothing,AbstractMatrix}=nothing,
                            strata::Union{Nothing,AbstractVector{Set{T}} where T}=nothing,
                            optimal_threshold_class::Union{Nothing,Integer}=nothing)

Return a plot and dictionary containing a battery of classifier performance
metrics that each compare `predicted_soft_labels` and/or `predicted_hard_labels`
agaist `elected_hard_labels`.

Where...

- `predicted_soft_labels` is a matrix of soft labels whose columns correspond to
  classes and whose rows correspond to samples in the evaluation set.

- `predicted_hard_labels` is a vector of hard labels where the `i`th element
  is the hard label predicted by the model for sample `i` in the evaulation set.

- `elected_hard_labels` is a vector of hard labels where the `i`th element
  is the hard label elected as "ground truth" for sample `i` in the evaulation set.

- `thresholds` are the range of thresholds used by metrics (e.g. PR curves) that
  are calculated on the `predicted_soft_labels` for a range of thresholds.

- `votes` is a matrix of hard labels whose columns correspond to voters and whose
  rows correspond to the samples in the test set that have been voted on. If
  `votes[sample, voter]` is not a valid hard label for `model`, then `voter` will
  simply be considered to have not assigned a hard label to `sample`.

- `strata` is a vector of sets of (arbitrarily typed) groups/strata for each sample
  in the evaluation set, or `nothing`. If not `nothing`, per-class and multiclass
  kappas will also be calculated per group/stratum.

- `optimal_threshold_class` is the class index (`1` or `2`) for which to calculate
  an optimal threshold for converting the `predicted_soft_labels` to
  `predicted_hard_labels`. If present, the input `predicted_hard_labels` will be
  ignored and new `predicted_hard_labels` will be recalculated from the new threshold.
  This is only a valid parameter when `length(classes) == 2`
"""
function evaluation_metrics_plot(predicted_hard_labels::AbstractVector,
                                 predicted_soft_labels::AbstractMatrix,
                                 elected_hard_labels::AbstractVector, classes, thresholds;
                                 votes::Union{Nothing,AbstractMatrix}=nothing,
                                 strata::Union{Nothing,AbstractVector{Set{T}} where T}=nothing,
                                 optimal_threshold_class::Union{Nothing,Integer}=nothing)
    _validate_threshold_class(optimal_threshold_class, classes)

    class_count = length(classes)
    class_vector = collect(classes) # Plots.jl expects this to be an `AbstractVector`
    class_labels = string.(class_vector)
    per_class_stats = per_class_confusion_statistics(predicted_soft_labels,
                                                     elected_hard_labels, thresholds)
    plot_dict = Dict()
    plot_dict["class_labels"] = class_labels
    plot_dict["thresholds"] = thresholds

    # ROC curves
    plot_dict["per_class_roc_curves"] = [(map(t -> t.false_positive_rate, stats),
                                          map(t -> t.true_positive_rate, stats))
                                         for stats in per_class_stats]
    plot_dict["per_class_roc_aucs"] = [area_under_curve(x, y)
                                       for (x, y) in plot_dict["per_class_roc_curves"]]

    # Optionally calculate optimal threshold
    if !isnothing(optimal_threshold_class)
        plot_dict["optimal_threshold_class"] = optimal_threshold_class
        threshold = nothing

        # If votes exist, calculate the threshold based on comparing against
        # vote probabilities. Otherwise, use the ROC curve.
        if !isnothing(votes)
            c = _calculate_optimal_threshold_from_discrimination_calibration(predicted_soft_labels,
                                                                             votes;
                                                                             thresholds=thresholds,
                                                                             class_of_interest_index=optimal_threshold_class)
            threshold = c.threshold
            plot_dict["discrimination_calibration_curve"] = c.plot_curve_data
            plot_dict["discrimination_calibration_score"] = c.mse

            expert_cal = _calculate_voter_discrimination_calibration(votes;
                                                                     class_of_interest_index=optimal_threshold_class)
            plot_dict["per_expert_discrimination_calibration_curves"] = expert_cal.plot_curve_data
            plot_dict["per_expert_discrimination_calibration_scores"] = expert_cal.mse
        else
            # ...based on ROC curve otherwise
            threshold = _get_optimal_threshold_from_ROC(plot_dict["per_class_roc_curves"];
                                                        thresholds=thresholds,
                                                        class_of_interest_index=optimal_threshold_class)
        end
        plot_dict["optimal_threshold"] = threshold

        # Recalculate `predicted_hard_labels` with this new threshold
        other_class = optimal_threshold_class == 1 ? 2 : 1
        for (i, row) in enumerate(eachrow(predicted_soft_labels))
            predicted_hard_labels[i] = row[optimal_threshold_class] .>= threshold ?
                                       optimal_threshold_class : other_class
        end
    end

    # PR curves
    plot_dict["per_class_pr_curves"] = [(map(t -> t.true_positive_rate, stats),
                                         map(t -> t.precision, stats))
                                        for stats in per_class_stats]

    # Cohen's kappa
    kappas = _calculate_ea_kappas(predicted_hard_labels, elected_hard_labels, class_count)
    plot_dict["per_class_kappas"] = kappas.per_class
    plot_dict["multiclass_kappa"] = kappas.multiclass
    ira = _calculate_ira_kappas(votes, classes)
    if !isnothing(ira)
        plot_dict["per_class_IRA_kappas"] = ira.per_class
        plot_dict["multiclass_IRA_kappas"] = ira.multiclass
    end

    # Stratified kappas
    if !isnothing(strata)
        plot_dict["stratified_kappas"] = _calculate_stratified_ea_kappas(predicted_hard_labels,
                                                                         elected_hard_labels,
                                                                         class_count,
                                                                         strata)
    end

    # Reliability calibration curves
    per_class_reliability_calibration_curves = map(1:class_count) do class_index
        class_probabilities = view(predicted_soft_labels, :, class_index)
        return calibration_curve(class_probabilities, elected_hard_labels .== class_index)
    end
    plot_dict["per_class_reliability_calibration_curves"] = map(x -> (mean.(x.bins),
                                                                      x.fractions),
                                                                per_class_reliability_calibration_curves)
    plot_dict["per_class_reliability_calibration_scores"] = map(x -> x.mean_squared_error,
                                                                per_class_reliability_calibration_curves)

    # Confusion matrix
    plot_dict["confusion_matrix"] = confusion_matrix(class_count,
                                                     zip(predicted_hard_labels,
                                                         elected_hard_labels))

    # Log Spearman correlation, iff this is a binary classification problem
    if length(classes) == 2 && !isnothing(votes)
        plot_dict["spearman_correlation"] = _calculate_spearman_correlation(predicted_soft_labels,
                                                                            votes, classes)
    end
    return evaluation_metrics_plot(plot_dict), plot_dict
end

function per_class_confusion_statistics(predicted_soft_labels::AbstractMatrix,
                                        elected_hard_labels::AbstractVector, thresholds)
    class_count = size(predicted_soft_labels, 2)
    confusions = [[confusion_matrix(2) for _ in 1:length(thresholds)]
                  for _ in 1:class_count]
    for class_index in 1:class_count
        for label_index in 1:length(elected_hard_labels)
            predicted_soft_label = predicted_soft_labels[label_index, class_index]
            elected = (elected_hard_labels[label_index] == class_index) + 1
            for (threshold_index, threshold) in enumerate(thresholds)
                predicted = (predicted_soft_label >= threshold) + 1
                confusions[class_index][threshold_index][predicted, elected] += 1
            end
        end
    end
    return [binary_statistics.(confusions[i], 2) for i in 1:class_count]
end

#####
##### `learn!`
#####

"""
    learn!(model::AbstractClassifier, logger,
           get_train_batches, get_test_batches, votes,
           elected=majority.(eachrow(votes), (1:length(classes(model)),));
           epoch_limit=100, post_epoch_callback=(_ -> nothing),
           optimal_threshold_class::Union{Nothing,Integer}=nothing,
           test_set_logger_prefix="test_set")

Return `model` after optimizing its parameters across multiple epochs of
training and test, logging Lighthouse's standardized suite of classifier
performance metrics to `logger` throughout the optimization process.

The following phases are executed at each epoch (note: in the below lists
of logged values, `\$resource` takes the values of the field names of
`Lighthouse.ResourceInfo`):

1. Train `model` by calling `train!(model, get_train_batches(), logger)`.
   The following quantities are logged to `logger` during this phase:
    - `train/loss_per_batch`
    - any additional quantities logged by the relevant model/framework-specific
      implementation of `train!`.

2. Compute `model`'s predictions on test set provided by `get_test_batches()`
   (see below for details). The following quantities are logged to `logger`
   during this phase:
    - `<test_set_logger_prefix>_prediction/loss_per_batch`
    - `<test_set_logger_prefix>_prediction/mean_loss_per_epoch`
    - `<test_set_logger_prefix>_prediction/\$resource_per_batch`

3. Compute a battery of metrics to evaluate `model`'s performance on the test
   set based on the test set prediction phase. The following quantities are
   logged to `logger` during this phase:
    - `<test_set_logger_prefix>_evaluation/metrics_per_epoch`
    - `<test_set_logger_prefix>_evaluation/\$resource_per_epoch`

4. Call `post_epoch_callback(current_epoch)`.

Where...

- `get_train_batches` is a zero-argument function that returns an iterable of
  training set batches. Internally, `learn!` uses this function when it calls
  `train!(model, get_train_batches(), logger)`.

- `get_test_batches` is a zero-argument function that returns an iterable
  of test set batches used during the current epoch's test phase. Each element of
  the iterable takes the form `(batch, votes_locations)`. Internally, `batch` is
  passed to [`loss_and_prediction`](@ref) as `loss_and_prediction(model, batch...)`,
  and `votes_locations[i]` is expected to yield the row index of `votes` that
  corresponds to the `i`th sample in `batch`.

- `votes` is a matrix of hard labels whose columns correspond to voters and whose
  rows correspond to the samples in the test set that have been voted on. If
  `votes[sample, voter]` is not a valid hard label for `model`, then `voter` will
  simply be considered to have not assigned a hard label to `sample`.

- `elected` is a vector of hard labels where the `i`th element is the hard label
  elected as "ground truth" out of `votes[i, :]`.

- `optimal_threshold_class` is the class index (`1` or `2`) for which to calculate
  an optimal threshold for converting `predicted_soft_labels` to `predicted_hard_labels`.
  This is only a valid parameter when `length(classes) == 2`. If `optimal_threshold_class`
  is present, test set evaluation will be based on predicted hard labels calculated
  with this threshold; if `optimal_threshold_class` is `nothing`, predicted hard labels
  will be calculated via `onecold(classifier, soft_label)`.
"""
function learn!(model::AbstractClassifier, logger, get_train_batches, get_test_batches,
                votes, elected=majority.(eachrow(votes), (1:length(classes(model)),));
                epoch_limit=100, post_epoch_callback=(_ -> nothing),
                optimal_threshold_class::Union{Nothing,Integer}=nothing,
                test_set_logger_prefix="test_set")
    # NOTE `votes` is currently unused except to construct `elected` by default,
    # but will be necessary for calculating multirater metrics e.g. Fleiss' kappa
    # later so we still required it in the API
    # TODO keeping `votes` alive might hog a lot of memory, so it would be convenient
    # to provide callers with a wrapper around `channel_unordered` (or at least example
    # code) for generating an iterator that batches `votes` rows to `soft labels` in
    # a manner that's respectful to throughput/memory constraints.
    # TODO is it better to pre-allocate these, or allocate them per-epoch? The
    # former ensures fixed memory usage, but the latter gives the GC a chance
    # to free up some RAM between epochs.
    _validate_threshold_class(optimal_threshold_class, classes(model))

    predicted = zeros(Float32, length(elected), length(classes(model)))
    log_event!(logger, "started learning")
    for current_epoch in 1:epoch_limit
        try
            train!(model, get_train_batches(), logger)
            predict!(model, predicted, get_test_batches(), logger;
                     logger_prefix="$(test_set_logger_prefix)_prediction")
            evaluate!(map(label -> onecold(model, label), eachrow(predicted)), predicted,
                      elected, classes(model), logger; logger_prefix="$(test_set_logger_prefix)_evaluation",
                      logger_suffix="_per_epoch", votes=votes,
                      optimal_threshold_class=optimal_threshold_class)
            post_epoch_callback(current_epoch)
            flush(logger)
        catch ex # support early stopping via exception handling
            if is_early_stopping_exception(model, ex)
                log_event!(logger, "`learn!` call stopped via $ex in epoch $current_epoch")
                break
            else
                log_event!(logger,
                           "`learn!` call encountered exception in epoch $current_epoch")
                rethrow(ex)
            end
        end
    end
    return model
end

#####
##### `post_epoch_callback` utilities
#####

"""
    upon(logged::Dict{String,Any}, field::AbstractString; condition, initial)

Return a closure that can be called to check the most recent state of
`logger.logged[field]` and trigger a caller-provided function when
`condition(recent_state, previously_chosen_state)` is `true`.

For example:

```
upon_loss_decrease = upon(logger, "test_set_prediction/mean_loss_per_epoch";
                          condition=<, initial=Inf)

save_upon_loss_decrease = _ -> begin
    upon_loss_decrease(new_lowest_loss -> save_my_model(model, new_lowest_loss),
                       consecutive_failures -> consecutive_failures > 10 && Flux.stop())
end

learn!(model, logger, get_train_batches, get_test_batches, votes;
       post_epoch_callback=save_upon_loss_decrease)
```

Specifically, the form of the returned closure is `f(on_true, on_false)` where
`on_true(state)` is called if `condition(state, previously_chosen_state)` is
`true`. Otherwise, `on_false(consecutive_falses)` is called where `consecutive_falses`
is the number of `condition` calls that have returned `false` since the last
`condition` call returned `true`.

Note that the returned closure is a no-op if `logger.logged[field]` has not
been updated since the most recent call.
"""
function upon(logged::Dict{String,Vector{Any}}, field::AbstractString; condition, initial)
    history = get!(() -> Any[], logged, field)
    previous_length = length(history)
    current = isempty(history) ? initial : last(history)
    consecutive_false_count = 0
    return (on_true, on_false=(_ -> nothing)) -> begin
        length(history) == previous_length && return nothing
        previous_length = length(history)
        candidate = last(history)
        if condition(candidate, current)
            consecutive_false_count = 0
            current = candidate
            on_true(current)
        else
            consecutive_false_count += 1
            on_false(consecutive_false_count)
        end
    end
end

function upon(logger::LearnLogger, field::AbstractString; condition, initial)
    return upon(logger.logged, field; condition=condition, initial=initial)
end


text:
"""
    GPU <: Processor

Abstract supertype for GPU processors.
"""
abstract type GPU <: Processor end

"""
    CPUGPU <: GPU

    CPUGPU()

Uses the CUDA GPU code on CPU using KernelAbstractions, to test it.
"""
struct CPUGPU{L} <: GPU
    spinlock::L
end
CPUGPU() = CPUGPU(Base.Threads.SpinLock())
Base.Threads.lock(opt::CPUGPU) = lock(opt.spinlock)
Base.Threads.unlock(opt::CPUGPU) = unlock(opt.spinlock)

kernel_setup(::CPUGPU) = KernelAbstractions.CPU(), 1

function maprule!(
    data::AbstractSimData, proc::GPU, opt, ruletype::Val{<:Rule}, rule, rkeys, wkeys
)
    kernel! = cu_cell_kernel!(kernel_setup(proc)...)
    kernel!(data, ruletype, rule, rkeys, wkeys; ndrange=gridsize(data)) |> wait
end
function maprule!(
    data::AbstractSimData, proc::GPU, opt,
    ruletype::Val{<:NeighborhoodRule}, rule, rkeys, wkeys
)
    kernel! = cu_neighborhood_kernel!(kernel_setup(proc)...)
    kernel!(data, opt, ruletype, rule, rkeys, wkeys, ndrange=gridsize(data)) |> wait
    return nothing
end

@kernel function cu_neighborhood_kernel!(
    data, opt, ruletype::Val{<:NeighborhoodRule}, rule, rkeys, wkeys
)
    I, J = @index(Global, NTuple)
    src = parent(source(_firstgrid(data, rkeys)))
    buf = _getwindow(src, neighborhood(rule), I, J)
    bufrule = _setbuffer(rule, buf)
    cell_kernel!(data, ruletype, bufrule, rkeys, wkeys, I, J)
    nothing
end

@kernel function cu_cell_kernel!(data, ruletype::Val, rule, rkeys, wkeys)
    i, j = @index(Global, NTuple)
    cell_kernel!(data, ruletype, rule, rkeys, wkeys, i, j)
    nothing
end

### Indexing. UNSAFE / LOCKS required

# This is not safe for general use.
# It can be used where only identical transformations of a cell
# can happen from any other cell, such as setting all 1s to 2.
@propagate_inbounds function _setindex!(d::WritableGridData, opt::GPU, x, I...)
    dest(d)[I...] = x
end

ext: jl


text:
mutable struct TSNode{T}
    value::T
    @atomic next::Union{TSNode{T},Nothing}
end

TSNode{T}(value::T) where {T} = TSNode{T}(value, nothing)

mutable struct ConcurrentStack{T}
    @atomic next::Union{TSNode{T},Nothing}
end

ConcurrentStack{T}() where {T} = ConcurrentStack{T}(nothing)

function Base.push!(stack::ConcurrentStack{T}, v) where {T}
    v = convert(T, v)
    node = TSNode{T}(v)

    next = @atomic stack.next
    while true
        @atomic node.next = next
        next, ok = @atomicreplace(stack.next, next => node)
        ok && break
    end

    return stack
end

function ConcurrentCollections.maybepop!(stack::ConcurrentStack)
    while true
        node = @atomic stack.next
        node === nothing && return nothing

        next = @atomic node.next
        next, ok = @atomicreplace(stack.next, node => next)
        if ok
            return Some(node.value)
        end
    end
end

function Base.pop!(stack::ConcurrentStack)
    r = maybepop!(stack)
    if r === nothing
        error("stack is empty")
    else
        return something(r)
    end
end


text:
################################################################################
# Submits test coverage for relevant code to Coveralls and Codecov
################################################################################

using Printf
using Coverage

# Add folders and jl files below if you want to exclude them from the test
# coverage reports. Folders are excluded recursively. All names are case
# insensitive, and relative to the src/ folder.
excluded_folders_and_files = map(f -> uppercase("src/" * f), [
    "gallery_extra/waveguide/waveguide_debug.jl",   # Some code to make larger verifications agains MATLAB (development phase)
    ])

coverage = process_folder()
unfiltered_count = length(coverage)

filter!(f -> !any(x -> startswith(uppercase(f.filename), x), excluded_folders_and_files), coverage)

covered_lines, total_lines = get_summary(coverage)
@printf("%d / %d lines covered by tests (%.2f %%), in %d source files (%d files excluded)\n",
    covered_lines, total_lines, covered_lines * 100.0 / total_lines, length(coverage), unfiltered_count - length(coverage))

#Coveralls.submit(coverage) # Coveralls submission is broken as of September 2018
Codecov.submit(coverage)


text:
using TableSchema
import TableSchema: read, validate

include("../src/DataPackage.jl")
using DataPackage
import DataPackage: Resource, add_resource, get_table, read

# Initialise a resource from a remote file
REMOTE_URL = "https://raw.githubusercontent.com/frictionlessdata/tableschema-jl/master/data/data_simple.csv"
res = Resource(REMOTE_URL, strict=false, name="sample")

# Add this resource to a blank data package
p = Package()
add_resource(p, res)

# Read the data from the table
data = read(res)
println( "The number of cells is ", length(data[:,1]) )
println( "Column 1 is called ", res.schema.fields[1].name )
println( join([ row for row in data[:,1] ], ", ", " and ") )

# Validate the table
table = get_table(res)
if validate(table); println("The table is valid according to the Schema"); end


text:
"""
    Index(exclude_decls_from_PCH, display_diagnostics)
Provide a shared context for creating translation units.

# Arguments
- `exclude_decls_from_PCH`: whether we only want to see "local" declarations (that did not come from a previous precompiled header). If false, we want to see all declarations.
- `display_diagnostics`: whether to display diagnostics.
"""
mutable struct Index
    ptr::CXIndex
    exclude_decls_from_PCH::Cint
    display_diagnostics::Cint
    function Index(exclude_decls_from_PCH, display_diagnostics)
        ptr = clang_createIndex(exclude_decls_from_PCH, display_diagnostics)
        @assert ptr != C_NULL
        obj = new(ptr, exclude_decls_from_PCH, display_diagnostics)
        finalizer(obj) do x
            if x.ptr != C_NULL
                clang_disposeIndex(x)
                x.ptr = C_NULL
            end
        end
        return obj
    end
end
Index(diagnostic::Bool) = Index(false, diagnostic)  #TODO: set `exclude_decls_from_PCH` on after adding PCH support
Index() = Index(true)

Base.unsafe_convert(::Type{CXIndex}, x::Index) = x.ptr

## TODO:
# clang_CXIndex_setGlobalOptions
# clang_CXIndex_getGlobalOptions
# clang_CXIndex_setInvocationEmissionPathOption


text:
using Test
using DexCall

@testset "DexCall" begin
    include("api.jl")
    include("evaluate.jl")
    include("native_function.jl")
end

max_stars_repo_name: justinchiu/dex-lang

max_issues_repo_name: justinchiu/dex-lang

max_forks_repo_name: justinchiu/dex-lang


text:
# code for Video 3_5:

## on command-line:
mongoimport -d iris -c lengths < iris.json
# 2015-11-30T15:10:38.525+0100    connected to: localhost
# 2015-11-30T15:10:38.661+0100    imported 150 documents

## in mongo shell:
> use iris
switched to db iris
> show collections
lengths
system.indexes
> db.lengths.count()
150
> db.lengths.findOne()
{
        "_id" : ObjectId("565c58de3449d12d21cefe41"),
        "sepalLength" : 5.1,
        "sepalWidth" : 3.5,
        "petalLength" : 1.4,
        "petalWidth" : 0.2,
        "species" : "setosa"
}

> db.lengths.find({species: "versicolor"}, {_id:0}).limit(5)
{ "sepalLength" : 7, "sepalWidth" : 3.2, "petalLength" : 4.7, "petalWidth" : 1.4,
 "species" : "versicolor" }
{ "sepalLength" : 6.4, "sepalWidth" : 3.2, "petalLength" : 4.5, "petalWidth" : 1.
5, "species" : "versicolor" }
{ "sepalLength" : 6.9, "sepalWidth" : 3.1, "petalLength" : 4.9, "petalWidth" : 1.
5, "species" : "versicolor" }
{ "sepalLength" : 5.5, "sepalWidth" : 2.3, "petalLength" : 4, "petalWidth" : 1.3,
 "species" : "versicolor" }
{ "sepalLength" : 6.5, "sepalWidth" : 2.8, "petalLength" : 4.6, "petalWidth" : 1.
5, "species" : "versicolor" }

Pkg.add("Mongo")
using Mongo, LibBSON

cl = MongoClient()
# 2015/11/30 19:54:37.0518: [ 9031]:    DEBUG:      cluster: Client initialized in direct mode.
# MongoClient(mongodb://localhost:27017/)

irisl = MongoCollection(cl, "iris", "lengths")
# MongoCollection(lengths)

spvir = BSONObject(Dict("species" => "virginica"))
# BSONObject({ "species" : "virginica" })

count(irisl, spvir) # 50

cur = find(irisl, spvir)
# Mongo.MongoCursor(Ptr{Void} @0x0000000003bd6620)

 for obj in cur
     println(obj["species"], ": ", obj["petalWidth"])
 end
# virginica: 2.5
# virginica: 1.9
# virginica: 2.1
# virginica: 1.8
# virginica: 2.2
# virginica: 2.1
# virginica: 1.7
# virginica: 1.8
# virginica: 1.8
# virginica: 2.5

new = BSONObject(Dict("species" => "virginica", "petalWidth" => 10.0))
# BSONObject({ "petalWidth" : 10.000000, "species" : "virginica" })

insert(irisl, new) # BSONOID(565d5f91f7294f259f22d551)
count(irisl, spvir) # 51

upd = BSONObject(Dict("species" => "virginica", "petalWidth" => 5.5))
# BSONObject({ "petalWidth" : 5.500000, "species" : "virginica" })
update(irisl, new, upd) # true

del = BSONObject(Dict("species" => "virginica", "petalWidth" => 5.5))
# BSONObject({ "petalWidth" : 5.500000, "species" : "virginica" })

delete(irisl, del) # true
count(irisl, spvir) # 50

## Redis:
Pkg.add("Redis")
using Redis

conn = RedisConnection()
# Redis.RedisConnection("127.0.0.1",6379,"",0,TCPSocket(open, 0 bytes waiting))

set(conn, "Iris species", "I. setosa") # true

keys(conn, *) # Set(Any["Iris species"])

println(get(conn, "Iris species")) # I. setosa
bgsave() # Ok

disconnect(conn)

max_stars_repo_path: Section 3/3_5.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 3/3_5.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 3/3_5.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
## R code 12.2
library(rethinking)
data(UCBadmit)
d <- UCBadmit
d$gid <- ifelse( d$applicant.gender=="male" , 1L , 2L )
dat <- list( A=d$admit , N=d$applications , gid=d$gid )
m12.1 <- ulam(
    alist(
        A ~ dbetabinom( N , pbar , theta ),
        logit(pbar) <- a[gid],
        a[gid] ~ dnorm( 0 , 1.5 ),
        transpars> theta <<- phi + 2.0,
        phi ~ dexp(1)
    ), data=dat , chains=4 )

## R code 12.6
library(rethinking)
data(Kline)
d <- Kline
d$P <- standardize( log(d$population) )
d$contact_id <- ifelse( d$contact=="high" , 2L , 1L )

dat2 <- list(
    T = d$total_tools,
    P = d$population,
    cid = d$contact_id )

m12.2 <- ulam(
    alist(
        T ~ dgampois( lambda , phi ),
        lambda <- exp(a[cid])*P^b[cid] / g,
        a[cid] ~ dnorm(1,1),
        b[cid] ~ dexp(1),
        g ~ dexp(1),
        phi ~ dexp(1)
    ), data=dat2 , chains=4 , log_lik=TRUE )

## R code 12.9
m12.3 <- ulam(
    alist(
        y ~ dzipois( p , lambda ),
        logit(p) <- ap,
        log(lambda) <- al,
        ap ~ dnorm( -1.5 , 1 ),
        al ~ dnorm( 1 , 0.5 )
    ) , data=list(y=y) , chains=4 )

## R code 12.11
m12.3_alt <- ulam(
    alist(
        y|y>0 ~ custom( log1m(p) + poisson_lpmf(y|lambda) ),
        y|y==0 ~ custom( log_mix( p , 0 , poisson_lpmf(0|lambda) ) ),
        logit(p) <- ap,
        log(lambda) <- al,
        ap ~ dnorm(-1.5,1),
        al ~ dnorm(1,0.5)
    ) , data=list(y=as.integer(y)) , chains=4 )

## R code 12.15
logit <- function(x) log(x/(1-x)) # convenience function
round( lco <- logit( cum_pr_k ) , 2 )

## R code 12.16
m12.4 <- ulam(
    alist(
        R ~ dordlogit( 0 , cutpoints ),
        cutpoints ~ dnorm( 0 , 1.5 )
    ) , data=list( R=d$response ), chains=4 , cores=4 )

## R code 12.17
m12.4q <- quap(
    alist(
        response ~ dordlogit( 0 , c(a1,a2,a3,a4,a5,a6) ),
        c(a1,a2,a3,a4,a5,a6) ~ dnorm( 0 , 1.5 )
    ) , data=d , start=list(a1=-2,a2=-1,a3=0,a4=1,a5=2,a6=2.5) )

## R code 12.24
dat <- list(
    R = d$response,
    A = d$action,
    I = d$intention,
    C = d$contact )
m12.5 <- ulam(
    alist(
        R ~ dordlogit( phi , cutpoints ),
        phi <- bA*A + bC*C + BI*I ,
        BI <- bI + bIA*A + bIC*C ,
        c(bA,bI,bC,bIA,bIC) ~ dnorm( 0 , 0.5 ),
        cutpoints ~ dnorm( 0 , 1.5 )
    ) , data=dat , chains=4 , cores=4 )

## R code 12.34
dat <- list(
    R = d$response ,
    action = d$action,
    intention = d$intention,
    contact = d$contact,
    E = as.integer( d$edu_new ),   # edu_new as an index
    alpha = rep( 2 , 7 ) )      # delta prior

m12.6 <- ulam(
    alist(
        R ~ ordered_logistic( phi , kappa ),
        phi <- bE*sum( delta_j[1:E] ) + bA*action + bI*intention + bC*contact,
        kappa ~ normal( 0 , 1.5 ),
        c(bA,bI,bC,bE) ~ normal( 0 , 1 ),
        vector[8]: delta_j <<- append_row( 0 , delta ),
        simplex[7]: delta ~ dirichlet( alpha )
    ), data=dat , chains=4 , cores=4 , coerce_int=FALSE )

## R code 12.37
dat$edu_norm <- normalize( d$edu_new )
m12.7 <- ulam(
    alist(
        R ~ ordered_logistic( mu , cutpoints ),
        mu <- bE*edu_norm + bA*action + bI*intention + bC*contact,
        c(bA,bI,bC,bE) ~ normal( 0 , 1 ),
        cutpoints ~ normal( 0 , 1.5 )
    ), data=dat , chains=4 , cores=4 )


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.DATA_DIR, "toy-data", "gdemo-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.MODELS_DIR, "toy-models", "gdemo.model.jl"))

@tbenchmark(HMC(20, 0.1, 3), simplegaussmodel, simplegaussstandata[1]["x"])
bench_res = @tbenchmark(HMC(2000, 0.1, 3), simplegaussmodel, simplegaussstandata[1]["x"])
logd = build_logd("Simple Gaussian Model", bench_res...)
logd["analytic"] = Dict("s" => 49/24, "m" => 7/6)

#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "gauss-stan.run.jl"))

#logd["stan"] = Dict("s" => mean(s_stan), "m" => mean(m_stan))
#logd["time_stan"] = sg_time

print_log(logd)
send_log(logd)


text:
using NullableArrays

srand(1)
N = 5_000_000
A = rand(N)
B = rand(Bool, N)
C = Array{Float64}(N)

X = NullableArray(A)
Y = NullableArray(A, B)
Z = NullableArray(Float64, N)

f(x::Float64) = 5 * x

function profile_map()

    println("Method: map!(f, dest, src) (0 missing entries)")
    print("  for Array{Float64}:          ")
    map!(f, C, A);
    @time map!(f, C, A);
    print("  for NullableArray{Float64}:  ")
    map!(f, Z, X);
    @time map!(f, Z, X);
    println()

    println("Method: map!(f, dest, src) (~half missing entries)")
    print("  for NullableArray{Float64}:  ")
    map!(f, Z, Y);
    @time map!(f, Z, Y);
    println()

    println("Method: map(f, src) (0 missing entries)")
    print("  for Array{Float64}:          ")
    map(f, A);
    @time map(f, A);
    print("  for NullableArray{Float64}:  ")
    map(f, X);
    @time map(f, X);
    println()

    println("Method: map(f, src) (~half missing entries)")
    print("  for NullableArray{Float64}:  ")
    map(f, Y);
    @time map(f, Y);
    println()
end


text:
using Flux, Gym, Printf
using Flux.Tracker: data, grad
using Flux.Optimise: _update_params!
using Statistics: mean
using DataStructures: CircularBuffer
using Distributions: sample

using CuArrays

#Load game environment

env = PendulumEnv()

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(reset!(env))
ACTION_SIZE = 1#length(env.actions)
ACTION_BOUND = 2#env.action_space.hi
MAX_EP = 50_000
MAX_EP_LENGTH = 1000

BATCH_SIZE = 64
MEM_SIZE = 100_000
MIN_EXP_SIZE = 50_000

γ = 99f-2     # discount rate

τ = 1f-3 # for running average while updating target networks
η_act = 1f-4   # Learning rate
η_crit = 1f-3
L2_DECAY = 1f-2

# Ornstein-Uhlenbeck Noise params
μ = 0f0
θ = 15f-2
σ = 2f-1

# --------------------------------- Memory ------------------------------------

memory = CircularBuffer{Any}(MEM_SIZE)

function getData(batch_size = BATCH_SIZE)
  # Getting data in shape
  minibatch = sample(memory, batch_size)
  x = hcat(minibatch...)

  s      =   hcat(x[1, :]...) |> gpu
  a      =   hcat(x[2, :]...) |> gpu
  r      =   hcat(x[3, :]...) |> gpu
  s′     =   hcat(x[4, :]...) |> gpu
  s_mask = .!hcat(x[5, :]...) |> gpu

  return s, a, r, s′, s_mask
end

# -------------------------------- Action Noise --------------------------------

struct OUNoise
  μ
  θ
  σ
  X
end

ou = OUNoise(μ, θ, σ, zeros(Float32, ACTION_SIZE) |> gpu)

function sample_noise(ou::OUNoise)
  dx     = ou.θ * (ou.μ .- ou.X)
  dx   .+= ou.σ * randn(Float32, length(ou.X)) |> gpu
  ou.X .+= dx
end

# Noise scale
τ_ = 25
ϵ  = exp(-1f0 / τ_)
noise_scale = 1f0 / ACTION_BOUND

# ----------------------------- Model Architecture -----------------------------

w_init(dims...) = 6f-3rand(Float32, dims...) .- 3f-3

actor = Chain(Dense(STATE_SIZE, 400, relu),
	      Dense(400, 300, relu),
              Dense(300, ACTION_SIZE, tanh, initW=w_init),
              x -> x * ACTION_BOUND) |> gpu
actor_target = deepcopy(actor)

struct crit
  state_crit
  act_crit
  sa_crit
end

Flux.@treelike crit

function (c::crit)(state, action)
  s = c.state_crit(state)
  a = c.act_crit(action)
  c.sa_crit(relu.(s .+ a))
end

Base.deepcopy(c::crit) = crit(deepcopy(c.state_crit),
                              deepcopy(c.act_crit),
			      deepcopy(c.sa_crit))

critic = crit(Chain(Dense(STATE_SIZE, 400, relu), Dense(400, 300)) |> gpu,
              Dense(ACTION_SIZE, 300) |> gpu,
	      Dense(300, 1, initW=w_init) |> gpu)
critic_target = deepcopy(critic)

# ------------------------------- Param Update Functions---------------------------------

function update_target!(target, model; τ = 1f0)
  for (p_t, p_m) in zip(params(target), params(model))
    p_t.data .= (1f0 - τ) * p_t.data .+ τ * p_m.data
  end
end

nullify_grad!(p) = p
nullify_grad!(p::TrackedArray) = (p.grad .= 0f0)

zero_grad!(model) = (model = mapleaves(nullify_grad!, model))

# ---------------------------------- Training ----------------------------------

function L2_loss(model)
  l2_loss = 0
  for p in params(model)
    l2_loss += sum(p .^ 2)
  end

  return L2_DECAY * l2_loss
end

opt_crit = ADAM(η_crit)
opt_act  = ADAM(η_act)

function replay()
  s, a, r, s′, s_mask = getData()

  # Update Critic
  a′ = data(actor_target(s′))
  v′ = data(critic_target(s′, a′))
  y = r .+ γ * v′ .* s_mask	# set v′ to 0 where s_ is terminal state

  v = critic(s, a)
  loss_crit = Flux.mse(y, v) + L2_loss(critic)

  zero_grad!(critic)
  Flux.back!(loss_crit)
  _update_params!(opt_crit, params(critic))

  # Update Actor
  actions = actor(s)
  crit_out = critic(s, actions)

  zero_grad!(actor)
  Flux.back!(-sum(crit_out))
  _update_params!(opt_act, params(actor))

  # Update Target models
  update_target!(actor_target, actor; τ = τ)
  update_target!(critic_target, critic; τ = τ)
end

# --------------------------- Helper Functions --------------------------------

# Stores tuple of state, action, reward, next_state, and done
remember(state, action, reward, next_state, done) =
  push!(memory, [data.((state, action, reward[1], next_state))..., done])

# Choose action according to policy PendulumPolicy
function action(state, train=true)
  state = reshape(data(state), size(state)..., 1)
  act_pred = data(actor(state |> gpu))
  if train
    act_pred .+= noise_scale * sample_noise(ou)
  end
  clamp.(act_pred, -ACTION_BOUND, ACTION_BOUND) # returns action
end

function episode!(env, train=true)
  total_reward = 0f0
  s = reset!(env)
  for ep=1:MAX_EP_LENGTH
    a = action(s, train)
    s′, r, done, _ = step!(env, a)
    total_reward += data(r)[1]
    if train
      remember(s, a, r, s′, done)
      replay()
    end
    s = s′
  end
  total_reward
end

# -------------------------------- Testing -------------------------------------

# Returns average score over 100 episodes
function test()
  score_mean = 0f0
  for _=1:100
    total_reward = episode!(env, false)
    score_mean += total_reward / 100
  end
  return score_mean
end

# ------------------------------ Training --------------------------------------

# Populate memory with random actions

s = reset!(env)
for e=1:MIN_EXP_SIZE
  global s
  a = 2rand(Float32) * ACTION_BOUND - ACTION_BOUND
  s′, r, done, _ = step!(env, a)
  remember(s, a, r, s′, done)
  s = s′
end

for e=1:MAX_EP
  global noise_scale
  total_reward = episode!(env)
  total_reward = @sprintf "%9.3f" total_reward
  print("Episode: $e | Score: $total_reward | ")
  score_mean = test()
  score_mean = @sprintf "%9.3f" score_mean
  println("Mean score over 100 test episodes: $score_mean")
  noise_scale *= ϵ
end


text: include("dqns/dqns.jl")


text:
##  Factor a  companion matrix
##  so  that  eigenvalue are the roots

## We have ps = [a0, a1, ..., an]
## we factor -(QR) so we use p(-x) coefficients here
function basic_decompose(ps::Vector{T}) where {T}

    #    ps, k = deflate_leading_zeros(ps)  ## ASSUMED
    n = length(ps) - 1
    p0, pNi = ps[1], 1/ps[end]
    qs = Vector{T}(undef, n+1)
    par = iseven(n) ? one(T) : -one(T)

    @inbounds for i in 1:n-1
        qs[i] =  par * (-1)^(i+1) * ps[i+1] * pNi
    end

    qs[n] =  par * p0 * pNi
    qs[n+1] = -one(T)

    qs

end

## A pencil function takes
## ps = [a0, a1, a2, ..., an] and splits into vs, ws with:
## v1 = a0,
## wn = an, and
## v_{i+1} + w_i = a_i for i-1,2,...,n-1.
##
## This gives matrices:
## V = [0 ....     -v1   W = [1  .... w1
##      1 0 .....  -v2        0 1 ... w2
##      0 1 .....  -v3        0 0 1 . w3
##          .....               .....
##      0 0 .... 1 -vn]       0 ....0 wn]
## V is Hessenberg, W upper triangular.
##
## The default, `basic_pencil`, is
## vs = [a0, a1, ..., a_{n-1}], ws = [0,0,..., an]
##
function basic_pencil(ps::Vector{T}) where {T}

    N = length(ps)-1

    qs = zeros(T, N)
    qs[N]  = ps[end]

    ps = ps[1:end-1]
    ps, qs

end

##################################################
##
## Factor R: (Z +  yt)
## xs are decompose(ps)

## Constructor for Real coefficients; no pencil
function r_factorization(xs::Vector{S}) where {S}
    N = length(xs) - 1
    T = S <: Real ? S : real(S)
    Ct = AscendingChain(Vector{Rotator{T,S}}(undef, N))
    B =  DescendingChain(Vector{Rotator{T,S}}(undef, N))
    D = SparseDiagonal(S, N+1)

    _r_factorization(xs, Ct, B, D)

    RFactorizationRankOne(Ct, B, D)
end

# constructor populating Ct, B, D for Real Case
function _r_factorization(xs::Vector{T}, Ct, B, D) where {T <: Real}
    N = length(xs) - 1
    c, s, tmp = givensrot(xs[N], -one(T))
    C = Rotator(c, s, N)
    Ct[1] = C' #Rotator(c, -s, N)
    B[N] = Rotator(s, -c, N)  # not the adjoint, so that Ct_1 * B_n = [0 -1; 1 0]

    @inbounds for i in (N-1):-1:1
        c,s,tmp = givensrot(xs[i], tmp)
        C = Rotator(c, s, i)
        Ct[N + 1 - i], B[i] = C', C
    end

    nothing
end

# constructor populating Ct, B, D for Complex Case (phase is impt)
function _r_factorization(xs::Vector{S}, Ct, B, D) where {S <: Complex}
    N = length(xs) - 1

    c, s, tmp = givensrot(xs[N], -one(S))

    nrm = norm(c)
    alpha = c/nrm

    C = Rotator(c, s, N)
    Ct[1] = C'
    B[N] = Rotator(s * alpha, -nrm, N)
    D[N] = conj(alpha)  # Ct*B*D gives [0 -1; 1 0] block at tail
    D[N+1] = alpha

    @inbounds for i in (N-1):-1:1
        c,s,tmp = givensrot(xs[i], tmp)
        C = Rotator(c, s, i)
        Ct[N + 1 - i], B[i] = C', C
    end

    nothing
end

## Factor Q using vector xs
function q_factorization(xs::Vector{S}) where {S}

    N = length(xs) - 1
    T = S <: Real ? S : real(S)

    Q =  DescendingChain(Vector{Rotator{T,S}}(undef, N-1))
    zt,ot,zs,os = zero(T), one(T), zero(S), one(S)

    @inbounds for ii = 1:(N-1)
        Q[ii] = Rotator(zero(S), one(T), ii)
    end

    q_factorization(Q)

end

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

"""

    amrvw(ps)

For a polynomial specified by `ps` computes a sparse factorization of its companion matrix.

"""
function amrvw(ps::Vector{S}) where {S}

    xs = basic_decompose(ps)

    QF = q_factorization(xs)
    RF = r_factorization(xs)

    QRFactorization(QF, RF)

end

## decomposition of ps into vs, ws needs to be expanded to work with framework
function adjust_pencil(vs::Vector{T}, ws::Vector{T}) where {T}

    ps = basic_decompose(vcat(vs, one(T)))

    ## adjust ws, as basic_decompose will use p(-x) for roots
    for i in length(ws)-1:-2:1
        ws[i] = -ws[i]
    end

    qs = vcat(ws, -one(T))  # or +1?

    ps, qs

end

"""

    amrvw(vs, ws)

Computes a sparse factorization of of the companion matrix of a polynomial specified througha  pencil decomposition.

A pencil decomposition of a polynomial, is a specification where
if `p = a0 + a1x^1 + ... + xn x^n`
then
`vs[1] = a0`,
`vs[i+1] + ws[i] = ai`, and
`ws[n] = an`.

"""
function amrvw(vs::Vector{S}, ws::Vector{S}) where {S}

    N = length(vs)

    ps, qs = adjust_pencil(vs, ws)

    QF = q_factorization(Vector{S}(undef, N+1)) # is this correct?
    ZF = pencil_factorization(ps, qs)

    QRFactorization(QF, ZF)

end

"""

    roots(ps)
    roots(vs, ws)

Use an algorithm of AMRVW to find roots of a polynomial over the reals of complex numbers.

ps: The coefficients, `[p0, p1, p2, ...., pn]`, of the polynomial `p0 + p1*x + p2*x^2 + ... + pn*x^n`

[vs, ws]: If a second set of coefficients is passed, a pencil decomposition is used.
A pencil decomposition satisfies `vs[1]=p0; vs[i+1]+ws[i] = pi and ws[n] = pn`.
A pencil can be used when there is a wide range in the coefficients of the polynomial, as though slower, it can be more stable.
The non-exported function `basic_pencil` implements the pencil which pulls off the leading coefficient `a_n`.

Returns a complex vector of roots. When the algorithm fails, a warning is issued about the number of non-identified roots.

Examples:

```
using AMRVW; const A = AMRVW
rs = rand(10)
A.roots(rs)      # uses Real double shift algorithm

rs = rand(Complex{Float64}, 10)
A.roots(rs)      # uses Complex Single Shift algorithm

rs = rand(10)
vs, ws = A.basic_pencil(rs)
A.roots(vs, ws)  # uses QZ pencil factorization
```

References:

* Fast and backward stable computation of the eigenvalues of matrix polynomials. By Aurentz, Jared & Mach, Thomas & Robol, Leonardo & Vandebril, Raf & Watkins, David. (2016). Mathematics of Computation. 88. DOI: 10.1090/mcom/3338.

* Fast and backward stable computation of roots of polynomials, Part II: backward error analysis; companion matrix and companion pencil
Jared L. Aurentz, Thomas Mach, Leonardo Robol, Raf Vandebril, David S. Watkins; arXiv:1611.02435

"""
function roots(ps::Vector{S}) where {S}
    # deflate 0s
    ps, K = deflate_leading_zeros(ps)

    if length(ps) <= 3
        rts = solve_simple_cases(ps)
        append!(rts, zeros(Complex{real(S)}, K))
        return rts
    end

    state = amrvw(ps)

    EIGS =  AMRVW_algorithm!(state)

    if K > 0 # put back in leading zeros
        append!(EIGS, zeros(eltype(EIGS), K))
    end

    EIGS

end

## A Pencil algorithm
## if p = a0 + a1x^1 + ... + xn x^n
## then
## vs[1] = a0
## vs[i+1] + ws[i] = ai
## ws[n] = an
## is a pencil decompostion.
## `basic_pencil` will break [a0, a1, ..., an] into vs=[a0,a1, ..., a_{n-1}], ws=[0,0,..., 0, an]
function roots(vs::Vector{S}, ws::Vector{S}) where {S}

    ps, qs, K = deflate_leading_zeros(vs, ws)

    ## solve simple cases
    N = length(ps)
    if N <= 2
        rs = vcat(ps, zero(S)) + vcat(zero(S), qs)
        rts = roots(rs)
        append!(rts, zeros(Complex{real(S)}, K))
        return rts
    end

    state = amrvw(ps, qs)
    EIGS =  AMRVW_algorithm!(state)

    if K > 0 # put back in leading zeros
        append!(EIGS, zeros(eltype(EIGS), K))
    end

    EIGS
end


text:
# ---
# title: 943. Find the Shortest Superstring
# id: problem943
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/find-the-shortest-superstring/description/>
# hidden: true
# ---
#
# Given an array A of strings, find any smallest string that contains each
# string in `A` as a substring.
#
# We may assume that no string in `A` is substring of another string in `A`.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["alex","loves","leetcode"]
#     Output: "alexlovesleetcode"
#     Explanation: All permutations of "alex","loves","leetcode" would also be accepted.
#
#
# **Example 2:**
#
#
#
#     Input: ["catg","ctaagt","gcta","ttca","atgcatc"]
#     Output: "gctaagttcatgcatc"
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 12`
#   2. `1 <= A[i].length <= 20`
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
challenge_str = let input = @__DIR__() * "/day16_input.txt"
    isfile(input) ? readline(input) : test_str1
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Parts 1 and 2

function read_packet(code, k)
    ver = parse(Int, code[k:k+2], base=2)
    id = parse(Int, code[k+3:k+5], base=2)
    k += 6
    if id == 4
        a = '1'
        b = ""
        while a == '1'
            a = code[k]
            b *= code[k+1:k+4]
            k += 5
        end
        val = parse(Int, b, base=2)
    else
        bl = code[k] == '1' ? 11 : 15
        cl = parse(Int, code[k+1:k+bl], base=2)
        k += bl + 1
        vals = Int[]
        if bl == 11 # cl contains the number of subpackets
            for _ in 1:cl
                loc_ver, loc_val, k = read_packet(code, k)
                ver += loc_ver
                push!(vals, loc_val)
            end
        else # cl contains the total length in bits of the subpackets
            k0 = k
            while k < k0 + cl
                loc_ver, loc_val, k = read_packet(code, k)
                ver += loc_ver
                push!(vals, loc_val)
            end
        end
        if id == 0
            val = sum(vals)
        elseif id == 1
            val = prod(vals)
        elseif id == 2
            val = minimum(vals)
        elseif id == 3
            val = maximum(vals)
        elseif id == 5
            val = first(vals) > last(vals) ? 1 : 0
        elseif id == 6
            val = first(vals) < last(vals) ? 1 : 0
        elseif id == 7
            val = first(vals) == last(vals) ? 1 : 0
        else
            @warn "id type not implemented"
        end
    end
    return ver, val, k
end

function solve(str) # solves both parts
    code = prod(
        string.(
            digits(
                parse(
                    BigInt,
                    str,
                    base=16
                ),
                base=2,
                pad=4*length(str)
            )[end:-1:1]'
        )
    )
    k = 1
    ver, val, = read_packet(code, k)
    return ver, val
end

testit("Test Part 1", solve("8A004A801A8002F478")[1], ==(16))
testit("Test Part 1", solve("620080001611562C8802118E34")[1], ==(12))
testit("Test Part 1", solve("C0015000016115A2E0802F182340")[1], ==(23))
testit("Test Part 1", solve("A0016C880162017C3686B18A3D4780")[1], ==(31))
testit("Challenge Part 1", solve(challenge_str)[1], ==(925))

testit("Test Part 2", solve("C200B40A82")[2], ==(3))
testit("Test Part 2", solve("04005AC33890")[2], ==(54))
testit("Test Part 2", solve("880086C3E88112")[2], ==(7))
testit("Test Part 2", solve("CE00C43D881120")[2], ==(9))
testit("Test Part 2", solve("D8005AC2A8F0")[2], ==(1))
testit("Test Part 2", solve("F600BC2D8F")[2], ==(0))
testit("Test Part 2", solve("9C005AC2F8F0")[2], ==(0))
testit("Test Part 2", solve("9C0141080250320F1802104A08")[2], ==(1))
testit("Challenge Part 2", solve(challenge_str)[2], ==(342997120375))

max_stars_repo_name: rmsrosa/adventofcode

max_issues_repo_name: rmsrosa/adventofcode

max_forks_repo_name: rmsrosa/adventofcode


text:
"""
    filled_img = imfill(img::AbstractArray{Bool}, interval; dims=coords_spatial(img))
    filled_img = imfill(img::AbstractArray{Bool}, interval, connectivity)

Connected components of an image is found using flood-fill algorithm and returns a copy of
the original image after filling objects that falls in the range of interval.
For filling objects, represent the holes (part to be filled) with `true` in your array.

Parameters:
-  img            = Input image (Boolean array type)
-  interval       = objects of size (# of voxels) in this range will be filled with `false`
-  connectivity   = a Boolean-valued connectivity pattern, see [`label_components`](@ref).

# Examples

```jldoctest; setup=:(using ImageMorphology)
julia> img = Bool[0 0 1 1 0 0;
                  0 1 0 1 1 0;
                  0 0 1 1 0 0]
3×6 Matrix{Bool}:
 0  0  1  1  0  0
 0  1  0  1  1  0
 0  0  1  1  0  0

julia> imfill(.!(img), 0:3)
3×6 BitMatrix:
 1  1  0  0  1  1
 1  0  0  0  0  1
 1  1  0  0  1  1

julia> .!(ans)
3×6 BitMatrix:
 0  0  1  1  0  0
 0  1  1  1  1  0
 0  0  1  1  0  0
```
"""
function imfill(img::AbstractArray{Bool}, interval::Tuple{Real,Real}, connectivity::AbstractArray{Bool})
    return _imfill(img, interval, label_components(img, connectivity))
end
function imfill(img::AbstractArray{Bool}, interval::Tuple{Real,Real}; dims=coords_spatial(img))
    return _imfill(img, interval, label_components(img; dims=dims))
end
function imfill(img::AbstractArray{Bool}, interval, args...; kwargs...)
    return imfill(img, (minimum(interval)::Real, maximum(interval)::Real), args...; kwargs...)
end

function _imfill(img::AbstractArray{Bool}, interval::Tuple{Real,Real}, labels)
    if interval[1] > interval[2] || interval[1] < 0 || interval[2] < 0
        msg = "Interval must be non-negative and in format (min_range,max_range)"
        throw(DomainError(interval, msg))
    end
    count_labels = component_lengths(labels)

    new_img = similar(img)
    for ind in eachindex(img)
        if img[ind] == true && interval[1] <= count_labels[labels[ind]] <= interval[2]
            new_img[ind] = false
        else
            new_img[ind] = img[ind]
        end
    end

    return new_img
end


text:
function parse_species(ft::RSSANetwork, fname)
    specs_ic = Dict{Symbol,Int}()

    open(fname, "r") do file
        for line in eachline(file)
            var,val = split(line, "=")
            specs_ic[Symbol(strip(var))] = parse(Int64, val)
        end
    end

    specs_ic
end

function parse_reactions(ft::RSSANetwork, fname)
    rxstrs  = Vector{String}()
    rxrates = Vector{Float64}()

    open(fname, "r") do file
        for line in eachline(file)
            rxstr,rate = split(line, ",")
            rxstr = replace(strip(rxstr), "-" => "--")
            rxstr = replace(strip(rxstr), "_" => "0")
            push!(rxstrs, strip(rxstr))
            push!(rxrates, parse(Float64, rate))
        end
    end

    rxstrs, rxrates
end

function build_rxnetwork(ft::RSSANetwork, networkname, rxstrs, rxrates; printrxs=false, kwargs...)

    # string representing the network
    rxiobuf = IOBuffer()
    write(rxiobuf, "@min_reaction_network $(networkname) begin\n")
    for (i,rxstr) in enumerate(rxstrs)
        rate = rxrates[i]
        write(rxiobuf, "\t $rate, $rxstr \n")
    end
    write(rxiobuf, "end\n")
    rnstr = String(take!(rxiobuf))

    if printrxs
        print(rnstr)
    end

    # build the network using DiffEqBiological
    rn = eval( Meta.parse(rnstr) )

    rn,rnstr
end

function get_init_condit(ft::RSSANetwork, rn, specs_ic)
    icvec = zeros(Int64, length(specs_ic))
    for (i,sym) in enumerate(rn.syms)
        icvec[i] = specs_ic[sym]
    end

    icvec
end

# for parsing the simple format from the book by Thanh et al.
function loadrxnetwork(ft::RSSANetwork, networkname, specs_ic_file, rxs_file; kwargs...)

    # parse initial conditions
    specs_ic = parse_species(ft, specs_ic_file)

    # parse reaction network
    rxstrs,rxrates = parse_reactions(ft, rxs_file)

    # build the DiffEqBiological representation of the network
    rn,rnstr = build_rxnetwork(ft, networkname, rxstrs, rxrates; kwargs...)
    initialpop = get_init_condit(ft, rn, specs_ic)

    ParsedReactionNetwork(rn, initialpop; rnstr=rnstr)
end

max_issues_repo_name: isaacsas/ParseRxNetwork


text:
@testset "1072.flip-columns-for-maximum-number-of-equal-rows.jl" begin
    @test max_equal_rows_after_flips(hcat([[0, 1], [1, 1]]...)) == 1
    @test max_equal_rows_after_flips(hcat([[0, 1], [1, 0]]...)) == 2
    @test max_equal_rows_after_flips(hcat([[0, 0, 0], [0, 0, 1], [1, 1, 0]]...)) == 2
end


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

"""
    PointSet(points)

A set of `points` (a.k.a. point cloud) seen as a single entity.

## Examples

All point sets below are the same and contain two points in R³:

```julia
julia> PointSet([Point(1,2,3), Point(4,5,6)])
julia> PointSet(Point(1,2,3), Point(4,5,6))
julia> PointSet([(1,2,3), (4,5,6)])
julia> PointSet((1,2,3), (4,5,6))
julia> PointSet([[1,2,3], [4,5,6]])
julia> PointSet([1,2,3], [4,5,6])
julia> PointSet([1 4; 2 5; 3 6])
```
"""
struct PointSet{Dim,T} <: Domain{Dim,T}
  points::Vector{Point{Dim,T}}
end

PointSet(points::Vararg{P}) where {P<:Point} = PointSet(collect(points))
PointSet(coords::AbstractVector{TP}) where {TP<:Tuple} = PointSet(Point.(coords))
PointSet(coords::Vararg{TP}) where {TP<:Tuple} = PointSet(collect(coords))
PointSet(coords::AbstractVector{V}) where {V<:AbstractVector} = PointSet(Point.(coords))
PointSet(coords::Vararg{V}) where {V<:AbstractVector} = PointSet(collect(coords))
PointSet(coords::AbstractMatrix) = PointSet(Point.(eachcol(coords)))

"""
    PointSet(domain)

A point set representation of the `domain`, i.e. a set of centroid points
for each element of the `domain`.
"""
PointSet(domain::Domain) = PointSet([centroid(domain, ind) for ind in 1:nelements(domain)])

# -----------------
# DOMAIN INTERFACE
# -----------------

Base.getindex(pset::PointSet, ind::Int) = getindex(pset.points, ind)

nelements(pset::PointSet) = length(pset.points)

centroid(pset::PointSet, ind::Int) = pset.points[ind]


text:
module Runtests

import ModiaResult
using  Test

@testset "Test ModiaPlot_GLMakie/test" begin
    ModiaResult.usePlotPackage("GLMakie")
    include("$(ModiaResult.path)/test/runtests_withPlot.jl")
    ModiaResult.usePreviousPlotPackage()
end

end


text:
function tournament_selection(pop, fitness, cfg)
    tournament = rand(1:length(fitness), cfg.tournament_size)
    winner = tournament[argmin(fitness[tournament])]
    return pop[winner]
end

random_selection(pop) = rand(pop)

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
using StochasticDiffEq, DiffEqBase, Random

function f(du,u,p,t)
  for i in 1:length(u)
    du[i] = (0.2/length(u))*u[i]
  end
end

function g(du,u,p,t)
  for i in 1:length(u)
    du[i] = (0.2/length(u))*u[i]
  end
end

function condition(u,t,integrator)
  1-maximum(u)
end

function callback_affect!(integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()
  u[maxidx] = Θ
  u[end] = 1-Θ
  nothing
end

callback = ContinuousCallback(condition,callback_affect!)

u0 = [0.2]
tspan = (0.0,100.0)
prob = SDEProblem(f,g,u0,tspan)

Random.seed!(3)
sol = solve(prob,SRIW1(),callback=callback)

#=
using Plots; pyplot()
p1 = plot(sol,vars=(0,1),plotdensity=10000,title="Amount of X in Cell 1")
scatter!(sol,denseplot=false)
p2 = plot(sol.t,map((x)->length(x),sol[:]),lw=3,
     ylabel="Number of Cells",xlabel="Time")
plot(p1,p2,layout=(2,1),size=(600,1000))
=#

Random.seed!(3)
sol = solve(prob,EM(),callback=callback,dt=1/4)
sol = solve(prob,RKMil(),callback=callback,dt=1/4)

function g(du,u,p,t)
  for i in 1:length(u)
    du[i] = (0.3/length(u))
  end
end
prob = SDEProblem(f,g,u0,tspan)

sol = solve(prob,SRA1(),callback=callback)


text:
using ReproductionNumbers

# setup
fig_name    =   "reported"
file_name   =   "../data/raw/reported-cases-data.csv" # this file is huge and is hence added to .gitignore
days_col    =   "Refdatum"
data_col    =   "AnzahlFall_sum"
ylabel_R    =   "Reproduction number R"
ylabel_N    =   "Reported cases"
# The reported cases by RKI are changed and updated daily
# That means, cases that were not submitted on time will be added the day after etc.
# This leads to today's and yesterday's values of R and N being rather useless
# Hence, the variable `pop_days` allows to go back as many days from today as we believe are credible
# Looking at the RKI dashboard we find that today - 2 days is fine.
pop_days    =   2

N, R, df_cases = main_reported(fig_name, file_name, days_col, data_col, ylabel_R, ylabel_N, "state", pop_days)

# # ## break down data for German states
# df = CSV.File(file_name) |> DataFrame

# df_filter = filter_to_relevant_cases(df)
# gdf = groupby(df_filter, :Bundesland)
# cases_dict = Dict{String, DataFrame}()

# for (i, df_state) in enumerate(gdf)
#     # get name of current state
#     state_name = unique(gdf[i].Bundesland) |> first
#     # aggregate state data by days
#     gdf_state = groupby(df_state, days_col)
#     # sum over the number of cases & sort
#     cases = combine(gdf_state, :AnzahlFall => sum) |> sort
#     # promote date column to Dates.DateFormat
#     cases[!, days_col] = Date.(cases[!, days_col], Dates.DateFormat("yyyy/mm/dd H:M:S"))
#     # standardize column names
#     df_case = get_reference_data(cases, days_col = days_col, data_col = data_col, kind = "cases")
#     df_case = df_case[df_case.days .>= Date("2020-03-15"), :]
#     cases_dict[state_name] = df_case
# end

# cases_dict = sort(cases_dict)


text:
importall DynamicDiscreteModels
#need to `importall` rather than `using` because we will extend `coef!()` and `dim()`
#in real applications this would typically be handled in a front-end package ToyModel.jl

import ForwardDiff

type ToyModel <: DynamicDiscreteModel

	#DynamicDiscreteModel fields
	m::Array{Float64,4}			  	#the transition matrix given as m[x,y,x',y']
	mu::Array{Float64,2}  			#initial distribution (dx,dy)

	#DynamicDiscreteModel's technical fields
	mjac::Array{Float64,5}			#jacobian
	rho::Array{Float64,1}
	phi::Array{Float64,1}
	psi::Array{Float64,1}
	rhojac::Array{Float64,1}
	phijac::Array{Float64,2}
	psijac::Array{Float64,2}
end

dx=2
dy=3

toymodel()=ToyModel(
	Array(Float64,dx,dy,dx,dy), 			#m
	fill(1/6,2,3), 							#mu
	Array(Float64,dx,dy,dx,dy,2),			#mjac
	Array(Float64,1), 						#rho
	Array(Float64,dx),						#phi
	Array(Float64,dx),						#psi
	Array(Float64,2), 						#rhojac
	Array(Float64,dx,2),					#phijac
	Array(Float64,dx,2)						#psijac
	)

function coef!(model::ToyModel,theta::Tuple)
	p1,p2=theta[1],theta[2]
	a=[p1 1-p1;1-p1 p1]
	p3=(1-p2)/2
	b=[p2 p3 p3; p3 p3 p2]
	hmm2ddm!(model,a,b)
end

theta2eta(theta::Tuple)=[log(theta[1]/(1-theta[1])),log(theta[2]/(1-theta[2]))]
eta2theta(eta::Array)=(exp(eta[1])/(1+exp(eta[1])),exp(eta[2])/(1+exp(eta[2])))
coef!(model::ToyModel,eta::Array)=coef!(model,eta2theta(eta))

dim(model::ToyModel)=2

function eta2ab(eta)
	theta=eta2theta(eta)
	a=[theta[1],1-theta[1],1-theta[1],theta[1]]
	p2=theta[2]
	p3=(1-p2)/2
	b=[p2, p3, p3, p3, p3, p2]
	(a,b)
end

#quick implementation of derivatives with numerical gradient
#explicit expressions should be easy to implement
function coef_jac!(model::ToyModel,eta)
	a,b=eta2ab(eta)
	fa(eta)=eta2ab(eta)[1]
	fb(eta)=eta2ab(eta)[2]
	ajac=reshape(ForwardDiff.jacobian(fa)(eta),2,2,2)
	bjac=reshape(ForwardDiff.jacobian(fb)(eta),2,3,2)
	hmm2ddm!(model::DynamicDiscreteModel,reshape(a,2,2),reshape(b,2,3),ajac,bjac)
end

# Plotting function used for the README plot

function lkcontour()
	theta0=(.65, .5)
	model=toymodel()
	coef!(model,theta0)
	data=rand(model,100,100)
	thetahat=eta2theta(mle(model,data))
	thetahat2=eta2theta(em(model,data))
	thetahat3=eta2theta(mle(model,data[1:10]))

	xx=linspace(.1,.9,50)
	yy=linspace(.1,.9,50)
	function ff(x,y,n)
		coef!(model,(x,y))
		loglikelihood(model,data[1:n])
	end
	gg1=[ff(x,y,10)::Float64 for x=xx,y=yy]
	gg2=[ff(x,y,100)::Float64 for x=xx,y=yy]

	qlevels=[0,.25,.5,.75,0.8,0.9,0.95,0.975,0.99,1]
	levels=quantile(vec(gg),qlevels)

	contour(xx,yy,gg,fill=true,levels=levels,color=ColorGradient(:heat,[0,0.98,1]))
	scatter!([theta0[1]],[theta0[2]],c=:green,leg=false)
	scatter!([thetahat[1]],[thetahat[2]],c=:blue,leg=false)
	contour(xx,yy,gg,fill=true,levels=levels,color=ColorGradient(:heat,[0,0.98,1]))
	scatter!([theta0[1]],[theta0[2]],c=:green,leg=false)
	scatter!([thetahat3[1]],[thetahat3[2]],c=:blue,leg=false)
end


text:
#=
Popis
  Funkcia vypocita pre kazdy vrchol hodnotu jeho load centrality.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String
  normalize : parameter pre aplikovanie normalizacie :: Bool

Navratova hodnota
  loadness : slovnik obsahujuci mapovanie vrcholov a k nim prisluchajucich hodnot
           closeness centrality (index_vrchola => hodnota_loadness_centrality) :: Dict{Int64, Float64}()

Poznamky
  Hodnota centrality moze byt normalizovana nastavenim parametra normalize na True/False.
=#

function centrality_load(g::GenericGraph; normalize::Bool=true, dist_key::AbstractString="weight")

  loadness = [v.index => 0.0 for v=vertices(g)]

  for s in keys(loadness)

    u_loadness = node_loadness(g, s, dist_key)

    for temp in keys(u_loadness)
      loadness[temp] += u_loadness[temp]
    end

  end

  if normalize

    vertex_count = num_vertices(g)

    if vertex_count <= 2
      return loadness
    end

    scale = 1.0 / ((vertex_count - 1) * (vertex_count - 2))

    for v in keys(loadness)
      loadness[v] *= scale
    end

  end

  return loadness

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
import Gadfly

function plot1D(gp::GP; clim::Tuple{Float64, Float64}=(minimum(gp.X), maximum(gp.X)), CI::Float64=1.96, res::Int=1000)

        sx = (clim[2]-clim[1])/(res-1)
        x=collect(clim[1]:sx:clim[2])
        mu, Sigma = predict(gp, x)
        conf = CI*sqrt(Sigma)
        u = mu + conf
        l = mu - conf

        Gadfly.plot(Gadfly.layer(x=x, y=mu, ymin=l, ymax=u, Gadfly.Geom.line, Gadfly.Geom.ribbon),
                Gadfly.layer(x=gp.X,y=gp.y,Gadfly.Geom.point))
end

function plot2D(gp::GP; clim::Tuple{Float64, Float64, Float64, Float64} = (minimum(gp.X[1,:]), maximum(gp.X[1,:]),
                                                                   minimum(gp.X[2,:]), maximum(gp.X[2,:])),
                res::Int=50)
        sx = (clim[2]-clim[1])/(res-1)
        sy = (clim[4]-clim[3])/(res-1)
        A = Array(Float64,2,res^2)
        for (i,a) in enumerate(clim[1]:sx:clim[2])
            for (j,b) in enumerate(clim[3]:sy:clim[4])
                A[:,(i-1)*res+j] = [a,b]
            end
        end

        mu = predict(gp,A)[1]
        z= reshape(mu,res,res)
        Gadfly.plot(z=z,x=collect(clim[1]:sx:clim[2]),y=collect(clim[3]:sy:clim[4]),Gadfly.Geom.contour)
end

function Gadfly.plot(gp::GP; kwargs...)
    d, n = size(gp.X)
    if d>2
        error("Can only plot one or two dimensional Gaussian processes")
    elseif d==1
        plot1D(gp; kwargs...)
    elseif d==2
        plot2D(gp; kwargs...)
    end
end


text:
using ValueHistories

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

tests = [
    "tst_history.jl"
    "tst_mvhistory.jl"
]

perf = [
    "bm_history.jl"
]

for t in tests
    @testset "[->] $t" begin
        include(t)
    end
end

for p in perf
    @testset "[->] $p" begin
        include(p)
    end
end


text:
using Base.LinAlg: axpy!, scale!

function update!(r::Reg; gscale=1, o...)
    initupdate(r)
    getp(r,:l1reg,0)!=0 && l1reg!(getp(r,:l1reg), r.out, r.dif)
    getp(r,:l2reg,0)!=0 && l2reg!(getp(r,:l2reg), r.out, r.dif)
    getp(r,:adagrad) && adagrad!(1e-8, getp(r,:ada), r.dif) # TODO: make 1e-8 a parameter
    getp(r,:rmsprop) && rmsprop!(1e-8, 0.9, getp(r,:rms), r.dif) # TODO: make 1e-8 a parameter
    getp(r,:momentum,0)!=0 && momentum!(getp(r,:momentum), getp(r,:mom), r.dif)
    getp(r,:nesterov,0)!=0 && nesterov!(getp(r,:nesterov), getp(r,:nes), r.dif)
    if getp(r, :adam)
        incp(r, :t)
        adam!(getp(r, :t), getp(r, :eps), getp(r, :b1), getp(r, :b2), getp(r, :fstm), getp(r, :scndm), r.dif)
    end
    scale = -1 * getp(r,:lr,1) * gscale # TODO: make scale a parameter for callbacks
    axpy!(scale, r.dif, r.out)
    getp(r,:average) && axpy!(1,r.out,getp(r,:avg))
end

function initupdate(r::Reg)
    isa(r.op, Par) || error("Update called on $(typeof(r.op))")
    if r.dif != nothing
        getp(r,:adagrad) && similarp(r, :ada, r.dif)
        getp(r,:rmsprop) && similarp(r, :rms, r.dif)
        getp(r,:momentum,0)!=0 && similarp(r, :mom, r.dif)
        getp(r,:nesterov,0)!=0 && similarp(r, :nes, r.dif)
        if getp(r,:adam)
            similarp(r, :fstm, r.dif)
            similarp(r, :scndm, r.dif)
            getp(r, :eps) == 0 && setp(r, :eps, 1e-8)
            getp(r, :b1) == 0 && setp(r, :b1, 0.9)
            getp(r, :b2) == 0 && setp(r, :b2, 0.999)
        end
    else
        error("r.dif==nothing in $r")
    end
    if r.out != nothing
        getp(r,:average) && similarp(r, :avg, r.out)
    else
        error("r.out==nothing in $r")
    end
end

function similarp(r,n,a)
    b = getp(r,n)
    issimilar(a,b) && return b
    b != false && warn("Changing $n")
    setp(r,n,fillsync!(similar(a),0))
end

l1reg!(l1, w, dw)=for i=1:length(dw); (w[i]>zero(w[i])) ? (dw[i]+=l1) : (w[i]<zero(w[i])) ? (dw[i]-=l1) : 0; end
l2reg!(l2, w, dw)=axpy!(l2, w, dw)
adagrad!(eps, dw2, dw)=for i=1:length(dw); dw2[i] += dw[i] * dw[i]; dw[i] /= sqrt(dw2[i] + eps); end
momentum!(m, dw2, dw)=(axpy!(m, dw2, dw); copysync!(dw2,dw))
nesterov!(m, dw2, dw)=(scale!(m, dw2); axpy!(1, dw, dw2); axpy!(m, dw2, dw))
rmsprop!(eps, rho, dw2, dw)=for i=1:length(dw); dw2[i] = dw2[i] * rho + (1 - rho) * dw[i] * dw[i]; dw[i] /= sqrt(dw2[i] + eps); end

function adam!(t, eps, b1, b2, fstm, scndm, dw)
    for i=1:length(dw)
        fstm[i] = b1*fstm[i] + (1-b1)*dw[i]
        scndm[i] = b2*scndm[i] + (1-b2)*(dw[i] *dw[i])
        dw[i] = (fstm[i] / (1 - b1 ^ t)) / (sqrt(scndm[i] / (1 - b2 ^ t)) + eps)
    end
end

@gpu adagrad!(eps, dw2::CudaArray{Float32}, dw::CudaArray{Float32})=ccall((:adagrad32,libknet),Void,(Cint,Cdouble,Ptr{Float32},Ptr{Float32}),length(dw),eps,dw2,dw)
@gpu adagrad!(eps, dw2::CudaArray{Float64}, dw::CudaArray{Float64})=ccall((:adagrad64,libknet),Void,(Cint,Cdouble,Ptr{Float64},Ptr{Float64}),length(dw),eps,dw2,dw)
@gpu l1reg!(l1, w::CudaArray{Float32}, dw::CudaArray{Float32})=ccall((:l1reg32,libknet),Void,(Cint,Cdouble,Ptr{Float32},Ptr{Float32}),length(dw),l1,w,dw)
@gpu l1reg!(l1, w::CudaArray{Float64}, dw::CudaArray{Float64})=ccall((:l1reg64,libknet),Void,(Cint,Cdouble,Ptr{Float64},Ptr{Float64}),length(dw),l1,w,dw)
@gpu rmsprop!(eps, rho, dw2::CudaArray{Float32}, dw::CudaArray{Float32})=ccall((:rmsprop32,libknet),Void,(Cint,Cdouble, Cdouble,Ptr{Float32},Ptr{Float32}),length(dw),eps, rho, dw2,dw)
@gpu rmsprop!(eps, rho, dw2::CudaArray{Float64}, dw::CudaArray{Float64})=ccall((:rmsprop64,libknet),Void,(Cint,Cdouble,Cdouble,Ptr{Float64},Ptr{Float64}),length(dw),eps, rho, dw2,dw)
@gpu adam!(t, eps, b1, b2, fstm::CudaArray{Float32}, scndm::CudaArray{Float32}, dw::CudaArray{Float32})=ccall((:adam32,libknet),Void,(Cint,Cint,Cdouble,Cdouble, Cdouble,Ptr{Float32},Ptr{Float32},Ptr{Float32}),length(dw),t,eps,b1,b2,fstm,scndm,dw)
@gpu adam!(t, eps, b1, b2, fstm::CudaArray{Float64}, scndm::CudaArray{Float64}, dw::CudaArray{Float64})=ccall((:adam64,libknet),Void,(Cint,Cint,Cdouble,Cdouble, Cdouble,Ptr{Float64},Ptr{Float64},Ptr{Float64}),length(dw),t,eps,b1,b2,fstm,scndm,dw)

# function maxnorm!(maxnorm, w)
#     error("Did not debug maxnorm yet.")
#     norms = sqrt(sum(w.^2, 2))
#     if any(norms > p.maxnorm)
#         scale = min(p.maxnorm ./ norms, 1)
#         w *= scale
#     end
# end

# nz(r,n,v)=(isdefined(r,n) && (r.(n) != v))


text:
module Utils

export @hascuda

const HAVE_CUDA = try
  using CUDAnative
  using CUDAdrv
  true
catch
  false
end
if HAVE_CUDA
  macro hascuda(ex)
    return :($(esc(ex)))
  end
else
  macro hascuda(ex)
    return :()
  end
end

end


text:
using BenchmarkTools

list_test_str = split(
"7,4,9,5,11,17,23,2,0,14,21,24,10,16,13,6,15,25,12,22,18,20,8,19,3,26,1

22 13 17 11  0
 8  2 23  4 24
21  9 14 16  7
 6 10  3 18  5
 1 12 20 15 19

 3 15  0  2 22
 9 18 13 17  5
19  8  7 25 23
20 11 10 24  4
14 21 16 12  6

14 21 17 24  4
10 16 15  9 19
18  8 23 26 20
22 11 13  6  5
 2  0 12  3  7",
'\n'
)

list_str = readlines("day4_input.txt")

# Part 1

list = copy(list_test_str)

function part1(list)
    numbers = parse.(Int, split(list[1], ','))
    num_boards = div(length(list) - 1, 6)
    boards = [Matrix{Int}(undef, 5, 5) for _ in 1:num_boards]
    masks = [BitMatrix(zeros(Bool, 5, 5)) for _ in 1:num_boards]

    for k in 1:num_boards
        for i in 1:5
            boards[k][i, :] .= parse.(Int, split(list[2 + 6 * (k - 1) + i]))
        end
    end

    round = 0
    winner = 0
    while winner == 0
        round += 1
        n = numbers[round]
        for k in 1:num_boards
            masks[k] .|= boards[k] .== n
            for i in 1:5
                if all(masks[k][i,:]) || all(masks[k][:,i])
                    winner = k
                end
            end
        end
    end

    return numbers[round] * sum((!).(masks[winner]) .* boards[winner])
end

#end

@info "Part 1"
@info "test:"
@show part1(list_test_str) == 4512
@info "challenge:"
@show part1(list_str) == 10374

# Part 2

function part2(list)
    numbers = parse.(Int, split(list[1], ','))
    num_boards = div(length(list) - 1, 6)
    boards = [Matrix{Int}(undef, 5, 5) for _ in 1:num_boards]
    masks = [BitMatrix(zeros(Bool, 5, 5)) for _ in 1:num_boards]

    for k in 1:num_boards
        for i in 1:5
            boards[k][i, :] .= parse.(Int, split(list[2 + 6 * (k - 1) + i]))
        end
    end

    round = 0
    winners = Int[]
    while length(winners) < num_boards
        round += 1
        n = numbers[round]
        for k in 1:num_boards
            masks[k] .|= boards[k] .== n
            if k ∉ winners
                for i in 1:5
                    if all(masks[k][i,:]) || all(masks[k][:,i])
                        push!(winners, k)
                        break
                    end
                end
            end
        end
    end

    return numbers[round] * sum((!).(masks[winners[end]]) .* boards[winners[end]])
end

@info "Part 2"
@info "test:"
@show part2(list_test_str) == 1924
@info "challenge:"
@show part2(list_str) == 24742

# Splitting up things

function init(list)
    numbers = parse.(Int, split(list[1], ','))
    num_boards = div(length(list) - 1, 6)
    boards = [Matrix{Int}(undef, 5, 5) for _ in 1:num_boards]
    masks = [BitMatrix(zeros(Bool, 5, 5)) for _ in 1:num_boards]

    for k in 1:num_boards
        for i in 1:5
            boards[k][i, :] .= parse.(Int, split(list[2 + 6 * (k - 1) + i]))
        end
    end
    return numbers, boards, masks
end

function part1!(numbers, boards, masks)
    num_boards = length(boards)
    round = 0
    winner = 0
    @inbounds while winner == 0
        round += 1
        n = numbers[round]
        for k in 1:num_boards
            masks[k] .|= boards[k] .== n
            for i in 1:5
                if all(@view(masks[k][i,:])) || all(@view(masks[k][:,i]))
                    winner = k
                end
            end
        end
    end

    return numbers[round] * sum((!).(masks[winner]) .* boards[winner])
end

function part2!(numbers, boards, masks)
    num_boards = length(boards)
    round = 0
    winners = Int[]
    @inbounds while length(winners) < num_boards
        round += 1
        n = numbers[round]
        for k in 1:num_boards
            masks[k] .|= boards[k] .== n
            if k ∉ winners
                for i in 1:5
                    if all(@view(masks[k][i,:])) || all(@view(masks[k][:,i]))
                        push!(winners, k)
                        break
                    end
                end
            end
        end
    end

    return numbers[round] * sum((!).(masks[winners[end]]) .* boards[winners[end]])
end

@info "Part 1 - splitted"
@info "test:"
numbers, boards, masks = init(list_test_str)
@show part1!(numbers, boards, masks) == 4512
@info "challenge:"
numbers, boards, masks = init(list_str)
@show part1!(numbers, boards, masks) == 10374

@info "Part 2 - splitted"
@info "test:"
numbers, boards, masks = init(list_test_str)
@show part2!(numbers, boards, masks) == 1924
@info "challenge:"
numbers, boards, masks = init(list_str)
@show part2!(numbers, boards, masks) == 24742

@info "Benchmarks"
@info "Init"
@btime init($list_str);
@info "Part 1"
@btime part1!($numbers, b, m) setup=(b=boards; m=masks) evals=1
@info "Part 2"
@btime part2!($numbers, b, m) setup=(b=boards; m=masks) evals=1


text:
using DifferentialEquations, Plots

𝞪 = 2.0
u0 = 1/2
f(t,u) = -𝞪*sin(u)
tspan = (0.0, 12.0)

prob = ODEProblem(f, u0, tspan)

sol1 = solve(prob, Euler(), dt=3.0)
sol2 = solve(prob, Euler(), dt=0.5)
sol3 = solve(prob, Euler(), dt=0.1)

#[t+2u for (t,u) in zip(sol1.t,sol.u)]

x = 0.0:0.05:12.0
p = plot(x, sol1(x))
plot!(p, x, sol2(x))
plot!(p, x, sol3(x))

gui()