text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -4.847398418942631 <= q <= 3.5461678827819525)

                     add_NL_constraint(m, :(log(1 + exp(-0.19745953189550702 + -0.698583212677272*log(1 + exp(-0.2024495957518173 + -0.1501654471281335*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + 0.7965619538518913*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.8656395810315112*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))) + 0.07925279857906231*log(1 + exp(-0.4774715210476401 + -0.18966514690787006*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + 0.1687755600264067*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.9882336248352352*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))) + -0.8658353105986021*log(1 + exp(-0.5967554189497073 + 0.7632066186026765*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + -0.767873838532454*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.02410968604346264*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + log(1 + exp(0.9047173820312868 + 0.09375520961098927*log(1 + exp(-0.2024495957518173 + -0.1501654471281335*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + 0.7965619538518913*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.8656395810315112*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))) + -0.9282166334997988*log(1 + exp(-0.4774715210476401 + -0.18966514690787006*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + 0.1687755600264067*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.9882336248352352*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))) + 0.3928411969784107*log(1 + exp(-0.5967554189497073 + 0.7632066186026765*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + -0.767873838532454*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.02410968604346264*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) + log(1 + exp(0.16193188153327842 + 0.4964335766578394*log(1 + exp(-0.2024495957518173 + -0.1501654471281335*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + 0.7965619538518913*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.8656395810315112*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))) + -0.6165233713312808*log(1 + exp(-0.4774715210476401 + -0.18966514690787006*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + 0.1687755600264067*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.9882336248352352*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))) + -0.4911670927879337*log(1 + exp(-0.5967554189497073 + 0.7632066186026765*log(1 + exp(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))) + -0.767873838532454*log(1 + exp(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))) + 0.02410968604346264*log(1 + exp(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module ImageMetadata

# The order here is designed to avoid an ambiguity warning in convert,
# see the top of ImageAxes
using ImageAxes
using ImageCore
using ColorVectorSpace   # for overriding math operations with Gray/RGB
import AxisArrays

import Base: +, -, *, /
import Base: permutedims

const ViewIndex = Union{Base.ViewIndex, Colon}

export
    # types
    ImageMeta,

    # functions
    copyproperties,
    arraydata,
    properties,
    shareproperties,
    spatialproperties

#### types and constructors ####

# Concrete types
"""
`ImageMeta` is an AbstractArray that can have metadata, stored in a dictionary.

Construct an image with `ImageMeta(A, props)` (for a properties dictionary
`props`), or with `ImageMeta(A, prop1=val1, prop2=val2, ...)`.
"""
struct ImageMeta{T,N,A<:AbstractArray,P<:AbstractDict{Symbol,Any}} <: AbstractArray{T,N}
    data::A
    properties::P
end

function ImageMeta(data::AbstractArray{T,N}, props::AbstractDict{Symbol,Any}) where {T,N}
    return ImageMeta{T,N,typeof(data),typeof(props)}(data,props)
end
function ImageMeta(data::AbstractArray{T,N}, props::Dict{Symbol,Any}) where {T,N}
    return ImageMeta{T,N,typeof(data),typeof(props)}(data,props)
end
function ImageMeta(data::AbstractArray, props::Dict{Symbol})
    return ImageMeta(data, convert(Dict{Symbol,Any}, props))
end
function ImageMeta(data::AbstractArray; kwargs...)
    return ImageMeta(data, kwargs2dict(kwargs))
end
ImageMeta(data::AbstractArray, props::AbstractDict) = throw(ArgumentError("properties must be an AbstractDict{Symbol,Any}"))

const ImageMetaArray{T,N,A<:Array} = ImageMeta{T,N,A}
const ImageMetaAxis{T,N,A<:AxisArray} = ImageMeta{T,N,A}

Base.size(A::ImageMeta) = size(arraydata(A))
Base.size(A::ImageMetaAxis, Ax::Axis) = size(arraydata(A), Ax)
Base.size(A::ImageMetaAxis, ::Type{Ax}) where {Ax<:Axis} = size(arraydata(A), Ax)
Base.axes(A::ImageMeta) = axes(arraydata(A))
Base.axes(A::ImageMetaAxis, Ax::Axis) = axes(arraydata(A), Ax)
Base.axes(A::ImageMetaAxis, ::Type{Ax}) where {Ax<:Axis} = axes(arraydata(A), Ax)

datatype(::Type{ImageMeta{T,N,A,P}}) where {T,N,A<:AbstractArray,P} = A

Base.IndexStyle(::Type{M}) where {M<:ImageMeta} = IndexStyle(datatype(M))

AxisArrays.HasAxes(A::ImageMetaAxis) = AxisArrays.HasAxes{true}()

# getindex and setindex!
for AType in (ImageMeta, ImageMetaAxis)
    @eval begin
        @inline function Base.getindex(img::$AType{T,1}, i::Int) where T
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds ret = arraydata(img)[i]
            ret
        end
        @inline function Base.getindex(img::$AType, i::Int)
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds ret = arraydata(img)[i]
            ret
        end
        @inline function Base.getindex(img::$AType{T,N}, I::Vararg{Int,N}) where {T,N}
            @boundscheck checkbounds(arraydata(img), I...)
            @inbounds ret = arraydata(img)[I...]
            ret
        end

        @inline function Base.setindex!(img::$AType{T,1}, val, i::Int) where T
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds arraydata(img)[i] = val
            val
        end
        @inline function Base.setindex!(img::$AType, val, i::Int)
            @boundscheck checkbounds(arraydata(img), i)
            @inbounds arraydata(img)[i] = val
            val
        end
        @inline function Base.setindex!(img::$AType{T,N}, val, I::Vararg{Int,N}) where {T,N}
            @boundscheck checkbounds(arraydata(img), I...)
            @inbounds arraydata(img)[I...] = val
            val
        end
    end
end

@inline function Base.getindex(img::ImageMetaAxis, ax::Axis, I...)
    result = arraydata(img)[ax, I...]
    maybe_wrap(img, result)
end
@inline function Base.getindex(img::ImageMetaAxis, i::Union{Integer,AbstractVector,Colon}, I...)
    result = arraydata(img)[i, I...]
    maybe_wrap(img, result)
end
maybe_wrap(img::ImageMeta{T}, result::T) where T = result
maybe_wrap(img::ImageMeta{T}, result::AbstractArray{T}) where T = copyproperties(img, result)

@inline function Base.setindex!(img::ImageMetaAxis, val, ax::Axis, I...)
    setindex!(arraydata(img), val, ax, I...)
end
@inline function Base.setindex!(img::ImageMetaAxis, val, i::Union{Integer,AbstractVector,Colon}, I...)
    setindex!(arraydata(img), val, i, I...)
end

Base.view(img::ImageMeta, ax::Axis, I...) = shareproperties(img, view(arraydata(img), ax, I...))
Base.view(img::ImageMeta{T,N}, I::Vararg{ViewIndex,N}) where {T,N} = shareproperties(img, view(arraydata(img), I...))
Base.view(img::ImageMeta, i::ViewIndex) = shareproperties(img, view(arraydata(img), i))
Base.view(img::ImageMeta, I::Vararg{ViewIndex,N}) where {N} = shareproperties(img, view(arraydata(img), I...))

function Base.getproperty(img::ImageMeta, propname::Symbol)
    # TODO remove these once deprecations are done
    if propname === :data
        Base.depwarn("img.data is deprecated use arraydata(img) to extract data img.", :arraydata)
        return arraydata(img)
    elseif propname === :properties
        Base.depwarn("img.properties is deprecated use properties(img) to extract data img.", :properties)
        return properties(img)
    else
        return properties(img)[propname]
    end
end

function Base.setproperty!(img::ImageMeta, propname::Symbol, val)
    # TODO remove these once deprecations are done
    if propname === :data
        Base.depwarn("Replacing the data field with is deprecated. Consider creating a new instance of ImageMeta instead.", :rawdata)
        setindex!(img, val, img)
    elseif propname === :properties
        Base.depwarn("Replacing the properties field is deprecated. Consider creating a new instance of ImageMeta instead.", :properties)
        for (k,v) in val
            setproperty!(img, k, v)
        end
    else
        properties(img)[propname] = val
    end
    return img
end

Base.propertynames(img::ImageMeta) = (keys(properties(img))...,)

Base.copy(img::ImageMeta) = ImageMeta(copy(arraydata(img)), deepcopy(properties(img)))

Base.convert(::Type{ImageMeta}, A::ImageMeta) = A
Base.convert(::Type{ImageMeta}, A::AbstractArray) = ImageMeta(A)
Base.convert(::Type{ImageMeta{T}}, A::ImageMeta{T}) where {T} = A
Base.convert(::Type{ImageMeta{T}}, A::ImageMeta) where {T} = shareproperties(A, convert(AbstractArray{T}, arraydata(A)))
Base.convert(::Type{ImageMeta{T}}, A::AbstractArray{T}) where {T} = ImageMeta(A)
Base.convert(::Type{ImageMeta{T}}, A::AbstractArray) where {T} = ImageMeta(convert(AbstractArray{T}, A))

# copy properties
function Base.copy!(imgdest::ImageMeta, imgsrc::ImageMeta, prop1::Symbol, props::Symbol...)
    setproperty!(imgdest, prop1, deepcopy(getproperty(imgsrc, prop1)))
    for p in props
        setproperty!(imgdest, p, deepcopy(getproperty(imgsrc, p)))
    end
    return imgdest
end

# similar
Base.similar(img::ImageMeta, ::Type{T}, shape::Dims) where {T} = ImageMeta(similar(arraydata(img), T, shape), copy(properties(img)))
Base.similar(img::ImageMetaAxis, ::Type{T}) where {T} = ImageMeta(similar(arraydata(img), T), copy(properties(img)))

"""
    copyproperties(img::ImageMeta, data) -> imgnew

Create a new "image," copying the properties dictionary of `img` but
using the data of the AbstractArray `data`. Note that changing the
properties of `imgnew` does not affect the properties of `img`.

See also: [`shareproperties`](@ref).
"""
copyproperties(img::ImageMeta, data::AbstractArray) =
    ImageMeta(data, copy(properties(img)))

"""
    shareproperties(img::ImageMeta, data) -> imgnew

Create a new "image," reusing the properties dictionary of `img` but
using the data of the AbstractArray `data`. The two images have
synchronized properties; modifying one also affects the other.

See also: [`copyproperties`](@ref).
"""
shareproperties(img::ImageMeta, data::AbstractArray) = ImageMeta(data, properties(img))

# Delete a property!
Base.delete!(img::ImageMeta, propname::Symbol) = delete!(properties(img), propname)

# Iteration
# Defer to the array object in case it has special iteration defined
Base.iterate(img::ImageMeta) = Base.iterate(arraydata(img))
Base.iterate(img::ImageMeta, s) = Base.iterate(arraydata(img), s)

# Show
const emptyset = Set()
function showim(io::IO, img::ImageMeta)
    IT = typeof(img)
    print(io, eltype(img).name.name, " ImageMeta with:\n  data: ", summary(arraydata(img)), "\n  properties:")
    showdictlines(io, properties(img), get(img, :suppress, emptyset))
end
Base.show(io::IO, img::ImageMeta) = showim(io, img)
Base.show(io::IO, ::MIME"text/plain", img::ImageMeta) = showim(io, img)

function Base.reinterpret(::Type{T}, img::ImageMeta) where {T}
    shareproperties(img, reinterpret(T, arraydata(img)))
end

"""
    arraydata(img::ImageMeta) -> array

Extract the data from `img`, omitting the properties
dictionary. `array` shares storage with `img`, so changes to one
affect the other.

See also: [`properties`](@ref).
"""
ImageAxes.arraydata(img::ImageMeta) = getfield(img, :data)

function ImageCore.permuteddimsview(A::ImageMeta, perm)
    ip = sortperm([perm...][[coords_spatial(A)...]])  # the inverse spatial permutation
    permutedims_props!(copyproperties(A, permuteddimsview(arraydata(A), perm)), ip)
end
ImageCore.channelview(A::ImageMeta) = shareproperties(A, channelview(arraydata(A)))
ImageCore.colorview(::Type{C}, A::ImageMeta{T,N}) where {C<:Colorant,T,N} = shareproperties(A, colorview(C, arraydata(A)))
ImageCore.colorview(::Type{ARGB32}, A::ImageMeta{T,N}) where {T,N} = shareproperties(A, colorview(ARGB32, arraydata(A)))
ImageCore.rawview(A::ImageMeta{T}) where {T<:Real} = shareproperties(A, rawview(arraydata(A)))
ImageCore.normedview(::Type{T}, A::ImageMeta{S}) where {T<:FixedPoint,S<:Unsigned} = shareproperties(A, normedview(T, arraydata(A)))

# AxisArrays functions
AxisArrays.axes(img::ImageMetaAxis) = AxisArrays.axes(arraydata(img))
AxisArrays.axes(img::ImageMetaAxis, d::Int) = AxisArrays.axes(arraydata(img), d)
AxisArrays.axes(img::ImageMetaAxis, Ax::Axis) = AxisArrays.axes(arraydata(img), Ax)
AxisArrays.axes(img::ImageMetaAxis, ::Type{Ax}) where {Ax<:Axis} = AxisArrays.axes(arraydata(img), Ax)
AxisArrays.axisdim(img::ImageMetaAxis, ax) = axisdim(arraydata(img), ax)
AxisArrays.axisnames(img::ImageMetaAxis) = axisnames(arraydata(img))
AxisArrays.axisvalues(img::ImageMetaAxis) = axisvalues(arraydata(img))

#### Properties ####

"""
    properties(imgmeta) -> props

Extract the properties dictionary `props` for `imgmeta`. `props`
shares storage with `img`, so changes to one affect the other.

See also: [`arraydata`](@ref).
"""
properties(img::ImageMeta) = getfield(img, :properties)

if isdefined(Base, :hasproperty) # or VERSION >= v"1.2"
    import Base: hasproperty
end
hasproperty(img::ImageMeta, k::Symbol) = haskey(properties(img), k)

Base.get(img::ImageMeta, k::Symbol, default) = get(properties(img), k, default)

# So that defaults don't have to be evaluated unless they are needed,
# we also define a @get macro (thanks Toivo Hennington):
struct IMNothing end   # to avoid confusion in the case where dict[key] === nothing
macro get(img, k, default)
    quote
        img, k = $(esc(img)), $(esc(k))
        val = get(properties(img), k, IMNothing())
        return isa(val, IMNothing) ? $(esc(default)) : val
    end
end

ImageAxes.timeaxis(img::ImageMetaAxis) = timeaxis(arraydata(img))
ImageAxes.timedim(img::ImageMetaAxis) = timedim(arraydata(img))
ImageAxes.colordim(img::ImageMetaAxis) = colordim(arraydata(img))

ImageCore.pixelspacing(img::ImageMeta) = pixelspacing(arraydata(img))

"""
    spacedirections(img)

Using ImageMetadata, you can set this property manually. For example, you
could indicate that a photograph was taken with the camera tilted
30-degree relative to vertical using

```
img["spacedirections"] = ((0.866025,-0.5),(0.5,0.866025))
```

If not specified, it will be computed from `pixelspacing(img)`, placing the
spacing along the "diagonal".  If desired, you can set this property in terms of
physical units, and each axis can have distinct units.
"""
ImageCore.spacedirections(img::ImageMeta) = @get img :spacedirections spacedirections(arraydata(img))

ImageCore.sdims(img::ImageMetaAxis) = sdims(arraydata(img))

ImageCore.coords_spatial(img::ImageMetaAxis) = coords_spatial(arraydata(img))

ImageCore.spatialorder(img::ImageMetaAxis) = spatialorder(arraydata(img))

ImageAxes.nimages(img::ImageMetaAxis) = nimages(arraydata(img))

ImageCore.size_spatial(img::ImageMetaAxis) = size_spatial(arraydata(img))

ImageCore.indices_spatial(img::ImageMetaAxis) = indices_spatial(arraydata(img))

ImageCore.assert_timedim_last(img::ImageMetaAxis) = assert_timedim_last(arraydata(img))

#### Permutations over dimensions ####

"""
    permutedims(img, perm, [spatialprops])

When permuting the dimensions of an ImageMeta, you can optionally
specify that certain properties are spatial and they will also be
permuted. `spatialprops` defaults to `spatialproperties(img)`.
"""
permutedims

function permutedims_props!(ret::ImageMeta, ip, spatialprops=spatialproperties(ret))
    if !isempty(spatialprops)
        for prop in spatialprops
            if hasproperty(ret, prop)
                a = getproperty(ret, prop)
                if isa(a, AbstractVector)
                    setproperty!(ret, prop, a[ip])
                elseif isa(a, Tuple)
                    setproperty!(ret, prop, a[ip])
                elseif isa(a, AbstractMatrix) && size(a,1) == size(a,2)
                    setproperty!(ret, prop, a[ip,ip])
                else
                    error("Do not know how to handle property ", prop)
                end
            end
        end
    end
    ret
end

function permutedims(img::ImageMetaAxis, perm)
    p = AxisArrays.permutation(perm, axisnames(arraydata(img)))
    ip = sortperm([p...][[coords_spatial(img)...]])
    permutedims_props!(copyproperties(img, permutedims(arraydata(img), p)), ip)
end
function permutedims(img::ImageMeta, perm)
    ip = sortperm([perm...][[coords_spatial(img)...]])
    permutedims_props!(copyproperties(img, permutedims(arraydata(img), perm)), ip)
end

# Note: `adjoint` does not recurse into ImageMeta properties.
function Base.adjoint(img::ImageMeta{T,2}) where {T<:Real}
    ip = sortperm([2,1][[coords_spatial(img)...]])
    permutedims_props!(copyproperties(img, adjoint(arraydata(img))), ip)
end

function Base.adjoint(img::ImageMeta{T,1}) where T<:Real
    check_empty_spatialproperties(img)
    copyproperties(img, arraydata(img)')
end

"""
    spatialproperties(img)

Return a vector of strings, containing the names of properties that
have been declared "spatial" and hence should be permuted when calling
`permutedims`.  Declare such properties like this:

    img[:spatialproperties] = [:spacedirections]
"""
spatialproperties(img::ImageMeta) = @get img :spatialproperties [:spacedirections]

function check_empty_spatialproperties(img)
    sp = spatialproperties(img)
    for prop in sp
        if hasproperty(img, prop)
            error("spatialproperties must be empty, have $prop")
        end
    end
    nothing
end

#### Low-level utilities ####
function showdictlines(io::IO, dict::AbstractDict, suppress::Set)
    for (k, v) in dict
        if k === :suppress
            continue
        end
        if !in(k, suppress)
            print(io, "\n    ", k, ": ")
            print(IOContext(io, :compact => true), v)
        else
            print(io, "\n    ", k, ": <suppressed>")
        end
    end
end
showdictlines(io::IO, dict::AbstractDict, prop::Symbol) = showdictlines(io, dict, Set([prop]))

# printdictval(io::IO, v) = print(io, v)
# function printdictval(io::IO, v::Vector)
#     for i = 1:length(v)
#         print(io, " ", v[i])
#     end
# end

# converts keyword argument to a dictionary
function kwargs2dict(kwargs)
    d = Dict{Symbol,Any}()
    for (k, v) in kwargs
        d[k] = v
    end
    return d
end

include("operators.jl")
include("deprecations.jl")

end # module

ext: jl


text:
@testset "Strings" begin
    @test LO(sqrt(d64"2")) < LO(nextfloat(sqrt(d64"2")))
end

@testset "Adjacent values" begin
    @test nextfloat(Double64(10))  > Double64(10)
    @test nextfloat(Double32(-10)) > Double32(-10)
    @test prevfloat(Double64(10))  < Double64(10)
    @test prevfloat(Double32(-10)) < Double32(-10)
    @test nextfloat(Double64(10),3)  > nextfloat(Double64(10),2)
    @test nextfloat(Double32(-10),3) > nextfloat(Double32(-10),2)
    @test prevfloat(Double64(10),3)  < prevfloat(Double64(10),2)
    @test prevfloat(Double32(-10),3) < prevfloat(Double32(-10),2)
end

@testset "Exponential functions"  begin
    @test Double64(10.0)^0 == Double64(1.0)
    @test Double64(10.0)^1 == Double64(10.0)
    @test_throws DomainError Double64(0.0)^0
end


text:
include("_common.jl")

using CategoricalArrays, StatsPlots

climate_fig_path = joinpath(FIG_PATH, "climate") * "/"
mkpath(climate_fig_path)

# Change plot theme for publication
theme(:default, grid=false)

dataset = copy(FULL_DATASET)
data_years = Dates.year.(dataset.Date)
start_year = minimum(data_years)
end_year = maximum(data_years)
half_decades = cut(data_years, start_year:5:end_year, extend=true)
annual = cut(data_years, start_year:1:end_year, extend=true)

dataset.half_decade_set = half_decades
dataset.annual = annual

half_decade_dataset = groupby(dataset, :half_decade_set)
annual_dataset = groupby(dataset, :annual)

# Plot historic rainfall
thresholds = quantile(FULL_DATASET[:, "410730_P"], [0.6, 0.95])
bar(FULL_DATASET.Date, FULL_DATASET[:, "410730_P"],
	title="Observed Rainfall",
	label="",
	xlabel="Year",
	ylabel="Rainfall [mm]")
hline!(thresholds, label="Bin thresholds")
savefig("$(climate_fig_path)rainfall_record.png")

rain_sym = Symbol("410730_P")
temp_sym = Symbol("410730_max_T")
flow_sym = Symbol("410730_Q")

climate_annual_stats = combine(annual_dataset,
						rain_sym => mean,
						rain_sym => std,
						rain_sym => maximum,
						rain_sym => sum,
						temp_sym => mean,
						temp_sym => std,
						temp_sym => maximum,
						temp_sym => sum)

# Plot annual rain
bar(start_year:end_year, climate_annual_stats[:, "410730_P_sum"],
	title="Observed Annual Rainfall",
	label="",
	xlabel="Year",
	ylabel="Rainfall [mm]")
savefig("$(climate_fig_path)rainfall.png")

climate_half_decade_stats = combine(half_decade_dataset,
							rain_sym => mean,
							rain_sym => std,
							rain_sym => maximum,
							rain_sym => sum)

flow_annual_stats = combine(annual_dataset,
						flow_sym => mean,
						flow_sym => std,
						flow_sym => maximum,
						flow_sym => sum)

flow_half_decade_stats = combine(half_decade_dataset,
							flow_sym => mean,
							flow_sym => std,
							flow_sym => maximum,
							flow_sym => sum)

# flow_half_decade_stats

bar(start_year+1:end_year, climate_annual_stats[:, "410730_P_sum"], label="",
	title="Total Annual Rainfall",
	xlabel="Year",
	ylabel="Rainfall [mm]")
savefig("$(climate_fig_path)rainfall_sum.png")

p_total = climate_annual_stats[:, "410730_P_sum"]
violin(["Calibration"], p_total[1:30], label="", title="Annual Rainfall", ylabel="Total [mm]")
violin!(["Validation"], p_total[31:end], label="")
dotplot!(["Calibration"], p_total[1:30],
		 marker=(:black, stroke(0)),
		 label="")
dotplot!(["Validation"], p_total[31:end],
		 marker=(:black, stroke(0)),
		 label="")
savefig("$(climate_fig_path)rainfall_comparison.png")

# Average temperatures
plot(start_year+1:end_year, climate_annual_stats[:, "410730_max_T_mean"], color="red",
	 label="",
	 title="Mean Annual Temperature",
	 xlabel="Year",
	 ylabel="Temperature [°C]")
savefig("$(climate_fig_path)average_max_T.png")

T_total = climate_annual_stats[:, "410730_max_T_mean"]
violin(["Calibration"], T_total[1:30], label="", title="Mean Max Temperature", ylabel="Degrees [°C]")
violin!(["Validation"], T_total[31:end], label="")
dotplot!(["Calibration"], T_total[1:30],
		 marker=(:black, stroke(0)),
		 label="")
dotplot!(["Validation"], T_total[31:end],
		 marker=(:black, stroke(0)),
		 label="")
savefig("$(climate_fig_path)temperature_comparison.png")

# Coefficient of variation
climate_annual_stats.rainfall_variability = climate_annual_stats[:, "410730_P_std"] ./ climate_annual_stats[:, "410730_P_mean"]

p_vari = climate_annual_stats.rainfall_variability
violin(["Calibration"], p_vari[1:30], label="", title="Rainfall Variability", ylabel="Variability [CV]")
violin!(["Validation"], p_vari[31:end], label="")
dotplot!(["Calibration"], p_vari[1:30],
		 marker=(:black, stroke(0)),
		 label="")
dotplot!(["Validation"], p_vari[31:end],
		 marker=(:black, stroke(0)),
		 label="")
savefig("$(climate_fig_path)rainfall_variability_comparison.png")

bar(start_year:end_year, flow_annual_stats[:, "410730_Q_sum"] ./ 1000,
		xlabel="Year", ylabel="Annual Flow [GL]", label="",
		title="Annual Flow [GL]"
)

savefig("$(climate_fig_path)annual_flow.png")

flow_annual_stats.flow_variability = (flow_annual_stats[:, "410730_Q_std"] ./ flow_annual_stats[:, "410730_Q_mean"])

bar(start_year:end_year, flow_annual_stats.flow_variability,
	xlabel="Year", ylabel="Flow variability [CV]", label="",
	title="Variability of Annual Streamflow"
)

# Qualitatively determined separation between calib start/end
c_end = start_year + 30
v_start = start_year + 31
vline!([v_start], label="")
annotate!((start_year+0.5, 1.6, text("Calibration", font(11), :left)))
annotate!((v_start+0.5, 1.6, text("Validation", font(11), :left)))
savefig("$(climate_fig_path)flow_variability.png")

flow_col = Symbol("flow_variability")
violin(["Calibration"], flow_annual_stats[1:30, flow_col],
		label="", ylabel="Variability of Flow [CV]",
		title="Annual Flow\nVariability")
violin!(["Validation"], flow_annual_stats[31:end, flow_col],
		label="")
dotplot!(["Calibration"], flow_annual_stats[1:30, flow_col],
		 marker=(:black, stroke(0)),
		 label="")
dotplot!(["Validation"], flow_annual_stats[31:end, flow_col],
		 marker=(:black, stroke(0)),
		 label="")
savefig("$(climate_fig_path)flow_variability_comparison.png")

"""
Köppen aridity index

See:
https://doi.org/10.1016/j.palaeo.2013.05.008

AIₖ : arid: < 40; humid: 40–160; perhumid: > 160
"""
AI_k(P, T) = P ./ (T .+ 33.0)

AP = climate_annual_stats[:, "410730_P_sum"]
AT = climate_annual_stats[:, "410730_max_T_mean"]
aridity = AI_k(AP, AT)

bar(start_year:end_year, aridity,
	label="",
	xlabel="Year",
	ylabel="\$AI_{k}\$",
	title="Köppen Aridity Index"
)

vline!([v_start], label="", linewidth=2.0)
annotate!((start_year, 34.0, text("Calibration", font(11), :left)))
annotate!((v_start+0.5, 34.0, text("Validation", font(11), :left)))
savefig("$(climate_fig_path)aridity.png")

violin(["Calibration"], aridity[1:30], label="",
		ylabel="\$AI_{k}\$",
		title="Aridity Characteristics")
violin!(["Validation"], aridity[31:end], label="")
dotplot!(["Calibration"], aridity[1:30],
		 marker=(:black, stroke(0)),
		 label="")
dotplot!(["Validation"], aridity[31:end],
		 marker=(:black, stroke(0)),
		 label="")
savefig("$(climate_fig_path)aridity_comparison.png")

@info "Timeframe:" start_year, end_year
@info minimum(dataset.Date), maximum(dataset.Date)

flow_col = Symbol("410730_Q_sum")
violin(["Calibration"], flow_annual_stats[1:30, flow_col] ./ 1000,
		label="", ylabel="Annual Flow [GL]",
		title="Annual Flow\nCharacteristics")
violin!(["Validation"], flow_annual_stats[31:end, flow_col] ./ 1000,
		label="")
dotplot!(["Calibration"], flow_annual_stats[1:30, flow_col] ./1000,
		 marker=(:black, stroke(0)),
		 label="")
dotplot!(["Validation"], flow_annual_stats[31:end, flow_col] ./1000,
		 marker=(:black, stroke(0)),
		 label="")
savefig("$(climate_fig_path)flow_comparison.png")

describe(flow_annual_stats[1:30, Symbol("410730_Q_sum")] ./ 1000)
# Summary Stats:
# Length:         30
# Missing Count:  0
# Mean:           48.433578
# Minimum:        6.393224
# 1st Quartile:   31.157348
# Median:         46.762521
# 3rd Quartile:   62.423117
# Maximum:        110.813315
# Type:           Float64

describe(flow_annual_stats[31:end, Symbol("410730_Q_sum")] ./ 1000)
# Summary Stats:
# Length:         21
# Missing Count:  0
# Mean:           35.346256
# Minimum:        8.673252
# 1st Quartile:   21.041727
# Median:         35.154885
# 3rd Quartile:   51.299297
# Maximum:        67.310798
# Type:           Float64

# Plot streamflow
plot(CALIB_DATES, FULL_DATASET[1:CALIB_LN, "410730_Q"], label="Calibration", color="darkblue", title="Historic Streamflow")
plot!(VALID_DATES, FULL_DATASET[CALIB_LN+1:end, "410730_Q"], label="Validation", color="lightblue", ylabel="Streamflow [ML]")
savefig("$(climate_fig_path)historic_flow.png")


text:
using Test
using JavisGraphs
using Graphs
using GraphPlot
using LaTeXStrings

include("graph.jl")
include("vertex.jl")
include("edge.jl")

@testset "Constructions" begin
    @test create_weighted_graph()
    @test create_reference_graph()
    @test check_graph_build()
    @test check_vertex_build()
    @test check_edge_build()
end

function create_dummy_graph()
    g = create_graph()
    v1 = GraphVertex(1, 1:2)
    v2 = GraphVertex(2, 1:2)
    e1 = GraphEdge(1, 2, 1:2)
    e2 = GraphEdge(2, 1, 1:2)
    return g
end

include("utils.jl")
@testset "Utils" begin
    @test check_get_edges()
    @test check_get_edges_unordered()
    @test check_get_edges_ordered()
    @test check_get_vertices()
    @test check_get_draw()
    @test check_graph_from_matrix()
    check_vertex_out_of_parent_range()
    check_edge_out_of_parent_range()
end

include("render.jl")
@testset "Render" begin
    @test check_graph_render()
end

@testset "Layout" begin
    check__global_layout()
end

@testset "Style utils" begin
    @test check_add_styles()
    @test check_register_style_opts()
end

include("styles.jl")
@testset "vertex styles" begin
    @test check_vertex_shape()
    @test check_vertex_text()
    @test check_vertex_fill()
    @test check_vertex_border()
end

@testset "Edge styles" begin
    @test check_edge_shape()
    @test check_edge_style()
    @test check_edge_arrow()
    @test check_edge_text()
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c7356460-055b-11eb-0d69-dd42369dc913
try using AddPackage catch; using Pkg; Pkg.add("AddPackage"); using AddPackage end

# ╔═╡ df4bc620-055b-11eb-39c7-e952028567da
begin
	using Base.Iterators: repeated
	@add using Colors
	@add using CUDA
	@add using Images
	@add using MLDatasets
	@add using Parameters
	@add using Flux
	using Flux.Data: DataLoader
	using Flux: @epochs, onehotbatch, onecold, logitcrossentropy, throttle
	using Statistics
end

# ╔═╡ d7a1f38e-055b-11eb-31ef-55de28cda23d
md"# Classifier for MNIST"

# ╔═╡ d8d88ce0-055d-11eb-32e8-31586d97a59c
md"## Training arguments"

# ╔═╡ 07f1cf70-055c-11eb-05ab-e77200ea8ea7
@with_kw mutable struct Args
	α::Float64 = 3e-4      # learning rate
	batchsize::Int = 1024  # batch size
	epochs::Int = 20       # number of epochs
	device::Function = gpu # set as gpu if available
	throttle::Int = 1      # throttle print every X seconds
end

# ╔═╡ f9bad300-055d-11eb-1e4d-f362cf07c6b5
md"## Data processing"

# ╔═╡ 37042240-055c-11eb-3eba-0da0d64627aa
function getdata(args)
	# load dataset
	xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
	xtest, ytest = MLDatasets.MNIST.testdata(Float32)

	# reshape data to flatten each image into a vector
	xtrain = Flux.flatten(xtrain)
	xtest = Flux.flatten(xtest)

	# one-hot encode the labels
	ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

	# split into batches
	traindata = DataLoader(xtrain, ytrain, batchsize=args.batchsize, shuffle=true)
	testdata = DataLoader(xtest, ytest, batchsize=args.batchsize)

	return traindata, testdata
end

# ╔═╡ ffa3bbb0-055d-11eb-2652-fd3b56dcf474
md"## Model"

# ╔═╡ 87cbe9b0-055c-11eb-1c34-9d9bedbe4ad5
function buildmodel(; imgsize=(28,28,1), nclasses=10)
	return Chain(
			Dense(prod(imgsize), 32, relu),
			Dense(32, nclasses))
end

# ╔═╡ 054402a2-055e-11eb-273e-1340d52bbed8
md"""
## Cost function
$$\begin{gather}
\mathcal{L}(\hat{y}, y) = -\frac{1}{n}\sum_{i=1}^n y \left(\hat{y} - \log\left(\sum e^{\hat{y}}\right)\right)\\
\mathcal{J}(\mathbf{\hat{y}}, \mathbf{y}) = \frac{1}{m}\sum \mathcal{L}(\hat{y}, y)
\end{gather}$$
"""

# ╔═╡ b391a990-055c-11eb-13b8-b3f9ad78cea7
cost(dataloader, model) = mean(logitcrossentropy(model(x), y) for (x,y) in dataloader)

# ╔═╡ ea7f0630-055e-11eb-0717-bde12a519dfc
md"## Accuracy"

# ╔═╡ db2160e0-055c-11eb-0c66-6b3cc667e18e
function accuracy(dataloader, model)
	acc = 0
	for (x,y) in dataloader
		acc += sum(onecold(cpu(model(x))) .== onecold(cpu(y))) / size(x,2)
	end
	return acc/length(dataloader)
end

# ╔═╡ ef2b0210-055e-11eb-244b-d18b38d83a71
md"## Model training"

# ╔═╡ 01bccd20-055d-11eb-07a0-73fa53c45c24
function train(; kwargs...)
	# initialize model parameters
	args = Args(; kwargs...)

	# load data
	traindata, testdata = getdata(args)

	# construct model
	m = buildmodel()
	traindata = args.device.(traindata)
	testdata = args.device.(testdata)
	m = args.device(m)
	loss(x,y) = logitcrossentropy(m(x), y)

	# callback, optimizer, and training
	callback = throttle(()->@show(cost(traindata, m)), args.throttle)
	opt = ADAM(args.α)

	@epochs args.epochs Flux.train!(loss, params(m), traindata, opt, cb=callback)

	@show accuracy(traindata, m)
	@show accuracy(testdata, m)

	return m, traindata, testdata
end

# ╔═╡ 6092c160-055d-11eb-3364-0df28086629e
m, traindata, testdata = train();

# ╔═╡ f89b15b0-055e-11eb-3a07-091d5e96e340
begin
	Markdown.parse(string("\$\$\\operatorname{Loss}_\\text{train} = ", 100round(accuracy(traindata, m), digits=5), "\\%\$\$", "\$\$\\operatorname{Loss}_\\text{test} = ", 100round(accuracy(testdata, m), digits=5), "\\%\$\$"))
end

# ╔═╡ Cell order:
# ╟─d7a1f38e-055b-11eb-31ef-55de28cda23d
# ╠═c7356460-055b-11eb-0d69-dd42369dc913
# ╠═df4bc620-055b-11eb-39c7-e952028567da
# ╟─d8d88ce0-055d-11eb-32e8-31586d97a59c
# ╠═07f1cf70-055c-11eb-05ab-e77200ea8ea7
# ╟─f9bad300-055d-11eb-1e4d-f362cf07c6b5
# ╠═37042240-055c-11eb-3eba-0da0d64627aa
# ╟─ffa3bbb0-055d-11eb-2652-fd3b56dcf474
# ╠═87cbe9b0-055c-11eb-1c34-9d9bedbe4ad5
# ╟─054402a2-055e-11eb-273e-1340d52bbed8
# ╠═b391a990-055c-11eb-13b8-b3f9ad78cea7
# ╟─ea7f0630-055e-11eb-0717-bde12a519dfc
# ╠═db2160e0-055c-11eb-0c66-6b3cc667e18e
# ╟─ef2b0210-055e-11eb-244b-d18b38d83a71
# ╠═01bccd20-055d-11eb-07a0-73fa53c45c24
# ╠═6092c160-055d-11eb-3364-0df28086629e
# ╟─f89b15b0-055e-11eb-3a07-091d5e96e340


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook <: SwaggerModel
    admissionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: admissionReviewVersions
    clientConfig::Any # spec type: Union{ Nothing, IoK8sApiAdmissionregistrationV1beta1WebhookClientConfig } # spec name: clientConfig
    failurePolicy::Any # spec type: Union{ Nothing, String } # spec name: failurePolicy
    matchPolicy::Any # spec type: Union{ Nothing, String } # spec name: matchPolicy
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespaceSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: namespaceSelector
    objectSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: objectSelector
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAdmissionregistrationV1beta1RuleWithOperations} } # spec name: rules
    sideEffects::Any # spec type: Union{ Nothing, String } # spec name: sideEffects
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook(;admissionReviewVersions=nothing, clientConfig=nothing, failurePolicy=nothing, matchPolicy=nothing, name=nothing, namespaceSelector=nothing, objectSelector=nothing, rules=nothing, sideEffects=nothing, timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("admissionReviewVersions"), admissionReviewVersions)
        setfield!(o, Symbol("admissionReviewVersions"), admissionReviewVersions)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("clientConfig"), clientConfig)
        setfield!(o, Symbol("clientConfig"), clientConfig)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("failurePolicy"), failurePolicy)
        setfield!(o, Symbol("failurePolicy"), failurePolicy)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("matchPolicy"), matchPolicy)
        setfield!(o, Symbol("matchPolicy"), matchPolicy)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("namespaceSelector"), namespaceSelector)
        setfield!(o, Symbol("namespaceSelector"), namespaceSelector)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("objectSelector"), objectSelector)
        setfield!(o, Symbol("objectSelector"), objectSelector)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("sideEffects"), sideEffects)
        setfield!(o, Symbol("sideEffects"), sideEffects)
        validate_property(IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook

const _property_map_IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook = Dict{Symbol,Symbol}(Symbol("admissionReviewVersions")=>Symbol("admissionReviewVersions"), Symbol("clientConfig")=>Symbol("clientConfig"), Symbol("failurePolicy")=>Symbol("failurePolicy"), Symbol("matchPolicy")=>Symbol("matchPolicy"), Symbol("name")=>Symbol("name"), Symbol("namespaceSelector")=>Symbol("namespaceSelector"), Symbol("objectSelector")=>Symbol("objectSelector"), Symbol("rules")=>Symbol("rules"), Symbol("sideEffects")=>Symbol("sideEffects"), Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook = Dict{Symbol,String}(Symbol("admissionReviewVersions")=>"Vector{String}", Symbol("clientConfig")=>"IoK8sApiAdmissionregistrationV1beta1WebhookClientConfig", Symbol("failurePolicy")=>"String", Symbol("matchPolicy")=>"String", Symbol("name")=>"String", Symbol("namespaceSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("objectSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("rules")=>"Vector{IoK8sApiAdmissionregistrationV1beta1RuleWithOperations}", Symbol("sideEffects")=>"String", Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook }) = collect(keys(_property_map_IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook))
Swagger.property_type(::Type{ IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook[name]))}
Swagger.field_name(::Type{ IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook }, property_name::Symbol) =  _property_map_IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook[property_name]

function check_required(o::IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook)
    (getproperty(o, Symbol("clientConfig")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAdmissionregistrationV1beta1ValidatingWebhook }, name::Symbol, val)
end


text:
"""
    config_diagnostics(driver_config, timeend; interval=nothing)

Returns the state and tendency diagnostic groups
"""
function config_diagnostics(driver_config, timeend; interval = nothing)
    FT = eltype(driver_config.grid)
    info = driver_config.config_info
    if interval == nothing
        interval = "$(cld(timeend, 2) + 10)ssecs"
        #interval = "10steps"
    end

    boundaries = [
        FT(0) FT(0) FT(0)
        FT(info.hmax) FT(info.hmax) FT(info.zmax)
    ]
    axes = (
        [FT(1)],
        [FT(1)],
        collect(range(boundaries[1, 3], boundaries[2, 3], step = FT(50)),),
    )
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries;
        axes = axes,
    )
    ds_dgngrp = setup_dump_state_diagnostics(
        SingleStackConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )
    dt_dgngrp = setup_dump_tendencies_diagnostics(
        SingleStackConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )
    return ClimateMachine.DiagnosticsConfiguration([ds_dgngrp, dt_dgngrp])
