text:
module OffsetArrays

using Base: tail, @propagate_inbounds
@static if !isdefined(Base, :IdentityUnitRange)
    const IdentityUnitRange = Base.Slice
else
    using Base: IdentityUnitRange
end

export OffsetArray, OffsetMatrix, OffsetVector

const IIUR = IdentityUnitRange{<:AbstractUnitRange{<:Integer}}

include("axes.jl")
include("utils.jl")
include("origin.jl")

# Technically we know the length of CartesianIndices but we need to convert it first, so here we
# don't put it in OffsetAxisKnownLength.
const OffsetAxisKnownLength = Union{Integer, AbstractUnitRange}
const OffsetAxis = Union{OffsetAxisKnownLength, Colon}
const ArrayInitializer = Union{UndefInitializer, Missing, Nothing}

## OffsetArray
"""
    OffsetArray(A, indices...)

Return an `AbstractArray` that shares element type and size with the first argument but
uses the supplied `indices` to infer its axes. If all the indices are `AbstractUnitRange`s then
these are directly used as the axis span along each dimension. Refer to the examples below for other
permissible types.

Alternatively it's possible to specify the coordinates of one corner of the array
and have the axes be computed automatically from the size of `A`.
This constructor makes it convenient to shift to
an arbitrary starting index along each axis, for example to a zero-based indexing scheme followed by
arrays in languages such as C and Python.
See [`Origin`](@ref) and the examples below for this usage.

# Example: offsets

There are two types of `indices`: integers and ranges-like types.

Integers are recognized as offsets, where `0` means no offsets are applied:

```jldoctest; setup=:(using OffsetArrays)
julia> A = OffsetArray(reshape(1:6, 2, 3), -1, -2)
2×3 OffsetArray(reshape(::UnitRange{$Int}, 2, 3), 0:1, -1:1) with eltype $Int with indices 0:1×-1:1:
 1  3  5
 2  4  6

julia> A[0, 1]
5
```

Examples of range-like types are: `UnitRange` (e.g, `-1:2`), `CartesianIndices`,
and `Colon()` (or concisely `:`). A `UnitRange` specifies the axis span along one particular dimension,
`CartesianIndices` specify the axis spans along multiple dimensions, and a `Colon` is a placeholder
that specifies that the `OffsetArray` shares its axis with its parent along that dimension.

```jldoctest; setup=:(using OffsetArrays)
julia> OffsetArray(reshape(1:6, 2, 3), 0:1, -1:1)
2×3 OffsetArray(reshape(::UnitRange{$Int}, 2, 3), 0:1, -1:1) with eltype $Int with indices 0:1×-1:1:
 1  3  5
 2  4  6

julia> OffsetArray(reshape(1:6, 2, 3), :, -1:1) # : as a placeholder to indicate that no offset is to be applied to the first dimension
2×3 OffsetArray(reshape(::UnitRange{$Int}, 2, 3), 1:2, -1:1) with eltype $Int with indices 1:2×-1:1:
 1  3  5
 2  4  6
```

Use `CartesianIndices` to specify the coordinates of two diagonally opposite corners:

```jldoctest; setup=:(using OffsetArrays)
julia> OffsetArray(reshape(1:6, 2, 3), CartesianIndex(0, -1):CartesianIndex(1, 1))
2×3 OffsetArray(reshape(::UnitRange{$Int}, 2, 3), 0:1, -1:1) with eltype $Int with indices 0:1×-1:1:
 1  3  5
 2  4  6
```

Integers and range-like types may not be combined in the same call:

```julia
julia> OffsetArray(reshape(1:6, 2, 3), 0, -1:1)
ERROR: [...]
```

# Example: origin

[`OffsetArrays.Origin`](@ref) may be used to specify the origin of the OffsetArray. The term origin here
refers to the corner with the lowest values of coordinates, such as the left edge for an `AbstractVector`,
the bottom left corner for an `AbstractMatrix` and so on. The coordinates of the origin sets the starting
index of the array along each dimension.

```jldoctest; setup=:(using OffsetArrays)
julia> a = [1 2; 3 4];

julia> OffsetArray(a, OffsetArrays.Origin(0, 1))
2×2 OffsetArray(::$(Array{Int,2}), 0:1, 1:2) with eltype $Int with indices 0:1×1:2:
 1  2
 3  4

julia> OffsetArray(a, OffsetArrays.Origin(0)) # set the origin to zero along each dimension
2×2 OffsetArray(::$(Array{Int, 2}), 0:1, 0:1) with eltype $Int with indices 0:1×0:1:
 1  2
 3  4
```

"""
struct OffsetArray{T,N,AA<:AbstractArray{T,N}} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
    @inline function OffsetArray{T, N, AA}(parent::AA, offsets::NTuple{N, Int}; checkoverflow = true) where {T, N, AA<:AbstractArray{T,N}}
        # allocation of `map` on tuple is optimized away
        checkoverflow && map(overflow_check, axes(parent), offsets)
        new{T, N, AA}(parent, offsets)
    end
end

"""
    OffsetVector(v, index)

Type alias and convenience constructor for one-dimensional [`OffsetArray`](@ref)s.
"""
const OffsetVector{T,AA<:AbstractVector{T}} = OffsetArray{T,1,AA}

"""
    OffsetMatrix(A, index1, index2)

Type alias and convenience constructor for two-dimensional [`OffsetArray`](@ref)s.
"""
const OffsetMatrix{T,AA<:AbstractMatrix{T}} = OffsetArray{T,2,AA}

# checks if the offset may be added to the range without overflowing
function overflow_check(r::AbstractUnitRange, offset::Integer)
    Base.hastypemax(eltype(r)) || return nothing
    # This gives some performance boost https://github.com/JuliaLang/julia/issues/33273
    throw_upper_overflow_error(val) = throw(OverflowError("offset should be <= $(typemax(Int) - val) corresponding to the axis $r, received an offset $offset"))
    throw_lower_overflow_error(val) = throw(OverflowError("offset should be >= $(typemin(Int) - val) corresponding to the axis $r, received an offset $offset"))

    # With ranges in the picture, first(r) might not necessarily be < last(r)
    # we therefore use the min and max of first(r) and last(r) to check for overflow
    firstlast_min, firstlast_max = minmax(first(r), last(r))

    if offset > 0 && firstlast_max > typemax(Int) - offset
        throw_upper_overflow_error(firstlast_max)
    elseif offset < 0 && firstlast_min < typemin(Int) - offset
        throw_lower_overflow_error(firstlast_min)
    end
    return nothing
end

# Tuples of integers are treated as offsets
# Empty Tuples are handled here
@inline function OffsetArray(A::AbstractArray, offsets::Tuple{Vararg{Integer}}; kw...)
    _checkindices(A, offsets, "offsets")
    OffsetArray{eltype(A), ndims(A), typeof(A)}(A, offsets; kw...)
end

# These methods are necessary to disallow incompatible dimensions for
# the OffsetVector and the OffsetMatrix constructors
for (FT, ND) in ((:OffsetVector, :1), (:OffsetMatrix, :2))
    @eval @inline function $FT(A::AbstractArray{<:Any,$ND}, offsets::Tuple{Vararg{Integer}}; kw...)
        _checkindices(A, offsets, "offsets")
        OffsetArray{eltype(A), $ND, typeof(A)}(A, offsets; kw...)
    end
    FTstr = string(FT)
    @eval @inline function $FT(A::AbstractArray, offsets::Tuple{Vararg{Integer}}; kw...)
        throw(ArgumentError($FTstr*" requires a "*string($ND)*"D array"))
    end
end

## OffsetArray constructors
for FT in (:OffsetArray, :OffsetVector, :OffsetMatrix)
    # Nested OffsetArrays may strip off the wrapper and collate the offsets
    # empty tuples are handled here
    @eval @inline function $FT(A::OffsetArray, offsets::Tuple{Vararg{Int}}; checkoverflow = true)
        _checkindices(A, offsets, "offsets")
        # ensure that the offsets may be added together without an overflow
        checkoverflow && map(overflow_check, axes(A), offsets)
        I = map(+, _offsets(A, parent(A)), offsets)
        $FT(parent(A), I, checkoverflow = false)
    end
    @eval @inline function $FT(A::OffsetArray, offsets::Tuple{Integer,Vararg{Integer}}; kw...)
        $FT(A, map(Int, offsets); kw...)
    end

    # In general, indices get converted to AbstractUnitRanges.
    # CartesianIndices{N} get converted to N ranges
    @eval @inline function $FT(A::AbstractArray, inds::Tuple{Any,Vararg{Any}}; kw...)
        $FT(A, _toAbstractUnitRanges(to_indices(A, axes(A), inds)); kw...)
    end

    # convert ranges to offsets
    @eval @inline function $FT(A::AbstractArray, inds::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}; kw...)
        _checkindices(A, inds, "indices")
        # Performance gain by wrapping the error in a function: see https://github.com/JuliaLang/julia/issues/37558
        throw_dimerr(lA, lI) = throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
        lA = size(A)
        lI = map(length, inds)
        lA == lI || throw_dimerr(lA, lI)
        $FT(A, map(_offset, axes(A), inds); kw...)
    end

    @eval @inline $FT(A::AbstractArray, inds::Vararg; kw...) = $FT(A, inds; kw...)
    @eval @inline $FT(A::AbstractArray; checkoverflow = false) = $FT(A, ntuple(zero, Val(ndims(A))), checkoverflow = checkoverflow)

    @eval @inline $FT(A::AbstractArray, origin::Origin; checkoverflow = true) = $FT(A, origin(A); checkoverflow = checkoverflow)
end

# conversion-related methods
@inline OffsetArray{T}(M::AbstractArray, I...; kw...) where {T} = OffsetArray{T,ndims(M)}(M, I...; kw...)

@inline function OffsetArray{T,N}(M::AbstractArray{<:Any,N}, I...; kw...) where {T,N}
    M2 = _of_eltype(T, M)
    OffsetArray{T,N}(M2, I...; kw...)
end
@inline OffsetArray{T,N}(M::OffsetArray{T,N}, I...; kw...) where {T,N} = OffsetArray(M, I...; kw...)
@inline OffsetArray{T,N}(M::AbstractArray{T,N}, I...; kw...) where {T,N} = OffsetArray{T,N,typeof(M)}(M, I...; kw...)

@inline OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::Vararg; kw...) where {T,N,A<:AbstractArray{T,N}} = OffsetArray{T,N,A}(M, I; kw...)
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::NTuple{N,Int}; checkoverflow = true) where {T,N,A<:AbstractArray{T,N}}
    checkoverflow && map(overflow_check, axes(M), I)
    Mv = no_offset_view(M)
    MvA = convert(A, Mv)::A
    Iof = map(+, _offsets(M), I)
    OffsetArray{T,N,A}(MvA, Iof, checkoverflow = false)
end
@inline function OffsetArray{T, N, AA}(parent::AbstractArray{<:Any,N}, offsets::NTuple{N, Integer}; kw...) where {T, N, AA<:AbstractArray{T,N}}
    OffsetArray{T, N, AA}(parent, map(Int, offsets)::NTuple{N,Int}; kw...)
end
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}; kw...) where {T,N,A<:AbstractArray{T,N}}
    _checkindices(M, I, "indices")
    # Performance gain by wrapping the error in a function: see https://github.com/JuliaLang/julia/issues/37558
    throw_dimerr(lA, lI) = throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
    lM = size(M)
    lI = map(length, I)
    lM == lI || throw_dimerr(lM, lI)
    OffsetArray{T,N,A}(M, map(_offset, axes(M), I); kw...)
end
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::Tuple; kw...) where {T,N,A<:AbstractArray{T,N}}
    OffsetArray{T,N,A}(M, _toAbstractUnitRanges(to_indices(M, axes(M), I)); kw...)
end
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}; kw...) where {T,N,A<:AbstractArray{T,N}}
    Mv = no_offset_view(M)
    MvA = convert(A, Mv)::A
    OffsetArray{T,N,A}(MvA, _offsets(M); kw...)
end
@inline OffsetArray{T,N,A}(M::A; checkoverflow = false) where {T,N,A<:AbstractArray{T,N}} = OffsetArray{T,N,A}(M, ntuple(zero, Val(N)); checkoverflow = checkoverflow)

Base.convert(::Type{T}, M::AbstractArray) where {T<:OffsetArray} = M isa T ? M : T(M)

@inline AbstractArray{T,N}(M::OffsetArray{S,N}) where {T,S,N} = OffsetArray{T}(M)

# array initialization
@inline function OffsetArray{T,N}(init::ArrayInitializer, inds::Tuple{Vararg{OffsetAxisKnownLength}}; kw...) where {T,N}
    _checkindices(N, inds, "indices")
    AA = Array{T,N}(init, map(_indexlength, inds))
    OffsetArray{T, N, typeof(AA)}(AA, map(_indexoffset, inds); kw...)
end
@inline function OffsetArray{T, N}(init::ArrayInitializer, inds::Tuple; kw...) where {T, N}
    OffsetArray{T, N}(init, _toAbstractUnitRanges(inds); kw...)
end
@inline OffsetArray{T,N}(init::ArrayInitializer, inds::Vararg; kw...) where {T,N} = OffsetArray{T,N}(init, inds; kw...)

