text:
"""
    Selector

Abstract supertype for all selectors.

Selectors are wrappers that indicate that passed values are not the array indices,
but values to be selected from the dimension index, such as `DateTime` objects for
a `Ti` dimension.

Selectors provided in DimensionalData are:

- [`At`](@ref)
- [`Between`](@ref)
- [`Near`](@ref)
- [`Where`](@ref)
- [`Contains`](@ref)

"""
abstract type Selector{T} end

val(sel::Selector) = sel.val
rebuild(sel::Selector, val) = basetypeof(sel)(val)

@inline maybeselector(I...) = maybeselector(I)
@inline maybeselector(I::Tuple) = map(maybeselector, I)
# Int AbstractArray and Colon do normal indexing
@inline maybeselector(i::StandardIndices) = i
# Selectors are allready selectors
@inline maybeselector(i::Selector) = i
# Anything else becomes `At`
@inline maybeselector(i) = At(i)

"""
    At <: Selector

    At(x, atol, rtol)
    At(x; atol=nothing, rtol=nothing)

Selector that exactly matches the value on the passed-in dimensions, or throws an error.
For ranges and arrays, every intermediate value must match an existing value -
not just the end points.

`x` can be any value or `Vector` of values.

`atol` and `rtol` are passed to `isapprox`.
For `Number` `rtol` will be set to `Base.rtoldefault`, otherwise `nothing`,
and wont be used.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(At(20)), Y(At(6))]

# output

5
```
"""
struct At{T,A,R} <: Selector{T}
    val::T
    atol::A
    rtol::R
end
At(val; atol=nothing, rtol=nothing) =
    At{typeof.((val, atol, rtol))...}(val, atol, rtol)

atol(sel::At) = sel.atol
rtol(sel::At) = sel.rtol

"""
    Near <: Selector

    Near(x)

Selector that selects the nearest index to `x`.

With [`Points`](@ref) this is simply the index values nearest to the `x`,
however with [`Intervals`](@ref) it is the interval _center_ nearest to `x`.
This will be offset from the index value for `Start` and
[`End`](@ref) loci.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(Near(23)), Y(Near(5.1))]

# output

4
```
"""
struct Near{T} <: Selector{T}
    val::T
end

"""
    Contains <: Selector

    Contains(x)

Selector that selects the interval the value is contained by. If the
interval is not present in the index, an error will be thrown.

Can only be used for [`Intervals`](@ref) or [`Categorical`](@ref).

## Example

```jldoctest
using DimensionalData

dims_ = X(10:10:20; mode=Sampled(sampling=Intervals())),
        Y(5:7; mode=Sampled(sampling=Intervals()))
A = DimArray([1 2 3; 4 5 6], dims_)
A[X(Contains(8)), Y(Contains(6.8))]

# output

3
```
"""
struct Contains{T} <: Selector{T}
    val::T
end

"""
    Between <: Selector

    Between(a, b)

Selector that retreive all indices located between 2 values,
evaluated with `>=` for the lower value, and `<` for the upper value.
This means the same value will not be counted twice in 2 adjacent
`Between` selections.

For [`Intervals`](@ref) the whole interval must be lie between the
values. For [`Points`](@ref) the points must fall between
the values. Different [`Sampling`](@ref) types may give different
results with the same input - this is the intended behaviour.

`Between` for [`Irregular`](@ref) intervals is a little complicated. The
interval is the distance between a value and the next (for `Start` locus)
or previous (for [`End`](@ref) locus) value.

For [`Center`](@ref), we take the mid point between two index values
as the start and end of each interval. This may or may not make sense for
the values in your indes, so use `Between` with `Irregular` `Intervals(Center())`
with caution.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(Between(15, 25)), Y(Between(4, 6.5))]

# output

1×2 DimArray{Int64,2} with dimensions:
  X: 20:10:20 (Sampled - Ordered Regular Points)
  Y: 5:6 (Sampled - Ordered Regular Points)
 4  5
```
"""
struct Between{T<:Union{Tuple{Any,Any},Nothing}} <: Selector{T}
    val::T
end
Between(args...) = Between(args)

Base.first(sel::Between) = first(val(sel))
Base.last(sel::Between) = last(val(sel))

"""
    Where <: Selector

    Where(f::Function)

Selector that filters a dimension by any function that accepts
a single value from the index and returns a `Bool`.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(19:21)))
A[X(Where(x -> x > 15)), Y(Where(x -> x in (19, 21)))]

# output

1×2 DimArray{Int64,2} with dimensions:
  X: Int64[20] (Sampled - Ordered Regular Points)
  Y: Int64[19, 21] (Sampled - Ordered Regular Points)
 4  6
```
"""
struct Where{T} <: Selector{T}
    f::T
end

val(sel::Where) = sel.f

# sel2indices ==========================================================================

# Converts Selectors to regular indices

@inline sel2indices(x, ls...) = sel2indices(dims(x), ls...)
@inline sel2indices(dims::Tuple, l1, ls...) = sel2indices(dims, (l1, ls...))
@inline sel2indices(dims::Tuple, lookup::Tuple) =
    map((d, l) -> sel2indices(d, l), dims, lookup)
@inline sel2indices(dims::Tuple, lookup::Tuple{}) = ()
@inline sel2indices(dim::Dimension, sel) = _sel2indices(dim, maybeselector(sel))

# First filter based on rough selector properties -----------------

# Standard indices are just returned.
@inline _sel2indices(::Dimension, sel::StandardIndices) = sel
# Vectors are mapped
@inline _sel2indices(dim::Dimension, sel::Selector{<:AbstractVector}) =
    [_sel2indices(mode(dim), dim, rebuild(sel, v)) for v in val(sel)]
@inline _sel2indices(dim::Dimension, sel::Selector) = _sel2indices(mode(dim), dim, sel)

# Where selector ==============================
# Yes this is everything. Where doesn't need mode specialisation
@inline _sel2indices(dim::Dimension, sel::Where) =
    [i for (i, v) in enumerate(index(dim)) if sel.f(v)]

# Then dispatch based on IndexMode -----------------
# Selectors can have varied behaviours depending on the index mode.

# Noindex Contains just converts the selector to standard indices. Implemented
# so the Selectors actually work, not because what they do is useful or interesting.
@inline _sel2indices(mode::NoIndex, dim::Dimension, sel::Union{At,Near,Contains}) = val(sel)
@inline _sel2indices(mode::NoIndex, dim::Dimension, sel::Union{Between}) =
    val(sel)[1]:val(sel)[2]
@inline _sel2indices(mode::Categorical, dim::Dimension, sel::Selector) =
    if sel isa Union{Contains,Near}
        _sel2indices(Points(), mode, dim, At(val(sel)))
    else
        _sel2indices(Points(), mode, dim, sel)
    end
@inline _sel2indices(mode::AbstractSampled, dim::Dimension, sel::Selector) =
    _sel2indices(sampling(mode), mode, dim, sel)

# For Sampled filter based on sampling type and selector -----------------

@inline _sel2indices(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::At) =
    at(sampling, mode, dim, sel)
@inline _sel2indices(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::Near) = begin
    span(mode) isa Irregular && locus(mode) isa Union{Start,End} && _nearirregularerror()
    near(sampling, mode, dim, sel)
end
@inline _sel2indices(sampling::Points, mode::IndexMode, dim::Dimension, sel::Contains) =
    _containspointserror()
@inline _sel2indices(sampling::Intervals, mode::IndexMode, dim::Dimension, sel::Contains) =
    contains(sampling, mode, dim, sel)
@inline _sel2indices(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::Between{<:Tuple}) =
    between(sampling, mode, dim, sel)

@noinline _nearirregularerror() = throw(ArgumentError("Near is not implemented for Irregular with Start or End loci. Use Contains"))
@noinline _containspointserror() = throw(ArgumentError("`Contains` has no meaning with `Points`. Use `Near`"))

# Unaligned IndexMode ------------------------------------------

# unalligned2indices is callled directly from dims2indices

# We use the transformation from the first Transformed dim.
# In practice the others could be empty.
@inline unalligned2indices(dims::DimTuple, sel::Tuple) = sel
@inline unalligned2indices(dims::DimTuple, sel::Tuple{<:Dimension,Vararg{<:Dimension}}) =
    unalligned2indices(dims, map(val, sel))
@inline unalligned2indices(dims::DimTuple, sel::Tuple{<:Selector,Vararg{<:Selector}}) = begin
    coords = [map(val, sel)...]
    transformed = transformfunc(mode(dims[1]))(coords)
    map(_transform2int, sel, transformed)
end

_transform2int(::At, x) = convert(Int, x)
_transform2int(::Near, x) = round(Int, x)

# Selector methods

# at =============================================================================

at(dim::Dimension, sel::At) = at(sampling(mode(dim)), mode(dim), dim, sel)
function at(::Sampling, mode::IndexMode, dim::Dimension, sel::At)
    relate(dim, at(indexorder(dim), val(dim), dim, val(sel), atol(sel), rtol(sel)))
end
#
function at(o::IndexOrder, a::AbstractArray{<:Union{Number,Dates.TimeType}}, dim::Dimension, selval, atol, rtol::Nothing)
    x = unwrap(selval)
    i = searchsortedlast(a, x; order=_ordering(o))
    _isat(x, a[i], atol, 1) || _selvalnotfound(dim, selval)
    return i
end
# range optimisation - we leave bounds checking to the array, but maybe this is not cool
function at(o::IndexOrder, a::AbstractRange{<:Number}, dim::Dimension, selval, atol, rtol::Nothing)
    x = unwrap(selval)
    f = (x - first(a)) / step(a)
    i = unsafe_trunc(Int, f)
    _isat(i, f, atol, step(a)) && checkbounds(Bool, a, i + 1) || _selvalnotfound(dim, selval)
    return i + 1
end

@inline _isat(i, f, atol::Nothing, scalar) = i == f
@inline _isat(i, f, atol, scalar) = abs(i - f) <= atol / scalar

# catch-all for an unordered or non-number index
function at(
    ::IndexOrder, a, dim::Dimension, selval, atol, rtol::Nothing
)
    i = findfirst(x -> x == unwrap(selval), index(dim))
    i == nothing && _selvalnotfound(dim, selval)
    return i
end
# compile-time indexing
@generated function at(
    ::IndexOrder, ::Val{Index}, dim::Dimension, selval::Val{X}, atol, rtol::Nothing
) where {Index,X}
    i = findfirst(x -> x == X, Index)
    if i == nothing
        :(_selvalnotfound(dim, selval))
    else
        return i
    end
end

@noinline _selvalnotfound(dim, selval) = throw(ArgumentError("$selval not found in $dim"))

# near ===========================================================================

# Finds the nearest point in the index, adjusting for locus if necessary.
# In Intevals we are finding the nearest point to the center of the interval.

near(dim::Dimension, sel::Near) = near(sampling(mode(dim)), mode(dim), dim, sel)
function near(::Sampling, mode::IndexMode, dim::Dimension, sel::Near)
    order = indexorder(dim)
    order isa UnorderedIndex && _nearunorderederror()
    locus = DD.locus(dim)

    v = _locus_adjust(locus, unwrap(val(sel)), dim)
    i = _inbounds(_searchorder(order)(order, dim, v), dim)
    i = if (order isa ForwardIndex ? (<=) : (>=))(i, _dimlower(order, dim))
        _dimlower(order, dim)
    else
        previ = _prevind(order, i)
        vl, vi = map(abs, (dim[previ] - v, dim[i] - v))
        # We have to use the right >/>= for Start/End locus
        _lt(locus)(vl, vi) ? previ : i
    end
    relate(dim, i)
end

@noinline _nearunorderederror() = throw(ArgumentError("`Near` has no meaning in an `Unordered` index"))

_locus_adjust(locus::Center, v, dim) = v
_locus_adjust(locus::Start, v, dim) = v - abs(step(dim)) / 2
_locus_adjust(locus::End, v, dim) = v + abs(step(dim)) / 2
_locus_adjust(locus::Start, v::DateTime, dim) = v - (v - (v - abs(step(dim)))) / 2
_locus_adjust(locus::End, v::DateTime, dim) = v + (v + abs(step(dim)) - v) / 2

# contains ================================================================================

