text:
#####################################
#     Basic interface functions     #
#####################################

"""
    cost,i1,i2 = dtw(seq1, seq2, [dist=SqEuclidean, postprocess=nothing])
    cost,i1,i2 = dtw(seq1, seq2, dist, i2min, i2max)

Perform dynamic-time warping to measure the distance between two sequences.

Find a set of indices (`i1`,`i2`) that align two series (`seq1`,`seq2`) by
dynamic axis warping. Also returns the distance (after warping) according to
the SemiMetric `dist`, which defaults to squared Euclidean distance (see
Distances.jl). If `seq1` and `seq2` are matrices, each column is considered
an observation.

If `i2min/max` are provided, do DTW to align `seq1` and `seq2` confined to a window. Vectors `i2min` and
`i2max` specify (inclusive) lower and upper bounds for `seq2` for each index in
`seq1`. Thus, `i2min` and `i2max` are required to be the same length as `seq1`.

If `filternernel::AbstractMatrix` is provided, it's used to filter the cost matrix. Create a suitable kerlen using, e.g., `ImageFiltering.Kernel.gaussian(3)`. The filtering of the cost matrix makes the warping smoother, effectively penalizing small-scale warping.

See also [`dtw_cost`](@ref) and [`dtwnn`](@ref).
"""
function dtw(args...; kwargs...)
    D = dtw_cost_matrix(args...; kwargs...)
    return trackback(D)
end

##############################
#  Cost matrix computations  #
##############################

Distances.pairwise(d::PreMetric, s1::AbstractVector, s2::AbstractVector; dims=2) = evaluate.(Ref(d), s1, transpose(s2))

function Distances.pairwise(d::PreMetric, s1::AbstractArray, s2::AbstractArray; dims=2)
    [evaluate(d, s1[!,i], s2[!,j]) for i in 1:lastlength(s1), j in lastlength(s2)]
end

@inbounds function dtw_cost_matrix(seq1::AbstractArray{T}, seq2::AbstractArray{T}, dist::SemiMetric = SqEuclidean();
    transportcost=1,
    postprocess = nothing) where T
    # Build the cost matrix
    m = lastlength(seq2)
    n = lastlength(seq1)

    # Initialize first column and first row
    D = pairwise(dist, seq2, seq1, dims=2)
    @assert size(D) == (m,n)

    for r=2:m
        D[r,1] += D[r-1,1]
    end
    for c=2:n
        D[1,c] += D[1,c-1]
    end

    # Complete the cost matrix
    for c = 2:n
        for r = 2:m
            best_neighbor_cost = min(transportcost*D[r-1, c], D[r-1, c-1], transportcost*D[r, c-1])
            D[r, c] += best_neighbor_cost
        end
    end

    if postprocess !== nothing
        D = postprocess(D)
    end

    return D
end

Base.@propagate_inbounds function dtw_cost_matrix(
    seq1::AbstractArray{T},
    seq2::AbstractArray{T},
    dist::SemiMetric,
    i2min::AbstractVector{U},
    i2max::AbstractVector{U};
    transportcost = 1
) where {T,U<:Integer}
    n = lastlength(seq1) # of columns in cost matrix
    m = lastlength(seq2) # of rows in cost matrix
    Base.@boundscheck begin
        n == length(i2min) || throw(ArgumentError("i2min does not match length of seq1."))
        n == length(i2max) || throw(ArgumentError("i2max does not match length of seq1."))
        1 == i2min[1]      || throw(ArgumentError("i2min must start at 1."))
        m == i2max[end]    || throw(ArgumentError("i2max must end at length(seq2), was $(i2max[end]) ≂̸ $(m)"))
    end

    # Build the (n x m) cost matrix into a WindowedMatrix, because it's ragged.
    # That type gives efficient storage with convenient [r,c] indexing and returns
    # Inf when accessed outside the window.
    D = WindowedMatrix(i2min, i2max, Inf)

    # First column first
    D[1, 1] = evaluate(dist, seq1[!, 1], seq2[!, 1])
    for r = 2:i2max[1]
        D[r, 1] = D[r-1, 1] + evaluate(dist, seq1[!, 1], seq2[!, r])
    end

    # Complete the cost matrix from columns 2 to m.
    for c = 2:n
        for r = i2min[c]:i2max[c]
            best_neighbor_cost = min(transportcost*D[r-1, c], D[r-1, c-1], transportcost*D[r, c-1])
            D[r, c] = best_neighbor_cost + evaluate(dist, seq1[!, c], seq2[!, r])
        end
    end

    return D
end

########################################
#  Find Best Path through Cost Matrix  #
########################################

"""
    cost,cols,rows = trackback(D::Matrix)

Given the cost matrix `D`, computes the optimal track from end to beginning.
Returns `cols` and `rows` which are vectors respectively holding the track.
"""
function trackback(D::AbstractMatrix{T}) where {T<:Number}

    # initialize trackback throught rows/columns
    r, c       = size(D)
    rows, cols = Int[r], Int[c]

    # estimate that we'll need N⋅logN elements
    N  = max(r, c)
    sz = 2 * N
    sizehint!(rows, sz)
    sizehint!(cols, sz)

    # do trackback
    @inbounds while r > 1 && c > 1
        tb, r, c = indmin3(D[r-1, c-1], D[r-1, c], D[r, c-1], r, c)
        push!(rows, r)
        push!(cols, c)
    end
    # Possibly either r>1 or c>1 at this point (but not both).
    # Add the unfinished part of the track to reach [1,1]
    for r = r-1:-1:1
        push!(rows, r)
        push!(cols, 1)
    end
    for c = c-1:-1:1
        push!(rows, 1)
        push!(cols, c)
    end
    return D[end, end], reverse(cols), reverse(rows)
end

"""
    dtw_cost(a::AbstractArray, b::AbstractArray, dist::Distances.SemiMetric, r::Int; best_so_far = Inf, cumulative_bound = Zeros(length(a)))

Perform dynamic time warping to measure the distance between two sequences.

Calculate the DTW cost between `a` and `b` with maximum warping radius `r`. You may provide values of `best_so_far` and `cumulative_bound` in order to enable early stopping.

# Keyword arguments:
- `best_so_far`: The best cost value obtained so far (optional)
- `cumulative_bound`: A vector the same length as a and b (optional)
- `s1`: Optional storage vector of length 2r+1, can be used to save allocations.
- `s2`: Optional storage vector of length 2r+1, can be used to save allocations.

Providing the two vectors `s1, s2` does not save very much time, but it makes the function completely allocation free. Can be useful in a threaded context.

See also [`dtw`](@ref) and [`dtwnn`](@ref).
This code was inspired by https://www.cs.ucr.edu/~eamonn/UCRsuite.html
"""
function dtw_cost(
    a::AbstractArray{QT},
    b::AbstractArray,
    dist::F,
    r::Int;
    transportcost = 1,
    best_so_far = typemax(floattype(QT)),
    cumulative_bound = Zeros(lastlength(a)),
    s1 = fill(typemax(floattype(QT)), 2r + 1),
    s2 = fill(typemax(floattype(QT)), 2r + 1),
    kwargs...
) where {QT, F <: Union{Distances.SemiMetric, Function}}

    T = floattype(QT)
    cost      = s1 # just change the name of these variables
    cost_prev = s2

    # Instead of using matrix of size O(m^2) or O(mr), we will reuse two array of size O(r).
    m, mb = lastlength(a), lastlength(b)
    mb == m || throw(ArgumentError("a and b must have the same length, got $m and $mb. To compare two series of different lengths, use function dtw"))
    length(cumulative_bound) == m || throw(ArgumentError("cumulative_bound and a must have the same length."))
    length(s1) == 2r+1 || throw(ArgumentError("s1 must be length 2r+1."))
    length(s2) == 2r+1 || throw(ArgumentError("s2 must be length 2r+1."))

    local k

    for i = 0:m-1
        k = max(0, r - i)
        min_cost = typemax(T)

        for j = max(0, i - r):min(m - 1, i + r)
            if j == 0 && i == 0
                cost[k+1] = dist(a[!,1], b[!,1]; kwargs...)
                min_cost = cost[k+1]
                k += 1
                continue
            end
            y = (j - 1 < 0) || (k - 1 < 0)     ? typemax(T) : cost[k]
            x = (i - 1 < 0) || (k + 1 > 2 * r) ? typemax(T) : transportcost*cost_prev[k+2]
            z = (i - 1 < 0) || (j - 1 < 0)     ? typemax(T) : transportcost*cost_prev[k+1]

            cost[k+1] = min(x, y, z) + dist(a[!,i+1], b[!,j+1]; kwargs...)

            # Find minimum cost in row for early stopping
            if cost[k+1] < min_cost
                min_cost = cost[k+1]
            end
            k += 1
        end

        # We can abandon early if the current cumulative distace with lower bound together are larger than best_so_far
        if ((i + r) < (m - 1)) && (min_cost + cumulative_bound[i+r+1] >= best_so_far)
            return min_cost + cumulative_bound[i+r+1]
        end

        cost_prev, cost = cost, cost_prev
    end

    # the DTW distance is in the last cell in the matrix of size O(m^2) or at the middle of our array.
    final_dtw = cost_prev[k]
    return T(final_dtw)
end

@inbounds function soft_dtw_cost_matrix(seq1::AbstractArray, seq2::AbstractArray, dist::SemiMetric = SqEuclidean(); γ = 1,
    transportcost=1, radius=nothing)
    # Build the cost matrix
    m = lastlength(seq2)
    n = lastlength(seq1)

    # Initialize first column and first row
    D = pairwise(dist, seq2, seq1, dims=2)
    @assert size(D) == (m,n)

    for r=2:m
        D[r,1] += D[r-1,1]
    end
    for c=2:n
        D[1,c] += D[1,c-1]
    end

    # Complete the cost matrix
    if radius === nothing
        for c = 2:n
            for r = 2:m
                D[r, c] += softmin(transportcost*D[r-1, c], D[r-1, c-1], transportcost*D[r, c-1], γ)
            end
        end
    else
        for c = 2:n
            for r = 2:m
                if abs(c-r) > radius
                    D[r, c] += 1/γ
                    # continue
                end
                D[r, c] += softmin(transportcost*D[r-1, c], D[r-1, c-1], transportcost*D[r, c-1], γ)
            end
        end
    end

    return D
end

"""
    soft_dtw_cost(args...; γ = 1, kwargs...)

Perform Soft DTW. This is a differentiable version of DTW. The "distance" returned by this function is quite far from a true distance and can be negative. A smaller value of `γ` makes the distance closer to the standard DTW distance.

To differentiate w.r.t. the first argument, try
```julia
using ReverseDiff
da = ReverseDiff.gradient(a->soft_dtw_cost(a,b), a)
```

Ref: "Soft-DTW: a Differentiable Loss Function for Time-Series" https://arxiv.org/pdf/1703.01541.pdf

#Arguments:
- `args`: same as for [`dtw`](@ref)
- `γ`: The smoothing factor. A small value means less smoothing and a result closer to [`dtw_cost`](@ref)
- `kwargs`: same as for [`dtw`](@ref)
"""
function soft_dtw_cost(args...; γ = 1, kwargs...)
    D = soft_dtw_cost_matrix(args...; γ = γ, kwargs...)
    D[end,end]
end

@fastmath @inline function softmin(a, b, c, γ)
    γ = -γ
    a,b,c = a/γ, b/γ, c/γ
    maxv = max(a,b,c)
    ae,be,ce = exp(a - maxv), exp(b - maxv), exp(c - maxv)
    γ*(log(ae+be+ce) + maxv)
end

const LVB = LoopVectorization.VectorizationBase
@inline function softmin(a::T, b::T, c::T, γ) where T <: Union{Float64, Float32}
    γ = -γ
    ninvγ = one(T) / γ
    v = LVB.Vec{4,T}(a, b, c, typemax(T))
    v = v * ninvγ
    @fastmath maxv = min(a,b,c) * ninvγ
    ve = exp(v - maxv) * LVB.Vec{4,T}(one(T), one(T), one(T), zero(T))
    γ*(log(LVB.vsum(ve)) + maxv)
end


text:
using EscapeAnalysis, Test
import Core: Argument, SSAValue, ReturnNode

isT(T) = (@nospecialize x) -> x === T
issubT(T) = (@nospecialize x) -> x <: T
isreturn(@nospecialize x) = isa(x, Core.ReturnNode) && isdefined(x, :val)
isthrow(@nospecialize x) = Meta.isexpr(x, :call) && Core.Compiler.is_throw_call(x)
isnew(@nospecialize x) = Meta.isexpr(x, :new)
isϕ(@nospecialize x) = isa(x, Core.PhiNode)
function with_normalized_name(@nospecialize(f), @nospecialize(x))
    if Meta.isexpr(x, :foreigncall)
        name = x.args[1]
        nn = EscapeAnalysis.normalize(name)
        return isa(nn, Symbol) && f(nn)
    end
    return false
end
isarrayalloc(@nospecialize x) = with_normalized_name(nn->!isnothing(Core.Compiler.alloc_array_ndims(nn)), x)
isarrayresize(@nospecialize x) = with_normalized_name(nn->!isnothing(EscapeAnalysis.array_resize_info(nn)), x)
isarraycopy(@nospecialize x) = with_normalized_name(nn->EscapeAnalysis.is_array_copy(nn), x)
import Core.Compiler: argextype, singleton_type
iscall(y) = @nospecialize(x) -> iscall(y, x)
function iscall((ir, f), @nospecialize(x))
    return iscall(x) do @nospecialize x
        singleton_type(Core.Compiler.argextype(x, ir, Any[])) === f
    end
end
iscall(pred::Function, @nospecialize(x)) = Meta.isexpr(x, :call) && pred(x.args[1])

# check if `x` is a statically-resolved call of a function whose name is `sym`
isinvoke(y) = @nospecialize(x) -> isinvoke(y, x)
isinvoke(sym::Symbol, @nospecialize(x)) = isinvoke(mi->mi.def.name===sym, x)
isinvoke(pred::Function, @nospecialize(x)) = Meta.isexpr(x, :invoke) && pred(x.args[1]::Core.MethodInstance)

"""
    is_load_forwardable(x::EscapeInfo) -> Bool

Queries if `x` is elibigle for store-to-load forwarding optimization.
"""
function is_load_forwardable(x::EscapeAnalysis.EscapeInfo)
    AliasInfo = x.AliasInfo
    AliasInfo === false && return true # allows this query to work for immutables since we don't impose escape on them
    # NOTE technically we also need to check `!has_thrown_escape(x)` here as well,
    # but we can also do equivalent check during forwarding
    return isa(AliasInfo, EscapeAnalysis.IndexableFields) || isa(AliasInfo, EscapeAnalysis.IndexableElements)
end

let setup_ex = quote
        mutable struct SafeRef{T}
            x::T
        end
        Base.getindex(s::SafeRef) = getfield(s, 1)
        Base.setindex!(s::SafeRef, x) = setfield!(s, 1, x)

        mutable struct SafeRefs{S,T}
            x1::S
            x2::T
        end
        Base.getindex(s::SafeRefs, idx::Int) = getfield(s, idx)
        Base.setindex!(s::SafeRefs, x, idx::Int) = setfield!(s, idx, x)
    end
    global function EATModule(setup_ex = setup_ex)
        M = Module()
        Core.eval(M, setup_ex)
        return M
    end
    Core.eval(@__MODULE__, setup_ex)
end

ext: jl


text:
"""
Sets how to connect to the Bloomberg API.

* `BLPAPI_CLIENTMODE_AUTO` tries to
connect to Desktop API, and falls back to Server API.

* `BLPAPI_CLIENTMODE_DAPI` connects to Desktop API.

* `BLPAPI_CLIENTMODE_SAPI` connects to Server API.

The default when creating `SessionOptions`
is `BLPAPI_CLIENTMODE_AUTO`.

See also [`Session`](@ref).
"""
@enum ClientMode::Cint begin
    BLPAPI_CLIENTMODE_AUTO        =  0
    BLPAPI_CLIENTMODE_DAPI        =  1
    BLPAPI_CLIENTMODE_SAPI        =  2
    BLPAPI_CLIENTMODE_COMPAT_33X  =  16
end

# as defined in blpapi_types.h
@enum BLPDataType::Cint begin
    BLPAPI_DATATYPE_BOOL           = 1  # Bool
    BLPAPI_DATATYPE_CHAR           = 2  # Char
    BLPAPI_DATATYPE_BYTE           = 3  # Unsigned 8 bit value
    BLPAPI_DATATYPE_INT32          = 4  # 32 bit Integer
    BLPAPI_DATATYPE_INT64          = 5  # 64 bit Integer
    BLPAPI_DATATYPE_FLOAT32        = 6  # 32 bit Floating point - IEEE
    BLPAPI_DATATYPE_FLOAT64        = 7  # 64 bit Floating point - IEEE
    BLPAPI_DATATYPE_STRING         = 8  # ASCIIZ string
    BLPAPI_DATATYPE_BYTEARRAY      = 9  # Opaque binary data
    BLPAPI_DATATYPE_DATE           = 10 # Date
    BLPAPI_DATATYPE_TIME           = 11 # Timestamp
    BLPAPI_DATATYPE_DECIMAL        = 12 #
    BLPAPI_DATATYPE_DATETIME       = 13 # Date and time
    BLPAPI_DATATYPE_ENUMERATION    = 14 # An opaque enumeration
    BLPAPI_DATATYPE_SEQUENCE       = 15 # Sequence type
    BLPAPI_DATATYPE_CHOICE         = 16 # Choice type
    BLPAPI_DATATYPE_CORRELATION_ID = 17 # Used for some internal
end

@enum SchemaStatus::Cint begin
    BLPAPI_STATUS_ACTIVE               = 0
    BLPAPI_STATUS_DEPRECATED           = 1
    BLPAPI_STATUS_INACTIVE             = 2
    BLPAPI_STATUS_PENDING_DEPRECATION  = 3
end

# valuetype field of CorrelationId
@enum CorrelationType::Cint begin
    BLPAPI_CORRELATION_TYPE_UNSET   = 0
    BLPAPI_CORRELATION_TYPE_INT     = 1
    BLPAPI_CORRELATION_TYPE_POINTER = 2
    BLPAPI_CORRELATION_TYPE_AUTOGEN = 3
end

# blpapi_defs.h
@enum EventType::Cint begin
    BLPAPI_EVENTTYPE_ADMIN                 = 1
    BLPAPI_EVENTTYPE_SESSION_STATUS        = 2
    BLPAPI_EVENTTYPE_SUBSCRIPTION_STATUS   = 3
    BLPAPI_EVENTTYPE_REQUEST_STATUS        = 4
    BLPAPI_EVENTTYPE_RESPONSE              = 5
    BLPAPI_EVENTTYPE_PARTIAL_RESPONSE      = 6
    BLPAPI_EVENTTYPE_SUBSCRIPTION_DATA     = 8
    BLPAPI_EVENTTYPE_SERVICE_STATUS        = 9
    BLPAPI_EVENTTYPE_TIMEOUT               = 10
    BLPAPI_EVENTTYPE_AUTHORIZATION_STATUS  = 11
    BLPAPI_EVENTTYPE_RESOLUTION_STATUS     = 12
    BLPAPI_EVENTTYPE_TOPIC_STATUS          = 13
    BLPAPI_EVENTTYPE_TOKEN_STATUS          = 14
    BLPAPI_EVENTTYPE_REQUEST               = 15
end

abstract type BLPResult end
abstract type BLPResultOk <: BLPResult end
abstract type BLPResultErr <: BLPResult end
abstract type BLPException <: Exception end

struct NamedTupleResult{T<:NamedTuple} <: BLPResultOk
    result::T
end

struct FieldDataVecResult{T} <: BLPResultOk
    field_exceptions::Dict{Symbol, Any}
    field_data_vec::Vector{T}
end

struct BLPTimeoutException <: BLPException
    timeout_milliseconds::UInt32
end

struct BLPResponseException <: BLPException
    msg::String
end

# Tells the user to file an issue on Github
struct BLPUnknownException <: BLPException
    cause::String
end

struct SecurityErr <: BLPResultErr
    security::String
    source::String
    code::Int32
    category::String
    subcategory::String
    message::String
end

struct FieldErr <: BLPResultErr
    security::String
    field::Symbol
    source::String
    code::Int32
    category::String
    subcategory::String
    message::String
end

struct BLPResultErrException{T<:BLPResultErr} <: BLPException
    err::T
end

"""
# Unwrap

If BLPResult holds an Err, panics.
If BLPResult holds Ok result, returns underlying result.

# NoUnwrap

Always returns the BLPResult itself.
"""
abstract type ErrorHandling end

struct Unwrap <: ErrorHandling
end

struct NoUnwrap <: ErrorHandling
end

"""
Represents the library version in use for Bloomberg API.

See [`get_version_info`](@ref).
"""
struct VersionInfo
    major::Int32
    minor::Int32
    patch::Int32
    build::Int32
end

mutable struct SessionOptions
    handle::Ptr{Cvoid}

    function SessionOptions(handle::Ptr{Cvoid})
        new_session_options = new(handle)
        finalizer(destroy!, new_session_options)
        return new_session_options
    end
end

function destroy!(opt::SessionOptions)
    if opt.handle != C_NULL
        blpapi_SessionOptions_destroy(opt.handle)
        opt.handle = C_NULL
    end
    nothing
end

abstract type AbstractEventSource end

mutable struct Session <: AbstractEventSource
    handle::Ptr{Cvoid}
    opened_services::Set{String}

    function Session(handle::Ptr{Cvoid}, opened_services::Set{String}, session_start_timeout_msecs::Integer, verbose::Bool)
        new_session = new(handle, opened_services)
        finalizer(destroy!, new_session)
        handle_session_start_events(new_session, session_start_timeout_msecs, verbose)
        return new_session
    end
end

function destroy!(session::Session)
    if session.handle != C_NULL
        blpapi_Session_destroy(session.handle)
        session.handle = C_NULL
    end
    nothing
end

struct CorrelationValue
    ptr::UInt
    user_data_array::NTuple{4, UInt}
    fun_ptr::UInt
end

mutable struct CorrelationId
    header::UInt32
    value::CorrelationValue
end

"""
A `Service` provides access to API data.
A service is obtained from a [`Session`](@ref)
and gives access to operations.

```julia
session = BLPData.Session()
service = session["//blp/refdata"]
println(BLPData.list_operation_names(service))
operation = service["HistoricalDataRequest"]
println(operation)
```
"""
mutable struct Service
    handle::Ptr{Cvoid}
    name::String
    session::Session
end

mutable struct BLPName
    handle::Ptr{Cvoid}
    symbol::Symbol

    function BLPName(handle::Ptr{Cvoid})
        ptr_check(handle, "Failed to create BLPName")
        sym = Symbol(unsafe_string(blpapi_Name_string(handle)))
        new_blp_name = new(handle, sym)
        finalizer(destroy!, new_blp_name)
        return new_blp_name
    end
end

function destroy!(name::BLPName)
    if name.handle != C_NULL
        blpapi_Name_destroy(name.handle)
        name.handle = C_NULL
    end
    nothing
end

# D is the value of datatype field
# T is the DataType for the Julia value
struct BLPConstant{D,T}
    name::BLPName
    description::String
    status::SchemaStatus
    datatype::BLPDataType
    value::T
end

struct BLPConstantList{D}
    name::BLPName
    description::String
    status::SchemaStatus
    datatype::BLPDataType
    list::Vector{BLPConstant}
end

"Wraps a blpapi_Datetime_t."
struct BLPDateTime
    parts::UInt8 # bitmask of date/time parts that are set
    hours::UInt8
    minutes::UInt8
    seconds::UInt8
    milliSeconds::UInt16
    month::UInt8
    day::UInt8
    year::UInt16
    offset::Int16 # (signed) minutes ahead of UTC
end

"""
    AbstractSchemaTypeDefinition{T}

Wraps a `blpapi_SchemaTypeDefinition_t` from BLPAPI.

The `T` type parameter is the value of the `datatype` field.

There are three concrete types for an `AbstractSchemaTypeDefinition`:

## Simple

A `SimpleSchemaTypeDefinition` is suitable when `datatype` implies a simple data type.
BLPAPI defines `blpapi_SchemaTypeDefinition_isSimple` to check for simple data type.

See also [`is_simple_datatype`](@ref).

## Complex

A `ComplexSchemaTypeDefinition` is suitable when `datatype` implies a compelx data type,
meaning either a `BLPAPI_DATATYPE_SEQUENCE` or `BLPAPI_DATATYPE_CHOICE`.
BLPAPI defines `blpapi_SchemaTypeDefinition_isComplex` to check for complex data type.

This struct has the same fields as for the Simple case,
with additional `elements::Vector{SchemaElementDefinition}` field.

See also [`is_complex_datatype`](@ref).

## Enumeration

A `EnumerationSchemaTypeDefinition` is suitable when `datatype` equals `BLPAPI_DATATYPE_ENUMERATION`.
BLPAPI defines `blpapi_SchemaTypeDefinition_isEnumeration` to check for enumeration data type.

An enumeration is also considered a simple datatype, so `blpapi_SchemaTypeDefinition_isSimple`
also returns a `true` value for this data type.

See also [`is_enumeration_datatype`](@ref).
"""
abstract type AbstractSchemaTypeDefinition{T} end

abstract type AbstractSchemaElementDefinition end

struct SchemaElementDefinition{T<:AbstractSchemaTypeDefinition} <: AbstractSchemaElementDefinition
    name::BLPName
    status::SchemaStatus
    schema_type::T
    alternate_names::Vector{BLPName}
    min_values::UInt64
    max_values::UInt64
end

"""
    SimpleSchemaTypeDefinition{T}

See docstring for [`AbstractSchemaTypeDefinition`](@ref).
"""
struct SimpleSchemaTypeDefinition{T} <: AbstractSchemaTypeDefinition{T}
    name::BLPName
    description::String
    status::SchemaStatus
    datatype::BLPDataType

    function SimpleSchemaTypeDefinition(name::BLPName, description::String, status::SchemaStatus, datatype::BLPDataType)
        @assert is_simple_datatype(datatype) "$datatype not a simple datatype."
        return new{datatype}(name, description, status, datatype)
    end
end

"""
    ComplexSchemaTypeDefinition{T}

See docstring for [`AbstractSchemaTypeDefinition`](@ref).
"""
struct ComplexSchemaTypeDefinition{T, D<:SchemaElementDefinition} <: AbstractSchemaTypeDefinition{T}
    name::BLPName
    description::String
    status::SchemaStatus
    datatype::BLPDataType
    elements::Vector{D}

    function ComplexSchemaTypeDefinition(name::BLPName, description::String, status::SchemaStatus, datatype::BLPDataType, elements::Vector{D}) where {D<:SchemaElementDefinition}
        @assert is_complex_datatype(datatype) "$datatype is not a complex datatype."
        return new{datatype, D}(name, description, status, datatype, elements)
    end
end

"""
    EnumerationSchemaTypeDefinition{T}

See docstring for [`AbstractSchemaTypeDefinition`](@ref).
"""
struct EnumerationSchemaTypeDefinition{T,L<:BLPConstantList} <: AbstractSchemaTypeDefinition{T}
    name::BLPName
    description::String
    status::SchemaStatus
    datatype::BLPDataType
    enumeration::L

    function EnumerationSchemaTypeDefinition(name::BLPName, description::String, status::SchemaStatus, datatype::BLPDataType, enumeration::L) where {L<:BLPConstantList}
        @assert datatype == BLPAPI_DATATYPE_ENUMERATION
        new{datatype, L}(name, description, status, datatype, enumeration)
    end
end

struct Operation{S<:SchemaElementDefinition}
    name::String
    request_definition::S
    response_definitions::Vector{AbstractSchemaElementDefinition}
end

mutable struct Request
    handle::Ptr{Cvoid}
    service::Service

    function Request(handle::Ptr{Cvoid}, service::Service)
        new_request = new(handle, service)
        finalizer(destroy!, new_request)
        return new_request
    end