@inline OffsetArray{T}(init::ArrayInitializer, inds::NTuple{N, OffsetAxisKnownLength}; kw...) where {T,N} = OffsetArray{T,N}(init, inds; kw...)
@inline function OffsetArray{T}(init::ArrayInitializer, inds::Tuple; kw...) where {T}
    OffsetArray{T}(init, _toAbstractUnitRanges(inds); kw...)
end
@inline OffsetArray{T}(init::ArrayInitializer, inds::Vararg; kw...) where {T} = OffsetArray{T}(init, inds; kw...)

Base.IndexStyle(::Type{OA}) where {OA<:OffsetArray} = IndexStyle(parenttype(OA))
parenttype(::Type{OffsetArray{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::OffsetArray) = parenttype(typeof(A))

Base.parent(A::OffsetArray) = A.parent

# TODO: Ideally we would delegate to the parent's broadcasting implementation, but that
#       is currently broken in sufficiently many implementation, namely RecursiveArrayTools, DistributedArrays
#       and StaticArrays, that it will take concentrated effort to get this working across the ecosystem.
#       The goal would be to have `OffsetArray(CuArray) .+ 1 == OffsetArray{CuArray}`.
# Base.Broadcast.BroadcastStyle(::Type{<:OffsetArray{<:Any, <:Any, AA}}) where AA = Base.Broadcast.BroadcastStyle(AA)

@inline Base.size(A::OffsetArray) = size(parent(A))

@inline Base.axes(A::OffsetArray) = map(IdOffsetRange, axes(parent(A)), A.offsets)
@inline Base.axes(A::OffsetArray, d) = d <= ndims(A) ? IdOffsetRange(axes(parent(A), d), A.offsets[d]) : IdOffsetRange(axes(parent(A), d))
@inline Base.axes1(A::OffsetArray{T,0}) where {T} = IdOffsetRange(axes(parent(A), 1))  # we only need to specialize this one

# Issue 128
# See https://github.com/JuliaLang/julia/issues/37274 for the issue reported in Base
# The fix https://github.com/JuliaLang/julia/pull/39404 should be available on v1.6
# The following method is added on older Julia versions to ensure correct behavior for OffsetVectors
if VERSION < v"1.6"
    @inline function Base.compute_linindex(A::OffsetVector, I::NTuple{N,Any}) where N
        IP = Base.fill_to_length(axes(A), Base.OneTo(1), Val(N))
        Base.compute_linindex(first(LinearIndices(A)), 1, IP, I)
    end
end

# Utils to translate a function to the parent while preserving offsets
unwrap(x) = x, identity
unwrap(x::OffsetArray) = parent(x), data -> OffsetArray(data, x.offsets, checkoverflow = false)
function parent_call(f, x)
    parent, wrap_offset = unwrap(x)
    wrap_offset(f(parent))
end

Base.similar(A::OffsetArray, ::Type{T}, dims::Dims) where T =
    similar(parent(A), T, dims)
function Base.similar(A::AbstractArray, ::Type{T}, shape::Tuple{OffsetAxisKnownLength,Vararg{OffsetAxisKnownLength}}) where T
    # strip IdOffsetRanges to extract the parent range and use it to generate the array
    new_shape = map(_strip_IdOffsetRange, shape)
    # route through _similar_axes_or_length to avoid a stack overflow if map(_strip_IdOffsetRange, shape) === shape
    # This tries to use new_shape directly in similar if similar(A, T, ::typeof(new_shape)) is defined
    # If this fails, it calls similar(A, T, map(_indexlength, new_shape)) to use the size along each axis
    # to generate the new array
    P = _similar_axes_or_length(A, T, new_shape, shape)
    return OffsetArray(P, map(_offset, axes(P), shape))
end
function Base.similar(::Type{T}, shape::Tuple{OffsetAxisKnownLength,Vararg{OffsetAxisKnownLength}}) where {T<:AbstractArray}
    new_shape = map(_strip_IdOffsetRange, shape)
    P = _similar_axes_or_length(T, new_shape, shape)
    OffsetArray(P, map(_offset, axes(P), shape))
end
# Try to use the axes to generate the parent array type
# This is useful if the axes have special meanings, such as with static arrays
# This method is hit if at least one axis provided to similar(A, T, axes) is an IdOffsetRange
# For example this is hit when similar(A::OffsetArray) is called,
# which expands to similar(A, eltype(A), axes(A))
_similar_axes_or_length(A, T, ax, ::Any) = similar(A, T, ax)
_similar_axes_or_length(AT, ax, ::Any) = similar(AT, ax)
# Handle the general case by resorting to lengths along each axis
# This is hit if none of the axes provided to similar(A, T, axes) are IdOffsetRanges,
# and if similar(A, T, axes::AX) is not defined for the type AX.
# In this case the best that we can do is to create a mutable array of the correct size
_similar_axes_or_length(A, T, ax::I, ::I) where {I} = similar(A, T, map(_indexlength, ax))
_similar_axes_or_length(AT, ax::I, ::I) where {I} = similar(AT, map(_indexlength, ax))

# reshape accepts a single colon
Base.reshape(A::AbstractArray, inds::OffsetAxis...) = reshape(A, inds)
function Base.reshape(A::AbstractArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}})
    AR = reshape(no_offset_view(A), map(_indexlength, inds))
    O = OffsetArray(AR, map(_offset, axes(AR), inds))
    return _popreshape(O, axes(AR), _filterreshapeinds(inds))
end

# Reshaping OffsetArrays can "pop" the original OffsetArray wrapper and return
# an OffsetArray(reshape(...)) instead of an OffsetArray(reshape(OffsetArray(...)))
# Short-circuit for AbstractVectors if the axes are compatible to get around the Base restriction
# to 1-based vectors
function _reshape(A::AbstractVector, inds::Tuple{OffsetAxis})
    @noinline throw_dimerr(ind::Integer) = throw(
        DimensionMismatch("parent has $(size(A,1)) elements, which is incompatible with length $ind"))
    @noinline throw_dimerr(ind) = throw(
        DimensionMismatch("parent has $(size(A,1)) elements, which is incompatible with indices $ind"))
    _checksize(first(inds), size(A,1)) || throw_dimerr(first(inds))
    A
end
_reshape(A, inds) = _reshape2(A, inds)
_reshape2(A, inds) = reshape(A, inds)
# avoid a stackoverflow by relegating to the parent if no_offset_view returns an offsetarray
_reshape2(A::OffsetArray, inds) = reshape(parent(A), inds)
_reshape_nov(A, inds) = _reshape(no_offset_view(A), inds)

Base.reshape(A::OffsetArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) =
    OffsetArray(_reshape(parent(A), inds), map(_toaxis, inds))
# And for non-offset axes, we can just return a reshape of the parent directly
Base.reshape(A::OffsetArray, inds::Tuple{Union{Integer,Base.OneTo},Vararg{Union{Integer,Base.OneTo}}}) = _reshape_nov(A, inds)
Base.reshape(A::OffsetArray, inds::Dims) = _reshape_nov(A, inds)
Base.reshape(A::OffsetVector, ::Colon) = A
Base.reshape(A::OffsetVector, ::Tuple{Colon}) = A
Base.reshape(A::OffsetArray, ::Colon) = reshape(A, (Colon(),))
Base.reshape(A::OffsetArray, inds::Union{Int,Colon}...) = reshape(A, inds)
Base.reshape(A::OffsetArray, inds::Tuple{Vararg{Union{Int,Colon}}}) = _reshape_nov(A, inds)

# permutedims in Base does not preserve axes, and can not be fixed in a non-breaking way
# This is a stopgap solution
Base.permutedims(v::OffsetVector) = reshape(v, (1, axes(v, 1)))