# Finds which interval contains a point

function contains(dim::Dimension, sel::Contains)
    contains(sampling(mode(dim)), mode(dim), dim, sel)
end
# Points --------------------------------------
@noinline function contains(::Points, ::IndexMode, dim::Dimension, sel::Contains)
    throw(ArgumentError("Points IndexMode cannot use 'Contains', use 'Near' instead."))
end
# Intervals -----------------------------------
function contains(sampling::Intervals, mode::IndexMode, dim::Dimension, sel::Contains)
    relate(dim, contains(span(mode), sampling, indexorder(mode), locus(mode), dim, sel))
end
# Regular Intervals ---------------------------
function contains(span::Regular, ::Intervals, order, locus, dim::Dimension, sel::Contains)
    v = val(sel); s = abs(val(span))
    _locus_checkbounds(locus, bounds(dim), v)
    i = _whichsearch(locus, order)(order, dim, _maybeaddhalf(locus, s, v))
    # Check the value is in this cell.
    # It is always for AbstractRange but might not be for Val tuple or Vector.
    if !(val(dim) isa AbstractRange)
        _lt(locus)(v, dim[i] + s) || _notcontainederror(v)
    end
    i
end
# Explicit Intervals ---------------------------
function contains(span::Explicit, ::Intervals, order, locus, dim::Dimension, sel::Contains)
    x = val(sel)
    i = searchsortedlast(view(val(span), 1, :), x; order=_ordering(order))
    if i <= 0 || val(span)[2, i] < x
        _notcontainederror(x)
    end
    i
end
# Irregular Intervals -------------------------
function contains(span::Irregular, ::Intervals, order::IndexOrder, locus::Locus, dim::Dimension, sel::Contains)
    _locus_checkbounds(locus, bounds(span), val(sel))
    _whichsearch(locus, order)(order, dim, val(sel))
end
function contains(span::Irregular, ::Intervals, order::IndexOrder, locus::Center, dim::Dimension, sel::Contains)
    v = val(sel)
    _locus_checkbounds(locus, bounds(span), v)
    i = _searchfirst(order, dim, v)
    i <= firstindex(dim) && return firstindex(dim)
    i > lastindex(dim) && return lastindex(dim)

    interval = abs(dim[i] - dim[i - 1])
    distance = abs(dim[i] - v)
    _order_lt(order)(interval / 2, distance) ? i - 1 : i
end

@noinline _notcontainederror(v) = throw(ArgumentError("No interval contains $(v)"))

_whichsearch(::Locus, ::ForwardIndex) = _searchlast
_whichsearch(::Locus, ::ReverseIndex) = _searchfirst
_whichsearch(::End, ::ForwardIndex) = _searchfirst
_whichsearch(::End, ::ReverseIndex) = _searchlast

_maybeaddhalf(::Locus, s, v) = v
_maybeaddhalf(::Center, s, v) = v + s / 2

_order_lt(::ForwardIndex) = (<)
_order_lt(::ReverseIndex) = (<=)

# between ================================================================================

# Finds all values between two points, adjusted for locus where necessary

struct _Upper end
struct _Lower end

between(dim::Dimension, sel::Between) = between(sampling(mode(dim)), mode(dim), dim, sel)
function between(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::Between)
    order = indexorder(dim)
    order isa UnorderedIndex && throw(ArgumentError("Cannot use `Between` with UnorderedIndex"))
    a, b = between(sampling, order, mode, dim, sel)
    relate(dim, a:b)
end
# Points ------------------------------------
function between(sampling::Points, o::IndexOrder, ::IndexMode, dim::Dimension, sel::Between)
    b1, b2 = _maybeflip(o, _sorttuple(sel))
    s1, s2 = _maybeflip(o, (_searchfirst, _searchlast))
    _inbounds((s1(o, dim, b1), s2(o, dim, b2)), dim)
end
# Intervals -------------------------
function between(sampling::Intervals, o::IndexOrder, mode::IndexMode, dim::Dimension, sel::Between)
    between(span(mode), sampling, o, mode, dim, sel)
end
# Regular Intervals -------------------------
function between(span::Regular, ::Intervals, o::IndexOrder, mode::IndexMode, dim::Dimension, sel::Between)
    b1, b2 = _maybeflip(o, _sorttuple(sel) .+ _locus_adjust(mode))
    _inbounds((_searchfirst(o, dim, b1), _searchlast(o, dim, b2)), dim)
end
# Explicit Intervals -------------------------
function between(span::Explicit, ::Intervals, o::IndexOrder, mode::IndexMode, dim::Dimension, sel::Between)
    _inbounds(
        (searchsortedfirst(view(val(span), 1, :), first(val(sel)); order=_ordering(o), lt=<),
         searchsortedlast(view(val(span), 2, :), last(val(sel)); order=_ordering(o), lt=<)),
        dim
    )
end
# Irregular Intervals -----------------------
function between(span::Irregular, ::Intervals, o::IndexOrder, mode::IndexMode, d::Dimension, sel::Between)
    l, h = _sorttuple(sel)
    bl, bh = bounds(span)
    a = l <= bl ? _dimlower(o, d) : between(_Lower(), locus(mode), o, d, l)
    b = h >= bh ? _dimupper(o, d) : between(_Upper(), locus(mode), o, d, h)
    _maybeflip(o, (a, b))
end
function between(x, locus::Union{Start,End}, o::IndexOrder, d::Dimension, v)
    _search(x, o, d, v) - _ordscalar(o) * (_locscalar(locus) + _endshift(x))
end
function between(x, locus::Center, o::IndexOrder, d::Dimension, v)
    r = _ordscalar(o); sh = _endshift(x)
    i = _search(x, o, d, v)
    interval = abs(d[i] - d[i-r])
    distance = abs(d[i] - v)
    # Use the right >/>= to match interval bounds
    _lt(x)(distance, (interval / 2)) ? i - sh * r : i - (1 + sh) * r
end

_locus_adjust(mode) = _locus_adjust(locus(mode), abs(step(span(mode))))
_locus_adjust(locus::Start, step) = zero(step), -step
_locus_adjust(locus::Center, step) = step/2, -step/2
_locus_adjust(locus::End, step) = step, zero(step)

_locscalar(::Start) = 1
_locscalar(::End) = 0
_endshift(::_Lower) = -1
_endshift(::_Upper) = 1
_ordscalar(::ForwardIndex) = 1
_ordscalar(::ReverseIndex) = -1

_search(x, order, dim, v) = _inbounds(_searchorder(order)(order, dim, v, _lt(x)), dim)

_lt(::_Lower) = (<)
_lt(::_Upper) = (<=)

# Shared utils ============================================================================

_searchlast(o::IndexOrder, dim::Dimension, v, lt=<) =
    searchsortedlast(index(dim), unwrap(v); order=_ordering(o), lt=lt)
_searchlast(o::IndexOrder, dim::Dimension{<:Val{Index}}, v, lt=<) where Index =
    searchsortedlast(Index, unwrap(v); order=_ordering(o), lt=lt)

_searchfirst(o::IndexOrder, dim::Dimension, v, lt=<) =
    searchsortedfirst(index(dim), unwrap(v); order=_ordering(o), lt=lt)
_searchfirst(o::IndexOrder, dim::Dimension{<:Val{Index}}, v, lt=<) where Index =
    searchsortedfirst(Index, unwrap(v); order=_ordering(o), lt=lt)

_asfunc(::Type{typeof(<)}) = <
_asfunc(::Type{typeof(<=)}) = <=

# Return an inbounds index
_inbounds(is::Tuple, dim::Dimension) = map(i -> _inbounds(i, dim), is)
_inbounds(i::Int, dim::Dimension) =
    if i > lastindex(dim)
        lastindex(dim)
    elseif i <= firstindex(dim)
        firstindex(dim)
    else
        i
    end

_sorttuple(sel::Between) = _sorttuple(val(sel))
_sorttuple((a, b)) = a < b ? (a, b) : (b, a)

_maybeflip(o::ForwardIndex, (a, b)) = (a, b)
_maybeflip(o::ReverseIndex, (a, b)) = (b, a)

_lt(::Locus) = (<)
_lt(::End) = (<=)
_gt(::Locus) = (>=)
_gt(::End) = (>)

_locus_checkbounds(loc, (l, h), v) =
    (_lt(loc)(v, l) || _gt(loc)(v, h)) && throw(BoundsError())

_prevind(::ForwardIndex, i) = i - 1
_prevind(::ReverseIndex, i) = i + 1

_dimlower(o::ForwardIndex, d) = firstindex(d)
_dimlower(o::ReverseIndex, d) = lastindex(d)
_dimupper(o::ForwardIndex, d) = lastindex(d)
_dimupper(o::ReverseIndex, d) = firstindex(d)

_searchorder(::ForwardIndex) = _searchfirst
_searchorder(::ReverseIndex) = _searchlast

_ordering(::ForwardIndex) = Base.Order.ForwardOrdering()
_ordering(::ReverseIndex) = Base.Order.ReverseOrdering()


text:
include("_common.jl")
include("_calib_funcs.jl")

using Streamfall: load_calibration

using StatsPlots, DataStructures
using Infiltrator
using HypothesisTests

ensemble_figs = joinpath(FIG_PATH, "ensemble_results", "cmd") * "/"
mkpath(ensemble_figs)

ensemble_res_path = joinpath(DATA_PATH, "ensemble_results") * "/"
mkpath(ensemble_res_path)

calib_ts = 3651
state = :storage
state_str = String(state)

Qo = FULL_DATASET[:, "410730_Q"]
Qo_calib = CALIB_OBS[calib_ts+1:CALIB_LN]
Qo_valid = VALID_OBS

# Load metric results
baseline_calib = joinpath(DATA_PATH, "baselines") * "/"
calib_param_path = joinpath(DATA_PATH, "calib_params") * "/"

if !@isdefined best_NmKGE_ensemble
    baseline_NmKGE_results = CSV.File(joinpath(ensemble_res_path, "overall_baseline_NmKGE_results.csv")) |> DataFrame
    state_NmKGE_results = CSV.File(joinpath(ensemble_res_path, "overall_state_NmKGE_results.csv")) |> DataFrame
    base_base_NmKGE_results = CSV.File(joinpath(ensemble_res_path, "overall_base-base_NmKGE_results.csv")) |> DataFrame
    base_state_NmKGE_results = CSV.File(joinpath(ensemble_res_path, "overall_base-state_NmKGE_results.csv")) |> DataFrame
    state_state_NmKGE_results = CSV.File(joinpath(ensemble_res_path, "overall_state-state_NmKGE_results.csv")) |> DataFrame

    baseline_A2k_results = CSV.File(joinpath(ensemble_res_path, "overall_baseline_A2k_results.csv")) |> DataFrame
    state_A2k_results = CSV.File(joinpath(ensemble_res_path, "overall_state_A2k_results.csv")) |> DataFrame
    base_base_A2k_results = CSV.File(joinpath(ensemble_res_path, "overall_base-base_A2k_results.csv")) |> DataFrame
    base_state_A2k_results = CSV.File(joinpath(ensemble_res_path, "overall_base-state_A2k_results.csv")) |> DataFrame
    state_state_A2k_results = CSV.File(joinpath(ensemble_res_path, "overall_state-state_A2k_results.csv")) |> DataFrame

    active_params = CSV.File(joinpath(ensemble_res_path, "base_NmKGE_state.csv")) |> DataFrame
    active_params = active_params[:, "NmKGE"]
end

function split_calib_valid_cols(df)
    calib_df = df[:, [:Overall_calib, :Wet_calib, :Mid_calib, :Dry_calib]]
    valid_df = df[:, [:Overall_valid, :Wet_valid, :Mid_valid, :Dry_valid]]
    rename!(calib_df, Dict(
        "Overall_calib"=>"Overall",
        "Wet_calib"=>"Wet",
        "Mid_calib"=>"Usual",
        "Dry_calib"=>"Dry",
    ))

    rename!(valid_df, Dict(
        "Overall_valid"=>"Overall",
        "Wet_valid"=>"Wet",
        "Mid_valid"=>"Usual",
        "Dry_valid"=>"Dry",
    ))

    return calib_df, valid_df
end