end

function destroy!(req::Request)
    if req.handle != C_NULL
        blpapi_Request_destroy(req.handle)
        req.handle = C_NULL
    end
    nothing
end

# A holds the boolean result of `blpapi_Element_isArray`.
# D holds the value of `datatype` field.
abstract type AbstractElement{A,D} end

# A holds the boolean result of `blpapi_Element_isArray`.
# D holds the value of `datatype` field.
# T holds the type of the `source` field.
mutable struct Element{A,D} <: AbstractElement{A,D}
    handle::Ptr{Cvoid}
    name::BLPName
    datatype::BLPDataType
    is_array::Bool
    source::Any # Request or any other parent

    function Element(handle::Ptr{Cvoid}, source)
        ptr_check(handle, "Failed to create Element")
        name = BLPName(blpapi_Element_name(handle))
        datatype = BLPDataType(blpapi_Element_datatype(handle))
        is_array = blpapi_Element_isArray(handle) != 0

        return new{is_array,datatype}(handle, name, datatype, is_array, source)
    end
end

mutable struct Event
    handle::Ptr{Cvoid}
    event_type::EventType
    source::AbstractEventSource

    function Event(handle::Ptr{Cvoid}, source::AbstractEventSource)
        ptr_check(handle, "Failed to create Event")
        event_type = EventType(blpapi_Event_eventType(handle))
        new_event = new(handle, event_type, source)
        finalizer(destroy!, new_event)
        return new_event
    end
end

function destroy!(event::Event)
    if event.handle != C_NULL
        err = blpapi_Event_release(event.handle)
        error_check(err, "Failed to release event")
        event.handle = C_NULL
    end
    nothing
end

mutable struct Message
    handle::Ptr{Cvoid}
    correlation_ids::Vector{CorrelationId}

    function Message(handle::Ptr{Cvoid})
        ptr_check(handle, "Failed to create Message")

        num_corr_ids = blpapi_Message_numCorrelationIds(handle)
        if num_corr_ids == 0
            correlation_ids = Vector{CorrelationId}()
        else
            @assert num_corr_ids > 0
            correlation_ids = [ blpapi_Message_correlationId(handle, i) for i in 0:(num_corr_ids-1) ]
        end

        return new(handle, correlation_ids)
    end
end

mutable struct MessageIterator
    handle::Ptr{Cvoid}

    function MessageIterator(handle::Ptr{Cvoid})
        ptr_check(handle, "Failed to create MessageIterator")
        new_message_iterator = new(handle)
        finalizer(destroy!, new_message_iterator)
        return new_message_iterator
    end
end

function destroy!(msg_iter::MessageIterator)
    if msg_iter.handle != C_NULL
        blpapi_MessageIterator_destroy(msg_iter.handle)
        msg_iter.handle = C_NULL
    end
    nothing
end

"""
An `EventQueue` can be used in `send_request`.
The application can then handle responses
in an async fashion as they arrive, or
handle all responses synchronously.

Use `EventQueue()` to create a new queue.
"""
mutable struct EventQueue <: AbstractEventSource
    handle::Ptr{Cvoid}

    function EventQueue(handle::Ptr{Cvoid})
        ptr_check(handle, "Failed to create EventQueue")
        new_event_queue = new(handle)
        finalizer(destroy!, new_event_queue)
        return new_event_queue
    end
end

function destroy!(queue::EventQueue)
    if queue.handle != C_NULL
        err = blpapi_EventQueue_destroy(queue.handle)
        error_check(err, "Failed to destroy EventQueue")
        queue.handle = C_NULL
    end
    nothing
end

"""
A list of `SubscriptionTopic`s.
This struct supports the basic vector API.

# Examples

```julia
list = BLPData.SubscriptionList() # creates an empty list
append!(list, ["//blp/mktdata/ticker/PETR4 BS Equity?fields=BID,ASK", "//blp/mktdata/ticker/VALE3 BS Equity?fields=BID,ASK"])

for topic in list
    println(topic)
end
```

See also [`SubscriptionTopic`](@ref).
"""
mutable struct SubscriptionList
    handle::Ptr{Cvoid}

    function SubscriptionList(handle::Ptr{Cvoid})
        ptr_check(handle, "Failed to create SubscriptionList")
        new_sublist = new(handle)
        finalizer(destroy!, new_sublist)
        return new_sublist
    end
end

function destroy!(sublist::SubscriptionList)
    if sublist.handle != C_NULL
        blpapi_SubscriptionList_destroy(sublist.handle)
        sublist.handle = C_NULL
    end
    nothing
end

"""
Represents a Topic related to the subscription API.

# Fields

* `correlation_id`: unique identifier for tracking events in the event stream related to this subscription.

* `topic`: a valid subscription string for the BLPAPI.

See also [`subscribe`](@ref).
"""
struct SubscriptionTopic
    correlation_id::CorrelationId
    topic::String
end

Base.:(==)(t1::SubscriptionTopic, t2::SubscriptionTopic) = t1.correlation_id == t2.correlation_id && t1.topic == t2.topic
Base.hash(t::SubscriptionTopic) = hash(t.correlation_id) + hash(t.topic)

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using Oceananigans.Advection
using Oceananigans.Coriolis
using Oceananigans.Operators
using Oceananigans.TurbulenceClosures: ∂ⱼ_2ν_Σ₁ⱼ, ∂ⱼ_2ν_Σ₂ⱼ, ∂ⱼ_2ν_Σ₃ⱼ, ∇_κ_∇c

@inline squared(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2

@inline gh2(i, j, k, grid, h, g) = @inbounds g/2 * h[i, j, k]^2

"""
Compute the tendency for the x-directional transport, uh
"""
@inline function uh_solution_tendency(i, j, k, grid,
                                      gravitational_acceleration,
                                      advection,
                                      coriolis,
                                      closure,
                                      bathymetry,
                                      solution,
                                      tracers,
                                      diffusivities,
                                      forcings,
                                      clock)

    g = gravitational_acceleration

    return ( - div_hUu(i, j, k, grid, advection, solution)
             - ∂xᶠᵃᵃ(i, j, k, grid, gh2, solution.h, gravitational_acceleration)
             - x_f_cross_U(i, j, k, grid, coriolis, solution)
             + forcings.uh(i, j, k, grid, clock, merge(solution, tracers)))
end

"""
Compute the tendency for the y-directional transport, vh.
"""
@inline function vh_solution_tendency(i, j, k, grid,
                                      gravitational_acceleration,
                                      advection,
                                      coriolis,
                                      closure,
                                      bathymetry,
                                      solution,
                                      tracers,
                                      diffusivities,
                                      forcings,
                                      clock)

     g = gravitational_acceleration

    return ( - div_hUv(i, j, k, grid, advection, solution)
             - ∂yᵃᶠᵃ(i, j, k, grid, gh2, solution.h, gravitational_acceleration)
             - y_f_cross_U(i, j, k, grid, coriolis, solution)
             + forcings.vh(i, j, k, grid, clock, merge(solution, tracers)))
end

"""
Compute the tendency for the height, h.
"""
@inline function h_solution_tendency(i, j, k, grid,
                                     gravitational_acceleration,
                                     coriolis,
                                     closure,
                                     bathymetry,
                                     solution,
                                     tracers,
                                     diffusivities,
                                     forcings,
                                     clock)

    return ( - div_Uh(i, j, k, grid, solution)
             + forcings.h(i, j, k, grid, clock, merge(solution, tracers)))
end

@inline function tracer_tendency(i, j, k, grid,
                                 val_tracer_index::Val{tracer_index},
                                 advection,
                                 closure,
                                 solution,
                                 tracers,
                                 diffusivities,
                                 forcing,
                                 clock) where tracer_index

    @inbounds c = tracers[tracer_index]

    return ( -  div_Uc(i, j, k, grid, advection, solution, c)
             + c_div_U(i, j, k, grid, solution, c)
             +  ∇_κ_∇c(i, j, k, grid, clock, closure, c, val_tracer_index, diffusivities, tracers, nothing)
             + forcing(i, j, k, grid, clock, merge(solution, tracers))
            )
end


text:
"""
	WattsStrogatzNetwork <: AbstractNetwork

A network contructed via the Watts-Strogatz method, possibly small-world.

# Fields
- N :: Integer -> Number of nodes
- k :: Integer -> Mean connectivity
- β :: Real -> Connection probability
- seed :: Integer -> Seed for the random creation
"""
mutable struct WattsStrogatzNetwork <: AbstractNetwork
	N :: Integer
	k :: Integer
	β :: Union{Real, Nothing}
	seed :: Integer

	_adjMat :: AbstractMatrix
	_props :: NetworkProperties
	_numShortcuts :: Integer
end

"""
	WattsStrogatzNetwork(N, k, β; directed=false, seed=-1)

Create an Watts-Strogatz network with `N` nodes, mean connectivity `k` and rewiring probability `β`.
"""
function WattsStrogatzNetwork(
	N::Integer,
	k::Integer,
	β::Real;
	directed::Bool=false,
	seed::Integer=-1
)
	(N <= 0) && throw(ArgumentError("Number of nodes must be a positive integer!"))

	mat = nothing
	numShortcuts = nothing
	(k % 2 == 1) && (k -= 1)

	# Deal with special cases
	if k < 0
		k = 0
		mat = adjMat(EmptyNetwork(N), sparse=true)
		numShortcuts = 0
		β = nothing
	elseif k >= N-1
		k = N-1
		mat = adjMat(GlobalNetwork(N))
		numShortcuts = 0
		β = nothing
	else
		(β < 0) && (β = 0.0)
		(β > 1) && (β = 1.0)
	end

	_seed = seed < 0 ? rand(1:999999999) : seed
	isnothing(mat) && ((mat, numShortcuts) = generateWSAdjMat(N, k, β; directed=directed, seed=_seed))

	WattsStrogatzNetwork(N, k, β,_seed, mat, NetworkProperties(N, directed), numShortcuts)
end

"""
	WattsStrogatzNetwork(N, k, numShortcuts; directed=false, seed=-1)

Create an Watts-Strogatz network with `N` nodes, mean connectivity `k` and `numShortcuts` rewired
shortcuts.
"""
function WattsStrogatzNetwork(
	N::Integer,
	k::Integer,
	numShortcuts::Integer;
	directed::Bool=false,
	seed::Integer=-1
)
	(N <= 0) && throw(ArgumentError("Number of nodes must be a positive integer!"))

	mat = nothing
	(k % 2 == 1) && (k -= 1)

	# Deal with special cases
	if k < 0
		k = 0
		mat = adjMat(EmptyNetwork(N))
		numShortcuts = 0
	elseif k >= N-1
		k = N-1
		mat = adjMat(GlobalNetwork(N))
		numShortcuts = 0
	else
		numConnections = k * N
		directed || (numConnections = Int(numConnections / 2))
		(numShortcuts < 0) && (numShortcuts = 0)
		(numShortcuts > numConnections) && (numShortcuts = numConnections)
	end

	_seed = seed < 0 ? rand(1:999999999) : seed
	isnothing(mat) && (mat = generateWSAdjMat(N, k, numShortcuts; directed=directed, seed=_seed))

	WattsStrogatzNetwork(N, k, nothing, _seed, mat, NetworkProperties(N, directed), numShortcuts)
end

function show(network::WattsStrogatzNetwork)
	print(isdirected(network) ? "Directed " : "Undirected ")
	println("Watts-Strogatz Network")
	println("- N = $(network.N)")
	println("- k = $(network.k)")
	println("- β = $(network.β)")
	println("- seed = $(network.seed)")
end


text:
"""
    Berlage(; <keyword arguments>)

Create a Berlage wavelet.

# Arguments
- `dt::AbstractFloat=0.002`: sampling interval in secs.
- `f0::Real=20.0`: central frequency in Hz.
- `m::Real=2`: exponential parameter of Berlage wavelet.
- `alpha::Real=180.0`: alpha parameter of Berlage wavelet in rad/secs.
- `phi0::Real=0.0`: phase rotation in radians.

# Example
```julia
julia> w = Berlage(); plot(w);
```
**Reference**
* Aldridge, David F., 1990, The berlage wavelet: GEOPHYSICS, 55, 1508--1511.
"""
function Berlage(; dt::AbstractFloat=0.002, f0::Real=20.0, m::Real=2, alpha::Real=180.0,
                 phi0::Real=0.0)

    nw = floor(Int, 2.0/(f0*dt))
    t = dt*collect(0:1:nw-1)
    w = [ (t[i] ^ m) * exp(-alpha*t[i]) * cos(2*pi*f0*t[i] + phi0) for i in eachindex(t) ]

    w = w/maximum(w)

end


text:
using Pkg
Pkg.activate("examples")
Pkg.instantiate()

import GPLinearODEMaker; GLOM = GPLinearODEMaker

kernel, n_kern_hyper = GLOM.include_kernel("pp")

n = 100
xs = 20 .* sort(rand(n))
y1 = sin.(xs) .+ 0.1 .* randn(n)
noise1 = 0.1 .* ones(n)
y2 = cos.(xs) .+ 0.2 .* randn(n)
noise2 = 0.2 .* ones(n)

ys = collect(Iterators.flatten(zip(y1, y2)))
noise = collect(Iterators.flatten(zip(noise1, noise2)))

prob_def = GLOM.GLO(kernel, n_kern_hyper, 2, 2, xs, ys; noise = noise, a0=[[1. 0];[0 1]])
total_hyperparameters = append!(collect(Iterators.flatten(prob_def.a0)), [10])
total_hyperparameters .*= (1 .+ 0.1 .* rand(length(total_hyperparameters)))
workspace = GLOM.nlogL_matrix_workspace(prob_def, total_hyperparameters)

using Optim

initial_x = GLOM.remove_zeros(total_hyperparameters)
f(non_zero_hyper::Vector{T} where T<:Real) = GLOM.nlogL_GLOM!(workspace, prob_def, non_zero_hyper)
function g!(G::Vector{T}, non_zero_hyper::Vector{T}) where T<:Real
    G[:] = GLOM.∇nlogL_GLOM!(workspace, prob_def, non_zero_hyper)
end
function h!(H::Matrix{T}, non_zero_hyper::Vector{T}) where T<:Real
    H[:, :] = GLOM.∇∇nlogL_GLOM!(workspace, prob_def, non_zero_hyper)
end

function optim_cb(x::OptimizationState)
    println()
    if x.iteration > 0
        println("Iteration:     ", x.iteration)
        println("Time so far:   ", x.metadata["time"], " s")
        println("Currnet score: ", x.value)
        println("Gradient norm: ", x.g_norm)
        println()
    end
    return false
end

# @time result = optimize(f, initial_x, NelderMead()) # 26s
# @time result = optimize(f, g!, initial_x, LBFGS()) # 40s
@time result = optimize(f, g!, h!, initial_x, NewtonTrustRegion(), Optim.Options(callback=optim_cb, g_tol=1e-6, iterations=50)) # 15s

fit_total_hyperparameters = GLOM.reconstruct_total_hyperparameters(prob_def, result.minimizer)

n_samp_points = convert(Int64, max(500, round(2 * sqrt(2) * length(prob_def.x_obs))))
x_samp = collect(range(minimum(prob_def.x_obs); stop=maximum(prob_def.x_obs), length=n_samp_points))
n_total_samp_points = n_samp_points * prob_def.n_out
n_meas = length(prob_def.x_obs)

mean_GP, σ, mean_GP_obs, Σ = GLOM.GP_posteriors(prob_def, x_samp, fit_total_hyperparameters; return_mean_obs=true)

#=
n_show = 5
show_curves = zeros(n_show, n_total_samp_points)
L = GLOM.ridge_chol(Σ).L
for i in 1:n_show
    show_curves[i, :] = L * randn(n_total_samp_points) + mean_GP
end

using Plots

function make_plot(output::Integer; show_draws::Bool=true)
    sample_output_indices = output:prob_def.n_out:n_total_samp_points
    obs_output_indices = output:prob_def.n_out:length(ys)
    p = scatter(xs, ys[obs_output_indices], yerror=noise1, leg=false)
    if show_draws
        for i in 1:n_show
            plot!(x_samp, show_curves[i, sample_output_indices], leg=false)
        end
    end
    plot!(x_samp, mean_GP[sample_output_indices]; ribbon=σ[sample_output_indices], alpha=0.3, leg=false)
    return p
end

plot(make_plot(1), make_plot(2), layout=(2,1), size=(1920,1080))
=#


text:
module GBQ

# external library dependancies
using JSON
using DataFrames

# functions to export
export gbq_large_results, gbq_show, gbq_head, gbq_list_tables,
gbq_list_datasets, gbq_list_projects, gbq_query, gbq_create_dataset,
gbq_copy_table, gbq_table_exists

# custom exception for GBQ package
struct GBQException <: Exception
    var::String
end

Base.showerror(io::IO, e::GBQException) = print(io, e.var);

# runs a query and saves results directly into a table
# query = the query to run
# destination = destination table (must not already exist)
# table = dataset to store table in
#
# returns DataFrame conatining a preview of the new table
function gbq_large_results(query, project, dataset, table)
  if gbq_table_exists(dataset, table)
    throw(GBQException("Destination Table Already Exists"))
  end
  response = JSON.parse(read(`bq --format=json --quiet=True query --destination_table=$project:$dataset.$table "$query"`, String))
  return _gbq_parse(response)
end

# check if table exists in dataset
#
# returns boolean
function gbq_table_exists(dataset, table)
    tables = gbq_list_tables(dataset)
    return table in tables[:tableId]
end

# internal function to parse json response returned from big query
#
# returns a dataframe
function _gbq_parse(response)
    cols = collect(keys(response[1]))
    values = Dict()
    for key in cols
        values[key] = []
    end
    for dict in response
        for key in cols
            push!(values[key], dict[key])
        end
    end
    return convert(DataFrame, values)
end

# Execute a query
#
# Returns a dataframe
function gbq_query(query; use_legacy_sql=false, quiet=true, max_rows=100000000)
  response = JSON.parse(read(`bq --format=json  --quiet="$quiet" query --use_legacy_sql="$use_legacy_sql" --max_rows="$max_rows" "$query"`, String))
  return _gbq_parse(response)
end

# list all projects associated with default billing account
#
# Returns a dataframe
function gbq_list_projects()
  response = JSON.parse(read(`bq ls --format=json -p`, String))
  return _gbq_parse(response)
end

# lists datasets in default project
#
# Returns a dataframe
function gbq_list_datasets()
  response = JSON.parse(read(`bq ls --format=json`, String))
  return _gbq_parse(response)
end

# lists all tables in a specified dataset
#
# Returns a dataframe
function gbq_list_tables(dataset)
  response = JSON.parse(read(`bq ls --format=json "$dataset"`, String))
  return _gbq_parse(response)
end

# preview the head of the table
#
# Returns a dataframe
function gbq_head(dataset, table, num_rows=10)
  response = JSON.parse(read(`bq --format=json head -n $num_rows $dataset.$table`, String))
  return _gbq_parse(response)
end

# examine the schema for a table
#
# Returns dict containing schema
function gbq_show(dataset, table)
  return JSON.parse(read(`bq --format=json show $dataset.$table`, String))
end

# create a new dataset in the default project
#
# returns string containing response from Google BigQuery
function gbq_create_dataset(dataset)
  return read(`bq mk $dataset`, String)
end

# copy a table
#
# returns string containing response from Google BigQuery
function gbq_copy_table(dataset1, table1, dataset2, table2)
  return read(`bq cp --quiet=True $dataset1.$table1 $dataset2.$table2`, String)
end

end

max_stars_repo_name: martineastwood/GBQ

max_issues_repo_name: martineastwood/GBQ

max_forks_repo_name: martineastwood/GBQ


text:
"""
    P(n, τ)

Returns the value of the Legendre polynomial of degree `n` at `τ`.
(Recursive function used in `../test`).
"""
function P(n::Integer, τ::T) where T <: AbstractFloat
    if n == 0
        return one(T)
    elseif n == 1
        return τ
    else
        return ( (2n-1) * τ * P(n-1, τ) - (n-1) * P(n-2, τ) ) / n
    end
end

"""
    legendre_polys!(P, τ)

If `τ` is a scalar, then `P` is a vector with `P[n+1]` equal to the value
of the Legendre polynomial of degree `n` at `τ`.

If `τ` is a vector, then `P` is a matrix with `P[n+1,j]` equal to the value
of the Legendre polynomial of degree `n` at `τ[j]`.
"""
function legendre_polys!(P::AbstractVector{T}, τ::T) where T <: AbstractFloat
    I = length(P)
    if I ≥ 1
        P[1] = one(T)
    end
    if I ≥ 2
        P[2] = τ
    end
    for j = 1:I, i = 1:I-2
        P[i+2] = ((2i+1) * τ * P[i+1] - i * P[i] ) / (i+1)
    end
end

function legendre_polys!(P::AbstractMatrix{T}, τ::AbstractVector{T}
                        ) where T <: AbstractFloat
    I = size(P, 1)
    J = size(P, 2)
    @argcheck length(τ) == J
    if I ≥ 1
        for j = 1:J
            P[1,j] = one(T)
        end
    end
    if I ≥ 2
        for j = 1:J
            P[2,j] = τ[j]
        end
    end
    for j = 1:J, i = 1:I-2
        P[i+2,j] = ((2i+1) * τ[j] * P[i+1,j] - i * P[i,j] ) / (i+1)
    end
end

"""
    dP(n, τ)

Returns the value of the derivative of the Legendre polynomial of degree `n`
at `τ`.  (Recursive function used in `../test`).
"""
function dP(n::Integer, τ::T) where T <: AbstractFloat
    if n == 0
        return zero(T)
    elseif n == 1
        return one(T)
    else
        return ( (2n-1) * τ * dP(n-1, τ) - n * dP(n-2, τ) ) / ( n - 1 )
    end
end

function deriv_legendre_polys!(dP::AbstractVector{T}, τ::T
                              ) where T <: AbstractFloat
    I = length(dP)
    if I ≥ 1
        dP[1] = zero(T)
    end
    if I ≥ 2
        dP[2] = one(T)
    end
    for i = 1:I-2
        dP[i+2] = ( (2i+1) * τ * dP[i+1] - (i+1) * dP[i] ) / i
    end
end

function deriv_legendre_polys!(dP::AbstractMatrix{T}, τ::AbstractVector{T}
                               ) where T <: AbstractFloat
    I = size(dP, 1)
    J = size(dP, 2)
    @argcheck length(τ) == J
    if I ≥ 1
        for j = 1:J
            dP[1,j] = zero(T)
        end
    end
    if I ≥ 2
        for j = 1:J
            dP[2,j] = one(T)
        end
    end
    for j = 1:J, i = 1:I-2
        dP[i+2,j] = ((2i+1) * τ[j] * dP[i+1,j] - (i+1) * dP[i,j] ) / i
    end
end

"""
    pcwise_t, pcwise_U = evaluate_pcwise_poly(U, t, ppI, store)
"""
function evaluate_pcwise_poly!(U::Vector{Vector{T}}, t::OffsetVector{T},
                               ppI::Integer, store::Store{T}
                              ) where { T <: AbstractFloat }
    τ = range(-one(T), stop=one(T), length=ppI)
    pcwise_t, pcwise_U = evaluate_pcwise_poly!(U, t, τ, store)
    return pcwise_t, pcwise_U
end

function evaluate_pcwise_poly!(U::Vector{Vector{T}}, t::OffsetVector{T},
                               τ::AbstractVector{T}, store::Store{T}
                              ) where T <: AbstractFloat
    N = length(U)
    ppI = length(τ)
    @argcheck ppI ≤ store.ppImax
    pcwise_t = Array{T}(undef, ppI, N)
    pcwise_U = similar(pcwise_t)
    evaluate_pcwise_poly!(pcwise_t, pcwise_U, U, t, τ, store)
    return pcwise_t, pcwise_U
end

function evaluate_pcwise_poly!(pcwise_t::Matrix{T}, pcwise_U::Matrix{T},
                               U::Vector{Vector{T}}, t::OffsetVector{T},
                               τ::AbstractVector{T}, store::Store{T}
                              ) where { T <: AbstractFloat }
    N = length(U)
    pts_per_interval = length(τ)
    @argcheck size(pcwise_t) == (pts_per_interval, N)
    @argcheck size(pcwise_U) == (pts_per_interval, N)
    @argcheck length(t) == N+1
    rmax = store.rmax
    Ψ = view(store.Ψ, 1:rmax, 1:pts_per_interval)
    legendre_polys!(Ψ, τ)
    for n = 1:N
        rn = length(U[n])
        for m = 1:pts_per_interval
            pcwise_t[m,n] = ( (1-τ[m])*t[n-1] + (1+τ[m])*t[n] ) / 2
            s = zero(T)
            for j = 1:rn
                s += U[n][j] * Ψ[j,m]
            end
            pcwise_U[m,n] = s
        end
    end
end

function Fourier_Legendre_coefs(u::Function, r::Integer, t::OffsetArray{T},
                                store::Store{T}) where T <: AbstractFloat
    N = length(t) - 1
    rmax = store.rmax
    a = Vector{Vector{T}}(undef, N)
    M = store.Mmax
    τ, wτ = rule(store.legendre[M])
    Ψ = view(store.Ψ, 1:r, 1:M)
    legendre_polys!(Ψ, τ)
    for n = 1:N
        a[n] = Vector{T}(undef, r)
        for j = 1:r
            s = zero(T)
            for m = 1:M
                tnm = ( (1-τ[m])*t[n-1] + (1+τ[m])*t[n] ) / 2
                s += wτ[m] * u(tnm) * Ψ[j,m]
            end
            a[n][j] = (2j-1) * s / 2
        end
    end
    return a
end

function reconstruction_pts(In::Tuple{T,T}, r::Integer) where T <: AbstractFloat
    tnm1, tn = In
    tstar = reconstruction_pts(T, r)
    tstar[0] = tnm1
    for m = 1:r-1
        tstar[m] = ( ( 1 - tstar[m] ) * tnm1 + ( tstar[m] + 1 ) * tn ) / 2
    end
    tstar[r] = tn
    return tstar
end

function reconstruction_pts(::Type{T}, r::Integer) where T <: AbstractFloat
    right = GaussQuadrature.right
    τ, wτ = GaussQuadrature.legendre(T, r, right)
    return OffsetArray([ -one(0); τ ], 0:r)
end

"""
    Uhat = reconstruction(U, u0, store)
"""
function reconstruction(U::Vector{Vector{T}}, u0::T,
                        store::Store{T}) where T <: AbstractFloat
    Uhat = reconstruction(U, u0, store.rmax)
    return Uhat
end

"""
    Uhat = reconstruction(U, u0, rmax)
"""
function reconstruction(U::Vector{Vector{T}}, u0::T,
	                rmax::Integer) where T <: AbstractFloat
    N = length(U)
    Uhat = Vector{Vector{T}}(undef, N)
    pow = OffsetArray{T}(undef, 0:rmax)
    pow[0] = one(T)
    for n = 1:rmax
        pow[n] = -pow[n-1]
    end
    r1 = length(U[1])
    @argcheck r1+1 ≤ rmax
    U_left = zero(T)   # U_left  = U(t[0] + 0) = value at the left endpoint.
    U_right = zero(T)  # U_right = U(t[1] - 0) = value at the right endpoint.
    for j = 1:r1
        U_left += pow[j-1] * U[1][j]
        U_right += U[1][j]
    end
    jumpU0 = U_left - u0
    Uhat[1] = Vector{T}(undef, r1+1)
    Uhat[1][1:r1] .= U[1][1:r1]
    Uhat[1][r1] += pow[r1] * jumpU0 / 2
    Uhat[1][r1+1] = - pow[r1] * jumpU0 / 2
    for n = 2:N
        rn = length(U[n])
        @argcheck rn + 1 ≤ rmax
        U_left = zero(T) # U_left = U(t[n-1] + 0 )
        for j = 1:rn
            U_left += pow[j-1] * U[n][j]
        end
        jumpUnm1 = U_left - U_right
        Uhat[n] = Vector{T}(undef, rn+1)
        Uhat[n][1:rn] .= U[n][1:rn]
        Uhat[n][rn] += pow[rn] * jumpUnm1 / 2
        Uhat[n][rn+1]  = -pow[rn] * jumpUnm1 / 2
        U_right = zero(T) # U_right = U(t[n] - 0 )
        for j = 1:rn
            U_right += U[n][j]
        end
    end
    return Uhat
end

function max_order(U::Vector{Vector{T}}) where T <: AbstractFloat
    N = length(U)
    r = length(U[1])
    for n = 2:N
        rn = length(U[n])
        r = max(r, rn)
    end
    return r
end

function jumps(U::Vector{Vector{T}}, t::OffsetVector{T},
	u0::T) where T <: AbstractFloat
    N = length(t) - 1
    r = max_order(U)
    JU = OffsetVector{T}(undef, 0:N-1)
    U_left = u0
    U_right = zero(T)
    pow = ones(T, r)
    pow[2:2:r] .= -one(T)
    for n = 1:N
	U_left = zero(T)
	for j = 1:r
	    U_left += pow[j] * U[n][j]
	end
	JU[n-1] = U_left
    end
    JU[0] -= u0
    for n = 2:N
	U_right = sum(U[n-1])
	JU[n-1] -= U_right
    end
    return JU
end

function dG_error_estimator(JU::OffsetVector{T}, t::OffsetVector{T},
        r::Integer, τ::AbstractVector{T}) where T <: AbstractFloat
    N = length(t) - 1
    pts_per_interval = length(τ)
    pcwise_t = Matrix{T}(undef, pts_per_interval, N)
    pcwise_approx_err = similar(pcwise_t)
    Ψ = Matrix{T}(undef, r+1, pts_per_interval)
    legendre_polys!(Ψ, τ)
    if r % 2 == 0
        sign = 1
    else
        sign = -1
    end
    for n = 1:N
        for m = 1:pts_per_interval
            pcwise_t[m,n] = ( (1-τ[m])*t[n-1] + (1+τ[m])*t[n] ) / 2
            pcwise_approx_err[m,n] = (sign/2) * JU[n-1] * ( Ψ[r+1,m] - Ψ[r,m] )
        end
    end
    return pcwise_t, pcwise_approx_err
end

function dG_error_estimator(JU::OffsetVector{T}, t::OffsetVector{T},
        r::Integer, pts_per_interval::Integer) where T <: AbstractFloat
    τ = range(-1.0, 1.0, length=pts_per_interval)
    return dG_error_estimator(JU, t, r, τ)
end


text:
struct Sphere <: Hitable
	center::Point3
	radius::Float64
	material::Material
end

export Sphere, Hit

function hit(s::Sphere, t, ray)
	p = at(ray, t)
	outward_normal = (p - s.center) / s.radius
	ff = dot(ray.direction, outward_normal)
	norm = ff < 0 ? outward_normal : -outward_normal
	Hit(p, norm, t, ff < 0)
end

function trace(sphere::Sphere, ray::Ray, t_min::Float64, t_max::Float64)::Float64
	oc = ray.origin - sphere.center
	a = magnitudesq(ray.direction)
	half_b = dot(oc, ray.direction)
	c = magnitudesq(oc) - sphere.radius^2
	discriminant = half_b^2 - a*c
	if discriminant < 0
		return -1.0
	end

	sqrtd = sqrt(discriminant)
	root = (-half_b - sqrtd) / a
	if t_max < root || root < t_min
		root = (-half_b + sqrtd) / a
		if t_max < root || root < t_min
			return -1.0
		end
	end

	root
end


text:
Nor(x;m=0.0,s=1.0) = 1.0/sqrt(2*pi)/s * exp( -0.5/(s^2)*(x-m)^2 )
Nor2(x;m1=0.0,m2=0.0,s=1.0) = 0.7/sqrt(2*pi)/s * exp( -0.5/(s^2)*(x-m1)^2 ) + 0.3/sqrt(2*pi)/s * exp( -0.5/(s^2)*(x-m2)^2 )
function simdata01_true(;N1=100,N2=100)
  # common struct to hold all simulated data
  data = DataGroundTruth()
  samples = SampleData(zeros(0,0), 1, N1+N2)

  temp1 = (randn(round(Int,0.7*N1)).-1.0)
  cluster1 = ([temp1; randn(N1-length(temp1)).-5.0]')
  cluster2 = randn(1,N2).+1

  data.clusters[1] = cluster1
  data.clustersizes[1] = N1
  data.clusters[2] = cluster2
  data.clustersizes[2] = N2
  # one dimensional, Array{Float64,1}, with N1+N2 elements [vector alias]
  sd = [cluster1';cluster2']'
  permu = randperm(N1+N2)
  samples.samples = sd[:,permu]

y = randn(10)

  data.assignment = ([ones(Int,N1);2*ones(Int,N2)])[permu]

  # [OPTIONAL] going to carry plot layers for convenience also
  drawN = 200
  pts1 = StatsBase.sample(cluster1[1,:][:], ProbabilityWeights(1.0/N1*ones(N1)), drawN, replace=false)
  yy1 = zeros(drawN)
  for i in 1:drawN  yy1[i] = Nor2(pts1[i],m1=-1.0,m2=-5.0) end
  data.plt_lyr_cluster[1] = layer(x=pts1, y=yy1, Geom.line, Theme(default_color=colorant"blue"))[1] #parse(Colorant,"red")
  data.plt_lyr_cluster_nocolor[1] = layer(x=pts1, y=yy1, Geom.line, Theme(default_color=colorant"gray"))[1] #parse(Colorant,"red")

  pts2 = StatsBase.sample(cluster2[1,:][:], ProbabilityWeights(1.0/N2*ones(N2)), drawN, replace=false)
  yy2 = zeros(drawN)
  for i in 1:drawN  yy2[i] = Nor(pts2[i],m=+1.0) end
  data.plt_lyr_cluster[2] = layer(x=pts2, y=yy2, Geom.line, Theme(default_color=colorant"red"))[1] #parse(Colorant,"red")
  data.plt_lyr_cluster_nocolor[2] = layer(x=pts2, y=yy2, Geom.line, Theme(default_color=colorant"gray"))[1] #parse(Colorant,"red")

  return samples, data
end

# two dimensional example
function simdata02_true(;N1=100,N2=100)
  # common struct to hold all simulated data
  data = DataGroundTruth()
  samples = SampleData(zeros(0,0), 1, N1+N2)

  cluster1 = rand(MvNormal([3.0;3.0],[2.5 -1.5; -1.5 2.5]),N1)
  cluster2 = rand(MvNormal([-2.0;-2.0],[4.0 3.0; 3.0 4.0]),N2)

  data.clusters[1] = cluster1
  data.clustersizes[1] = N1
  data.clusters[2] = cluster2
  data.clustersizes[2] = N2
  # two rows, N1+N2 columns (Array{Float64,2})
  # data.samples = [cluster1';cluster2']'
  # TODO: not sure if sd should be a vector or matrix
  sd = [cluster1';cluster2']'
  permu = randperm(N1+N2)
  samples.samples = sd[:,permu]

  data.assignment = ([ones(Int,N1);2*ones(Int,N2)])[permu]

  drawN = 200

  pts1 = StatsBase.sample(cluster1[1,:][:], StatsBase.ProbabilityWeights(1.0/N1*ones(N1)), drawN, replace=false)
  data.plt_lyr_cluster[1] = layer(x=cluster1[1,:], y=cluster1[2,:], Geom.point, Theme(default_color=colorant"blue"))
  data.plt_lyr_cluster_nocolor[1] = layer(x=cluster1[1,:], y=cluster1[2,:], Geom.point, Theme(default_color=colorant"gray"))
  # data.plt_lyr_cluster[1] = layer(x=pts1, y=yy1, Geom.point, Theme(default_color=colorant"blue"))[1] #parse(Colorant,"red")
  # data.plt_lyr_cluster_nocolor[1] = layer(x=pts1, y=yy1, Geom.point, Theme(default_color=colorant"gray"))[1] #parse(Colorant,"red")

  pts2 = StatsBase.sample(cluster2[1,:][:], StatsBase.ProbabilityWeights(1.0/N2*ones(N2)), drawN, replace=false)
  data.plt_lyr_cluster[2] = layer(x=cluster2[1,:], y=cluster2[2,:], Geom.point, Theme(default_color=colorant"red"))
  data.plt_lyr_cluster_nocolor[2] = layer(x=cluster2[1,:], y=cluster2[2,:], Geom.point, Theme(default_color=colorant"gray"))
  # data.plt_lyr_cluster[2] = layer(x=pts2, y=yy2, Geom.point, Theme(default_color=colorant"red"))[1] #parse(Colorant,"red")
  # data.plt_lyr_cluster_nocolor[2] = layer(x=pts2, y=yy2, Geom.point, Theme(default_color=colorant"gray"))[1] #parse(Colorant,"red")

  return samples, data
end


text:
module TagBot

using Dates: Day, UTC, now

using GitHub: GitHub
using JSON: JSON

const GH = GitHub

const AUTH = Ref{GH.OAuth2}()
const ISSUE_TITLE = "TagBot trigger issue"
const ISSUE_BODY = """
This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your `TagBot.yml` to include issue comment triggers.
Please see [this post on Discourse](https://discourse.julialang.org/t/ann-required-updates-to-tagbot-yml/49249) for instructions and more details.
"""

include("cron.jl")
include("pull_request.jl")

function main()
    AUTH[] = GH.authenticate(ENV["GITHUB_TOKEN"])
    event = JSON.parse(read(ENV["GITHUB_EVENT_PATH"], String))
    if is_merged_pull_request(event)
        handle_merged_pull_request(event)
    elseif is_cron(event)
        handle_cron(event)
    end
end

function repo_and_version_of_pull_request_body(body)
    if occursin("JLL package", body)
        @info "Skipping JLL package registration"
        return nothing, nothing
    end
    m = match(r"Repository: .*github\.com[:/](.*)", body)
    repo = m === nothing ? nothing : strip(m[1])
    repo !== nothing && endswith(repo, ".git") && (repo = repo[1:end-4])
    m = match(r"Version: (.*)", body)
    version = m === nothing ? nothing : strip(m[1])
    return repo, version
end

function clone_repo(repo)
    dir = mktempdir()
    run(`git clone --depth=1 https://github.com/$repo $dir`)
    return dir
end

function is_tagbot_enabled(repo)
    # TODO: Traversing the file tree should be possible via GitHub API,
    # but GitHub.jl doesn't seem capable.
    dir = clone_repo(repo)
    workflows = joinpath(dir, ".github", "workflows")
    isdir(workflows) || return false
    for workflow in readdir(workflows)
        contents = read(joinpath(workflows, workflow), String)
        occursin("JuliaRegistries/TagBot", contents) && return true
    end
    return false
end

function get_repo_notification_issue(repo)
    # TODO: Get the authenticated user (how?) and use it as `creator`.
    params = (; creator="JuliaTagBot", state="closed")
    issues, _ = GH.issues(repo; auth=AUTH[], params=params)
    filter!(x -> x.pull_request === nothing, issues)
    return if isempty(issues)
        @info "Creating new notification issue"
        params = (; title=ISSUE_TITLE, body=ISSUE_BODY)
        issue = GH.create_issue(repo; auth=AUTH[], params=params)
        GH.edit_issue(repo, issue; auth=AUTH[], params=(; state="closed"))
        issue
    else
        @info "Found existing notification issue"
        issues[1]
    end
end

function notification_body(event)
    url = get(get(event, "pull_request", Dict()), "html_url", "")
    body = "Triggering TagBot for merged registry pull request"
    return isempty(url) ? body : "$body: $url"
end

function notify(repo, issue, body)
    return GH.create_comment(repo, issue, :issue; auth=AUTH[], params=(; body=body))
end

function tag_exists(repo, version)
    return try
        GH.tag(repo, version; auth=AUTH[])
        true
    catch e
        if !occursin("404", e.msg)
            @warn "Unknown error when checking for existing tag" ex=(e, catch_backtrace())
        end
        false
    end
end

function maybe_notify(event, repo, version; check_tag=false)
    @info "Processing version $version of $repo"
    if !is_tagbot_enabled(repo)
        @info "TagBot is not enabled on $repo"
        return
    end
    if check_tag && tag_exists(repo, version)
        @info "Tag $version already exists for $repo"
        return
    end
    issue = get_repo_notification_issue(repo)
    body = notification_body(event)
    notify(repo, issue, body)
end

end


text:
# This file was generated, do not modify it. # hide
coerce!(data, :Wind_gen => Continuous)
coerce!(data, :Solar_gen => Continuous)

schema(data)


text:
# using Calculus
# define types first
struct Substance
    name::String
    density::Function
    enthalpy::Function
    heat_capacity::Function
    viscosity::Function
    thermal_conductivity::Function
    alt_heat_capacity::Function
    visc_kin::Function
    exp_coeff::Function
    therm_diff::Function
    stens::Function
    Pr::Function
end

function Substance(name, x::Vector{Function})
    dens, enth, Cp, visc, thermCond, stens = x
    function alt_heat_capacity(p, T)
        return ForwardDiff.derivative(T_ -> enth(p, T_), T)
    end
    function visc_kin(p, T)
        return visc(p, T)/dens(p, T)
    end
    function exp_coeff(p, T)
        dRdT = ForwardDiff.derivative(T_ -> dens(p, T_), T)
        return (-1/dens(p, T))*dRdT
    end
    function thermDiff(p, T)
        if Cp == nothing
            Cp_ = alt_heat_capacity(p, T)
        else
            Cp_ = Cp(p, T)
        end
        return thermCond(p, T) / (dens(p,T)*Cp_)
    end
    function Pr(p, T)
        if Cp == nothing
            Cp_ = alt_heat_capacity(p, T)
        else
            Cp_ = Cp(p, T)
        end
        return Cp_ * visc(p, T) / thermCond(p, T)
    end
    return Substance(name, dens, enth, Cp, visc, thermCond, alt_heat_capacity, visc_kin, exp_coeff, thermDiff, stens, Pr)
end

function na_enthalpy(p, T)
    """This function takes the temperature of liquid sodium in K and returns the enthalpy in kJ/kg.
    The reference state is taken as solid sodium at 298 K. The correlation takes this form:
    H(l, T) - H(s, 298.15) = 365.77 + 1.6582 T - 4.2395 × 10^-4 * T^2 + 1.4847 × 10^-7 * T^3 + 2992.6 * T^-1
    This correlation is recommended by a report prepared by Argonne National Lab that can be found here:
    http://www.ne.anl.gov/eda/ANL-RE-95-2.pdf
    The correlation is valid between 371 and 2000 K."""
    a = -365.77
    b = 1.6582
    c = 4.2395E-4
    d = 1.4847E-7
    e = 2992.6
    return (a + b * T - c * T ^ 2 + d * T ^ 3 + e / T) * 1000
end

function na_Cp(p, T)
    """This function takes the temperature of liquid sodium and returns the heat
    capacity at constant pressure in kJ/kg*K. The Argonne report gives a very
    complicated correlation for Cp but notes that it's almost equal to the
    ForwardDiff.derivative of enthalpy with temperature up to about 1800 K. So I used the ForwardDiff.derivative here."""
    a = -365.77
    b = 1.6582
    c = 4.2395E-4
    d = 1.4847E-7
    e = 2992.6
    return (b - 2*c * T + 3 * d * T ^ 2 - e / (T^2)) * 1000
    #return ForwardDiff.derivative(enthalpy, T, dx = 1E-6)
end

function na_density(p, T)
    """This function takes the temperature of liquid sodium in K and returns the density in kg/m^3.
    The correlation takes this form:
    rho = rho_c + f * (1 - Tr) + g * (1 - Tr)^h
    where rho_c is the critical density of sodium, Tr is the reduced temperature T/T_c,
    and f g and h are constants given in the function.
    This correlation is recommended by a report prepared by Argonne National Lab that can be found here:
    http://www.ne.anl.gov/eda/ANL-RE-95-2.pdf
    The correlation is valid between 371 and ~2000 K."""
    rho_c = 219
    f = 275.32
    g = 511.58
    h = 0.5
    Tc = 2503.7
    Tr = T / Tc
    return rho_c + f * (1 - Tr) + g * (1 - Tr) ^ h
end

function na_thermal_conductivity(p, T)
    """This function takes the temperature of liquid sodium in K and returns the thermal conductivity in W/m*K.
    The reference state is taken as solid sodium at 298 K. The correlation takes this form:
    k = 124.67 - 0.11381 * T + 5.5226 × 10^-5 * T^2 - 1.1842 × 10^-8 * T^3
    This correlation is recommended by a report prepared by Argonne National Lab that can be found here:
    http://www.ne.anl.gov/eda/ANL-RE-95-2.pdf
    The correlation is valid between 371 and 2000 K."""
    a = 124.67
    b = 0.11381
    c = 5.5226E-5
    d = 1.1842E-8
    return a - b * T + c * T ^ 2 - d * T ^ 3
end

function na_dynamic_viscosity(p, T)
    """This function takes the temperature of liquid sodium in K and returns the dynamic viscosity in Pa * s.
    The reference state is taken as solid sodium at 298 K. The correlation takes this form:
    ln eta = - 6.4406 - 0.3958 ln (T) + 556.835 / T
    This correlation is recommended by a report prepared by Argonne National Lab that can be found here:
    http://www.ne.anl.gov/eda/ANL-RE-95-2.pdf
    The correlation is valid between 371 and 2500 K."""
    a = -6.4406
    b = -0.3958
    c = 556.835
    return exp(a + b * log(T) + c / T)
end


text:
"""
    checkresult(result::PSResult;
        logger::Union{Nothing, SimpleLogger} = nothing)->(Union{Nothing,T} where T<:Any)

Checks the value of the argument `result` which is of type PSResult for a value or an error.
Returns the value stored in result, or throws the error if the value stored in `result` is an Exception type.
If the optional argument `logger` is not Nothing, the error message is written to a log file.
"""
function checkresult(result::PSResult;
    logger::Union{Nothing,SimpleLogger}=nothing)::(Union{Nothing,T} where T <: Any)

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

    # Error case -
    if (isa(result.value, Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object)

        # log -
        if (isnothing(logger) == false)
            with_logger(logger) do
                @error(error_message)
            end
        end

        # throw -
        rethrow(result.value)
    end

    # default -
    return result.value
end

"""
    checkresult(result::PSResult)->(Union{Nothing,T} where T<:Any)

Checks the value of the argument `result` which is of type PSResult for a value or an error.
Returns the value stored in result, or rethrows the error if the value stored in `result` is an Exception type.
"""
function checkresult(result::PSResult)::(Union{Nothing,T} where T <: Any)

    # ok, so check, do we have an error object?
    # Yes: rethrow the error.
    # No: return the result.value

    # Error case -
    if (isa(result.value, Exception) == true)

        # throw -
        rethrow(result.value)
    end

    # default -
    return result.value
end


text:
using SACBP
using POMDPs
using MCTS
import Base.isequal

# Localiztion Task

# include("gradient_greedy.jl")

struct AugMDP <: POMDPs.MDP{AugState{<:PhysPos,<:VecBelMvNormal{T} where T<:Real},Vector{Float64}}
    dtc::Float64
    dto::Float64
    u_param_min::Vector{Float64}
    u_param_max::Vector{Float64}
    numActions::Int64
    Cu::Matrix{Float64}
    Q::Matrix{Float64}
end
function AugMDP(simPosUKF::SimulatePosRangeLocalization2D,
                u_param_min::Vector{Float64}, u_param_max::Vector{Float64},
                numActions::Int64)
    return AugMDP(simPosUKF.dtc, simPosUKF.dto,
                  u_param_min, u_param_max, numActions,
                  simPosUKF.Cu, simPosUKF.Q)
end

function POMDPs.actions(mdp::AugMDP)
    u_x = range(mdp.u_param_min[1],stop=mdp.u_param_max[1],length=mdp.numActions);
    u_y = range(mdp.u_param_min[2],stop=mdp.u_param_max[2],length=mdp.numActions);
    return [[u_x[ii],u_y[jj]] for ii = 1:length(u_x) for jj = 1:length(u_y)];
end
#=
struct Greedy <: POMDPs.Policy
    u_param_min::Vector{Float64}
    u_param_max::Vector{Float64}
    Q::Matrix{Float64}
    dto::Float64
end

function POMDPs.action(p::Greedy,s::AugState{PhysPos{Float64},VecBelMvNormal{Float64}})
    return gradientControlVal(UKFPosRange(),s,p.u_param_min,p.u_param_max,p.Q,p.dto)
end
=#

function POMDPs.gen(::DDNOut{(:sp, :r)},
                    mdp::AugMDP,
                    s::AugState{<:PhysPos,<:VecBelMvNormal{T} where T<:Real},
                    u::Vector{Float64},
                    rng::AbstractRNG)
    pathCost = 0.;
    SArray = Vector{typeof(s)}(undef,Int64(mdp.dto/mdp.dtc)+1);
    SArray[1] = s;
    for ii = 1:length(SArray)-1
        if SArray[ii].t < round(s.t + mdp.dto,digits=5)
            u_current = PosControl(round(s.t + mdp.dtc*(ii-1),digits=5),u);
            p = trans(TransModel_Pos(),SArray[ii].p,u_current,mdp.dtc);
            bVec = similar(SArray[ii].b);
            for jj = 1:length(bVec)
                bVec[jj] = BelMvNormal(p.t,SArray[ii].b[jj].params);
            end
            SArray[ii+1] = AugState(p,bVec);
            pathCost += instCost(CostPosRangeLocalization(),SArray[ii+1],u_current,mdp.Cu)*mdp.dtc;
        else
            u_current = PosControl(round(s.t + mdp.dtc*(ii-1),5),zeros(2));
            p = trans(TransModel_Pos(),SArray[ii].p,u_current,mdp.dtc);
            bVec = similar(SArray[ii].b);
            for jj = 1:length(bVec)
                bVec[jj] = BelMvNormal(p.t,SArray[ii].b[jj].params);
            end
            SArray[ii+1] = AugState(p,bVec);
            pathCost += instCost(CostPosRangeLocalization(),SArray[ii+1],u_current,mdp.Cu)*mdp.dtc;
        end
    end
    yVec,SArray[end] = trans(UKFPosRange(),SArray[end],mdp.dto,mdp.Q,rng)
    termCostDifference = termCost(CostPosRangeLocalization(),SArray[end]) - termCost(CostPosRangeLocalization(),s); # This is equivalent to having the terminal cost.
    reward = -(pathCost + termCostDifference);
    return SArray[end],reward
end;

function isequal(s1::AugState{SACBP.PhysPos{Float64},Array{SACBP.BelMvNormal{Float64},1}},
                 s2::AugState{SACBP.PhysPos{Float64},Array{SACBP.BelMvNormal{Float64},1}})
    return s1.t == s2.t && s1.p.pos == s2.p.pos && all([s1.b[ii].params == s2.b[ii].params for ii = 1:length(s1.b)])
end

POMDPs.discount(::AugMDP) = 1.0;

function mctsControlUpdate(policy::MCTS.DPWPlanner,
                           s::AugState{PhysPos{Float64},VecBelMvNormal{Float64}},
                           UArray::Vector{<:PosControl}, dto::Float64)
    u_val = action(policy,s)
    UArray_new = copy(UArray)
    for ii = 1:length(UArray_new)
        if UArray_new[ii].t < round(s.t + dto,digits=5)
            UArray_new[ii] = PosControl(UArray_new[ii].t,u_val);
        end
    end
    return UArray_new
end

struct MCTSParams
    numActions::Int64
    # treeDepth::Int64
    numSamples::Int64
    expConst::Float64
    kAct::Float64
    αAct::Float64
    kState::Float64
    αState::Float64
end

# Manipulation Task

include("p_control.jl")

struct BelMDP <: POMDPs.MDP{BelMvNormal{T} where T <: Real, Vector{Float64}}
    dtc::Float64
    dto::Float64
    u_param_min::Vector{Float64}
    u_param_max::Vector{Float64}
    numActions::Int64
    Q::Matrix{Float64}
    R::Matrix{Float64}
    Cs::Matrix{Float64}
    Cu::Matrix{Float64}
end
function BelMDP(simModel::SimulateManipulate2D,
                u_param_min::Vector{Float64}, u_param_max::Vector{Float64},
                numActions::Int64)
    return BelMDP(simModel.dtc, simModel.dto,
                  u_param_min, u_param_max, numActions,
                  simModel.Q, simModel.R, simModel.Cs, simModel.Cu)
end

function POMDPs.actions(mdp::BelMDP)
    u_fx = range(mdp.u_param_min[1],stop=mdp.u_param_max[1],length=mdp.numActions);
    u_fy = range(mdp.u_param_min[2],stop=mdp.u_param_max[2],length=mdp.numActions);
    u_tr = range(mdp.u_param_min[3],stop=mdp.u_param_max[3],length=mdp.numActions);
    return [[u_fx[ii],u_fy[jj],u_tr[kk]] for ii = 1:length(u_fx) for jj = 1:length(u_fy) for kk = 1:length(u_tr)];
end

struct DEKFManipulate2D <: GaussianFilter
# Discrete EKF model for belief transition in planar manipulation problem.
end

function trans_mcts(model::DEKFManipulate2D, b::BelMvNormal{T}, u::MControl2D{<:Real},
                    dt::Real, Q::Matrix{Float64}, R::Matrix{Float64}, rng::AbstractRNG) where T <: Real
    # EKF Predict and then Update for MCTS Transition
    if b.t != u.t
        error(ArgumentError("BelState and Control have inconsistent time parameters."));
    else
        transModel = TransModel_Manipulate2D();
        μ,Σ = b.params[1:b.dim],reshape(b.params[b.dim+1:end],b.dim,b.dim);

        # Predict.
        μMP = PhysManipulate2D(b.t,μ);
        A = trans_jacobi(transModel,μMP,u)*dt + Matrix(1.0I, b.dim, b.dim); # euler approximation
        Γ = A*Σ*A' + Q;
        Γ = Symmetric(Γ)
        try
            W = cholesky(Γ).L;
        catch
            println(μ)
            println(u)
            println(Σ)
        end
        W = cholesky(Γ).L;
        μ = vec(trans(transModel,PhysManipulate2D(b.t,μ),u,dt));

        # Update. (square-root form)
        μMP = PhysManipulate2D(b.t,μ);
        observeModel = ObserveModel_Manipulate2D();
        C = observe_jacobi(observeModel,μMP,u);
        Z = W'*C';
        H = Z'*Z + R
        H = Symmetric(H);
        try
            U = cholesky(H).L;
        catch
            println(μ)
            println(u)
            println(W)
        end
        U = cholesky(H).L;
        V = cholesky(R).L;
        K = Γ*C'/H;
        y = observe(observeModel,μMP,u,Matrix(H),rng);
        y_pred_mean = observe(observeModel,μMP,u);
        if y[3] - y_pred_mean[3] > pi   # Correct angular difference.
            y[3] -= 2*pi;
        elseif y[3] - y_pred_mean[3] < -pi
            y[3] += 2*pi;
        end
        μ += K*(y - y_pred_mean);
        μ[3] = mod2pi(μ[3]);
        Ω = W*(Matrix(1.0I, b.dim, b.dim) - Z/(U')/(U+V)*Z'); # See Andrews 1968
        g = [μ;vec(Symmetric(Ω*Ω'))];
        return BelMvNormal(round(b.t+dt,digits=5),g)
    end
end

function POMDPs.gen(::DDNOut{(:sp, :r)},
                    mdp::BelMDP,
                    b::BelMvNormal{<:Real},
                    u::Vector{Float64},
                    rng::AbstractRNG)
    b_new = trans_mcts(DEKFManipulate2D(),b,MControl2D(b.t,u),mdp.dto,mdp.Q,mdp.R,rng)
    reward = -instCost(CostManipulate2D(),b,MControl2D(b.t,u),mdp.Cs,mdp.Cu)*mdp.dto;
    return b_new,reward
end;

function isequal(b1::BelMvNormal{Float64},
                 b2::BelMvNormal{Float64})
    return b1.t == b2.t && b1.params == b2.params
end;

POMDPs.discount(::BelMDP) = 1.0;

function mctsControlUpdate(policy::MCTS.DPWPlanner,
                           b::BelMvNormal{Float64},
                           UArray::Vector{<:MControl2D},
                           dt::Float64)
    u_val = action(policy,b)
    UArray_new = copy(UArray)
    for ii = 1:length(UArray_new)
        if UArray_new[ii].t < round(b.t + dt,digits=5)
            UArray_new[ii] = MControl2D(UArray_new[ii].t,u_val);
        end
    end
    return UArray_new
end

struct PControl <: POMDPs.Policy
    u_param_min::Vector{Float64}
    u_param_max::Vector{Float64}
    posGain::Real
    rotGain::Real
end;

function POMDPs.action(p::PControl,b::BelMvNormal{Float64})
    return pControlVal(b,p.u_param_min,p.u_param_max,p.posGain,p.rotGain)
end


text:
###################################################
#
# Copyright © Akamai Technologies. All rights reserved.
#
# File: Tenant.jl
#
# Functions to communicate with the mPulse Repository REST API regarding Tenant Objects.
# This file MUST be `include()`d from `mPulseAPI.jl`
#
###################################################

export
    getRepositoryTenant

"""
Fetches a Tenant object from the mPulse repository

At least one of `tenantID` or `name` must be passed in to identify the tenant.

The tenant will be cached in memory for 1 hour, so subsequent calls using a matching `tenantID`, or `name` return
quickly without calling out to the API.  This can be a problem if the tenant changes in the repository.
You can clear the cache for this tenant using [`mPulseAPI.clearTenantCache`](@ref) and passing in one of `tenantID` or `name`.

#### Arguments
`token::AbstractString`
:    The Repository authentication token fetched by calling [`getRepositoryToken`](@ref)

#### Keyword Arguments
`tenantID::Int64`
:    The ID of the tenant to fetch.  This is the fastest method, but it can be hard to figure out a tenant's ID

`name::AbstractString`
:    The Tenant name in mPulse.  This is available from the mPulse tenant list.

#### Returns
`{Dict}` The `tenant` object with the following fields:

`name::AbstractString`
:    The tenant's name

`id::Int64`
:    The tenant's ID

`body::XMLElement`
:    An XML object representing the tenant's XML definition or an empty node if you do not have permission to see the full tenant

`parentID::Int64`
:    The ID of the parent folder that this tenant is in

`parentType::AbstractString`
:    The type of parent object (typically `tenantFolder`)

`path::AbstractString`
:    The folder path that this tenant is in

`description::AbstractString`
:    The description of this tenant entered into mPulse

`created::DateTime`
:    The timestamp when this object was created

`lastModified::DateTime`
:    The timestamp when this object was created

`attributes::Dict`
:    A `Dict` of attributes for this tenant

`dswbUrls::Array{AbstractString}`
:    An array of DSWB URLs that are valid auth redirect targets for this tenant

#### Throws
`ArgumentError`
:    if token is empty or tenantID and name are both empty

`mPulseAPIException`
:    if API access failed for some reason

"""
function getRepositoryTenant(token::AbstractString; tenantID::Int64=0, name::AbstractString="")
    tenant = getRepositoryObject(
                token,
                "tenant",
                Dict{Symbol, Any}(:id => tenantID, :name => name)
        )

    # If the object came out of cache, then it already contains these fields
    if !haskey(tenant, "dswbUrls")
        try
            tenant["dswbUrls"] = filter(u -> u != "", split(getNodeContent(tenant, "DSWBURLs", ""), ','))
        catch ex
            # If this is an Exception that we are not prepared to deal with
            if !isa(ex, LightXML.XMLParseError)
                rethrow()
            end
        end
    end

    tenant["dswb_dsn_name"] = "tenant_$(tenant["id"])"

    return tenant
end


text:
module PooksoftOptionsKit

# include -
include("Include.jl")

# export methods -
export compute_profit_loss_value
export compute_breakeven_price
export compute_option_profit_and_loss_at_expiration
export build_simulation_lattice_data_structure
export build_binary_lattice_data_structure
export build_ternary_lattice_data_structure
export build_simulation_contract_set
export build_simulation_price_array
export compute_underlying_price_distribution

# methods for pricing options (uses binary tree)
export option_contract_price
export longstaff_option_contract_price

# implied volatility -
export estimate_implied_volatility
export compute_weighted_volatility

# the greeks -
export delta
export theta
export gamma
export vega
export rho

# export types -
export PSCallOptionContract
export PSPutOptionContract
export PSEquityAsset

# Parameters for greek calculation -
export PSOptionKitPricingParameters
export PSBinaryLatticeModel
export PSTernaryLatticeModel

end # module


text:
using Zygote, IndexFunArrays
f(of,s) = sum(xx((10,),offset=of, scale=s))  # Float64(idx[1]+idx[2])*
f(2.2,(3.3,))
gradient(f, (1,), (2,))

kdsjfls
using View5D, IndexFunArrays

sz = (50,50)
N = 1
offsets = 1 .+ (sz.-1) .* rand(2,N)
# sigmas = 2.0 .*(0.3 .+rand(2,N))
mygaussians(of) = gaussian(sz,offset = of, weight=1.0, sigma=2.0);   # , sigma=sigmas
data =  mygaussians(offsets);
@vv data  # use "e" to advance spectral channels, "n" in the bottom right (element) panel to toggle normalization
@ve  mygaussians(offsets)

loss(of) = sum(abs2.(data .- mygaussians(of)))  # a loss function
loss(offsets)  # is zero
gradient(loss,(41,37))

using Zygote, IndexFunArrays

# f(idx,of,sc) = (idx[1].*sc - of[1] + of[2])^2
function c(of,sc)
    f(idx) = (idx[1].*sc - of[1] + of[2])^2
    IndexFunArray(f,(10,))
end
c((2.2,1.1),1.0)
loss(of,sc) = sum(c(of,sc))
gradient(loss,(2.2,1.3),1.2)

using Zygote, IndexFunArrays
f(of,s) = sum(xx((10,),offset=of, scale=s))  # Float64(idx[1]+idx[2])*
f(2.2,(3.3,))
gradient(f, 1, (2,))

f(of) = sum(xx((10,),offset=of))  # Float64(idx[1]+idx[2])*
f(2.2)
gradient(f, 1)

ntuple(i -> i ∈ dims ? scale[i] : zero(scale[1]), N)

# by "r" "g" "b" you can choose colors for elements. "v" toggles the color in and out of overlay,
# color overlay is accessed by "C" (shift-"c")


text:
using LoopFieldCalc

# Define grid on which we compute magnetic field strength.
# This is a 2D grid in the y-z plane
xs = 0.0:0.0
ys = -0.4:0.005:0.4
zs = -0.3:0.005:0.3

# Define number of wire turns per coil
nturns = 10

# Define geometry of inner coil
inner_coil_radius = 0.1
inner_coil_width = 0.15
inner_coil_offset = -0.1
inner_coil_current = 1.0

# Generate loops for inner coil
inner_coil = [
    LoopFieldCalc.CurrentLoop(
        radius = inner_coil_radius,
        current = inner_coil_current / nturns,
        center = LoopFieldCalc.CylindricalPoint(inner_coil_offset + z, 0.0, 0.0)
    )
    for z in LinRange(0.0, inner_coil_width, nturns)
]

# Define geometry for outer coil
outer_coil_radius = 0.25
outer_coil_width = 0.15
outer_coil_offset = -0.15
outer_coil_current = 0.5 * inner_coil_current

# Generate loops for outer coil
outer_coil = [
    LoopFieldCalc.CurrentLoop(
        radius = outer_coil_radius,
        current = outer_coil_current / nturns,
        center = LoopFieldCalc.CylindricalPoint(outer_coil_offset + z, 0.0, 0.0)
    )
    for z in LinRange(0.0, outer_coil_width, nturns)
]

# Combine inner and outer coil loops into a single vector and compute field
loops = [outer_coil; inner_coil]
Bx, By, Bz = LoopFieldCalc.field_on_grid(loops, xs, ys, zs)

# Write output. By reversing z and y and transposing the magnetic field matrices, we
# can rotate the output so that the loop axis is aligned with the x axis
LoopFieldCalc.write_field("examples/example_2.dat", zs, ys, Bz[1, :, :]', By[1, :, :]')


text:
# 2D arrays are confusing, wonder if I can do this better?
room = open("a11.txt") do file
  hcat([collect(line) for line in eachline(file)]...)
end

cartesian = CartesianIndices(room)
linear = LinearIndices(room)
Ifirst = first(cartesian)
Ilast = last(cartesian)
I1 = CartesianIndex(1,1)

seats = linear[findall((c) -> c == 'L', room)]

function simulate(affecting, tolerance)
  occupied = BitSet()
  while true
    change = BitSet()
    for seat in seats
      neighbours = length(intersect(affecting[seat], occupied))
      if in(seat, occupied)
        neighbours >= tolerance && push!(change, seat)
      else
        neighbours == 0 && push!(change, seat)
      end
    end
    length(change) == 0 && break
    symdiff!(occupied, change)
  end
  length(occupied)
end

function part1()
  adjacent = fill(BitSet(), axes(room))
  for I in cartesian
    if room[I] != 'L'
      continue
    end
    adjacent[I] = BitSet([linear[J] for J in max(Ifirst, I-I1):min(Ilast, I+I1) if I != J && room[J] == 'L'])
  end
  simulate(adjacent, 4)
end

using BenchmarkTools: @btime
ans1 = @btime part1()
println("$ans1")

function part2()
  visible = fill(BitSet(), axes(room))
  directions = [
    CartesianIndex(x,y) for x in -1:1 for y in -1:1 if (x,y) != (0,0)
  ]
  for I in cartesian
    if room[I] != 'L'
      continue
    end
    visibles = []
    for D in directions
      J = I
      while true
        J += D
        checkbounds(Bool, room, J) || break
        room[J] == 'L' || continue
        push!(visibles, linear[J])
        break
      end
    end
    visible[I] = BitSet(visibles)
  end
  simulate(visible, 5)
end

ans2 = @btime part2()
println("$ans2")


text:
for (mod, fname) in ((:Base, :sum), (:Base, :prod), (:Base, :maximum), (:Base, :minimum), (:Statistics, :mean))
    _fname = Symbol('_', fname)
    @eval begin
        @inline ($mod.$fname)(A::AbDimArray) = ($mod.$fname)(parent(A))
        @inline ($mod.$_fname)(A::AbstractArray, dims::AllDimensions) =
            rebuild(A, ($mod.$_fname)(parent(A), dimnum(A, dims)), reducedims(A, dims))
        @inline ($mod.$_fname)(f, A::AbstractArray, dims::AllDimensions) =
            rebuild(A, ($mod.$_fname)(f, parent(A), dimnum(A, dims)), reducedims(A, dims))
        @inline ($mod.$_fname)(A::AbDimArray, dims::Union{Int,Base.Dims}) =
            rebuild(A, ($mod.$_fname)(parent(A), dims), reducedims(A, dims))
        @inline ($mod.$_fname)(f, A::AbDimArray, dims::Union{Int,Base.Dims}) =
            rebuild(A, ($mod.$_fname)(f, parent(A), dims), reducedims(A, dims))
    end
end

for fname in (:std, :var)
    _fname = Symbol('_', fname)
    @eval begin
        @inline (Statistics.$fname)(A::AbDimArray) = (Statistics.$fname)(parent(A))
        @inline (Statistics.$_fname)(A::AbstractArray, corrected::Bool, mean, dims::AllDimensions) =
            rebuild(A, (Statistics.$_fname)(A, corrected, mean, dimnum(A, dims)), reducedims(A, dims))
        @inline (Statistics.$_fname)(A::AbDimArray, corrected::Bool, mean, dims::Union{Int,Base.Dims}) =
            rebuild(A, (Statistics.$_fname)(parent(A), corrected, mean, dims), reducedims(A, dims))
    end
end

Statistics.median(A::AbDimArray) = Statistics.median(parent(A))
Statistics._median(A::AbstractArray, dims::AllDimensions) =
    rebuild(A, Statistics._median(parent(A), dimnum(A, dims)), reducedims(A, dims))
Statistics._median(A::AbDimArray, dims::Union{Int,Base.Dims}) =
    rebuild(A, Statistics._median(parent(A), dims), reducedims(A, dims))

Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractArray, dims::AllDimensions) =
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbDimArray, dims::Union{Int,Base.Dims}) =
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbDimArray, dims::Colon) =
    Base._mapreduce_dim(f, op, nt, parent(A), dims)
# Unfortunately Base/accumulate.jl kwargs methods all force dims to be Integer.
# accumulate wont work unless that is relaxed, or we copy half of the file here.
Base._accumulate!(op, B, A, dims::AllDimensions, init::Union{Nothing, Some}) =
    Base._accumulate!(op, B, A, dimnum(A, dims), init)

Base._dropdims(A::AbstractArray, dim::Union{AbDim,Type{<:AbDim}}) =
    rebuildsliced(A, Base._dropdims(A, dimnum(A, dim)), dims2indices(A, basetypeof(dim)(1)))
Base._dropdims(A::AbstractArray, dims::AbDimTuple) =
    rebuildsliced(A, Base._dropdims(A, dimnum(A, dims)),
                  dims2indices(A, Tuple((basetypeof(d)(1) for d in dims))))

@inline Base.map(f, A::AbDimArray) = rebuild(A, map(f, parent(A)), dims(A))

# TODO cov, cor mapslices, eachslice, reverse, sort and sort! need _methods without kwargs in base so
# we can dispatch on dims. Instead we dispatch on array type for now, which means
# these aren't usefull unless you inherit from AbDimArray.

Base.mapslices(f, A::AbDimArray; dims=1, kwargs...) = begin
    dimnums = dimnum(A, dims)
    data = mapslices(f, parent(A); dims=dimnums, kwargs...)
    rebuild(A, data, reducedims(A, DimensionalData.dims(A, dimnums)))
end

# This is copied from base as we can't efficiently wrap this function
# through the kwarg with a rebuild in the generator. Doing it this way
# wierdly makes it faster to use a dim than an integer.
if VERSION > v"1.1-"
    Base.eachslice(A::AbDimArray; dims=1, kwargs...) = begin
        if dims isa Tuple && length(dims) == 1
            throw(ArgumentError("only single dimensions are supported"))
        end
        dim = first(dimnum(A, dims))
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
        return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
    end

end

for fname in (:cor, :cov)
    @eval Statistics.$fname(A::AbDimArray{T,2}; dims=1, kwargs...) where T = begin
        newdata = Statistics.$fname(parent(A); dims=dimnum(A, dims), kwargs...)
        I = dims2indices(A, dims, 1)
        newdims, newrefdims = slicedims(A, I)
        rebuild(A, newdata, (newdims[1], newdims[1]), newrefdims)
    end
end

@inline Base.reverse(A::AbDimArray{T,N}; dims=1) where {T,N} = begin
    dnum = dimnum(A, dims)
    # Reverse the dimension. TODO: make this type stable
    newdims = revdims(DimensionalData.dims(A), dnum)
    # Reverse the data
    newdata = reverse(parent(A); dims=dnum)
    rebuild(A, newdata, newdims, refdims(A))
end

# TODO change order after reverse
@inline revdims(dimstorev::Tuple, dnum) = begin
    dim = dimstorev[end]
    if length(dimstorev) == dnum
        dim = rebuild(dim, reverse(val(dim)))
    end
    (revdims(Base.front(dimstorev), dnum)..., dim)
end
@inline revdims(dims::Tuple{}, i) = ()

for (pkg, fname) in [(:Base, :permutedims), (:Base, :adjoint),
                     (:Base, :transpose), (:LinearAlgebra, :Transpose)]
    @eval begin
        @inline $pkg.$fname(A::AbDimArray{T,2}) where T =
            rebuild(A, $fname(parent(A)), reverse(dims(A)), refdims(A))
    end
end

for fname in [:permutedims, :PermutedDimsArray]
    @eval begin
        @inline Base.$fname(A::AbDimArray{T,N}, perm) where {T,N} =
            rebuild(A, $fname(parent(A), dimnum(A, perm)), permutedims(dims(A), perm))
    end
end


text:
"""
    song(n::Integer)

Repeat lyrics to song starting at `n` and counting down to 1.
Creadit: Rosetta Code
"""

function song(start::Integer)
    lines(i) = i == 0? "No more lines" : i == 1 ? "1 line" : "$i lines"
    for i in start:-1:1
        PRINTLN("""
                    $(lines(i)) of Code on the Screen,
                    $(lines(i)) of Code...
                    \tControl the Flow,
                    \tDelete one Row,
                    $(lines(i - 1)) of Code on the Screen!!!
                """)
    end
end

max_stars_repo_path: 99 Lines of Code/song.jl

max_issues_repo_path: 99 Lines of Code/song.jl

max_forks_repo_path: 99 Lines of Code/song.jl


text:
"""
    center_mode(X::AbstractArray, I::CartesianIndex, offsets)
    center_mode(X::IntegralArray, I::CartesianIndex, bo::CartesianIndex, offsets)

Compute the local binary pattern using center pixel as the block mode. This is the
version that original LBP is used.

See also [`average_mode`](@ref) for mLBP (modified LBP) version.
```
"""
Base.@propagate_inbounds center_mode(X, I, offsets) = X[I]
Base.@propagate_inbounds function center_mode(X::IntegralArray, I, bo, offsets)
    Rp = I:min(I+bo, last(CartesianIndices(X)))
    X[first(Rp)..last(Rp)]/length(Rp)
end

"""
    average_mode(X, I, offsets) -> value
    local_binary_pattern(average_mode, X, args...; kwargs...)

Compute the local binary pattern using the average mode of the block. This is also called
mLBP (modified LBP) in the literature [1].

Original local binary pattern compares the neighbors with the center value `X[I]`, this
modified version instead uses the mean value of the block.

See also [`center_mode`](@ref) for original LBP version.

# References

- [1] Wikipedia contributors. "Local binary patterns." _Wikipedia, The Free Encyclopedia_. Wikipedia, The Free Encyclopedia, 17 Nov. 2021. Web. 1 Jan. 2022.
"""
Base.@propagate_inbounds function average_mode(X, I::CartesianIndex, offsets)
    v = X[I]
    ax = axes(X)
    rst = v + mapreduce(+, offsets) do o
        p = I.I .+ o
        inbounds = map(in, p, ax)
        all(inbounds) ? _inbounds_getindex(X, p) : v
    end
    return rst / (length(offsets) + 1)
end

Base.@propagate_inbounds function average_mode(X::IntegralArray, I::CartesianIndex, bo::CartesianIndex, offsets)
    # multi-block version
    R = CartesianIndices(X)
    Rfirst, Rlast = first(R), last(R)
    Rp = I:min(I+bo, Rlast)
    v = X[I..last(Rp)]/length(Rp)
    rst = v + mapreduce(+, offsets) do o
        p = CartesianIndex(I.I .+ o)
        Rp = max(p, Rfirst):min(p+bo, Rlast)
        !isempty(Rp) ? X[first(Rp)..last(Rp)]/length(Rp) : v
    end
    return rst / (length(offsets) + 1)
end


text:
using LinearAlgebraicRepresentation
using Base.Test

function BoxCalculation(Vertices)
	Minx=minimum(Vertices[1,:])
	Maxx=maximum(Vertices[1,:])
	Miny=minimum(Vertices[2,:])
	Maxy=maximum(Vertices[2,:])
	dx=Maxx-Minx
	dy=Maxy-Miny
	Box=dx*dy
	if size(Vertices,1)==3
		Minz=minimum(Vertices[3,:])
		Maxz=maximum(Vertices[3,:])
		dz=Maxz-Minz
		Box=Box*dz
	end
	return Box
end

@testset "circle" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.circle()()[1])==4
	@test BoxCalculation(LinearAlgebraicRepresentation.circle(2., 2*pi)()[1])==16
	@test BoxCalculation(LinearAlgebraicRepresentation.circle(3, 2*pi)()[1])==36
	@test BoxCalculation(LinearAlgebraicRepresentation.circle(5, pi/2)()[1])==25
	@test size(LinearAlgebraicRepresentation.circle(3,2*pi)(60)[1],2)==60
	@test length(LinearAlgebraicRepresentation.circle(3,2*pi)(60)[2])==60
end

@testset "helix" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.helix()()[1])==8
	@test BoxCalculation(LinearAlgebraicRepresentation.helix(2, 2, 2)()[1])==64
	@test BoxCalculation(LinearAlgebraicRepresentation.helix(1, 2, 5)()[1])==40
	@test BoxCalculation(LinearAlgebraicRepresentation.helix(1, 1, 3)()[1])==12
	#=
	1, 1, 2, is contained
	in a box that has volume (radius*2)^2*pitch*nturns
	=#
	@test size(LinearAlgebraicRepresentation.helix(5, 7, 9)()[1],2)==325
	@test length(LinearAlgebraicRepresentation.helix(5, 7, 9)()[2])==324
end

@testset "disk" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.disk(1., 2*pi)([36, 1])[1])==4
	@test BoxCalculation(LinearAlgebraicRepresentation.disk(2., 2*pi)()[1])==16
	@test BoxCalculation(LinearAlgebraicRepresentation.disk(2., pi)()[1])==8
	@test BoxCalculation(LinearAlgebraicRepresentation.disk(2., pi/2)()[1])==4
	@test size(LinearAlgebraicRepresentation.disk(10, pi/7)()[1],2)==75
	@test length(LinearAlgebraicRepresentation.disk(10, pi/7)()[2])==108
end

@testset "helicoid" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.helicoid()()[1])==8
	@test BoxCalculation(LinearAlgebraicRepresentation.helicoid(2,1,2,2)()[1])==64
	@test BoxCalculation(LinearAlgebraicRepresentation.helicoid(1,.5,2,5)()[1])==40
	@test BoxCalculation(LinearAlgebraicRepresentation.helicoid(1,.3,1,3)()[1])==12
	#=
	1, 1, 2, is contained in a
	box that has volume (radius*2)^2*pitch*nturns
	=#
	@test size(LinearAlgebraicRepresentation.helicoid(1.3,.7,1,3)()[1],2)==327
	@test length(LinearAlgebraicRepresentation.helicoid(1.3,.7,1,3)()[2])==432
end

@testset "ring" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.ring(1.,3.,2*pi)()[1])==36
	@test BoxCalculation(LinearAlgebraicRepresentation.ring(1.,2,pi)()[1])==8
	@test BoxCalculation(LinearAlgebraicRepresentation.ring(1.,5,pi/2)()[1])==25
	#(radius*2)^2
	@test size(LinearAlgebraicRepresentation.ring(1,5,pi/2)()[1],2)==74
	@test length(LinearAlgebraicRepresentation.ring(5,10,pi/6)()[2])==36
end

@testset "cylinder" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.cylinder(1,5,2*pi)()[1])==20
	@test BoxCalculation(LinearAlgebraicRepresentation.cylinder(2,2,pi)()[1])==16
	@test BoxCalculation(LinearAlgebraicRepresentation.cylinder(1,4,pi)()[1])==8
	#((radius*2)^2)*height
	@test size(LinearAlgebraicRepresentation.cylinder(3.4,20,pi/7)()[1],2)==74
	@test length(LinearAlgebraicRepresentation.cylinder(3.4,20,pi/7)()[2])==36
end

@testset "sphere" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.sphere(2,pi,2*pi)()[1])==64
	@test BoxCalculation(LinearAlgebraicRepresentation.sphere(6,pi,pi)()[1])==864
	@test BoxCalculation(LinearAlgebraicRepresentation.sphere(4,pi,2*pi)()[1])==8^3
	@test size(LinearAlgebraicRepresentation.sphere(2.5,pi/3,pi/5)()[1],2)==703
	@test length(LinearAlgebraicRepresentation.sphere(2.5,pi/3,pi/5)()[2])==1296
end

@testset "toroidal" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.toroidal(1,3,2*pi,2*pi)()[1])==128
	@test BoxCalculation(LinearAlgebraicRepresentation.toroidal(2,3,2*pi,2*pi)()[1])==400
	#(((R+r)*2)^2)*(r*2)
	@test size(LinearAlgebraicRepresentation.toroidal(1.3,4.6,pi/4,pi/7)()[1],2)==925
	@test length(LinearAlgebraicRepresentation.toroidal(1.3,4.6,pi/4,pi/7)()[2])==1728
end

@testset "crown" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.crown(1., 3., 2*pi)()[1])==128
	@test BoxCalculation(LinearAlgebraicRepresentation.crown(2., 3., 2*pi)()[1])==400
	#(((R+r)*2)^2)*(r*2)
	@test size(LinearAlgebraicRepresentation.crown(1.5, 5.6, pi/8)()[1],2)==481
	@test length(LinearAlgebraicRepresentation.crown(1.5, 5.6, pi/8)()[2])==864
end

@testset "ball" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.ball(1, pi, 2*pi)()[1])==8
	@test BoxCalculation(LinearAlgebraicRepresentation.ball(6, pi, pi)()[1])==864
	@test BoxCalculation(LinearAlgebraicRepresentation.ball(3, pi, 2*pi)()[1])==6^3
	@test size(LinearAlgebraicRepresentation.ball(2.6, pi/5, pi/9)()[1],2)==2813
	@test length(LinearAlgebraicRepresentation.ball(2.6, pi/5, pi/9)()[2])==2592
end

@testset "rod" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.rod(1,5,2*pi)()[1])==20
	@test BoxCalculation(LinearAlgebraicRepresentation.rod(2,2,pi)()[1])==16
	@test BoxCalculation(LinearAlgebraicRepresentation.rod(1,4,pi)()[1])==8
	#((radius*2)^2)*height
	@test size(LinearAlgebraicRepresentation.rod(3.7,8.9,pi/9)()[1],2)==74
	@test length(LinearAlgebraicRepresentation.rod(3.7,8.9,pi/9)()[2])==1
end

@testset "hollowCyl" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.hollowCyl(0,1.,5,2*pi)()[1])==20
	@test BoxCalculation(LinearAlgebraicRepresentation.hollowCyl(1,2.,4,pi)()[1])==32
	@test BoxCalculation(LinearAlgebraicRepresentation.hollowCyl(1,4,3,pi/2)()[1])==48
	#((radius*2)^2)*height
	@test size(LinearAlgebraicRepresentation.hollowCyl(3,4.,7.8,pi/5)()[1],2)==148
	@test length(LinearAlgebraicRepresentation.hollowCyl(3,4.,7.8,pi/5)()[2])==36
end

@testset "hollowBall" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.hollowBall(1,2,pi,2*pi)([36,36,1])[1]) ==64
	@test BoxCalculation(LinearAlgebraicRepresentation.hollowBall(1,6,pi,pi)([36,36,1])[1])== 864
	@test BoxCalculation(LinearAlgebraicRepresentation.hollowBall(2,4,pi,2*pi)([36,36,1])[1])==8^3
	#(radius*2)^3
	@test size(LinearAlgebraicRepresentation.hollowBall(1.5,6.7,pi/3,2*pi/3)()[1],2)==3700
	@test length(LinearAlgebraicRepresentation.hollowBall(1.5,6.7,pi/3,2*pi/3)()[2])==2592
end

@testset "torus" begin
	@test BoxCalculation(LinearAlgebraicRepresentation.torus(1.,2.,.5,2*pi,2*pi)()[1])==147
	@test BoxCalculation(LinearAlgebraicRepresentation.torus(2,3,.5,2*pi,2*pi)()[1])==605.0
	#(((R+r)*2)^2)*(r*2)
	@test size(LinearAlgebraicRepresentation.torus(5.2,7,.5,pi/3,pi/4)()[1],2)==4625
	@test length(LinearAlgebraicRepresentation.torus(5.2,7,.5,pi/3,pi/4)()[2])==3456
end