end


text:
mutable struct ScoreContext{P <: Parameters} <: ExecutionContext
    select::ConstrainedSelection
    weight::Float64
    visited::Visitor
    params::P
    function Score(obs::Vector{Tuple{K, P}}) where {P, K <: Union{Symbol, Pair}}
        c_sel = selection(obs)
        new{EmptyParameters}(c_sel, 0.0, Parameters())
    end
    ScoreContext(obs::K, params::P) where {K <: ConstrainedSelection, P <: Parameters} = new{P}(obs, 0.0, Visitor(), params)
end
Score(obs::Vector) = ScoreContext(selection(obs))
Score(obs::ConstrainedSelection) = ScoreContext(obs, Parameters())
Score(obs::ConstrainedSelection, params) = ScoreContext(obs, params)

# ------------ Convenience ------------ #

function score(sel::L, fn::Function, args...) where L <: ConstrainedSelection
    ctx = Score(sel)
    ret = ctx(fn, args...)
    b, missed = compare(sel.query, ctx.visited)
    b || error("ScoreError: did not visit all constraints in selection.\nDid not visit: $(missed).")
    return ret, ctx.weight
end

function score(sel::L, params, fn::Function, args...) where L <: ConstrainedSelection
    ctx = Score(sel, params)
    ret = ctx(fn, args...)
    b, missed = compare(sel.query, ctx.visited)
    b || error("ScoreError: did not visit all constraints in selection.\nDid not visit: $(missed).")
    return ret, ctx.weight
end

function score(sel::L, fn::typeof(rand), d::Distribution{K}) where {L <: ConstrainedSelection, K}
    ctx = Score(sel)
    addr = gensym()
    ret = ctx(fn, addr, d)
    b, missed = compare(sel.query, ctx.visited)
    b || error("ScoreError: did not visit all constraints in selection.\nDid not visit: $(missed).")
    return ret, ctx.weight
end

# TODO: fix for dispatch on params.
function score(sel::L, fn::typeof(markov), call::Function, len::Int, args...; params = Parameters()) where L <: ConstrainedSelection
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = Score(v_sel, params)
    ret = ctx(fn, addr, call, len, args...)
    b, missed = compare(sel.query, ctx.visited)
    b || error("ScoreError: did not visit all constraints in selection.\nDid not visit: $(missed).")
    return ret, ctx.weight
end

function score(sel::L, fn::typeof(plate), call::Function, args::Vector; params = Parameters()) where L <: ConstrainedSelection
    ctx = Score(sel, params)
    addr = gensym()
    ret = ctx(fn, addr, call, args)
    b, missed = compare(sel.query, ctx.visited)
    b || error("ScoreError: did not visit all constraints in selection.\nDid not visit: $(missed).")
    return ret, ctx.weight
end

function score(sel::L, fn::typeof(plate), d::Distribution{K}, len::Int; params = Parameters()) where {L <: ConstrainedSelection, K}
    addr = gensym()
    v_sel = selection(addr => sel)
    ctx = Score(v_sel, params)
    ret = ctx(fn, addr, d, len)
    b, missed = compare(sel.query, ctx.visited)
    b || error("ScoreError: did not visit all constraints in selection.\nDid not visit: $(missed).")
    return ret, ctx.weight