function create_violin(calib, valid, title)
    vplot = violin(title=title, ylabel="NKGE'", label=["Calibration", "Validation"], legend=:bottomleft)
    for (idx, n) in enumerate(names(calib))
        lab = idx > 1 ? "" : "Calibration"
        violin!([n], calib[:, n], side=:left, color="lightblue", label=lab)
        dotplot!([n], calib[:, n], side=:left, marker=(:black,stroke(0)), markersize=3, label="")

        # don't plot if no state in validation period
        if !all(ismissing.(valid[:, n]))
            lab = idx > 1 ? "" : "Validation"
            violin!([n], replace(valid[:, n], missing=>0.0), side=:right, color="orange", label=lab)
            dotplot!([n], replace(valid[:, n], missing=>0.0), side=:right, marker=(:black,stroke(0)), markersize=3, label="")
        end
    end

    return vplot
end

function relabel_approaches(dfs; cols=false)
    replacements = [
        "mean_NmKGE"=>"NμKGE'⁻¹",
        "NmKGE"=>"NKGE'",
        "NnpKGE"=>"NKGEnp"
    ]

    for df in dfs
        if !cols
            for idx in 1:nrow(df)
                app_name = df[idx, :Approach]
                for w in replacements
                    app_name = replace(app_name, w)
                end
                app_name = replace(app_name, "_"=>" ")
                df[idx, :Approach] = app_name
            end
        else
            for w in replacements
                rename!(df, replace.(names(df), w))
            end
            rename!(df, replace.(names(df), "_"=>" "))
        end
    end
end

df_names = ["Baseline", "State-based", "Base-Base", "Base-State", "State-State"]
NmKGE_dfs = [baseline_NmKGE_results, state_NmKGE_results, base_base_NmKGE_results, base_state_NmKGE_results, state_state_NmKGE_results]
A2k_dfs = [baseline_A2k_results, state_A2k_results, base_base_A2k_results, base_state_A2k_results, state_state_A2k_results]

relabel_approaches(A2k_dfs)
relabel_approaches(NmKGE_dfs)

# Find top performing ensemble
constituent_scores = DataFrame(
    :model => [],
    :overall_calib => [],
    :wet_calib => [],
    :usual_calib => [],
    :dry_calib => [],
    :overall_valid => [],
    :wet_valid => [],
    :usual_valid => [],
    :dry_valid => [],
)

function sort_scores(df; col=:Overall_valid, rev=false)
    df = deepcopy(df)

    # Do not consider any model that indicates
    # that drier conditions are not experienced
    filter!(row -> !(ismissing(row.Dry_valid)),  df)

    df[:, :Mean_valid] = mean.(skipmissing.(eachrow(df[:, [:Wet_valid, :Mid_valid, :Dry_valid]])))

    # :Mean_valid,
    sort!(df, [:Mean_valid, col], rev=rev)
end

function create_readable_name(text)
    if occursin(",", text)
        text = join(eval(Meta.parse(text)), " - ")
    end

    return text
end

function create_NmKGE_label(label, Qo, Qm, calib_pos=nothing, valid_pos=nothing)
    local Qo_calib = Qo[calib_ts+1:CALIB_LN]
    local Qo_valid = Qo[CALIB_LN+1:end]
    local Qm_calib = Qm[calib_ts+1:CALIB_LN]
    local Qm_valid = Qm[CALIB_LN+1:end]

    if !isnothing(calib_pos)
        calib_score = round(Streamfall.NmKGE(Qo_calib[calib_pos], Qm_calib[calib_pos]), digits=2)
        valid_score = round(Streamfall.NmKGE(Qo_valid[valid_pos], Qm_valid[valid_pos]), digits=2)
    else
        calib_score = round(Streamfall.NmKGE(Qo_calib, Qm_calib), digits=2)
        valid_score = round(Streamfall.NmKGE(Qo_valid, Qm_valid), digits=2)
    end

    label = "$(label) (C: $(calib_score); V: $(valid_score))"
    return label
end

sorted_NmKGEs = map((df) -> sort_scores(df, col=:Overall_valid, rev=true), NmKGE_dfs)
sorted_A2ks = map((df) -> sort_scores(df, col=:Overall_valid), A2k_dfs)

best_A2k_ensemble = nothing
best_A2k_score = 9999.9
best_A2k_id = nothing
for (df_id, s_A2k) in enumerate(sorted_A2ks)
    score = mean(skipmissing(s_A2k[1, :Mean_valid]))

    println(df_names[df_id], " - Best A2k score: ", score)

    global best_A2k_score
    global best_A2k_id
    global best_A2k_ensemble
    if score < best_A2k_score
        best_A2k_score = score
        best_A2k_ensemble = s_A2k[1, :Approach]
        best_A2k_id = df_id
    end
end

best_performing_A2k_ensemble = sorted_A2ks[best_A2k_id][1, :]
best_performing_A2k_approach = df_names[best_A2k_id]
println("Best Performing: ", best_performing_A2k_approach, " ", best_performing_A2k_ensemble[:Approach], " Score: ", best_A2k_score)

best_A2k_col = create_readable_name(best_A2k_ensemble)

best_NmKGE_ensemble = nothing
best_NmKGE_score = -1.0
best_NmKGE_id = nothing
for (df_id, s_NmKGE) in enumerate(sorted_NmKGEs)
    score = mean(skipmissing(s_NmKGE[1, :Mean_valid]))

    println(df_names[df_id], " - Best NmKGE: ", score)

    global best_NmKGE_score
    global best_NmKGE_id
    if score > best_NmKGE_score
        best_NmKGE_score = score
        best_NmKGE_ensemble = s_NmKGE[1, :Approach]
        best_NmKGE_id = df_id
    end
end

best_performing_baseline = sorted_NmKGEs[1][1, :]
best_performing_NmKGE_ensemble = sorted_NmKGEs[best_NmKGE_id][1, :]
best_performing_NmKGE_name = df_names[best_NmKGE_id]
best_NmKGE_ensemble = best_performing_NmKGE_ensemble[:Approach]
println("Best Performing: ", best_performing_NmKGE_name, " ", best_NmKGE_ensemble, " Score: ", best_NmKGE_score)

base_Qe = CSV.File(joinpath(ensemble_res_path, "baseline_outflows.csv")) |> DataFrame
state_Qe = CSV.File(joinpath(ensemble_res_path, "state_outflows.csv")) |> DataFrame
base_base_Qe = CSV.File(joinpath(ensemble_res_path, "base_base_outflows.csv")) |> DataFrame
base_state_Qe = CSV.File(joinpath(ensemble_res_path, "base_state_outflows.csv")) |> DataFrame
state_state_Qe = CSV.File(joinpath(ensemble_res_path, "state_state_outflows.csv")) |> DataFrame

outflows = [base_Qe, state_Qe, base_base_Qe, base_state_Qe, state_state_Qe]
relabel_approaches(outflows, cols=true)

best_bl_col = create_readable_name(best_performing_baseline[:Approach])
best_bl_display = "Baseline ($best_bl_col)"

best_NmKGE_col = create_readable_name(best_NmKGE_ensemble)
best_NmKGE_display = "$(best_performing_NmKGE_name) ($(best_NmKGE_col))"

base_NmKGE_Q = base_Qe[:, "NKGE'"]
base_NmKGE_calib = base_NmKGE_Q[calib_ts+1:CALIB_LN]
base_NmKGE_valid = base_NmKGE_Q[CALIB_LN+1:end]

best_bl_Q = base_Qe[:, best_bl_col]
best_bl_calib = best_bl_Q[calib_ts+1:CALIB_LN]
best_bl_valid = best_bl_Q[CALIB_LN+1:end]

best_NmKGE_Q = outflows[best_NmKGE_id][:, best_NmKGE_col]
best_NmKGE_Q_calib = best_NmKGE_Q[calib_ts+1:CALIB_LN]
best_NmKGE_Q_valid = best_NmKGE_Q[CALIB_LN+1:end]

base_NmKGE_label = create_NmKGE_label("Baseline NKGE'", Qo, base_NmKGE_Q)
best_bl_label = create_NmKGE_label(best_bl_display, Qo, best_bl_Q)
best_NmKGE_label = create_NmKGE_label(best_NmKGE_display, Qo, best_NmKGE_Q)

res_qq = qqplot(Qo_valid, base_NmKGE_valid, label=base_NmKGE_label, xlabel="Observed [ML]", ylabel="Modeled Flow [ML]")
qqplot!(Qo_valid, best_bl_valid, label=best_bl_label, xlabel="Observed [ML]", ylabel="Modeled Flow [ML]", alpha=0.5)
qqplot!(Qo_valid, best_NmKGE_Q_valid, alpha=0.4, label=best_NmKGE_label)
savefig("$(ensemble_figs)comparison_baseline_best_ensembles_valid_qqplot.png")

res_scatter = scatter(Qo_valid, base_NmKGE_valid .- Qo_valid, label=base_NmKGE_label, xlabel="Observed [ML]", ylabel="Modeled - Observed [ML]")
scatter!(Qo_valid, best_bl_valid - Qo_valid, alpha=0.5, label=best_bl_label)
scatter!(Qo_valid, best_NmKGE_Q_valid, alpha=0.4, label=best_NmKGE_label)
savefig("$(ensemble_figs)comparison_baseline_best_ensembles_valid_residuals.png")

#######

# Read in state params
function separate_state_periods(activity)
    low = findall(activity .== 1)
    mid = findall(activity .== 2)
    high = findall(activity .== 3)

    return low, mid, high
end

function separate_state_periods(activity, t_min, t_max=nothing)
    if !isnothing(t_max)
        active = activity[t_min:t_max]
    else
        active = activity[t_min:end]
    end

    return separate_state_periods(active)
end

# Use baseline NmKGE catchment conditions
calib_low, calib_mid, calib_high = separate_state_periods(active_params, calib_ts+1, CALIB_LN)
valid_low, valid_mid, valid_high = separate_state_periods(active_params, CALIB_LN+1)

# Show each state
base_dry_valid = base_NmKGE_valid[valid_high]
bl_dry_valid = best_bl_valid[valid_high]
NmKGE_dry_valid = best_NmKGE_Q_valid[valid_high]
Qo_dry_valid = Qo_valid[valid_high]
dry_scatter = scatter(Qo_dry_valid, base_dry_valid .- Qo_dry_valid, label=base_NmKGE_label, xlabel="Observed [ML]", ylabel="Modeled [ML]")

base_NmKGE_mid_label = create_NmKGE_label("Baseline NKGE'", Qo, base_NmKGE_Q, calib_mid, valid_mid)
best_bl_mid_label = create_NmKGE_label(best_bl_display, Qo, best_bl_Q, calib_mid, valid_mid)
best_NmKGE_mid_label = create_NmKGE_label(best_NmKGE_display, Qo, best_NmKGE_Q, calib_mid, valid_mid)
base_mid_qq = qqplot(Qo_valid[valid_mid], base_NmKGE_valid[valid_mid], label=base_NmKGE_mid_label, xlabel="Observed [ML]", ylabel="Modeled Flow [ML]")
bl_mid_qq = qqplot!(Qo_valid[valid_mid], best_bl_valid[valid_mid], label=best_bl_mid_label, alpha=0.3)
NmKGE_mid_qq = qqplot!(Qo_valid[valid_mid], best_NmKGE_Q_valid[valid_mid], label=best_NmKGE_mid_label, alpha=0.3)
savefig("$(ensemble_figs)comparison_usual_valid_qq.png")

base_NmKGE_dry_label = create_NmKGE_label("Baseline NKGE'", Qo, base_NmKGE_Q, calib_high, valid_high)
bl_NmKGE_dry_label = create_NmKGE_label(best_bl_display, Qo, best_bl_Q, calib_high, valid_high)
best_NmKGE_dry_label = create_NmKGE_label(best_NmKGE_display, Qo, best_NmKGE_Q, calib_high, valid_high)
base_dry_qq = qqplot(Qo_dry_valid, base_dry_valid, label=base_NmKGE_dry_label, xlabel="Observed [ML]", ylabel="Modeled Flow [ML]", legend=:topleft, log=true)
bl_dry_qq = qqplot!(Qo_dry_valid, bl_dry_valid, label=bl_NmKGE_dry_label, alpha=0.3)
NmKGE_dry_qq = qqplot!(Qo_dry_valid, NmKGE_dry_valid, label=best_NmKGE_dry_label, alpha=0.3)
savefig("$(ensemble_figs)comparison_dry_valid_qq.png")