Base.fill(v, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(Array{typeof(v)}, inds), v)
Base.zeros(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(similar(Array{T}, inds), zero(T))
Base.ones(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(similar(Array{T}, inds), one(T))
Base.trues(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(BitArray, inds), true)
Base.falses(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(BitArray, inds), false)

Base.zero(A::OffsetArray) = parent_call(zero, A)
Base.fill!(A::OffsetArray, x) = parent_call(Ap -> fill!(Ap, x), A)

## Indexing

# Note this gets the index of the parent *array*, not the index of the parent *range*
# Here's how one can think about this:
#   Δi = i - first(r)
#   i′ = first(r.parent) + Δi
# and one obtains the result below.
parentindex(r::IdOffsetRange, i) = i - r.offset

@propagate_inbounds Base.getindex(A::OffsetArray{<:Any,0})  = A.parent[]

@inline function Base.getindex(A::OffsetArray{<:Any,N}, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(A, I...)
    J = map(parentindex, axes(A), I)
    @inbounds parent(A)[J...]
end

@propagate_inbounds Base.getindex(A::OffsetArray{<:Any,N}, c::Vararg{Colon,N}) where N =
    parent_call(x -> getindex(x, c...), A)

# With one Colon we use linear indexing.
# In this case we may forward the index to the parent, as the information about the axes is lost
# The exception to this is with OffsetVectors where the axis information is preserved,
# but that case is handled by getindex(::OffsetArray{<:Any,N}, ::Vararg{Colon,N})
@propagate_inbounds Base.getindex(A::OffsetArray, c::Colon) = A.parent[:]

@inline function Base.getindex(A::OffsetVector, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[parentindex(Base.axes1(A), i)]
end
@propagate_inbounds Base.getindex(A::OffsetArray, i::Int)  = parent(A)[i]

@inline function Base.setindex!(A::OffsetArray{T,N}, val, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    J = map(parentindex, axes(A), I)
    @inbounds parent(A)[J...] = val
    A
end

@inline function Base.setindex!(A::OffsetVector, val, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[parentindex(Base.axes1(A), i)] = val
    A
end
@propagate_inbounds function Base.setindex!(A::OffsetArray, val, i::Int)
    parent(A)[i] = val
    A
end

@inline Base.iterate(a::OffsetArray, i...) = iterate(parent(a), i...)

Base.in(x, A::OffsetArray) = in(x, parent(A))
Base.copy(A::OffsetArray) = parent_call(copy, A)

Base.strides(A::OffsetArray) = strides(parent(A))
Base.elsize(::Type{OffsetArray{T,N,A}}) where {T,N,A} = Base.elsize(A)
@inline Base.unsafe_convert(::Type{Ptr{T}}, A::OffsetArray{T}) where {T} = Base.unsafe_convert(Ptr{T}, parent(A))

# For fast broadcasting: ref https://discourse.julialang.org/t/why-is-there-a-performance-hit-on-broadcasting-with-offsetarrays/32194
Base.dataids(A::OffsetArray) = Base.dataids(parent(A))
Broadcast.broadcast_unalias(dest::OffsetArray, src::OffsetArray) = parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src)

### Special handling for AbstractRange
const OffsetRange{T} = OffsetVector{T,<:AbstractRange{T}}
const OffsetUnitRange{T} = OffsetVector{T,<:AbstractUnitRange{T}}

Base.step(a::OffsetRange) = step(parent(a))

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, or::OffsetRange) = Base.checkindex(Bool, inds, parent(or))

# Certain special methods for linear indexing with integer ranges (or OffsetRanges)
# These may bypass the default getindex(A, I...) pathway if the parent types permit this
# For example AbstractUnitRanges and Arrays have special linear indexing behavior defined

# If both the arguments are offset, we may unwrap the indices to call (::OffsetArray)[::AbstractRange{Int}]
@propagate_inbounds function Base.getindex(A::OffsetArray, r::OffsetRange{Int})
    _indexedby(A[parent(r)], axes(r))
end
# If the indices are offset, we may unwrap them and pass the parent to getindex
@propagate_inbounds function Base.getindex(A::AbstractRange, r::OffsetRange{Int})
    _indexedby(A[parent(r)], axes(r))
end

# An OffsetUnitRange might use the rapid getindex(::Array, ::AbstractUnitRange{Int}) for contiguous indexing
@propagate_inbounds function Base.getindex(A::Array, r::OffsetUnitRange{Int})
    B = A[_contiguousindexingtype(parent(r))]
    OffsetArray(B, axes(r), checkoverflow = false)
end

# avoid hitting the slow method getindex(::Array, ::AbstractRange{Int})
# instead use the faster getindex(::Array, ::UnitRange{Int})
if VERSION <= v"1.7.0-DEV.1039"
    @propagate_inbounds function Base.getindex(A::Array, r::Union{IdOffsetRange, IIUR})
        B = A[_contiguousindexingtype(r)]
        _indexedby(B, axes(r))
    end
end

# Linear Indexing of OffsetArrays with AbstractUnitRanges may use the faster contiguous indexing methods
@inline function Base.getindex(A::OffsetArray, r::AbstractUnitRange{Int})
    @boundscheck checkbounds(A, r)
    # nD OffsetArrays do not have their linear indices shifted, so we may forward the indices provided to the parent
    @inbounds B = parent(A)[_contiguousindexingtype(r)]
    _indexedby(B, axes(r))
end
@inline function Base.getindex(A::OffsetVector, r::AbstractUnitRange{Int})
    @boundscheck checkbounds(A, r)
    # OffsetVectors may have their linear indices shifted, so we subtract the offset from the indices provided
    @inbounds B = parent(A)[_subtractoffset(r, A.offsets[1])]
    _indexedby(B, axes(r))
end

# This method added mainly to index an OffsetRange with another range
@inline function Base.getindex(A::OffsetVector, r::AbstractRange{Int})
    @boundscheck checkbounds(A, r)
    @inbounds B = parent(A)[_subtractoffset(r, A.offsets[1])]
    _indexedby(B, axes(r))
end

# In general we would pass through getindex(A, I...) which calls to_indices(A, I) and finally to_index(I)
# An OffsetUnitRange{Int} has an equivalent IdOffsetRange with the same values and axes,
# something similar also holds for OffsetUnitRange{BigInt}
# We may replace the former with the latter in an indexing operation to obtain a performance boost
@inline function Base.to_index(r::OffsetUnitRange{<:Union{Int,BigInt}})
    of = first(axes(r,1)) - 1
    IdOffsetRange(_subtractoffset(parent(r), of), of)
end

@inline function _boundscheck_index_retaining_axes(r, s)
    @boundscheck checkbounds(r, s)
    @inbounds pr = r[UnitRange(s)]
    _indexedby(pr, axes(s))
end
@inline _boundscheck_return(r, s) = (@boundscheck checkbounds(r, s); s)

for OR in [:IIUR, :IdOffsetRange]
    for R in [:StepRange, :StepRangeLen, :LinRange, :UnitRange]
        @eval @inline Base.getindex(r::$R, s::$OR) = _boundscheck_index_retaining_axes(r, s)
    end

    # this method is needed for ambiguity resolution
    @eval @inline function Base.getindex(r::StepRangeLen{T,<:Base.TwicePrecision,<:Base.TwicePrecision}, s::$OR) where T
        _boundscheck_index_retaining_axes(r, s)
    end
end
Base.getindex(r::Base.OneTo, s::IdOffsetRange) = _boundscheck_index_retaining_axes(r, s)
if VERSION < v"1.7.0-beta2"
    Base.getindex(r::Base.OneTo, s::IIUR) = _boundscheck_index_retaining_axes(r, s)
end

# These methods are added to avoid ambiguities with Base.
# The ones involving Base types should be ported to Base and version-limited here
@inline Base.getindex(r::IdentityUnitRange, s::IIUR) = _boundscheck_return(r, s)
@inline Base.getindex(r::IdentityUnitRange, s::IdOffsetRange) = _boundscheck_return(r, s)
if IdentityUnitRange !== Base.Slice
    @inline Base.getindex(r::Base.Slice, s::IIUR) = _boundscheck_return(r, s)
    @inline Base.getindex(r::Base.Slice, s::IdOffsetRange) = _boundscheck_return(r, s)
end

# eltype conversion
# This may use specialized map methods for the parent
Base.map(::Type{T}, O::OffsetArray) where {T} = parent_call(x -> map(T, x), O)
Base.map(::Type{T}, r::IdOffsetRange) where {T<:Real} = _indexedby(map(T, UnitRange(r)), axes(r))
if eltype(IIUR) === Int
    # This is type-piracy, but there is no way to convert an IdentityUnitRange to a non-Int type in Base
    Base.map(::Type{T}, r::IdentityUnitRange) where {T<:Real} = _indexedby(map(T, UnitRange(r)), axes(r))
end

# mapreduce is faster with an IdOffsetRange than with an OffsetUnitRange
# We therefore convert OffsetUnitRanges to IdOffsetRanges with the same values and axes
function Base.mapreduce(f, op, As::OffsetUnitRange{<:Integer}...; kw...)
    ofs = map(A -> first(axes(A,1)) - 1, As)
    AIds = map((A, of) -> IdOffsetRange(_subtractoffset(parent(A), of), of), As, ofs)
    mapreduce(f, op, AIds...; kw...)
end

# Optimize certain reductions that treat an OffsetVector as a list
for f in [:minimum, :maximum, :extrema, :sum]
    @eval Base.$f(r::OffsetRange) = $f(parent(r))
end

function Base.show(io::IO, r::OffsetRange)
    show(io, r.parent)
    print(io, " with indices ", UnitRange(axes(r, 1)))
end
Base.show(io::IO, ::MIME"text/plain", r::OffsetRange) = show(io, r)

### Some mutating functions defined only for OffsetVector ###

Base.resize!(A::OffsetVector, nl::Integer) = (resize!(A.parent, nl); A)
Base.push!(A::OffsetVector, x...) = (push!(A.parent, x...); A)
Base.pop!(A::OffsetVector) = pop!(A.parent)
Base.append!(A::OffsetVector, items) = (append!(A.parent, items); A)
Base.empty!(A::OffsetVector) = (empty!(A.parent); A)

# These functions keep the summary compact
function Base.inds2string(inds::Tuple{Vararg{Union{IdOffsetRange, IdentityUnitRange{<:IdOffsetRange}}}})
    Base.inds2string(map(UnitRange, inds))
end
Base.showindices(io::IO, ind1::IdOffsetRange, inds::IdOffsetRange...) = Base.showindices(io, map(UnitRange, (ind1, inds...))...)

function Base.showarg(io::IO, @nospecialize(a::OffsetArray), toplevel)
    print(io, "OffsetArray(")
    Base.showarg(io, parent(a), false)
    Base.showindices(io, axes(a)...)
    print(io, ')')
    if toplevel
        print(io, " with eltype ", eltype(a))
    end
end

function Base.replace_in_print_matrix(A::OffsetArray{<:Any,2}, i::Integer, j::Integer, s::AbstractString)
    J = map(parentindex, axes(A), (i,j))
    Base.replace_in_print_matrix(parent(A), J..., s)
end
function Base.replace_in_print_matrix(A::OffsetArray{<:Any,1}, i::Integer, j::Integer, s::AbstractString)
    ip = parentindex(axes(A,1), i)
    Base.replace_in_print_matrix(parent(A), ip, j, s)
end

"""
    no_offset_view(A)

Return an `AbstractArray` that shares structure and underlying data with the argument,
but uses 1-based indexing. May just return the argument when applicable.
Not exported.

The default implementation uses `OffsetArrays`, but other types should use something more
specific to remove a level of indirection when applicable.

```jldoctest; setup=:(using OffsetArrays)
julia> A = [1 3 5; 2 4 6];

julia> O = OffsetArray(A, 0:1, -1:1)
2×3 OffsetArray(::$(Matrix{Int}), 0:1, -1:1) with eltype $Int with indices 0:1×-1:1:
 1  3  5
 2  4  6

julia> OffsetArrays.no_offset_view(O)[1,1] = -9
-9

julia> A
2×3 $(Matrix{Int}):
 -9  3  5
  2  4  6
```
"""
no_offset_view(A::OffsetArray) = no_offset_view(parent(A))
if isdefined(Base, :IdentityUnitRange)
    # valid only if Slice is distinguished from IdentityUnitRange
    no_offset_view(a::Base.Slice{<:Base.OneTo}) = a
    no_offset_view(a::Base.Slice) = Base.Slice(UnitRange(a))
    no_offset_view(S::SubArray) = view(parent(S), map(no_offset_view, parentindices(S))...)
end
no_offset_view(a::Array) = a
no_offset_view(i::Number) = i
no_offset_view(A::AbstractArray) = _no_offset_view(axes(A), A)
_no_offset_view(::Tuple{}, A::AbstractArray{T,0}) where T = A
_no_offset_view(::Tuple{Base.OneTo, Vararg{Base.OneTo}}, A::AbstractArray) = A
# the following method is needed for ambiguity resolution
_no_offset_view(::Tuple{Base.OneTo, Vararg{Base.OneTo}}, A::AbstractUnitRange) = A
_no_offset_view(::Any, A::AbstractArray) = OffsetArray(A, Origin(1))
_no_offset_view(::Any, A::AbstractUnitRange) = UnitRange(A)

#####
# center/centered
# These two helpers are deliberately not exported; their meaning can be very different in
# other scenarios and will be very likely to cause name conflicts if exported.
#####
"""
    center(A, [r::RoundingMode=RoundDown])::Dims

Return the center coordinate of given array `A`. If `size(A, k)` is even,
a rounding procedure will be applied with mode `r`.

!!! compat "OffsetArrays 1.9"
    This method requires at least OffsetArrays 1.9.

# Examples

```jldoctest; setup=:(using OffsetArrays)
julia> A = reshape(collect(1:9), 3, 3)
3×3 $(Matrix{Int}):
 1  4  7
 2  5  8
 3  6  9

julia> c = OffsetArrays.center(A)
(2, 2)

julia> A[c...]
5

julia> Ao = OffsetArray(A, -2, -2); # axes (-1:1, -1:1)

julia> c = OffsetArrays.center(Ao)
(0, 0)

julia> Ao[c...]
5
```

To shift the center coordinate of the given array to `(0, 0, ...)`, you
can use [`centered`](@ref OffsetArrays.centered).
"""
function center(A::AbstractArray, r::RoundingMode=RoundDown)
    map(axes(A)) do inds
        round(Int, (length(inds)-1)/2, r) + first(inds)
    end
end

"""
    centered(A, cp=center(A)) -> Ao

Shift the center coordinate/point `cp` of array `A` to `(0, 0, ..., 0)`. Internally, this is
equivalent to `OffsetArray(A, .-cp)`.

!!! compat "OffsetArrays 1.9"
    This method requires at least OffsetArrays 1.9.

# Examples

```jldoctest; setup=:(using OffsetArrays)
julia> A = reshape(collect(1:9), 3, 3)
3×3 $(Matrix{Int}):
 1  4  7
 2  5  8
 3  6  9

julia> Ao = OffsetArrays.centered(A); # axes (-1:1, -1:1)

julia> Ao[0, 0]
5

julia> Ao = OffsetArray(A, OffsetArrays.Origin(0)); # axes (0:2, 0:2)

julia> Aoo = OffsetArrays.centered(Ao); # axes (-1:1, -1:1)

julia> Aoo[0, 0]
5
```

Users are allowed to pass `cp` to change how "center point" is interpreted, but the meaning of the
output array should be reinterpreted as well. For instance, if `cp = map(last, axes(A))` then this
function no longer shifts the center point but instead the bottom-right point to `(0, 0, ..., 0)`.
A commonly usage of `cp` is to change the rounding behavior when the array is of even size at some
dimension:

```jldoctest; setup=:(using OffsetArrays)
julia> A = reshape(collect(1:4), 2, 2) # Ideally the center should be (1.5, 1.5) but OffsetArrays only support integer offsets
2×2 $(Matrix{Int}):
 1  3
 2  4

julia> OffsetArrays.centered(A, OffsetArrays.center(A, RoundUp)) # set (2, 2) as the center point
2×2 OffsetArray(::$(Matrix{Int}), -1:0, -1:0) with eltype $(Int) with indices -1:0×-1:0:
 1  3
 2  4

julia> OffsetArrays.centered(A, OffsetArrays.center(A, RoundDown)) # set (1, 1) as the center point
2×2 OffsetArray(::$(Matrix{Int}), 0:1, 0:1) with eltype $(Int) with indices 0:1×0:1:
 1  3
 2  4
```

See also [`center`](@ref OffsetArrays.center).
"""
centered(A::AbstractArray, cp::Dims=center(A)) = OffsetArray(A, .-cp)

centered(A::AbstractArray, i::CartesianIndex) = centered(A, Tuple(i))

####
# work around for segfault in searchsorted*
#  https://github.com/JuliaLang/julia/issues/33977
####

function _safe_searchsorted(v::OffsetArray, x, ilo::T, ihi::T, o::Base.Ordering) where T<:Integer
    u = T(1)
    lo = ilo - u
    hi = ihi + u
    @inbounds while lo < hi - u
        m = (lo + hi) ÷ 2
        if Base.lt(o, v[m], x)
            lo = m
        elseif Base.lt(o, x, v[m])
            hi = m
        else
            a = searchsortedfirst(v, x, max(lo,ilo), m, o)
            b = searchsortedlast(v, x, m, min(hi,ihi), o)
            return a : b
        end
    end
    return (lo + 1) : (hi - 1)
end
function _safe_searchsortedfirst(v::OffsetArray, x, lo::T, hi::T, o::Base.Ordering) where T<:Integer
    u = T(1)
    lo = lo - u
    hi = hi + u
    @inbounds while lo < hi - u
        m = (lo + hi) ÷ 2
        if Base.lt(o, v[m], x)
            lo = m
        else
            hi = m
        end
    end
    return hi
end
function _safe_searchsortedlast(v::OffsetArray, x, lo::T, hi::T, o::Base.Ordering) where T<:Integer
    u = T(1)
    lo = lo - u
    hi = hi + u
    @inbounds while lo < hi - u
        m = (lo + hi) ÷ 2
        if Base.lt(o, x, v[m])
            hi = m
        else
            lo = m
        end
    end
    return lo
end

if VERSION ≤ v"1.2"
    # ambiguity warnings in earlier versions
    Base.searchsorted(v::OffsetArray, x, ilo::Int, ihi::Int, o::Base.Ordering) =
        _safe_searchsorted(v, x, ilo, ihi, o)
    Base.searchsortedfirst(v::OffsetArray, x, lo::Int, hi::Int, o::Base.Ordering) =
        _safe_searchsortedfirst(v, x, lo, hi, o)
    Base.searchsortedlast(v::OffsetArray, x, lo::Int, hi::Int, o::Base.Ordering) =
        _safe_searchsortedlast(v, x, lo, hi, o)
end

Base.searchsorted(v::OffsetArray, x, ilo::T, ihi::T, o::Base.Ordering) where T<:Integer =
    _safe_searchsorted(v, x, ilo, ihi, o)
Base.searchsortedfirst(v::OffsetArray, x, lo::T, hi::T, o::Base.Ordering) where T<:Integer =
    _safe_searchsortedfirst(v, x, lo, hi, o)
Base.searchsortedlast(v::OffsetArray, x, lo::T, hi::T, o::Base.Ordering) where T<:Integer =
    _safe_searchsortedlast(v, x, lo, hi, o)

if VERSION < v"1.1.0-DEV.783"
    Base.copyfirst!(dest::OffsetArray, src::OffsetArray) = (maximum!(parent(dest), parent(src)); return dest)
end

if VERSION <= v"1.7.0-DEV.400"
    # https://github.com/JuliaLang/julia/pull/39393
    # index for zero-argument getindex should be first linear index instead of 1 (#194)
    Base._to_linear_index(A::OffsetArray) = first(LinearIndices(A))
end

##
# Adapt allows for automatic conversion of CPU OffsetArrays to GPU OffsetArrays
##
import Adapt
Adapt.adapt_structure(to, O::OffsetArray) = parent_call(x -> Adapt.adapt(to, x), O)

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

##
# Deprecations
##

# This is a bad API design as it introduces counter intuitive results (#250)
@deprecate centered(A::AbstractArray, r::RoundingMode) OffsetArray(A, .-center(A, r)) false

end # module


text:
function computeStress(mesh::FluidMesh,fluidsol::FluidSolution,param::P) where P<:AbstractConstantParameter
    Nnodes = length(mesh.xy)
    Nel    = length(mesh.cm)
    σxx = zeros(Float64,Nnodes)
    σxy = zeros(Float64,Nnodes)
    σyy = zeros(Float64,Nnodes)

    xN = zeros(Float64,9)
    yN = zeros(Float64,9)

    nodeCounter = zeros(Int,Nnodes)

    # compute pressure at Q2 mesh
    xyArr = NodesToArray(mesh.xy)
    xypArr = NodesToArray(mesh.xyp)
    cmpArr = ElementsToArray(mesh.cmp)
    p = meshTransform(xypArr,cmpArr,fluidsol.p,xyArr)

    s = [-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,-1.0,0.0]
    t = [-1.0,-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,0.0]

    for el=1:Nel
        c = mesh.cm[el].NodeList
        xN = [mesh.xy[i].x for i in c]
        yN = [mesh.xy[i].y for i in c]

        for ti=1:9
            phi,dphidx,dphidy,jac = derivShape2D(s[ti],t[ti],xN,yN,2)

            # interpolate values at gauss points
            dudxg = shapeEval(fluidsol.u[c],dphidx)/jac
            dudyg = shapeEval(fluidsol.u[c],dphidy)/jac
            dvdxg = shapeEval(fluidsol.v[c],dphidx)/jac
            dvdyg = shapeEval(fluidsol.v[c],dphidy)/jac

            σxx[c[ti]] += -p[c[ti]] + 2*param.μ*dudxg
            σxy[c[ti]] += param.μ*(dudyg + dvdxg)
            σyy[c[ti]] += -p[c[ti]] + 2*param.μ*dvdyg

            nodeCounter[c[ti]] += 1
        end
    end

    # divide through by number of times contributed (finds average at each node from surrounding elements)
    for ti=1:Nnodes
        σxx[ti] = σxx[ti]/nodeCounter[ti]
        σxy[ti] = σxy[ti]/nodeCounter[ti]
        σyy[ti] = σyy[ti]/nodeCounter[ti]
    end

    return σxx, σxy, σyy
end

function computeStress(mesh::FluidMesh,fluidsol::FluidSolution,param::P) where P<:AbstractVariableParameter
    Nnodes = length(mesh.xy)
    Nel    = length(mesh.cm)
    σxx = zeros(Float64,Nnodes)
    σxy = zeros(Float64,Nnodes)
    σyy = zeros(Float64,Nnodes)

    xN = zeros(Float64,9)
    yN = zeros(Float64,9)

    nodeCounter = zeros(Int,Nnodes)

    # compute pressure at Q2 mesh
    xyArr = NodesToArray(mesh.xy)
    xypArr = NodesToArray(mesh.xyp)
    cmpArr = ElementsToArray(mesh.cmp)
    p = meshTransform(xypArr,cmpArr,fluidsol.p,xyArr)

    s = [-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,-1.0,0.0]
    t = [-1.0,-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,0.0]

    for el=1:Nel
        c = mesh.cm[el].NodeList
        xN = [mesh.xy[i].x for i in c]
        yN = [mesh.xy[i].y for i in c]

        for ti=1:9
            phi,dphidx,dphidy,jac = derivShape2D(s[ti],t[ti],xN,yN,2)

            # interpolate values at gauss points
            dudxg = shapeEval(fluidsol.u[c],dphidx)/jac
            dudyg = shapeEval(fluidsol.u[c],dphidy)/jac
            dvdxg = shapeEval(fluidsol.v[c],dphidx)/jac
            dvdyg = shapeEval(fluidsol.v[c],dphidy)/jac
            μg    = shapeEval(param.μ[c],phi)

            σxx[c[ti]] += -p[c[ti]] + 2*μg*dudxg
            σxy[c[ti]] += μg*(dudyg + dvdxg)
            σyy[c[ti]] += -p[c[ti]] + 2*μg*dvdyg

            nodeCounter[c[ti]] += 1
        end
    end

    # divide through by number of times contributed (finds average at each node from surrounding elements)
    for ti=1:Nnodes
        σxx[ti] = σxx[ti]/nodeCounter[ti]
        σxy[ti] = σxy[ti]/nodeCounter[ti]
        σyy[ti] = σyy[ti]/nodeCounter[ti]
    end

    return σxx, σxy, σyy
end

function computeDerivative(mesh::M,u::Vector{Float64},var::Symbol) where M<:AbstractMesh
    Nnodes = length(mesh.xy)
    Nel    = length(mesh.cm)
    dudi = zeros(Float64,Nnodes)

    xN = zeros(Float64,9)
    yN = zeros(Float64,9)

    nodeCounter = zeros(Int,Nnodes)

    s = [-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,-1.0,0.0]
    t = [-1.0,-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,0.0]

    for el=1:Nel
        c = mesh.cm[el].NodeList
        xN = [mesh.xy[i].x for i in c]
        yN = [mesh.xy[i].y for i in c]

        for ti=1:9
            phi,dphidx,dphidy,jac = derivShape2D(s[ti],t[ti],xN,yN,2)

            # interpolate values at gauss points
            if var == :dx
                dudig = shapeEval(u[c],dphidx)/jac
            elseif var == :dy
                dudig = shapeEval(u[c],dphidy)/jac
            else
                error("variable ain't right")
            end

            dudi[c[ti]] += dudig
            nodeCounter[c[ti]] += 1
        end
    end

    # divide through by number of times contributed (finds average at each node from surrounding elements)
    for ti=1:Nnodes
        dudi[ti] = dudi[ti]/nodeCounter[ti]
    end

    return dudi
end

function computeExtension(mesh::FluidMesh,fluidsol::FluidSolution,param)
    nNodes = length(mesh.xy)

    σxx,σxy,σyy = computeStress(mesh,fluidsol,param)

    λField = zeros(Float64,nNodes)
    evecXField = zeros(Float64,nNodes)
    evecYField = zeros(Float64,nNodes)

    for ti=1:nNodes
        # construct stress tensor
        T = [σxx[ti] σxy[ti]; σxy[ti] σyy[ti]]

        # find eigenvalues and eigenvectors of stress tensor T
        F = eigfact(T)

        # find maximum eigenvalue
        λField[ti]     = abs(maximum(abs.(F[:values])))
    end

    return λField
end

function computeCompressibility(mesh::M,fSol::FluidSolution) where M<:AbstractMesh
    u = fSol.u; v=fSol.v

    dudx = computeDerivative(mesh,fSol.u,:dx)
    dvdy = computeDerivative(mesh,fSol.v,:dy)

    compress = dudx .+ dvdy

    return compress
end

ext: jl


text:
# Test the precision of the series expansion
function test_impact_parameter_accuracy(n)
    # Setup the simulation
    BJD = 7250.0; t0_ic = 7.0
    tmax = 2.0
    ic = setup_ICs(n,BJD,t0_ic)
    intr = setup_integrator(ic, tmax)
    tt = compute_transit_times(ic, intr)
    pd = compute_pd(ic, tt, intr)

    # Normalize points to stellar radius
    rstar = 0.00465047 * 0.1192 # Trappist-1 (Rstar/AU)
    normalize_points!(pd.points, rstar)

    # Get radius ratios
    k = get_radius_ratios_trappist(n)

    for it in 1:length(pd.times)
        # Compute the transit duration
        ib = pd.bodies[it]
        t0 = pd.times[it];
        s = State(ic)
        els = NbodyGradient.get_orbital_elements(s, ic)
        b0 = compute_impact_parameter(t0, t0, 2e-2, pd.points[ib,it,:,:])
        a = els[ib].a / rstar

        # Compute impact parameter over transit duration
        resolution = 0.0001
        bound = els[ib].P / pi * asin(sqrt((1 + k[ib - 1])^2 + b0^2) / (a * sin(els[ib].I))) / 2
        tc = collect(pd.times[it] - bound:resolution:pd.times[it] + bound)
        xc = components(@views(pd.points[ib,it,:,1]), 2e-2)
        yc = components(@views(pd.points[ib,it,:,2]), 2e-2)
        b_series = compute_impact_parameter.(tc, t0, Ref(xc), Ref(yc))

        # Now output the actual sky-separation from NbodyGradient
        positions = zeros(length(tc), 2)
        for (i, t) in enumerate(tc)
            s = State(ic)
            intr(s, t,;grad=false)
            positions[i,:] .= s.x[1:2,ib] ./ rstar
        end
        b_actual = sqrt.(positions[:,1].^2 .+ positions[:,2].^2)

        @test isapprox(b_actual, b_series, norm=x -> maximum(abs.(x)))
    end
end

function test_impact_parameter_derivatives(n)
    # Test the derivatives of the impact parameter with respect
    # to the initial cartesian coordinates and masses

    # Setup the simulation
    BJD = 7250.0; t0_ic = 7.0
    tmax = 2.0
    ic = setup_ICs(n, BJD, t0_ic)
    intr = setup_integrator(ic, tmax)
    tt = compute_transit_times(ic, intr)
    pd = compute_pd(ic, tt, intr, grad=true)

    # Now run BigFloat precision simulations
    ic_copy = setup_ICs(n, big(BJD), big(t0_ic))
    intr_copy = setup_integrator(ic_copy, big(tmax))
    s_copy = State(ic_copy)
    tt_copy = compute_transit_times(ic_copy, intr_copy)

    # Compute the impact parameter
    # Used in finite difference derivatives
    function compute_b(coords, s_copy, tt, intr, ic, it)
        T = eltype(coords)
        h::T = 2e-2
        s = deepcopy(s_copy)
        coords = reshape(coords, 7, n)
        for (i, col) in enumerate(eachcol(coords))
            s.x[:,i] .= col[1:3]
            s.v[:,i] .= col[4:6]
            s.m[i]    = col[end]
        end

        pd = compute_pd(s, deepcopy(ic), deepcopy(tt), intr, grad=false)

        ib = pd.bodies[it]
        t0 = pd.times[it]
        points = pd.points[ib,it,:,:]
        xc = components(points[:,1], h)
        yc = components(points[:,2], h)
        return compute_impact_parameter(t0, t0, xc, yc)
    end

    # Compute the impact parameter and analytic gradient
    function compute_grad_b(coords, s_copy, tt, intr, ic, it)
        h = 2e-2
        s = deepcopy(s_copy)
        coords = reshape(coords, 7, n)
        for (i, col) in enumerate(eachcol(coords))
            s.x[:,i] .= col[1:3]
            s.v[:,i] .= col[4:6]
            s.m[i]    = col[end]
        end

        pd = compute_pd(s, deepcopy(ic), deepcopy(tt), intr, grad=true)

        ib = pd.bodies[it]
        t0 = pd.times[it]
        points = pd.points[ib,it,:,:]
        xc = components(points[:,1], h)
        yc = components(points[:,2], h)
        N = length(coords)
        grad = zeros(N)
        inds = [1,2,3,4,5,6,7]
        dxc = [components(pd.dpoints[ib,it,:,1,k,i], h) for i in inds, k in 1:ic.nbody][:]
        dyc = [components(pd.dpoints[ib,it,:,2,k,i], h) for i in inds, k in 1:ic.nbody][:]
        compute_impact_parameter!(t0, t0, xc, yc, dxc, dyc, grad)
        return grad
    end

    # Compare numerical and analytic derivative wrt the inital cartesian coordinates and masses
    coords = vcat([[x...,v...,m] for (x, v, m) in zip(eachcol(s_copy.x), eachcol(s_copy.v), s_copy.m)]...)
    for it in 1:length(pd.times)
        grad_num = grad(central_fdm(5, 1), c -> compute_b(c, s_copy, tt_copy, intr_copy, ic_copy, it), big.(coords))[1]
        grad_analytic = compute_grad_b(Float64.(coords), State(ic), tt, intr, ic, it)
        # println(maximum(abs.((Float64.(grad_num) .- grad_analytic) ./ grad_num)))
        @test isapprox(asinh.(Float64.(grad_num)), asinh.(grad_analytic), norm=x -> maximum(abs.(x)))
    end
end

@testset "Impact Parameter" begin
    test_impact_parameter_accuracy(8)
    test_impact_parameter_derivatives(3) # Takes a long time...
end


text:
#####
##### record diagnostics of plankton processes at each time step
#####

function gpu_diags_proc_kernel!(diags_proc, proc, ac, x, y, z)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:size(ac,1)
        @inbounds CUDA.@atomic diags_proc[x[i], y[i], z[i]] += proc[i] * ac[i]
    end
end
function diags_proc!(diags_proc, proc, ac, x, y, z, ::GPU)
    @cuda threads=256 blocks=ceil(Int, size(ac,1)/256) gpu_diags_proc_kernel!(diags_proc, proc, ac, x, y, z)
    return nothing
end
function diags_proc!(diags_proc, proc, ac, x, y, z, ::CPU)
    for i in 1:size(ac,1)
        @inbounds diags_proc[x[i], y[i], z[i]] += proc[i] * ac[i]
    end
end

function diags_spcs!(diags_sp, plank, ac, x, y, z, mode::AbstractMode, arch::Architecture)
    diags = (:PS, :resp, :Bm, :Chl)
    if isa(mode, CarbonMode)
        nothing
    elseif isa(mode, QuotaMode)
        diags = (:PS, :BS, :VDOC, :VNH4, :VNO3, :VPO4, :resp, :exu, :Bm, :Cq, :Nq, :Pq, :Chl)
    elseif isa(mode, MacroMolecularMode)
        nothing
    end

    for diag in keys(diags_sp)
        if diag in (:num, :graz, :mort, :dvid)
            nothing
        elseif diag in diags
            diags_proc!(diags_sp[diag], getproperty(plank, diag), ac, x, y, z, arch)
        end
    end
end


text:
using XDF, Downloads, Test

@testset "Minimal XDF file" begin
    url = "https://github.com/xdf-modules/example-files/blob/master/minimal.xdf?raw=true"

    streams = read_xdf(Downloads.download(url))
    @test 0 in keys(streams)
    @test streams[0]["nchannels"] == 3
    @test streams[0]["name"] == "SendDataC"
    @test streams[0]["dtype"] == Int16
    @test streams[0]["srate"] == 10.0
    @test streams[0]["type"] == "EEG"
    @test size(streams[0]["data"]) == (9, 3)

    @test 46202862 in keys(streams)
    @test streams[46202862]["nchannels"] == 1
    @test streams[46202862]["name"] == "SendDataString"
    @test streams[46202862]["dtype"] == String
    @test streams[46202862]["srate"] == 10.0
    @test streams[46202862]["type"] == "StringMarker"
    @test size(streams[46202862]["data"]) == (9, 1)
    @test startswith(streams[46202862]["data"][1, 1], "<?xml version")
    @test endswith(streams[46202862]["data"][1, 1], "</info>")
    @test streams[46202862]["data"][2:end, 1] == ["Hello", "World", "from", "LSL", "Hello",
                                                  "World", "from", "LSL"]
end

@testset "XDF file with clock resets" begin
    url = "https://github.com/xdf-modules/example-files/blob/master/data_with_clock_resets.xdf?raw=true"

    streams = read_xdf(Downloads.download(url))
    @test 1 in keys(streams)
    @test streams[1]["nchannels"] == 1
    @test streams[1]["name"] == "MyMarkerStream"
    @test streams[1]["dtype"] == String
    @test streams[1]["srate"] == 0.0
    @test streams[1]["type"] == "Markers"
    @test size(streams[1]["data"]) == (175, 1)

    @test 2 in keys(streams)
    @test streams[2]["nchannels"] == 8
    @test streams[2]["name"] == "BioSemi"
    @test streams[2]["dtype"] == Float32
    @test streams[2]["srate"] == 100.0
    @test streams[2]["type"] == "EEG"
    @test size(streams[2]["data"]) == (27815, 8)
end


text:
module Lexer

include("Tokens.jl")

function lex(code::String)::Vector{Tokens.AbstractToken}
    result::Vector{Tokens.AbstractToken} = Vector{Tokens.AbstractToken}()
    pos::Int = 1

    while length(code) >= pos
        c::Char = code[pos]

        if isdigit(c)
            number::Vector{Char} = Vector{Char}()
            push!(number, c)

            haspoint::Bool = false
            pos += 1
            while length(code) >= pos
                c = code[pos]
                if isdigit(c)
                    push!(number, c)
                elseif c == '.'
                    if haspoint
                        throw(ErrorException("two decimal points in a number literal not allowed!"))
                    else
                        push!(number, c)
                        haspoint = true
                    end
                else
                    break
                end
                pos += 1
            end

            if haspoint
                push!(result, Tokens.Num(parse(Float64, join(number))))
            else
                push!(result, Tokens.Num(parse(Int64, join(number))))
            end
            pos -= 1

        elseif c == '+'
            push!(result, Tokens.Plus())
        elseif c == '-'
            push!(result, Tokens.Minus())
        elseif c == '*'
            push!(result, Tokens.Mul())
        elseif c == '/'
            push!(result, Tokens.Div())
        elseif c == '^'
            push!(result, Tokens.Pow())

        elseif c == '('
            push!(result, Tokens.RBrace())
        elseif c == ')'
            push!(result, Tokens.LBrace())
        else
            throw(ErrorException("unexpected token '$c'"))
        end

        pos += 1
    end

    return result
end

end

max_stars_repo_name: MhmCats/AsciiMath

max_issues_repo_name: MhmCats/AsciiMath

max_forks_repo_name: MhmCats/AsciiMath


text:
const MemoryData = NTuple{5,Vector{Int}}

abstract type AgentMemory end

Base.@kwdef struct CircularBufferMemory <: AgentMemory
    data::CircularBuffer = CircularBuffer{MemoryData}(100_000)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 7c860884-5afd-4b4a-b358-4f2bce8bd395
using Pkg; Pkg.activate("../..")

# ╔═╡ 38c4e220-d708-11ec-3968-fbbd41c26155
using PlutoUI, Plots, DifferentialEquations, NonLinearDynamicsCourse

# ╔═╡ 2ac364c2-cbdd-49b4-9f26-9fe89382be5e
TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)

# ╔═╡ 6b6a5e19-5298-4969-b1c9-699aa1cb2996
html"""

"""

# ╔═╡ 62cba86d-4406-4d16-8715-b29bee7d3178
md"""
# Logistic Equation

$\dot{x}=Rx\left(1-\frac{x}{K}\right)$
"""

# ╔═╡ 6d64cc27-1059-4dd7-8f17-1c6fd9a2eca0
html"""

"""

# ╔═╡ 4dc32b5a-5794-4d2d-844e-cdb80bec2e7b
logistic(x,p,t)=p[1]*x*(1.0-x/p[2])

# ╔═╡ 77766c99-4142-4acb-a855-ab133e67aa66
@bind pars2 (
	PlutoUI.combine() do bind
		md"""
		R: $(bind(Slider(0:0.02:2.0,default=0.1;show_value=true))) \
		K: $(bind(Slider(0.01:0.02:2.0,default=0.1;show_value=true))) \
		x0: $(bind(Slider(0:0.02:2.0,default=0.1;show_value=true)))
		"""
	end
)

# ╔═╡ 16a156b5-92ab-4e1c-ab70-c3b348186c57
flux1D(logistic,pars2[3],100.0,pars2;xlims=[-0.1,2.0],title="Logistic")

# ╔═╡ ce4da9bf-316b-45d2-8504-f2a63cdda247
md"""
# Logistic Equation with Harvest

$\dot{x}=Rx\left(1-\frac{x}{K}\right)-H$
"""

# ╔═╡ 29474ca0-839b-45de-a7bf-7d91c5ea59aa
logharvest1(x,p,t)=p[1]*x*(1.0-x/p[2])-p[3]

# ╔═╡ 8fc407f4-a242-4cdc-b7e0-062c8f5782d1
@bind pars_harvest (
	PlutoUI.combine() do bind
		md"""
		R: $(bind(Slider(0:0.02:2.0,default=0.1;show_value=true))) \
		K: $(bind(Slider(0.01:0.02:2.0,default=0.1;show_value=true))) \
		H: $(bind(Slider(0.01:0.02:2.0,default=0.1;show_value=true))) \
		x0: $(bind(Slider(0:0.02:2.0,default=0.1;show_value=true)))
		"""
	end
)

# ╔═╡ 794936ae-d9ae-4319-994f-4282c99c4238
flux1D(logharvest1,pars_harvest[4],300.0,pars_harvest,(u)->(u<0);xlims=[0.0,2.0],title="Logistic with Harvest")

# ╔═╡ b13d261e-1c1b-41b0-ae86-e20cd2304727
md"""
## Critical Slowing Down
"""

# ╔═╡ c3757398-060c-412c-997b-58331b2dee04
@bind pars_csd (
	PlutoUI.combine() do bind
		md"""
		H: $(bind(Slider(0.2:0.002:0.25,default=0.1;show_value=true))) \
		S: $(bind(Slider(0:0.001:0.1,default=0.0;show_value=true)))
		"""
	end
)

# ╔═╡ f1f04dd4-0447-4e48-9bab-f34058a6a0f1
flux1D(logharvest1,0.5+sqrt(0.25-pars_csd[1]),200.0,[1.0,1.0,pars_csd[1]],10.0,pars_csd[2],(u)->(u<0);xlims=[0.0,1.0],title="Log whith Harvest perturbed")

# ╔═╡ c803ed90-958c-4d0e-84d5-3e8cc9c8fbfe
md"""
# Consumer Equation

$\dot{x} = Rx\left(1-\frac{x}{K}\right) - Px$
"""

# ╔═╡ b4ea8364-a582-443c-bdb1-46e75a5c4d4e
# Consumer Equation
consumer(x,p,t)=p[1]*x*(1.0-x/p[2])-p[3]*x

# ╔═╡ 5afa3aff-af15-429b-a0fb-9a0dfcad740d
@bind pars_consumer (
	PlutoUI.combine() do bind
		md"""
		R: $(bind(Slider(0:0.02:1.0,default=0.1;show_value=true))) \
		K: $(bind(Slider(0.01:0.02:2.0,default=0.1;show_value=true))) \
		P: $(bind(Slider(0.0:0.01:0.5,default=0.1;show_value=true))) \
		x0: $(bind(Slider(0:0.02:2.0,default=0.1;show_value=true)))
		"""
	end
)

# ╔═╡ 5018d342-0a51-4275-8c0f-1b5d253c2ec0
flux1D(consumer,pars_consumer[4],300.0,pars_consumer,xlims=[0.0,2.0],title="Consumer Equation")

# ╔═╡ 84a596c1-3a0a-4edd-a7ea-f3fcd8839c2a
md"""
# Logistic Equation with Outbreak

$\dot{x} = Rx\left(1-\displaystyle\frac{x}{K}\right)-P\displaystyle\frac{x^2}{1+x^2}$
"""

# ╔═╡ 91b2718a-a921-4bae-99c0-91ec9c7e6479
html"""

"""

# ╔═╡ 5a0b156e-b0b2-43b8-b15e-bcdb65c154bf
logoutbreak(x,p,t)=p[1]*x*(1.0-x/p[2])-p[3]*x*x/(1+x*x)

# ╔═╡ 6528501f-3bb2-48a7-b24e-d32400033538
@bind pars_outbreak (
	PlutoUI.combine() do bind
		md"""
		R: $(bind(Slider(0:0.02:2.0,default=0.1;show_value=true))) \
		K: $(bind(Slider(0.01:0.02:10.0,default=0.1;show_value=true))) \
		P: $(bind(Slider(0.0:0.01:1.0,default=0.1;show_value=true))) \
		x0: $(bind(Slider(0:0.02:8.0,default=0.1;show_value=true)))
		"""
	end
)

# ╔═╡ a3270d8b-f167-4fb9-bcc5-c62c91b39649
flux1D(logoutbreak,pars_outbreak[4],300.0,pars_outbreak;xlims=[-0.2,8.0],title="Log with Outbreak")

# ╔═╡ 0d0d6d2f-56b2-4e01-b9a4-6e469eac90ad
html"""

"""

# ╔═╡ 76cfd5c8-6298-4ae2-9d2d-9e60deb12f49
html"""

"""

# ╔═╡ 0bf40552-b627-4fa7-b3ba-1d8c7ed4a568
html"""

"""

# ╔═╡ Cell order:
# ╟─7c860884-5afd-4b4a-b358-4f2bce8bd395
# ╠═38c4e220-d708-11ec-3968-fbbd41c26155
# ╟─2ac364c2-cbdd-49b4-9f26-9fe89382be5e
# ╟─6b6a5e19-5298-4969-b1c9-699aa1cb2996
# ╟─62cba86d-4406-4d16-8715-b29bee7d3178
# ╟─6d64cc27-1059-4dd7-8f17-1c6fd9a2eca0
# ╠═4dc32b5a-5794-4d2d-844e-cdb80bec2e7b
# ╟─16a156b5-92ab-4e1c-ab70-c3b348186c57
# ╟─77766c99-4142-4acb-a855-ab133e67aa66
# ╟─ce4da9bf-316b-45d2-8504-f2a63cdda247
# ╠═29474ca0-839b-45de-a7bf-7d91c5ea59aa
# ╟─794936ae-d9ae-4319-994f-4282c99c4238
# ╟─8fc407f4-a242-4cdc-b7e0-062c8f5782d1
# ╟─b13d261e-1c1b-41b0-ae86-e20cd2304727
# ╟─f1f04dd4-0447-4e48-9bab-f34058a6a0f1
# ╟─c3757398-060c-412c-997b-58331b2dee04
# ╟─c803ed90-958c-4d0e-84d5-3e8cc9c8fbfe
# ╠═b4ea8364-a582-443c-bdb1-46e75a5c4d4e
# ╟─5018d342-0a51-4275-8c0f-1b5d253c2ec0
# ╟─5afa3aff-af15-429b-a0fb-9a0dfcad740d
# ╟─84a596c1-3a0a-4edd-a7ea-f3fcd8839c2a
# ╟─91b2718a-a921-4bae-99c0-91ec9c7e6479
# ╠═5a0b156e-b0b2-43b8-b15e-bcdb65c154bf
# ╟─a3270d8b-f167-4fb9-bcc5-c62c91b39649
# ╟─6528501f-3bb2-48a7-b24e-d32400033538
# ╟─0d0d6d2f-56b2-4e01-b9a4-6e469eac90ad
# ╟─76cfd5c8-6298-4ae2-9d2d-9e60deb12f49
# ╟─0bf40552-b627-4fa7-b3ba-1d8c7ed4a568


text:
# ENV["JULIA_DEBUG"] = "all"

include("debugger.jl")

try
    VSCodeDebugger.startdebug(VSCodeDebugger.clean_up_ARGS_in_launch_mode())
catch err
    Base.display_error(err, catch_backtrace())
end

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode


text:
const _ENABLED = Ref(true)

"""
Calling `SmartAsserts.set_enabled(false)` will make all future `@smart_assert`
be compiled into no-ops. Hence, simply call this at the beginning of your module
to disable all `@smart_asserts` in your project.
(You might also want to call `SmartAsserts.set_enabled(true)` at the end of your
module to not accidentally turn off others' `@smart_asserts`).
"""
function set_enabled(value::Bool)
    _ENABLED[] = value
end

"""
Evaluate a boolean expression, and if it is false, also try to show additional information
to explain why it's false.

Returns `(value::Bool, reason::Union{Nothing, String})`
"""
macro show_assertion(ex)
    _show_assertion(ex)
end

function _show_assertion(ex)
    ex_type = if Meta.isexpr(ex, :(<:))
        args = ex.args
        :type_assert
    elseif ex.head == :call
        op, pos_args, kw_args = group_function_args(ex.args)
        args = tuple(pos_args..., getindex.(kw_args, 2)...)
        kw_names = getindex.(kw_args, 1)
        :func_call
    elseif ex.head == :comparison
        args = Tuple(e for (i, e) in enumerate(ex.args) if isodd(i))
        ops = Tuple(e for (i, e) in enumerate(ex.args) if iseven(i))
        :comparison
    else
        :other
    end

    ex_q = QuoteNode(ex)
    if ex_type == :other
        :($(esc(ex)), "Condition `$($ex_q)` failed.")
    else
        args_q = Expr(:tuple, QuoteNode.(args)...)
        arg_names = [gensym("arg$i") for i in 1:length(args)]

        cond_ex = if ex_type == :type_assert
            Expr(:(<:), arg_names...)
        elseif ex_type == :func_call
            n_pos_args = length(pos_args)
            Expr(
                :call,
                esc(op),
                Expr(
                    :parameters,
                    (
                        Expr(:kw, k, v) for
                        (k, v) in zip(kw_names, arg_names[(n_pos_args + 1):end])
                    )...,
                ),
                arg_names[1:n_pos_args]...,
            )
        elseif ex_type == :comparison
            comp_args = (
                isodd(i) ? arg_names[1 + i ÷ 2] : ops[i ÷ 2] for i in eachindex(ex.args)
            )
            Expr(:comparison, comp_args...)
        else
            error("This should not be reached.")
        end

        assigns = Expr(
            :block, (Expr(:(=), n, esc(e)) for (n, e) in zip(arg_names, args))...
        )
        args_tuple = Expr(:tuple, arg_names...)

        quote
            $assigns
            if !$(cond_ex)
                eval_string = join(
                    [
                        "\t`$ex` evaluates to $val" for
                        (ex, val) in zip($args_q, $args_tuple)
                        if string(ex) != string(val)
                    ],
                    "\n",
                )
                reason_text = "Condition `$($ex_q)` failed due to:\n" * eval_string
                (false, reason_text)
            else
                (true, nothing)
            end
        end
    end
end

"""
Group the args part of a function call Expr into 3 parts: the operator,
the positional args, and the keyword args.
"""
function group_function_args(args)
    op = args[1]
    has_kw_params = length(args) > 1 && args[2] isa Expr && args[2].head == :parameters

    pos_args = []
    kw_args = Pair{Symbol}[]
    rest_args = has_kw_params ? args[3:end] : args[2:end]
    for a in rest_args
        if a isa Expr && a.head == :kw
            push!(kw_args, (a.args[1] => a.args[2]))
        else
            push!(pos_args, a)
        end
    end
    if has_kw_params
        for a in args[2].args
            pair = if a isa Symbol
                a => a
            else
                @assert a isa Expr
                if a.head == :kw
                    a.args[1] => a.args[2]
                elseif a.head == :.
                    @assert a.args[2] isa QuoteNode
                    a.args[2].value => a
                else
                    error("Unsupported expression in keyword parameters: $a")
                end
            end
            push!(kw_args, pair)
        end
    end
    op, pos_args, kw_args
end

"""
Like @assert, but try to also print out additional information about the arguments.
Note that each argument is only evaluated once, so there is no extra overhead compared
to a normal assert.

Currently, additional information will be returned for the following types of expressions:

- function calls: e.g., `a <= b`, `bool_f(a,k1=...; k2...)`

- type asserts: e.g., `Type1 <: Type2`

- comparison expressions: e.g., `a == b <= c + d <= e`

See also [`SmartAsserts.set_enabled`](@ref) on how to disable `@smart_assert`s
at compile-time.

## Examples
```julia
julia> let a = 1.0, rtol=0.1
           @smart_assert isapprox(a, sin(a), atol=0.05; rtol)
       end
ERROR: AssertionError: Condition `isapprox(a, sin(a), atol = 0.05; rtol)` failed due to:
        `a` evaluates to 1.0
        `sin(a)` evaluates to 0.8414709848078965
        `0.05` evaluates to 0.05
        `rtol` evaluates to 0.1
Stacktrace:
 ...
```
"""
macro smart_assert(ex, msg=nothing)
    if !_ENABLED[]
        return :(nothing)
    end

    has_msg = msg !== nothing
    r_ex = _show_assertion(ex)

    quote
        (result, reason) = $r_ex
        error_msg = $has_msg ? "$($(esc(msg)))\nCaused by $reason" : reason
        if !result
            throw(AssertionError(error_msg))
        end
    end |> Base.remove_linenums!
end


text:
using CellwiseEditingDifferentiation
using Distributions
using Lora

import Base: rand

immutable MvBeta <: ContinuousMultivariateDistribution
  betas::Vector{Beta}
end

Base.rand(d::MvBeta) = [rand(d.betas[i]) for i in 1:length(d.betas)]

ncells = 20

λ = Hyperparameter(:λ)

μ = Data(:μ)
coverage = Data(:coverage)
edited = Data(:edited)
rate = Data(:rate)

function wlogtarget(w::Float64, s::Vector)
  v = s[2]*(1-s[2])*inv_logit(w)
  a, b = beta_pars_from_mv(s[2], v)
  mapreduce(x -> logpdf(Beta(a, b), x), +, s[7])+vpcprior(s[2], v, a, b, s[1])+logdvdw(w, s[2])
end

w = BasicContUnvParameter(:w, logtarget=wlogtarget, nkeys=7)

function psetpdf(p::Vector{Float64}, s::Vector)
  a, b = beta_pars_from_mv(s[2], s[2]*(1-s[2])*inv_logit(s[6]))
  MvBeta([Beta(a+s[4][i], b+s[3][i]-s[4][i]) for i in 1:ncells])
end

p = BasicContMuvParameter(:p, setpdf=psetpdf, nkeys=7)

model = GenericModel(
  [λ, μ, coverage, edited, rate, w, p],
  [λ w; μ w; μ p; coverage p; edited p; w p],
  isindexed=false
)

init = Dict(
  :λ => 10.,
  :μ => 0.5,
  :coverage => fill(10.0, 20),
  :edited => rand([1., 2., 3., 8., 9.], 20),
  :w => 0.
)
init[:rate] = float(init[:edited]./init[:coverage])
init[:p] = rand(psetpdf(fill(NaN, ncells), [init[k] for k in [:λ, :μ, :coverage, :edited, :rate, :w]]))

wjob = BasicMCJob(
  model,
  RAM(),
  BasicMCRange(nsteps=100, burnin=10),
  init;
  pindex=6,
  resetpstate=false,
  verbose=false
)

mcrange = BasicMCRange(nsteps=1100, burnin=100)

outopts = [Dict(:monitor => [:value]), Dict(:monitor => [:value])]

gibbsjob = GibbsJob(model, Dict(:w => wjob), mcrange, init)

@time run(gibbsjob)

output(gibbsjob)

chains = Dict(gibbsjob)

v = map(x -> init[:μ]*(1-init[:μ])*inv_logit(x), chains[:w].value)

writecsv("pchains.csv", chains[:p].value)
writecsv("vchain.csv", v)


text:
# https://www.codewars.com/kata/514b92a657cdc65150000006

module Solution

  triOrQuin(n) = n % 3 == 0 || n % 5 == 0

  function sumofmultiples(n)
    sum(filter(triOrQuin, 1:n-1))
  end

end

max_stars_repo_name: hiljusti/codewars-solutions

max_issues_repo_name: hiljusti/codewars-solutions

max_forks_repo_name: hiljusti/codewars-solutions


text:
using Distributions: MultivariateNormal

function optimaltransport_safe(M::AbstractMatrix, a::AbstractVector, b::AbstractVector;
                        λ::Real=1.0, maxitter::Int=10, ϵ::Real=1e-5)
    Q = exp.(λ * M)
    iter = 0
    while iter < maxitter && maximum(abs.(sum(Q, dims=2)[:] - a)) > ϵ
        iter += 1
        Q = a .* (sum(Q, dims=2) .\ Q)
        Q = Q .* transpose(b) ./ sum(Q, dims=1)
    end
    return Q
end

function laplace_posterior(M, P; γ=1e-2,
                maxitter::Int=20, tol=1e-4)
    a, b = marginals(P)
    Q = M -> Q = optimaltransport_safe(exp.(M), a, b)
    l = M -> relative_entropy(P, Q(M)) + γ * sum(M .* M)
    hess = Zygote.hessian(l, M)
    return MultivariateNormal(vec(exp.(M)), Symmetric(hess+tol*I)^-1)
end

#TODO: test and explore this.

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
const zerosVector3f = [0f0,0,0]
const onesVector3f = [1f0,1,1]

#-------------------------------------------------------------------------------------------------
# Base implementations for StaticArray

Base.isless{N,T<:Number}(a::StaticArray{Tuple{N},T,1}, b::StaticArray{Tuple{N},T,1}) = begin
  for i=N:-1:1
    if a[i] < b[i] return true
    elseif a[i] > b[i] return false
    end
  end
  false
end

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

VectorType{N, T} = FieldVector{N, T}

struct Vector1{T} <: VectorType{1, T}
  x::T
end

struct Vector2{T} <: VectorType{2, T}
  x::T
  y::T
end

struct Vector3{T} <: VectorType{3, T}
  x::T
  y::T
  z::T
end

struct Vector4{T} <: VectorType{4, T}
  x::T
  y::T
  z::T
  w::T
end

struct Vector{N, T}
  Vector{N,T}(x) where {N,T} = (N==4 ? Vector4{T} : N==3 ? Vector3{T} : N==2 ? Vector2{T} : VectorType{N, T})(x)
end

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

# set
Vec2{T}(x::T) = Vector2{T}(x,x)
Vec3{T}(x::T) = Vector2{T}(x,x,x)
Vec4{T}(x::T) = Vector2{T}(x,x,x,x)

# copy
Vec2{T}(v::Vector2{T}) = Vector2{T}(v.x,v.y)
Vec3{T}(v::Vector3{T}) = Vector3{T}(v.x,v.y,v.z)
Vec4{T}(v::Vector4{T}) = Vector4{T}(v.x,v.y,v.z,v.w)

# reduce
Vec2{T}(v::Union{Vector3{T}, Vector4{T}}) = Vector2{T}(v.x,v.y)
Vec3{T}(v::Vector4{T}) = Vector3{T}(v.x,v.y,v.z)

# extend
Vec3{T<:Number}(v::Vector2{T}) = Vector3{T}(v.x,v.y,0)
Vec4{T<:Number}(v::Vector2{T}) = Vector4{T}(v.x,v.y,0,0)
Vec4{T<:Number}(v::Vector3{T}) = Vector4{T}(v.x,v.y,v.z,0)

# combine
Vec3{T<:Number}(v::Vector2{T},z::Number) = Vector3{T}(v.x,v.y,z)
Vec4{T<:Number}(v::Vector3{T},w::Number) = Vector4{T}(v.x,v.y,v.z,w)
Vec4{T}(v1::Vector2{T},v2::Vector2{T}) = Vector4{T}(v1.x,v1.y,v2.x,v2.y)

# operate
Base.:+{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:+.(v1,v2))
Base.:-{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:-.(v1,v2))
Base.:*{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:*.(v1,v2))
Base.:/{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:/.(v1,v2))
Base.:\{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:\.(v1,v2))
Base.:^{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:^.(v1,v2))
Base.:%{N, T<:Number}(v1::VectorType{N, T}, v2::VectorType{N, T}) = Vector{N, T}(Base.:%.(v1,v2))

#normalize
#normalize{N, T<:Number}(v::VectorType{N, T}) = Vector{N, T}(normalize(v))

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

const Vector2f = Vector2{Float32}
const Vector3f = Vector3{Float32}
const Vector4f = Vector4{Float32}

const Vector2i = Vector2{Int32}
const Vector3i = Vector3{Int32}
const Vector4i = Vector4{Int32}

const Vector2u = Vector2{UInt32}
const Vector3u = Vector3{UInt32}
const Vector4u = Vector4{UInt32}

const Vector2d = Vector2{Integer}
const Vector3d = Vector3{Integer}
const Vector4d = Vector4{Integer}

const Vec2f = Vector2f
const Vec3f = Vector3f
const Vec4f = Vector4f

const Vec2i = Vector2i
const Vec3i = Vector3i
const Vec4i = Vector4i

const Vec2u = Vector2u
const Vec3u = Vector3u
const Vec4u = Vector4u

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

Vec2f() = zeros(Vec2f)
Vec3f() = zeros(Vec3f)
Vec4f() = zeros(Vec4f)

Vec2i() = zeros(Vec2i)
Vec3i() = zeros(Vec3i)
Vec4i() = zeros(Vec4i)

Vec2u() = zeros(Vec2u)
Vec3u() = zeros(Vec3u)
Vec4u() = zeros(Vec4u)

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

Vec3f{T<:Number}(v::Vector2{T},z::Number) = Vec3(v,z)
Vec4f{T<:Number}(v::Vector3{T},w::Number) = Vec4(v,w)
Vec4f{T<:Number}(v1::Vector2{T},v2::Vector2{T}) = Vec4(v1,v2)

Vec3i{T<:Number}(v::Vector2{T},z::Number) = Vec3(v,z)
Vec4i{T<:Number}(v::Vector3{T},w::Number) = Vec4(v,w)
Vec4i{T<:Number}(v1::Vector2{T},v2::Vector2{T}) =Vec4(v1,v2)

Vec3u{T<:Number}(v::Vector2{T},z::Number) = Vec3(v,z)
Vec4u{T<:Number}(v::Vector3{T},w::Number) = Vec4(v,w)
Vec4u{T<:Number}(v1::Vector2{T},v2::Vector2{T}) = Vec4(v1,v2)

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
@testset "MST" begin
    using STMO
    using STMO.MST

    edges = [(1, :A, :D),
            (2, :A, :B),
            (3, :D, :C),
            (4, :B, :D),
            (5, :B, :C)]

   vertices = [:A, :B, :C, :D]

   mst_prim, cprim = prim(vertices, edges, :A)
   @test cprim == 6
   @test (4, :B, :D) ∉ mst_prim && (4, :D, :B) ∉ mst_prim

   @test prim(vertices, edges)[2] == 6

   mst_kruskal, ckruskal = kruskal(vertices, edges)
   @test ckruskal == 6
   @test (4, :B, :D) ∉ mst_kruskal && (4, :D, :B) ∉ mst_kruskal
end


text:
using Test
using MIPVerify
using MIPVerify: LInfNormBoundedPerturbationFamily
using MIPVerify: get_example_network_params, read_datasets, get_image, get_label
@isdefined(TestHelpers) || include("../../../TestHelpers.jl")

@testset "MNIST.WK17a_linf0.1_authors.jl" begin
    nn = get_example_network_params("MNIST.WK17a_linf0.1_authors")
    mnist = read_datasets("mnist")

    test_cases = [(1, NaN), (2, NaN), (9, 0.0940014207), (248, 0)]

    for test_case in test_cases
        (index, expected_objective_value) = test_case
        @testset "Sample $index (1-indexed) with expected objective value $expected_objective_value" begin
            input = get_image(mnist.test.images, index)
            label = get_label(mnist.test.labels, index)
            TestHelpers.test_find_adversarial_example(
                nn,
                input,
                label + 1,
                LInfNormBoundedPerturbationFamily(0.1),
                Inf,
                0,
                expected_objective_value,
                invert_target_selection = true,
            )
        end
    end
end


text:
"""
This file contains the code for the work by

Y. Chandak, S. Shankar, and P. S. Thomas.
High Confidence Off-Policy (or Counterfactual) Variance Estimation.
In Proceedings of the Thirty-Fifth AAAI Conference on Artificial Intelligence, 2021

Code repo: https://github.com/yashchandak/UnO
"""

using Statistics
using Random
using Bootstrap

include("./HCOPE.jl")

######################
# Variance Estimator
######################

function CDIS_fulltraj(data, props, ξG2=0.0)
    """
    Estimator for E[ρG^2] using full-traj IS

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    ξG2: Control variate to ensure G^2 is always non-positive or non-negative
    """

    # Meta-properties
    n, horizon, d = size(data)
    γ = props["gamma"]
    total = fill(0.0, n)

    for i in 1:1:n
        returns = 0.0
        ρ = 1.0
        for j in 1:1:horizon
            ρ = ρ * data[i,j,1]
            returns += γ^(j-1) * data[i,j,2]
        end;
        total[i] = ρ * (returns^2 - ξG2)
    end;
    return total
end;

function CDIS(data, props, ξr2=0.0)
    """
    Estimator for E[ρG^2] using
    coupled decision importance sampling (CDIS) [Thm 3 by Chandak et al. (2021)]

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    ξr2: Control variate to ensure r^2 is always non-positive or non-negative
    """

    n, horizon, d = size(data)
    γ = props["gamma"]
    total = fill(0.0, n)

    for i in 1:1:n
        val = 0.0
        ρ_outer = 1.0

        ############
        # Alternate form
        # ∑_{j=0}^{T} ∑_{k=0}^T ρ(0,max(j, k)) γ^{j+k} R_j R_k
        # = ∑_{j=0}^T [ ρ(0, j)γ^{2j} R_j^2  + 2 ∑_{k=j+1}^T ρ(0,k) γ^{j+k} R_j R_k ]
        ############

        for j in 1:1:horizon
            ρ_outer = ρ_outer * data[i,j,1]
            R_j = data[i,j,2]
            val = val + ρ_outer * γ^(2*j-2) * (R_j^2 - ξr2)

            ρ_inner = ρ_outer
            for k in j+1:1:horizon
                ρ_inner = ρ_inner * data[i,k,1]
                R_k = data[i,k,2]
                val = val + 2 * ρ_inner * γ^(j+k-2) * (R_j * R_k - ξr2)
            end;
        end;
        total[i] = val
    end;
    return total
end;

function var_estimator_naive1(data, props)
    """
    This is the ̂σ^{!!} estimator

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    """

    n, horizon, d = size(data)
    IS_returns = OPE_PDIS(data, props)

    return var(IS_returns)
end;

function var_estimator_naive2(data, props)
    """
    This is the ̂σ^{!} estimator

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    """

    n, horizon, d = size(data)
    IS_returns = OPE_PDIS(data, props)

    μ = mean(IS_returns)

    T, γ = props["horizon"], props["gamma"]
    total = fill(0.0, n)
    for i in 1:1:n
        returns = 0.0
        ρ = 1.0
        for j in 1:1:horizon
            ρ = ρ * data[i,j,1]
            returns += γ^(j-1) * data[i,j,2]
        end;
        total[i] = ρ * (returns - μ)^2
    end;
    return sum(total)/(n-1)

end;

function var_estimator(data, props)
    """
    This is the ̂σ_n estimator

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    """
    n, horizon, d = size(data)

    # Per/coupled decision variants
    term1 = mean(CDIS(data, props))                         # E[ρG^2]
    term2 = mean(OPE_PDIS(data[1:(n÷2), :, :], props))      # E[ρG] using first half of data
    term3 = mean(OPE_PDIS(data[(n÷2)+1:n, :, :], props))    # E[ρG] using second half of data

    return term1 - term2*term3                              # E[ρG^2] - E[ρG]E[ρG]
end;

function var_estimator_wo_CDIS(data, props)
    n, horizon, d = size(data)

    # No per-step CDIS
    term1 = mean(CDIS_fulltraj(data, props))
    term2 = mean(OPE_PDIS(data[1:(n÷2), :, :], props))
    term3 = mean(OPE_PDIS(data[(n÷2)+1:n, :, :], props))

    return term1 - term2*term3
end;

##########################################
#  High confidence Bounds
#######################################

function global_bounds(val, props)
    """
    Intersection of the proposed upper/lower values with
    lowest possible value and the higheest possible value (from Popoviciu's ineqality)

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    """

    rmin, rmax = props["rmin"], props["rmax"]
    gmin, gmax =  props["Gmin"], props["Gmax"]
    T, γ = props["horizon"], props["gamma"]

    geom_sum = (1 - γ^T)/(1 - γ)
    min_g = max(geom_sum*rmin, gmin)   # Maximum of the two lower possibilities
    max_g = min(geom_sum*rmax, gmax)   # Minimum of thw two upper possibilties

    max_var = ((max_g - min_g)^2) /4.0  # Maximum possible variance
    min_var = 0.0                       # Minimum possible variance

    return clamp.(val, min_var, max_var)
end;

function CDIS_lower(data, props, δ=0.05)
    """
    Lower bound on E[ρG^2]
    This uses the bound by Thomas et al. (2015) as a sub-routine

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    n, horizon, d = size(data)

    rmin, rmax = props["rmin"], props["rmax"]
    gmin, gmax =  props["Gmin"], props["Gmax"]
    T, γ = props["horizon"], props["gamma"]

    geom_sum = (1 - γ^T)/(1 - γ)
    min_g = max(geom_sum*rmin, gmin)       # Maximum of the two lower possibilities
    max_g = min(geom_sum*rmax, gmax)       # Minimum of thw two upper possibilties

    # minimum G^2 value possible
    min_g2 = 0.0
    if !(min_g < 0 && max_g > 0)
        # If the sign of min_g and max_g are the same
        min(min_g^2, max_g^2)
    end;

    # Get the CDIS returns
    # Note: ρG^2 is always positive
    # Therefore, no need to use control variates to subtract anything.
    estimates = CDIS(data, props, 0.0)

    if any(x -> x<0, estimates)
        println("CDIS Bound will not work, items in array less than 0")
    end;

    function lower_bound(D, c, n_post)
        Y = clamp.(D, 0.0, c)
        term1 = mean(Y)
        term2 = (7*c*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        lb = term1 - term2 - term3                      # Thm1 from Thomas (2015)
        # Note: No control variate shift is done here, as nothing was subtracted in CDIS

        return lb
    end;

    #########################################################
    # Prepare data splits into dpre and dpost
    # Size based on the split-value suggested by Thomas (2015)
    n_pre = Int(floor(0.05 * n))
    n_post = n - n_pre

    perm_idxs = shuffle(collect(1:n))
    D_pre = estimates[perm_idxs[1:n_pre]]
    D_post = estimates[perm_idxs[n_pre + 1:n]]

    ##########################################################
    # Search for best truncation point using basic random search in 1D
    min_obs, max_obs = minimum(D_pre), maximum(D_pre)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_max, c_max = -Inf, cs[1]

    # Search for c that maximizes the lower bound
    for c in cs
        val = lower_bound(D_pre, c, n_post)
        if val > v_max
            v_max, c_max = val, c
        end;
    end;

    ########################################################
    # Bound by Thomas(2015)
    lb = lower_bound(D_post, c_max, n_post)

    # If lower bound is lower than minimum G^2 possible, then clip.
    lb = max(lb, min_g2)

    return lb
end;

function CDIS_upper(data, props, δ=0.05)
    """
    Upper bound on E[ρG^2]
    This uses the bound by Thomas et al. (2015) as a sub-routine

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    n, horizon, d = size(data)
    rmin, rmax = props["rmin"], props["rmax"]
    gmin, gmax =  props["Gmin"], props["Gmax"]
    T, γ = props["horizon"], props["gamma"]

    geom_sum = (1 - γ^T)/(1 - γ)
    ξr2 = max(rmin^2, rmax^2)               # maximum possible r^2 value

    min_g = max(geom_sum*rmin, gmin)        # Maximum of the two lower possibilities
    max_g = min(geom_sum*rmax, gmax)        # Minimum of thw two upper possibilties
    ξG2 = max(min_g^2, max_g^2)             # maximum G^2 value possible

    # Get the importance sampled returns ** with control variate **
    # Choose the control varaite smartly
    # if the full trajectory control is smaller, then use CDIS_fulltraj
    # Otherwise use CDIS with coupled decision control
    # Note: ξr2*(c^2) is never smaller than ξG2; it can only be equal to it at best.
    control = 0.0
    if ξG2 < ξr2*(geom_sum^2)
        estimates = CDIS_fulltraj(data, props, ξG2)
        control = ξG2
    else
        estimates = CDIS(data, props, ξr2)
        control = ξr2*(geom_sum^2)
    end;

    if any(x -> x>0, estimates)
        println("CDIS Bound will not work, items in array more than 0")
    end;

    function upper_bound(D, c, n_post)
        Y = clamp.(D, c, 0.0)
        term1 = mean(Y)
        term2 = (7*abs(c)*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        ub = term1 + term2 + term3                      # Thm1 from Thomas (2015)
        ub = ub + control      # add back Gmax^2 to the bound
        return ub
    end;

    #########################################################
    # Prepare data splits into dpre and dpost
    # Size based on the split-value suggested by Thomas (2015)
    n_pre = Int(floor(0.05 * n))
    n_post = n - n_pre

    perm_idxs = shuffle(collect(1:n))
    D_pre = estimates[perm_idxs[1:n_pre]]
    D_post = estimates[perm_idxs[n_pre + 1:n]]

    ##########################################################
    # Search for best truncation point using basic random search in 1D
    min_obs, max_obs = minimum(D_pre), maximum(D_pre)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_min, c_min = +Inf, cs[1]

    # Search for c that minimizes the upper bound
    for c in cs
        val = upper_bound(D_pre, c, n_post)
        if val < v_min
            v_min, c_min = val, c
        end;
    end;

    ########################################################
    # Bound by Thomas(2015)
    ub = upper_bound(D_post, c_min, n_post)

    # If upper bound is higher than maximum G^2 possible, then clip.
    ub = min(ub, ξG2)

    return ub
end;

function HCOVE_upper(data, props, δ=0.05)
    """
    Upper bound on the variance (Equation 6)

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    n, horizon, d = size(data)

    # δ/2 upper bound on E[ρG^2]
    term1 = CDIS_upper(data, props, δ/2)

    # δ/2 lower bound on E[ρG]^2
    # This requires two-sided bound on E[ρG]
    upper = CI_upper(data, props, δ/4)
    lower = CI_lower(data, props, δ/4)

    # Obtain lower bound on E[ρG]^2 from E[ρG] using bound propagation
    if lower <= 0.0 && upper >= 0.0
        term2 = 0.0
    else
        term2 = min(lower^2, upper^2)
    end;

    ub = term1 - term2

    # Bound by Popoviciu
    ub = global_bounds(ub, props)

    return ub
end;

function HCOVE_lower(data, props, δ=0.05)
    """
    Lower bound on the variance (Equation 5)

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    n, horizon, d = size(data)

    # δ/2 lower bound on E[ρG^2]
    term1 = CDIS_lower(data, props, δ/2)

    # δ/2 upper bound on E[ρG]^2
    # This requires two-sided bound on E[ρG]
    upper = CI_upper(data, props, δ/4)
    lower = CI_lower(data, props, δ/4)

    # Obtain upper bound on E[ρG]^2 from E[ρG] using bound propagation
    term2 = max(lower^2, upper^2)

    lb = term1 - term2

    # Bound by Popoviciu
    lb = global_bounds(lb, props)

    return lb
end;

function HCOVE_Bootstrap(data, props)
    """
    Bootstrap based bounds for the variance

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate (0.05)
    """

    full_size, _, _ = size(data)
    idxs = collect(1:full_size)

    bs = bootstrap(x -> var_estimator(data[x, :, :], props), idxs, BasicSampling(100))   # Bootstrap mean's distribution
    # bci = confint(bs, BCaConfInt(0.95))               # Bootstrap CI on the mean
    bci = confint(bs, PercentileConfInt(0.95))          # Bootstrap CI on the mean

    estimate, lb, ub = bci[1]
    lb, ub = global_bounds([lb, ub], props)             # Intersect with the global possibility
    return [estimate, lb, ub]
end;

max_stars_repo_name: yashchandak/UnO

max_issues_repo_name: yashchandak/UnO

max_forks_repo_name: yashchandak/UnO


text:
export remove_response!

function remove_response!(data::AbstractArray, freq_response::AbstractArray)
    data[1] = 0.0 + 0.0 * im # set DC frequency to zero
    for i in 2:length(data)
        data[i] /= freq_response[i]
    end
end


text:
module Tst
    include("tmp.jl")
    using .Tmp

    # add test code below
    a = get_baseline_points()
    println(size(a))

end

max_stars_repo_name: mtmcrobie/ultra-proc

max_issues_repo_name: mtmcrobie/ultra-proc

max_forks_repo_name: mtmcrobie/ultra-proc


text:
# Update position of walker
function update(w::Walker1D, s::Step1D;width::Int=1,length::Int=1)
    w+s
end

function update(w::Walker2D, s::Step2D,width::Int=1,length::Int=1)
    w+s
end

# Update position of walker on periodic XY domain discretised such that
# width and length are integers indices
function updateperiod(w::Walker2D, s::Step2D,width::Int,length::Int)
    pos = update(w, s)
    w₊ = (1/2)*width
    w₋ = -(1/2)*width
    l₊ = (1/2)*length
    l₋ = -(1/2)*length

    # Check x boundary
    if pos.x > w₊
        x = w₋
    elseif pos.x < w₋
        x = w₊
    else
        x = pos.x
    end

    # Check y boundary
    if pos.y > l₊
        y = l₋
    elseif pos.y < l₋
        y =  l₊
    else
        y = pos.y
    end

    return Walker2D(Int(x),Int(y))
end

# Update where domain is periodic in one direction and
# midline reflective in the other direction
# Rectangle version of a capsule
# Not ready

function updatecell(w::Walker2D, s::Step2D,width::Int,length::Int)
    # x width say 0 to 100
    pos = update(w, s)
    w₊ = (1/2)*width
    w₋ = -(1/2)*width
    l₊ = (1/2)*length
    l₋ = -(1/2)*length

    # Check x boundary
    if pos.x > w₊
        x = w₋
    elseif pos.x < w₋
        x = w₊
    else
        x = pos.x
    end

    # Check y boundary
    if pos.y > l₊ # north
        x = rand(1:width)
        y = w.y
    elseif pos.y < l₋ #
        x = rand(1:width)
        y = w.y
    else
        y = pos.y
    end

    Walker2D(Int(x), Int(y))
end

# Update position of walker on periodic XY domain discretised such that
# width and length are integers indices
function updatenoflux(w::Walker2D, s::Step2D,sₙ::Step2D, sₛ::Step2D,width::Int,length::Int)
    pos = update(w, s)
    w₊ = (1/2)*width
    w₋ = -(1/2)*width
    l₊ = (1/2)*length
    l₋ = -(1/2)*length

    # Check x boundary
    if pos.x > w₊
        x = w₋
    elseif pos.x < w₋
        x = w₊
    else
        x = pos.x
    end

    # Check y boundary
    if pos.y > l₊ # north
        y = update(w, sₙ).y
    elseif pos.y < l₋ #
        y = update(w, sₛ).y
    else
        y = pos.y
    end

    return Walker2D(Int(x),Int(y))
end

function updatenoflux(state::Int64,w::Walker2D, rate::Float64,width::Int,length::Int)
    pos = update(w, stepper(state,stepping,rate))
    w₊ = (1/2)*width
    w₋ = -(1/2)*width
    l₊ = (1/2)*length
    l₋ = -(1/2)*length

    # Check x boundary
    if pos.x > w₊
        x = w₋
    elseif pos.x < w₋
        x = w₊
    else
        x = pos.x
    end

    # Check y boundary
    if pos.y > l₊ # north
        y = update(w, stepper(state,stepping_no_north,rate)).y
    elseif pos.y < l₋ #
        y = update(w, stepper(state,stepping_no_south,rate)).y
    else
        y = pos.y
    end

    return Walker2D(Int(x),Int(y))
end


text:
function is_undirected(n::Int, cptr::Vector, rval::Vector)
    for col = OneTo(n)
        l1 = cptr[col+1]-1
        for i = 0 : (l1 - cptr[col])
            if rval[l1-i] <= col
                break
            end
            return false
        end
    end
    return true
end

@noinline function _check_type(n::Integer, T::Type)
    if 0 ≤ n && (!isbitstype(Ti) || n ≤ typemax(Ti))
        throw(ArgumentError("number of vertices (nv = $n) does not fit in eltype = $(T)"))
    end
    return nothing
end

# stronger check for debugging purposes
# && all(issorted(@view rowval[colptr[i]:colptr[i+1]-1]) for i=1:n)
_goodbuffers(g::AdjacencyMap) = _goodbuffers(nvertices(g), getcolptr(g), rowvals(g), metadata(g))
function _goodbuffers(n, cptr, rval, ::NoMetadata)
    if !(length(cptr) == n + 1 && cptr[end] - 1 == length(rval))
        throw(ArgumentError("Illegal buffers for AdjacencyMap construction $n $colptr $rowval $nzval"))
    end
    return nothing
end
function _goodbuffers(n, cptr, rval, m)
    if !(length(cptr) == n + 1 && cptr[end] - 1 == length(rval) == length(m))
        throw(ArgumentError("Illegal buffers for AdjacencyMap construction $n $colptr $rowval $nzval"))
    end
    return nothing
end

function _check_length(rowstr, rowval, minlen, Ti)
    throwmin(len, minlen, rowstr) = throw(ArgumentError("$len == length($rowstr) < $minlen"))
    throwmax(len, max, rowstr) = throw(ArgumentError("$len == length($rowstr) >= $max"))

    len = length(rowval)
    len >= minlen || throwmin(len, minlen, rowstr)
    !isbitstype(Ti) || len < typemax(Ti) || throwmax(len, typemax(Ti), rowstr)
    return nothing
end

_check_meta_length(ne, ::NoMetadata) = nothing
function _check_meta_length(ne, m)
    if ne != length(m)
        throw(ArgumentError("""
        cannot assign AdjacencyMap list of metadata with length different than number of edges:
        number of edges = $(Int(ne))
        length metadata list = $(length(m))
        """))
    end
    return nothing
end
# This is a bad but it's also the only thing the stops me from writing the print code for
# MetaAdjView here.
SparseArrays._checkbuffers(m::MetaAdjView) = _goodbuffers(parent(m))


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module Collage

export collage, collage_compact

using Compat
import Compat.ASCIIString
using ..DefineSave
using ..Visualize
using ..Visualize.TikzUtils
using ..SaveHelpers
using RLESUtils, FileUtils

using TikzPictures
import PGFPlots
import PGFPlots: Plots, Axis, GroupPlot

"""
Make a collage out of the vertical profile plots of the first N json.gz files
Expects clusterdirs to be in cluster number order
"""
function collage{T<:AbstractString}(outfileroot::T, clusterdirs::Vector{T},
    ncols::Int64, nrows::Int64=1)

    td = TikzDocument()
    enc_ids = Int64[]
    N = nrows * ncols
    for (i, dir) in enumerate(clusterdirs)
        fs = readdirGZs(dir)
        NN = min(N, length(fs))
        fs = fs[1:NN]
        g = GroupPlot(ncols, nrows, groupStyle="horizontal sep=0.5cm, vertical sep=0.5cm")
        empty!(enc_ids)
        for f in fs
            sav = trajLoad(f)
            ax = pgfplot_alt(sav)
            push!(g, ax)
            push!(enc_ids, sv_encounter_id(sav)[1])
        end
        for j = 1:N-NN #fill remaining with empty axis for scaling
            push!(g, Axis())
        end
        tp = PGFPlots.plot(g)
        use_geometry_package!(tp, landscape=true)
        use_aircraftshapes_package!(tp)
        @compat cap = "Cluster $i: encounters=$(join(string.(enc_ids),","))"
        push!(td, tp; caption=cap)
    end
    TikzPictures.save(TEX(outfileroot), td)
    wrap_tikzpicture!("$(outfileroot).tex", "\\resizebox{\\textwidth}{!}{%\n", "\n}")
    lualatex_compile("$(outfileroot).tex")
    td
end

"""
Make a collage out of the vertical profile plots of the first N json.gz files
compact version
"""
function collage_compact{T<:AbstractString}(outfileroot::T, clusterdirs::Vector{T},
    N::Int64)

    td = TikzDocument()
    nclusters = length(clusterdirs)
    g = GroupPlot(N, nclusters, groupStyle="horizontal sep=1.75cm, vertical sep=1.75cm")
    enc_ids = Int64[]
    caps = ASCIIString[]
    for (i, dir) in enumerate(clusterdirs)
        fs = readdirGZs(dir)
        NN = min(N, length(fs))
        fs = fs[1:NN]
        empty!(enc_ids)
        for f in fs
            sav = trajLoad(f)
            ax = pgfplot_alt(sav)
            push!(g, ax)
            push!(enc_ids, sv_encounter_id(sav)[1])
        end
        for j = 1:N-NN #fill remaining with empty axis for scaling
            push!(g, Axis())
        end
        @compat push!(caps, "(Cluster $i: encounters=$(join(string.(enc_ids),",")))")
    end
    tp = PGFPlots.plot(g)
    use_geometry_package!(tp, landscape=true)
    use_aircraftshapes_package!(tp)
    push!(td, tp)
    TikzPictures.save(TEX(outfileroot), td)
    wrap_tikzpicture!("$(outfileroot).tex", "\\resizebox{!}{0.95\\textheight}{%\n", "\n}")
    lualatex_compile("$(outfileroot).tex")

    #output encounter numbers to text file
    textfile("$(outfileroot).txt", caps...)
    td
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS