text:
export
    AbstractCoordinate,
    LatLonAlt,
    ECEF,
    UTM,
    ENU,
    NED,

    WGS_84,
    INTERNATIONAL,

    check_is_in_radians,
    ensure_lon_between_pies,
    deg_min_sec_to_degrees,
    degrees_to_deg_min_sec,
    get_earth_radius,

    GeodeticDatum,
    eccentricity

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

"""
Defines a reference ellipsoid used to map a position on the Earth
relative to said ellipsoid
"""
struct GeodeticDatum
    a::Float64 # equitorial radius
    b::Float64 # polar radius
end
const WGS_84        = GeodeticDatum(6378137.0, 6356752.314245)
const INTERNATIONAL = GeodeticDatum(6378388.0, 6356911.94613)

eccentricity_sq(datum::GeodeticDatum) = 1.0 - datum.b^2 / datum.a^2
eccentricity(datum::GeodeticDatum) = sqrt(eccentricity_sq(datum))
function radius_of_curvature_in_the_meridian(ϕ::Float64, datum::GeodeticDatum)

    # ϕ = meridian angle (latitude)

    a = datum.a
    e = eccentricity(datum)
    e² = e*e
    sin²ϕ = sin(ϕ)^2
    ρ = a*(1-e²) / ((1-e²*sin²ϕ)^1.5)

    ρ
end
function radius_of_curvature_in_the_prime_vertical(ϕ::Float64, datum::GeodeticDatum=WGS_84)

    # ϕ = meridian angle (latitude)

    a = datum.a
    e = eccentricity(datum)
    e² = e*e
    sin²ϕ = sin(ϕ)^2
    ν = a / sqrt(1 - e²*sin²ϕ)

    ν
end

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

abstract type AbstractCoordinate end

"""
The latitude, longitude, altitude coordinate system
"""
struct LatLonAlt <: AbstractCoordinate
    lat::Float64 # [rad] ϕ, angle between the equatorial plane and the straight line that passes through the point and through the center of the Earth
    lon::Float64 # [rad] λ, angle (typically east) from a reference meridian to another meridian that passes through that point
    alt::Float64 # [m] above reference ellipsoid
end

Base.print(io::IO, lla::LatLonAlt) = @printf(io, "%6.2f %7.2f %10.2f", rad2deg(lla.lat), rad2deg(lla.lon), lla.alt)

function check_is_in_radians(lla::LatLonAlt)
    if abs(lla.lat) > π/2 || abs(lla.lon) > π
        warn("lat and lon may not be in radians!")
    end
    nothing
end
function ensure_lon_between_pies(lla::LatLonAlt)

    # lan ∈ [-π/2, π/2]
    # lon ∈ [-π, π]

    lon = lla.lon

    while lon ≤ -π
        lon += 2π
    end
    while lon ≥ π
        lon -= 2π
    end

    LatLonAlt(lla.lat, lon, lla.alt)
end
deg_min_sec_to_degrees(deg::Real, min::Real, sec::Real) = convert(Float64, deg + min/60.0 + sec/3600.0)
function degrees_to_deg_min_sec(degrees::Real)
    degrees = convert(Float64, degrees)
    deg, remainder = divrem(degrees, 1.0)
    min, remainder = divrem(remainder*60.0, 1.0)
    sec = remainder * 60.0
    (convert(Int, deg), convert(Int, min), sec)
end
function get_earth_radius(lat::Float64, datum::GeodeticDatum=WGS_84)

    a, b = datum.a, datum.b
    c, s = cos(lat), sin(lat)
    return sqrt(((a^2*c)^2 + (b^2*s)^2)/((a*c)^2 + (b*s)^2))
end

Base.convert(::Type{VecE3}, lla::LatLonAlt) = VecE3(lla.lat, lla.lon, lla.alt)
Base.convert(::Type{LatLonAlt}, p::VecE3) = LatLonAlt(p.x, p.y, p.z)

"""
The Earth-centered, Earth-fixed coordinate system
"""
struct ECEF <: AbstractCoordinate
    x::Float64 # [m]
    y::Float64 # [m]
    z::Float64 # [m]
end
Base.convert(::Type{VecE3}, p::ECEF) = VecE3(p.x, p.y, p.z)
Base.convert(::Type{ECEF}, p::VecE3) = ECEF(p.x, p.y, p.z)

"""
Universal Transverse Mercator coordinate system
"""

const UTM_LATITUDE_LIMIT_NORTH = deg2rad(84)
const UTM_LATITUDE_LIMIT_SOUTH = deg2rad(-80)

struct UTM <: AbstractCoordinate
    e::Float64 # [m]
    n::Float64 # [m]
    u::Float64 # [m] above reference ellipsoid
    zone::Int # UTM zone (∈ 1:60)
end
Base.convert(::Type{VecE3}, p::UTM) = VecE3(p.e, p.n, p.u)
function Base.:+(a::UTM, b::UTM)
    @assert(a.zone == b.zone)
    UTM(a.e+b.e, a.n+b.n, a.u+b.u, a.zone)
end
function Base.:-(a::UTM, b::UTM)
    @assert(a.zone == b.zone)
    UTM(a.e-b.e, a.n-b.n, a.u-b.u, a.zone)
end

"""
East North Up coordinate system
relative to some point on the surface of the earth
"""
struct ENU <: AbstractCoordinate
    e::Float64 # [m]
    n::Float64 # [m]
    u::Float64 # [m] above reference location
end
Base.convert(::Type{VecE3}, p::ENU) = VecE3(p.e, p.n, p.u)
Base.convert(::Type{ENU}, p::VecE3) = ENU(p.x, p.y, p.z)
Base.:+(a::ENU, b::ENU) = ENU(a.e+b.e, a.n+b.n, a.u+b.u)
Base.:-(a::ENU, b::ENU) = ENU(a.e-b.e, a.n-b.n, a.u-b.u)

"""
North East Down coordinate system
"""
struct NED <: AbstractCoordinate
    n::Float64 # [m]
    e::Float64 # [m]
    d::Float64 # [m] below reference location
end
Base.convert(::Type{VecE3}, p::NED) = VecE3(p.n, p.e, p.d)
Base.:+(a::NED, b::NED) = NED(a.n+b.n, a.e+b.e, a.d+b.d)
Base.:-(a::NED, b::NED) = NED(a.n-b.n, a.e-b.e, a.d-b.d)

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

function Base.convert(::Type{ECEF}, lla::LatLonAlt, datum::GeodeticDatum=WGS_84)

    ϕ = lla.lat
    λ = lla.lon
    h = lla.alt

    e² = eccentricity_sq(datum)

    sϕ = sin(ϕ)
    cϕ = cos(ϕ)
    sλ = sin(λ)
    cλ = cos(λ)

    Rₙ = radius_of_curvature_in_the_prime_vertical(ϕ, datum)

    x = (Rₙ         + h)*cϕ*cλ
    y = (Rₙ         + h)*cϕ*sλ
    z = (Rₙ*(1.0-e²)+ h)*sϕ

    ECEF(x,y,z)
end
function Base.convert(::Type{LatLonAlt}, ecef::ECEF, datum::GeodeticDatum=WGS_84)

    # https://microem.ru/files/2012/08/GPS.G1-X-00006.pdf

    #=
    Convert Earth-Centered-Earth-Fixed (ECEF) to lat, Lon, Altitude
    Input is x, y, z in meters
    Returns lat and lon in radians, and altitude in meters

    correction for altitude near poles left out.
    =#

    a = datum.a
    b = datum.b
    e² = eccentricity_sq(datum)

    x = ecef.x
    y = ecef.y
    z = ecef.z

    λ = atan2(y, x)
    p = sqrt(x*x + y*y)
    θ = atan2(z*a, (p*b))

    h = 0.0
    ϕ = atan2(z,p*(1.0-e²))

    # the equation diverges at the poles
    if isapprox(abs(θ), π/2, atol=0.1)
        Rₙ = radius_of_curvature_in_the_prime_vertical(ϕ, datum)
        ϕₛ = sin(ϕ)
        L = z + e²*Rₙ*ϕₛ
        h = L / ϕₛ - Rₙ
    else
        Δh = Inf
        while Δh > 1e-4
            N = a / sqrt(1-e²*sin(ϕ)^2)
            h₂ = p/cos(ϕ) - N
            ϕ = atan2(z, p*(1-e²*N/(N + h₂)))
            Δh = abs(h - h₂)
            h = h₂
        end
    end

    ensure_lon_between_pies(LatLonAlt(ϕ, λ, h))
end

function Base.convert(::Type{UTM}, lla::LatLonAlt, datum::GeodeticDatum=WGS_84, zone::Integer=-1)
    # see DMATM 8358.2 by the Army
    # code verified using their examples

    check_is_in_radians(lla)
    lla = ensure_lon_between_pies(lla)

    lat = lla.lat
    lon = lla.lon

    if lat > UTM_LATITUDE_LIMIT_NORTH || lat < UTM_LATITUDE_LIMIT_SOUTH
        error("latitude $(rad2deg(lat)) is out of limits!")
    end

    const FN = 0.0      # false northing, zero in the northern hemisphere
    const FE = 500000.0 # false easting
    const ko = 0.9996   # central scale factor

    zone_centers = -177.0*pi/180 + 6.0*pi/180*collect(0:59) # longitudes of the zone centers
    if zone == -1
        zone = indmin(map(x->abs(lon - x), zone_centers)) # index of min zone center
    end
    long0 = zone_centers[zone]

    s  = sin(lat)
    c  = cos(lat)
    t  = tan(lat)

    a = datum.a
    b = datum.b

    n  = (a-b)/(a+b)
    e₁  = sqrt((a^2-b^2)/a^2) # first eccentricity
    ep = sqrt((a^2-b^2)/b^2) # second eccentricity
    nu = a/(1-e₁^2*s^2)^0.5   # radius of curvature in the prime vertical
    dh = lon - long0         # longitudinal distance from central meridian

    A  = a*(1 - n + 5.0/4*(n^2-n^3) + 81.0/64*(n^4-n^5))
    B  = 3.0/2*a*(n-n^2 + 7.0/8*(n^3-n^4) + 55.0/64*n^5)
    C  = 15.0/16*a*(n^2-n^3 + 3.0/4*(n^4-n^5))
    D  = 35.0/48*a*(n^3-n^4 + 11.0/16.0*n^5)
    E  = 315.0/512*a*(n^4-n^5)
    S  = A*lat - B*sin(2*lat) + C*sin(4*lat) - D*sin(6*lat) + E*sin(8*lat)

    T1 = S*ko
    T2 = nu*s*c*ko/2
    T3 = nu*s*c^3*ko/24  * (5  -    t^2 + 9*ep^2*c^2 + 4*ep^4*c^4)
    T4 = nu*s*c^5*ko/720 * (61 - 58*t^2 + t^4 + 270*ep^2*c^2 - 330*(t*ep*c)^2
                               + 445*ep^4*c^4 + 324*ep^6*c^6 - 680*t^2*ep^4*c^4
                               + 88*ep^8*c^8  - 600*t^2*ep^6*c^6 - 192*t^2*ep^8*c^8)
    T5 = nu*s*c^7*ko/40320 * (1385 - 3111*t^2 + 543*t^4 - t^6)
    T6 = nu*c*ko
    T7 = nu*c^3*ko/6   * (1 - t^2 + ep^2*c^2)
    T8 = nu*c^5*ko/120 * (5 - 18*t^2 + t^4 + 14*ep^2*c^2 - 58*t^2*ep^2*c^2 + 13*ep^4*c^4
                            + 4*ep^6*c^6 - 64*t^2*ep^4*c^4 - 24*t^2*ep^6*c^6)
    T9 = nu*c^7*ko/5040 * (61 - 479*t^2 + 179*t^4 - t^6)

    E = FE + (dh*T6 + dh^3*T7 + dh^5*T8 + dh^7*T9)        # easting  [m]
    N = FN + (T1 + dh^2*T2 + dh^4*T3 + dh^6*T4 + dh^8*T5) # northing [m]
    U = lla.alt

    UTM(E, N, U, zone)