combined_state_qq = plot(
    plot(base_mid_qq, title="Usual", legend=:bottomright),
    plot(base_dry_qq, title="Dry", legend=:bottomright),
    layout=(1, 2),
    margin=5Plots.mm,
    size=(1000,400),
    markerstrokewidth=0,
    xaxis=:log,
    yaxis=:log
)
savefig(combined_state_qq, "$(ensemble_figs)combined_state_qq.png")
#####

base_NmKGE_label = create_NmKGE_label("Baseline NKGE'", Qo, base_NmKGE_Q)
best_NmKGE_label = create_NmKGE_label(best_NmKGE_display, Qo, best_NmKGE_Q)
res_scatter2 = scatter(Qo_valid, base_NmKGE_valid - Qo_valid, label=base_NmKGE_label,
                       xlabel="Observed [ML]", ylabel="Modeled - Observed [ML]")
scatter!(Qo_valid, best_bl_valid - Qo_valid, alpha=0.5, label=best_bl_label)
scatter!(Qo_valid, best_NmKGE_Q_valid - Qo_valid, alpha=0.4, label=best_NmKGE_label)

combined_res = plot(
    plot(res_qq, xaxis=:log, yaxis=:log, legend=:bottomright),
    plot(res_scatter2, legend=false),
    layout=(1, 2),
    margin=5Plots.mm,
    size=(1000,400),
    #link=:y,
    markerstrokewidth=0
)
savefig(combined_res, "$(ensemble_figs)combined_residuals.png")

layout_4panel = @layout [a b; c]
combined_4panel = plot(
    plot(res_qq, title="Q-Q Plot", xaxis=:log, yaxis=:log, legend=:bottomright),
    plot(res_scatter2, title="Residuals", legend=false),
    combined_state_qq,
    layout=layout_4panel,
    margin=5Plots.mm,
    size=(1000,860),
    markerstrokewidth=0
)
savefig(combined_4panel, "$(ensemble_figs)combined_residuals_4panel.png")

# QQ of Best performing instance out of each approach
base_NmKGE_label = create_NmKGE_label("Baseline NKGE'", Qo, base_NmKGE_Q)
best_NmKGE_qq = qqplot(Qo_valid, base_Qe[:, "NKGE'"][CALIB_LN+1:end], label=base_NmKGE_label,
                       legend=:bottomright, xaxis=:log, yaxis=:log)
flow_data = [base_Qe, state_Qe, base_base_Qe, base_state_Qe, state_state_Qe]
for (df_name, flow, df) in zip(df_names, flow_data, sorted_NmKGEs)
    best_name = create_readable_name(df[1, :Approach])

    if (df_name == "Baseline") && (best_name == "NKGE'")
        continue
    end

    best_NmKGE_label = create_NmKGE_label(best_name, Qo, flow[:, best_name])
    qqplot!(Qo_valid, flow[CALIB_LN+1:end, best_name], label="$(df_name) [$(best_NmKGE_label)]",
            alpha=0.3, legend=:bottomright, markerstrokewidth=0)
end
savefig(best_NmKGE_qq, "$(ensemble_figs)best_NmKGE_qq.png")

# QQ based on A2k score
# best_A2k_qq = qqplot(Qo_valid, base_Qe[:, "NKGE'"][CALIB_LN+1:end], label="Baseline NKGE'", legend=:topleft)
# for (df_name, flow, df) in zip(df_names, flow_data, sorted_A2ks)
#     best_name = create_readable_name(df[1, :Approach])

#     if (df_name == "Baseline") && (best_name == "NKGE'")
#         continue
#     end

#     ME_calib = round(mean(Qo_calib - flow[calib_ts+1:CALIB_LN, best_name]), digits=2)
#     ME_valid = round(mean(Qo_valid - flow[CALIB_LN+1:end, best_name]), digits=2)
#     qqplot!(Qo_valid, flow[CALIB_LN+1:end, best_name], label="$(df_name) ($(best_name); C: $(ME_calib) | V: $(ME_valid))",
#             alpha=0.3, legend=:topleft, markerstrokewidth=0)
# end
# savefig(best_A2k_qq, "$(ensemble_figs)best_A2k_qq.png")

# Create violin plots
bl_NmKGE_calib, bl_NmKGE_valid = split_calib_valid_cols(sorted_NmKGEs[1])  # baseline NmKGEs
state_NmKGE_calib, state_NmKGE_valid = split_calib_valid_cols(sorted_NmKGEs[2])  # state-based NmKGEs
base_base_NmKGE_calib, base_base_NmKGE_valid = split_calib_valid_cols(sorted_NmKGEs[3])  # base-base NmKGEs
base_state_NmKGE_calib, base_state_NmKGE_valid = split_calib_valid_cols(sorted_NmKGEs[4])  # base-state NmKGEs
state_state_NmKGE_calib, state_state_NmKGE_valid = split_calib_valid_cols(sorted_NmKGEs[5])  # state-state NmKGEs

baseline_vplot = create_violin(bl_NmKGE_calib, bl_NmKGE_valid, "Baseline")
savefig(baseline_vplot, "$(ensemble_figs)baseline_NmKGE_results.png")

state_based_vplot = create_violin(state_NmKGE_calib, state_NmKGE_valid, "State-based")
savefig(state_based_vplot, "$(ensemble_figs)state-based_NmKGE_results.png")

base_base_vplot = create_violin(base_base_NmKGE_calib, base_base_NmKGE_valid, "Base-Base")
savefig(base_base_vplot, "$(ensemble_figs)base-base_NmKGE_results.png")

base_state_vplot = create_violin(base_state_NmKGE_calib, base_state_NmKGE_valid, "Base-State")
savefig(base_state_vplot, "$(ensemble_figs)base-state_NmKGE_results.png")

state_state_vplot = create_violin(state_state_NmKGE_calib, state_state_NmKGE_valid, "State-State")
savefig(state_state_vplot, "$(ensemble_figs)state-state_NmKGE_results.png")

combined_violin = plot(
    baseline_vplot,
    plot(state_based_vplot, ylabel="", legend=false),
    plot(state_state_vplot, ylabel="", legend=false),
    layout=(1, 3),
    margin=5Plots.mm,
    size=(1200,400),
    link=:y
)
savefig(combined_violin, "$(ensemble_figs)combined_violin.png")

combined_violin_all = plot(
    baseline_vplot,
    plot(state_based_vplot, ylabel="", legend=false),
    plot(base_base_vplot, ylabel="", legend=false),
    plot(base_state_vplot, ylabel="", legend=false),
    plot(state_state_vplot, ylabel="", legend=false),
    layout=(1, 5),
    margin=5Plots.mm,
    size=(2000,400),
    link=:y
)
savefig(combined_violin_all, "$(ensemble_figs)combined_violin_all.png")

ext: jl


text:
module ApplicationBuilder

using Glob, PackageCompiler

export build_app_bundle

@static if Sys.isapple()
    include("sign_mac_app.jl")
    include("mac_commandline_app.jl")
end

"""
    build_app_bundle(juliaprog_main;
        appname, builddir, binary_name, resources, libraries, verbose, bundle_identifier,
        app_version, icns_file, certificate, entitlements_file, snoopfile, autosnoop,
        commandline_app, cpu_target)

Compile `juliaprog_main` into an executable, and bundle it together with all its
`resources` and `libraries` into an App called `appname`.

juliaprog_main: Path to a ".jl" file that defines the function `julia_main()`

# Examples
```julia-repl
julia> build_app_bundle("main.jl", appname="MyApp", resources=["img.jpg"],
           libraries=[MyPackage._libp])
```
"""
function build_app_bundle(juliaprog_main;
        appname = splitext(basename(juliaprog_main))[1], builddir = "builddir",
        binary_name = splitext(basename(juliaprog_main))[1],
        resources = String[], libraries = String[], verbose = false,
        bundle_identifier = nothing, app_version = "0.1", icns_file = nothing,
        certificate = nothing, entitlements_file = nothing,
        snoopfile = nothing, autosnoop = false, cpu_target="x86-64",
        create_installer = false, commandline_app = false,
    )

    # ----------- Input sanity checking --------------

    builddir = abspath(builddir)
    @static if Sys.isapple()
        if bundle_identifier == nothing
            bundle_identifier = make_bundle_identifier(appname)
            println("  Using calculated bundle_identifier: '$bundle_identifier'")
        end
        # Bundle identifier requirements: https://apple.stackexchange.com/a/238381/52530
        if occursin(r"\s", bundle_identifier) throw(ArgumentError("Bundle identifier must not contain whitespace.")) end
        if occursin(r"[^A-Za-z0-9-.]", bundle_identifier) throw(ArgumentError("Bundle identifier must contain only alphanumeric characters (A-Z,a-z,0-9), hyphen (-), and period (.).")) end

    else
        if commandline_app
            @warn "Ignore `commandline_app=true` on non-macOS system."
            commandline_app = false
        end
    end

    @static if !Sys.iswindows()
        if create_installer
            @warn "Ignoring `create_installer=true` on non-Windows system."
            create_installer = false
        end
    end

    if !isfile(juliaprog_main) throw(ArgumentError("Cannot build application. No such file '$juliaprog_main'")) end

    # ----------- Initialize App ---------------------
    appext = @static if Sys.isapple() ".app" else "" end
    appbundle = joinpath(builddir, appname * appext)
    println("~~~~~~ Creating App bundle in \"$appbundle\" ~~~~~~~")

    app_dir, launcher_dir, resources_dir, libs_dir = @static if Sys.isapple()
        ("$appbundle/Contents",
        "$appbundle/Contents/MacOS",
        "$appbundle/Contents/Resources",
        "$appbundle/Contents/Libraries",)
    else
        (appbundle,
        joinpath(appbundle, "bin"),
        joinpath(appbundle, "res"),
        joinpath(appbundle, "lib"),)
    end

    applet_name = nothing
    if commandline_app  # MacOS only
        # TODO: What if the user specifies Resources that could overwrite
        #  applet resources? (ie Scripts/ or applet.rsrc)
        applet_name = build_commandline_app_bundle(builddir, binary_name, appname, verbose)
    end

    function has_files(files)
        return !isempty(files) && !all(isempty, [strip(f) for f in files])
    end

    mkpath(launcher_dir)
    mkpath(resources_dir)  # Always create resources directory (For .icns file)
    has_files(libraries) && mkpath(libs_dir)

    # ----------- Copy user libs & assets -------------
    run_verbose(verbose, cmd::Cmd) = (verbose && println("    julia> run($cmd)") ; run(cmd))

    if has_files(libraries) || has_files(resources)
        # Copy assets and libs early so they're available during compilation.
        #  This is mostly relevant if you have .dylibs in your assets, which the compiler needs to look at.

        println("~~~~~~ Copying user-specified libraries & resources to bundle... ~~~~~~~")

        println("  Resources:")
        for res in resources
            res = clean_file_pattern(res, "-R")
            print("    - $res ...")
            copy_file_dir_or_glob(res, resources_dir)
            println("............ done")
        end
        println("  Libraries:")
        for lib in libraries
            lib = clean_file_pattern(lib, "-L")
            print("    - $lib ...")
            copy_file_dir_or_glob(lib, libs_dir)
            println("............ done")
        end
    end

    # ----------- Compile a binary ---------------------
    # Compile the binary right into the app.
    println("~~~~~~ Compiling a binary from '$juliaprog_main'... ~~~~~~~")

    if autosnoop
        if snoopfile != nothing
            println("WARNING: autosnoop is overwriting user-specified snoopfile.")
        end
        snoopfile = "$builddir/$appname-snoopfile.jl"
        write(snoopfile, """Base.include(@__MODULE__, raw"$(abspath("$juliaprog_main"))");  julia_main([""]); """)
    end

    custom_program_c = "$(@__DIR__)/program.c"
    cc_flags = Sys.isapple() ? `-mmacosx-version-min=10.10 -headerpad_max_install_names` : nothing
    # Provide an environment variable telling the code it's being compiled into a mac bundle.
    withenv("LD_LIBRARY_PATH"=>"$libs_dir:$libs_dir/julia",
            "COMPILING_APPLE_BUNDLE"=>"true") do
        verbose && println("  PackageCompiler.static_julia(...)")
        # Compile executable and copy julia libs to $launcher_dir.
        PackageCompiler.build_executable(juliaprog_main, binary_name, custom_program_c;
                builddir=launcher_dir, verbose=verbose, optimize="3",
                snoopfile=snoopfile, debug="0", cpu_target=cpu_target,
                compiled_modules="yes",
                cc_flags=cc_flags)
    end

    @static if Sys.isapple()
        for b in ["$launcher_dir/$binary_name", "$launcher_dir/$binary_name.dylib"]
            run_verbose(verbose, `install_name_tool -add_rpath "@executable_path/../Frameworks/" $b`)
            run_verbose(verbose, `install_name_tool -add_rpath "@executable_path/../Libraries/" $b`)
        end

        # In order to pass Apple's GateKeeper, the App must not reference any external Julia libs.
        for binary_file in glob("*", launcher_dir)
            try
                #  an example output line from otool: "         path /Applications/Dev Apps/Julia-0.6.app/Contents/Resources/julia/lib (offset 12)"
                external_julia_deps = readlines(pipeline(`otool -l $binary_file`,
                     `grep $(dirname(Sys.BINDIR))`,  # filter julia lib deps
                     `sed 's/\s*path//'`, # remove leading "  path"
                     `sed 's/(.*)$//'`)) # remove trailing parens
                for line in external_julia_deps
                    path = strip(line)
                    run_verbose(verbose, `install_name_tool -delete_rpath "$path" $binary_file`)
                end
            catch
            end
            # Also need to strip any non-x86 architectures to make Apple happy.
            # (It looks like this only affects libgcc_s.1.dylib.)
            try
                if success(pipeline(`file $binary_file`, `grep 'i386'`))
                    run_verbose(verbose, `lipo $binary_file -thin x86_64 -output $binary_file`)
                end
            catch
            end
        end

        # ---------- Create Info.plist to tell it where to find stuff! ---------
        # This lets you have a different app name from your juliaprog_main.
        println("~~~~~~ Generating 'Info.plist' for '$bundle_identifier'... ~~~~~~~")

        info_plist() = """
            <?xml version="1.0" encoding="UTF-8"?>
            <!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
            <plist version="1.0">
            <dict>
            	<key>CFBundleAllowMixedLocalizations</key>
            	<true/>
            	<key>CFBundleDevelopmentRegion</key>
            	<string>en</string>
            	<key>CFBundleDisplayName</key>
            	<string>$appname</string>
            	<key>CFBundleExecutable</key>
            	<string>$(commandline_app ? applet_name : binary_name)</string>
            	<key>CFBundleIconFile</key>
            	<string>$appname.icns</string>
            	<key>CFBundleIdentifier</key>
            	<string>$bundle_identifier</string>
                <key>CFBundleInfoDictionaryVersion</key>
                <string>6.0</string>
            	<key>CFBundleName</key>
            	<string>$appname</string>
            	<key>CFBundlePackageType</key>
            	<string>APPL</string>
            	<key>CFBundleShortVersionString</key>
            	<string>$app_version</string>
            	<key>CFBundleVersion</key>
            	<string>$app_version</string>
            	<key>NSHighResolutionCapable</key>
                <string>YES</string>
            	<key>LSMinimumSystemVersionByArchitecture</key>
            	<dict>
            		<key>x86_64</key>
            		<string>10.10</string>
            	</dict>
            	<key>LSRequiresCarbon</key>
            	<true/>
            	<key>NSHumanReadableCopyright</key>
            	<string>© 2018 $bundle_identifier </string>
            </dict>
            </plist>
            """

        verbose && println(info_plist())
        write("$app_dir/Info.plist", info_plist());

        # Copy Julia icons
        julia_app_resources_dir() = joinpath(Sys.BINDIR, "..","..")
        if (icns_file == nothing)
            icns_file = joinpath(julia_app_resources_dir(),"julia.icns")
            verbose && println("Attempting to copy default icons from Julia.app: $icns_file")
        end
        if isfile(icns_file)
            cp(icns_file, "$resources_dir/$appname.icns", force=true);
        else
            @warn "Skipping nonexistent icons file: '$icns_file'"
        end

        # --------------- CLEAN UP before distributing ---------------
        println("~~~~~~ Cleaning up temporary files... ~~~~~~~")

        # Delete the tmp build files
        function delete_if_present(file, path)
            files = glob(file, path)
            if !isempty(files) run(`rm -r $(files)`) end
        end
        delete_if_present("*.ji",launcher_dir)
        delete_if_present("*.o",launcher_dir)

        # Remove debug .dylib libraries and any precompiled .ji's
        delete_if_present("*.dSYM",libs_dir)
        delete_if_present("*.dSYM","$libs_dir/julia")
        delete_if_present("*.backup","$libs_dir/julia")
        delete_if_present("*.ji","$libs_dir/julia")
        delete_if_present("*.o","$libs_dir/julia")

        if certificate != nothing
            println("~~~~~~ Signing the binary and all libraries ~~~~~~~")
            sign_application_libs(launcher_dir, certificate)
            if entitlements_file != nothing
                set_entitlements("$launcher_dir/$binary_name", certificate, entitlements_file)
            end
        end
    elseif Sys.iswindows()
        create_installer && win_installer(builddir, name = appname)
    end

    println("~~~~~~ Done building '$appbundle'! ~~~~~~~")
    return 0
end

function make_bundle_identifier(appname)
    cleanregex = r"[^a-zA-Z0-9]"
    cleanuser = replace(lowercase(ENV["USER"]), cleanregex => "")
    cleanapp = replace(lowercase(appname), cleanregex => "")
    "com.$cleanuser.$cleanapp"
end

function copy_file_dir_or_glob(pattern, destdir)
    if isfile(pattern)  # Copy the file to destdir
        cp(pattern, joinpath(destdir, basename(pattern)), force=true)
    elseif isdir(pattern)  # Copy the entire *dir* (not its contents) to destdir.
        cp(pattern, joinpath(destdir, basename(pattern)), force=true)
    # if not a file or a directory, assume it's a glob.
    elseif pattern[1] == '/'  # Glob.glob can't handle globs starting with '/'
        files = glob(pattern[2:end], pattern[1:1])
        for file in files
            cp(file, joinpath(destdir, basename(file)), force=true)
        end
    elseif !isempty(pattern)  # Copy the specified glob pattern file[s] to destdir.
        for file in glob(pattern)
            cp(file, joinpath(destdir, basename(file)), force=true)
        end
    else
        @warn "Skipping unknown file '$pattern'!"
    end
end

function clean_file_pattern(pattern, errorMsg_fileCmd)
    clean_pattern = strip(pattern)
    if isempty(clean_pattern)
        throw(ArgumentError("ERROR: $errorMsg_fileCmd '$pattern' must not be empty."))
    elseif clean_pattern[1] == '~'
        clean_pattern = homedir() * clean_pattern[2:end]
    end
    return clean_pattern
end

end  # module


text:
function compute_char(board_array, row, col)
    board_array[row, col] && return '*'
    rows, cols = size(board_array)
    row_range = max(1, row-1):min(rows, row+1)
    col_range = max(1, col-1):min(cols, col+1)
    count = sum(board_array[row_range, col_range])
    count > 0 ? '0' + count : ' '
end

function annotate(board)
    isempty(board) && return []
    isempty(board[1]) && return [""]
    board_array = map(c -> c == '*', mapreduce(collect, hcat, board))'
    board_array = map(i -> compute_char(board_array, i[1], i[2]), CartesianIndices(board_array))
    map(row -> join(board_array[row,:]), 1:length(board))
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
using CompactBasisHelFEM
using Test

@testset "CompactBasisHelFEM.jl" begin
    # Write your tests here.
end


text:
using Statistics: mean

@testset "softmax integer input" begin
    @test softmax(Int[0, 0]) == [0.5, 0.5]
end

@testset "softmax on different dims" begin
    xs = rand(fill(2, 5)...)
    out = similar(xs)
    for (fn!, fn) in [(softmax!, softmax), (logsoftmax!, logsoftmax)], i = 1:ndims(xs)
        @test fn!(out, xs; dims = i) == fn(xs; dims = i)
    end
end

@testset "softmax" begin
    xs = rand(5, 5)
    @test all(sum(softmax(xs), dims = 1) .≈ 1)
    @test all(sum(softmax(xs; dims = 2), dims = 2) .≈ 1)
    @test sum(softmax(vec(xs))) ≈ 1
    @test log.(softmax(xs; dims = 2)) ≈ logsoftmax(xs; dims = 2)

    xs = [-100_000.0, -100_000.0]
    @test softmax(xs) ≈ [0.5, 0.5]
    @test logsoftmax(xs) ≈ log.([0.5, 0.5])

    xs = rand(5)
    @test softmax(xs) ≈ exp.(xs) ./ sum(exp.(xs))
    @test logsoftmax(xs) ≈ log.(softmax(xs))

    xs = Float32[1, 2, 3000.0]
    @test logsoftmax(xs) ≈ [-2999, -2998, 0]

    xs = Float32[1 2 3; 1000 2000 3000]
    @test logsoftmax(xs) ≈ [-999 -1998 -2997; 0 0 0.0]

    y = logsoftmax(xs)
    @test ∇logsoftmax(ones(Float32, size(xs)), xs, y) ≈ Float32[1 1 1; -1 -1 -1]

    y = softmax(xs)
    @test ∇softmax(ones(Float32, size(xs)), xs, y) ≈ zeros(Float32, size(xs))

    # These values precalculated using PyTorch's nn.LogSoftmax
    xs = [
        -0.238639 0.748142 -0.283194 -0.525461 -1.5348 -0.797842
        0.690384 0.211427 0.254794 -0.213572 -0.314174 -0.372663
        -1.146370 -0.577988 0.718952 0.919720 -0.620773 0.929977
    ]
    ys = [
        0.237703 -0.621474 0.448193 0.546047 0.564185 0.632273
        -0.930163 0.0519798 0.0549979 0.3799 -0.477112 0.437428
        0.69246 0.569494 -0.503191 -0.925947 -0.0870738 -1.0697
    ]

    y = logsoftmax(xs)
    @test ∇logsoftmax(ones(size(xs)), xs, y) ≈ ys rtol = 1e-6

    y = softmax(xs)
    @test ∇softmax(ones(size(xs)), xs, y) ≈ zeros(size(xs)) atol = 1e-6
end

@testset "softmax with Inf, NaN" begin
    @test softmax(Float32[1 2; 3 Inf]) ≈    Float32[0.11920292 0.0; 0.880797 1.0]
    @test softmax(Float32[1 -Inf; 3 Inf]) ≈ Float32[0.11920292 0.0; 0.880797 1.0]
    @test softmax(Float32[1 Inf; 3 Inf]) ≈  Float32[0.11920292 0.5; 0.880797 0.5]

    @test softmax(Float32[1 2; 3 NaN]) ≈    Float32[0.11920292 NaN; 0.880797 NaN] nans=true
    @test softmax(Float32[1 2; 3 Inf]; dims=2) ≈ Float32[0.26894143 0.7310586; 0.0 1.0]
    @test softmax(Float32[1 2; 3 Inf]; dims=(:)) ≈ Float32[0.0 0.0; 0.0 1.0]
    @test softmax(Float32[1 2; 3 Inf]; dims=(1,2)) ≈ Float32[0.0 0.0; 0.0 1.0]

    @test exp.(logsoftmax(Float32[1 2; 3 Inf])) ≈ softmax(Float32[1 2; 3 Inf])
    @test exp.(logsoftmax(Float32[1 -Inf; 3 Inf])) ≈ softmax(Float32[1 -Inf; 3 Inf])
    @test exp.(logsoftmax(Float32[1 Inf; 3 Inf])) ≈ softmax(Float32[1 Inf; 3 Inf])
end