end

# ------------ includes ------------ #

include("hierarchical/score.jl")
include("plate/score.jl")
include("markov/score.jl")
include("factor/score.jl")

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct ScoreContext{P <: Parameters} <: ExecutionContext
    select::ConstrainedSelection
    weight::Float64
    params::P
end
```

The `ScoreContext` is used to score selections according to a model function. For computation in the `ScoreContext` to execute successfully, the `select` selection must provide a choice for every address visited in the model function, and the model function must allow the context to visit every constraints expressed in `select`.

Inner constructors:

```julia
function Score(obs::Vector{Tuple{K, P}}) where {P, K <: Union{Symbol, Pair}}
    c_sel = selection(obs)
    new{EmptyParameters}(c_sel, 0.0, Parameters())
end
```

Outer constructors:

```julia
ScoreContext(obs::K, params) where {K <: ConstrainedSelection} = new(obs, 0.0, params)
Score(obs::Vector) = ScoreContext(selection(obs))
Score(obs::ConstrainedSelection) = ScoreContext(obs, Parameters())
Score(obs::ConstrainedSelection, params) = ScoreContext(obs, params)
```
""", ScoreContext)

@doc(
"""
```julia
ret, w = score(sel::L, fn::Function, args...; params = Parameters()) where L <: ConstrainedSelection
ret, w = score(sel::L, fn::typeof(rand), d::Distribution{K}; params = Parameters()) where {L <: ConstrainedSelection, K}
ret, w = score(sel::L, fn::typeof(markov), call::Function, len::Int, args...; params = Parameters()) where L <: ConstrainedSelection
ret, w = score(sel::L, fn::typeof(plate), call::Function, args::Vector; params = Parameters()) where L <: ConstrainedSelection
ret, w = score(sel::L, fn::typeof(plate), d::Distribution{K}, len::Int; params = Parameters()) where {L <: ConstrainedSelection, K}
```

`score` provides an API to the `ScoreContext` execution context. You can use this function on any of the matching signatures above - it will return the return value `ret`, and the likelihood weight score of the user-provided selection `sel`. The selection should satisfy the following requirement:

1. At any random choice in any branch traveled according to the constraints of `sel`, `sel` must provide a constraint for that choice.

Simply put, this just means you need to provide a constraint for each `ChoiceSite` you encounter.
""", score)


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# This module introduces an N-dimensional iterator for use in for-loops with
# arbitrary number of variables. Each variable will be incremented by one.
#
# Currently, we implement an iterator interface for the iteration over the
# elements of the iterator, not to advance the iterator.
#
# TODO: Implement iterator interface defined by Julia documentation. The object
# that the iterator returns is the array of indices.

module NDIterators

export NDIterator
export advance

struct NDIterator{T} # <: AbstractArray{T,1}
    idx::Array{T,1}
    idxmin::Array{T,1}
    idxmax::Array{T,1}
    i::Array{Int,0}
end

function NDIterator(idxmin::AbstractArray, idxmax::AbstractArray)
    @assert length(idxmin) == length(idxmax)
    idx = deepcopy(idxmin)
    idx[end] -= 1
    ndi = NDIterator(idx, idxmin, idxmax, fill(length(idx)))
    for i=1:length(idx)
        if idxmax[i] < idxmin[i]
            ndi.i[] = 0
        end
    end
    return ndi
end
NDIterator(idxmin::Integer, idxmax::Integer; N=1) = NDIterator(fill(idxmin,N), fill(idxmax,N))
NDIterator(idxmin::Integer, idxmax) = NDIterator(fill(idxmin, length(idxmax)), convert(Array, idxmax))
NDIterator(idxmin, idxmax::Integer) = NDIterator(convert(Array, idxmin), fill(idxmax, length(idxmin)))
NDIterator(idxmin, idxmax) = NDIterator(convert(Array, idxmin), convert(Array, idxmax))

Base.getindex(ndi::NDIterator, i::Int) = ndi.idx[i]
Base.:-(ndi::NDIterator) = - ndi.idx
Base.convert(::Type{T}, ndi::NDIterator) where {T<:AbstractArray} = T(ndi.idx)

Base.firstindex(ndi::NDIterator) = 1
Base.lastindex(ndi::NDIterator) = length(ndi.idx)

Base.show(io::IO, ndi::NDIterator) = show(io, ndi.idx)

# for @__dot__ syntax
Base.length(ndi::NDIterator) = length(ndi.idx)
Base.iterate(ndi::NDIterator, i=1) = i > length(ndi) ? nothing : (ndi.idx[i], i+1)
Base.eltype(::NDIterator{T}) where {T} = T

function advance(ndi::NDIterator)
    i = ndi.i[]
    N = length(ndi.idx)
    while true
        if i <= 0
            ndi.i[] = i
            return false
        end
        ndi.idx[i] += 1
        if ndi.idx[i] > ndi.idxmax[i]
            ndi.idx[i] = ndi.idxmin[i] - 1
            i -= 1
        elseif i < N
            i += 1
        else
            break
        end
    end
    ndi.i[] = i
    return true
end

end

# vim: set sw=4 et sts=4 :


text:
@testset "thresholdexceedance.jl" begin
    n = 10000

    pd = GeneralizedPareto(1.0, 0.1)
    y = rand(pd, n)

    x = rand(n)

    ev = [Variable("x", x)]

    smodel = Extremes.ThresholdExceedance(Variable("y", y))
    nsmodel = Extremes.ThresholdExceedance(Variable("y", y), logscalecov = ev, shapecov = ev)

    @testset "ThresholdExceedance(exceedances; logscalecov, shapecov)" begin
        ev10 = Extremes.Variable("t", collect(1:n+10))

        # Build with logscalecov length != y length
        @test_throws AssertionError Extremes.ThresholdExceedance(Variable("y", y), logscalecov = [ev10])

        # Build with shapecov length != y length
        @test_throws AssertionError Extremes.ThresholdExceedance(Variable("y", y), shapecov = [ev10])

        # Build with all optional parameters set
        @test nsmodel.data.value == y
        @test nsmodel.logscale.covariate == ev
        @test nsmodel.shape.covariate == ev
    end

    @testset "paramindex(model)" begin
        # model with stationary and non-stationary parameters
        model = Extremes.ThresholdExceedance(Variable("y", y), logscalecov = ev)

        paramin = Extremes.paramindex(model)

        @test length(paramin) == 2
        @test paramin[:ϕ] == [1,2]
        @test paramin[:ξ] == [3]
    end

    @testset "getcovariatenumber(model)" begin
        # stationary
        @test Extremes.getcovariatenumber(smodel) == 0

        # non-stationary
        @test Extremes.getcovariatenumber(nsmodel) == 2
    end

    @testset "nparameter(model)" begin
        # stationary
        @test Extremes.nparameter(smodel) == 2

        # non-stationary
        @test Extremes.nparameter(nsmodel) == 4
    end

    @testset "getdistribution(model, θ)" begin
        # length(θ) != nparameter(model) throws
        @test_throws AssertionError Extremes.getdistribution(smodel, [1.0])

        # stationary
        n = 100

        σ = 1.0
        ξ = .1
        ϕ = log(σ)

        θ = [ϕ ; ξ]

        pd = GeneralizedPareto(σ, ξ)
        y = rand(pd,n)

        model = ThresholdExceedance(Variable("y", y))

        fd = Extremes.getdistribution(model, θ)[]

        @test fd == pd

        # non-stationary
        n = 100

        x₁ = randn(n)/3
        x₂ = randn(n)/3
        x₃ = randn(n)/10

        θ = [-.5 ; 1.0 ; 1.0 ; 0 ; 1.0]

        ϕ = θ[1] .+ θ[2] * x₁ .+ θ[3] * x₂
        ξ = θ[4] .+ θ[5] * x₃

        σ = exp.(ϕ)

        pd = GeneralizedPareto.(σ, ξ)
        y = rand.(pd)

        model = ThresholdExceedance(Variable("y", y), logscalecov = [Variable("x₁", x₁), Variable("x₂", x₂)], shapecov = [Variable("x₃", x₃)])

        fd = Extremes.getdistribution(model, θ)

        @test fd == pd

    end

    @testset "getinitialvalue(::Type{GeneralizedPareto},y)" begin
        # Test with valid pwm GPD estimates
        y = [0.0, 1.0, 2.0]
        ini = Extremes.getinitialvalue(GeneralizedPareto,y)
        @test !isapprox(ini[2],0)
        pd = GeneralizedPareto(exp(ini[1]), ini[2])
        @test all(insupport.(pd,y))

        # Test with invalid pwm GPD estimates
        y = [0.0 , 1.0, 1.0, 1.0, 3.0]
        ini = Extremes.getinitialvalue(GeneralizedPareto,y)
        @test ini[2] ≈ 0
        pd = GeneralizedPareto(exp(ini[1]), ini[2])
        @test all(insupport.(pd,y))

    end

    @testset "getinitialvalue(model)" begin
        # Test with valid pwm GPD estimates
        y = [0.0, 1.0, 2.0]
        model = ThresholdExceedance(Variable("y", y))
        ini = Extremes.getinitialvalue(model)
        @test ini[1] ≈ .-.693 atol = .001
        @test ini[2] ≈ .5 atol = .001

        # Test with invalid pwm GPD estimates
        y = [0.0 , 1.0, 1.0, 1.0, 3.0]
        model = ThresholdExceedance(Variable("y", y))
        ini = Extremes.getinitialvalue(model)
        @test ini[1] ≈ .182 atol = .001
        @test ini[2] ≈ 0 atol = .001

    end

    @testset "showEVA(io, obj; prefix)" begin
        # print does not throw
        buffer = IOBuffer()
        @test_logs Extremes.showEVA(buffer, smodel, prefix = "\t")

    end

end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export GallopingNeighborhood

mutable struct GallopingNeighborhood <: NeighborhoodAlgorithm
    increase_factor::Float32
    decrease_factor::Float32
    lower::Float32
    upper::Float32
    neighborhood::Int
end

function GallopingNeighborhood()
    return GallopingNeighborhood(1.5)
end

function GallopingNeighborhood(factor)
    return GallopingNeighborhood(factor, 1f0/factor, 0.95, 1.05, 1)
end

function optimize_neighborhood!(algo::GallopingNeighborhood, index::SearchGraph{T}, dist::Function, perf, recall) where {T}
    # restarts, beam_size, montecarlo_size = index.restarts, index.beam_size, index.montecarlo_size
    # index.restarts, index.beam_size, index.montecarlo_size = 1, 1, 1
    SearchType = typeof(index.search_algo)
    tmp, index.search_algo = index.search_algo, SearchType()
    p = probe(perf, index, dist)
    if p.recall <= recall * algo.lower
        # the recall is too low, it needs larger neighborhoods
        algo.neighborhood = algo.increase_factor * algo.neighborhood |> ceil |> Int
    elseif p.recall >= recall * algo.upper
        # the recall is too high, it needs smaller neighborhoods
        algo.neighborhood = algo.decrease_factor * algo.neighborhood |> ceil |> Int
    end
    n = length(index.db)
    algo.neighborhood = min(max(1, algo.neighborhood), ceil(Int, log2(n)^2))
    index.search_algo = tmp
    # index.restarts, index.beam_size, index.montecarlo_size = restarts, beam_size, montecarlo_size
end

function neighborhood(algo::GallopingNeighborhood, index::SearchGraph{T}, dist::Function, item::T) where {T}
    k = algo.neighborhood
    nbuffer = Vector{Int32}(undef, 0)
    knn = search(index, dist, item, KnnResult(k))

    for p in knn
        push!(nbuffer, p.objID)
    end

    return knn, nbuffer
end


text:
# Autogenerated wrapper script for fplll_jll for x86_64-linux-gnu
export fplll, latticegen, latsieve, libfplll

using CompilerSupportLibraries_jll
using GMP_jll
using MPFR_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `fplll`
const fplll_splitpath = ["bin", "fplll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
fplll_path = ""

# fplll-specific global declaration
function fplll(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(fplll_path)
    end
end

# Relative path to `latticegen`
const latticegen_splitpath = ["bin", "latticegen"]

# This will be filled out by __init__() for all products, as it must be done at runtime
latticegen_path = ""

# latticegen-specific global declaration
function latticegen(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(latticegen_path)
    end
end

# Relative path to `latsieve`
const latsieve_splitpath = ["bin", "latsieve"]

# This will be filled out by __init__() for all products, as it must be done at runtime
latsieve_path = ""

# latsieve-specific global declaration
function latsieve(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(latsieve_path)
    end
end

# Relative path to `libfplll`
const libfplll_splitpath = ["lib", "libfplll.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libfplll_path = ""

# libfplll-specific global declaration
# This will be filled out by __init__()
libfplll_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libfplll = "libfplll.so.6"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"fplll")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (CompilerSupportLibraries_jll.PATH_list, GMP_jll.PATH_list, MPFR_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (CompilerSupportLibraries_jll.LIBPATH_list, GMP_jll.LIBPATH_list, MPFR_jll.LIBPATH_list,))

    global fplll_path = normpath(joinpath(artifact_dir, fplll_splitpath...))

    push!(PATH_list, dirname(fplll_path))
    global latticegen_path = normpath(joinpath(artifact_dir, latticegen_splitpath...))

    push!(PATH_list, dirname(latticegen_path))
    global latsieve_path = normpath(joinpath(artifact_dir, latsieve_splitpath...))

    push!(PATH_list, dirname(latsieve_path))
    global libfplll_path = normpath(joinpath(artifact_dir, libfplll_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libfplll_handle = dlopen(libfplll_path)
    push!(LIBPATH_list, dirname(libfplll_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
using Printf

function generate_descriptor(Q, τ)
    τstring = @sprintf("%.2f", τ)
    Qstring = string(round(Int, Q))
    descriptor = "_τ" * τstring * "_Q" * Qstring
    return descriptor
end

function channel_looper(Q, τ)
    # Architecture
    CUDA.allowscalar(true)
    arch = GPU()
    FT   = Float64

    # File IO
    ic_load = true

    write_slices = false
    write_zonal  = true
    slice_output_interval = 48hour
    zonal_output_interval = 5*365day

    time_avg_window =  zonal_output_interval / 2  # needs to be a float
    checkpoint_interval = 365 * 10 *  day

    resolution = 16
    descriptor = string(resolution)
    descriptor2 = generate_descriptor(Q, τ)
    location = "/storage6/MesoscaleRuns/"
    filename = location * "Abernathy_" * descriptor * descriptor2

    if ic_load
        loadfilename = "Abernathy_" * "16"
        filepath = pwd() * "/" * getlatest(loadfilename)
    end
    ## Domain
    Lx = 1000.0kilometer
    Ly = 2000.0kilometer
    Lz = 2985
    # Discretization
    maxΔt = 1100.0 * 16 / resolution  # [s]
    Δt =  ic_load ? maxΔt : Δt = 300.0 * 16 / resolution # [s]

    end_time = 100 * 365day
    advection   = WENO5()
    timestepper = :RungeKutta3
    # Rough target resolution
    Δx = Δy = 5kilometer # 5km
    Δz = 100meter
    # Multiple of 16 gridpoints
    Nx = round(Int, Lx / Δx / 16) * resolution
    Ny = round(Int, Ly / Δy / 16) * resolution
    Nz = round(Int, Lz / Δz / 16) * 16
    # Create Grid
    topology = (Periodic, Bounded, Bounded)
    grid = RegularCartesianGrid(topology=topology,
                                size=(Nx, Ny, Nz),
                                x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

    # Parameters
    f = -1e-4
    β = 1 * 10^(-11)
    coriolis = FPlane(FT, f=f)
    coriolis = BetaPlane(FT, f₀ = f, β = β)

    α  = 2e-4     # [K⁻¹] Thermal expansion coefficient
    g  = 9.8061   # [m/s²] gravitational constant
    cᵖ = 3994.0   # [J/K]  heat capacity
    ρ  = 999.8    # [kg/m³] density
    h = 1000.0     # [m] e-folding length scale for northern sponge
    ΔB = 8 * α * g # [m/s²] total change in buoyancy from surface to bottom
    eos = LinearEquationOfState(FT, α=α, β=0)
    buoyancy = BuoyancyTracer()

    κh = 0.5e-5 # [m²/s] horizontal diffusivity
    νh = 12.0   # [m²/s] horizontal viscocity
    κv = 0.5e-5 # [m²/s] vertical diffusivity
    νv = 3e-4   # [m²/s] vertical viscocity

    closure = AnisotropicDiffusivity(νx = νh, νy = νh, νz =νv,
                                    κx = κh, κy = κh, κz=κv)

    parameters = (
        Ly = Ly,                   # y-domain length
        τ = τ,                   # [N m⁻²] Zonal stress
        ρ = ρ,                     # [kg / m³]
        μ = 1.1e-3,                # [m/s]  linear drag
        H = Lz,                    # [m]
        h = h,                     # [m]    relexaction profile scale
        ΔB = ΔB,                   # [m/s²] buoyancy jump
        Lz = Lz,                   # [m]
        Lsponge = 1900kilometer,   # [m]
        λᵗ = 7.0day,               # [s]
        Qᵇ = Q/(ρ * cᵖ) * α * g,  # [m² / s³]
        Qᵇ_cutoff = Ly * 5/6.      # [m]
    )

    # Momentum Boundary Conditions
    @inline windshape(y, p) = sin( π * y / p.Ly)
    @inline windstress(x, y, t, p) = - p.τ / p.ρ * windshape(y, p)
    @inline ulineardrag(i, j, grid, clock, state, p) = @inbounds - p.μ * state.u[i, j, 1]
    @inline vlineardrag(i, j, grid, clock, state, p) = @inbounds - p.μ * state.v[i, j, 1]
    # Zonal Velocity
    top_u_bc = BoundaryCondition(Flux, windstress, parameters = parameters)
    bottom_u_bc =  BoundaryCondition(Flux, ulineardrag, discrete_form = true, parameters = parameters)
    u_bcs = UVelocityBoundaryConditions(grid, top = top_u_bc, bottom = bottom_u_bc)
    # Meridional Velocity
    bottom_v_bc =  BoundaryCondition(Flux, vlineardrag, discrete_form = true, parameters = parameters)
    v_bcs = VVelocityBoundaryConditions(grid, bottom = bottom_v_bc)

    # Buoyancy Boundary Conditions Forcing. Note: Flux convention opposite of Abernathy
    @inline cutoff(j, grid, p ) = grid.yC[j] > p.Qᵇ_cutoff ? -0.0 : 1.0
    @inline surface_flux(j, grid, p) = p.Qᵇ * cos(3π * grid.yC[j] / p.Ly) * cutoff(j, grid, p)
    @inline relaxation(i, j, grid, clock, state, p) = @inbounds surface_flux(j, grid, p)
    top_b_bc = BoundaryCondition(Flux, relaxation, discrete_form = true, parameters = parameters)
    b_bcs = TracerBoundaryConditions(grid, top = top_b_bc)

    # Save boundary conditions as named tuple
    bcs = (b = b_bcs,  u = u_bcs, v = v_bcs,)

    # Forcing Functions
    # Sponge layers
    relu(y) = (abs(y) + y) * 0.5
    # Northern Wall Relaxation
    @inline relaxation_profile_north(k, grid, p) = p.ΔB * ( exp(grid.zC[k]/p.h) - exp(-p.Lz/p.h) ) / (1 - exp(-p.Lz/p.h))
    function Fb_function(i, j, k, grid, clock, state, p)
        return @inbounds - (1/p.λᵗ)  * (state.b[i,j,k]
            - relaxation_profile_north(k, grid, p)) * relu( (grid.yC[j]-p.Lsponge) / (p.Ly - p.Lsponge))
    end
    Fb = Forcing(Fb_function, parameters = parameters, discrete_form = true)

    # Record forcings
    forcings = (b = Fb, )

    # Model Setup
    model = IncompressibleModel(
            architecture = arch,
                float_type = FT,
                    grid = grid,
                coriolis = coriolis,
                buoyancy = buoyancy,
                    closure = closure,
                    tracers = (:b,),
        boundary_conditions = bcs,
                    forcing = forcings,
                advection = advection,
                timestepper = timestepper,
    )

    # Initial Conditions
    if ic_load
        println("loading from " * filepath)
        ic!(model, filepath, ArrayType = archarray(arch))
    else
        ε(σ) = σ * randn()
        B₀(x, y, z) = ΔB * ( exp(z/h) - exp(-Lz/h) ) / (1 - exp(-Lz/h)) + ε(1e-8)
        set!(model, b=B₀)
    end

    # Input Output Details
    ## Progress printer
Ni = 1000

function print_progress(simulation)
    model = simulation.model
    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (model.clock.time / end_time)

    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)
    meanb = mean(view(model.tracers.b.data, 2:Nx-1, 2:Ny-1, 2:Nz-1))

    @printf("[%05.2f%%] i: %d, t: %.2e days, umax: (%6.3e, %6.3e, %6.3e) m/s, CFL: %6.4e, next Δt: %.1e s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), Δt_wizard.Δt)
    println(" ")
    @printf("The mean value of b is %.2e", meanb)
    println(" ")
end

## Diagnostics
if write_slices
fields = Dict(
    "u" => model.velocities.u,
    "v" => model.velocities.v,
    "w" => model.velocities.w,
    "b" => model.tracers.b
)

surface_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_surface.nc",
			           schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(k = Nz))

middepth_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_middepth.nc",
                       schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(k = Int(floor(Nz/2))))

zonal_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_zonal.nc",
                       schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(j = Int(floor(Ny/2))))

meridional_output_writer =
    NetCDFOutputWriter(model, fields, filepath = filename * "_meridional.nc",
                       schedule=TimeInterval(slice_output_interval), field_slicer = FieldSlicer(i = Int(floor(Nx/2))))
function debug_f(debug)
    if debug
        close(surface_output_writer)
        close(middepth_output_writer)
        close(zonal_output_writer)
        close(meridional_output_writer)
    end
    return nothing
end
end

## Zonal and Time Averages
if write_zonal
u, v, w = model.velocities
b = model.tracers.b

u_scratch = XFaceField(model.architecture, model.grid)
v_scratch = YFaceField(model.architecture, model.grid)
w_scratch = ZFaceField(model.architecture, model.grid)
b_scratch =  CellField(model.architecture, model.grid)

zonal_fields = Dict(
    :u => AveragedField( u, dims= 1),
    :v => AveragedField( v, dims= 1),
    :w => AveragedField( w, dims= 1),
    :b => AveragedField( b, dims= 1),
    :uu => AveragedField(u * u, dims= 1),
    :vv => AveragedField(v * v, dims= 1),
    :ww => AveragedField(w * w, dims= 1),
    :uv => AveragedField(u * v, dims= 1),
    :vw => AveragedField(w * v, dims= 1),
    :uw => AveragedField(w * u, dims= 1),
    :ub => AveragedField(u * b, dims= 1),
    :vb => AveragedField(v * b, dims= 1),
    :wb => AveragedField(w * b, dims= 1),
)
#
# AveragedTimeInterval(time_avg_window, window=zonal_output_interval, stride=1)
zonalschedule = AveragedTimeInterval(zonal_output_interval, window=time_avg_window, stride=5) # TimeInterval(checkpoint_interval)
zonal_statistics = JLD2OutputWriter(model, zonal_fields,
                                    schedule = zonalschedule,
                                    prefix = filename * "_zonal_averages", force = true)
end
## Checkpointer
checkpointer = Checkpointer(model, prefix = filename * "_checkpoint",
                            schedule = TimeInterval(checkpoint_interval),
                            force = true)

    # Create Simulation
    Δt_wizard = TimeStepWizard(cfl = 1.0, Δt = Δt, max_change = 1.05, max_Δt = maxΔt)
    cfl = AdvectiveCFL(Δt_wizard)
    simulation = Simulation(model, Δt=Δt_wizard,
                            stop_time=end_time,
                            iteration_interval=Ni,
                            progress=print_progress)
    # IO
    if write_slices
        simulation.output_writers[:surface] = surface_output_writer
        simulation.output_writers[:middepth] = middepth_output_writer
        simulation.output_writers[:zonal] = zonal_output_writer
        simulation.output_writers[:meridional] = meridional_output_writer
    end
    if write_zonal
        simulation.output_writers[:statistics] = zonal_statistics
    end
    simulation.output_writers[:checkpoint] = checkpointer

    ## add output for Xiaozhou
    #=
    u, v, w = model.velocities
    b = model.tracers.b
    Xiaozhou_fields = Dict(
        :v => v,
        :b => b,
    )
    Xiaozhou_schedule = TimeInterval(5day)
    location = "/storage6/MesoscaleRuns/"
    Xiaozhou_output = JLD2OutputWriter(model, Xiaozhou_fields,
                                        schedule = Xiaozhou_schedule,
                                        prefix = location * filename * "_moc_data", force = true)
    simulation.output_writers[:Xiaozhou] = Xiaozhou_output
    =#
    ## Run
    run!(simulation)
    return nothing
end

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
using IntUtils

function main()
    if length(ARGS) > 0
        n = parse(Int64, ARGS[1])
        primes = sieve_of_eratosthenes(1_000_000)
        num_primes = 0
        nth = -1
        for i in 1:length(primes)
            if primes[i]
                num_primes += 1
                if num_primes == n
                    nth = i
                    break
                end
            end
        end

        @printf("The %dth prime is %d\n", n, nth)
    else
        error("this program needs a single numeric argument")
    end
end

main()


text:
export Transform, IdentityTransform, ScaleTransform, LowRankTransform, FunctionTransform, ChainTransform
export transform

"""
```julia
    transform(t::Transform, X::AbstractMatrix)
    transform(k::Kernel, X::AbstractMatrix)