end
function Base.convert(::Type{LatLonAlt}, utm::UTM, datum::GeodeticDatum=WGS_84)
    # see DMATM 8358.2 by the Army
    # code verified using their examples

    zone = utm.zone
    @assert(1 ≤ zone ≤ 60)
    λₒ = deg2rad(6.0*zone - 183) # longitude of the origin of the projection

    kₒ = 0.9996 # central scale factor
    a = datum.a
    b = datum.b

    e = eccentricity(datum)
    e² = e*e
    e⁴ = e²*e²
    e⁶ = e⁴*e²

    eₚ = sqrt((a^2-b^2)/b^2) # second eccentricity
    eₚ² = eₚ*eₚ
    eₚ⁴ = eₚ²*eₚ²
    eₚ⁶ = eₚ⁴*eₚ²
    eₚ⁸ = eₚ⁴*eₚ⁴

    N = utm.n
    arc = N / kₒ
    μ = arc / (a * (1 - e²/4 - 3*e⁴/64 - 5*e⁶/256))
    ei = (1 - sqrt(1 - e²)) / (1 + sqrt(1 - e²))
    ca = 3*ei/2 - 27*ei^3 / 32.0
    cb = 21 * ei^2 / 16 - 55*ei^4 / 32
    cc = 151 * ei^3 / 96
    cd = 1097 * ei^4 / 512
    ϕₚ = μ + ca * sin(2μ) + cb*sin(4μ) + cc*sin(6μ) + cd*sin(8μ)

    # const FN = 0.0      # false northing, zero in the northern hemisphere
    E = utm.e
    FE = 500000.0 # false easting
    ΔE = E - FE
    ΔE² = ΔE*ΔE
    ΔE³ = ΔE²*ΔE
    ΔE⁴ = ΔE²*ΔE²
    ΔE⁵ = ΔE³*ΔE²
    ΔE⁶ = ΔE³*ΔE³
    ΔE⁷ = ΔE⁴*ΔE³
    ΔE⁸ = ΔE⁴*ΔE⁴

    tanϕₚ = tan(ϕₚ)
    tan²ϕₚ = tanϕₚ*tanϕₚ
    tan⁴ϕₚ = tan²ϕₚ*tan²ϕₚ
    tan⁶ϕₚ = tan⁴ϕₚ*tan²ϕₚ

    cosϕₚ = cos(ϕₚ)
    cos²ϕₚ = cosϕₚ*cosϕₚ
    cos⁴ϕₚ = cos²ϕₚ*cos²ϕₚ
    cos⁶ϕₚ = cos⁴ϕₚ*cos²ϕₚ
    cos⁸ϕₚ = cos⁴ϕₚ*cos⁴ϕₚ

    kₒ² = kₒ*kₒ
    kₒ³ = kₒ²*kₒ
    kₒ⁴ = kₒ²*kₒ²
    kₒ⁵ = kₒ³*kₒ²
    kₒ⁶ = kₒ³*kₒ³
    kₒ⁷ = kₒ⁴*kₒ³
    kₒ⁸ = kₒ⁴*kₒ⁴

    ν = radius_of_curvature_in_the_prime_vertical(ϕₚ, datum)
    ν² = ν*ν
    ν³ = ν²*ν
    ν⁴ = ν²*ν²
    ν⁵ = ν³*ν²
    ν⁶ = ν³*ν³
    ν⁷ = ν⁴*ν³
    ν⁸ = ν⁴*ν⁴

    ρ = radius_of_curvature_in_the_meridian(ϕₚ, datum)

    T10 = tanϕₚ/(2ρ*ν*kₒ²)
    T11 = tanϕₚ/(24ρ*ν³*kₒ⁴) * (5.0 + 3*tan²ϕₚ + eₚ²*cos²ϕₚ - 4*eₚ⁴*cosϕₚ - 9*tan²ϕₚ*eₚ²*cos²ϕₚ)
    T12 = tanϕₚ/(720ρ*ν⁵*kₒ⁶) * (61.0 + 90*tan²ϕₚ + 46.0*eₚ²*cos²ϕₚ + 45*tan⁴ϕₚ - 252*tan²ϕₚ*eₚ²*cos²ϕₚ -
                                   3*eₚ⁴*cos⁴ϕₚ + 100*eₚ⁶*cos⁶ϕₚ - 66*tan²ϕₚ*eₚ⁴*cos⁴ϕₚ -
                                   90*tan⁴ϕₚ*eₚ²*cos²ϕₚ + 88*eₚ⁸*cos⁸ϕₚ + 225*tan⁴ϕₚ*eₚ⁴*cos⁴ϕₚ +
                                   84*tan²ϕₚ*eₚ⁶*cos⁶ϕₚ - 192*tan²ϕₚ*eₚ⁸*cos⁸ϕₚ)
    T13 = tanϕₚ/(40320ρ*ν⁷*kₒ⁸) * (1385.0 + 3633*tan²ϕₚ + 4095*tan⁴ϕₚ + 1575*tan⁶ϕₚ)

    T14 = 1.0/(ν*cosϕₚ*kₒ)
    T15 = (1.0 + 2*tan²ϕₚ + eₚ²*cos²ϕₚ) / (6*(ν³*cosϕₚ*kₒ³))
    T16 = (5.0 + 6*eₚ²*cos²ϕₚ + 28*tan²ϕₚ - 3*eₚ⁴*cos⁴ϕₚ + 8*tan²ϕₚ*eₚ²*cos²ϕₚ +
           24*tan⁴ϕₚ - 4*eₚ⁶*cos⁶ϕₚ + 4*tan²ϕₚ*eₚ⁴*cos⁴ϕₚ + 24*tan²ϕₚ*eₚ⁶*cos⁶ϕₚ) /
            (120*(ν⁵*cosϕₚ*kₒ⁵))
    T17 = (61.0 + 662*tan²ϕₚ + 1320*tan⁴ϕₚ + 720*tan⁶ϕₚ)/(5040*(ν⁷*cosϕₚ*kₒ⁷))

    ϕ = ϕₚ - ΔE²*T10 + ΔE⁴*T11 - ΔE⁶*T12 + ΔE⁸*T13
    λ = λₒ + ΔE*T14 - ΔE³*T15 + ΔE⁵*T16 - ΔE⁷*T17
    alt = utm.u

    LatLonAlt(ϕ, λ, alt)
end

function Base.convert(::Type{ECEF}, enu::ENU, reference::LatLonAlt)

    # http://www.navipedia.net/index.php/Transformations_between_ECEF_and_ENU_coordinates

    ϕ = reference.lat
    λ = reference.lon

    sϕ, cϕ = sin(ϕ), cos(ϕ)
    sλ, cλ = sin(λ), cos(λ)

    R = [-sλ  -cλ*sϕ  cλ*cϕ;
          cλ  -sλ*sϕ  sλ*cϕ;
         0.0      cϕ     sϕ]

    p = [enu.e, enu.n, enu.u]
    ecef = R*p

    ref_ecef = convert(ECEF, reference)

    ECEF(ecef[1] + ref_ecef.x, ecef[2] + ref_ecef.y, ecef[3] +  ref_ecef.z)
end
function Base.convert(::Type{ENU}, ecef::ECEF, refLLA::LatLonAlt, refECEF::ECEF=convert(ECEF, refLLA))

    # http://www.navipedia.net/index.php/Transformations_between_ECEF_and_ENU_coordinates

    ϕ = refLLA.lat
    λ = refLLA.lon
    a = refLLA.alt

    sϕ, cϕ = sin(ϕ), cos(ϕ)
    sλ, cλ = sin(λ), cos(λ)

    R = [-sλ      cλ     0.0;
         -cλ*sϕ  -sλ*sϕ  cϕ;
          cλ*cϕ   sλ*cϕ  sϕ]

    p = [ecef.x - refECEF.x, ecef.y - refECEF.y, ecef.z - refECEF.z]
    enu = R*p

    ENU(enu[1], enu[2], enu[3])
end

Base.convert(::Type{NED}, enu::ENU) = NED(enu.n, enu.e, -enu.u)
Base.convert(::Type{ENU}, ned::NED) = ENU(ned.e, ned.n, -ned.d)


text:
"""
	autoaccelerate(a::AbstractVector) -> AcceleratedVector

Generate an `AcceleratedVector` from any `AbstractVector` by choosing an appropriate index
based on the type of the vector.
This function is a no-op for `AcceleratedVector`s.

```jldoctest; setup = :(using IndexedDims: autoaccelerate)
julia> autoaccelerate(1:3)
3-element UnitRange{Int64} + UniqueSortIndex:
 1
 2
 3

julia> autoaccelerate(["foo", "bar", "baz", "bar"])
4-element Array{String,1} + HashIndex (3 unique elements):
 "foo"
 "bar"
 "baz"
 "bar"

julia> autoaccelerate(["foo", "bar", "baz"])
3-element Array{String,1} + UniqueHashIndex:
 "foo"
 "bar"
 "baz"
```
"""
function autoaccelerate end

autoaccelerate(v::AcceleratedVector) = v

autoaccelerate(r::OrdinalRange) = accelerate(r, UniqueSortIndex)
function autoaccelerate(v::AbstractVector)
    sorted = issorted(v)
    uniq = allunique(v)

    if sorted && uniq
        return accelerate(v, UniqueSortIndex)
    elseif sorted
        return accelerate(v, SortIndex)
    elseif uniq
        return accelerate(v, UniqueHashIndex)
    else
        return accelerate(v, HashIndex)
    end
end

for index_type in (UniqueSortIndex, SortIndex, UniqueHashIndex)
    @eval _reindex_type(::Type{<:$index_type}, index::AbstractUnitRange) = $index_type
end

_reindex_type(::Type{<:Union{HashIndex, UniqueHashIndex}}, index) = HashIndex

function reindex(v::AcceleratedVector{<:Any, <:Any, I}, index) where I
    return accelerate(parent(v)[index], _reindex_type(I, index))
end

reindex(v::AcceleratedVector, ::Colon) = v

# no way to know whether a given AbstractArray is unique or sorted
function reindex(v::AcceleratedVector{<:Any, <:Any, I}, index::AbstractArray) where I
    return accelerate(parent(v)[index], HashIndex)
end

ext: jl


text:
type LinearCombination{T}<:AbstractLinearMap{T}
    maps::Vector{AbstractLinearMap}
    coeffs::Vector{T}
    function LinearCombination(maps::Vector{AbstractLinearMap},coeffs::Vector{T})
        N=length(maps)
        N==length(coeffs) || error("Number of coefficients doesn't match number of terms")
        sz=size(maps[1])
        for n=1:N
            size(maps[n])==sz || throw(DimensionMismatch("LinearCombination"))
            promote_type(T,eltype(maps[n]))==T || throw(InexactError())
        end
        new(maps,coeffs)
    end
end

# basic methods
Base.size(A::LinearCombination,n)=size(A.maps[1],n)
Base.size(A::LinearCombination)=size(A.maps[1])
Base.isreal(A::LinearCombination)=all(isreal,A.maps) && all(isreal,A.coeffs) # sufficient but not necessary
Base.issym(A::LinearCombination)=all(issym,A.maps) # sufficient but not necessary
Base.ishermitian(A::LinearCombination)=all(ishermitian,A.maps) && all(isreal,A.coeffs) # sufficient but not necessary
Base.isposdef(A::LinearCombination)=all(isposdef,A.maps) && all(isposdef,A.coeffs) # sufficient but not necessary

# adding linear maps
function +(A1::LinearCombination,A2::LinearCombination)
    size(A1)==size(A2) || throw(DimensionMismatch("+"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1.maps...,A2.maps...],T[A1.coeffs...,A2.coeffs...])
end
function +(A1::AbstractLinearMap,A2::LinearCombination)
    size(A1)==size(A2) || throw(DimensionMismatch("+"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1,A2.maps...],T[one(T),A2.coeffs...])
end
+(A1::LinearCombination,A2::AbstractLinearMap)=+(A2,A1)
function +(A1::AbstractLinearMap,A2::AbstractLinearMap)
    size(A1)==size(A2) || throw(DimensionMismatch("+"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1,A2],T[one(T),one(T)])
end
function -(A1::LinearCombination,A2::LinearCombination)
    size(A1)==size(A2) || throw(DimensionMismatch("-"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1.maps...,A2.maps...],T[A1.coeffs...,map(-,A2.coeffs)...])
end
function -(A1::AbstractLinearMap,A2::LinearCombination)
    size(A1)==size(A2) || throw(DimensionMismatch("-"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1,A2.maps...],T[one(T),map(-,A2.coeffs)...])
end
function -(A1::LinearCombination,A2::AbstractLinearMap)
    size(A1)==size(A2) || throw(DimensionMismatch("-"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1.maps...,A2],T[A1.coeffs...,-one(T)])
end
function -(A1::AbstractLinearMap,A2::AbstractLinearMap)
    size(A1)==size(A2) || throw(DimensionMismatch("-"))
    T=promote_type(eltype(A1),eltype(A2))
    return LinearCombination{T}(AbstractLinearMap[A1,A2],T[one(T),-one(T)])
end

# scalar multiplication
-(A::AbstractLinearMap)=LinearCombination{eltype(A)}(AbstractLinearMap[A],[-one(eltype(A))])
-(A::LinearCombination)=LinearCombination{eltype(A)}(A.maps,map(-,A.coeffs))

function *(alpha::Number,A::AbstractLinearMap)
    T=promote_type(eltype(alpha),eltype(A))
    return LinearCombination{T}(AbstractLinearMap[A],T[alpha])
end
*(A::AbstractLinearMap,alpha::Number)=*(alpha,A)
function *(alpha::Number,A::LinearCombination)
    T=promote_type(eltype(alpha),eltype(A))
    return LinearCombination{T}(A.maps,alpha*A.coeffs)
end
*(A::LinearCombination,alpha::Number)=*(alpha,A)

function \(alpha::Number,A::AbstractLinearMap)
    T=promote_type(eltype(alpha),eltype(A))
    return LinearCombination{T}(AbstractLinearMap[A],T[one(T)/alpha])
end
/(A::AbstractLinearMap,alpha::Number)=\(alpha,A)
function \(alpha::Number,A::LinearCombination)
    T=promote_type(eltype(alpha),eltype(A))
    return LinearCombination{T}(A.maps,A.coeffs/alpha)
end
/(A::LinearCombination,alpha::Number)=\(alpha,A)

# comparison of LinearCombination objects
==(A::LinearCombination,B::LinearCombination)=(eltype(A)==eltype(B) && A.maps==B.maps && A.coeffs==B.coeffs)

# special transposition behavior
Base.transpose(A::LinearCombination)=LinearCombination{eltype(A)}(AbstractLinearMap[transpose(l) for l in A.maps],A.coeffs)
Base.ctranspose(A::LinearCombination)=LinearCombination{eltype(A)}(AbstractLinearMap[ctranspose(l) for l in A.maps],conj(A.coeffs))

# multiplication with vectors
function Base.A_mul_B!(y::AbstractVector,A::LinearCombination,x::AbstractVector)
    # no size checking, will be done by individual maps
    Base.A_mul_B!(y,A.maps[1],x)
    scale!(A.coeffs[1],y)
    z=similar(y)
    for n=2:length(A.maps)
        Base.A_mul_B!(z,A.maps[n],x)
        Base.axpy!(A.coeffs[n],z,y)
    end
    return y
end
function Base.At_mul_B!(y::AbstractVector,A::LinearCombination,x::AbstractVector)
    # no size checking, will be done by individual maps
    Base.At_mul_B!(y,A.maps[1],x)
    scale!(A.coeffs[1],y)
    z=similar(y)
    for n=2:length(A.maps)
        Base.At_mul_B!(z,A.maps[n],x)
        Base.axpy!(A.coeffs[n],z,y)
    end
    return y
end
function Base.Ac_mul_B!(y::AbstractVector,A::LinearCombination,x::AbstractVector)
    # no size checking, will be done by individual maps
    Base.Ac_mul_B!(y,A.maps[1],x)
    scale!(conj(A.coeffs[1]),y)
    z=similar(y)
    for n=2:length(A.maps)
        Base.Ac_mul_B!(z,A.maps[n],x)
        Base.axpy!(conj(A.coeffs[n]),z,y)
    end
    return y
end


text:
module NURBStoolbox_mixedtest02

using Test;
using MagMechFEM_Matlab2Julia.NURBStoolbox;

@testset "Mixed Tests: 02" begin

  crv = nrbtestcrv();
  evalpoints = vec([0.0 0.25 0.5 0.75 1.0]);
  evalresult = [0.5 3.0 3.75 7.5 8.5;
                3.0 5.5 3.5  1.5 4.5;
                0.0 0.0 0.0  0.0 0.0];

  @test nrbeval(crv,evalpoints) ≈ evalresult;
  @test nrbeval(nrbkntins(crv,[0.5;0.75]),evalpoints) ≈ evalresult;
  @test nrbeval(nrbdegelev(crv,2),evalpoints) ≈ evalresult;
  @test nrbeval(nrbtform(crv,vectrans([2.0;-3.0;5.0])),
                evalpoints) ≈ evalresult .+ [2.0;-3.0;5.0];

end # testset

end # module


text:
######################################################################
# Voidables UTs
# -----
# Licensed under Apache License 2.0
module TestVoidables
using Test
using Pia.Voidables

negate(x) = !x

@testset "Voidables" begin
    @testset "fundamentals" begin
        let voidable = Voidable{Int}()
            @test voidable_param(voidable) === Int
            @test isvoided(voidable)
            @test_throws VoidError voidable[]

            @test (voidable[] = 42) == 42
            @test !isvoided(voidable) && voidable[] == 42
        end

        let voidable = Voidable(42)
            @test voidable_param(voidable) === Int
            @test !isvoided(voidable) && voidable[] == 42

            @test void!(voidable) === void
            @test isvoided(voidable)
            @test_throws VoidError voidable[]
        end
    end

    @testset "multi-ref" begin
        let base = Voidable{Int}(42), ref1 = base, ref2 = base, ref3 = ref1
            @test all(negate.(isvoided.((base, ref1, ref2, ref3))))
            @test void!(base) === void
            @test all(isvoided.((base, ref1, ref2, ref3)))
        end
    end
end

end # module TestVoidables


text:
#=
Suppose an arithmetic expression is given as a binary tree. Each leaf is an integer and each internal node is one of '+', '−', '∗', or '/'.

Given the root to such a tree, write a function to evaluate it.

For example, given the following tree:

    *
   / \
  +    +
 / \  / \
3  2  4  5

You should return 45, as it is (3 + 2) * (4 + 5).
=#

struct Node
    key::Union{Int, Char}
    left::Union{Node, Int}
    right::Union{Node, Int}
end

function evaluate(tree::Node)
    out = 0
    if typeof(tree.left) == Int && typeof(tree.right) == Int
        if tree.key == '+'
            out = tree.left + tree.right
        elseif tree.key == '-'
            out = tree.left - tree.right
        elseif tree.key == '*'
            out = tree.left * tree.right
        elseif tree.key == '/'
            out = tree.left / tree.right
        end
    elseif typeof(tree.left.key) == Char && typeof(tree.right.key) == Char
        if tree.key == '+'
            out = evaluate(tree.left) + evaluate(tree.right)
        elseif tree.key == '-'
            out = evaluate(tree.left) - evaluate(tree.right)
        elseif tree.key == '*'
            out = evaluate(tree.left) * evaluate(tree.right)
        elseif tree.key == '/'
            out = evaluate(tree.left) / evaluate(tree.right)
        end
    end
    return out
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
function pay_cashflows!(investors::Array{T, 1} where {T<:AbstractInvestor}, assets::Array{T, 1} where {T<:AbstractAsset})
    rates = Dict{AbstractAsset, Float64}()
    for asset in assets
        rate = get_rate(asset)
        for investor in investors
            n_assets = investor.positions[asset]
            n_asset_new = n_assets * (1 + rate)
            n_asset_new = rate > 0 ? floor(n_asset_new) : ceil(n_asset_new)
            investor.positions[asset] = Int64(n_asset_new)
        end
        rates[asset] = rate
    end
    rates
end

function pay_cashflow(asset::AbstractAsset, investor::AbstractInvestor)
    rate = get_rate(asset)
    n_assets = investor.positions[asset]
    n_asset_new = n_stocks * (1 - rate)
    n_asset_new = rate > 0 ? floor(n_asset_new) : ceil(n_asset_new)
    inv[payment_currency] += int(n_asset_new)
    return rate
end

"Place holder function to be defined"
function get_rate(asset::AbstractAsset)
    return 0
end


text:
using TextAnalysis
using MutableStrings

const nstr = 10^7
const lstr = 6

function makestr_ascii()
    iob = IOBuffer()
    for idx in 1:nstr
        write(iob, randstring(lstr))
        write(iob, ' ')
    end
    takebuf_string(iob)
end

function makestr_utf8()
    iob = IOBuffer()
    utf8chars = convert(Array{Char,1}, "abcdefΑΒΓΔΕΖΗΘαβγδεζηθ")
    for idx in 1:nstr
        if (idx % 10) == 0
            shuffle!(utf8chars)
            write(iob, convert(UTF8String, utf8chars[1:lstr]))
        else
            write(iob, randstring(lstr))
        end
        write(iob, ' ')
    end
    takebuf_string(iob)
end

function warmup_ascii()
    _mstr = MutableASCIIString("asdsdsf");
    l = length(_mstr)
    lowercase!(_mstr);
    uppercase!(_mstr);
    reverse!(_mstr);
    ucfirst!(_mstr);
    lcfirst!(_mstr);
    replace!(_mstr, "abc", ' ');
    remove_case!(StringDocument("ssdfsdfdsf"));
    TextAnalysis.remove_patterns!(StringDocument("ssdfsdfdsf"), r"abc");
    TextAnalysis.remove_patterns!(StringDocument(MutableASCIIString("ssdfsdfdsf")), r"abc");
    TextAnalysis.stem!(StringDocument("ssdfsdfdsf"));
    TextAnalysis.stem!(StringDocument(MutableASCIIString("ssdfsdfdsf")));
end

function warmup_utf8()
    _mstr = MutableUTF8String("abcdefΑΒΓΔΕΖΗΘαβγδεζηθ");
    l = length(_mstr)
    lowercase!(_mstr);
    uppercase!(_mstr);
    reverse!(_mstr);
    ucfirst!(_mstr);
    lcfirst!(_mstr);
    replace!(_mstr, "abc", ' ');
    remove_case!(StringDocument("abcdefΑΒΓΔΕΖΗΘαβγδεζηθ"));
    TextAnalysis.remove_patterns!(StringDocument("abcdefΑΒΓΔΕΖΗΘαβγδεζηθ"), r"abc");
    TextAnalysis.remove_patterns!(StringDocument(MutableUTF8String("abcdefΑΒΓΔΕΖΗΘαβγδεζηθ")), r"abc");
    TextAnalysis.stem!(StringDocument("abcdefΑΒΓΔΕΖΗΘαβγδεζηθ"));
    TextAnalysis.stem!(StringDocument(MutableUTF8String("abcdefΑΒΓΔΕΖΗΘαβγδεζηθ")));
end

function test_ascii()
    str = makestr_ascii();
    mstr = MutableASCIIString(copy(str.data));

    strdoc = StringDocument(str);
    mstrdoc = StringDocument(mstr);

    @printf("%20s%30s%40s\n", "", "string", "mutable string")
    @printf("%20s%30s%40s\n", "length", "$(length(str))", "$(length(mstr))")
    @printf("%20s%20s%20s%20s%20s\n", "", "time", "bytes", "time", "bytes")

    gc(); _ret, t, b = @timed remove_case!(strdoc);
    gc(); _mret, mt, mb = @timed remove_case!(mstrdoc)
    @printf("%20s%20s%20s%20s%20s\n", "remove_case!", "$t", "$b", "$mt", "$mb")

    strdoc = StringDocument(str)
    mstr = MutableASCIIString(copy(str.data));
    mstrdoc = StringDocument(mstr)

    gc(); _ret, t, b = @timed prepare!(strdoc, strip_case)
    gc(); _mret, mt, mb = @timed prepare!(mstrdoc, strip_case)
    @printf("%20s%20s%20s%20s%20s\n", "prepare!->strip_case", "$t", "$b", "$mt", "$mb")

    rx = r"a[a-z]+\s"
    gc(); _ret, t, b = @timed TextAnalysis.remove_patterns!(strdoc, rx)
    gc(); _mret, mt, mb = @timed TextAnalysis.remove_patterns!(mstrdoc, rx)
    @printf("%20s%20s%20s%20s%20s\n", "remove_patterns!", "$t", "$b", "$mt", "$mb")

    stm = Stemmer("english")
    str = "running skipping " ^ 10^5
    mstr = MutableASCIIString(copy(str.data))
    strdoc = StringDocument(str)
    mstrdoc = StringDocument(mstr)
    gc(); _ret, t, b = @timed TextAnalysis.stem!(stm, strdoc)
    gc(); _mret, mt, mb = @timed TextAnalysis.stem!(stm, mstrdoc)
    @printf("%20s%20s%20s%20s%20s\n", "stem!", "$t", "$b", "$mt", "$mb")
end

function test_utf8()
    str = makestr_utf8();
    mstr = MutableUTF8String(copy(str.data));

    strdoc = StringDocument(str);
    mstrdoc = StringDocument(mstr);

    @printf("%20s%30s%40s\n", "", "string", "mutable string")
    @printf("%20s%30s%40s\n", "length", "$(length(str))", "$(length(mstr))")
    @printf("%20s%20s%20s%20s%20s\n", "", "time", "bytes", "time", "bytes")

    gc(); _ret, t, b = @timed remove_case!(strdoc);
    gc(); _mret, mt, mb = @timed remove_case!(mstrdoc)
    @printf("%20s%20s%20s%20s%20s\n", "remove_case!", "$t", "$b", "$mt", "$mb")

    strdoc = StringDocument(str)
    mstr = MutableUTF8String(copy(str.data));
    mstrdoc = StringDocument(mstr)

    gc(); _ret, t, b = @timed prepare!(strdoc, strip_case)
    gc(); _mret, mt, mb = @timed prepare!(mstrdoc, strip_case)
    @printf("%20s%20s%20s%20s%20s\n", "prepare!->strip_case", "$t", "$b", "$mt", "$mb")

    rx = r"a[a-z]+\s"
    gc(); _ret, t, b = @timed TextAnalysis.remove_patterns!(strdoc, rx)
    gc(); _mret, mt, mb = @timed TextAnalysis.remove_patterns!(mstrdoc, rx)
    @printf("%20s%20s%20s%20s%20s\n", "remove_patterns!", "$t", "$b", "$mt", "$mb")

    stm = Stemmer("english")
    str = "running skipping " ^ 10^5
    mstr = MutableUTF8String(copy(str.data))
    strdoc = StringDocument(str)
    mstrdoc = StringDocument(mstr)
    gc(); _ret, t, b = @timed TextAnalysis.stem!(stm, strdoc)
    gc(); _mret, mt, mb = @timed TextAnalysis.stem!(stm, mstrdoc)
    @printf("%20s%20s%20s%20s%20s\n", "stem!", "$t", "$b", "$mt", "$mb")
end

warmup_ascii()
test_ascii()

warmup_utf8()
test_utf8()


text:
# ---
# title: 1201. Ugly Number III
# id: problem1201
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/ugly-number-iii/description/>
# hidden: true
# ---
#
# Write a program to find the `n`-th ugly number.
#
# Ugly numbers are **  positive integers** which are divisible by `a`  **or**
# `b`  **or** `c`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3, a = 2, b = 3, c = 5
#     Output: 4
#     Explanation: The ugly numbers are 2, 3, 4, 5, 6, 8, 9, 10... The 3rd is 4.
#
# **Example 2:**
#
#
#
#     Input: n = 4, a = 2, b = 3, c = 4
#     Output: 6
#     Explanation: The ugly numbers are 2, 3, 4, 6, 8, 9, 10, 12... The 4th is 6.
#
#
# **Example 3:**
#
#
#
#     Input: n = 5, a = 2, b = 11, c = 13
#     Output: 10
#     Explanation: The ugly numbers are 2, 4, 6, 8, 10, 11, 12, 13... The 5th is 10.
#
#
# **Example 4:**
#
#
#
#     Input: n = 1000000000, a = 2, b = 217983653, c = 336916467
#     Output: 1999999984
#
#
#
#
# **Constraints:**
#
#   * `1 <= n, a, b, c <= 10^9`
#   * `1 <= a * b * c <= 10^18`
#   * It's guaranteed that the result will be in range `[1, 2 * 10^9]`
#
#
## @lc code=start
using LeetCode

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


text:
# Contiguous on-device arrays

export
    CuDeviceArray, CuBoundsError

## construction

"""
    CuDeviceArray(dims, ptr)
    CuDeviceArray{T}(dims, ptr)
    CuDeviceArray{T,N}(dims, ptr)

Construct an `N`-dimensional dense CUDA device array with element type `T` wrapping a
pointer, where `N` is determined from the length of `dims` and `T` is determined from the
type of `ptr`. `dims` may be a single scalar, or a tuple of integers corresponding to the
lengths in each dimension). If the rank `N` is supplied explicitly as in `Array{T,N}(dims)`,
then it must match the length of `dims`. The same applies to the element type `T`, which
should match the type of the pointer `ptr`.
"""
CuDeviceArray

# NOTE: we can't support the typical `tuple or series of integer` style construction,
#       because we're currently requiring a trailing pointer argument.

struct CuDeviceArray{T,N} <: AbstractArray{T,N}
    shape::NTuple{N,Int}
    ptr::Ptr{T}

    # inner constructors (exact types, ie. Int not <:Integer)
    CuDeviceArray{T,N}(shape::NTuple{N,Int}, ptr::Ptr{T}) where {T,N} = new(shape, ptr)
end

# outer constructors, non-parameterized
CuDeviceArray(dims::NTuple{N,<:Integer}, p::Ptr{T})                where {T,N} = CuDeviceArray{T,N}(dims, p)
CuDeviceArray(len::Integer, p::Ptr{T})                             where {T}   = CuDeviceArray{T,1}((len,), p)

# outer constructors, partially parameterized
(::Type{CuDeviceArray{T}})(dims::NTuple{N,<:Integer}, p::Ptr{T})   where {T,N} = CuDeviceArray{T,N}(dims, p)
(::Type{CuDeviceArray{T}})(len::Integer, p::Ptr{T})                where {T}   = CuDeviceArray{T,1}((len,), p)

# outer constructors, fully parameterized
(::Type{CuDeviceArray{T,N}})(dims::NTuple{N,<:Integer}, p::Ptr{T}) where {T,N} = CuDeviceArray{T,N}(Int.(dims), p)
(::Type{CuDeviceArray{T,1}})(len::Integer, p::Ptr{T})              where {T}   = CuDeviceArray{T,1}((Int(len),), p)

Base.convert(::Type{CuDeviceArray{T,N}}, a::CuArray{T,N}) where {T,N} =
    CuDeviceArray{T,N}(a.shape, Base.unsafe_convert(Ptr{T}, a.devptr))

Base.unsafe_convert(::Type{Ptr{T}}, a::CuDeviceArray{T}) where {T} = a.ptr::Ptr{T}

## array interface

Base.size(g::CuDeviceArray) = g.shape
Base.length(g::CuDeviceArray) = prod(g.shape)

@inline function Base.getindex(A::CuDeviceArray{T}, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    Base.pointerref(Base.unsafe_convert(Ptr{T}, A), index, 8)::T
end

@inline function Base.setindex!(A::CuDeviceArray{T}, x, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    Base.pointerset(Base.unsafe_convert(Ptr{T}, A), convert(T, x)::T, index, 8)
end

Base.IndexStyle(::Type{<:CuDeviceArray}) = Base.IndexLinear()

Base.show(io::IO, a::CuDeviceArray{T,1}) where {T} =
    print(io, "$(length(a))-element device array at $(pointer(a))")
Base.show(io::IO, a::CuDeviceArray{T,N}) where {T,N} =
    print(io, "$(join(a.shape, '×')) device array at $(pointer(a))")

## quirks

# bounds checking is currently broken due to a PTX assembler issue (see #4)
Base.checkbounds(::CuDeviceArray, I...) = nothing

# replace boundserror-with-arguments to a non-allocating, argumentless version
# TODO: can this be fixed by stack-allocating immutables with heap references?
struct CuBoundsError <: Exception end
@inline Base.throw_boundserror(A::CuDeviceArray, I) =
    (Base.@_noinline_meta; throw(CuBoundsError()))

# idem
function Base.unsafe_view(A::CuDeviceArray{T,1}, I::Vararg{Base.ViewIndex,1}) where {T}
    Base.@_inline_meta
    ptr = Base.unsafe_convert(Ptr{T}, A) + (I[1].start-1)*sizeof(T)
    len = I[1].stop - I[1].start + 1
    return CuDeviceArray(len, ptr)
end


text:
@testset "447.number-of-boomerangs.jl" begin
    @test number_of_boomerangs([[0, 0], [1, 0], [2, 0]]) == 2
    @test number_of_boomerangs([[1, 1], [2, 2], [3, 3]]) == 2
    @test number_of_boomerangs([[1, 1]]) == 0
end


text:
using fci
using JLD2

@load "_testdata_h4_triplet_integrals.jld2"
@load "_testdata_h4_triplet_alphaconfigs.jld2"
configs = alpha_configs
int1e = one
int2e = two
y_matrix = ya
norbs = 4
nalpha = 3
ndets = factorial(norbs)÷(factorial(nalpha)*factorial(norbs-nalpha))
index_table = index_table_a

@testset "ss terms" begin
    Ha_test = zeros(size(H_alpha))
    nelecs = size(configs[1].config)[1]

    for I in configs #Ia or Ib, configs=list of all possible determinants
        I_idx = fci.get_index(I.config, y_matrix, I.norbs)
        F = zeros(ndets)
        orbs = [1:I.norbs;]
        vir = filter!(x->!(x in I.config), orbs)

        #single excitation
        for k in I.config
            for l in vir
                #annihlate electron in orb k
                config_single, sorted_config, sign_s = fci.excit_config(deepcopy(I.config), [k,l])
                config_single_idx = index_table[k,l,I.label]
                F[abs(config_single_idx)] += sign_s*int1e[k,l]
                for m in I.config
                    if m!=k
                        F[abs(config_single_idx)] += sign_s*(int2e[k,l,m,m]-int2e[k,m,l,m])
                    end
                end
            end
        end

        #double excitation
        for k in I.config
            for i in I.config
                if i>k
                    for l in vir
                        for j in vir
                            if j>l
                                single, sorted_s, sign_s = fci.excit_config(deepcopy(I.config), [k,l])
                                double, sorted_d, sign_d = fci.excit_config(deepcopy(sorted_s), [i,j])
                                idx = get_index(double, y_matrix, I.norbs)
                                if sign_d == sign_s
                                    F[abs(idx)] += (int2e[i,j,k,l] - int2e[i,l,j,k]) #one that works for closed

                                else
                                    F[abs(idx)] -= (int2e[i,j,k,l] - int2e[i,l,j,k]) #one that works for closed
                                end
                            end
                        end
                    end
                end
            end
        end

        Ha_test[:,I_idx] .= F
    end

    @test isapprox(H_alpha, Ha_test, atol=0.05)
end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci


text:
module Events

struct Emitter
    listeners::Dict{Symbol,Vector{Function}}
end

Emitter() = Emitter(Dict())

add!(emitter, event::Symbol, f::Function) = begin
    if !haskey(emitter.listeners, event)
        emitter.listeners[event] = []
    end
    push!(emitter.listeners[event], f)
    emitter
end

listeners(emitter, event::Symbol) = haskey(emitter.listeners, event) ?
    emitter.listeners[event] :
    []

on!(f::Function, emitter, event::Symbol) = begin
    add!(emitter, event, f)
    () -> off!(emitter, event, f)
end

off!(emitter, event::Symbol, f::Function) = begin
    if haskey(emitter.listeners, event)
        filter!(emitter.listeners[event]) do f′
            f′ ≠ f
        end
    end
    emitter
end

off!(emitter, event::Symbol) = begin
    delete!(emitter.listeners, event)
    emitter
end

off!(emitter) = begin
    empty!(emitter.listeners)
    emitter
end

emit(emitter, event::Symbol, args...) = begin
    for listener ∈ listeners(emitter, event)
        listener(args...)
    end
end

export
    Emitter,
    add!,
    listeners,
    on!,
    off!,
    emit

end


text:
############################
# Types & Constructors
############################
Base.@kwdef mutable struct DefaultStatus <: Status
    feasible::Bool = false;
    optimal::Bool  = false;
    objfunction::Float64 = 0.0;
    constraintLhsConsumption::Array{Float64,1};
end

"""
    constructStatus(nConstraints)

Constructs a `DefaultStatus`. It is necessary to specify the number of constraints of the problem.

# Example
```jldoctest
julia> OREnvironment.constructStatus(5);
```
"""
function constructStatus(numberConstraints::Int)
    return DefaultStatus(constraintLhsConsumption = zeros(Float64, numberConstraints));
end

############################
# General methods for Status
############################
"""
    is_feasible(s)

Returns the feasibility value of status `s`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.is_feasible(s)
false
```
"""
is_feasible(s::Status)::Bool = s.feasible;

"""
    set_feasible!(s, val)

Sets the feasibility value of status `s` to `val`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.set_feasible!(s, true);
julia> OREnvironment.is_feasible(s)
true
```
"""
set_feasible!(s::Status, value::Bool) = s.feasible = value;

"""
    is_optimal(s)

Returns the optimality value of status `s`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.is_optimal(s)
false
```
"""
is_optimal(s::Status)::Bool = s.optimal;

"""
    set_optimal!(s, val)

Sets the optimality value of status `s` to `val`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.set_optimal!(s, true);
julia> OREnvironment.is_optimal(s)
true
```
"""
set_optimal!(s::Status, value::Bool) = s.optimal = value;

"""
    get_objfunction(s)

Returns the objective function value of status `s`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.get_objfunction(s)
0.0
```
"""
get_objfunction(s::Status)::Float64 = s.objfunction;

"""
    set_objfunction!(s, val)

Sets the objective function value of status `s` to `val`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.set_objfunction!(s, 19.8);
julia> OREnvironment.get_objfunction(s)
19.8
```
"""
set_objfunction!(s::Status, value::Float64) = s.objfunction = value;

"""
    get_constraint_consumption(s, idx)

Returns the consumption of the `idx`-th constraint in status `s`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.get_constraint_consumption(s, 1)
0.0
```
"""
@inline function get_constraint_consumption(s::Status, idxConstraint::Int)::Float64
    @inbounds s.constraintLhsConsumption[idxConstraint];
end

"""
    set_constraint_consumption!(s, idx, val)

Sets the consumption of the `idx`-th constraint to `val` in status `s`.

# Example
```jldoctest
julia> s = OREnvironment.constructStatus(5);
julia> OREnvironment.set_constraint_consumption!(s, 2, 3.4);
julia> OREnvironment.get_constraint_consumption(s, 2)
3.4
```
"""
@inline function set_constraint_consumption!(s::Status, idxConstraint::Int, value::Float64)
    @inbounds s.constraintLhsConsumption[idxConstraint] = value;
end

"""
    worst_value(objSense)

Returns the worst value of the objetive function according to the objective sense of the problem.

# Example
```jldoctest
julia> OREnvironment.worst_value(:max)
-Inf
julia> OREnvironment.worst_value(:min)
Inf
```
"""
worst_value(objsense::Symbol)::Float64 = worst_value(Val(objsense));
worst_value(::Val{:max})::Float64 = typemin(Float64);
worst_value(::Val{:min})::Float64 = typemax(Float64);

"""
    is_first_status_better(s1, s2, objSense, feasibilityRequiered)

Returns if status `s1` is better than status `s2`. For the comparaison is required the objective function sense and if feasibility is a must for the comparaison.

# Example
```jldoctest
julia> s1 = OREnvironment.constructStatus(5);
julia> OREnvironment.set_feasible!(s1, true);
julia> OREnvironment.set_objfunction!(s1, 19.8);
julia> s2 = OREnvironment.constructStatus(5);
julia> OREnvironment.set_feasible!(s2, false);
julia> OREnvironment.set_objfunction!(s2, 100.0);
julia> OREnvironment.is_first_status_better(s1, s2, :max, true)
true
julia> OREnvironment.is_first_status_better(s1, s2, :max, false)
false
julia> OREnvironment.set_feasible!(s2, true);
julia> OREnvironment.is_first_status_better(s1, s2, :max, true)
false
julia> OREnvironment.is_first_status_better(s1, s2, :min, true)
true
```
"""
function is_first_status_better(s1::Status,
                                s2::Status,
                                objSense::Symbol,
                                feasibilityRequiered::Bool)::Bool
    if feasibilityRequiered
        if is_feasible(s1) && !is_feasible(s2)
            return true;
        elseif !is_feasible(s1) && is_feasible(s2)
            return false;
        end
    end
    # In any other situation we compare the objective function
    return is_first_obj_function_better(s1, s2, Val(objSense));
end

@inline function is_first_obj_function_better(s1::Status, s2::Status, ::Val{:max})::Bool
    get_objfunction(s1) > get_objfunction(s2);
end

@inline function is_first_obj_function_better(s1::Status, s2::Status, ::Val{:min})::Bool
    get_objfunction(s1) < get_objfunction(s2);
end

############################
# Specific methods for Status
############################
"""
    update_status!(s1, s2)

Copies info of status `s2` to `s1`.

# Example
```jldoctest
julia> s1 = OREnvironment.constructStatus(5);
julia> s2 = OREnvironment.constructStatus(5);
julia> OREnvironment.set_feasible!(s2, true);
julia> OREnvironment.set_optimal!(s2, true);
julia> OREnvironment.set_objfunction!(s2, 100.0);
julia> OREnvironment.is_feasible(s1)
false
julia> OREnvironment.is_optimal(s1)
false
julia> OREnvironment.get_objfunction(s1)
0.0
julia> OREnvironment.update_status!(s1, s2);
julia> OREnvironment.is_feasible(s1)
true
julia> OREnvironment.is_optimal(s1)
true
julia> OREnvironment.get_objfunction(s1)
100.0
```
"""
function update_status!(oldStatus::DefaultStatus,
                        newStatus::DefaultStatus)
    set_objfunction!(oldStatus, get_objfunction(newStatus));
    set_feasible!(oldStatus, is_feasible(newStatus));
    set_optimal!(oldStatus, is_optimal(newStatus));
    copy!(oldStatus.constraintLhsConsumption, newStatus.constraintLhsConsumption);
    return nothing;
end


text:
using Test, VarBCFiles

@testset "Only version 6" begin
   # Check that we only allow version 6 files
   @test_throws AssertionError("version == \"VARBC_cycle.version006\"") read("data/VARBC.cycle_v5",VarBC)
end

v1 = read("data/VARBC.cycle1",VarBC)
v2 = read("data/VARBC.cycle2",VarBC)

# Check that we write back file exactly as we read them
@testset "write" begin
     file, = mktemp()
     write(file,v2)
     @test read(file) == read("data/VARBC.cycle2") # this reads files as Array{UInt8,1}
     rm(file)
end

@testset "filter!" begin
    v = read("MEPSprod/VARBC.cycle",VarBC)
    @test length(filter!(p-> p.second.class=="sfcobs",v)) == 95
end

@testset "merge!" begin
    v1 = read("data/VARBC.cycle1_merge",VarBC)
    v2 = read("data/VARBC.cycle2_merge",VarBC)
    v3 = read("data/VARBC.cycle_merge2into1",VarBC)
    merge!(v1,v2)
    @test v3.records == v1.records
end

@testset "Equality" begin
    @test v1 != v2
end

@testset "Conversions" begin
    # Check that rmdi2nan works
    key="D0MZ8zJj 16199110       1     6"
    @test isnan(v2.records[key].param0[1])
end


text:
using JLD2
using Plots
using Printf
using Oceananigans
using Oceananigans.Grids

filename = "oceananigans_unstable_bickley_jet_Nh512_WENO5_fields.jld2"

file = jldopen(filename)

iterations = parse.(Int, keys(file["timeseries/t"]))
grid = file["serialized/grid"]

xc, yc, zc = nodes((Cell, Cell, Cell), grid)
xω, yω, zω = nodes((Face, Face, Cell), grid)

snapshots = []

for iteration in iterations[1:6:54]

    t = file["timeseries/t/$iteration"]
    ω = file["timeseries/ω/$iteration"][:, :, 1]
    c = file["timeseries/c/$iteration"][:, :, 1]

    kwargs = Dict(
                  :aspectratio => 1,
                  :linewidth => 0,
                  :colorbar => false,
                  :xticks => false,
                  :yticks => false,
                  :clims => (-1, 1),
                  :xlims => (-grid.Lx/2, grid.Lx/2),
                  :ylims => (-grid.Ly/2, grid.Ly/2),
                  :title => @sprintf("t = %.1f", t)
                 )

    #snap = heatmap(xc, yc, clamp.(c, -1, 1)'; color = :thermal, kwargs...)
    snap = heatmap(xω, yω, clamp.(ω, -1, 1)'; color = :balance, kwargs...)

    push!(snapshots, snap)
end

close(file)

plt = plot(snapshots..., size = (1800, 1600))

display(plt)

savefig("vorticity.png")
#savefig("dye.png")


text:
"""
Constructs a Differentiable Optimizer model from a MOI Optimizer.
Supports `forward` and `backward` methods for solving and differentiating the model respectectively.

## Note
Currently supports differentiating linear and quadratic programs only.
"""

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

const SUPPORTED_OBJECTIVES = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64},
    MOI.ScalarQuadraticFunction{Float64}
}

const SUPPORTED_SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64}
}

const SUPPORTED_SCALAR_FUNCTIONS = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64}
}

const SUPPORTED_VECTOR_FUNCTIONS = Union{
    MOI.VectorOfVariables,
    MOI.VectorAffineFunction{Float64},
}

const SUPPORTED_VECTOR_SETS = Union{
    MOI.Zeros,
    MOI.Nonpositives,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    MOI.PositiveSemidefiniteConeTriangle,
}

"""
    diff_optimizer(optimizer_constructor)::Optimizer

Creates a `DiffOpt.Optimizer`, which is an MOI layer with an internal optimizer
and other utility methods. Results (primal, dual and slack values) are obtained
by querying the internal optimizer instantiated using the
`optimizer_constructor`. These values are required for find jacobians with respect to problem data.

One define a differentiable model by using any solver of choice. Example:

```julia
julia> using DiffOpt, GLPK

julia> model = diff_optimizer(GLPK.Optimizer)
julia> model.add_variable(x)
julia> model.add_constraint(...)

julia> backward!(model)  # for convex quadratic models

julia> backward!(model)  # for convex conic models
```
"""
function diff_optimizer(optimizer_constructor)::Optimizer
    return Optimizer(MOI.instantiate(optimizer_constructor, with_bridge_type=Float64))
end

mutable struct Optimizer{OT <: MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT
    primal_optimal::Vector{Float64}  # solution
    dual_optimal::Vector{Union{Vector{Float64}, Float64}}  # refer - https://github.com/jump-dev/DiffOpt.jl/issues/21#issuecomment-651130485
    var_idx::Vector{VI}
    con_idx::Vector{CI}

    function Optimizer(optimizer_constructor::OT) where {OT <: MOI.ModelLike}
        new{OT}(
            optimizer_constructor,
            zeros(0),
            zeros(0),
            Vector{VI}(),
            Vector{CI}(),
        )
    end
end

function MOI.add_variable(model::Optimizer)
    vi = MOI.add_variable(model.optimizer)
    push!(model.var_idx, vi)
    return vi
end

function MOI.add_variables(model::Optimizer, N::Int)
    return VI[MOI.add_variable(model) for i in 1:N]
end

function MOI.add_constraint(model::Optimizer, f::SUPPORTED_SCALAR_FUNCTIONS, s::SUPPORTED_SCALAR_SETS)
    ci = MOI.add_constraint(model.optimizer, f, s)
    push!(model.con_idx, ci)
    return ci
end

function MOI.add_constraint(model::Optimizer, vf::SUPPORTED_VECTOR_FUNCTIONS, s::SUPPORTED_VECTOR_SETS)
    ci = MOI.add_constraint(model.optimizer, vf, s)
    push!(model.con_idx, ci)   # adding it as a whole here; need to define a method to extract matrices
    return ci
end

function MOI.add_constraints(model::Optimizer, f::AbstractVector{F}, s::AbstractVector{S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return CI{F, S}[MOI.add_constraint(model, f[i], s[i]) for i in eachindex(f)]
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES}, f::SUPPORTED_OBJECTIVES)
    MOI.set(model.optimizer, attr, f)
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    return MOI.set(model.optimizer, attr, sense)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.set(model.optimizer,attr,ci,s)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.set(model.optimizer,attr,ci,f)
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.set(model.optimizer,attr,ci,s)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.set(model.optimizer, attr, ci, f)
end

function MOI.optimize!(model::Optimizer)
    MOI.optimize!(model.optimizer)

    # do not fail. interferes with MOI.Tests.linear12test
    if !in(MOI.get(model.optimizer, MOI.TerminationStatus()),  (MOI.LOCALLY_SOLVED, MOI.OPTIMAL))
        @warn "problem status: $(MOI.get(model.optimizer, MOI.TerminationStatus()))"
        return
    end
    # save the solution
    model.primal_optimal = MOI.get(model.optimizer, MOI.VariablePrimal(), model.var_idx)
    model.dual_optimal = MOI.get.(model.optimizer, MOI.ConstraintDual(), model.con_idx)
    return
end

# """
#     Method to differentiate and obtain gradients/jacobians
#     of z, λ, ν  with respect to the parameters specified in
#     in argument
# """
# function backward(params)
#     grads = []
#     LHS = create_LHS_matrix(z, λ, Q, G, h, A)

#     # compute the jacobian of (z, λ, ν) with respect to each
#     # of the parameters recieved in the method argument
#     # for instance, to get the jacobians w.r.t vector `b`
#     # substitute db = I and set all other differential terms
#     # in the right hand side to zero. For more info refer
#     # equation (6) of https://arxiv.org/pdf/1703.00443.pdf
#     for param in params
#         if param == "Q"
#             RHS = create_RHS_matrix(z, ones(nz, nz), zeros(nz, 1),
#                                     λ, zeros(nineq, nz), zeros(nineq,1),
#                                     ν, zeros(neq, nz), zeros(neq, 1))
#             push!(grads, LHS \ RHS)
#         elseif param == "q"
#             RHS = create_RHS_matrix(z, zeros(nz, nz), ones(nz, 1),
#                                     λ, zeros(nineq, nz), zeros(nineq,1),
#                                     ν, zeros(neq, nz), zeros(neq, 1))
#             push!(grads, LHS \ RHS)
#         elseif param == "G"
#             RHS = create_RHS_matrix(z, zeros(nz, nz), zeros(nz, 1),
#                                     λ, ones(nineq, nz), zeros(nineq,1),
#                                     ν, zeros(neq, nz), zeros(neq, 1))
#             push!(grads, LHS \ RHS)
#         elseif param == "h"
#             RHS = create_RHS_matrix(z, zeros(nz, nz), zeros(nz, 1),
#                                     λ, zeros(nineq, nz), ones(nineq,1),
#                                     ν, zeros(neq, nz), zeros(neq, 1))
#             push!(grads, LHS \ RHS)
#         elseif param == "A"
#             RHS = create_RHS_matrix(z, zeros(nz, nz), zeros(nz, 1),
#                                     λ, zeros(nineq, nz), zeros(nineq,1),
#                                     ν, ones(neq, nz), zeros(neq, 1))
#             push!(grads, LHS \ RHS)
#         elseif param == "b"
#             RHS = create_RHS_matrix(z, zeros(nz, nz), zeros(nz, 1),
#                                     λ, zeros(nineq, nz), zeros(nineq,1),
#                                     ν, zeros(neq, nz), ones(neq, 1))
#             push!(grads, LHS \ RHS)
#         else
#             push!(grads, [])
#         end
#     end
#     return grads
# end

"""
    backward!(model::Optimizer, params::Array{String}, dl_dz::Array{Float64})

Method to differentiate optimal solution `z` and return
product of jacobian matrices (`dz / dQ`, `dz / dq`, etc) with
the backward pass vector `dl / dz`

The method computes the product of
1. jacobian of problem solution `z*` with respect to
    problem parameters `params` recieved as method arguments
2. a backward pass vector `dl / dz`, where `l` can be a loss function

Note that this method *does not returns* the actual jacobians.

For more info refer eqn(7) and eqn(8) of https://arxiv.org/pdf/1703.00443.pdf
"""
function backward!(model::Optimizer, params::Vector{String}, dl_dz::Vector{Float64})
    Q, q, G, h, A, b, nz, var_idx, nineq, ineq_con_idx, neq, eq_con_idx = get_problem_data(model.optimizer)

    z = model.primal_optimal

    # separate λ, ν
    λ = filter(con -> !is_equality(MOI.get(model.optimizer, MOI.ConstraintSet(), con)), model.con_idx)
    ν = filter(con -> is_equality(MOI.get(model.optimizer, MOI.ConstraintSet(), con)), model.con_idx)

    λ = [MOI.get(model.optimizer, MOI.ConstraintDual(), con) for con in λ]
    ν = [MOI.get(model.optimizer, MOI.ConstraintDual(), con) for con in ν]

    grads = []
    LHS = create_LHS_matrix(z, λ, Q, G, h, A)
    RHS = [dl_dz; zeros(neq+nineq)]

    partial_grads = -(LHS \ RHS)

    dz = partial_grads[1:nz]
    if nineq > 0
        dλ = partial_grads[nz+1:nz+nineq]
    end
    if neq > 0
        dν = partial_grads[nz+nineq+1:nz+nineq+neq]
    end

    for param in params
        if param == "Q"
            push!(grads, 0.5 * (dz * z' + z * dz'))
        elseif param == "q"
            push!(grads, dz)
        elseif param == "G"
            push!(grads, Diagonal(λ) * dλ * z' - λ * dz')
        elseif param == "h"
            push!(grads, -Diagonal(λ) * dλ)
        elseif param == "A"
            push!(grads, dν * z' - ν * dz')
        elseif param == "b"
            push!(grads, -dν)
        else
            push!(grads, [])
        end
    end
    return grads
end

# `MOI.supports` methods

function MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}})
    return true
end

function MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}})
    return true
end

function MOI.supports(::Optimizer, ::MOI.AbstractModelAttribute)
    return true
end

function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction)
    return false
end

function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES})
    return true
end

MOI.supports_constraint(::Optimizer, ::Type{<: SUPPORTED_SCALAR_FUNCTIONS}, ::Type{<: SUPPORTED_SCALAR_SETS}) = true
MOI.supports_constraint(::Optimizer, ::Type{<: SUPPORTED_VECTOR_FUNCTIONS}, ::Type{<: SUPPORTED_VECTOR_SETS}) = true
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS} = true
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS} = true

MOI.get(model::Optimizer, attr::MOI.SolveTime) = MOI.get(model.optimizer, attr)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.optimizer)
    empty!(model.primal_optimal)
    empty!(model.dual_optimal)
    empty!(model.var_idx)
    empty!(model.con_idx)
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           isempty(model.primal_optimal) &&
           isempty(model.dual_optimal) &&
           isempty(model.var_idx) &&
           isempty(model.con_idx)
end

# now supports name too
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = true #!copy_names

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    return MOIU.default_copy_to(model.optimizer, src, copy_names)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return MOI.get(model.optimizer, MOI.TerminationStatus())
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::VI, value::Float64)
    MOI.set(model.optimizer, MOI.VariablePrimalStart(), vi, value)
end

function MOI.supports(model::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return MOI.supports(model.optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(model.optimizer, attr)
    return MOI.get(model.optimizer, attr, vi)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    filter!(x -> x ≠ ci, model.con_idx)
    MOI.delete(model.optimizer, ci)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    filter!(x -> x ≠ ci, model.con_idx)
    MOI.delete(model.optimizer, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.is_valid(model::Optimizer, v::VI)
    return v in model.var_idx
end

function MOI.is_valid(model::Optimizer, con::CI)
    return con in model.con_idx
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintBasisStatus, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintBasisStatus(), ci)
end

# helper methods to check if a constraint contains a Variable
function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.SingleVariable, S}) where {S <: SUPPORTED_SCALAR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v == func.variable
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.ScalarAffineFunction{Float64}, S}) where {S <: SUPPORTED_SCALAR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.variable_index, func.terms)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorOfVariables, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v in func.variables
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorAffineFunction{Float64}, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.scalar_term.variable_index, func.terms)
end

function MOI.delete(model::Optimizer, v::VI)
    # remove those constraints that depend on this Variable
    filter!(ci -> !_constraint_contains(model, v, ci), model.con_idx)

    # delete in inner solver
    MOI.delete(model.optimizer, v)

    # delete from var_idx
    filter!(x -> x ≠ v, model.var_idx)
end

# for array deletion
function MOI.delete(model::Optimizer, indices::Vector{VI})
    for i in indices
        MOI.delete(model, i)
    end
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.AbstractFunctionModification
)
    MOI.modify(
        model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        chg
    )
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:MOI.VectorAffineFunction{Float64}, S <: SUPPORTED_VECTOR_SETS}
    MOI.modify(model.optimizer, ci, chg)
end

"""
    π(cones::Array{<: MOI.AbstractVectorSet}, v::Vector{Vector{Float64}})

Find projection of vectors in `v` on product of `cones`.
For more info, refer https://github.com/matbesancon/MathOptSetDistances.jl
"""
π(cones, v) = MOSD.projection_on_set(MOSD.DefaultDistance(), v, cones)

"""
    Dπ(cones::Array{<: MOI.AbstractVectorSet}, v::Vector{Vector{Float64}})

Find gradient of projection of vectors in `v` on product of `cones`.
For more info, refer https://github.com/matbesancon/MathOptSetDistances.jl
"""
Dπ(cones, v) = MOSD.projection_gradient_on_set(MOSD.DefaultDistance(), v, cones)

"""
    backward_conic!(model::Optimizer, dA::Array{Float64,2}, db::Array{Float64}, dc::Array{Float64})

Method to differentiate optimal solution `x`, `y`, `s` given perturbations related to
conic program parameters `A`, `b`, `c`. This is similar to [`backward!`](@ref) method
but it this *does returns* the actual jacobians.

For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043
"""
function backward_conic!(model::Optimizer, dA::Matrix{Float64}, db::Vector{Float64}, dc::Vector{Float64})
    if !in(
            MOI.get(model, MOI.TerminationStatus()), (MOI.LOCALLY_SOLVED, MOI.OPTIMAL)
        )
        error("problem status: ", MOI.get(model.optimizer, MOI.TerminationStatus()))
    end

    # fetch matrices from MatrixOptInterface
    conic_form = MatOI.get_conic_form(Float64, model.optimizer, model.con_idx)

    A = conic_form.A
    A = CSRToCSC(A)
    b = conic_form.b
    c = conic_form.c

    # programs in tests were cross-checked against `diffcp`, which follows SCS format
    # hence, some arrays saved during `MOI.optimize!` are not same across all optimizers
    # specifically there's an extra preprocessing step for `PositiveSemidefiniteConeTriangle` constraint for SCS/Mosek

    # get x,y,s
    x = model.primal_optimal
    s = MOI.get(model, MOI.ConstraintPrimal(), model.con_idx)
    y = model.dual_optimal

    # pre-compute quantities for the derivative
    m = A.m
    n = A.n
    N = m + n + 1
    cones = MOI.get(model, MOI.ConstraintSet(), model.con_idx)
    # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
    (u, v, w) = (x, y - s, 1.0)

    # find gradient of projections on dual of the cones
    Dπv = Dπ(MOI.dual_set.(cones), v)

    # Q = [
    #      0   A'   c;
    #     -A   0    b;
    #     -c' -b'   0;
    # ]
    # M = (Q- I) * B + I
    # with B =
    # [
    #  I    .   .
    #  .  Dπv   .
    # .    .   I
    # ]

    # NOTE: double transpose because Dπv is BlockDiagonal
    # https://github.com/invenia/BlockDiagonals.jl/issues/16
    M = [
          spzeros(n,n)  (A' * Dπv)    c
         -A               -Dπv + I    b
         -c'           -(Dπv' * b)'   0
    ]
    # find projections on dual of the cones
    vp = π(MOI.dual_set.(cones), v)

    # dQ * [u, vp, max(0, w)]
    RHS = [dA' * vp + dc; -dA * u + db; -dc ⋅ u  - db ⋅ vp]

    dz = if norm(RHS) <= 1e-4
        RHS .= 0
    else
        lsqr(M, RHS)
    end

    du, dv, dw = dz[1:n], dz[n+1:n+m], dz[n+m+1]
    dx = du - x * dw
    dy = Dπv * dv - collect(Iterators.flatten(y)) * dw
    ds = Dπv * dv - dv - collect(Iterators.flatten(s)) * dw
    return -dx, -dy, -ds
end

MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{MOI.VariableIndex}) = true

function MOI.get(model::Optimizer, ::MOI.VariableName, v::VI)
    return MOI.get(model.optimizer, MOI.VariableName(), v)
end

function MOI.set(model::Optimizer, ::MOI.VariableName, v::VI, name::String)
    MOI.set(model.optimizer, MOI.VariableName(), v, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    return MOI.get(model.optimizer, MOI.VariableIndex, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, con::CI, name::String)
    MOI.set(model.optimizer, MOI.ConstraintName(), con, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    return MOI.get(model.optimizer, MOI.ConstraintIndex, name)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, con::CI)
    return MOI.get(model.optimizer, MOI.ConstraintName(), con)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{F,S})
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{VF, VS}}) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{VF,VS})
end

function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    MOI.set(model.optimizer, MOI.Name(), name)
end

function MOI.get(model::Optimizer, ::Type{CI{F, S}}, name::String) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, CI{F,S}, name)
end

function MOI.get(model::Optimizer, ::Type{CI{VF, VS}}, name::String) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, CI{VF,VS}, name)
end

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Real, Nothing})
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
end
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

MOI.supports(model::Optimizer, ::MOI.Silent) = MOI.supports(model.optimizer, MOI.Silent())

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    MOI.set(model.optimizer, MOI.Silent(), value)
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return MOI.get(model.optimizer, MOI.Silent())
end

function MOI.get(model::Optimizer, ::MOI.SolverName)
    return MOI.get(model.optimizer, MOI.SolverName())
end


text:
function check_same_dict(dict1::Dict, dict2::Dict, key; atol = 1e-3, ignore = ["source_version", "source_type", "name", "source_id"])
    if key in ignore
        return true
    end
    v1, v2 = dict1[key], dict2[key]
    if isa(v1, Dict)
        if isa(v2, Dict)
            return check_same_dict(v1,v2)
        else
            println("$v2 is not a dict")
            return false
        end
    elseif isa(v1, Real)
        if isa(v2, Real)
            if !isapprox(v1, v2, atol=atol)
                println("$v1 != $v2")
                return false
            end
        else
            println("$v2 is not a number")
            return false
        end
    else
        if v1 != v2
            println("$v1 != $v2")
            return false
        end
    end
    return true
end

function check_same_dict(dict1::Dict, dict2::Dict; atol = 1e-3, ignore = ["source_version", "source_type", "name", "source_id"])
    bools = Bool[]
    for (k,v1) in dict1
        if !(k in ignore)
            res = check_same_dict(dict1, dict2, k, ignore = ignore)
            push!(bools, res)
            if res
            else
            end
        end
    end
    return !(false in bools)
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSStopMaxIterations

using Test, RLinearAlgebra

@testset "LS Max Iteration Stop -- Procedural" begin

    # Verify appropriate super type
    @test supertype(LSStopMaxIterations) == LinSysStopCriterion

    # Verify check_stop_criterion functionality
    log = LSLogFull()
    stop = LSStopMaxIterations(10)

    log.iterations = 0
    @test RLinearAlgebra.check_stop_criterion(log, stop) == false

    log.iterations = 10
    @test RLinearAlgebra.check_stop_criterion(log, stop) == true

    log.iterations = 11
    @test RLinearAlgebra.check_stop_criterion(log, stop) == false
end

end # End module


text:
"""
$(TYPEDEF)
"""
struct AnalyticalProblem{uType,tType,isinplace,P,F,K} <: AbstractAnalyticalProblem{uType,tType,isinplace}
  f::F
  u0::uType
  tspan::tType
  p::P
  kwargs::K
  @add_kwonly function AnalyticalProblem{iip}(f,u0,tspan,p=NullParameters();
           kwargs...) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(_tspan),iip,typeof(p),
        typeof(f),typeof(kwargs)}(f,u0,_tspan,p,kwargs)
  end
end

function AnalyticalProblem(f,u0,tspan,p=NullParameters();kwargs...)
  iip = DiffEqBase.isinplace(f,4)
  AnalyticalProblem{iip}(f,u0,tspan,p;kwargs...)
end

export AnalyticalProblem, AbstractAnalyticalProblem


text:
reload("SimpleNeuralNets"); S=SimpleNeuralNets

X = vars'
Y = res'
# subset=idx; minibatch=100; steps=100; freq=0; noise=0.

nn = S.NN([19,10,10,1], [S.Relu(), S.Sigmoid()])
# copy!(nn.ws[1], ls[1])
# copy!(nn.ws[2], ls[2])
# copy!(nn.ws[3], ls[3])
λ0, μ = 1e-2, 0.99
S.sgd(nn, λ0, μ, X, Y, subset=idx, freq=200, steps=1000)

##############
reload("SimpleNeuralNets"); S=SimpleNeuralNets
nn = S.NN([19,10,10,1], [S.Relu(), S.Sigmoid()])

@time S.calc(nn, X) # 3.776
@time nn.ws[1] * X # 1.48 s

tst = convert(Array{Float64}, nn.ws[1])
@time tst * X # 0.005 s

@time S.calc2(nn, X) # 3.776
Profile.clear()
@profile S.calc(nn,X)
Profile.print()

Profile.print()
@profile S.sgd(nn, λ0, μ, X, Y, subset=idx, freq=200, steps=100)
Profile.print()

isa(X, StridedArray)
methods(S.cycle!, (S.TrainState, S.NN, SubArray{Float64,2,Array{Float64,2},Tuple{Colon,UnitRange{Int64}},2}, SubArray{Float64}))
methods(S.cycle!, (S.TrainState, S.NN, StridedArray{Float64}, StridedArray))

methods(S.cycle!, (SimpleNeuralNets.TrainState, SimpleNeuralNets.NN, SubArray{Float64,2,Array{Float64,2},Tuple{Colon,Array{Int64,1}},1}, SubArray))

typeof(sub(X,:,1:100))
cycle!(::SimpleNeuralNets.TrainState{3,100}, ::SimpleNeuralNets.NN{3}, ::SubArray{Float64,2,Array{Float64,2},Tuple{Colon,Array{Int64,1}},1}, ::SubArray{Float64,2,Array{Float64,2},Tuple{Colon,Array{Int64,1}},1})

T = sub(vars,120:320,:)
isa(T, SubArray)
isa(T, StridedArray)
methods(A_mul_B!, (DenseArray, SubArray, DenseArray))

@time S.sgd(nn, λ0, μ, X, Y, subset=idx, steps=10000) # 1.4 680Mb
@time S.sgd(nn, λ0, μ, X, Y, minibatch=10, subset=idx, steps=100000) # 10s 5400Mb
@time S.sgd(nn, λ0, μ, X, Y, minibatch=1000, subset=idx, steps=1000) # 0.52s 210Mb

# nn = S.NN((randn(10,19), randn(10,10), randn(1,10)), (S.Relu(), S.Relu(), S.Relu()))
# nn = S.NN((ls...,), (S.Relu(), S.Relu(), S.Relu()))

# sgd(nn, λ0, μ, X, Y, subset=idx, freq=200, steps=10000)

# @time sgd(nn, λ0, μ, X, Y, freq=200, steps=1000)
# @time sgd(nn, λ0, μ, X, Y, freq=0, steps=1000)

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
# authors: Carlo Baldassi and Carlo Lucibello
"""
    type NewtonMethod <: AbstractRootsMethod
        dx::Float64
        maxiters::Int
        verb::Int
        atol::Float64
    end

Type containg the parameters for Newton's root finding algorithm.
The default parameters are:

    NewtonMethod(dx=1e-7, maxiters=1000, verb=0, atol=1e-10)

"""
mutable struct NewtonMethod <: AbstractRootsMethod
    dx::Float64
    maxiters::Int
    verb::Int
    atol::Float64
end

NewtonMethod(; dx=1e-7, maxiters=1000, verb =0, atol=1e-10) =
                                    NewtonMethod(dx, maxiters, verb, atol)

function ∇!(∂f::Matrix, f::Function, x0, δ, f0, x1)
    n = length(x0)
    copy!(x1, x0)
    for i = 1:n
        x1[i] += δ
        ∂f[:,i] = (f(x1) - f0) / δ
        x1[i] = x0[i]
    end
    #=cf = copy(∂f)=#
    #=@time ∂f[:,:] = @parallel hcat for i = 1:n
        x1[i] += δ
        d = (f(x1) - f0) / δ
        x1[i] = x0[i]
        d
    end=#
    #@assert cf == ∂f
end

∇(f::Function, x0::Real, δ::Real, f0::Real) = (f(x0 + δ) - f0) / δ

"""
    newton(f::Function, x₀; pars=NewtonParameters())

Apply Newton's method with parameters `pars` to find a zero of `f` starting from the point
`x₀`.
The derivative of `f` is computed by numerical discretization. Multivariate
functions are supported.

Returns a tuple `(ok, x, it, normf)`.

**Usage Example**
ok, x, it, normf = newton(x->exp(x)-x^4, 1.)
ok || normf < 1e-10 || @warn("Newton Failed")

"""
function newton(f::Function, x₀::Float64, m::NewtonMethod)
    η = 1.0
    ∂f = 0.0
    x = x₀
    x1 = 0.0

    f0 = f(x)
    @assert isa(f0, Real)
    normf0 = abs(f0)
    it = 0
    while normf0 ≥ m.atol
        it > m.maxiters && return (false, x, it, normf0)
        it += 1
        if m.verb > 1
            println("(𝔫) it=$it")
            println("(𝔫)   x=$x")
            println("(𝔫)   f(x)=$f0")
            println("(𝔫)   normf=$(abs(f0))")
            println("(𝔫)   η=$η")
        end
        δ = m.dx
        while true
            try
                ∂f = ∇(f, x, δ, f0)
                break
            catch err
                @warn("newton: catched error:")
                Base.display_error(err, catch_backtrace())
                δ /= 2
                @warn("new δ = $δ")
            end
            if δ < 1e-15
                normf0 ≥ m.atol && @warn("newton:  δ=$δ")
                return (false, x, it, normf0)
            end
        end
        Δx = -f0 / ∂f
        m.verb > 1 && println("(𝔫)  Δx=$Δx")
        while true
            x1 = x + Δx * η
            local new_f0, new_normf0
            try
                new_f0 = f(x1)
                new_normf0 = abs(new_f0)
            catch err
                @warn("newton: catched error:")
                Base.display_error(err, catch_backtrace())
                new_normf0 = Inf
            end
            if new_normf0 < normf0
                η = min(1.0, η * 1.1)
                f0 = new_f0
                normf0 = new_normf0
                x = x1
                break
            end
            η /= 2
            η < 1e-15 && return (false, x, it, normf0)
        end
    end
    return true, x, it, normf0
end

function newton(f::Function, x₀, m::NewtonMethod)

    η = 1.0
    n = length(x₀)
    ∂f = Matrix{Float64}(n, n)
    x = Float64[x₀[i] for i = 1:n]
    x1 = Vector{Float64}(n)

    f0 = f(x)
    @assert length(f0) == n
    normf0 = vecnorm(f0)
    it = 0
    while normf0 ≥ m.atol
        it > m.maxiters && return (false, x, it, normf0)
        it += 1
        if m.verb > 1
            println("(𝔫) it=$it")
            println("(𝔫)   x=$x")
            println("(𝔫)   f0=$f0")
            println("(𝔫)   normf=$(vecnorm(f0))")
            println("(𝔫)   η=$η")
        end
        δ = m.dx
        while true
            try
                ∇!(∂f, f, x, δ, f0, x1)
                break
            catch
                δ /= 2
            end
            if δ < 1e-15
                normf0 ≥ m.atol && @warn("newton:  δ=$δ")
                return (false, x, it, normf0)
            end
        end
        if isa(f0, Vector)
            Δx = -∂f \ f0
        else
            Δx = -f0 / ∂f[1,1]
        end
        m.verb > 1 && println("(𝔫)  Δx=$Δx")
        while true
            for i = 1:n
                x1[i] = x[i] + Δx[i] * η
            end
            local new_f0, new_normf0
            try
                new_f0 = f(x1)
                new_normf0 = vecnorm(new_f0)
            catch
                new_normf0 = Inf
            end
            if new_normf0 < normf0
                η = min(1.0, η * 1.1)
                if isa(f0, Vector)
                    copy!(f0, new_f0)
                else
                    f0 = new_f0
                end
                normf0 = new_normf0
                copy!(x, x1)
                break
            end
            η /= 2
            η < 1e-15 && return (false, x, it, normf0)
        end
    end
    return true, x, it, normf0
end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
@testset "1648.sell-diminishing-valued-colored-balls.jl" begin
    @test max_profit([497978859, 167261111, 483575207, 591815159], 836556809) == 373219333
    @test max_profit([2, 5], 4) == 14
    @test max_profit([3, 5], 6) == 19
    @test max_profit([2, 8, 4, 10, 6], 20) == 110
    @test max_profit([1000000000], 1000000000) == 21
end


text:
"""
    Step{T, S<:Integer}(start, decay, step_sizes)
    Step(;λ, γ, step_sizes)

A step schedule decays exponentially by `γ` every step
in `step_sizes`.
The output conforms to
```text
λ * γ^{i - 1}
```
where `sum(step_sizes[1:(i - 1)]) < t <= sum(step_sizes[1:i])`

# Arguments:
- `start`/`λ`: the starting value
- `decay`/`γ`: the decay rate
- `step_sizes::Union{<:Integer, <:Vector}`: the step sizes
"""
struct Step{T, S<:Integer}
    start::T
    decay::T
    step_sizes::Vector{S}
end
Step(λ, γ, step_size::Integer) = Step(λ, γ, [step_size])
Step(;λ, γ, step_sizes) = Step(λ, γ, step_sizes)

function (schedule::Step)(t)
    i = findlast(x -> t > x, cumsum(schedule.step_sizes))
    i = isnothing(i) ? 0 :
            (i >= length(schedule.step_sizes)) ? length(schedule.step_sizes) - 1 : i

    return schedule.start * schedule.decay^i
end

Base.eltype(::Type{<:Step{T}}) where T = T
Base.IteratorSize(::Type{<:Step}) = Base.IsInfinite()

function Base.iterate(schedule::Step, state = (1, 1, 1))
    t, i, t0 = state
    if (i < length(schedule.step_sizes)) && (t >= t0 + schedule.step_sizes[i])
        # move onto next step range
        i += 1
        t0 = t
    end

    return schedule.start * schedule.decay^(i - 1), (t + 1, i, t0)
end

Base.axes(::Step) = (OneToInf(),)

"""
    Exp{T}(start, decay)
    Exp(;λ, γ)

A exponential decay schedule at rate `γ`.
The output conforms to
```text
λ * γ^{t - 1}
```

# Arguments:
- `start`/`λ`: the base value
- `decay`/`γ`: the decay rate
"""
struct Exp{T}
    start::T
    decay::T
end
Exp(;λ, γ) = Exp(λ, γ)

(schedule::Exp)(t) = schedule.start * schedule.decay^(t - 1)

Base.eltype(::Type{<:Exp{T}}) where T = T
Base.IteratorSize(::Type{<:Exp}) = Base.IsInfinite()

Base.iterate(schedule::Exp, t = 1) = schedule(t), t + 1

Base.axes(::Exp) = (OneToInf(),)

"""
    Poly{T, S<:Integer}(start, degree, max_iter)
    Poly(;λ, p, max_iter)

A polynomial schedule decays with degree `p`.
The output conforms to
```text
λ / (1 - (t - 1) / max_iter)^p
```

# Arguments
- `start`/`λ`: the base value
- `degree`/`p::Integer`: the degree of the polynomial
- `max_iter::Integer`: the total number of iterations
"""
struct Poly{T, S<:Integer}
    start::T
    degree::S
    max_iter::S
end
Poly(;λ, p, max_iter) = Poly(λ, p, max_iter)

function (schedule::Poly)(t)
    (t <= length(schedule)) || throw(BoundsError("Cannot index Poly for t > max_iter"))
    return schedule.start * (1 - (t - 1) / schedule.max_iter)^schedule.degree
end

Base.eltype(::Type{<:Poly{T}}) where T = T
Base.IteratorSize(::Type{<:Poly}) = Base.HasLength()
Base.length(schedule::Poly) = schedule.max_iter

Base.iterate(schedule::Poly, t = 1) = schedule(t), t + 1

Base.axes(schedule::Poly) = 1:length(schedule)

"""
    Inv{T, S<:Integer}(start, decay, degree)
    Inv(;λ, γ, p)

A decay schedule that inversely decays with rate `γ`.
The output conforms to
```text
λ / (1 + (t - 1) * γ)^p
```

# Arguments
- `start`/`λ`: the base value
- `decay`/`γ`: the decay rate
- `degree`/`p::Integer`: the degree of decay
"""
struct Inv{T<:Number, S<:Integer}
    start::T
    decay::T
    degree::S
end
Inv(;λ, γ, p) = Inv(λ, γ, p)

(schedule::Inv)(t) = schedule.start / (1 + (t - 1) * schedule.decay)^schedule.degree

Base.eltype(::Type{<:Inv{T}}) where T = T
Base.IteratorSize(::Type{<:Inv}) = Base.IsInfinite()

Base.iterate(schedule::Inv, t = 1) = schedule(t), t + 1

Base.axes(::Inv) = (OneToInf(),)


text:
using Homebrew

# First time there is an error. Works the 2nd time though. Not sure why, but
# seems to be a problem with GnuTLS.
try
  using GnuTLS
catch
end

using GnuTLS

using Requests

r = get("http://rubygems.org/profiles/feldt")

m = match(r"
([0-9,]+) all time
", r.data)

if m == nothing

  println("ERROR! Could not identify a dl count!")
  exit(-1)

else

  count = replace(m.captures[1], ",", "")
  num_dls = int(count)
  rounded = int( floor(num_dls/1000) * 1000 )
  datestr = strftime("%B %e %Y", time())
  print("Feldt's Ruby gems has been downloaded ($(datestr)) more than $(rounded) times ($(num_dls))")

end


text:
using OceanParameterizations
using WindMixing
using OceanTurb
using BenchmarkTools
using OrdinaryDiffEq

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis",
            "weak_wind_strong_cooling", "strong_wind_weak_cooling", "strong_wind_weak_heating"]

test_file = "free_convection"

𝒟test = WindMixing.data(test_file)
𝒟train = 𝒟test # just for Benchmarking purposes
les = read_les_output(test_file)

# KPP Parameterization (no training)
#
# parameters = KPP.Parameters() # default parameters
# KPP_model = closure_kpp_full_evolution(parameters, 𝒟test.T.coarse[:,1], les)
#
# # Time building the model
# @btime closure_kpp_full_evolution(parameters, 𝒟test.T.coarse[:,1], les) # 86.448 μs (168 allocations: 19.03 KiB)
#
# # Time running the model
# @btime KPP_model() # 25.762 ms (4774 allocations: 242.36 KiB)
#
# ## TKE Parameterization (no training; use default parameters)
#
# # les = read_les_output(test_file)
# parameters = TKEMassFlux.TKEParameters() # default parameters
# TKE_model = closure_tke_full_evolution(parameters, 𝒟test.T.coarse[:,1], les)
#
# # Time building the model
# @btime closure_tke_full_evolution(parameters, 𝒟test.T.coarse[:,1], les) # 92.100 μs (213 allocations: 23.08 KiB)
#
# # Time running the model
# @btime TKE_model() # 2.710 ms (5894 allocations: 259.86 KiB)

## GPR
function build()
    # Set the kernel manually
    uw_kernel = get_kernel(1,0.1,0.0,euclidean_distance)
    vw_kernel = get_kernel(1,0.1,0.0,euclidean_distance)
    wT_kernel = get_kernel(1,0.1,0.0,euclidean_distance)

    # Trained GP models
    uw_GP_model = gp_model(𝒟train.uw, uw_kernel)
    vw_GP_model = gp_model(𝒟train.vw, vw_kernel)
    wT_GP_model = gp_model(𝒟train.wT, wT_kernel)

    return (uw_GP_model, vw_GP_model, wT_GP_model)
end

println("Time to build 3 GP models:")
@btime build() # 42.818 ms (148 allocations: 10.07 MiB)

uw_GP_model, vw_GP_model, wT_GP_model = build()
uw_GP = predict(𝒟test.uw, uw_GP_model)
vw_GP = predict(𝒟test.vw, vw_GP_model)
wT_GP = predict(𝒟test.wT, wT_GP_model)

@time predict(𝒟test.wT, wT_GP_model)
println("Time to predict wT:")
println(@time predict(𝒟test.wT, wT_GP_model))

function run_model(uw_GP_model, vw_GP_model, wT_GP_model)
    # GP predictions on test data
    # uw_GP = predict(𝒟test.uw, uw_GP_model)
    # vw_GP = predict(𝒟test.vw, vw_GP_model)
    # wT_GP = predict(𝒟test.wT, wT_GP_model)
    t = 𝒟test.t
    uvT₀      = 𝒟test.uvT_scaled[:,1]
    zF_coarse = 𝒟test.uw.z
    zC_coarse = 𝒟test.u.z
    f⁰        = les.f⁰
    Nz        = 32

    ∂z(vec) = (vec[1:Nz] .- vec[2:Nz+1]) ./ diff(zF_coarse)
    function f(dx, x, p, t)
        u = x[1:Nz]
        v = x[Nz+1:2*Nz]
        dx[1:Nz] .= -∂z(uw_GP_model(x)) .+ f⁰ .* v
        dx[Nz+1:2*Nz] .= -∂z(vw_GP_model(x)) .- f⁰ .* u
        dx[2*Nz+1:end] .= -∂z(wT_GP_model(x))
    end

    prob = ODEProblem(f, uvT₀, (t[1],t[289]), saveat=t)
    sol = solve(prob, ROCK4())

    # split_array(uvT) = (uvT[1:Nz,:], uvT[Nz+1:2*Nz,:], uvT[2*Nz+1:end,:])
    # u_pred, v_pred, T_pred = split_array(sol)
    return sol
end

println("Time to run 3 GP models + solve:")
println(@time run_model(uw_GP_model, vw_GP_model, wT_GP_model))


text:
using Oceananigans
using Oceananigans.Architectures: arch_array
using Oceananigans: AbstractModel
using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel

function get_parameter(filename, group, parameter_name, default=nothing)
    parameter = default

    jldopen(filename) do file
        if parameter_name ∈ keys(file["$group"])
            parameter = file["$group/$parameter_name"]
        end
    end

    return parameter
end

"""
    OneDimensionalEnsembleModel(observations::OneDimensionalTimeSeriesBatch; architecture = CPU(), ensemble_size = 50, kwargs...)

Build an Oceananigans `HydrostaticFreeSurfaceModel` with many independent
columns. The model grid is given by the data in `observations`, and the
dynamics in each column is attached to its own `CATKEVerticalDiffusivity` closure stored in
an `(Nx, Ny)` Matrix of closures. `ensemble_size = 50` is the
desired count of ensemble members for calibration with Ensemble Kalman Inversion (EKI), and the
remaining keyword arguments `kwargs` define the default closure across all columns.

In the "many columns" configuration, we run the model on a 3D grid with `(Flat, Flat, Bounded)` boundary
conditions so that many independent columns can be evolved at once with much of the computational overhead
split among the columns. The `Nx` rows of vertical columns are each reserved for an "ensemble" member
whose attached parameter value (updated at each iteration of EKI) sets the diffusivity closure
used to predict the model solution for the `Ny` physical scenarios described by the simulation-specific
`OneDimensionalTimeSeries` objects in `observations`.
"""
function OneDimensionalEnsembleModel(observations::OneDimensionalTimeSeriesBatch;
                       architecture = CPU(),
                       ensemble_size = 50,
                       closure = CATKEVerticalDiffusivity(Float64; warning=false, kwargs...)
                      )

    #
    # Harvest observation file metadata
    #

    datapath = observations.file_path

    α = get_parameter(datapath, "buoyancy", "equation_of_state/α", 2e-4)
    g = get_parameter(datapath, "buoyancy", "gravitational_acceleration", 9.81)
    αg = α * g

    f = get_parameter(datapath, "coriolis", "f", 0.0)

    # Surface fluxes
    Qᵘ = get_parameter(datapath, "parameters", "boundary_condition_u_top", 0.0)
    Qᵛ = get_parameter(datapath, "parameters", "boundary_conditions_v_top", 0.0)
    Qᶿ = get_parameter(datapath, "parameters", "boundary_condition_θ_top", 0.0)
    Qᵇ = Qᶿ * αg

    # Bottom gradients
    dudz_bottom = get_parameter(datapath, "parameters", "boundary_condition_u_bottom", 0.0)
    dθdz_bottom = get_parameter(datapath, "parameters", "boundary_condition_θ_bottom", 0.0)
    dbdz_bottom = dθdz_bottom * αg

    #
    # Build model using metadata
    #

    grid = OneDimensionalEnsembleGrid(observations.grid; size=(ensemble_size, length(observations), data_grid.Nz))

    closure = [closure for i=1:ensemble_size, j=1:length(observations)]

    coriolis = [FPlane(f=td.constants[:f]) for i=1:ensemble_size, td in observations]
    coriolis = arch_array(architecture, coriolis)

    bc_matrix(x) = [x for i = 1:ensemble_size, d in observations]
    Qᵇ = bc_matrix(Qᵇ)
    Qᵘ = bc_matrix(Qᵘ)
    Qᵛ = bc_matrix(Qᵛ)
    dbdz_bottom = bc_matrix(dbdz_bottom)
    dudz_bottom = bc_matrix(dudz_bottom)

    # Convert to CuArray if necessary
    closure = arch_array(architecture, closure)
    Qᵇ = arch_array(architecture, Qᵇ)
    Qᵘ = arch_array(architecture, Qᵘ)
    Qᵛ = arch_array(architecture, Qᵛ)
    dbdz_bottom = arch_array(architecture, dbdz_bottom)
    dudz_bottom = arch_array(architecture, dudz_bottom)

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ),
                                    bottom = GradientBoundaryCondition(dudz_bottom))
    v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),
                                    bottom = GradientBoundaryCondition(dbdz_bottom))

    model = HydrostaticFreeSurfaceModel(architecture = architecture,
                                         grid = grid,
                                         tracers = (:b, :e),
                                         buoyancy = BuoyancyTracer(),
                                         coriolis = coriolis,
                                         boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                         closure = closure)

    return model
end


text:
export Client, open!, wait, exit_session_on_err, request, set_event_handler!

mutable struct Client
    id::String
    token::String # Authorization
    intents::Int # https://discord.com/developers/docs/topics/gateway#list-of-intents
    conn::VersionedWebSocket
    geteway_url::String
    request_func::Function
    ready::Bool # True after connectiong to geteway
    hb::Heartbeat
    session_id::String # got from Ready
    seq_num::Int # sequence number of the last event received
    event_h::EventHandler

    function Client(id::String, token::String, intents::Int)
        token_field = "Bot $token"
        request_func = create_request_func("$URL/v$API_VERSION/", Dict(
                "Content-Type" => "application/json",
                "Authorization" => token_field,
            )
        )
        event_h = EventHandler()
        set!(event_h, "READY", (client, data) -> client.session_id = data[:session_id])
        new(
            id,
            token_field,
            intents,
            VersionedWebSocket(0, nothing),
            "",
            request_func,
            false,
            Heartbeat(),
            "",
            0,
            event_h
         )
    end
end

function request(c::Client, method::Symbol, dest::String; payload="", json::Bool=true)
    return try
        !isempty(payload) && json && (payload = JSON.json(payload))
        c.request_func(method, dest; body=payload), nothing
    catch e
        nothing, e
    end
end

Base.isopen(c::Client) = c.ready && isopen(c.conn)
set_event_handler!(c::Client, event_name::String, handler::Function) = set!(c.event_h, event_name, handler)

function open!(c::Client; resume::Bool=false, reconn::Bool=false)
    c.ready || ErrorException("connection with the gateway has already been established")
    if isempty(c.geteway_url)
        c.geteway_url, err = get_gateway_url()
        isnothing(err) || return err
        @info "Received the geteway url."
    end

    if reconn || !isopen(c.conn)
        @info "Establishing a WebSocket connection..." geteway = c.geteway_url

        err = open!(c.conn, c.geteway_url)
        isnothing(err) || return err
        reconn && version_update!(c.conn)

        data, err = read_json(c.conn.io)
        isnothing(err) || return err
        GOPCODES(data[:op]) == HELLO || return ErrorException("on getting unexpected opcode: $(data[:op])")
        @info "HELLO opcode was received."
        set_interval!(c.hb, data[:d][:heartbeat_interval])
    end

    err = if resume && !isempty(c.session_id)
        @info "Attempting to resume the lost session..."
        err = resume_conn(c)
        isnothing(err) && @info begin
            session_id = c.session_id
            seq_num = c.seq_num
            "Connection was resumed successfully."
        end
        err
    else
        @info "Attempting to authorize with a token..."
        hello(c)
    end
    isnothing(err) || return err
    c.ready = true
    @async event_listener!(c)
    @async heartbeat_loop!(c)
    nothing
end

function close!(c::Client)
    c.ready = false
    isopen(c.conn) && close!(c.conn)
end

const RECONNECT_WAIT_TIME = 30
function Base.wait(c::Client)
    @info "Started waiting on socket connection..."
    while isopen(c)
        wait(c.conn)
        sleep(RECONNECT_WAIT_TIME)
    end
    @info "Stopped waiting." client = c
end

function exit_session_on_err(c::Client, err)
    isnothing(err) && return
    close!(c)
    @error err client = c
    # exit(1)
end

function event_listener!(c::Client)
    socket_v = c.conn.v
    @info "Started listening for new events..."
    try
        while isopen(c) && socket_v == c.conn.v
            data, err = read_json(c.conn.io)
            throw_on_err(err)
            err = handle(data[:op])(c, data)
            isnothing(err) || @error err
        end
    catch e
        @error "Error has occured while listening for events." err = e
        @info "Trying to reconnect..."
        exit_session_on_err(c, reconnect!(c))
        @info "Reconnected successfully." socket_version = c.conn.v
    end
    @info "Stopped listening for events..." socket_version = socket_v
end

function heartbeat_loop!(c::Client)
    socket_v = c.conn.v
    @info "Started sending the heartbeat signals every $(c.hb.interval)s..." socket_version = socket_v
    sleep(rand(1:c.hb.interval))
    try
        while isopen(c) && socket_v == c.conn.v
            # Handle a zombie connection
            if is_zombie(c.hb)
                # TODO: RECONNECT FALLBACK
                throw(ErrorException("on missing a heartbeat-ack message from the geteway"))
            end
            throw_on_err(heartbeat(c))
            update!(c.hb)
            sleep(c.hb.interval)
        end
    catch e
        @error "Error has occured while sending a heartbeat." err = e
    end
end


text:
using BoxChopper
using Base.Test

const Float = typeof(0.0)

@testset "BoxChopper" begin
    @testset "triangular cylinder" begin
        box = ((0,1), (0,1), (0,1))
        nout = (1,1,0)

        @test (r₀ = (0.5,0,0); volfrac(box, nout, r₀) ≈ 0.125)
        @test (r₀ = (0.5,0,0); volfrac(box, (-nout[1],-nout[2],-nout[3]), r₀) ≈ 0.875)
        @test (r₀ = (1,0,0); volfrac(box, nout, r₀) ≈ 0.5)
        @test (r₀ = (1,0,0); volfrac(box, (-nout[1],-nout[2],-nout[3]), r₀) ≈ 0.5)
        @test (r₀ = (1,0,0); nout = (1,2,0); volfrac(box, nout, r₀) ≈ 0.25)
        @test (r₀ = (1,0,0); nout = (1,2,0); volfrac(box, (-nout[1],-nout[2],-nout[3]), r₀) ≈ 0.75)
    end

    @testset "quadrangular cylinder" begin
        @test begin
            result = true
            for i = 1:100
                box = ((-rand(),rand()), (-rand(),rand()), (-rand(),rand()))
                r₀ = ((box[1][1]+box[1][2])/2, (box[2][1]+box[2][2])/2, (box[3][1]+box[3][2])/2)
                nout = randn(3)
                nout[rand(1:3)] = 0
                nout = (nout...)
                result &= volfrac(box, nout, r₀) ≈ 0.5
            end
            result
        end
    end

    @testset "general cases" begin
        # Test random cases.
        @test begin
            result = true
            for i = 1:100
                box = ((-rand(),rand()), (-rand(),rand()), (-rand(),rand()))
                r₀ = ((box[1][1]+box[1][2])/2, (box[2][1]+box[2][2])/2, (box[3][1]+box[3][2])/2)
                nout = (randn(), randn(), randn())
                result &= volfrac(box, nout, r₀) ≈ 0.5
            end
            result
        end

        @test begin
            result = true
            for i = 1:100
                box = ((-rand(),rand()), (-rand(),rand()), (-rand(),rand()))
                r₀ = (randn(), randn(), randn())
                nout = (randn(), randn(), randn())
                result &= (volfrac(box, nout, r₀) + volfrac(box, (-nout[1],-nout[2],-nout[3]), r₀) ≈ 1.)
            end
            result
        end

        # Test boundary cases.
        box = ((0,1), (0,1), (0,1))
        r₀ = (0,0,0)
        @test (nout = (1,0,0); volfrac(box, nout, r₀) ≈ 0.)  # completely outside
        @test (nout = (-1,0,0); volfrac(box, nout, r₀) ≈ 1.)  # completely inside
        @test (nout = (-1,1,0); volfrac(box, nout, r₀) ≈ 0.5)  # rvol_tricyl()
        @test (nout = (-1,-1,1); volfrac(box, nout, r₀) ≈ 5/6)  # rvol_gensect()
        @test (nout = (1,-2,1); volfrac(box, nout, r₀) ≈ 0.5)  # rvol_quadsect()
        r₀ = (0.5,0.5,0)
        @test (nout = (-2,1,0); volfrac(box, nout, r₀) ≈ 0.5)  # rvol_quadcyl()

        # Test rvol_quadsect() for nontrivial cases.
        box = ((0,1), (0,1), (0,2))
        r₀ = (0.5, 0.5, 0.5)
        @test begin
            result = true
            for i = 1:100
                nout = (randn()/20, randn()/20, 1)
                result &= volfrac(box, nout, r₀) ≈ 0.5/2
            end
            result
        end
    end

end  # @testset "BoxChopper"