@testset "mutating softmax" begin
    map([
        Float64[1 2 3; 5 6 7],
        Float64[
            -0.238639 0.748142 -0.283194 -0.525461 -1.5348 -0.797842
            0.690384 0.211427 0.254794 -0.213572 -0.314174 -0.372663
            -1.146370 -0.577988 0.718952 0.919720 -0.620773 0.929977
        ],
    ]) do xs
        out = similar(xs)
        softmax!(out, xs)
        @test out ≈ softmax(xs) rtol = 1e-6
        logsoftmax!(out, xs)
        @test out ≈ logsoftmax(xs) rtol = 1e-6

        @testset "$fn(Float64, $(size(xs)))" for fn in [zeros, ones, rand]
            Δ = fn(Float64, size(xs))
            y = softmax(xs)
            ∇softmax!(out, Δ, xs, y)
            @test out ≈ ∇softmax(Δ, xs, y)  rtol = 1e-6

            y = logsoftmax(xs)
            ∇logsoftmax!(out, Δ, xs, y)
            @test out ≈ ∇logsoftmax(Δ, xs, y)  rtol = 1e-6
        end
    end
end

@testset "logsumexp" begin
    flogsoft(x; dims) = mean(x .- logsoftmax(x; dims = dims), dims = dims)

    x = rand(3, 4)
    @test logsumexp(x) ≈ flogsoft(x, dims = :)
    @test logsumexp(x; dims = 1) ≈ flogsoft(x, dims = 1)
end

@testset "AutoDiff" begin
    for f in (softmax, logsoftmax), d in (:, 1, 2)
        gradtest(f, (3,4); fkwargs=(; dims=d), check_rrule=true)
    end
    gradtest(x -> softmax(x).*(1:3), 3)
    gradtest(x -> softmax(x).*(1:3), (3,5), atol=1e-4)
    gradtest(x -> softmax(x, dims=2).*(1:3), (3,5), atol=1e-4)
    gradtest(x -> logsoftmax(x).*(1:3), 3)
    gradtest(x -> logsoftmax(x).*(1:3), (3,5))
    gradtest(x -> logsoftmax(x, dims=2).*(1:3), (3,5))

    for d  in (:, 1, 2)
        gradtest(logsumexp, (3,4), fkwargs=(; dims=d))
    end
end


text:
using BenchmarkTools
using JultraDark
using LinearAlgebra
using MPI

if ~MPI.Initialized()
    MPI.Init()
end

n = 64

suite = BenchmarkGroup()

function fft_broadcast!(grids)
    grids.ψk .= grids.fft_plan * grids.ψx
end

function fft_inplace!(grids)
    mul!(grids.ψk, grids.fft_plan, grids.ψx)
end

function fft_inv_broadcast!(grids)
    grids.ψx .= grids.fft_plan \ grids.ψk
end

function fft_inv_inplace!(grids)
    ldiv!(grids.ψx, grids.fft_plan, grids.ψk)
end

for (label, grids) in [("Grids", Grids(1.0, n)), ("PencilGrids", PencilGrids(1.0, n))]

    suite[label] = BenchmarkGroup()

    suite[label]["fft_broadcast!"] = @benchmarkable fft_broadcast!($grids)
    suite[label]["fft_inplace!"] = @benchmarkable fft_inplace!($grids)
    suite[label]["fft_inv_broadcast!"] = @benchmarkable fft_inv_broadcast!($grids)
    suite[label]["fft_inv_inplace!"] = @benchmarkable fft_inv_inplace!($grids)

end

@show tune!(suite)

results = run(suite)

@show res_min = minimum(results)
@show res_med = median(results)

judge(res_med["Grids"], res_med["PencilGrids"])


text:
using NiLang
using ReversibleSeismic
using Test

@testset "integrate" begin
    k = 20
    nx = ny = 100
    N = 1999

    upre = zeros(nx+2, ny+2)
    u = zeros(nx+2, ny+2)
    φ = zeros(nx+2, ny+2)
    ψ = zeros(nx+2, ny+2)

    param = AcousticPropagatorParams(nx=nx, ny=ny,
        Rcoef=0.2, dx=20.0, dy=20.0, dt=0.05, nstep=N)
    c = 1000*ones(nx+2, ny+2)
    src = (nx ÷ 2, ny ÷ 2)
    srcv = Ricker(param, 100.0, 500.0)

    x = ReversibleSeismic.SeismicState(upre, u, φ, ψ, Ref(0))

    x_last = NiLang.direct_emulate(bennett_step!, copy(x), param, src, srcv, copy(c); N=N)
    x_last_b = bennett(bennett_step!, zero(x), copy(x), param, src, srcv, copy(c); k=k, N=N)[2]
    @test isapprox(x_last.u, x_last_b.u; atol=1e-5)
    @test sum(abs2, x_last.u) ≈ 6.234873084873294e-5
    @test isapprox(x_last.ψ, x_last_b.ψ; atol=1e-5)
    @test isapprox(x_last.φ, x_last_b.φ; atol=1e-5)
    @test isapprox(x_last.upre, x_last_b.upre; atol=1e-5)
    @test x_last.step[] == x_last_b.step[]

    @i function loss(out, step, y, x, param, src, srcv, c; kwargs...)
        bennett((@skip! step), y, x, param, src, srcv, c; kwargs...)
        out += y.u[src...]
    end
    _,_,_,gx,_,_,gsrcv,gc = NiLang.AD.gradient(loss, (0.0, bennett_step!, zero(x), copy(x), param, src, srcv, copy(c)); iloss=1, k=k, N=N)
    x_last_2 = NiLang.direct_emulate(bennett_step!, (x2=copy(x); x2.u[src...]+=1e-5; x2), param, src, srcv, copy(c); N=N)
    @test gx.u[src...] ≈ (x_last_2.u[src...] - x_last.u[src...])/1e-5
    x_last_3 = NiLang.direct_emulate(bennett_step!, copy(x), param, src, srcv, (c2=copy(c); c2[src...]+=1e-5; c2); N=N)
    @test isapprox(gc[src...], (x_last_3.u[src...] - x_last.u[src...])/1e-5; rtol=1e-3)

    gn = ReversibleSeismic.SeismicState(Float64, nx, ny)
    gn.u[src...] = 1
    gx2, gv, gc2 = bennett_solve(x, gn; param=param, src, srcv=srcv, c=c, N=N, k=k)
    @test gx.u ≈ gx2.u
    @test gc ≈ gc2
end


text:
function generate_data(n, seed, density)
  function random_conflicts()
    E = Vector{Set}()
    D = 0
    while D < density
      D = 2*length(E)/(n*(n-1))
      new_edge = Set(rand(1:n,2))
      if !∈(new_edge,E) && length(new_edge) == 2
        push!(E,new_edge)
      end
    end
    return E
  end

  Random.seed!(seed)
  to_print = []
  push!(to_print,"InstanceID")
  push!(to_print,"$(seed)_$(n)")
  push!(to_print,"NumObjects")
  push!(to_print,n)
  push!(to_print,"Probabilities")
  p = rand(0.5:0.01:1.0,n)
  for i in 1:n
    push!(to_print,p[i])
  end
  push!(to_print,"Rewards")
  c = rand(1:1:1000,n)
  for i in 1:n
    push!(to_print,-c[i])
  end
  push!(to_print,"Conflicts")
  conflicts = random_conflicts()
  for (i,j) in conflicts
    push!(to_print,"$i $j")
  end
  writedlm(string(dirname(@__FILE__),"/../data/Smuggler/$(seed)_$(n)_$(density).txt"),to_print)
end

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

function read_data(instance::AbstractString)
  file = readdlm(string(dirname(@__FILE__),"/../data/Smuggler/",instance))
  @info "read instance $(file[2,1])"
  data = Data()
  data.constraint_uncertainty = false
  data.n = file[4,1]
  line = findfirst(file.=="Probabilities")[1]
  data.p = file[line+1:line+data.n,1]
  line = findfirst(file.=="Rewards")[1]
  data.hc = file[line+1:line+data.n,1]
  line = findfirst(file.=="Conflicts")[1]
  data.conflicts = Set.([file[i,1:2] for i in (line+1):size(file)[1]])
  data.compatibility = trues(data.n,data.n)
  for (i,j) in data.conflicts
    data.compatibility[i,j] = false
    data.compatibility[j,i] = false
  end
  data.E = Vector{Vector{Int64}}(undef, data.n)
  for i in 1:data.n
    data.E[i] = findall(x->x!=1,data.compatibility[i,:])
  end
  return data
end

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

mutable struct Data
  n::Int
  M::Int
  E::Vector{Vector{Int64}}
  hc::Vector{Float64}           # nominal rewards
  p::Vector{Float64}            # probabilities
  conflicts::Vector{Set}
  compatibility::Matrix{Bool}
  constraint_uncertainty::Bool
  nX_before_dom::Int
  heuristic_ub::Float64
  Data() = new()
end

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

function h(x)
  return prod(data.p[l] for l in x)
end

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

function generate_solutions(ub)
  function gen_sol(x::BitArray,i::Int64,ℓ::Int64,admissible::BitArray)
    time_elapsed = time() - starting_time
    if time_elapsed < TIME_LIM
      selected = findall(x)
      c = sum((1+ρ)*data.hc[selected])
      p = prod(data.p[selected])
      if p*c ≤ ub # weak and quick optimality check
        worstcase = create_model(TIME_LIM)
        @timeit to "build model" begin
          @variable(worstcase, (1+ρ)*data.hc[i] ≤ ξ[i in selected] ≤ (1-ρ)*data.hc[i] )
          @constraint(worstcase, sum( ((ξ[i] - data.hc[i])/data.hc[i])^2 for i in selected) ≤ Ω)
          @objective(worstcase, Min, sum(ξ[i] for i in selected)) # we are being optimistic here
        end
        @timeit to "optimize model" optimize!(worstcase)
        if p*objective_value(worstcase) ≤ ub # filter by optimality
          push!(X,copy(x))
        end
      end
      if i < data.n
        i = i+1
        x[i] = false
        gen_sol(x,i,ℓ,admissible)
        if ℓ < B
          x[i] = true
          if admissible[i] # filter by conflicts constraints
            ℓ += 1
            p *= data.p[i]
            c += (1+ρ)*data.hc[i]
            next_admissible = copy(admissible)
            next_admissible[i+1:end] = admissible[i+1:end].* data.compatibility[i,i+1:end]
            gen_sol(x,i,ℓ,next_admissible)
          end
        end
      end
    end
  end
  starting_time = time()
  X = Vector{BitArray{1}}() # set of all feasible vectors for the knapsack
  x = falses(data.n)
  i = 0
  ℓ = 0
  admissible = trues(data.n)
  println("Computing X")
  gen_sol(x,i,ℓ,admissible)
  println("|X|=$(length(X))")
  data.nX_before_dom = length(X)

  # Next we perform dominance over a full ellipsoid set because the comparison can be made much faster
  non_dom = trues(length(X))
  @timeit to "dominance" begin
    p = zeros(length(X))
    μ = Vector{Dict{Int64,Float64}}(undef,length(X))
    μ_square = Vector{Dict{Int64,Float64}}(undef,length(X))
    sum_μ = Vector{Float64}(undef,length(X))
    sqrtΩ = sqrt(Ω)
    for i in 1:length(X)
      μ[i] = Dict{Int64,Float64}()
      μ_square[i] = Dict{Int64,Float64}()
      x_i = findall(X[i])
      p[i] = prod(data.p[x_i])
      for l in x_i
        μ[i][l] = p[i]*data.hc[l]
        μ_square[i][l] = (μ[i][l])^2
      end
      sum_μ[i] = sum(μ[i][l] for l in x_i)
    end
    for i in 1:length(X)
      x_i = findall(X[i])
      for j in 1:length(X)
        sumsquares = 0
        @timeit to "get xj" x_j = findall(X[j])
        @timeit to "define sumsquare" begin
          for l in x_j
            if X[i] == 1
              sumsquares += (μ[i][l] - μ[j][l])^2
            else
              sumsquares += μ_square[j][l]
            end
          end
          for l in x_i
            if X[j] == 0
              sumsquares += μ_square[i][l]
            end
          end
        end
        # now we maximize μ^T ξ over the ellipsoid, performaning the change of variable η_i = (ξ_i - ĉ_i)/ĉ_i
        constant_term = sum_μ[j] - sum_μ[i]
        norm_term = sqrtΩ*sqrt(sumsquares)
        # We remove solution i if there exists a solution j that is cheaper than solution i over the full elliposid.
        if constant_term + norm_term < - ϵ
          non_dom[i] = false
          break
        end
      end
    end
  end
  X = X[findall(non_dom)]
  println("After dominance check: |X|=$(length(X))")
  return X
end

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

function compute_basic_ub(x)
  worstcase = create_model(TIME_LIM)
  selected = findall(x .> 0.99)
  @variable(worstcase, (1+ρ)*data.hc[i] ≤ ξ[i in selected] ≤ (1-ρ)*data.hc[i] )
  @constraint(worstcase, sum( ((ξ[i] - data.hc[i])/data.hc[i])^2 for i in selected) ≤ Ω )
  @objective(worstcase, Max, sum(ξ[i] for i in selected))
  optimize!(worstcase)
  cost = objective_value(worstcase)
  for i in selected
    cost *= data.p[i]
  end
  return cost
end

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

function generate_neighbors(sol,aa)
  X = Vector{BitArray{1}}()
  for k in 1:size(sol,1)
    x = falses(data.n)
    for j in N
      if sol[k,j]>0 x[j]=true end
    end
    push!(X,x)
    sol_items=findall(x)
    for i in sol_items
      admissible=trues(data.n)
      others = setdiff(sol_items,i)
      for o in others
        admissible[findall(x->x==false,data.compatibility[o,:])].=false
      end
      for iprime in setdiff(findall(admissible),sol_items)
        y = copy(x)
        y[i]=false
        y[iprime]=true
        push!(X,y)
      end
    end
  end
  return X
end

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

"Add the constraints forming X"
function add_constraints_X(model::Model, x::Array{VariableRef,1})
  @constraint(model,[i in 1:data.n],length(data.E[i])*x[i]+sum(x[j] for j in data.E[i])≤ length(data.E[i]))
  @constraint(model,sum(x[i] for i in 1:data.n) ≤ B)
end

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

"Create aa and cc with k rows. Notice aa is empty for this application."
function populate_scenarios()
  cc = Array{Float64,2}(undef,1,data.n)
  aa = Vector{Array{Float64,2}}(undef,1)
  aa[1] = Array{Float64,2}(undef,1,data.n)
  for l in N cc[1,l] = data.hc[l] end
  # for i = 1:k-1
  #   new_cost = [data.hc[l] for l in N] #This might be improved by adding other scenarios
  #   cc = [cc; new_cost']
  # end
  # Notice the first argument returned is aa, which is empty for this application
  return [], cc
end

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

"This function is non-trivial when the uncertainty set is defined through an extended formulation."
function add_scenario(ξ)
  return ξ, []
end

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

function build_and_solve_separation(x,v)
  time_remaining = TIME_LIM - (time()-STARTTIME)
  # Compute the worst-case cost:
  # max{η : ξ \in Ξ, η ≤ f(ξ,x^k), ∀ k in [k]}
  worstcase = create_model(TIME_LIM)
  @variable(worstcase, η)
  @variable(worstcase, (1+ρ)*data.hc[i] ≤ ξ[i in N] ≤ (1-ρ)*data.hc[i])
  @constraint(worstcase, sum(((ξ[i] - data.hc[i])/data.hc[i])^2 for i in N) ≤ Ω )
  @constraint(worstcase, sep_cost[k in K], η ≤ sum(ξ[i]*x[k,i] for i in N)*prod(data.p[i]^x[k,i] for i in N))
  @objective(worstcase, Max, η)
  @timeit to "separation" optimize!(worstcase)
  if termination_status(worstcase) == MOI.OPTIMAL
    return value.(ξ),objective_value(worstcase)
  else
    return [],-Inf
  end
end

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

"The heuristic constraint generation requires a starting upper bound. We provide a simple greedy bound instead
of the true static solution, which would require to solve a MINLP."

function solve_static()
  function callback(cb_data)
      selected = findall([ callback_value(cb_data,x[i]) > 0.9 for i in 1:data.n])
      # compute the worst ξ ∈ Ξ
      worstcase = create_model(TIME_LIM)
      @variable(worstcase, (1+ρ)*data.hc[i] ≤ ξ[i in selected] ≤ (1-ρ)*data.hc[i] )
      @constraint(worstcase, sum( ((ξ[i] - data.hc[i])/data.hc[i])^2 for i in selected) ≤ Ω )
      @objective(worstcase, Max, sum(ξ[i] for i in selected))
      optimize!(worstcase)

      #@warn worstcase

      #@warn "$(objective_value(worstcase)) ? > ? $(callback_value(cb_data,ω) + ϵ)"

      if objective_value(worstcase) > callback_value(cb_data,ω) + ϵ
        # add a cut related to an optimal scenario ξ
        con = @build_constraint( ω ≥ sum(value(ξ[i]) * x[i] for i in selected) + sum(data.hc[i] * x[i] for i in setdiff(1:data.n,selected)) )
        MOI.submit(static, MOI.LazyConstraint(cb_data), con)
        ncuts += 1
      end
   end

  ncuts = 0
  static = create_model(TIME_LIM)
  @variable(static, x[i in N], Bin,container=Array)
  @variable(static, sum(data.hc) ≤ ω ≤ 0 )
  add_constraints_X(static,x)
  @objective(static, Min, ω)
  MOI.set(static, MOI.LazyConstraintCallback(), callback)
  optimize!(static)
  probability = prod(data.p[findall(value.(x) .> 0.9)])
  #@warn probability
  @info "static value is $(probability*objective_value(static)), obtained by generating $ncuts cuts"

  "
  selected = []
  for index in sortperm(data.hc)
    length(selected) == L && break
    # check if index′ compatible with all selected items
    if all([ data.compatibility[index, index′] for index′ in selected ])
      push!(selected,index)
    end
  end

  probability = prod(data.p[selected])

  # compute the worst ξ ∈ Ξ
  worstcase = create_model(TIME_LIM)
  @variable(worstcase, (1-ρ)*data.hc[i] ≤ ξ[i in selected] ≤ (1+ρ)*data.hc[i] )
  @constraint(worstcase, sum( (ξ[i] - data.hc[i])^2 for i in selected) ≤ Ω )
  @objective(worstcase, Max, sum(ξ[i] for i in selected))
  optimize!(worstcase)"

  return probability*objective_value(static)
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
using CLIMAParameters.Planet: cv_d, T_0

export InitStateBC

export AtmosBC,
    Impenetrable,
    FreeSlip,
    NoSlip,
    DragLaw,
    Insulating,
    PrescribedTemperature,
    PrescribedEnergyFlux,
    BulkFormulaEnergy,
    Impermeable,
    OutflowPrecipitation,
    ImpermeableTracer,
    PrescribedMoistureFlux,
    BulkFormulaMoisture,
    PrescribedTracerFlux

export average_density_sfc_int

"""
    AtmosBC(momentum = Impenetrable(FreeSlip())
            energy   = Insulating()
            moisture = Impermeable()
            precipitation = OutflowPrecipitation()
            tracer  = ImpermeableTracer())

The standard boundary condition for [`AtmosModel`](@ref). The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct AtmosBC{M, E, Q, P, TR, TC}
    momentum::M = Impenetrable(FreeSlip())
    energy::E = Insulating()
    moisture::Q = Impermeable()
    precipitation::P = OutflowPrecipitation()
    tracer::TR = ImpermeableTracer()
    turbconv::TC = NoTurbConvBC()
end

boundary_conditions(atmos::AtmosModel) = atmos.problem.boundaryconditions

function boundary_state!(
    nf,
    bc::AtmosBC,
    atmos::AtmosModel,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    t,
    args...,
)
    atmos_boundary_state!(
        nf,
        bc,
        atmos,
        state⁺,
        aux⁺,
        n,
        state⁻,
        aux⁻,
        t,
        args...,
    )
    # update moisture auxiliary variables (perform saturation adjustment, if necessary)
    # to make thermodynamic quantities consistent with the boundary state
    atmos_nodal_update_auxiliary_state!(atmos.moisture, atmos, state⁺, aux⁺, t)
end

function boundary_state!(
    nf::Union{CentralNumericalFluxHigherOrder, CentralNumericalFluxDivergence},
    bc::AtmosBC,
    m::AtmosModel,
    x...,
)
    nothing
end

function atmos_boundary_state!(nf, bc::AtmosBC, atmos, args...)
    atmos_momentum_boundary_state!(nf, bc.momentum, atmos, args...)
    atmos_energy_boundary_state!(nf, bc.energy, atmos, args...)
    atmos_moisture_boundary_state!(nf, bc.moisture, atmos, args...)
    atmos_precipitation_boundary_state!(nf, bc.precipitation, atmos, args...)
    atmos_tracer_boundary_state!(nf, bc.tracer, atmos, args...)
    turbconv_boundary_state!(nf, bc.turbconv, atmos, args...)
end

function normal_boundary_flux_second_order!(
    nf,
    bc::AtmosBC,
    atmos::AtmosModel,
    fluxᵀn::Vars{S},
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    t,
    args...,
) where {S}
    atmos_normal_boundary_flux_second_order!(
        nf,
        bc,
        atmos,
        fluxᵀn,
        n⁻,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        t,
        args...,
    )
end

function atmos_normal_boundary_flux_second_order!(
    nf,
    bc::AtmosBC,
    atmos::AtmosModel,
    args...,
)
    atmos_momentum_normal_boundary_flux_second_order!(
        nf,
        bc.momentum,
        atmos,
        args...,
    )
    atmos_energy_normal_boundary_flux_second_order!(
        nf,
        bc.energy,
        atmos,
        args...,
    )
    atmos_moisture_normal_boundary_flux_second_order!(
        nf,
        bc.moisture,
        atmos,
        args...,
    )
    atmos_precipitation_normal_boundary_flux_second_order!(
        nf,
        bc.precipitation,
        atmos,
        args...,
    )
    atmos_tracer_normal_boundary_flux_second_order!(
        nf,
        bc.tracer,
        atmos,
        args...,
    )
    turbconv_normal_boundary_flux_second_order!(nf, bc.turbconv, atmos, args...)
end

"""
    average_density(ρ_sfc, ρ_int)

Average density between the surface and the interior point, given
 - `ρ_sfc` density at the surface
 - `ρ_int` density at the interior point
"""
function average_density(ρ_sfc::FT, ρ_int::FT) where {FT <: Real}
    return FT(0.5) * (ρ_sfc + ρ_int)
end

include("bc_momentum.jl")
include("bc_energy.jl")
include("bc_moisture.jl")
include("bc_precipitation.jl")
include("bc_initstate.jl")
include("bc_tracer.jl")

max_forks_repo_name: jm-c/CLIMA


text:
module MomentArithmetic

using IntervalArithmetic, ProbabilityBoundsAnalysis

import ProbabilityBoundsAnalysis: plot, makepbox
import IntervalArithmetic: issubset, intersect

using ProbabilityBoundsAnalysis: issubset

import Base.+
import Base.-
import Base.*
import Base./
import Base.^

import Base.issubset, Base.intersect
import Base: exp, sqrt, ^, log, +, -, *, /, abs, min, max, inv

global Nsub = [20]        # Number of sub-intervals

abstract type AbstractMoment end

export

    # From Moments.jl
    Moments, intervalM,
    make_consistent,
    makepbox, plot,

    # From Arithmetic.jl
    rowe, roweNoSub, rowevar, rowevarNoSub,
    +, -, *, /, reciprocal, exp, ^, sqrt, ln, log,
    abs, inv,
    sumIndep, subIndep, multIndep, divIndep, minIndep, maxIndep,
    sumFrechet, subFrechet, multFrechet, divFrechet, minFrechet,
    maxFrechet,

    sumPerfect, subPerfect, sumOpposite, subOpposite,

    sumCor, subCor, multCor, divCor, sumCov, subCov, multCov,

    # From Utils.jl
    env, left, right, split, issubset, intersect, hasIntMean, hasIntVar, hasIntMoments,

    area, dist,

    cov,

    # From IntervalArithmetic.jl
    interval, Interval,

    # From ProbabilityBoundsAnalysis.jl
    pbox

include("Moments.jl")
include("Utils.jl")
include("Arithmetic.jl")

end


text:
mutable struct Engine
    position_iterations
    velocity_iterations
    constraint_iterations
    timestep
    timestamp
    world
end

function Engine()
    position_iterations = 6
    velocity_iterations = 4
    constraint_iterations = 2
    world = World()
    return Engine(
        position_iterations,
        velocity_iterations,
        constraint_iterations,
        1,
        0,
        world)
end

# TODO
function update!(engine, dt)

    world = engine.world
    broadphase = engine.broadphase

    # update timestamp
    engine.timestamp += dt * engine.timestep

    # apply gravity to all bodies in world
    apply!(world.bodies, world.gravity)

    # ...

    # narrowphase pass
    collisions = detect!(broadphase, broadphase_pairs, engine)
    update!(engine.pairs, collisions, engine.timestamp)
    remove_old!(engine.pairs, engine.timestamp)

    # ...

    clear_forces!(engine.bodies)
end

function _clear_forces!(bodies)
    for body in bodies
        body.force = Vector(0., 0.)
        body.torque = 0.
    end
end

function _update!(bodies, dt)
    for body in bodies
        update!(body, dt)
    end
end


text:
"""
The macro `@iteratelist` writes a `for` loop over the `list` but avoiding `getfield` calls
inside the loop. The first argument of the macro is the `PairwiseListMatrix` that is going
to be iterated and the second is the body of the loop.
In the body `list` will be the list field of the `PairwiseListMatrix` and `k` the index
over that list. Other variables are expanded in the current scope. You must not modify the
value of `k`.

```jldoctest
julia> using PairwiseListMatrices

julia> PLM = PairwiseListMatrix([1,2,3], false)
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 0  1  2
 1  0  3
 2  3  0

julia> @iteratelist PLM println(list[k])
1
2
3

```
"""
macro iteratelist(plm, exp)
    esc(quote
        let list = $plm.list
            for k in 1:length(list)
                $exp
            end
        end
    end)
end

"""
The macro `@iteratediag` writes a `for` loop over the `diag` field of a
`PairwiseListMatrix{T,false,VT}` but avoiding calls to `getfield` inside the loop. The first
argument of the macro is the `PairwiseListMatrix` that is going to be iterated and the
second is the body of the loop. In the body `diag` will be the diag field of the
`PairwiseListMatrix` and `k` the index over that vector.
Other variables are expanded in the current scope. You must not modify the value of `k`.

```jldoctest
julia> using PairwiseListMatrices

julia> PLM = PairwiseListMatrix([1,2,3], false)
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 0  1  2
 1  0  3
 2  3  0

julia> @iteratediag PLM diag[k] += 10k

julia> PLM
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 10   1   2
  1  20   3
  2   3  30

```
"""
macro iteratediag(plm, exp)
    esc(quote
        if !$hasdiagonal($plm)
            local diag = $plm.diag
            for k in 1:length(diag)
                $exp
            end
        end
    end)
end

"""
The macro `@iterateupper` iterates over the upper triangular part of the
`PairwiseListMatrix` that is given as first argument. The second argument should be
`true` if the diagonal need to be included in the iteration or `false` otherwise.
The last argument is the body of the loop, where `list` is the list and diag fields of
the `PairwiseListMatrix` and `k` is the index over that `list`.
You can also use the respective `i` and `j` indexes for that position `k` in the upper
triangular part of the matrix. Other variables are expanded in the current scope.
You must not modify the values of `i`, `j` or `k`.

```jldoctest
julia> using PairwiseListMatrices

julia> PLM = PairwiseListMatrix([1,2,3], true)
2×2 PairwiseListMatrix{Int64,true,Array{Int64,1}}:
 1  2
 2  3

julia> mat = zeros(Int, 2, 2)
2×2 Array{Int64,2}:
 0  0
 0  0

julia> let mat = mat # To avoid using global
           @iterateupper PLM true mat[i,j] = list[k]
       end

julia> mat
2×2 Array{Int64,2}:
 1  2
 0  3

```
"""
macro iterateupper(plm, use_diag, exp)
    esc(quote
        let N = $plm.nelements
            local k, diag, list
            if $hasdiagonal($plm)
                if $use_diag
                    k = 0
                    list = $plm.list
                    for i in 1:N
                        for j in i:N
                            k += 1
                            $exp
                        end
                    end
                else
                    k = 0
                    list = $plm.list
                    for i in 1:N
                        for j in i:N
                            k += 1
                            if i != j
                                $exp
                            end
                        end
                    end
                end
            else
                if $use_diag
                    k = 0
                    diag = $plm.diag
                    list = $plm.list
                    for i in 1:N
                        for j in i:N
                            if i != j
                                k += 1
                                $exp
                            else
                                let list = diag, k = i
                                    $exp
                                end
                            end
                        end
                    end
                else
                    k = 0
                    list = $plm.list
                    for i in 1:(N-1)
                        for j in (i+1):N
                            k += 1
                            $exp
                        end
                    end
                end
            end
        end
    end)
end


text:
"""
    invariantset(embedding, tolerance)

Determine whether an embedding forms an invariant set.
"""
function invariantset(embedding)
    lastpoint = embedding[end, :]
    dim = size(embedding, 2)

    # Triangulate the embedding using all points but the last
    triangulation = Simplices.Delaunay.delaunayn(embedding[1:end-1, :])

    points = embedding[1:end-1, :]
    simplex_indices = triangulation

    # Centroids and radii of simplices in the triangulation
    centroids, radii = centroids_radii2(points, simplex_indices)

    lastpoint_matrix = repmat(lastpoint', size(centroids, 1), 1)

    # Find simplices that can contain the last point (not all can)
    dists_lastpoint_and_centroids = sum((lastpoint_matrix - centroids).^2, 2)
    distdifferences = radii.^2 - dists_lastpoint_and_centroids

    # Find the row indices of the simplices that possibly contain the last point (meaning that
    # dist(simplex_i, lastpoint) <= radius(simplex), so the circumsphere of the simplex
    # contains the last point.
    valid_simplex_indices = find(heaviside0(distdifferences) .* collect(1:size(triangulation[2], 1)))

    n_validsimplices = size(valid_simplex_indices, 1)

    # Loop over valid simplices and check whether the corresponding simplex actually
    # contains the last point.
    i = 1
    lastpoint_contained = false

    # Continue checking if simplices contain the last point until some simplex contains it.
    # Then the point must necessarily be inside the convex hull of the triangulation formed
    # by those simplices.
    while i <= n_validsimplices && !lastpoint_contained
        # Subsample the valid simplex and compute its orientation
        simplex = embedding[simplex_indices[valid_simplex_indices[i], :], :]
        orientation_simplex = det([ones(dim + 1, 1) simplex])

        beta = 1 # Convex expansion coefficient

        j = 1
        while j <= dim + 1 && beta >= 0
            tmp = copy(simplex)
            tmp[j, :] = copy(lastpoint)
            beta = det([ones(dim + 1, 1) tmp]) * sign(orientation_simplex)
            j = j + 1
        end

        # If the last convex expansion coefficient is positive, the last point is contained
        # in the triangulation (because all the previous coefficients must have been nonnegative)
        if beta >= 0
            lastpoint_contained = true
        end

        i = i + 1
    end
    # If the last point is contained in the triangulation, the set is invariant.
    return lastpoint_contained
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
function train!(model::Network, x, y)
    values = []
    J = @diff loss(model, x, y; loss_values=values)
    for par in params(model)
        g = grad(J, par)
        update!(value(par), g, par.opt)
    end
    return (sum(values), values[1], values[2])
end

function epoch!(model::Network, data)
    atype = typeof(value(model.qnetwork.mu_layer.w))
    Lx = Lz = 0.0f0
    iter = 0
    @time for (x, y) in data
        J1, J2 = train!(model, atype(reshape(x, 784, size(x,4))), y)
        Lx += J1
        Lz += J2
        iter += 1
    end
    lossval = Lx+Lz
    return lossval/iter, Lx/iter, Lz/iter
end

function validate(model::Network, data)
    atype = typeof(value(model.qnetwork.mu_layer.w))
    Lx = Lz = 0.0f0
    iter = 0
    values = []
    for (x, y) in data
        loss(model, atype(reshape(x, 784, size(x,4))), y; loss_values=values)
        J1, J2 = values
        Lx += J1
        Lz += J2
        iter += 1
        empty!(values)
    end
    lossval = Lx + Lz
    return lossval/iter, Lx/iter, Lz/iter
end

function parse_options(args)
    s = ArgParseSettings()
    s.description = "DRAW model on MNIST."

    @add_arg_table s begin
        # ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}");
        #  help="array and float type to use")
        ("--batchsize"; arg_type=Int; default=64; help="batch size")
        ("--zdim"; arg_type=Int; default=10; help="noise dimension")
        ("--encoder"; arg_type=Int; default=256; help="hidden units")
        ("--decoder"; arg_type=Int; default=256; help="hidden units")
        ("--embed"; arg_type=Int; default=100; help="condition dimension")
        ("--epochs"; arg_type=Int; default=20; help="# of training epochs")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
        ("--gridsize"; arg_type=Int; nargs=2; default=[9,9])
        ("--gridscale"; arg_type=Float32; default=2.0f0)
        ("--optim"; default="Adam(;beta1=0.5f0, gclip=5.0f0)")
        ("--loadfile"; default=nothing; help="file to load trained models")
        ("--outdir"; default=nothing; help="output dir for models/generations")
        ("--A"; arg_type=Int; default=28)
        ("--B"; arg_type=Int; default=28)
        ("--N"; arg_type=Int; default=5)
        ("--T"; arg_type=Int; default=10)
        ("--valid"; action=:store_true)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    # o[:atype] = eval(parse(o[:atype]))
    if o[:outdir] != nothing
        o[:outdir] = abspath(o[:outdir])
    end
    return o
end

function train(args)
    println("script started"); flush(stdout)
    o = parse_options(args)
    o[:seed] > 0 && Knet.seed!(o[:seed])
    println("options parsed"); flush(stdout)
    display(o)

    model = Network(
        o[:A], o[:B], o[:N], o[:T], o[:encoder],
        o[:decoder], o[:zdim]; embed=o[:embed])
    println("model initialized"); flush(stdout)
    init_opt!(model, o[:optim])
    println("optimization parameters initiazlied"); flush(stdout)
    dtrn, dtst = mnistdata(xtype=Array{Float32})
    println("data loaded"); flush(stdout)

    bestloss = Inf
    for epoch = 1:o[:epochs]
        trnloss = epoch!(model, dtrn)
        tstloss = o[:valid] ? validate(model, dtst) : (.0f0, .0f0, .0f0)
        report(epoch, trnloss, tstloss)
        if tstloss[1] < bestloss
            bestloss = tstloss[1]
        end
    end
    return model
end

max_stars_repo_name: ilkerkesen/draw

max_issues_repo_name: ilkerkesen/draw

max_forks_repo_name: ilkerkesen/draw


text:
# Random Fourier Features
# from: Rahimi,Recht, "Random Features for Large-Scale Kernel Machines", NIPS 2007
#
# Danny Perry (dperry@cs.utah.edu)
# Oct 2014

module RandomFourierFeatures

export RandomFourierFunctions, GenerateFunctionsGaussian, GenerateFunctionsGaussianSinCos, GenerateFunctionsGaussianBiased, RFFProject, RFFProjectSinCos, FastFood, GenerateFastFood, FFProject

include("GaussianKernel.jl")

# A set of random fourier functions
type RandomMaclaurinFunctions
	directions::Array{Float64,2} # random direction
	offsets::Array{Float64,1} # random shift in [0,2*pi]
end

# project the input data into the approximate kernel space
#
# @param X - input data
# @param RFFs - the random fourier functions
function RMFProject(X::Array{Float64,2}, RFFs::RandomMaclaurinFunctions)
	n = size(X,1)
	z = sqrt(2) * cos( (X * RFFs.directions') .+ ones(n,1) * RFFs.offsets' ) * (1/sqrt(length(RFFs.offsets)))
end

include("PolynomialKernel.jl")

end # module