```
Apply the transfomration `t` or `k.transform` on the input `X`
"""
transform

include("scaletransform.jl")
include("lowranktransform.jl")
include("functiontransform.jl")
include("selecttransform.jl")

"""
Chain a series of transform, here `t1` will be called first
```
    t1 = ScaleTransform()
    t2 = LowRankTransform(rand(3,4))
    ct = ChainTransform([t1,t2]) #t1 will be called first
    ct == t2∘t1
```
"""
struct ChainTransform <: Transform
    transforms::Vector{Transform}
end

Base.length(t::ChainTransform) = length(t.transforms) #TODO Add test

function ChainTransform(v::AbstractVector{<:Transform})
    ChainTransform(v)
end

function ChainTransform(v::AbstractVector{<:Type{<:Transform}},θ::AbstractVector)
    ChainTransform(v.(θ))
end

function transform(t::ChainTransform,X::T,obsdim::Int=defaultobs) where {T}
    Xtr = copy(X)
    for tr in t.transforms
        Xtr = transform(tr,Xtr,obsdim)
    end
    return Xtr
end

set_params!(t::ChainTransform,θ) = set_params!(t.transforms,θ)
params(t::ChainTransform) = (base_transform(t.transforms),params.(t.transforms))

Base.:∘(t₁::Transform,t₂::Transform) = ChainTransform([t₂,t₁])
Base.:∘(t::Transform,tc::ChainTransform) = ChainTransform(vcat(tc.transforms,t)) #TODO add test
Base.:∘(tc::ChainTransform,t::Transform) = ChainTransform(vcat(t,tc.transforms))

"""
IdentityTransform
Return exactly the input
"""
struct IdentityTransform <: Transform end

params(t::IdentityTransform) = nothing

transform(t::IdentityTransform,x::AbstractArray,obsdim::Int=defaultobs) = x #TODO add test

### TODO Maybe defining adjoints could help but so far it's not working

# @adjoint function ScaleTransform(s::T) where {T<:Real}
#     @check_args(ScaleTransform, s, s > zero(T), "s > 0")
#     ScaleTransform{T}(s),Δ->ScaleTransform{T}(Δ)
# end
#
# @adjoint function ScaleTransform(s::A) where {A<:AbstractVector{<:Real}}
#     @check_args(ScaleTransform, s, all(s.>zero(eltype(A))), "s > 0")
#     ScaleTransform{A}(s),Δ->begin; @show Δ,size(Δ); ScaleTransform{A}(Δ); end
# end

# @adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},x::AbstractVector{<:Real}) = transform(t,x),Δ->(ScaleTransform(nothing),t.s.*Δ)
#
#     @adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int) = transform(t,X,obsdim),Δ->begin
#     @show Δ,size(Δ);
#     return (obsdim == 1 ? ScaleTransform()Δ'.*X : ScaleTransform()Δ.*X,transform(t,Δ,obsdim),nothing)
#     end
#
# @adjoint transform(t::ScaleTransform{T},x::AbstractVecOrMat,obsdim::Int) where {T<:Real} = transform(t,x), Δ->(ScaleTransform(one(T)),t.s.*Δ,nothing)


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using RELOG, JSON, GZip

@testset "Reports" begin
    @testset "from solve" begin
        solution = RELOG.solve("$(pwd())/../instances/s1.json")
        tmp_filename = tempname()
        # The following should not crash
        RELOG.write_plant_emissions_report(solution, tmp_filename)
        RELOG.write_plant_outputs_report(solution, tmp_filename)
        RELOG.write_plants_report(solution, tmp_filename)
        RELOG.write_products_report(solution, tmp_filename)
        RELOG.write_transportation_emissions_report(solution, tmp_filename)
        RELOG.write_transportation_report(solution, tmp_filename)
    end
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
struct Polytope{T<:Integer}
    points::Matrix{Rational{T}}
    rays::Matrix{T}
    bidrays::Matrix{T}
end

function find_generators(A::Matrix{T}, b::Vector{T}) where T<:Rational
    mat = [A b]
    imat = intmat(mat)
    Aint, bint = imat[:,1:end-1], imat[:,end]
    find_generators(Aint, bint)
end

"""
# Description
Find generators for the polytope
  P = { x : Ax ≥ b }
that is, calculate the following representation
  P = { ∑ᵢ λᵢUᵢ + ∑ⱼ μⱼVⱼ + ∑ₖ νₖWₖ : ∑ᵢλᵢ = 1, μ ≥ 0 }

This is done through the use of the generalised Chernikova
algorithm on the cone homogenization of the above system.

# Returns
(points, rays, bidrays):
  `points::Matrix{Rational{T}}`: matrix whose columns correspond extremal points of polytope
  `rays::Matrix{T}`: matrix whose columns correspond to unidirectional rays of polytope
  `bidrays::Matrix{T}`: matrix whose columns correspond to bidirectional rays of polytopes
"""
function find_generators(A::Matrix{T}, b::Vector{T}) where T<:Integer
    # Create homogenized system
    m, n = size(A)
    Ab = [[A -b]; zeros(T, 1, n+1)]
    Ab[m+1,n+1] = one(T)
    bid, uni = chernikova_general(Ab)

    n_bid = size(bid, 2)
    n_uni = size(uni, 2)

    # Count number of trailing zeros (rays of polytope)
    # Note that all bidirectional rays in homogenized system
    # must have a trailing zero
    # (due to positivity constraint on this coordinate)
    n_zeros_uni = 0
    for j in 1:n_uni
        if uni[n+1,j] == 0 n_zeros_uni += 1 end
    end

    # If no unidirectional rays intersect with {x : x[n+1] == 0}
    # then polytope is empty
    if n_zeros_uni == n_uni
        return Array{Rational{T}}(undef, n, 0), Array{T}(undef, n, 0), Array{T}(undef, n, 0)
    end

    # Allocate matrices and get points.
    # The vertices of the polytope correspond to the unidirectional rays with non-zero trailing coordinate
    points = Array{Rational{T}}(undef, n, (n_uni - n_zeros_uni))
    rays = Array{T}(undef, n, n_zeros_uni)
    p_count = 1
    r_count = 1
    for j in 1:n_uni
        if uni[n+1, j] != 0
            points[:, p_count] = uni[1:n,j]//uni[n+1,j]
            p_count+=1
        else
            rays[:,r_count] = uni[1:n,j]
            r_count+=1
        end
    end

    bidrays = bid[1:n,:]
    return points, rays, bidrays
end

"""
# Description
Calculates the minimum and maximum values of
the projection of the polytope onto
specified vectors:

  min_{x ∈ P} Tᵢx
  max_{x ∈ P} Tᵢx

where Tᵢ are the columns of a matrix T.
# Arguments
* `P::Polytope{I}`: Polytope to project
* `T::Matrix{I}`: Matrix onto whose columns we project polytope
# Returns
`(mins::Array{Rational{T}}, maxs::Array{Rational{T}})`
"""
function min_max_projections(P::Polytope{S}, T::Matrix{Rational{S}} = Array{Rational{S}}(I, length(P), length(P))) where S<:Integer
    n = size(T, 2)
    n_points, n_rays, n_bidrays = size(P.points, 2), size(P.rays, 2), size(P.bidrays, 2)
    mins, maxs = fill(typemax(Rational{S}), n), fill(typemin(Rational{S}), n)
    for i in 1:n
        for j in 1:n_points
            tp = dot(T[:,i], P.points[:,j])
            if tp < mins[i]
                mins[i] = tp
            end
            if tp > maxs[i]
                maxs[i] = tp
            end
        end
        for j in 1:n_rays
            tr = dot(T[:,i], P.rays[:,j])
            if tr < 0
                mins[i] = -Inf
                break
            elseif tr > 0
                maxs[i] = Inf
                break
            end
        end
        for j in 1:n_bidrays
            if dot(T[:,i], P.bidrays[:,j]) != 0
                maxs[i] = Inf
                mins[i] = -Inf
                break
            end
        end
    end
    return mins, maxs
end

length(P::Polytope{T} where T<:Integer) = size(P.points, 1)


text:
using Colors, Test, FixedPointNumbers
using Colors: _de2000_t, _de2000_rot

# Test the colordiff function against example input and output from:
#
#   Sharma, G., Wu, W., & Dalal, E. N. (2005). The CIEDE2000 color‐difference
#   formula: Implementation notes, supplementary test data, and mathematical
#   observations. Color Research & Application, 30(1), 21–30. doi:10.1002/col
#

@testset "ColorDiff" begin
    abds = [
        ((50.0000,   2.6772, -79.7751), (50.0000,   0.0000, -82.7485),  2.0425),
        ((50.0000,   3.1571, -77.2803), (50.0000,   0.0000, -82.7485),  2.8615),
        ((50.0000,   2.8361, -74.0200), (50.0000,   0.0000, -82.7485),  3.4412),
        ((50.0000,  -1.3802, -84.2814), (50.0000,   0.0000, -82.7485),  1.0000),
        ((50.0000,  -1.1848, -84.8006), (50.0000,   0.0000, -82.7485),  1.0000),
        ((50.0000,  -0.9009, -85.5211), (50.0000,   0.0000, -82.7485),  1.0000),
        ((50.0000,   0.0000,   0.0000), (50.0000,  -1.0000,   2.0000),  2.3669),
        ((50.0000,  -1.0000,   2.0000), (50.0000,   0.0000,   0.0000),  2.3669),
        ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0009),  7.1792),
        ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0010),  7.1792),
        ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0011),  7.2195),
        ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0012),  7.2195),
        ((50.0000,  -0.0010,   2.4900), (50.0000,   0.0009,  -2.4900),  4.8045),
        ((50.0000,  -0.0010,   2.4900), (50.0000,   0.0010,  -2.4900),  4.8045),
        ((50.0000,  -0.0010,   2.4900), (50.0000,   0.0011,  -2.4900),  4.7461),
        ((50.0000,   2.5000,   0.0000), (50.0000,   0.0000,  -2.5000),  4.3065),
        ((50.0000,   2.5000,   0.0000), (73.0000,  25.0000, -18.0000), 27.1492),
        ((50.0000,   2.5000,   0.0000), (61.0000,  -5.0000,  29.0000), 22.8977),
        ((50.0000,   2.5000,   0.0000), (56.0000, -27.0000,  -3.0000), 31.9030),
        ((50.0000,   2.5000,   0.0000), (58.0000,  24.0000,  15.0000), 19.4535),
        ((50.0000,   2.5000,   0.0000), (50.0000,   3.1736,   0.5854),  1.0000),
        ((50.0000,   2.5000,   0.0000), (50.0000,   3.2972,   0.0000),  1.0000),
        ((50.0000,   2.5000,   0.0000), (50.0000,   1.8634,   0.5757),  1.0000),
        ((50.0000,   2.5000,   0.0000), (50.0000,   3.2592,   0.3350),  1.0000),
        ((60.2574, -34.0099,  36.2677), (60.4626, -34.1751,  39.4387),  1.2644),
        ((63.0109, -31.0961,  -5.8663), (62.8187, -29.7946,  -4.0864),  1.2630),
        ((61.2901,   3.7196,  -5.3901), (61.4292,   2.2480,  -4.9620),  1.8731),
        ((35.0831, -44.1164,   3.7933), (35.0232, -40.0716,   1.5901),  1.8645),
        ((22.7233,  20.0904, -46.6940), (23.0331,  14.9730, -42.5619),  2.0373),
        ((36.4612,  47.8580,  18.3852), (36.2715,  50.5065,  21.2231),  1.4146),
        ((90.8027,  -2.0831,   1.4410), (91.1528,  -1.6435,   0.0447),  1.4441),
        ((90.9257,  -0.5406,  -0.9208), (88.6381,  -0.8985,  -0.7239),  1.5381),
        (( 6.7747,  -0.2908,  -2.4247), ( 5.8714,  -0.0985,  -2.2286),  0.6377),
        ((2.0776,    0.0795,  -1.1350), ( 0.9033,  -0.0636,  -0.5514),  0.9082)]

    eps_cdiff = 0.00005

    @testset "CIEDE2000" begin
        metric = DE_2000()
        for (a, b, dexpect) in abds
            a64, b64 = Lab(a...), Lab(b...)
            @test colordiff(a64, b64; metric=metric) ≈ dexpect atol=eps_cdiff
            @test colordiff(b64, a64; metric=metric) ≈ dexpect atol=eps_cdiff
            a32, b32 = Lab{Float32}(a64), Lab{Float32}(b64)
            @test colordiff(a32, b32; metric=metric) ≈ dexpect atol=eps_cdiff
            @test colordiff(b32, a32; metric=metric) ≈ dexpect atol=eps_cdiff
        end

        h64 = 0.0:0.1:360.0
        h32 = 0.0f0:0.1f0:360.0f0
        @test all(h -> isapprox(_de2000_t(h), _de2000_t(big(h)), atol=3eps(Float64)), h64)
        @test all(h -> isapprox(_de2000_t(h), _de2000_t(big(h)), atol=2eps(Float32)), h32)
        @test all(h -> isapprox(_de2000_rot(h), _de2000_rot(big(h)), atol=eps(Float64)), h64)
        @test all(h -> isapprox(_de2000_rot(h), _de2000_rot(big(h)), atol=eps(Float32)), h32)
    end

    jl_red    =    RGB{N0f8}(Colors.JULIA_LOGO_COLORS.red)
    jl_green  = RGB{Float32}(Colors.JULIA_LOGO_COLORS.green)
    jl_blue   = HSV{Float64}(Colors.JULIA_LOGO_COLORS.blue)
    jl_purple = Lab{Float32}(Colors.JULIA_LOGO_COLORS.purple)
    colors = (jl_red, jl_green, jl_blue, jl_purple)
    pairs = ((a, b) for a in colors for b in Iterators.filter(c -> c != a,colors))
    tri = (colorant"navy", colorant"peru", colorant"gold")

    @testset "properties of metrics" begin
        metrics = (DE_2000(), DE_94(), DE_JPC79(), DE_CMC(), DE_BFD(),
                   DE_AB(), DE_DIN99(), DE_DIN99d(), DE_DIN99o())

        @testset "$metric" for metric in metrics
            # identity of indiscernibles
            @test all(c -> colordiff(c, c; metric=metric) == 0, colors)
            # positivity
            @test all(p -> colordiff(p[1], p[2]; metric=metric) > 0, pairs)
            # symmetry
            if metric isa DE_94 || metric isa DE_CMC # quasimetric
                @test all(p -> abs(colordiff(p[1], p[2]; metric=metric) -
                                   colordiff(p[2], p[1]; metric=metric)) > eps_cdiff, pairs)
            else
                @test all(p -> abs(colordiff(p[1], p[2]; metric=metric) -
                                   colordiff(p[2], p[1]; metric=metric)) < eps_cdiff, pairs)
            end
            # triangle inequality
            @test all(p -> colordiff(p[1],   p[2]; metric=metric) <=
                           colordiff(p[1], jl_red; metric=metric) +
                           colordiff(jl_red, p[2]; metric=metric), pairs)
            if any(m -> metric isa m, (DE_2000, DE_94, DE_CMC)) # semimetric
                @test colordiff(tri[1], tri[3]; metric=metric) >
                      colordiff(tri[1], tri[2]; metric=metric) +
                      colordiff(tri[2], tri[3]; metric=metric)
            end
        end
    end

    @testset "colordiff with grays" begin
        a, b = rand(100), rand(100)
        @test all(@. colordiff(a, b) == colordiff(Gray(a), b) == colordiff(a, Gray(b)) == colordiff(Gray(a), Gray(b)))
    end

    @testset "colordiff with transparent colors" begin
        @test colordiff(RGBA(1,0,0,1), RGB(1,0,0)) == 0
        @test colordiff(RGB(1,0,0), RGBA(1,0,0,1), metric=DE_AB()) == 0
        @test colordiff(RGBA(1,0,0,1), RGBA(1,0,0,1)) == 0
        @test_throws ArgumentError colordiff(RGBA(1,0,0,0.5), RGB(1,0,0))
        @test_throws ArgumentError colordiff(RGBA(1,0,0,0.5), RGBA(1,0,0,0.5))
    end
end


text:
#=
# Gradients with Zygote

It is possible to calculate gradients with Zygote.
Let's compare performance to explicit gradients.
=#

import MKL
import MLDatasets
import Makie
import CairoMakie
import RestrictedBoltzmannMachines as RBMs
using ValueHistories: MVHistory
using RestrictedBoltzmannMachines: BinaryRBM, initialize!, cd!, cdad!
nothing #hide

# Setup

Float = Float32
epochs = 100
batchsize = 128
nothing #hide

# Load MNIST

train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float}(train_x[:, :, train_y .== 0] .≥ 0.5)
nothing #hide

# Train using explicit gradients

rbm_∂s = BinaryRBM(Float, (28,28), 128)
@time cd!(rbm_∂s, train_x) # warm-up run so as to not consider pre-compilation times
initialize!(rbm_∂s, train_x)
history_∂s = MVHistory()
time_0 = time()
@time cd!(
    rbm_∂s, train_x; epochs, batchsize,
    callback = function(@nospecialize(args...); @nospecialize(kw...))
        push!(history_∂s, :t, time() - time_0)
    end
)
nothing #hide

# Train using Zygote gradients

rbm_ad = BinaryRBM(Float, (28,28), 128)
@time cdad!(rbm_ad, train_x) # warm-up run so as to not consider pre-compilation times
initialize!(rbm_ad, train_x)
history_ad = MVHistory()
time_0 = time()
@time cdad!(
    rbm_ad, train_x; epochs, batchsize,
    callback = function(@nospecialize(args...); @nospecialize(kw...))
        push!(history_ad, :t, time() - time_0)
    end
)
nothing #hide

# Compare timings

fig = Makie.Figure(resolution=(600, 400))
ax = Makie.Axis(fig[1,1], xlabel="batch", ylabel="seconds")
Makie.lines!(ax, get(history_∂s, :t)..., label="manual")
Makie.lines!(ax, get(history_ad, :t)..., label="zygote")
Makie.axislegend(ax, position=:rt)
fig


text:
module BasicTest

using Printf

include("../src/Dubins3D.jl")
using .Dubins3D

include("instances3D.jl")

function print_results(data, name)
  println("NAME ", name , " -------------  ")
  for i in 1:length(data)
    d = data[i]
    qi = d[1]
    qf = d[2]
    dubins = DubinsManeuver3D(qi, qf, rhomin, pitchmax)

    t = @timed DubinsManeuver3D(qi, qf, rhomin, pitchmax)

    lb = getLowerBound(qi, qf, rhomin, pitchmax)
    ub = getUpperBound(qi, qf, rhomin, pitchmax)
    best_lb = lb.length

    print(join([
      "$(@sprintf("%6s %d", name, i))",
      "$(@sprintf(" L[m]:%9.3f", dubins.length))",
      "$(@sprintf(" T[ms]: %5.3f", (t[2]) * 1000))",
      "$(@sprintf(" LB[m]: %8.3f", lb.length))",
      "$(@sprintf(" GAP[%%]: %6.3f", 100.0*(dubins.length - best_lb)/dubins.length))",
      "\n",
    ]," "))
  end
end

print_results(LONG, "Long")
print_results(SHORT, "Short")
print_results(BEARD, "Beard")
print_results(OTHERS, "Others")
print_results(HOTA, "Hota")

COUNT = 100

print("END -------------  \n\n")

function bold(best, text)
  if best
    return "\\textbf{" * text * "}"
  else
    return "        " * text * " "
  end
end

function text_value(value)
  if value > 100000000.0
    return "    --- "
  else
    return "$(@sprintf("%8.2f", value))"
  end
end

function text_value_gap(value)
  if isnan(value)
    return "    --- "
  else
    return "$(@sprintf("%8.3f", value))"
  end
end

end


text:
using Test
using OpenQuantumSystems
using SparseArrays, LinearAlgebra

@testset "metrics" begin

    b1 = SpinBasis(1 // 2)
    b2 = FockBasis(6)

    psi1 = spinup(b1) ⊗ coherentstate(b2, 0.1)
    psi2 = spindown(b1) ⊗ fockstate(b2, 2)

    rho = tensor(psi1, dagger(psi1))
    sigma = tensor(psi2, dagger(psi2))
    sigma_s = Commutator(sigma)
    rho_s = Commutator(rho)

    # tracenorm
    @test tracenorm(0 * rho_s) ≈ 0.0
    @test tracenorm_h(0 * rho_s) ≈ 0.0
    @test tracenorm_nh(0 * rho_s) ≈ 0.0

    @test tracenorm(rho_s) ≈ 26.0 # 1 * 6 * 4
    @test tracenorm_h(rho_s) ≈ 26.0
    @test tracenorm_nh(rho_s) ≈ 26.0

    # tracedistance
    @test tracedistance(rho_s, sigma_s) ≈ 26.0
    @test tracedistance_h(rho_s, sigma_s) ≈ 26.0
    @test tracedistance_nh(rho_s, sigma_s) ≈ 26.0

    @test tracedistance(rho_s, rho_s) ≈ 0.0
    @test tracedistance_h(rho_s, rho_s) ≈ 0.0
    @test tracedistance_nh(rho_s, rho_s) ≈ 0.0

    @test tracedistance(sigma_s, sigma_s) ≈ 0.0
    @test tracedistance_h(sigma_s, sigma_s) ≈ 0.0
    @test tracedistance_nh(sigma_s, sigma_s) ≈ 0.0

    # tracedistance
    @test tracedistance(rho_s, sigma_s) ≈ 26.0
    @test tracedistance(rho_s, rho_s) ≈ 0.0
    @test tracedistance(sigma_s, sigma_s) ≈ 0.0

    @test tracedistance_nh(rho_s, rho_s) ≈ 0.0

end # testset

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems


text:
# Managing save data, and creating a unified approach to the results.jld2 that is relatively extensible.

function save(x)
end


text:
using OperatorFlux
using Test
using Flux
using LinearAlgebra

@testset "SpectralConv - FourierTransform" begin
    # auxiliary data
    M = 32

    # test constructor
    trafo = FourierTransform(modes = (10, 23))
    ch = 3 => 13
    sc = SpectralConv(trafo, ch)
    @test ndims(sc.weights) == 4

    # test convolution
    a, b = -π, π
    x = (b - a) .* collect(0:(M - 1)) / M .+ a
    y = copy(x)
    z = sin.(x) * sin.(y)'
    z = z[:, :, :, :] # emulate channels and batches
    z = permutedims(z, (3, 1, 2, 4))
    modes = (12, 16)
    trafo = FourierTransform(modes = modes)
    ch = 1 => 1
    sc = SpectralConv(trafo, ch)
    sc.weights .= ones(eltype(z), 1, size(trafo)..., 1)
    @test norm(sc(z) - z) ≤ eps(norm(z))
    @test ndims(sc.weights) == 4

    for M in [31, 32]
        even = Bool(mod(M, 2)) ? false : true

        # test more than one out_channel
        ch = 1 => 2
        trafo = FourierTransform(modes = (16, 16), even = even)
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 1, M, M, 2)
        @test size(sc(x)) == (2, M, M, 2)

        # test different dimensionalities
        # 1d
        ch = 1 => 2
        trafo = FourierTransform(modes = (16,), even = even)
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 1, M, 2)
        @test size(sc(x)) == (2, M, 2)

        # 3d
        ch = 13 => 3
        trafo = FourierTransform(modes = (16, 16, 14), even = even)
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 13, M, M, M, 2)
        @test size(sc(x)) == (3, M, M, M, 2)

        # 4d
        ch = 13 => 7
        trafo = FourierTransform(modes = (16, 16, 14, 3), even = even)
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 13, M, M, M, 63, 2)
        @test size(sc(x)) == (7, M, M, M, 63, 2)
    end

    # test Base extensions
    trafo = FourierTransform(modes = (10,))
    ch = 3 => 13
    sc = SpectralConv(trafo, ch)
    @test ndims(sc) == 1
end

@testset "SpectralConv - ChebyshevTransformTransform" begin
    # auxiliary data
    M = 32

    # test constructor
    trafo = ChebyshevTransform(modes = (10, 23))
    ch = 3 => 13
    sc = SpectralConv(trafo, ch)
    @test ndims(sc.weights) == 4

    # test convolution
    a, b = -π, π
    x = @. (cos(pi * (0:M) / M) + 1) / 2
    x = (b - a) .* x .+ a
    y = copy(x)
    z = x * y'
    z = z[:, :, :, :] # emulate channels and batches
    z = permutedims(z, (3, 1, 2, 4))
    modes = (12, 17)
    trafo = ChebyshevTransform(modes = modes)
    ch = 1 => 1
    sc = SpectralConv(trafo, ch)
    sc.weights .= ones(eltype(z), 1, collect(modes)..., 1)
    @test norm(sc(z) - z) ≤ 2 * eps(norm(z))
    @test ndims(sc.weights) == 4

    for M in [31, 32]
        # test more than one out_channel
        ch = 1 => 2
        trafo = ChebyshevTransform(modes = (16, 16))
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 1, M, M, 2)
        @test size(sc(x)) == (2, M, M, 2)

        # test different dimensionalities
        # 1d
        ch = 1 => 2
        trafo = ChebyshevTransform(modes = (16,))
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 1, M, 2)
        @test size(sc(x)) == (2, M, 2)

        # 3d
        ch = 13 => 3
        trafo = ChebyshevTransform(modes = (16, 16, 14))
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 13, M, M, M, 2)
        @test size(sc(x)) == (3, M, M, M, 2)

        # 4d
        ch = 13 => 7
        trafo = ChebyshevTransform(modes = (16, 16, 14, 3))
        sc = SpectralConv(trafo, ch)
        x = rand(Float32, 13, M, M, M, 63, 2)
        @test size(sc(x)) == (7, M, M, M, 63, 2)
    end

    # test Base extensions
    trafo = ChebyshevTransform(modes = (10,))
    ch = 3 => 13
    sc = SpectralConv(trafo, ch)
    @test ndims(sc) == 1
end

@testset "SpectralConv - LegendreTransform" begin
    # auxiliary data / legendre polynomials
    P₀(x) = 1
    P₁(x) = x
    P₂(x) = 0.5 * (3 * x^2 - 1)
    P₃(x) = 0.5 * (5 * x^3 - 3 * x)
    P₄(x) = 0.125 * (35 * x^4 - 30 * x^2 + 3)
    P₅(x) = 0.125 * (63 * x^5 - 70 * x^3 + 15 * x)
    polynomials = [P₀, P₁, P₂, P₃, P₄, P₅]

    # test constructors
    for mode in [(4,), (4, 5), (4, 5, 6)]
        trafo = LegendreTransform(modes = mode)
        ch = 3 => 13
        sc = SpectralConv(trafo, ch)
        @test ndims(sc.weights) == length(mode) + 2
    end

    # rountrip integration test
    for mode in [(5,), (5, 6), (5, 7, 9)]
        z = rand(1, mode..., 5)
        trafo = LegendreTransform(modes = mode)
        ch = 1 => 1
        sc = SpectralConv(trafo, ch)
        sc.weights .= ones(eltype(z), 1, collect(mode)..., 1)
        @test norm(sc(z) - z) ≤ 10 * eps(norm(z))
    end

    # test more than one out_channel
    for mode in [(5,), (5, 6), (5, 7, 9)]
        for c in [2, 3]
            z = rand(1, mode..., 5)
            trafo = LegendreTransform(modes = mode)
            ch = 1 => c
            sc = SpectralConv(trafo, ch)
            sc.weights .= ones(eltype(z), 1, collect(mode)..., 1)
            @test size(sc(z)) == (c, mode..., 5)
        end
    end

    # test Base extensions
    trafo = LegendreTransform(modes = (10,))
    ch = 3 => 13
    sc = SpectralConv(trafo, ch)
    @test ndims(sc) == 1
end


text:
using Colors, DifferentialEquations, LaTeXStrings, PGFPlotsX

function sir!(du, u, p, t)
    du[1] = -(p[1]/p[2] * u[2] * u[1])
    du[2] = (p[1]/p[2] * u[2] * u[1]) - u[2]/p[2]
    du[3] = u[2] / p[2]
end

function sirplot(r0, d, t, i0 = 1e-6)
    u0 = [1-i0, i0, 0.0]
    solv = solve(ODEProblem(sir!, u0, t, [r0, d]), saveat=1)
    r0lab = latexstring(replace("\\mathcal{R}_0=$(r0)", "."=>"{,}"))
    i0lab = latexstring(replace("I_0=$i0", "."=>"{,}"))
    dlab = latexstring("D=$d")
    p = @pgf Axis({"yticklabel style"={"/pgf/number format/use comma"},
       legend_pos="outer north east",
       title="SIR för $(r0lab), $(dlab), $(i0lab)", xlabel=L"t",
       ylabel=L"p",
       xmajorgrids, ymajorgrids})
    cs = ["S","I","R"]
    plotcolors = distinguishable_colors(length(cs)+1, [RGB(1,1,1)])[2:end]
    for (i, c) in enumerate(cs)
       @pgf push!(p, PlotInc({no_markers, color = plotcolors[i]},
           Table([solv.t, map(r -> r[i], solv.u)])),
           LegendEntry(latexstring(c)))
    end
    @pgf push!(p, PlotInc({no_markers, style = {"dashed"}, color = plotcolors[1]},
       Coordinates([(t[1], 1/r0), (t[2], 1/r0)])),
       LegendEntry(L"1/\mathcal{R}_0"))
    p
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust