text:
module LULESH

using CEnum
using MPI
using Parameters
using Printf
using Random
using StaticArrays
using LinearAlgebra

abstract type AbstractDomain end
# Device vector types
# const VD{T} = Union{Vector{T}, CuVector{T}} where T
const VD{T} = Vector{T} where T

# Index type
const IndexT = Int
# MPI Message Tags
const MSG_COMM_SBN      = 1024
const MSG_SYNC_POS_VEL  = 2048
const MSG_MONOQ         = 3072
const MAX_FIELDS_PER_MPI_COMM = 6
# Assume 128 byte coherence
# Assume Float64 is an "integral power of 2" bytes wide
const CACHE_COHERENCE_PAD_REAL = div(128, sizeof(Float64))

# Only supported configuration
const SEDOV_SYNC_POS_VEL_EARLY = true
const ALLOW_UNPACKED_PLANE = false

struct LuleshProblem
    num_iters::Int
    structured::Bool
    nx::Int
    nr::Int
    balance::Int
    cost::Int
    col::Int
    row::Int
    plane::Int
    side::Int
    devicetype
    floattype
    comm::Union{MPI.Comm, Nothing}
    function LuleshProblem(num_iters, structured, nx, nr, balance, cost, devicetype, floattype, comm)
        col, row, plane, side = initMeshDecomp(comm)
        return new(num_iters, structured, nx, nr, balance, cost, col, row, plane, side, devicetype, floattype, comm)
    end
end

include("bc.jl")
include("types.jl")
include("domain.jl")
include("mpi.jl")
include("comm.jl")
include("utils.jl")

export printUsage, IndexT, Domain, LuleshProblem, getMyRank, getNumRanks, getWtime,
       lagrangeLeapFrog, comm_max, timeIncrement!, nodalMass, commRecv, MSG_COMM_SBN, verifyAndWriteFinalOutput
export commSend, commRecv, commSBN

function initMeshDecomp(comm)
    # Assume cube processor layout for now
    numRanks = getNumRanks(comm)
    myRank = getMyRank(comm)
    testProcs = floor(cbrt(numRanks+0.5))
    if (testProcs*testProcs*testProcs != numRanks)
        error("Num processors must be a cube of an integer (1, 8, 27, ...)")
    end
    # TODO This is not good for the padding
    if (MAX_FIELDS_PER_MPI_COMM > getCacheCoherencePad(testProcs))
        error("corner element comm buffers too small. MAX_FIELDS_PER_MPI_COMM > CACHE_COHERENCE_PAD_REAL ($MAX_FIELDS_PER_MPI_COMM > $(getCacheCoherencePad(testProcs))")
    end
    testProcs = convert(Int, testProcs)
    dx = convert(Int, testProcs)
    dy = convert(Int, testProcs)
    dz = convert(Int, testProcs)

    # temporary test
    if dx*dy*dz != numRanks
        error("error -- must have as many domains as procs\n") ;
    end
    remainder = (dx*dy*dz) % numRanks
    myDom = Int(0)
    if (myRank < remainder)
        myDom = myRank*( 1+ (div(dx*dy*dz, numRanks)))
    else
      myDom = remainder*( 1+ (div(dx*dy*dz, numRanks))) +
         (myRank - remainder)*(div(dx*dy*dz, numRanks))
    end
    col = myDom % dx
    row = div(myDom, dx) % dy
    plane = div(myDom, (dx*dy))
    side = testProcs

    return col, row, plane, side
end
export InitMeshDecomp

function printNormField(domain, fields)
    if getMyRank(domain.comm) == 0
        for field in fields
            # println("$field[", length(getfield(domain, field)), "]: ", norm(getfield(domain, field)))
        end
    end
end

function printField(domain, fields)
    if getMyRank(domain.comm) == 0
        for field in fields
            # println("$field[", length(getfield(domain, field)), "]: ", getfield(domain, field))
        end
    end
end

function printNormAllFields(domain, location="")
    if getMyRank(domain.comm) == 0
        # println("Location: ", location)
    end
    # fields = [:x, :xd, :fx, :nodalMass, :symmX, :dxx, :delv_xi]
    # fields = [:x, :commDataSend, :commDataRecv]
    # printNormField(domain, fields)
end

function printAllFields(domain, location="")
    if getMyRank(domain.comm) == 0
        # println("Location: ", location)
    end
    # fields = [:x, :xd, :fx, :nodalMass, :symmX, :nodelist]
    # fields = [:commDataSend, :commDataRecv]
    # printField(domain, fields)
end

export printField, printNormField, printAllFields, printNormAllFields
end # module


text:
# ---------------------------------------------------
## CLASSIFICATION METRICS (PROBABILISTIC PREDICTIONS)
# >> CrossEntropy
# >> BriersScore

# -----------------------------------------------------
# cross entropy

struct CrossEntropy{R} <: Measure where R <: AbstractFloat
    eps::R
end
CrossEntropy(;eps=eps()) = CrossEntropy(eps)

metadata_measure(CrossEntropy;
    name                     = "cross_entropy",
    target_scitype           = Vec{<:Finite},
    prediction_type          = :probabilistic,
    orientation              = :loss,
    reports_each_observation = true,
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "Cross entropy loss with probabilities " *
                               "clamped between eps and 1-eps; aliases: " *
                               "`cross_entropy`.",
    distribution_type        = UnivariateFinite)

"""
    cross_entropy

$(docstring(CrossEntropy()))

    CrossEntropy(; eps=eps())
    ce(ŷ, y)

Given an abstract vector of distributions `ŷ` and an abstract vector
of true observations `y`, return the corresponding Cross-Entropy
loss (aka log loss) scores.

Since the score is undefined in the case of the true observation has
predicted probability zero, probablities are clipped between `eps` and
`1-eps` where `eps` can be specified.

If `sᵢ` is the predicted probability for the true class `yᵢ` then
the score for that example is given by

    -log(clamp(sᵢ, eps, 1-eps))

For more information, run `info(cross_entropy)`.
"""
cross_entropy = CrossEntropy()

# for single observation:
_cross_entropy(d, y, eps) = -log(clamp(pdf(d, y), eps, 1 - eps))

function (c::CrossEntropy)(ŷ::Vec{<:UnivariateFinite},
                           y::Vec{<:CategoricalElement})
    check_dimensions(ŷ, y)
    check_pools(ŷ, y)
    return broadcast(_cross_entropy, ŷ, y, c.eps)
end

# -----------------------------------------------------
# brier score

# TODO: support many distributions/samplers D below:

struct BrierScore{D} <: Measure end

# As this measure is parametric, the use of `metadata_measure` is not
# appropriate.

MLJModelInterface.name(::Type{<:BrierScore{D}}) where D =
    "BrierScore{$(string(D))}"
MLJModelInterface.docstring(::Type{<:BrierScore{D}}) where D =
    "Brier proper scoring rule for distributions of type $D; "*
    "aliases: `BrierScore{$D}`"
MLJModelInterface.docstring(::Type{<:BrierScore{<:UnivariateFinite}}) =
    "Brier proper scoring rule for `MultiClass` or `OrderedFactor` data; "*
    "aliases: `BrierScore()`, `BrierScore{UnivariateFinite}`"
MLJModelInterface.docstring(::Type{BrierScore}) =
    "Brier proper scoring rule for various distribution types; " *
    "use `brier_score` for `BrierScore{UnivariateFinite}` "*
"(`Multiclass` or `OrderedFactor` targets)."
MLJModelInterface.target_scitype(::Type{<:BrierScore{D}}) where D =
    AbstractVector{<:Finite}
MLJModelInterface.prediction_type(::Type{<:BrierScore}) = :probabilistic
orientation(::Type{<:BrierScore}) = :score
reports_each_observation(::Type{<:BrierScore}) = true
is_feature_dependent(::Type{<:BrierScore}) = false
MLJModelInterface.supports_weights(::Type{<:BrierScore}) = true
distribution_type(::Type{<:BrierScore{D}}) where D =
    UnivariateFinite

"""
    BrierScore(; distribution=UnivariateFinite)(ŷ, y)

Given an abstract vector of distributions `ŷ` of type `distribution`,
and an abstract vector of true observations `y`, return the
corresponding Brier (aka quadratic) scores.

Currently only `distribution=UnivariateFinite` is supported, which is
applicable to superivised models with `Finite` target scitype. In this
case, if `p(y)` is the predicted probability for a *single*
observation `y`, and `C` all possible classes, then the corresponding
Brier score for that observation is given by

``2p(y) - \\left(\\sum_{η ∈ C} p(η)^2\\right) - 1``

Note that `BrierScore()=BrierScore{UnivariateFinite}` has the alias
`brier_score`.

*Warning.* `BrierScore` defines a true *score* (bigger is better). In
Brier's original 1950 paper, and many other places, it is defined as a
loss, despite the name. Moreover, the present implementation does not
treat the binary case as special, so that the score may differ, in
that case, by a factor of two from usage elsewhere.

For more information, run `info(BrierScore)`.

"""
function BrierScore(; distribution=UnivariateFinite)
    distribution == UnivariateFinite ||
        error("Only `UnivariateFinite` Brier scores currently supported. ")
    return BrierScore{distribution}()
end

# For single observations (no checks):

# UnivariateFinite:
function _brier_score(d::UnivariateFinite, y)
    levels = classes(d)
    pvec = broadcast(pdf, d, levels)
    offset = 1 + sum(pvec.^2)
    return 2 * pdf(d, y) - offset
end

# For multiple observations:

# UnivariateFinite:
function (::BrierScore{<:UnivariateFinite})(ŷ::Vec{<:UnivariateFinite},
                                            y::Vec{<:CategoricalElement})
    check_dimensions(ŷ, y)
    check_pools(ŷ, y)
    return broadcast(_brier_score, ŷ, y)
end

function (score::BrierScore{<:UnivariateFinite})(ŷ, y, w::Vec{<:Real})
    check_dimensions(y, w)
    return w .* score(ŷ, y) ./ (sum(w)/length(y))
end

const brier_score = BrierScore()

# ============================================================
## CLASSIFICATION METRICS (DETERMINISTIC PREDICTIONS)

const INVARIANT_LABEL_MULTICLASS = "This metric is invariant to class labelling and can be used for multiclass classification."
const INVARIANT_LABEL_BINARY = "This metric is invariant to class labelling and can be used only for binary classification."
const VARIANT_LABEL_BINARY = "This metric is labelling-dependent and can only be used for binary classification."

# ==============================================================
## MULTICLASS

# ---------------------------------------------------
# misclassification rate

struct MisclassificationRate <: Measure end

metadata_measure(MisclassificationRate;
    name                     = "misclassification_rate",
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "misclassification rate; aliases: " *
                               "`misclassification_rate`, `mcr`.")

"""
    misclassification_rate

$(docstring(MisclassificationRate()))

    misclassification_rate(ŷ, y)
    misclassification_rate(ŷ, y, w)
    misclassification_rate(conf_mat)

Returns the rate of misclassification of the (point) predictions `ŷ`,
given true observations `y`, optionally weighted by the weights
`w`. All three arguments must be abstract vectors of the same length.
A confusion matrix can also be passed as argument.
$INVARIANT_LABEL_MULTICLASS

For more information, run `info(misclassification_rate)`.

"""
const misclassification_rate = MisclassificationRate()
const mcr = misclassification_rate
const MCR = MisclassificationRate

(::MCR)(ŷ::Vec{<:CategoricalElement},
        y::Vec{<:CategoricalElement}) = mean(y .!= ŷ)

(::MCR)(ŷ::Vec{<:CategoricalElement},
        y::Vec{<:CategoricalElement},
        w::Vec{<:Real}) = sum((y .!= ŷ) .* w) / sum(w)

(::MCR)(cm::ConfusionMatrix) = 1.0 - sum(diag(cm.mat)) / sum(cm.mat)

# -------------------------------------------------------------
# accuracy

struct Accuracy <: Measure end

metadata_measure(Accuracy;
    name                     = "accuracy",
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :score,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "Classification accuracy; aliases: `accuracy`.")

"""
    accuracy

$(docstring(Accuracy()))

    accuracy(ŷ, y)
    accuracy(ŷ, y, w)
    accuracy(conf_mat)

Returns the accuracy of the (point) predictions `ŷ`,
given true observations `y`, optionally weighted by the weights
`w`. All three arguments must be abstract vectors of the same length.
$INVARIANT_LABEL_MULTICLASS

For more information, run `info(accuracy)`.

"""
const accuracy = Accuracy()

(::Accuracy)(args...) = 1.0 - misclassification_rate(args...)
(::Accuracy)(m::ConfusionMatrix) = sum(diag(m.mat)) / sum(m.mat)

# -----------------------------------------------------------
# balanced accuracy

struct BalancedAccuracy <: Measure end

metadata_measure(BalancedAccuracy;
    name                     = "balanced_accuracy",
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :score,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = true,
    docstring                = "Balanced classification accuracy; aliases: "*
                               "`balanced_accuracy`, `bacc`, `bac`.")

"""
    balanced_accuracy

$(docstring(BalancedAccuracy()))

    balanced_accuracy(ŷ, y [, w])
    balanced_accuracy(conf_mat)

Return the balanced accuracy of the point prediction `ŷ`, given true
observations `y`, optionally weighted by `w`. The balanced accuracy takes
into consideration class imbalance.
All  three arguments must have the same length.
$INVARIANT_LABEL_MULTICLASS

For more information, run `info(balanced_accuracy)`.

"""
const balanced_accuracy = BalancedAccuracy()
const bacc = balanced_accuracy
const bac  = bacc
const BACC = BalancedAccuracy

function (::BACC)(ŷ::Vec{<:CategoricalElement},
                  y::Vec{<:CategoricalElement})
    class_count = Dist.countmap(y)
    ŵ = 1.0 ./ [class_count[yi] for yi in y]
    return sum( (ŷ .== y) .* ŵ ) / sum(ŵ)
end

function (::BACC)(ŷ::Vec{<:CategoricalElement},
                  y::Vec{<:CategoricalElement},
                  w::Vec{<:Real})
    levels_ = levels(y)
    ŵ = similar(w)
    @inbounds for i in eachindex(w)
        ŵ[i] = w[i] / sum(w .* (y .== y[i]))
    end
    return sum( (ŷ .== y) .* ŵ ) / sum(ŵ)
end

# ==================================================================
## BINARY AND ORDER-INDEPENDENT

# ------------------------------------------------------------------
# Matthew's correlation

struct MatthewsCorrelation <: Measure end

metadata_measure(MatthewsCorrelation;
    name                     = "matthews_correlation",
    target_scitype           = Vec{<:Finite{2}},
    prediction_type          = :deterministic,
    orientation              = :score,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "Matthew's correlation; aliases: " *
                               "`matthews_correlation`, `mcc`")

"""
    matthews_correlation

$(docstring(MatthewsCorrelation()))

    matthews_correlation(ŷ, y)
    matthews_correlation(conf_mat)

Return Matthews' correlation coefficient corresponding to the point
prediction `ŷ`, given true observations `y`.
$INVARIANT_LABEL_MULTICLASS

For more information, run `info(matthews_correlation)`.

"""
const matthews_correlation = MatthewsCorrelation()
const mcc = matthews_correlation
const MCC = MatthewsCorrelation

function (::MCC)(cm::ConfusionMatrix{C}) where C
    # http://rk.kvl.dk/introduction/index.html
    # NOTE: this is O(C^3), there may be a clever way to
    # speed this up though in general this is only used for low  C
    num = 0
    @inbounds for k in 1:C, l in 1:C, m in 1:C
        num += cm[k,k] * cm[l,m] - cm[k,l] * cm[m,k]
    end
    den1 = 0
    den2 = 0
    @inbounds for k in 1:C
        a = sum(cm[k, :])
        b = sum(cm[setdiff(1:C, k), :])
        den1 += a * b
        a = sum(cm[:, k])
        b = sum(cm[:, setdiff(1:C, k)])
        den2 += a * b
    end
    mcc = num / sqrt(den1 * den2)

    isnan(mcc) && return 0
    return mcc
end

(m::MCC)(ŷ::Vec{<:CategoricalElement},
         y::Vec{<:CategoricalElement}) =
             confmat(ŷ, y, warn=false) |> m

# ---------------------------------------------------------
# area under the ROC curve

struct AUC <: Measure end

metadata_measure(AUC;
    name                     = "area_under_curve",
    target_scitype           = Vec{<:Finite},
    prediction_type          = :probabilistic,
    orientation              = :score,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = false,
    docstring                = "Area under the ROC curve; "*
       "aliases: `area_under_curve`, `auc`")

"""
    area_under_curve

$(docstring(AUC()))

    area_under_curve(ŷ, y)

Return the area under the receiver operator characteristic (curve),
for probabilistic predictions `ŷ`, given ground truth `y`.
$INVARIANT_LABEL_BINARY

For more information, run `info(area_under_curve)`.

"""
const area_under_curve = AUC()
const auc = AUC()

function (::AUC)(ŷ::Vec{<:UnivariateFinite},
                 y::Vec{<:CategoricalElement})
    # implementation drawn from https://www.ibm.com/developerworks/community/blogs/jfp/entry/Fast_Computation_of_AUC_ROC_score?lang=en
    lab_pos = levels(y)[2]         # 'positive' label
    scores  = pdf.(ŷ, lab_pos)     # associated scores
    y_sort  = y[sortperm(scores)]  # sort by scores
    n       = length(y)
    n_neg   = 0  # to keep of the number of negative preds
    auc     = 0
    @inbounds for i in 1:n
        # y[i] == lab_p --> it's a positive label in the ground truth
        # in that case increase the auc by the cumulative sum
        # otherwise increase the number of negatives by 1
        δ_auc, δ_neg = ifelse(y_sort[i] == lab_pos, (n_neg, 0), (0, 1))
        auc   += δ_auc
        n_neg += δ_neg
    end
    n_pos = n - n_neg
    return auc / (n_neg * n_pos)
end

# ==========================================================================
## BINARY AND ORDER DEPENDENT

const CM2 = ConfusionMatrix{2}

# --------------------------------------------------------------------------
# F_β-Score

"""
    FScore{β}(rev=nothing)

One-parameter generalization, ``F_β``, of the F-measure or balanced F-score.

[Wikipedia entry](https://en.wikipedia.org/wiki/F1_score)

    FScore{β}(ŷ, y)

Evaluate ``F_β`` score on observations ,`ŷ`, given ground truth values, `y`.

By default, the second element of `levels(y)` is designated as
`true`. To reverse roles, use `FScore{β}(rev=true)` instead of
`FScore{β}`.

For more information, run `info(FScore)`.

"""
struct FScore{β} <: Measure rev::Union{Nothing,Bool} end

FScore{β}(; rev=nothing) where β = FScore{β}(rev)

metadata_measure(FScore;
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :score,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = false)

MLJModelInterface.name(::Type{<:FScore{β}}) where β = "FScore{$β}"
MLJModelInterface.name(::Type{FScore})            = "FScore" # for registry
MLJModelInterface.docstring(::Type{<:FScore})       = "F_β score; aliases: " *
                                        "`FScore{β}`, `f1score=FScore{1}`"
const f1score      = FScore{1}()

# -------------------------------------------------------------------------
# truepositive, true_negative, etc

for M in (:TruePositive, :TrueNegative, :FalsePositive, :FalseNegative,
          :TruePositiveRate, :TrueNegativeRate, :FalsePositiveRate,
          :FalseNegativeRate, :FalseDiscoveryRate, :Precision, :NPV)
    ex = quote
        struct $M <: Measure rev::Union{Nothing,Bool} end
        $M(; rev=nothing) = $M(rev)
    end
    eval(ex)
end

metadata_measure.((FalsePositive, FalseNegative);
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    aggregation              = Sum(),
    is_feature_dependent     = false,
    supports_weights         = false)

metadata_measure.((FalsePositiveRate, FalseNegativeRate, FalseDiscoveryRate);
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :loss,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = false)

metadata_measure.((TruePositive, TrueNegative);
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :score,
    reports_each_observation = false,
    aggregation              = Sum(),
    is_feature_dependent     = false,
    supports_weights         = false)

metadata_measure.((TruePositiveRate, TrueNegativeRate, Precision, NPV);
    target_scitype           = Vec{<:Finite},
    prediction_type          = :deterministic,
    orientation              = :score,
    reports_each_observation = false,
    is_feature_dependent     = false,
    supports_weights         = false)

# adjustments
MMI.name(::Type{<:TruePositive})       = "true_positive"
MMI.docstring(::Type{<:TruePositive})  = "Number of true positives; " *
                                         "aliases: `true_positive`, `truepositive`."
MMI.name(::Type{<:TrueNegative})       = "true_negative"
MMI.docstring(::Type{<:TrueNegative})  = "Number of true negatives; " *
                                         "aliases: `true_negative`, `truenegative`."
MMI.name(::Type{<:FalsePositive})      = "false_positive"
MMI.docstring(::Type{<:FalsePositive}) = "Number of false positives; " *
                                         "aliases: `false_positive`, `falsepositive`."
MMI.name(::Type{<:FalseNegative})      = "false_negative"
MMI.docstring(::Type{<:FalseNegative}) = "Number of false negatives; " *
                                         "aliases: `false_negative`, `falsenegative`."

MMI.name(::Type{<:TruePositiveRate})      = "true_positive_rate"
MMI.docstring(::Type{<:TruePositiveRate}) = "True positive rate; aliases: " *
                               "`true_positive_rate`, `truepositive_rate`, `tpr`, `sensitivity`, " *
                               "`recall`, `hit_rate`."
MMI.name(::Type{<:TrueNegativeRate})      = "true_negative_rate"
MMI.docstring(::Type{<:TrueNegativeRate}) = "true negative rate; aliases: " *
                               "`true_negative_rate`, `truenegative_rate`, `tnr`, `specificity`, " *
                               "`selectivity`."
MMI.name(::Type{<:FalsePositiveRate})      = "false_positive_rate"
MMI.docstring(::Type{<:FalsePositiveRate}) = "false positive rate; aliases: " *
                               "`false_positive_rate`, `falsepositive_rate`, `fpr`, `fallout`."
MMI.name(::Type{<:FalseNegativeRate})      = "false_negative_rate"
MMI.docstring(::Type{<:FalseNegativeRate}) = "false negative rate; aliases: " *
                               "`false_negative_rate`, `falsenegative_rate`, `fnr`, `miss_rate`."
MMI.name(::Type{<:FalseDiscoveryRate})      = "false_discovery_rate"
MMI.docstring(::Type{<:FalseDiscoveryRate}) = "false discovery rate; "*
                               "aliases: `false_discovery_rate`, `falsediscovery_rate`, `fdr`."
MMI.name(::Type{<:NPV})      = "negative_predictive_value"
MMI.docstring(::Type{<:NPV}) = "negative predictive value; aliases: " *
                               "`negative_predictive_value`, `negativepredictive_value`, `npv`."

MMI.name(::Type{<:Precision})         = "positive_predictive_value"
MMI.docstring(::Type{<:Precision})    = "positive predictive value "*
  "(aka precision); aliases: `positive_predictive_value`, `ppv`, `Precision()`, `positivepredictive_value`. "

"""
    true_positive

$(docstring(TruePositive()))

    true_positive(ŷ, y)

Number of true positives for observations `ŷ` and ground truth
`y`. Assigns `false` to first element of `levels(y)`. To reverse roles,
use `TruePositive(rev=true)` instead of `true_positive`.

For more information, run `info(true_positive)`.

"""
const true_positive = TruePositive()
const tp = TruePositive()
const truepositive  = TruePositive()

"""
    true_negative

$(docstring(TrueNegative()))

    true_negative(ŷ, y)

Number of true negatives for observations `ŷ` and ground truth
`y`. Assigns `false` to first element of `levels(y)`. To reverse roles,
use `TrueNegative(rev=true)` instead of `true_negative`.

For more information, run `info(true_negative)`.

"""
const true_negative = TrueNegative()
const tn = TrueNegative()
const truenegative  = TrueNegative()

"""
    false_positive

$(docstring(FalsePositive()))

    false_positive(ŷ, y)

Number of false positives for observations `ŷ` and ground truth
`y`. Assigns `false` to first element of `levels(y)`. To reverse roles,
use `FalsePositive(rev=true)` instead of `false_positive`.

For more information, run `info(false_positive)`.

"""
const false_positive = FalsePositive()
const fp = FalsePositive()
const falsepositive = FalsePositive()

"""
    false_negative

$(docstring(FalseNegative()))

    false_negative(ŷ, y)

Number of false positives for observations `ŷ` and ground truth
`y`. Assigns `false` to first element of `levels(y)`. To reverse roles,
use `FalseNegative(rev=true)` instead of `false_negative`.

For more information, run `info(false_negative)`.

"""
const false_negative = FalseNegative()
const fn = FalseNegative()
const falsenegative = FalseNegative()

"""
    true_positive_rate

$(docstring(TruePositiveRate()))

    true_positive_rate(ŷ, y)

True positive rate for observations `ŷ` and ground truth `y`. Assigns
`false` to first element of `levels(y)`. To reverse roles, use
`TruePositiveRate(rev=true)` instead of `true_positive_rate`.

For more information, run `info(true_positive_rate)`.

"""
const true_positive_rate = TruePositiveRate()
const tpr = TruePositiveRate()
const TPR = TruePositiveRate
const truepositive_rate  = TPR()
const recall = TPR()
const Recall = TPR
const sensitivity  = recall
const hit_rate     = recall

"""
    true_negative_rate

$(docstring(TrueNegativeRate()))

    true_negative_rate(ŷ, y)

True negative rate for observations `ŷ` and ground truth `y`. Assigns
`false` to first element of `levels(y)`. To reverse roles, use
`TrueNegativeRate(rev=true)` instead of `true_negative_rate`.

For more information, run `info(true_negative_rate)`.

"""
const true_negative_rate = TrueNegativeRate()
const tnr = TrueNegativeRate()
const TNR = TrueNegativeRate
const truenegative_rate  = TNR()
const Specificity = TNR
const specificity  = truenegative_rate
const selectivity  = specificity

"""
    false_positive_rate

$(docstring(FalsePositiveRate()))

    false_positive_rate(ŷ, y)

False positive rate for observations `ŷ` and ground truth `y`. Assigns
`false` to first element of `levels(y)`. To reverse roles, use
`FalsePositiveRate(rev=true)` instead of `false_positive_rate`.

For more information, run `info(false_positive_rate)`.

"""
const false_positive_rate = FalsePositiveRate()
const fpr = FalsePositiveRate()
const FPR = FalsePositiveRate
const falsepositive_rate = FPR()
const fallout      = falsepositive_rate

"""
    false_negative_rate

$(docstring(FalseNegativeRate()))

    false_negative_rate(ŷ, y)

False negative rate for observations `ŷ` and ground truth `y`. Assigns
`false` to first element of `levels(y)`. To reverse roles, use
`FalseNegativeRate(rev=true)` instead of `false_negative_rate`.

For more information, run `info(false_negative_rate)`.

"""
const false_negative_rate = FalseNegativeRate()
const fnr = FalseNegativeRate()
const FNR = FalseNegativeRate
const falsenegative_rate = FNR()

"""
    false_discovery_rate

$(docstring(FalseDiscoveryRate()))

    false_discovery_rate(ŷ, y)

False discovery rate for observations `ŷ` and ground truth `y`. Assigns
`false` to first element of `levels(y)`. To reverse roles, use
`FalseDiscoveryRate(rev=true)` instead of `false_discovery_rate`.

For more information, run `info(false_discovery_rate)`.

"""
const false_discovery_rate = FalseDiscoveryRate()
const fdr = FalseDiscoveryRate()
const FDR = FalseDiscoveryRate
const falsediscovery_rate = FDR()
const miss_rate    = falsenegative_rate

"""
    negative_predictive_value

$(docstring(NPV()))

    negative_predictive_value(ŷ, y)

Negative predictive value for observations `ŷ` and ground truth
`y`. Assigns `false` to first element of `levels(y)`. To reverse roles,
use `NPV(rev=true)` instead of `negative_predictive_value`.

For more information, run `info(negative_predictive_value)`.

"""
const negative_predictive_value = NPV()
const npv = NPV()
const negativepredictive_value = NPV()

"""
    positive_predictive_value

$(docstring(Precision()))

    positive_predictive_value(ŷ, y)

Positive predictive value for observations `ŷ` and ground truth
`y`. Assigns `false` to first element of `levels(y)`. To reverse roles,
use `Precision(rev=true)` instead of `positive_predictive_value`.

For more information, run `info(positive_predictive_value)`.

"""
const positive_predictive_value = Precision()
const ppv = Precision()
const positivepredictive_value = Precision()
const PPV = Precision

## INTERNAL FUNCTIONS ON CONFUSION MATRIX

_tp(m::CM2) = m[2,2]
_tn(m::CM2) = m[1,1]
_fp(m::CM2) = m[2,1]
_fn(m::CM2) = m[1,2]

_tpr(m::CM2) = tp(m) / (tp(m) + fn(m))
_tnr(m::CM2) = tn(m) / (tn(m) + fp(m))
_fpr(m::CM2) = 1 - _tnr(m)
_fnr(m::CM2) = 1 - _tpr(m)

_fdr(m::CM2) = fp(m) / (tp(m) + fp(m))
_npv(m::CM2) = tn(m) / (tn(m) + fn(m))

## Callables on CM2
# NOTE: here we assume the CM was constructed a priori with the
# proper ordering so the field `rev` in the measure is ignored

(::TruePositive)(m::CM2)  = _tp(m)
(::TrueNegative)(m::CM2)  = _tn(m)
(::FalsePositive)(m::CM2) = _fp(m)
(::FalseNegative)(m::CM2) = _fn(m)

(::TPR)(m::CM2) = _tpr(m)
(::TNR)(m::CM2) = _tnr(m)
(::FPR)(m::CM2) = _fpr(m)
(::FNR)(m::CM2) = _fnr(m)

(::FDR)(m::CM2) = _fdr(m)
(::NPV)(m::CM2) = _npv(m)

(::Precision)(m::CM2) = 1.0 - _fdr(m)

function (::FScore{β})(m::CM2) where β
    β2   = β^2
    prec = precision(m)
    rec  = recall(m)
    return (1 + β2) * (prec * rec) / (β2 * prec + rec)
end

## Callables on vectors

for M in (TruePositive, TrueNegative, FalsePositive, FalseNegative,
          TPR, TNR, FPR, FNR,
          FDR, Precision, NPV, FScore)
    (m::M)(ŷ, y) = confmat(ŷ, y; rev=m.rev) |> m
end

# specify this as `precision` is in Base and so is ambiguous
Base.precision(m::CM2) = m |> Precision()
Base.precision(ŷ, y)   = confmat(ŷ, y) |> Precision()

## ROC COMPUTATION

"""
    _idx_unique_sorted(v)

Internal function to return the index of unique elements in `v` under the
assumption that the vector `v` is sorted in decreasing order.
"""
function _idx_unique_sorted(v::Vec{<:Real})
    n    = length(v)
    idx  = ones(Int, n)
    p, h = 1, 1
    cur  = v[1]
    @inbounds while h < n
        h     += 1                  # head position
        cand   = v[h]               # candidate value
        cand   < cur || continue    # is it new? otherwise skip
        p     += 1                  # if new store it
        idx[p] = h
        cur    = cand               # and update the last seen value
    end
    p < n && deleteat!(idx, p+1:n)
    return idx
end

"""
    tprs, fprs, ts = roc_curve(ŷ, y) = roc(ŷ, y)

Return the ROC curve for a two-class probabilistic prediction `ŷ` given the
ground  truth `y`. The true positive rates, false positive rates over a range
of thresholds `ts` are returned. Note that if there are `k` unique scores,
there are correspondingly  `k` thresholds and `k+1` "bins" over which the FPR
and TPR are constant:

* `[0.0 - thresh[1]]`
* `[thresh[1] - thresh[2]]`
* ...
* `[thresh[k] - 1]`

consequently, `tprs` and `fprs` are of length `k+1` if `ts` is of length `k`.

To draw the curve using your favorite plotting backend, do `plot(fprs, tprs)`.
"""
function roc_curve(ŷ::Vec{<:UnivariateFinite},
                   y::Vec{<:CategoricalElement})

    n       = length(y)
    lab_pos = levels(y)[2]
    scores  = pdf.(ŷ, lab_pos)
    ranking = sortperm(scores, rev=true)

    scores_sort = scores[ranking]
    y_sort_bin  = (y[ranking] .== lab_pos)

    idx_unique = _idx_unique_sorted(scores_sort)
    thresholds = scores_sort[idx_unique]

    # detailed computations with example:
    # y = [  1   0   0   1   0   0   1]
    # s = [0.5 0.5 0.2 0.2 0.1 0.1 0.1] thresh are 0.5 0.2, 0.1 // idx [1, 3, 5]
    # ŷ = [  0   0   0   0   0   0   0] (0.5 - 1.0] # no pos pred
    # ŷ = [  1   1   0   0   0   0   0] (0.2 - 0.5] # 2 pos pred
    # ŷ = [  1   1   1   1   0   0   0] (0.1 - 0.2] # 4 pos pred
    # ŷ = [  1   1   1   1   1   1   1] [0.0 - 0.1] # all pos pre

    idx_unique_2 = idx_unique[2:end]   # [3, 5]
    n_ŷ_pos      = idx_unique_2 .- 1   # [2, 4] implicit [0, 2, 4, 7]

    cs   = cumsum(y_sort_bin)          # [1, 1, 1, 2, 2, 2, 3]
    n_tp = cs[n_ŷ_pos]                 # [1, 2] implicit [0, 1, 2, 3]
    n_fp = n_ŷ_pos .- n_tp             # [1, 2] implicit [0, 1, 2, 4]

    # add end points
    P = sum(y_sort_bin) # total number of true positives
    N = n - P           # total number of true negatives

    n_tp = [0, n_tp..., P] # [0, 1, 2, 3]
    n_fp = [0, n_fp..., N] # [0, 1, 2, 4]

    tprs = n_tp ./ P  # [0/3, 1/3, 2/3, 1]
    fprs = n_fp ./ N  # [0/4, 1/4, 2/4, 1]

    return fprs, tprs, thresholds
end

const roc = roc_curve

ext: jl


text:
using BitTwiddle
using BitTwiddle: ConstructorPlan, ADTPlan, emit
using MacroTools


text:
# Works for Nosc={1, 2, 3}
"""
    plt = plotunitary(us, params, guardlev)

Plot the evolution of the state vector.

# Arguments
- `us:: Array{Complex{Float64},3})`: State vector history for each timestep
- `param:: objparams`: Struct with problem definition
- `guardlev:: Array{Bool,1})`: Boolean array indicating if a certain level is a guard level
"""
function plotunitary(us, params, guardlev)
    nsteps = length(us[1,1,:])
    Ntot = length(us[:,1,1])
    N =  length(us[1,:,1])
    t = range(0, stop = params.T, length = nsteps)
    if nsteps < 10000
        stride = 1
    elseif nsteps < 20000
        stride = 2
    else
        stride = 4 # only plot every 4th data point
    end
    rg = 1:stride:nsteps # range object

  # one figure for the response of each basis vector
    plotarray = Array{Plots.Plot}(undef, N) #empty array for separate plots

    for ii in 1:N # N = Ne[1] * Ne[2], ordered as
        if params.Nosc == 1
            statestr = string(ii-1)
        elseif params.Nosc == 2
            # Example: Ne[2] = 3, Ne[1]=4
            #   0,   1,  2,   3,   4,  5,   6,   7,   8,  9,  10, 11  (ii-1)
            # 00, 01, 02, 03, 10, 11, 12, 13, 20, 21, 22, 23
            statestr = string( div((ii-1), params.Ne[1]), mod(ii-1, params.Ne[1]))
        elseif params.Nosc == 3
            # Example: Ne[3] = 2, Ne[2] = 3, Ne[1]=4
            #     0,    1,     2,    3,    4,     5,     6,    7,     8,    9,   10,   11,   12,   13,   14,  15,   16,   17,   18,  19,   20,   21,   22,   23  (ii-1)
            # 000, 001, 002, 003, 010, 011, 012, 013, 020, 021, 022, 023,  100, 101, 102, 103, 110, 111, 112, 113, 120, 121, 122, 123
            s3 = div((ii-1), params.Ne[1]*params.Ne[2])
            s12 = (ii-1) % (params.Ne[1]*params.Ne[2])
            s2 = div(s12, params.Ne[1])
            s1 = s12 %  params.Ne[1]
            statestr = string( s3, s2, s1 )
        end
        titlestr = raw"Evolution from state $|" * statestr * raw"\rangle$"
#        h = plot(title = titlestr, size=(650, 400), legend= :outerright)
        h = plot(title = titlestr, legend= :outerright)
        for jj in 1:Ntot
            # Is jj an essential level?
            if !guardlev[jj]
                if params.Nosc == 1
                    labstr = string(jj-1) # "" for no labels
                elseif params.Nosc == 2
                    #   for s2 = 0:Nt[2]-1
                    #     for s1 = 0:Nt[1]-1
                    #        jj = s2*Nt[1] + s1 + 1
                    labstr = string( div((jj-1),params.Nt[1]), mod(jj-1, params.Nt[1]) )
                elseif params.Nosc == 3
                    # for s3 = 0:Nt[3]-1
                    #   for s2 = 0:Nt[2]-1
                    #     for s1 = 0:Nt[1]-1
                    #        jj = s3*Nt[1]*Nt[2] + s2*Nt[1] + s1 + 1
                    s3 = div((jj-1), params.Nt[1]*params.Nt[2])
                    s12 = jj -1 - s3 * params.Nt[1]*params.Nt[2]
                    s2 = div(s12, params.Nt[1])
                    s1 = s12 - s2 * params.Nt[1]
                    labstr = string( s3, s2, s1 )
                end
                plot!(t[rg], abs.(us[jj,ii,rg]).^2, lab = labstr, xlabel = "Time [ns]", ylabel="Population")
            end
        end
        plotarray[ii] = h
    end
    if N <= 2
        plt = plot(plotarray..., layout = (N,1))
    else
        plt = plot(plotarray..., layout = N, size=(650, 400))
    end

    return plt
end

# Works for Nosc={1, 2, 3}
"""
    plt = plotspecified(us, params, guardlev, specifiedlev)

Plot the evolution of the state vector for specified levels.

# Arguments
- `us:: Array{Complex{Float64},3})`: State vector history for each timestep
- `param:: objparams`: Struct with problem definition
- `us:: Array{Bool,1})`: Boolean array indicating if a certain level is a guard level
- `specifiedlev:: Array{Bool,1}`: Boolean array indicating which levels to be plotted
"""
function plotspecified(us, params, guardlev::Array{Bool,1}, specifiedlev::Array{Bool,1})
    nsteps = length(us[1,1,:])
    Ntot = length(us[:,1,1])
    N =  length(us[1,:,1])
    t = range(0, stop = params.T, length = nsteps)
    if nsteps < 10000
        stride = 1
    elseif nsteps < 20000
        stride = 2
    else
        stride = 4 # only plot every 4th data point
    end
    rg = 1:stride:nsteps # range object

    # How many specified levels are there?
    nForb = 0
    gLev=-1
    for jj in 1:Ntot
        if specifiedlev[jj]
            nForb += 1
            gLev = jj-1
        end
    end

    # plot rows corresponding to specified guard levels = specified levels
#    plotarray = Array{Plots.Plot}(undef, N) #empty array for separate plots

    if nForb == 1
        titlestr = raw"Population of state $|" * string(gLev) * raw"\rangle$"
    else
        titlestr = "Population of guard levels"
    end
    h = plot(title = titlestr, size=(700, 350), legend= :outerright) # make it big to fit the legend
#    h = plot(title = titlestr, size=(700, 350)) # put legend inside plot

    iplot = 0
    for col in 1:Ntot
        if !guardlev[col] # col is an essential level: plot this column in the 'us' array
            iplot += 1
            if params.Nosc == 1
                statestr = string(col-1)
            elseif params.Nosc == 2
                s2 = div((col-1), params.Nt[1])
                s1 = (col-1) % params.Nt[1]
                statestr = string( s2, s1 )
            elseif params.Nosc == 3
                s3 = div((col-1), params.Nt[1]*params.Nt[2])
                s12 = (col-1) % (params.Nt[1]*params.Nt[2])
                s2 = div(s12, params.Nt[1])
                s1 = s12 %  params.Nt[1]
                statestr = string( s3, s2, s1 )
            end
            # tmp
            # println("col=", col, " guardlev=", guardlev[col], " statestr= ", statestr)
            # end tmp
            for row in 1:Ntot
                # Is row a specified level?
                if specifiedlev[row] # only plot the specified = forbidden levels (rows in 'us')
                    if params.Nosc == 1
                        labstr = string(row-1)
                    elseif params.Nosc == 2
                        s2 = div((row-1), params.Nt[1])
                        s1 = (row-1) % params.Nt[1]
                        labstr = string( s2, s1 )
                    elseif params.Nosc == 3
                        s3 = div((row-1), params.Nt[1]*params.Nt[2])
                        s12 = row -1 - s3 * params.Nt[1]*params.Nt[2]
                        s2 = div(s12, params.Nt[1])
                        s1 = s12 - s2 * params.Nt[1]
                        labstr = string( s3, s2, s1 )
                    end
                    # Not enough room for too many labels
                    if nForb > 16
                        labstr = ""
                    end
                    plot!(t[rg], abs.(us[row,iplot,rg]).^2, lab = labstr * " from " * statestr * " state", xlabel = "Time [ns]", ylabel="Population")
                end
            end # for row
#            plotarray[iplot] = h
        end
    end
    plt3 = h
    # if N <= 2
    #     plt3 = plot(plotarray..., layout = (N,1))
    # else
    #     plt3 = plot(plotarray..., layout = N)
    # end

    return plt3
end

function plot_forward(us, T)
    nsteps = length(us[1,1,:])
    Ntot = length(us[:,1,1])
    N =  length(us[1,:,1])

    U0 = us[:,1,1]
    dum, ii = findmax(abs.(U0)) # Which basis did we start from?
    # test
    println("plot_forward: initial condition corresponds to element # ", ii-1, " (vectorized representation)")

    t = range(0, stop = T, length = nsteps)
    rg = 1:nsteps # range object  # plot all data points

#    ii = 1;
    titlestr = raw"Evolution from state $|" * string(ii-1) * raw"\rangle$"
    plt = plot(xlabel = "Time [ns]", ylabel = "Population", title = titlestr, size=(750, 400), legend=:outerright) # abs^2
    for jj = 1:Ntot
#        if jj != ii # don't plot the evolution of the initial state
            labstr = string("State ", jj-1)
            plot!(t[rg], abs.(us[jj,1,:]).^2, lab=labstr) # abs^2
#        end
    end
    return plt
end

# Evaluate the control functions on a grid in time in units of GHz
"""
    pj [, qj] = evalctrl(params, pcof0, td, func)

Evaluate the control function with index `func` at an array of time levels `td`.

NOTE: the control function index `func` is 1-based.

NOTE: The return value(s) depend on `func`. For `func∈[1,Ncoupled]`, `pj, qj` are returned. Otherwise,
only `pj` is returned, corresponding to control number `func`.

# Arguments
- `params:: objparams`: Struct with problem definition
- `pcof0:: Array{Float64,1})`: Vector of parameter values
- `td:: Array{Float64,1})`: Time values control is to be evaluated
- `jHam:: Int64`: Index of the control signal desired
"""
function evalctrl(params::objparams, pcof0:: Array{Float64, 1}, td:: Array{Float64, 1}, jHam:: Int64)
    if params.pFidType == 3
        nCoeff = length(pcof0)-1
    else
        nCoeff = length(pcof0)
    end
    pcof = pcof0[1:nCoeff]
    D1 = div(nCoeff, (2*params.Ncoupled + params.Nunc)*params.Nfreq)  # number of B-spline coeff per control function

    if (params.use_bcarrier)
        # B-splines with carrier waves
        splinepar = bcparams(params.T, D1, params.Ncoupled, params.Nunc, params.Cfreq, pcof)
    else
        # regular B-splines
        splinepar = splinepar(params.T, D1, 2*params.Ncoupled + params.Nunc, pcof)
    end

    # define inline function to enable vectorization over t
    controlplot(t, splinefunc) = controlfunc(t, splinepar, splinefunc)

    fact = 1.0/(2*pi) # conversion factor to GHz
    fact = 1.0 # conversion factor to rad/ns

    if(jHam <= params.Ncoupled && params.Ncoupled > 0)
        # coupled controls
        qs = (jHam-1)*2
        qa = qs+1

        pj = fact.*controlplot.(td, qs)
        qj = fact.*controlplot.(td, qa)
        return pj, qj
    else
        qu = jHam
        pj = controlplot.(td, qu)
        return pj
    end

end

"""
    guardlev = identify_guard_levels(params[, custom = 0])

Build a Bool array indicating if a given energy level is a guard
level in the simulation.

# Arguments
- `params:: objparams`: Struct with problem definition
- `custom:: Int64`: A nonzero value gives a special stirap pulses case
"""
function identify_guard_levels(params::Juqbox.objparams, custom:: Int64 = 0)
    # identify all guard levels
    Ntot = params.N+params.Nguard
    guardlev = fill(false, Ntot)

    if params.Nosc == 1
        if custom == 0
            guardlev[params.N+1:Ntot] .= true
        else # special case for stirap pulses
            guardlev[2] = true
            guardlev[4] = true
        end
    elseif params.Nosc == 2
        for q2 in 1:params.Nt[2]
            for q1 in 1:params.Nt[1]
                if q1 > params.Ne[1] || q2 > params.Ne[2]
                    guardlev[(q2-1)*params.Nt[1] + q1] = true
                end
            end
        end
    elseif params.Nosc == 3
        for q3 in 1:params.Nt[3]
            for q2 in 1:params.Nt[2]
                for q1 in 1:params.Nt[1]
                    if q1 > params.Ne[1] || q2 > params.Ne[2] || q3 > params.Ne[3]
                        guardlev[(q3-1)*params.Nt[1]*params.Nt[2] + (q2-1)*params.Nt[1] + q1] = true
                    end
                end
            end
        end
    end
    return guardlev
end #identify_guard_levels

"""
    forbiddenlev = identify_guard_levels(params[, custom = 0])

Build a Bool array indicating which energy levels are forbidden levels in the state vector. The
forbidden levels in a state vector are defined as thos corresponding to the highest energy level in
at least one of its subsystems.

# Arguments
- `params:: objparams`: Struct with problem definition
- `custom:: Int64`: For nonzero value special stirap pulses case
"""
function identify_forbidden_levels(params:: Juqbox.objparams, custom::Int64 = 0)
    # identify all forbidden levels
    Ntot = params.N+params.Nguard
    forbiddenlev = fill(false, Ntot)

    if params.Nosc == 1
        if custom == 0
            forbiddenlev[Ntot] = true
        else # Special case for stirap pulses
            forbiddenlev[2] = true
            forbiddenlev[4] = true
        end
    elseif params.Nosc == 2
        k=0
        for q2 in 1:params.Nt[2]
            for q1 in 1:params.Nt[1]
                k += 1
                if q1 == params.Nt[1] || q2 == params.Nt[2]
                    forbiddenlev[k] = true
                end
            end
        end
    elseif params.Nosc == 3
        k=0
        for q3 in 1:params.Nt[3]
            for q2 in 1:params.Nt[2]
                for q1 in 1:params.Nt[1]
                    k += 1
                    if q1 == params.Nt[1] || q2 == params.Nt[2] || q3 == params.Nt[3]
                        forbiddenlev[k] = true
                    end
                end
            end
        end
    end
    return forbiddenlev
end #identify_forbidden_levels

function specify_level3(params:: Juqbox.objparams, Nl3:: Int64) # Nl3 is 0-based
    # identify all  levels
    Ntot = params.N+params.Nguard
    specifiedlev = fill(false, Ntot)

    if params.Nosc == 3
        k=0
        for q3 in 1:params.Nt[3]
            for q2 in 1:params.Nt[2]
                for q1 in 1:params.Nt[1]
                    k += 1
                    if q3 == Nl3+1
                        specifiedlev[k] = true
                    end
                end
            end
        end
    end
    return specifiedlev
end #specify_level3

"""
    marg_prob = marginalize3(params, unitaryhist)

Evaluate marginalized probabilities for the case of 3 subsystems.

# Arguments
- `param:: objparams`: Struct with problem definition
- `unitaryhist:: Array{Complex{Float64},3})`: State vector history for each timestep
"""
function marginalize3(params:: Juqbox.objparams, unitaryhist:: Array{Complex{Float64},3})
    nsteps1 = size(unitaryhist,3)
    if params.Nosc == 3
        marg_prob = zeros(params.Nt[3], params.N, nsteps1)

        k=0
        for q3 in 1:params.Nt[3]
            for q2 in 1:params.Nt[2]
                for q1 in 1:params.Nt[1]
                    #                    offset = (q3-1)*params.Nt[1]*params.Nt[2]
                    k += 1
                    marg_prob[q3, :, :] += abs.(unitaryhist[k, :, :]).^2
                end # for q1
            end # for q2
        end # for q3

        return marg_prob
    end # if Nosc = 3
end #marginalize3

"""
    pconv = plot_conv_hist(params [, convname:: String=""])

Plot the optimization convergence history, including history of
the different terms in the objective function and the norm of the gradient.

# Arguments
- `param:: objparams`: Struct with problem definition
- `convname:: String`: Name of plot file to be generated
"""
function plot_conv_hist(params:: Juqbox.objparams, convname:: String="")
    pconv = Plots.plot(xlabel="Iteration", title="Convergence history", size=(400, 300))
    if params.saveConvHist && length(params.objHist)>0
        nIter = length(params.objHist)
        # note that primaryHist and secondaryHist have one additional elements corresponding to iter=0

        if minimum(params.objHist) > 0 && minimum(params.primaryHist) > 0 &&
            (params.Nguard == 0 || minimum(params.secondaryHist) > 0)
            Plots.plot!(yscale=:log10)
        end
        Plots.plot!(1:nIter, params.objHist, lab=L"{\cal G}") #, markershape=:utriangle, markersize=5)
        Plots.plot!(1:nIter, params.primaryHist, lab=L"{\cal J}_1", style=:dash)
        if (params.Nguard > 0)
            Plots.plot!(1:nIter, params.secondaryHist, lab=L"{\cal J}_2")
        end
        Plots.plot!(1:nIter, params.dualInfidelityHist, lab=L"\|\nabla{\cal G} - z\|_\infty") # dual infeasibility
        Plots.xlims!((0, nIter+1))

        if length(convname)>0
            Plots.savefig(pconv, convname)
            println("Saved convergence history plot on file '", convname, "'")
        end
    else
        println("Warning: plot_conv_hist: no convergence history to plot")
    end
    return pconv
end #plot_conv_hist


text:
module test_RotationMatrix

import ModiaMath

#  Desired:
#    using Test
#    using LinearAlgebra
#    using StaticArrays
#    using Unitful
#
#  In order that these packages need not to be defined in the user environment, they are included via ModiaMath:
using ModiaMath.Test
using ModiaMath.LinearAlgebra
using ModiaMath.StaticArrays
using ModiaMath.Unitful

angle1 = pi / 2
angle2 = 90u"°"

R1a = ModiaMath.rot1(angle1)
R1b = ModiaMath.rot1(angle2)
R1c = ModiaMath.rot_e([1,0,0], angle1)
R1d = ModiaMath.rot_e([1,0,0], angle2)
R1e = ModiaMath.rot_nxy([1,0,0], [0,0,1])

R2a = ModiaMath.rot2(angle1)
R2b = ModiaMath.rot2(angle2)
R2c = ModiaMath.rot_e([0,1,0], angle1)
R2d = ModiaMath.rot_e([0,1,0], angle2)

R3a = ModiaMath.rot3(angle1)
R3b = ModiaMath.rot3(angle2)
R3c = ModiaMath.rot_e([0,0,1], angle1)
R3d = ModiaMath.rot_e([0,0,1], angle2)

R4a = ModiaMath.rot123(angle1, angle1, angle1)
R4b = ModiaMath.rot123(angle2, angle2, angle2)
R4c = ModiaMath.rot_e([1,0,0], -angle1) *
      ModiaMath.rot_e([0,1,0], -angle1) *
      ModiaMath.rot_e([0,0,1], -angle1)

R5  = ModiaMath.rot_nxy([1,0,0], [0,0,1])

v1  = SVector{3,Float64}(2.0, -1.0, 4.0)

v2a = ModiaMath.resolve2(R1a, v1)
v1a = ModiaMath.resolve1(R1a, v2a)

R6a = ModiaMath.absoluteRotation(R1a, R2a)
R6b = ModiaMath.rot123(angle1, angle1, 0.0)
R7  = ModiaMath.relativeRotation(R1a, R6b)

@testset "ModiaMath.Frames: test RotationMatrix" begin
    @test isapprox(R1b, R1a)
    @test isapprox(R1c, R1a)
    @test isapprox(R1d, R1a)
    @test isapprox(R1e, R1a)

    @test isapprox(R2b, R2a)
    @test isapprox(R2c, R2a)
    @test isapprox(R2d, R2a)

    @test isapprox(R3b, R3a)
    @test isapprox(R3c, R3a)
    @test isapprox(R3d, R3a)

    @test isapprox(R4b, R4a)
    @test isapprox(R4c, R4a)

    @test isapprox(R5, R1a)

    @test isapprox(v1, v1a)

    @test isapprox(R6a, R6b)
    @test isapprox(R7, R2a)

    R4d = ModiaMath.inverseRotation(R4a)
    R4e = ModiaMath.inverseRotation(R4d)
    @test isapprox(R4a, R4e)

    angle4 = 45u"°"
    e      = normalize([1.0, 1.0, 1.0])
    R      = ModiaMath.rot_e(e, angle4)
    v3     = [1.0, 2.0, 3.0]
    v4     = ModiaMath.resolve2(R, v3)
    angle5 = ModiaMath.planarRotationAngle(e, v3, v4)
    isapprox(angle4, angle5)

    isapprox(-ModiaMath.eAxis(3), ModiaMath.eAxis(-3))
end

end


text:
s = """
module test; if true; end
"""
msgs = lintstr(s)
@test msgs[1].code == :E112
@test occursin("incomplete:", msgs[1].message)

s = """
2 +
"""
msgs = lintstr(s)
@test msgs[1].code == :E112
@test occursin("incomplete:", msgs[1].message)


text:
# TODO:
"""
- load Contours.jl similar to DataFrames
- method to build grid from x/y/z vectors
- method to wrap contours creation
- method to plot contours as custom shapes (TODO: create Stroke and Fill types and add markerstroke/markerfill args)

"""

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

# immutable Vertex
#   x::Float64
#   y::Float64
#   z::Float64
# end

# immutable Edge
#   v::Vertex
#   u::Vertex
# end

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

# # one rectangle's z-values and the center vertex
# # z is ordered: topleft, topright, bottomright, bottomleft
# immutable GridRect
#   z::Vector{Float64}
#   center::Vertex
#   data::Vector{Vertex}
# end

# type Grid
#   xs::Vector{Float64}
#   ys::Vector{Float64}
#   rects::Matrix{GridRect}
# end

# function splitDataEvenly(v::AbstractVector{Float64}, n::Int)
#   vs = sort(v)

# end

# # the goal here is to create the vertical and horizontal partitions
# # which define the grid, so that the data is somewhat evenly split
# function bucketData(x, y, z)

# end

# function buildGrid(x, y, z)
#   # create
# end


text:
# Tolerance
function truncated_svd(A::AbstractMatrix{T}, truncation::Real) where T <: Number
    truncation = min(truncation, one(T))
    U, S, V = svd(A)
    r = vec(S .> truncation*maximum(S))
    U = U[:, r]
    S = S[r]
    V = V[:, r]
    return U, S, V
end

# Explicit rank
function truncated_svd(A::AbstractMatrix{T}, truncation::Int) where T <: Number
    U, S, V = svd(A)
    r = [((i <= truncation && S[i] > zero(T)) ? true : false) for i in 1:length(S)]
    U = U[:, r]
    S = S[r]
    V = V[:, r]
    return U, S, V
end

## Cost functions
# Pareto functions
G(opt::AbstractKoopmanAlgorithm) = g(f) = first(f)

# Evaluate F
function F(opt::AbstractKoopmanAlgorithm)
    f(x, A, y::AbstractArray) = [norm(y .- A*x, 2)] # explicit
    return f
end

# General method with inputs
function (x::AbstractKoopmanAlgorithm)(X::AbstractArray, Y::AbstractArray, U::AbstractArray, B::AbstractArray)
    K, _ = x(X, Y-B*U)
    return (K, B)
end

"""
$(TYPEDEF)

Approximates the Koopman operator `K` based on

```julia
K = Y / X
```

where `Y` and `X` are data matrices. Returns a  `Eigen` factorization of the operator.

# Fields
$(FIELDS)

# Signatures
$(SIGNATURES)
"""
mutable struct DMDPINV <: AbstractKoopmanAlgorithm end;

# Fast but more allocations
function (x::DMDPINV)(X::AbstractArray, Y::AbstractArray)
    K = Y / X
    return (eigen(K),[])
 end

# DMDC
function (x::DMDPINV)(X::AbstractArray, Y::AbstractArray, U::AbstractArray)
    nx, m = size(X)
    nu, m = size(U)

    K̃ = Y / [X;U]
    K = K̃[:, 1:nx]
    B = K̃[:, nx+1:end]

    return (eigen(K), B)
end

"""
$(TYPEDEF)

Approximates the Koopman operator `K` based on the singular value decomposition
of `X` such that:

```julia
K = Y*V*Σ*U'
```

where `Y` and `X = U*Σ*V'` are data matrices. The singular value decomposition is truncated via
the `truncation` parameter, which can either be an `Int` indiciating an index based truncation or a `Real`
indiciating a tolerance based truncation. Returns a `Eigen` factorization of the operator.

# Fields
$(FIELDS)

# Signatures
$(SIGNATURES)
"""
mutable struct DMDSVD{T} <: AbstractKoopmanAlgorithm where T <: Number
    """Indiciates the truncation"""
    truncation::T
end;

DMDSVD() = DMDSVD(0.0)

# Slower but fewer allocations
function (x::DMDSVD{T})(X::AbstractArray, Y::AbstractArray) where T <: Real
    U, S, V = truncated_svd(X, x.truncation)
    xone = one(eltype(X))
    # Computed the reduced operator
    Sinv = Diagonal(xone ./ S)
    B = Y*V*Sinv
    Ã = U'B
    # Compute the modes
    λ, ω = eigen(Ã)
    φ = B*ω
    return (Eigen(λ, φ),[])
end

# DMDc
function (x::DMDSVD{T})(X::AbstractArray, Y::AbstractArray, U::AbstractArray) where T <: Real
    nx, m = size(X)
    nu, m = size(U)
    # Input space svd
    Ũ, S̃, Ṽ = truncated_svd([X;U], x.truncation)
    # Output space svd
    Û, Ŝ, V̂ = svd(Y)

    # Split the svd
    U₁, U₂ = Ũ[1:nx,:], Ũ[nx+1:end,:]

    xone = one(eltype(X))
    # Computed the reduced operator
    C = Y*Ṽ*Diagonal(xone ./ S̃ ) # Common submatrix
    # We do not project onto a reduced subspace here.
    # This would mess up our initial conditions, since sometimes we have
    # x1->x2, x2->x1
    Ã = Û'C*U₁'Û
    B̃ = C*U₂'
    # Compute the modes
    λ, ω = eigen(Ã)
    φ = C*U₁'Û*ω
    return (Eigen(λ, φ), B̃)
end

"""
$(TYPEDEF)

Approximates the Koopman operator `K` based on the singular value decomposition
of `X` and `Y` such that:

```julia
K₁ = Y*V₁*Σ₁*U₁'
K₂ = X*V₂*Σ₂*U₂'
K = (K₁*inv(K₂))^(0.5)
```

where `Y = U₂*Σ₂*V₂'` and `X = U₁*Σ₁*V₁'` are data matrices. The singular value decomposition is truncated via
the `truncation` parameter, which can either be an `Int` indiciating an index based truncation or a `Real`
indiciating a tolerance based truncation. Returns a `Eigen` factorization of the operator.

# Fields
$(FIELDS)

# Signatures
$(SIGNATURES)
"""
mutable struct fbDMD{T} <: AbstractKoopmanAlgorithm where T <: Number
    """Indiciates the truncation"""
    truncation::T
end;
fbDMD() = fbDMD(0.0)
function (x::fbDMD{T})(X::AbstractArray, Y::AbstractArray) where T <: Real
    U₁, S₁, V₁ = truncated_svd(X, x.truncation)
    xone = one(eltype(X))
    # Computed the reduced operator
    Sinv₁ = Diagonal(xone ./ S₁)
    B₁ = Y*V₁*Sinv₁
    A₁ = U₁'B₁  # forward DMD

    U₂, S₂, V₂ = truncated_svd(Y, x.truncation)
    yone = one(eltype(Y))
    # Computed the reduced operator
    Sinv₂ = Diagonal(yone ./ S₂)
    B₂ = X*V₂*Sinv₂
    A₂ = U₂'B₂  # Backward DMD

    Ã = sqrt(A₁/A₂)
    # Compute the modes
    λ, ω = eigen(Ã)
    φ = B₁*ω
    return Eigen(λ, φ)
end
"""
$(TYPEDEF)

Approximates the Koopman operator `K` with the algorithm `alg` over the rank-reduced data
matrices `Xᵣ = X Qᵣ` and `Yᵣ = Y Qᵣ`, where `Qᵣ` originates from the singular value decomposition of
the joint data `Z = [X; Y]`. Based on [this paper](http://cwrowley.princeton.edu/papers/Hemati-2017a.pdf).

If `rtol` ∈ (0, 1) is given, the singular value decomposition is reduced to include only
entries bigger than `rtol*maximum(Σ)`. If `rtol` is an integer, the reduced SVD up to `rtol` is used
for computation.

# Fields
$(FIELDS)

# Signatures
$(SIGNATURES)
"""
mutable struct TOTALDMD{R, A} <: AbstractKoopmanAlgorithm where {R <: Number, A <: AbstractKoopmanAlgorithm}
    truncation::R
    alg::A
end

TOTALDMD() = TOTALDMD(0.0, DMDPINV())

function (x::TOTALDMD)(X::AbstractArray, Y::AbstractArray)
    _ , _, Q = truncated_svd([X; Y], x.truncation)
    return x.alg(X*Q, Y*Q)
end

function (x::TOTALDMD)(X::AbstractArray, Y::AbstractArray, U::AbstractArray)
    _ , _, Q = truncated_svd([X; Y], x.truncation)
    return x.alg(X*Q, Y*Q, U*Q)
end

function (x::TOTALDMD)(X::AbstractArray, Y::AbstractArray, U::AbstractArray, B::AbstractArray)
    _ , _, Q = truncated_svd([X; Y], x.truncation)
    K, _ = x.alg(X*Q, (Y-B*U)*Q)
    return (K, B)
end


text:
module CompensatedArithmetic

!isdefined(Float) && typealias Float AbstractFloat;

export csdDot, csdDet2x2, csdCross3D, csdHorner

using ErrorfreeArithmetic

include("arith.jl")
include("det.jl")
include("horner.jl")
include("dotproduct.jl")

end # module


text:
# Copyright 2017-2021 Nicolae Cindea
# See accompanying license file.

module SymFEL
__precompile__(false)

using SymPy

# polynomial variable

include("lagrange.jl")
include("hermite.jl")
include("assemble1d.jl")
include("assemble2d.jl")
include("assemble3d.jl")

export x
export y
export z
export h
export hx
export hy
export hz
export xa, xb

"""
    x = SymPy.symbols("x")

Symbolic variable x.

This variable is exported.
"""
x = SymPy.symbols("x")

"""
    y = SymPy.symbols("y")

Symbolic variable y.

This variable is exported.
"""
y = SymPy.symbols("y")

"""
    z = SymPy.symbols("z")

Symbolic variable z.

This variable is exported.
"""
z = SymPy.symbols("z")

"""
    h = SymPy.symbols("h")

Symbolic variable h.

This variable is exported.
"""
h = SymPy.symbols("h")

"""
    hx = SymPy.symbols("hx")

Symbolic variable h.

This variable is exported.
"""
hx = SymPy.symbols("hx")

"""
    hy = SymPy.symbols("hy")

Symbolic variable hy.

This variable is exported.
"""
hy = SymPy.symbols("hy")

"""
    hz = SymPy.symbols("hz")

Symbolic variable hz.

This variable is exported.
"""
hz = SymPy.symbols("hz")

"""
    xa, xb = SymPy.symbols("xa xb")

Symbolic vaiables xa and xb.

This variables are used for variable coefficients problems. This variables are exported.
"""
xa, xb = SymPy.symbols("xa xb")

"""
    get_em(deg1=1, deg2=1, der1=0, der2=0; fe1="Lagrange", fe2="Lagrange")

Gives 1d elementary matrices for Lagrange and Hermite elements.

The elementary matrices are computed for elements of length `h`.
"""
function get_em(deg1=1, deg2=1, der1=0, der2=0; fe1="Lagrange", fe2="Lagrange")
    global x
    global h
    if fe1 == "Lagrange"
        if der1 ∈ collect(0:deg1)
            p1 = get_lagrange_basis(deg1)
        else
            error("The third argument should verify der1 ∈ [0, …, deg1]...")
        end
    end
    if fe2 == "Lagrange"
        if der2 ∈ collect(0:deg2)
            p2 = get_lagrange_basis(deg2)
        else
            error("The fourth parameter should verify der2 ∈ [0, …, deg2]...")
        end
    end
    if fe1 == "Hermite"
        if der1 ∈ collect(0:floor(Int, deg1 / 2 + 1))
            p1 = get_hermite_basis(deg1)
        else
            error("The third parameter should verify der1 ∈ [0, …,int(deg1 / 2 + 1)]...")
        end
    end
    if fe2 == "Hermite"
        if der2 ∈ collect(0:floor(Int, deg2 / 2 + 1))
            p2 = get_hermite_basis(deg2)
        else
            error("The fourth parameter should verify der2 ∈ [0, …,int(deg2 / 2 + 1)]...")
        end
    end
    if (fe1 != "Lagrange" && fe1 != "Hermite") || (fe2 != "Lagrange" && fe2 != "Hermite")
        error("Only Lagrange and Hermite FEM are implemented.")
    end

    l1 = length(p1)
    l2 = length(p2)
    M = Array{SymPy.Sym}(undef, l1, l2)
    for i = 1:l1
        for j = 1:l2
            M[i, j] = simplify(integrate(diff(p1[i], x, der1) * diff(p2[j], x, der2), (x, 0, h)))
        end
    end
    M
end

"""
    get_square_em(Mx::Array{SymPy.Sym, 2},
                  My::Array{SymPy.Sym, 2},
                  nc::Tuple{Array{Int64,1},Array{Int64,1}},
                  nr::Tuple{Array{Int64,1},Array{Int64,1}})

Get elementary matrices for a squared or a rectangular element.

# Arguments
  * `Mx` : elementary matrix for the x variable.
  * `My` : elementary matrix for the y variable
  * `nc` : order of nodes for column
  * `nr` : order of nodes for row

"""
function get_square_em(Mx::Array{SymPy.Sym, 2},
                       My::Array{SymPy.Sym, 2},
                       nc::Tuple{Array{Int64,1},Array{Int64,1}},
                       nr::Tuple{Array{Int64,1},Array{Int64,1}})

    pxc = size(Mx, 2)
    pxr = size(Mx, 1)

    pyc = size(My, 2)
    pyr = size(My, 1)

    pc = pxc * pyc
    pr = pxr * pyr

    M = Array{SymPy.Sym}(undef, pr, pc)

    for i = 1:pr
        for j = 1:pc
            M[i,j] = Mx[nr[1][i], nc[1][j]] * My[nr[2][i], nc[2][j]]
        end
    end
    M
end

"""
    get_cube_em(Mx::Array{SymPy.Sym, 2},
                My::Array{SymPy.Sym, 2},
                Mz::Array{SymPy.Sym, 2},
                nc::Tuple{Array{Int64,1},Array{Int64,1},Array{Int64,1}},
                nr::Tuple{Array{Int64,1},Array{Int64,1}},Array{Int64,1})

Get elementary matrices for a cube or a rectangular hexahedron element.

# Arguments
  * `Mx` : elementary matrix for the x variable.
  * `My` : elementary matrix for the y variable
  * `Mz` : elementary matrix for the z variable
  * `nc` : order of nodes for column
  * `nr` : order of nodes for row

"""
function get_cube_em(Mx::Array{SymPy.Sym, 2},
                     My::Array{SymPy.Sym, 2},
                     Mz::Array{SymPy.Sym, 2},
                     nc::Tuple{Array{Int64,1},Array{Int64,1},Array{Int64,1}},
                     nr::Tuple{Array{Int64,1},Array{Int64,1},Array{Int64,1}})

    pxc = size(Mx, 2)
    pxr = size(Mx, 1)

    pyc = size(My, 2)
    pyr = size(My, 1)

    pzc = size(Mz, 2)
    pzr = size(Mz, 1)

    pc = pxc * pyc * pzc
    pr = pxr * pyr * pzr

    M = Array{SymPy.Sym}(undef, pr, pc)

    for i = 1:pr
        for j = 1:pc
            M[i,j] = Mx[nr[1][i], nc[1][j]] * My[nr[2][i], nc[2][j]] * Mz[nr[3][i], nc[3][j]]
        end
    end
    M
end

end # module


text:
module ImageIO

function readimage() end
function writeimage() end

include(joinpath("tiff", "tiff.jl"))
include(joinpath("jpeg", "jpeg.jl"))
include(joinpath("png", "png.jl"))

export ModTIFF, ModJPEG, ModPNG

end # module


text:
module HPack

struct DecodeError <: Exception
    message::AbstractString
end

# package code goes here
bytearr(a::Vector{UInt8}) = a
bytearr(cs::Base.CodeUnits{UInt8,String}) = convert(Vector{UInt8}, cs)
bytearr(s::String) = bytearr(codeunits(s))
include("table.jl")
include("huffman.jl")
include("encode.jl")
include("decode.jl")

export encode
export decode
export DecodeError

end # module


text:
using DelaunayTriangulations

@testset "Triangulation" begin
    vertices = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
                (0.0, 0.5), (0.5, 0.5), (1.0, 0.5),
                (0.0, 1.0), (0.5, 1.0), (1.0, 1.0)]
    triangulation = Triangulation(vertices)
    @testset "Triangulation Construction" begin
        add_triangle!(triangulation, (1, 2, 5))
        @test ((1, 2) => 5) in triangulation.edges
        @test ((2, 5) => 1) in triangulation.edges
        @test ((5, 1) => 2) in triangulation.edges
        @test length(triangulation.edges) == 3
        @test adjacent(triangulation, (1, 2)) == 5
        @test adjacent(triangulation, (2, 5)) == 1
        @test adjacent(triangulation, (5, 1)) == 2
        @test adjacent(triangulation, (1, 5)) == 0 # No triangle containing edge
    end
    @testset "Add Triangles" begin
        # If (u,v,w) is an existing triangle, triangles (u,v,x),
        # (x,v,w), and (u,x,w) are rejected for any x
        add_triangle!(triangulation, (1, 2, 6))
        add_triangle!(triangulation, (3, 2, 5))
        add_triangle!(triangulation, (1, 4, 5))
        @test ((1, 2) => 5) in triangulation.edges
        @test ((2, 5) => 1) in triangulation.edges
        @test ((5, 1) => 2) in triangulation.edges
        @test length(triangulation.edges) == 3
        add_triangle!(triangulation, (1, 5, 4))
        add_triangle!(triangulation, (2, 3, 6))
        add_triangle!(triangulation, (2, 6, 5))
        add_triangle!(triangulation, (4, 5, 7))
        add_triangle!(triangulation, (5, 8, 7))
        add_triangle!(triangulation, (5, 6, 8))
        add_triangle!(triangulation, (8, 6, 9))
        @test length(triangulation.edges) == 24
        @test adjacent_to_vertex(triangulation, 1) in Set([(2, 5), (5, 4)])
        @test adjacent_to_vertex(triangulation, 2) in Set([(5, 1), (6, 5)])
        @test adjacent_to_vertex(triangulation, 3) == (6, 2)
        @test adjacent_to_vertex(triangulation, 9) == (8, 6)
    end
    @testset "Geometry Predicates" begin
        push!(triangulation.vertices, (0.25, 0.0001))
        push!(triangulation.vertices, (0.25, 0.0))
        @testset "incircle" begin
            @test incircle(triangulation, 1, 2, 5, 10) > 0.0
            @test incircle(triangulation, 10, 1, 2) > 0.0
            @test incircle(triangulation, 10, 2, 1) < 0.0
            @test incircle(triangulation, 3, 1, 2) == 0.0
            @test incircle(triangulation, 11, 1, 2) == 1.0
        end
        @testset "orient" begin
            @test orient(triangulation, 3, 1, 2) == 0.0
            @test orient(triangulation, 5, 1, 2) > 0.0
            @test orient(triangulation, 4, 5, 2) < 0.0
        end
    end
    @testset "DeleteTriangles" begin
        delete_triangle!(triangulation, (5, 2, 1)) # Triangle not in triangulation, nothing should change
        @test length(triangulation.edges) == 24
        @test adjacent(triangulation, (5, 1)) == 2
        delete_triangle!(triangulation, (2, 5, 1))
        @test length(triangulation.edges) == 21
        @test adjacent(triangulation, (5, 1)) == 0
        delete_triangle!(triangulation, (5, 2, 6))
        @test adjacent_to_vertex(triangulation, 2) == (3, 6)
        delete_triangle!(triangulation, (3, 6, 2))
        @test adjacent_to_vertex(triangulation, 2) == (0, 0)
        delete_triangle!(triangulation, (5, 4, 1))
        delete_triangle!(triangulation, (9, 8, 6))
        delete_triangle!(triangulation, (8, 5, 6))
        delete_triangle!(triangulation, (4, 5, 7))
        delete_triangle!(triangulation, (8, 7, 5))
        @test length(triangulation.edges) == 0
    end
end

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
"""
    runsim(config, seed, prerun)

Performs a simulation run using configuration file `config`, random seed `seed`
and other settings provided via commandline, configuration file or the defaults.
"""
function runsim(config::String = "", seed::Integer = 0)
    initsettings(defaultSettings()) #needed for log calls during `getsettings()`
    initsettings(getsettings(config, seed))
    Random.seed!(setting("seed"))
    setupdatadir()
    world = Patch[]
    timesteps = 0
    timeoffset = 0
    correctmode!()
    for m in 1:length(setting("maps"))
        timeoffset += timesteps
        timesteps, maptable = readmapfile(setting("maps")[m])
        if m == 1
            world = createworld(maptable)
            writedata(world, timeoffset)
            setting("fasta") != "off" && writefasta(world, timeoffset)
        else
            world = updateworld(world, maptable)
        end
        simulate!(world, timesteps, timeoffset)
    end
    world
end

"""
    rungemm(config, seed)

Wrapper for `runsim()`
Runs a simulation using configuration file `config`, random seed `seed`
and other settings provided via commandline, configuration file or the defaults.
"""
function rungemm(config::String = "", seed::Integer = 0)
    @time world = runsim(config, seed)
end

"""
    correctmode!()

A heuristic function to keep backward compatibility. (The `mode` setting was
only introduced after several different experiments had already been carried
out with GeMM, so old config files don't include it.)
"""
function correctmode!()
    if setting("mode") == "default"
        # detect invasion mode
        if setting("global-species-pool") > 0 || setting("propagule-pressure") > 0
            simlog("Detected implicit invasion mode, updated settings.", 'w')
            updatesetting("mode", "invasion")
        end
        # XXX further modes?
    end
    if setting("precrange") != 0
        updatesetting("maxprec", setting("precrange"))
        simlog("`precrange` is deprecated, use `maxprec` instead.", 'w')
    end
end

let settings::Dict{String, Any}
    """
        initsettings(newsettings)

    Define a new settings dict.
    """
    global function initsettings(newsettings::Dict{String, Any})
        settings = newsettings
    end

    """
        setting(param)

    Return a configuration parameter from the settings.
    """
    global function setting(param::String)
        settings[param]
    end

    """
        updatesetting(param, value)

    Change the value of an individual config parameter.
    Use with caution!
    """
    global function updatesetting(param::String, value::Any)
        settings[param] = value
    end

    """
        settingkeys()

    Return all keys in the settings dict.
    """
    global function settingkeys()
        keys(settings)
    end
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
# Download latest binary shared library of CoolProp project
import JSON
const destpathbase = abspath(joinpath(@__FILE__,"..","..","lib"));
println(destpathbase);
const srctablepath = abspath(joinpath(@__FILE__,"..","..","src","Tables"));
println(srctablepath);
m1=0
m2=0
m3=0
function canconvert(ty::Type,matrix::Array,database,sizearray)
  global m1,m2,m3
  j=1
  si=Int(0)
  r::Tuple
  r=(NaN,)
  for(j in 1:size(matrix)[1])
    try
      r=ty(tuple(matrix[j,:]...))
    catch
      return matrix[j,:],false
    end
    if (sizearray!=nothing)
      m1>0 && (r2=rpad(r[2],m1,' '))
      m2>0 && (r3=rpad(r[3],m2,' '))
      m3>0 && (r4=rpad(r[4],m3,' '))
      m1>0 && (r=(r[1],r2,r3,r4,r[5],r[6]))
      si+=writetup(database,r)
    else
      b=[(r[i]==matrix[j,i] || isnan(r[i])) for i in 1:length(r)]
      if(sum(b)!=length(r))
        println(matrix[j,:])
        println(r)
        return false,matrix[j,:]
      end
    end
  end
  (sizearray!=nothing) && push!(sizearray,(si,j))
  return true
end
function writetup(database,t::Tuple)
  i=0
  i+=write(database,t...)
  return i;
end
function buildbinarydatabase(database,si)
  ret=true
  global m1,m2,m3
  propnames=["Criticals",
              "LiquidsDensities",
              "FormationEnergy",
              "Compounds",
              "Cp",
              "LiquidsCp",
              "VaporizHeat",
              "LiquidThermal",
              "VaporPressure",
              "LiquidViscos",
              "VaporThermal",
              "VaporViscos"]

  tablenames=["perryCriticals_Table2_141.table",
              "perryDensities_Table2_32.table",
              "perryEnergiesOfFormation_Table2_179.table",
              "perryFormulaCompounds.table",
              "perryHeatCapIdealGas_Tables156_155.table",
              "perryHeatCapLiquids_Table2_153.table",
              "perryHeatsofVaporization_Table2_150.table",
              "perryLiqidThermalConductivity_Table2_315.table",
              "perryLiquidsVaporPressure_Table2_8.table",
              "perryLiquidViscosity_Table2_313.table",
              "perryVaporThermalConductivity_Table2_314.table",
              "perryVaporViscosity_Table2_312.table"]

  rowtypes=[Tuple{[UInt16;[Float64 for i=1:6]]...},
            Tuple{[UInt16;[Float64 for i=1:10];UInt8]...},
            Tuple{[UInt16;[Float64 for i=1:5]]...},
            Tuple{[UInt16,ASCIIString,ASCIIString,ASCIIString,Float64,Float64]...},
            Tuple{[UInt16;[Float64 for i=1:10];UInt8]...},
            Tuple{[UInt16;[Float64 for i=1:10];UInt8]...},
            Tuple{[UInt16;[Float64 for i=1:10]]...},
            Tuple{[UInt16;[Float64 for i=1:10]]...},
            Tuple{[UInt16;[Float64 for i=1:10]]...},
            Tuple{[UInt16;[Float64 for i=1:10]]...},
            Tuple{[UInt16;[Float64 for i=1:9]]...},
            Tuple{[UInt16;[Float64 for i=1:9]]...}]

  for i=1:length(tablenames)
    tablename=tablenames[i]
    file=open(joinpath(srctablepath, tablename))
    matrix=readdlm(file,';',header=false)
    close(file)
    if tablename=="perryFormulaCompounds.table"
      m1=(maximum(Int[length(i) for i in matrix[:,2]]))
      m2=(maximum(Int[length(i) for i in matrix[:,3]]))
      m3=(maximum(Int[length(i) for i in matrix[:,4]]))
    else
      m1=0;m2=0;m3=0;
    end
    ty=rowtypes[i]
    if canconvert(ty,matrix,database,si)
      database==nothing && println("$tablename convertion simulated.")
    else
      println("$tablename convertion failed.")
      ret=false
    end
  end
  if si!=nothing
    databaseaddress=open(joinpath(destpathbase,"perryanalytic.addresses"),"w")
    println("Writing addresses .....")
    startaddress=0;
    for i=1:length(tablenames)
      s=si[i]
      println(databaseaddress,propnames[i],';',startaddress,';',s[2],';',round(Int,s[1]/s[2]),';',rowtypes[i])
      startaddress+=s[1]
    end
    close(databaseaddress)
  end
  return ret
end

try
  mkdir(destpathbase)
  const OS_ARCH_CoolProp = (WORD_SIZE == 64) ? "64bit" : "32bit__cdecl";
  const OS_ARCH_FreeSteam = (WORD_SIZE == 64) ? "win64" : "win32";
  const latestVersion_CoolProp = JSON.parse(readall(download("http://sourceforge.net/projects/coolprop/best_release.json")))["release"]["filename"][11:15];
  const latestVersion_FreeSteam = "2.1"
  println("CoolProp latestVersion = $latestVersion_CoolProp, by default I am going to install it...")

  @windows_only begin
      urlbase = "http://netix.dl.sourceforge.net/project/coolprop/CoolProp/$latestVersion_CoolProp/shared_library/Windows/$OS_ARCH_CoolProp/"
      download(joinpath(urlbase,"CoolProp.dll"),joinpath(destpathbase,"CoolProp.dll"))
      download(joinpath(urlbase,"CoolProp.lib"),joinpath(destpathbase,"CoolProp.lib"))
      download(joinpath(urlbase,"exports.txt"),joinpath(destpathbase,"exports.txt"))
      println("downloaded => lib/CoolProp.dll")
      # FreeSteam
      println("I'm Getting FreeSteam Binaries...")
      urlbase = "http://cdn.rawgit.com/DANA-Laboratory/FreeSteamBinary/master/$latestVersion_FreeSteam/$OS_ARCH_FreeSteam/"
      download(joinpath(urlbase,"freesteam.dll"),joinpath(destpathbase,"freesteam.dll"))
      println("downloaded => lib/freesteam.dll")
  end
  @linux_only begin
      # CoolProp
      urlbase = "http://netix.dl.sourceforge.net/project/coolprop/CoolProp/$latestVersion_CoolProp/shared_library/Linux/64bit/libCoolProp.so.$latestVersion_CoolProp"
      download(urlbase,joinpath(destpathbase,"libCoolProp.so"))
      println("downloaded => lib/libCoolProp.so")
      # FreeSteam
      println("I'm Getting FreeSteam Binaries...")
      urlbase = "http://cdn.rawgit.com/DANA-Laboratory/FreeSteamBinary/master/"
      download(joinpath(urlbase,"$latestVersion_FreeSteam","linux","libfreesteam.so.1.0"),joinpath(destpathbase,"libfreesteam.so"))
      println("downloaded => lib/libfreesteam.so")
      download(joinpath(urlbase,"libgsl.so.0"),joinpath(destpathbase,"libgsl.so"))
      println("downloaded => lib/libgsl.so")
      download(joinpath(urlbase,"libgslcblas.so.0"),joinpath(destpathbase,"libgslcblas.so"))
      println("downloaded => lib/libgslcblas.so")
  end
catch err
  println("$err\n => If lib folder exists, and you want to refresh existing files, remove lib folder and retry")
end
if buildbinarydatabase(nothing, nothing) # Simulate convertion
  database=open(joinpath(destpathbase,"perryanalytic.binary"),"w")
  si=Vector{Tuple{Int,Int}}()
  println("********** Simulation done, I'm writing perryanalytic.binary **********")
  buildbinarydatabase(database,si)
  close(database)
end


text:
using SuperTB, PyPlot, ProgressMeter, LinearAlgebra
close("all")

# set up a model (amorphous system with AII-class hopping)
tb = TBModel(hopAgarwala)
tb.periodicity = trues(2) # periodic boundary conditions
tb.basisdim    = 2

# run over this parameter range
M = range(-4.5, 4.5, length=50)
sites = 12:12:24^2
Nsamples = 1 # number of samples in ensemble

β = zeros(length(M), length(sites), Nsamples)

@showprogress 2 "Building phase diagram (#sites-progress)..." for s = eachindex(sites)
    for g = 1:Nsamples
        tb.pos = rand(sites[s],2)

        for m = eachindex(M)
            tb.hopfun(ri,rj,per) = hopAgarwala(ri, rj, per,
                                               1.0/6,  # R
                                               M[m],   # M
                                               0.25,   # t₂
                                               0.5,    # λ
                                               1.0/24) # a
            # solve the system
            H = Hermitian(Matrix(setHopping!(tb)))
            F = eigen(H)

            # compute Bott indices
            β[m,s,g] = bott(F.vectors,tb,sites[s])[1]
        end
    end
end
βᵉ = dropdims(sum(β,3),dims=3)./Nsamples
ρ = sites./24^2

figure()
surf(ρ,M,βᵉ)

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
include("euler/euler.jl")
using .Calculus: len_collatz_chain
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result, max_chain, hashmap = 0, 0, Dict(1 => 1)
    for i ∈ 2:n-1
        if len_collatz_chain(i, hashmap) > max_chain result, max_chain = i, hashmap[i] end
    end
    return result
end

compute(1_000_000)

@benchmark compute(1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using GraphMatFun, LsqFit, LinearAlgebra
include("exp_reset_all.jl");
#include("simulationtools.jl");
include("simulationtools.jl");

f=exp
m=4
target_n=200 # Measure error wrt this discr
n=50;        # Optimize wrt this discr
rho=6.9e-1;  # SID rho m=4

its=5;
droptol0=1e-10;
opt_kwargs=Dict(:logger=> 0,:γ0 => 0.5,:droptol => droptol0,
                :linlsqr => :real_svd,:maxit => its)
base = init_state_mult(f,rho,n,m; eltype=Complex{BigFloat})
base.params[:n]=n;
base.params[:target_n]=target_n;
base.params[:opt_kwargs]=opt_kwargs;

# First the optimization free
sid_org=init_state_mult!(deepcopy(base),:sid,m,showmeta=true)

bbc_org=init_state_mult!(deepcopy(base),:bbc,m,showmeta=true)

sastre_org=init_state_mult!(deepcopy(base),:sastre,m,showmeta=true)

ps_org=init_state_mult!(deepcopy(base),:ps,m,showmeta=true)

mono_org=init_state_mult!(deepcopy(base),:mono,m,showmeta=true)

(mono,simlist,commandlist)=
     run_sequence(mono_org,"sGGGssddsddsddsddsggggggssssssGGGssssssssddsddsdsgggggggsdssssssssssssssssGGGsGGsGGsssssssdddssssssssssssszq");

(ps,simlist,commandlist)=
     run_sequence(ps_org,
     "sssssssssssddsssddssssdssGgggssdsssssssssssssDdgGGssssddssssssssssssssssddsssssssssssddssssssssssssddsssssssddsssssssssssssddgggssssssssssGssGGGssssssgggsdsssssssssssssssssssssssssddsssssssssssssssssssssssssssssssssddsssssssssssssssssssssssddssssddssssssddssddDssdggggsssssGssGGGGGGGsGsssssdsgggggggsddggsssssGGGGGGGsssssssssssssdsdszq")

(sastre,simlist,commandlist)=
     run_sequence(sastre_org,
                  "sssssddssssddsssssssssssdddssdddssssdddssdddsssNsssdddsssssssszq");

(sid,simlist,commandlist)=
     run_sequence(sid_org,
     "sssssddssssddsssssssssssdddssdddssssdddssdddsssNsssdddsssssssssdddssssddsdsdssddsddsddsddsddsddsddsddsNssdszq")

(bbc,simlist,commandlist)=
     run_sequence(bbc_org,
     "sssSSddddddddddddSSSSSSGGGSSddggggggggsssSSSSSSSSSGGGGGGGGSSzq")

include("exp_print_all.jl");
println("Run include(exp_save_all.jl) if you want to save");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using Catalyst, Graphviz_jll

rn = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
    S^2, R --> 0
end α β

# check can make a graph
gr = Graph(rn)

# check can save a graph
fname = Base.Filesystem.tempname()
savegraph(gr, fname, "svg")


text:
abstract type AbstractCurveSolution end

"""
    mutable struct MDCSolution{S, F} <: AbstractCurve
holds solution for minimally disruptive curve.
fields: sol, N, cost
sol is a conventional ODESolution structure
"""
mutable struct MDCSolution{S,F} <: AbstractCurveSolution
    sol::S
    N::Int64
    cost::F
end

"""
    function MDCSolution(sol, costf=nothing)
returns an mdc::MDCSolution out of a sol::ODESolution
"""
function MDCSolution(sol, costf=nothing)
    N = length(sol.prob.u0) ÷ 2
    return MDCSolution(sol, N, costf)
end

"""
    function (mdc::MDCSolution)(t::N) where N <: Number
returns array view of states and costates of curve at distance t
"""
function (mdc::MDCSolution)(t::N) where N <: Number
    return (states = (@view mdc.sol(t)[1:mdc.N]), costates = (@view mdc.sol(t)[mdc.N + 1:end]))
end

"""
    function (mdc::MDCSolution)(ts::A) where A <: Array
returns array of states and costates of curve at distances t
"""
function (mdc::MDCSolution)(ts::A) where A <: Array
    states_ = Array{Float64,2}(undef, mdc.N, length(ts))
    costates_ = Array{Float64,2}(undef, mdc.N, length(ts))
    for (i, el) in enumerate(ts)
        states_[:,i] = mdc(el)[:states]
        costates_[:,i] = mdc(el)[:costates]
    end
    return (states = states_, costates = costates_)
end

trajectory(mdc::MDCSolution) = Array(mdc.sol)[1:mdc.N, :]
trajectory(mdc::MDCSolution, ts) = mdc(ts)[:states]

costate_trajectory(mdc::MDCSolution) = Array(mdc.sol)[mdc.N + 1:end,:]
costate_trajectory(mdc::MDCSolution, ts) = mdc(ts)[:costates]

distances(mdc::MDCSolution) = mdc.sol.t
Δ(mdc::MDCSolution) = trajectory(mdc) .- mdc(0.)[:states]
Δ(mdc::MDCSolution, ts) = mdc(ts)[:states] .- mdc(0.)[:states]

"""
    cost_trajectory(mdc::MDCSolution, ts)
calculates cost on mdc curve at each point in the array/range ts
"""
function cost_trajectory(mdc::MDCSolution, ts)
    if mdc.cost === nothing
        @warn "MDCSolution struct has no cost function. You need to run mdc = add_cost(mdc, cost)"
        return
    else
        return [mdc.cost(el) for el in eachcol(mdc(ts)[:states])]
    end
end


text:
################
### old code ###
################

# this works but we will not use it
# instead use omp! as OMP function
# then omp() will compute entire path of betas
function omp_working_old{T <: Float}(
    Φ :: DenseMatrix{T},
    v :: DenseVector{T},
    m :: Int;
    quiet::Bool = true
)
    n,p = size(Φ);
    β = zeros(T, p);
    r = copy(v);
    Λ = Int[];
    a = zeros(T, n)
    x = [zero(T)]
    dot_p = zeros(T, p)

    for iter = 1:m
        # dot_p = Φ' * r;
        BLAS.gemv!('T', one(T), Φ, r, zero(T), dot_p)

#        λ = indmax(abs(dot_p));
        λ = indmax_abs(dot_p);
#        Λ = [Λ i];
        push!(Λ, λ)
#        sort!(Λ)
        temp = Φ[:, Λ];

        # x = argmin_{y} ( norm(v - temp*y) )
#        x = pinv(temp) * v;
#        x = temp \ v
        x, = lsqr(temp, v)

        # a = temp*x
        BLAS.gemv!('N', 1.0, temp, x, 0.0, a)

        # r = v - a;
        copy!(r, v)
        BLAS.axpy!(-1.0, a, r)

        # output if desired
        quiet || @printf("Iter %d, residual: %f\n", iter, norm(r));
    end
    β[Λ] = x
    return β
end


text:
"""
    query_intervals(spdx::AbstractString, created::AbstractVector{<:Interval{ZonedDateTime}})

Return count of search results based on the license for each created interval.
"""
function query_intervals(created::Vector{Interval{DateTime,Closed,Open}})
    @unpack pat, spdx = PARALLELENABLER
    subsquery = join([ string("_$idx:search(query:\"is:public fork:false mirror:false archived:false license:$spdx created:",
                       format(created[idx].first, "yyyy-mm-ddTHH:MM:SS\\Z"),
                       "..",
                       format(created[idx].last, "yyyy-mm-ddTHH:MM:SS\\Z"),
                       "\",type:REPOSITORY){...A}") for idx in eachindex(created)]);
    query = "fragment A on SearchResultItemConnection{repositoryCount} query Search{$subsquery}" |>
        (obj -> replace(obj, r"\s+" => " ")) |>
        strip |>
        string;
    result = graphql(query)
    json = JSON3.read(result.Data)
    data = json.data
    DataFrame(created = created, count = getproperty.(values(data), :repositoryCount))
end
"""
    query_intervals(created::AbstractVector{<:AbstractVector{Interval{ZonedDateTime}}})::DataFrame

Returns a
"""
function query_intervals(created::Vector{Vector{Interval{DateTime,Closed,Open}}})
    output = DataFrame()
    graphqlremaining =
        ( fetch(@spawnat(
            w,
            (proc = w,
             remaining = GHOST.PARALLELENABLER.pat.limits.graphql.remaining,
             resetat = GHOST.PARALLELENABLER.pat.limits.graphql.reset)
            )) for w in workers()
        ) |>
        DataFrame
    maptovalidprocs = sort!(graphqlremaining, [order(2, rev = true), order(3)])[!,1][1:min(length(READY.x), length(created))] .- 1
    for w in maptovalidprocs
        READY.x[w] = remotecall(GHOST.query_intervals, w + 1, popfirst!(created))
    end
    while !isempty(created)
        # local w
        w = findfirst(isready, @view(READY.x[maptovalidprocs]))
        if isnothing(w)
            sleep(3)
        else
            append!(output, fetch(READY.x[maptovalidprocs][w]))
            READY.x[maptovalidprocs[w]] = remotecall(GHOST.query_intervals, maptovalidprocs[w] + 1, popfirst!(created))
        end
    end
    while any(!isready, READY.x)
        sleep(3)
    end
    for f in @view(READY.x[maptovalidprocs][isready.(READY.x[maptovalidprocs])])
        append!(output, fetch(f))
    end
    sort!(output)
end
"""
    cleanintervals(row)

Returns the input if the count is 1,000 records or fewer.
If there are more than a 1,000 it splits them based on the ratio of the count.
"""
function cleanintervals(row)
    created = row.created
    cnt = row.count
    if ismissing(cnt) || cnt ≤ 1_000
        DataFrame(row)
    else
        new_step = ceil((created.last - created.first) ÷ (cnt ÷ 1_000 + 1), Second)
        new_date_period = range(created.first, created.last, step = new_step)
        DataFrame([ (created = Interval{Closed, Open}(ds, de), count = missing) for (ds, de) in zip(new_date_period[1:end - 1], new_date_period[2:end]) ])
    end
end
"""
    prune(data)

Prune the intervals based on the created and count values.
"""
function prune(data)
    isempty(data) && return data
    created = data.created
    cnt = data.count
    output = similar(data, 0)
    running_cnt = 0
    date_start = first(created).first
    for (created, cnt) in zip(created, cnt)
        if (running_cnt + cnt) ≥ 1_000
            push!(output, (created = Interval{Closed, Open}(date_start, created.first), count = running_cnt))
            date_start = created.first
            running_cnt = cnt
        else
            running_cnt += cnt
        end
    end
    push!(output, (created = Interval{Closed, Open}(date_start, created[end].last), count = running_cnt))
    output
end
"""
    format_tsrange(obj::Interval{ZonedDateTime})

Return the Postgres compatible form.
"""
format_tsrange(obj::Interval{DateTime}) = replace(string(obj), " .." => ",")
"""
    queries(conn::Connection,
            spdx::AbstractString,
            schema::AbstractString = "gh_2007_\$(Dates.year(floor(now(), Year) - Day(1)))")

This will upload the queries to the database with:
- spdx::text NOT NULL
- created::tsrange NOT NULL
- count::smallint NOT NULL
- asof::time
- done::bool NOT NULL
"""
function queries(spdx::AbstractString)
    # spdx = "mit"
    # schema = "gh_2007_$(Dates.year(floor(now(), Year) - Day(1)))"
    @unpack conn, schema = PARALLELENABLER
    @everywhere GHOST.PARALLELENABLER.spdx = $spdx
    calendaryear = parse(Int, schema[end - 3:end])
    created = vcat(floor(GH_FIRST_REPO_TS, Day),
                   DateTime("2009-01-01"):Month(6):DateTime("2010-01-01"),
                   DateTime("2010-01-01"):Month(3):DateTime("2011-01-01"),
                   DateTime("2011-01-01"):Month(2):DateTime("2012-01-01"),
                   DateTime("2012-01-01"):Month(1):DateTime("2013-01-01"),
                   DateTime("2013-01-01"):Week(1):DateTime("2013-12-30"),
                   DateTime("2013-12-30"):Day(1):DateTime(floor(now(utc_tz), Year), UTC)) |>
        unique
    created = [ Interval{Closed, Open}(start, stop) for (start, stop) in zip(@view(created[1:end - 1]), @view(created[2:end])) ]
    created = [ created[start:stop] for (start, stop) in zip(1:185:length(created), vcat((0:185:length(created))[2:end], length(created))) ]
    data = GHOST.query_intervals(created)
    data = GHOST.prune(data)
    data = reduce(vcat, cleanintervals(row) for row in eachrow(data))
    toreplace = ismissing.(data.count)
    created = data.created[toreplace]
    while !isempty(created)
        created = [ created[start:stop] for (start, stop) in zip(1:185:length(created), vcat((0:185:length(created))[2:end], length(created))) ]
        vals = query_intervals(created)
        data.count[toreplace] .= get.(Ref(Dict(zip(vals.created, vals.count))), data.created[toreplace], missing)
        data = reduce(vcat, cleanintervals(row) for row in eachrow(data))
        toreplace = ismissing.(data.count)
        created = data.created[toreplace]
    end
    data = prune(data)
    # Dates.canonicalize(Dates.CompoundPeriod(minimum(elem.last - elem.first for elem in data.created)))
    data[!,:created] .= format_tsrange.(data.created)
    data[!,:spdx] .= spdx
    sort!(data, :created)
    execute(conn, "BEGIN;")
    load!(data[!,[:spdx,:created,:count]],
          conn,
          """
          INSERT INTO $schema.queries (spdx, created, count) VALUES ($(join(("\$$i" for i in 1:3), ',')))
          ON CONFLICT ON CONSTRAINT nonoverlappingqueries DO NOTHING;
          """,
          )
    execute(conn, "COMMIT;")
    nothing
end
"""
    find_repo_count_for_intervals(spdx::AbstractString, created::AbstractVector{<:Interval{DateTime}})
"""
function find_repo_count_for_intervals(spdx::AbstractString, created::AbstractVector{<:Interval{DateTime}})
    repositoryCount = zeros(UInt32, length(created))
    indices = range(firstindex(created), lastindex(created), step = 286)
    for idx₀ in indices
        vals = idx₀:min(idx₀ + step(indices) - 1, lastindex(created))
        subsquery = join([ string("_$idx:search(query:\"is:public fork:false mirror:false archived:false license:$spdx created:",
                                  format(created[idx].first, "yyyy-mm-ddTHH:MM:SS+00"),
                                  "..",
                                  format(created[idx].last, "yyyy-mm-ddTHH:MM:SS+00"),
                                  "\",type:REPOSITORY){...A}") for idx in vals
                         ]);
        query = "fragment A on SearchResultItemConnection{repositoryCount} query Search{$subsquery}" |>
            (obj -> replace(obj, r"\s+" => " ")) |>
            strip |>
            string;
        try
            sleep(0.75)
            result = graphql(query)
            json = JSON3.read(result.Data)
            repositoryCount[vals] .= getproperty.(values(json.data), :repositoryCount)
        catch err
            println(vals)
            println(result)
            println(err)
        end
    end
    data = DataFrame(created = created, count = repositoryCount)
    pruned = prune(data)
    cleaned_prune = reduce(vcat, cleanintervals(row) for row in eachrow(pruned))
end
"""
    fill_missing_intervals(spdx::AbstractString, data::AbstractDataFrame)
"""
function fill_missing_intervals(spdx::AbstractString, data::AbstractDataFrame)
    new_data = find_repo_count_for_intervals(spdx, data.created[ismissing.(data.count)])
    new_data = join(data, new_data, on = :created, makeunique = true)
    new_data[!,:count] .= coalesce.(new_data.count_1, new_data.count)
    pruned = prune(new_data[!,[:created, :count]])
    cleaned_prune = reduce(vcat, cleanintervals(row) for row in eachrow(pruned))
end
"""
    find_queries(spdx::AbstractString)
"""
function find_queries(spdx::AbstractString)
    created = vcat(DateTime("2009-01-01") - GH_FIRST_REPO_TS,
                   fill(Month(6), 2),
                   fill(Month(4), 3),
                   fill(Month(2), 6),
                   fill(Month(1), 12),
                   fill(Day(5), 73))
    created = vcat(GH_FIRST_REPO_TS, GH_FIRST_REPO_TS .+ cumsum(created))
    created = vcat(created, created[end]:Day(1):DateTime(year(floor(now(utc_tz), Year))))
    created = Interval{Closed, Open}.(created[begin:end - 1], created[nextind(created, firstindex(created)):end])
    data = find_repo_count_for_intervals(spdx, created)
    while any(ismissing, data.count)
        data = fill_missing_intervals(spdx, data)
    end
    data[!,:spdx] .= spdx
    data[!,[:spdx, :created, :count]]
end


text:
s = open("input") do file
    readlines(file)
end

nums = map(x->parse(Int64,x), s);

lookingFor = Set();

for num in nums
    if(in(num, lookingFor))
        println(num * (2020-num));
        break;
    end

    push!(lookingFor, 2020-num);
end


text:
module RNSpectra
using LsqFit, QuadGK, Plots

function toyozawa(En,α,p0=[1,0.08,5.1,0.2,0.2])
    αf(E,p) = ( p[1].*sqrt( abs(p[2]) )./(E.*(2π) )).*(sum(m-> (2*p[2]/(m^3))*(1/p[4])*exp.(-((E .- (p[3] - p[2]/(m^2))).^2)/(2*(p[4]^2))),1:1000  ) .+ (1/p[5]).*quadgk(x-> exp.(-(x .- E).^2/(2*p[5].^2))./(1 .- exp.(-2π*sqrt( abs(p[2])./(abs(x .- p[3]) ))) ),p[3],Inf,rtol=1e-6)[1]  )

    αₙₒᵣₘ = α/maximum(α)

    #p0=[1,0.08,5.1,0.2,0.2] # [c0 R Eg Γm Γc]
    fit=curve_fit(αf,vec(En),αₙₒᵣₘ,p0)
    cf=coef(fit)
    J=fit.jacobian
    R=fit.resid

    return vec(cf), J, R, vec(αf(En,cf))

end

#function ElliottPlot(En,α,kn)
#        αₙₒᵣₘ = α/maximum(α)
#        plotly()
#        plot()
#        p1=plot(vec(En),vec(αₙₒᵣₘ)*maximum(α))
#        p2=scatter!(vec(En),vec(kn*maximum(α)),dpi=1200,legend=:false)
#        return p2
#
#end

export toyozawa#, ElliottPlot, αk

end