text:
"""
    axes2reduce(z, x) -> axes::Vector{Int}
axes need to be reduced, `z` and `x` comes from `z = broadcast(::typeof(+-*/...), x, y)`\n

# Example
    axes2reduce(rand(3,4,5),rand(1,4))    -> (1,3)
    axes2reduce(rand(3,4,5),rand(1,4,1))  -> (1,3)
"""
function axes2reduce(z, x)
    a = Int[]
    for i = 1:ndims(z)
        if size(x, i) == 1
            push!(a, i)
        end
    end
    return a
end

"""
    unbcast(δx::AbstractArray, x::AbstractArray) -> ∇x

reduced `δx` to `∇x` according to shape difference from `x` and `δx`

# Params
`x`  : comes from `z = broadcast(::typeof(+-*/...), x, y)`\n
`δx` : unreduced gradient, i.e. `δx = δz .* ∂z/∂x`\n
`∇x` : reduced gradient, i.e. ⤓⤓⤓\n
       Δx = sum(δx, dims=axes2reduce(δx, x)) # reduced but still has redundant dimensions\n
       ∇x = reshape(Δx, size(x))
"""
function unbcast(δx::AbstractArray, x::AbstractArray)
    if size(δx) == size(x)
        return δx
    elseif length(δx) == length(x)
        return reshape(δx, size(x))
    else
        Δx = sum(δx, dims=axes2reduce(δx,x))
        return reshape(Δx, size(x))
    end
end

# z = x .+ y
function Base.Broadcast.broadcasted(::typeof(+), x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) .+ ᵛ(y), backprop)
    if backprop
        function DotAddBackward()
            if need2computeδ!(x)
                δx = δ(z)
                δ(x) .+= unbcast(δx, ᵛ(x))
            end
            if need2computeδ!(y)
                δy = δ(z)
                δ(y) .+= unbcast(δy, ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, DotAddBackward)
    end
    return z
end

# z = x .- y
function Base.Broadcast.broadcasted(::typeof(-), x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) .- ᵛ(y), backprop)
    if backprop
        function DotMinusBackward()
            if need2computeδ!(x)
                δx = δ(z)
                δ(x) .+= unbcast(δx, ᵛ(x))
            end
            if need2computeδ!(y)
                δy = - δ(z)
                δ(y) .+= unbcast(δy, ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, DotMinusBackward)
    end
    return z
end

# z = x .* y
function Base.Broadcast.broadcasted(::typeof(*), x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) .* ᵛ(y), backprop)
    if backprop
        function DotMulBackward()
            if need2computeδ!(x)
                δx = δ(z) .* ᵛ(y)
                δ(x) .+= unbcast(δx, ᵛ(x))
            end
            if need2computeδ!(y)
                δy = δ(z) .* ᵛ(x)
                δ(y) .+= unbcast(δy, ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, DotMulBackward)
    end
    return z
end

# z = x ./ y
function Base.Broadcast.broadcasted(::typeof(/), x::Variable{T1}, y::Variable{T2}) where {T1,T2}
    @assert T1 <: T2 || T1 >: T2
    T = T1 <: T2 ? T1 : T2
    backprop = (x.backprop || y.backprop)
    z = Variable{T}(ᵛ(x) ./ ᵛ(y), backprop)
    if backprop
        function DotDivBackward()
            δx = δ(z) ./ ᵛ(y)
            if need2computeδ!(x)
                δ(x) .+= unbcast(δx, ᵛ(x))
            end
            if need2computeδ!(y)
                δy = - δx .* ᵛ(z)
                δ(y) .+= unbcast(δy, ᵛ(y))
            end
            ifNotKeepδThenFreeδ!(z);
        end
        push!(graph.backward, DotDivBackward)
    end
    return z
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# see below for imfilter docstring

# Step 1: if necessary, determine the output's element type
@inline function imfilter(img::AbstractArray, kernel, args...)
    imfilter(filter_type(img, kernel), img, kernel, args...)
end

# Step 2: if necessary, put the kernel into cannonical (factored) form
@inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, args...)
    imfilter(T, img, factorkernel(kernel), args...)
end

# Step 3: if necessary, fill in the default border
function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, args...)
    imfilter(T, img, kernel, "replicate", args...)
end

function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::AbstractString, args...)
    if border ∈ valid_borders
        return imfilter(T, img, kernel, Pad(Symbol(border)), args...)
    elseif border == "inner"
        error("specifying Inner as a string is deprecated, use `imfilter(img, kern, Inner())` instead")
    else
        throw(ArgumentError("$border not a recognized border"))
    end
end

# Step 4: if necessary, allocate the ouput
@inline function imfilter{T}(::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny, args...)
    imfilter!(allocate_output(T, img, kernel, border), img, kernel, border, args...)
end

# Now do the same steps for the case where the user supplies a Resource
@inline function imfilter(r::AbstractResource, img::AbstractArray, kernel, args...)
    imfilter(r, filter_type(img, kernel), img, kernel, args...)
end

@inline function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ArrayLike, args...)
    imfilter(r, T, img, factorkernel(kernel), args...)
end

# For steps 3 & 4, we make args... explicit as a means to prevent
# specifying both r and an algorithm
function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel)
    imfilter(r, T, img, kernel, Pad(:replicate))  # supply the default border
end

function imfilter{T}(r::AbstractResource, ::Type{T}, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny)
    imfilter!(r, allocate_output(T, img, kernel, border), img, kernel, border)
end

"""
    imfilter([T], img, kernel, [border="replicate"], [alg]) --> imgfilt
    imfilter([r], img, kernel, [border="replicate"], [alg]) --> imgfilt
    imfilter(r, T, img, kernel, [border="replicate"], [alg]) --> imgfilt

Filter an array `img` with kernel `kernel` by computing their
correlation.

`kernel[0,0,..]` corresponds to the origin (zero displacement) of the
kernel; you can use `centered` to place the origin at the array
center, or use the OffsetArrays package to set `kernel`'s indices
manually. For example, to filter with a random *centered* 3x3 kernel,
you could use either of the following:

    kernel = centered(rand(3,3))
    kernel = OffsetArray(rand(3,3), -1:1, -1:1)

`kernel` can be specified as an array or as a "factored kernel," a
tuple `(filt1, filt2, ...)` of filters to apply along each axis of the
image. In cases where you know your kernel is separable, this format
can speed processing.  Each of these should have the same
dimensionality as the image itself, and be shaped in a manner that
indicates the filtering axis, e.g., a 3x1 filter for filtering the
first dimension and a 1x3 filter for filtering the second
dimension. In two dimensions, any kernel passed as a single matrix is
checked for separability; if you want to eliminate that check, pass
the kernel as a single-element tuple, `(kernel,)`.

Optionally specify the `border`, as one of `Fill(value)`,
`"replicate"`, `"circular"`, `"symmetric"`, `"reflect"`, `NA()`, or
`Inner()`. The default is `"replicate"`. These choices specify the
boundary conditions, and therefore affect the result at the edges of
the image. See `padarray` for more information.

`alg` allows you to choose the particular algorithm: `FIR()` (finite
impulse response, aka traditional digital filtering) or `FFT()`
(Fourier-based filtering). If no choice is specified, one will be
chosen based on the size of the image and kernel in a way that strives
to deliver good performance. Alternatively you can use a custom filter
type, like `IIRGaussian`.

Optionally, you can control the element type of the output image by
passing in a type `T` as the first argument.

You can also dispatch to different implementations by passing in a
resource `r` as defined by the ComputationalResources package.  For
example,

    imfilter(ArrayFire(), img, kernel)

would request that the computation be performed on the GPU using the
ArrayFire libraries.

See also: imfilter!, centered, padarray, Pad, Fill, Inner, IIRGaussian.
"""
imfilter

# see below for imfilter! docstring

# imfilter! can be called directly, so we take steps 2&3 here too. We
# have to be a little more cautious to make sure that later methods
# don't inadvertently call back to these: in methods that take an
# AbstractResource argument, exclude `NoPad()` as a border option.
function imfilter!(out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, args...)
    imfilter!(out, img, factorkernel(kernel), args...)
end

function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian})
    imfilter!(r, out, img, factorkernel(kernel))
end

function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, border::AbstractString)
    imfilter!(r, out, img, kernel, Pad(Symbol(border)))
end

function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::Union{ArrayLike,Laplacian}, border::BorderSpec)
    imfilter!(r, out, img, factorkernel(kernel), border)
end

function imfilter!(out::AbstractArray, img::AbstractArray, kernel::ProcessedKernel)
    imfilter!(out, img, kernel, Pad(:replicate))
end

function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::ProcessedKernel)
    imfilter!(r, out, img, kernel, Pad(:replicate))
end

function imfilter!(r::AbstractResource, out::AbstractArray, img::AbstractArray, kernel::ProcessedKernel, border::AbstractString)
    imfilter!(r, out, img, kernel, Pad(Symbol(border)))
end

# Step 5: if necessary, pick an algorithm
function imfilter!(out::AbstractArray, img::AbstractArray, kernel::ProcessedKernel, border::AbstractString)
    imfilter!(out, img, kernel, Pad(Symbol(border)))
end

function imfilter!(out::AbstractArray, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny)
    imfilter!(out, img, kernel, border, filter_algorithm(out, img, kernel))
end

function imfilter!(out::AbstractArray, img::AbstractArray, kernel::ProcessedKernel, border::BorderSpecAny, alg::Alg)
    imfilter!(default_resource(alg), out, img, kernel, border)
end

"""
    imfilter!(imgfilt, img, kernel, [border="replicate"], [alg])
    imfilter!(r, imgfilt, img, kernel, border, [inds])
    imfilter!(r, imgfilt, img, kernel, border::NoPad, [inds=indices(imgfilt)])

Filter an array `img` with kernel `kernel` by computing their
correlation, storing the result in `imgfilt`.

The indices of `imgfilt` determine the region over which the filtered
image is computed---you can use this fact to select just a specific
region of interest, although be aware that the input `img` might still
get padded.  Alteratively, explicitly provide the indices `inds` of
`imgfilt` that you want to calculate, and use `NoPad` boundary
conditions. In such cases, you are responsible for supplying
appropriate padding: `img` must be indexable for all of the locations
needed for calculating the output. This syntax is best-supported for
FIR filtering; in particular, that that IIR filtering can lead to
results that are inconsistent with respect to filtering the entire
array.

See also: imfilter.
"""
imfilter!

# Step 6: pad the input
# NA "padding": normalizing by the number of available values (similar to nanmean)
function imfilter!{T,S,N}(r::AbstractResource,
                          out::AbstractArray{S,N},
                          img::AbstractArray{T,N},
                          kernel::ProcessedKernel,
                          border::NA{0})
    _imfilter_na!(r, out, img, kernel, border)
end

function _imfilter_na!{T,S,N}(r::AbstractResource,
                              out::AbstractArray{S,N},
                              img::AbstractArray{T,N},
                              kernel::ProcessedKernel,
                              border::NA{0})
    nanflag = isnan.(img)
    hasnans = any(nanflag)
    if hasnans || !isseparable(kernel)
        imfilter_na_inseparable!(r, out, img, nanflag, kernel)
    else
        imfilter_na_separable!(r, out, img, kernel)
    end
    out
end

# for types that can't have NaNs, we can skip the isnan check
function _imfilter_na!{S,T<:Union{Integer,FixedColorant},N}(r::AbstractResource,
                                                            out::AbstractArray{S,N},
                                                            img::AbstractArray{T,N},
                                                            kernel::ProcessedKernel,
                                                            border::NA{0})
    if isseparable(kernel)
        imfilter_na_separable!(r, out, img, kernel)
    else
        nanflag = similar(dims->falses(dims), indices(img))
        imfilter_na_inseparable!(r, out, img, nanflag, kernel)
    end
end

# Any other kind of padding
function imfilter!{S,T,N}(r::AbstractResource,
                          out::AbstractArray{S,N},
                          img::AbstractArray{T,N},
                          kernel::ProcessedKernel,
                          border::BorderSpec)
    bord = border(kernel, img, Alg(r))  # if it's FFT, the size of img is also relevant
    A = padarray(S, img, bord)
    # By specifying NoPad(), we ensure that dispatch will never
    # accidentally "go back" to an earlier routine and apply more
    # padding
    imfilter!(r, out, A, kernel, NoPad(border))
end

# An optimized case that performs only "virtual padding"
function imfilter!{S,T,N}(r::AbstractCPU{FIR},
                          out::AbstractArray{S,N},
                          img::AbstractArray{T,N},
                          kernel::ProcessedKernel,
                          border::Pad{0})
    # The fast path: handle the points that don't need padding
    iinds = map(intersect, interior(img, kernel), indices(out))
    imfilter!(r, out, img, kernel, NoPad(border), iinds)
    # The not-so-fast path: handle the edges
    padded = view(img, padindices(img, border(kernel))...)
    pkernel = kernelconv(kernel...)
    _imfilter_iter!(r, out, padded, pkernel, EdgeIterator(indices(out), iinds))
end

### "Scheduler" methods (all with NoPad)

# These methods handle much of what Halide calls "the schedule."
# Together they handle the order-of-operations for separable and/or
# cascaded kernels, and even implement multithreadable tiling for FIR
# filtering.

# Trivial kernel (a copy operation)
function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{}, ::NoPad, inds::Indices=indices(out))
    R = CartesianRange(inds)
    copy!(out, R, A, R)
end

# A single kernel
function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any}, border::NoPad, inds::Indices=indices(out))
    kern = kernel[1]
    iscopy(kern) && return imfilter!(r, out, A, (), border, inds)
    imfilter!(r, out, A, samedims(out, kern), border, inds)
end

const fillbuf_nan = Ref(false)  # used only for testing purposes

# A filter cascade (2 or more filters)
function imfilter!(r::AbstractResource, out::AbstractArray, A::AbstractArray, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds=indices(out))
    kern = kernel[1]
    iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border, inds)
    # For multiple stages of filtering, we introduce a second buffer
    # and swap them at each stage. The first of the two is the one
    # that holds the most recent result.
    A2 = similar(A)  # for type-stability, let's hope it's really the *same* type...
    if fillbuf_nan[]
        fill!(A2, NaN)  # for testing purposes
    end
    indsstep = shrink(expand(inds, calculate_padding(kernel)), kern)
    _imfilter!(r, out, A, A2, kernel, border, indsstep)
    return out
end

### When everything is FIR, we can instead use a tiled algorithm for the cascaded case
function imfilter!{S,T,N}(r::AbstractCPU{FIR}, out::AbstractArray{S,N}, A::AbstractArray{T,N}, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, inds=indices(out))
    kern = kernel[1]
    iscopy(kern) && return imfilter!(r, out, A, tail(kernel), border, inds)
    tmp = tile_allocate(filter_type(A, kernel), kernel)
    _imfilter_tiled!(r, out, A, kernel, border, tmp, inds)
    out
end

### Scheduler support methods

## No tiling, filter cascade

# For these (internal) methods, `indsstep` refers to `inds` for this
# step of filtering, not the indices of `out` that we want to finally
# target.

# When `kernel` is (originally) a tuple that has both TriggsSdika and
# FIR filters, the overall padding gets doubled, yet we only trim off
# the minimum at each stage. Consequently, `indsstep` might be
# optimistic about the range available in `out`; therefore we use
# `intersect`.
function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{}, border::NoPad, indsstep::Indices)
    imfilter!(r, out, A1, kernel, border, map(intersect, indsstep, indices(out)))
end

function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any}, border::NoPad, indsstep::Indices)
    imfilter!(r, out, A1, kernel[1], border, map(intersect, indsstep, indices(out)))
end

# For IIR, it's important to filter over the whole passed-in range,
# and then copy! to out
function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{AnyIIR}, border::NoPad, indsstep::Indices)
    if indsstep != indices(out)
        imfilter!(r, A2, A1, kernel[1], border, indsstep)
        R = CartesianRange(map(intersect, indsstep, indices(out)))
        return copy!(out, R, A2, R)
    end
    imfilter!(r, out, A1, kernel[1], border, indsstep)
end

function _imfilter!(r, out::AbstractArray, A1, A2, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, indsstep::Indices)
    kern = kernel[1]
    iscopy(kern) && return _imfilter!(r, out, A1, A2, tail(kernel), border, indsstep)
    kernN = samedims(A2, kern)
    imfilter!(r, A2, A1, kernN, border, indsstep)  # store result in A2
    kernelt = tail(kernel)
    newinds = next_shrink(indsstep, kernelt)
    _imfilter!(r, out, A2, A1, tail(kernel), border, newinds)          # swap the buffers
end

# Single-threaded, pair of kernels (with only one temporary buffer required)
function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}, indsout)
    k1, k2 = kernel
    tile = tiles[1]
    indsk1, indstile = indices(k1), indices(tile)
    sz = map(length, indstile)
    chunksz = map(length, shrink(indstile, indsk1))
    inds = expand(indsout, k2)
    for tileinds in TileIterator(inds, chunksz)
        tileb = TileBuffer(tile, tileinds)
        imfilter!(r, tileb, A, samedims(tileb, k1), border, tileinds)
        imfilter!(r, out, tileb, samedims(out, k2), border, shrink(tileinds, k2))
    end
    out
end

# Multithreaded, pair of kernels
function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any}, border::NoPad, tiles::Vector{AA}, indsout)
    k1, k2 = kernel
    tile = tiles[1]
    indsk1, indstile = indices(k1), indices(tile)
    sz = map(length, indstile)
    chunksz = map(length, shrink(indstile, indsk1))
    tileinds_all = collect(TileIterator(expand(indsout, k2), chunksz))
    _imfilter_tiled_threads!(CPU1(r), out, A, samedims(out, k1), samedims(out, k2), border, tileinds_all, tiles)
end
# This must be in a separate function due to #15276
@noinline function _imfilter_tiled_threads!{AA<:AbstractArray}(r1, out, A, k1, k2, border, tileinds_all, tile::Vector{AA})
    Threads.@threads for i = 1:length(tileinds_all)
        id = Threads.threadid()
        tileinds = tileinds_all[i]
        tileb = TileBuffer(tile[id], tileinds)
        imfilter!(r1, tileb, A, k1, border, tileinds)
        imfilter!(r1, out, tileb, k2, border, shrink(tileinds, k2))
    end
    out
end

# Single-threaded, multiple kernels (requires two tile buffers, swapping on each iteration)
function _imfilter_tiled!{AA<:AbstractArray}(r::CPU1, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}, indsout)
    k1, kt = kernel[1], tail(kernel)
    tilepair = tiles[1]
    indsk1, indstile = indices(k1), indices(tilepair[1])
    sz = map(length, indstile)
    chunksz = map(length, shrink(indstile, indsk1))
    inds = expand(indsout, kt)
    for tileinds in TileIterator(inds, chunksz)
        tileb1 = TileBuffer(tilepair[1], tileinds)
        imfilter!(r, tileb1, A, samedims(tileb1, k1), border, tileinds)
        _imfilter_tiled_swap!(r, out, kt, border, (tileb1, tilepair[2]))
    end
end

# Multithreaded, multiple kernels
function _imfilter_tiled!{AA<:AbstractArray}(r::CPUThreads, out, A, kernel::Tuple{Any,Any,Vararg{Any}}, border::NoPad, tiles::Vector{Tuple{AA,AA}}, indsout)
    k1, kt = kernel[1], tail(kernel)
    tilepair = tiles[1]
    indsk1, indstile = indices(k1), indices(tilepair[1])
    sz = map(length, indstile)
    chunksz = map(length, shrink(indstile, indsk1))
    tileinds_all = collect(TileIterator(expand(indsout, kt), chunksz))
    _imfilter_tiled_threads!(CPU1(r), out, A, samedims(out, k1), kt, border, tileinds_all, tiles)
end
# This must be in a separate function due to #15276
@noinline function _imfilter_tiled_threads!{AA<:AbstractArray}(r1, out, A, k1, kt, border, tileinds_all, tiles::Vector{Tuple{AA,AA}})
    Threads.@threads for i = 1:length(tileinds_all)
        tileinds = tileinds_all[i]
        id = Threads.threadid()
        tile1, tile2 = tiles[id]
        tileb1 = TileBuffer(tile1, tileinds)
        imfilter!(r1, tileb1, A, k1, border, tileinds)
        _imfilter_tiled_swap!(r1, out, kt, border, (tileb1, tile2))
    end
    out
end

# The first of the pair in `tmp` has the current data. We also make
# the second a plain array so there's no doubt about who's holding the
# proper indices.
function _imfilter_tiled_swap!(r, out, kernel::Tuple{Any,Any,Vararg{Any}}, border, tmp::Tuple{TileBuffer,Array})
    tileb1, tile2 = tmp
    k1, kt = kernel[1], tail(kernel)
    parentinds = indices(tileb1)
    tileinds = shrink(parentinds, k1)
    tileb2 = TileBuffer(tile2, tileinds)
    imfilter!(r, tileb2, tileb1, samedims(tileb2, k1), border, tileinds)
    _imfilter_tiled_swap!(r, out, kt, border, (tileb2, parent(tileb1)))
end

# on the last call we write to `out` instead of one of the buffers
function _imfilter_tiled_swap!(r, out, kernel::Tuple{Any}, border, tmp::Tuple{TileBuffer,Array})
    tileb1 = tmp[1]
    k1 = kernel[1]
    parentinds = indices(tileb1)
    tileinds = shrink(parentinds, k1)
    imfilter!(r, out, tileb1, samedims(out, k1), border, tileinds)
end

### FIR filtering

"""
    imfilter!(::AbstractResource, imgfilt, img, kernel, NoPad(), [inds=indices(imgfilt)])

Filter an array `img` with kernel `kernel` by computing their
correlation, storing the result in `imgfilt`, defaulting to a finite-impulse
response (FIR) algorithm. Any necessary padding must have already been
supplied to `img`. If you want padding applied, instead call

    imfilter!([r::AbstractResource,] imgfilt, img, kernel, border)

with a specific `border`, or use

    imfilter!(imgfilt, img, kernel, [Algorithm.FIR()])

for default padding.

If `inds` is supplied, only the elements of `imgfilt` with indices in
the domain of `inds` will be calculated. This can be particularly
useful for "cascaded FIR filters" where you pad over a larger area and
then calculate the result over just the necessary/well-defined region
at each successive stage.

See also: imfilter.
"""
function imfilter!{S,T,N}(r::AbstractResource,
                          out::AbstractArray{S,N},
                          A::AbstractArray{T,N},
                          kern::NDimKernel{N},
                          border::NoPad,
                          inds::Indices{N}=indices(out))
    (isempty(A) || isempty(kern)) && return out
    indso, indsA, indsk = indices(out), indices(A), indices(kern)
    if iscopy(kern)
        R = CartesianRange(inds)
        return copy!(out, R, A, R)
    end
    for i = 1:N
        # Check that inds is inbounds for out
        indsi, indsoi, indsAi, indski = inds[i], indso[i], indsA[i], indsk[i]
        if first(indsi) < first(indsoi) || last(indsi) > last(indsoi)
            throw(DimensionMismatch("output indices $indso disagrees with requested indices $inds"))
        end
        # Check that input A is big enough not to throw a BoundsError
        if      first(indsAi) > first(indsi) + first(indski) ||
                last(indsA[i])  < last(indsi)  + last(indski)
            throw(DimensionMismatch("requested indices $inds and kernel indices $indsk do not agree with indices of padded input, $indsA"))
        end
    end
    _imfilter_inbounds!(r, out, A, kern, border, inds)
end

function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, kern::ReshapedIIR, border::NoPad, inds)
    indspre, ind, indspost = iterdims(inds, kern)
    _imfilter_dim!(r, out, A, kern.data, CartesianRange(indspre), ind, CartesianRange(indspost), border[])
end

function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, kern, border::NoPad, inds)
    indsk = indices(kern)
    R, Rk = CartesianRange(inds), CartesianRange(indsk)
    p = A[first(R)+first(Rk)] * first(kern)
    TT = typeof(p+p)
    for I in R
        tmp = zero(TT)
        @unsafe for J in Rk
            tmp += A[I+J]*kern[J]
        end
        @unsafe out[I] = tmp
    end
    out
end

function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, kern::ReshapedOneD, border::NoPad, inds)
    Rpre, ind, Rpost = iterdims(inds, kern)
    k = kern.data
    R, Rk = CartesianRange(inds), CartesianRange(indices(kern))
    p = A[first(R)+first(Rk)] * first(k)
    TT = typeof(p+p)
    _imfilter_inbounds(r, TT, out, A, k, Rpre, ind, Rpost)
end

function _imfilter_inbounds{TT}(r::AbstractResource, ::Type{TT}, out, A::AbstractArray, k::AbstractVector, Rpre::CartesianRange, ind, Rpost::CartesianRange)
    indsk = indices(k, 1)
    for Ipost in Rpost
        for Ipre in Rpre
            for i in ind
                tmp = zero(TT)
                @unsafe for j in indsk
                    tmp += A[Ipre,i+j,Ipost]*k[j]
                end
                @unsafe out[Ipre,i,Ipost] = tmp
            end
        end
    end
    out
end

function _imfilter_iter!(r::AbstractResource, out, padded, kernel::AbstractArray, iter)
    p = first(padded) * first(kernel)
    TT = typeof(p+p)
    Rk = CartesianRange(indices(kernel))
    for I in iter
        tmp = zero(TT)
        for J in Rk
            tmp += padded[I+J]*kernel[J]
        end
        out[I] = tmp
    end
    out
end

### FFT filtering

"""
    imfilter!(::AbstractResource{FFT}, imgfilt, img, kernel, NoPad())

Filter an array `img` with kernel `kernel` by computing their
correlation, storing the result in `imgfilt`, using a fast Fourier
transform (FFT) algorithm. Any necessary padding must have already
been applied to `img`. If you want padding applied, instead call

    imfilter!(::AbstractResource{FFT}, imgfilt, img, kernel, border)

with a specific `border`, or use

    imfilter!(imgfilt, img, kernel, Algorithm.FFT())

for default padding.

See also: imfilter.
"""
function imfilter!{S,T,K,N}(r::AbstractCPU{FFT},
                            out::AbstractArray{S,N},
                            img::AbstractArray{T,N},
                            kernel::AbstractArray{K,N},
                            border::NoPad)
    imfilter!(r, out, img, (kernel,), border)
end

function imfilter!{S,T,N}(r::AbstractCPU{FFT},
                          out::AbstractArray{S,N},
                          A::AbstractArray{T,N},
                          kernel::Tuple{AbstractArray},
                          border::NoPad)
    _imfilter_fft!(r, out, A, kernel, border)  # ambiguity resolution
end
function imfilter!{S,T,N}(r::AbstractCPU{FFT},
                          out::AbstractArray{S,N},
                          A::AbstractArray{T,N},
                          kernel::Tuple{AbstractArray,Vararg{AbstractArray}},
                          border::NoPad)
    _imfilter_fft!(r, out, A, kernel, border)
end

function _imfilter_fft!{S,T,N}(r::AbstractCPU{FFT},
                          out::AbstractArray{S,N},
                          A::AbstractArray{T,N},
                          kernel::Tuple{AbstractArray,Vararg{AbstractArray}},
                          border::NoPad)
    kern = samedims(A, kernelconv(kernel...))
    krn = FFTView(zeros(eltype(kern), map(length, indices(A))))
    for I in CartesianRange(indices(kern))
        krn[I] = kern[I]
    end
    Af = filtfft(A, krn)
    if map(first, indices(out)) == map(first, indices(Af))
        R = CartesianRange(indices(out))
        copy!(out, R, Af, R)
    else
        # Exploit the periodic boundary conditions of FFTView
        dest = FFTView(out)
        src = OffsetArray(view(FFTView(Af), indices(dest)...), indices(dest))
        copy!(dest, src)
    end
    out
end

filtfft(A, krn) = irfft(rfft(A).*conj(rfft(krn)), length(indices(A,1)))
function filtfft{C<:Colorant}(A::AbstractArray{C}, krn)
    Av, dims = channelview_dims(A)
    kernrs = kreshape(C, krn)
    Avf = irfft(rfft(Av, dims).*conj(rfft(kernrs, dims)), length(indices(Av, dims[1])), dims)
    colorview(base_colorant_type(C){eltype(Avf)}, Avf)
end
channelview_dims{C<:Colorant,N}(A::AbstractArray{C,N}) = channelview(A), ntuple(d->d+1, Val{N})
if ImageCore.squeeze1
    channelview_dims{C<:ImageCore.Color1,N}(A::AbstractArray{C,N}) = channelview(A), ntuple(identity, Val{N})
end

function kreshape{C<:Colorant}(::Type{C}, krn::FFTView)
    kern = parent(krn)
    kernrs = FFTView(reshape(kern, 1, size(kern)...))
end
if ImageCore.squeeze1
    kreshape{C<:ImageCore.Color1}(::Type{C}, krn::FFTView) = krn
end

### Triggs-Sdika (modified Young-van Vliet) recursive filtering
# B. Triggs and M. Sdika, "Boundary conditions for Young-van Vliet
# recursive filtering". IEEE Trans. on Sig. Proc. 54: 2365-2367
# (2006).

# Note this is safe for inplace use, i.e., out === img

function imfilter!{S,T,N}(r::AbstractResource{IIR},
                          out::AbstractArray{S,N},
                          img::AbstractArray{T,N},
                          kernel::Tuple{TriggsSdika, Vararg{TriggsSdika}},
                          border::BorderSpec)
    isa(border, Pad) && border.symbol != :replicate && throw(ArgumentError("only \"replicate\" is supported"))
    length(kernel) <= N || throw(DimensionMismatch("cannot have more kernels than dimensions"))
    inds = indices(img)
    _imfilter_inplace_tuple!(r, out, img, kernel, CartesianRange(()), inds, CartesianRange(tail(inds)), border)
end

"""
    imfilter!(r::AbstractResource, imgfilt, img, kernel::Tuple{TriggsSdika...}, border)
    imfilter!(r::AbstractResource, imgfilt, img, kernel::TriggsSdika, dim::Integer, border)

Filter an array `img` with a Triggs-Sdika infinite impulse response
(IIR) `kernel`, storing the result in `imgfilt`. Unlike the `FIR` and
`FFT` algorithms, this version is safe for inplace operations, i.e.,
`imgfilt` can be the same array as `img`.

Either specify one kernel per dimension (as a tuple), or a particular
dimension `dim` along which to filter. If you exhaust `kernel`s before
you run out of array dimensions, the remaining dimension(s) will not
be filtered.

With Triggs-Sdika filtering, the only border options are `NA()`,
`"replicate"`, or `Fill(value)`.

See also: imfilter, TriggsSdika, IIRGaussian.
"""
function imfilter!(r::AbstractResource, out, img, kernel::TriggsSdika, dim::Integer, border::BorderSpec)
    inds = indices(img)
    k, l = length(kernel.a), length(kernel.b)
    # This next part is not type-stable, which is why _imfilter_dim! has a @noinline
    Rbegin = CartesianRange(inds[1:dim-1])
    Rend   = CartesianRange(inds[dim+1:end])
    _imfilter_dim!(r, out, img, kernel, Rbegin, inds[dim], Rend, border)
end
function imfilter!(r::AbstractResource, out, img, kernel::TriggsSdika, dim::Integer, border::AbstractString)
    imfilter!(r, out, img, kernel, dim, Pad(Symbol(border)))
end

function imfilter!(r::AbstractResource, out, A::AbstractVector, kern::TriggsSdika, border::NoPad, inds::Indices=indices(out))
    indspre, ind, indspost = iterdims(inds, kern)
    _imfilter_dim!(r, out, A, kern, CartesianRange(indspre), ind, CartesianRange(indspost), border[])
end

# Lispy and type-stable inplace (currently just Triggs-Sdika) filtering over each dimension
function _imfilter_inplace_tuple!(r, out, img, kernel, Rbegin, inds, Rend, border)
    ind = first(inds)
    _imfilter_dim!(r, out, img, first(kernel), Rbegin, ind, Rend, border)
    _imfilter_inplace_tuple!(r,
                             out,
                             out,
                             tail(kernel),
                             CartesianRange(CartesianIndex((Rbegin.start.I..., first(ind))),
                                            CartesianIndex((Rbegin.stop.I...,  last(ind)))),
                             tail(inds),
                             _tail(Rend),
                             border)
end
# When the final kernel has been used, return the output
_imfilter_inplace_tuple!(r, out, img, ::Tuple{}, Rbegin, inds, Rend, border) = out

# This is the "workhorse" function that performs Triggs-Sdika IIR
# filtering along a particular dimension. The "pre" dimensions are
# encoded in Rbegin, the "post" dimensions in Rend, and the dimension
# we're filtering is sandwiched between these. This design is
# type-stable and cache-friendly for any dimension---we update values
# in memory-order rather than along the chosen dimension. Nor does it
# require that the arrays have efficient linear indexing. For more
# information, see http://julialang.org/blog/2016/02/iteration.
@noinline function _imfilter_dim!{T,k,l}(r::AbstractResource,
                                         out, img, kernel::TriggsSdika{T,k,l},
                                         Rbegin::CartesianRange, ind::AbstractUnitRange,
                                         Rend::CartesianRange, border::BorderSpec)
    if iscopy(kernel)
        if !(out === img)
            copy!(out, img)
        end
        return out
    end
    if length(ind) <= max(k, l)
        throw(DimensionMismatch("size of img along dimension $dim $(length(inds[dim])) is too small for filtering with IIR kernel of length $(max(k,l))"))
    end
    for Iend in Rend
        # Initialize the left border
        indleft = ind[1:k]
        for Ibegin in Rbegin
            leftborder!(out, img, kernel, Ibegin, indleft, Iend, border)
        end
        # Propagate forwards. We omit the final point in case border
        # is "replicate", so that the original value is still
        # available. rightborder! will handle that point.
        for i = ind[k]+1:ind[end-1]
            @unsafe for Ibegin in Rbegin
                tmp = one(T)*img[Ibegin, i, Iend]
                for j = 1:k
                    tmp += kernel.a[j]*out[Ibegin, i-j, Iend]
                end
                out[Ibegin, i, Iend] = tmp
            end
        end
        # Initialize the right border
        indright = ind[end-l+1:end]
        for Ibegin in Rbegin
            rightborder!(out, img, kernel, Ibegin, indright, Iend, border)
        end
        # Propagate backwards
        for i = ind[end-l]:-1:ind[1]
            @unsafe for Ibegin in Rbegin
                tmp = one(T)*out[Ibegin, i, Iend]
                for j = 1:l
                    tmp += kernel.b[j]*out[Ibegin, i+j, Iend]
                end
                out[Ibegin, i, Iend] = tmp
            end
        end
        # Final scaling
        for i in ind
            @unsafe for Ibegin in Rbegin
                out[Ibegin, i, Iend] *= kernel.scale
            end
        end
    end
    out
end

# Implements the initialization in the first paragraph of Triggs & Sdika, section II
function leftborder!(out, img, kernel, Ibegin, indleft, Iend, border::Fill)
    _leftborder!(out, img, kernel, Ibegin, indleft, Iend, convert(eltype(img), border.value))
end
function leftborder!(out, img, kernel, Ibegin, indleft, Iend, border::Pad)
    _leftborder!(out, img, kernel, Ibegin, indleft, Iend, img[Ibegin, indleft[1], Iend])
end
function _leftborder!{T,k,l}(out, img, kernel::TriggsSdika{T,k,l}, Ibegin, indleft, Iend, iminus)
    uminus = iminus/(1-kernel.asum)
    n = 0
    for i in indleft
        n += 1
        tmp = one(T)*img[Ibegin, i, Iend]
        for j = 1:n-1
            tmp += kernel.a[j]*out[Ibegin, i-j, Iend]
        end
        for j = n:k
            tmp += kernel.a[j]*uminus
        end
        out[Ibegin, i, Iend] = tmp
    end
    out
end

# Implements Triggs & Sdika, Eqs 14-15
function rightborder!(out, img, kernel, Ibegin, indright, Iend, border::Fill)
    _rightborder!(out, img, kernel, Ibegin, indright, Iend, convert(eltype(img), border.value))
end
function rightborder!(out, img, kernel, Ibegin, indright, Iend, border::Pad)
    _rightborder!(out, img, kernel, Ibegin, indright, Iend, img[Ibegin, indright[end], Iend])
end
function _rightborder!{T,k,l}(out, img, kernel::TriggsSdika{T,k,l}, Ibegin, indright, Iend, iplus)
    # The final value from forward-filtering was not calculated, so do that here
    i = last(indright)
    tmp = one(T)*img[Ibegin, i, Iend]
    for j = 1:k
        tmp += kernel.a[j]*out[Ibegin, i-j, Iend]
    end
    out[Ibegin, i, Iend] = tmp
    # Initialize the v values at and beyond the right edge
    uplus = iplus/(1-kernel.asum)
    vplus = uplus/(1-kernel.bsum)
    vright = kernel.M * rightΔu(out, uplus, Ibegin, last(indright), Iend, kernel) + vplus
    out[Ibegin, last(indright), Iend] = vright[1]
    # Propagate inward
    n = 1
    for i in last(indright)-1:-1:first(indright)
        n += 1
        tmp = one(T)*out[Ibegin, i, Iend]
        for j = 1:n-1
            tmp += kernel.b[j]*out[Ibegin, i+j, Iend]
        end
        for j = n:l
            tmp += kernel.b[j]*vright[j-n+2]
        end
        out[Ibegin, i, Iend] = tmp
    end
    out
end

# Part of Triggs & Sdika, Eq. 14
function rightΔu{T,l}(img, uplus, Ibegin, i, Iend, kernel::TriggsSdika{T,3,l})
    @inbounds ret = SVector(img[Ibegin, i,   Iend]-uplus,
                            img[Ibegin, i-1, Iend]-uplus,
                            img[Ibegin, i-2, Iend]-uplus)
    ret
end

### NA boundary conditions

function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple{Vararg{TriggsSdika}})
    fc, fn = Fill(zero(T)), Fill(zero(eltype(T)))  # color, numeric
    copy!(out, img)
    out[nanflag] = zero(T)
    validpixels = copy!(similar(Array{eltype(T)}, indices(img)), mappedarray(x->!x, nanflag))
    # TriggsSdika is safe for inplace operations
    imfilter!(r, out, out, kernel, fc)
    imfilter!(r, validpixels, validpixels, kernel, fn)
    for I in eachindex(out)
        out[I] /= validpixels[I]
    end
    out[nanflag] = nan(T)
    out
end

function imfilter_na_inseparable!{T}(r, out::AbstractArray{T}, img, nanflag, kernel::Tuple)
    fc, fn = Fill(zero(T)), Fill(zero(eltype(T)))  # color, numeric
    imgtmp = copy!(similar(out, indices(img)), img)
    imgtmp[nanflag] = zero(T)
    validpixels = copy!(similar(Array{eltype(T)}, indices(img)), mappedarray(x->!x, nanflag))
    imfilter!(r, out, imgtmp, kernel, fc)
    vp = imfilter(r, validpixels, kernel, fn)
    for I in eachindex(out)
        out[I] /= vp[I]
    end
    out[nanflag] = nan(T)
    out
end

function imfilter_na_separable!{T}(r, out::AbstractArray{T}, img, kernel::Tuple)
    fc, fn = Fill(zero(T)), Fill(zero(eltype(T)))  # color, numeric
    imfilter!(r, out, img, kernel, fc)
    normalize_separable!(r, out, kernel, fn)
end

### Utilities

filter_type{S}(img::AbstractArray{S}, kernel) = filter_type(S, kernel)

filter_type{S,T}(::Type{S}, kernel::ArrayLike{T}) = typeof(zero(S)*zero(T) + zero(S)*zero(T))
filter_type{S<:Union{UFixed,FixedColorant}}(::Type{S}, ::Laplacian) = float32(S)
filter_type{S<:Colorant}(::Type{S}, kernel::Laplacian) = S
filter_type{S<:AbstractFloat}(::Type{S}, ::Laplacian) = S
filter_type{S<:Signed}(::Type{S}, ::Laplacian) = S
filter_type{S<:Unsigned}(::Type{S}, ::Laplacian) = signed_type(S)
filter_type(::Type{Bool}, ::Laplacian) = Int8

signed_type(::Type{UInt8})  = Int16
signed_type(::Type{UInt16}) = Int32
signed_type(::Type{UInt32}) = Int64

@inline function filter_type{S}(::Type{S}, kernel::Tuple{Any,Vararg{Any}})
    T = filter_type(S, kernel[1])
    filter_type(T, S, tail(kernel))
end
@inline function filter_type{T,S}(::Type{T}, ::Type{S}, kernel::Tuple{Any,Vararg{Any}})
    Tnew = promote_type(T, filter_type(S, kernel[1]))
    filter_type(Tnew, S, tail(kernel))
end
filter_type{T,S}(::Type{T}, ::Type{S}, kernel::Tuple{}) = T

factorkernel(kernel::AbstractArray) = (kernelshift(indices(kernel), kernel),)
factorkernel(L::Laplacian) = (L,)

function factorkernel{T}(kernel::AbstractMatrix{T})
    inds = indices(kernel)
    m, n = map(length, inds)
    kern = Array{T}(m, n)
    copy!(kern, 1:m, 1:n, kernel, inds[1], inds[2])
    factorstridedkernel(inds, kern)
end

function factorstridedkernel(inds, kernel::StridedMatrix)
    SVD = svdfact(kernel)
    U, S, Vt = SVD[:U], SVD[:S], SVD[:Vt]
    separable = true
    EPS = sqrt(eps(eltype(S)))
    for i = 2:length(S)
        separable &= (abs(S[i]) < EPS)
    end
    if !separable
        ks = kernelshift(inds, kernel)
        return (dummykernel(indices(ks)), ks)
    end
    s = S[1]
    u, v = U[:,1:1], Vt[1:1,:]
    ss = sqrt(s)
    (kernelshift((inds[1], dummyind(inds[1])), ss*u),
     kernelshift((dummyind(inds[2]), inds[2]), ss*v))
end

prod_kernel(kern::AbstractArray) = kern
prod_kernel(kern::AbstractArray, kern1, kerns...) = broadcast(*, kern, kern1, kerns...) #prod_kernel(kern.*kern1, kerns...)
prod_kernel{N}(::Type{Val{N}}, args...) = prod_kernel(Val{N}, prod_kernel(args...))
prod_kernel{_,N}(::Type{Val{N}}, kernel::AbstractArray{_,N}) = kernel
function prod_kernel{N}(::Type{Val{N}}, kernel::AbstractArray)
    inds = indices(kernel)
    newinds = fill_to_length(inds, oftype(inds[1], 0:0), Val{N})
    reshape(kernel, newinds)
end

kernelshift(::Tuple{}, A::StridedArray) = A
kernelshift(::Tuple{}, A) = A
kernelshift{N}(inds::NTuple{N,Base.OneTo}, A::StridedArray) = _kernelshift(inds, A)
kernelshift{N}(inds::NTuple{N,Base.OneTo}, A) = _kernelshift(inds, A)
function _kernelshift(inds, A)
    Base.depwarn("assuming that the origin is at the center of the kernel; to avoid this warning, call `centered(kernel)` or use an OffsetArray", :_kernelshift)  # this may be necessary long-term?
    centered(A)
end
kernelshift(inds::Indices, A::StridedArray) = OffsetArray(A, inds...)
function kernelshift(inds::Indices, A)
    @assert indices(A) == inds
    A
end

# Note this is not type-stable. Fortunately, all the outputs are
# allocated by the time this gets called.
function filter_algorithm(out, img, kernel::Union{ArrayType,Tuple{Vararg{ArrayType}}})
    L = maxlen(kernel)
    if L > 30
        return FFT()
    end
    FIR()
end
filter_algorithm(out, img, kernel::Tuple{AnyIIR,Vararg{AnyIIR}}) = IIR()
filter_algorithm(out, img, kernel) = Mixed()

maxlen(A::AbstractArray) = _length(A)
@inline maxlen(kernel::Tuple) = _maxlen(0, kernel...)
_maxlen(len, kernel1, kernel...) = _maxlen(max(len, _length(kernel1)), kernel...)
_maxlen(len) = len

_length(A::AbstractArray) = length(linearindices(A))
_length(A::ReshapedVector) = _length(A.data)

isseparable(kernels::Tuple{Vararg{AnyIIR}}) = true
isseparable(kernels::Tuple) = all(x->nextendeddims(x)==1, kernels)

normalize_separable!(r::AbstractResource, A, ::Tuple{}, border) = error("this shouldn't happen")
function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,TriggsSdika}, border)
    inds = indices(A)
    function imfilter_inplace!(r, a, kern, border)
        imfilter!(r, a, a, (kern,), border)
    end
    filtdims = ntuple(d->imfilter_inplace!(r, similar(dims->ones(dims), inds[d]), kernels[d], border), Val{N})
    normalize_dims!(A, filtdims)
end
function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N,ReshapedIIR}, border)
    normalize_separable!(r, A, map(_vec, kernels), border)
end

function normalize_separable!{N}(r::AbstractResource, A, kernels::NTuple{N}, border)
    inds = indices(A)
    filtdims = ntuple(d->imfilter(r, similar(dims->ones(dims), inds[d]), _vec(kernels[d]), border), Val{N})
    normalize_dims!(A, filtdims)
end

function normalize_dims!{T,N}(A::AbstractArray{T,N}, factors::NTuple{N})
    for I in CartesianRange(indices(A))
        tmp = A[I]/factors[1][I[1]]
        for d = 2:N
            tmp /= factors[d][I[d]]
        end
        A[I] = tmp
    end
    A
end

iscopy(kernel::AbstractArray) = all(x->x==0:0, indices(kernel)) && first(kernel) == 1
iscopy(kernel::Laplacian) = false
iscopy(kernel::TriggsSdika) = all(x->x==0, kernel.a) && all(x->x==0, kernel.b) && kernel.scale == 1
iscopy(kernel::ReshapedOneD) = iscopy(kernel.data)

kernelconv(kernel) = kernel
function kernelconv(k1, k2, kernels...)
    out = similar(Array{filter_type(eltype(k1), k2)}, calculate_padding((k1, k2)))
    fill!(out, zero(eltype(out)))
    k1N, k2N = samedims(out, k1), samedims(out, k2)
    R1, R2 = CartesianRange(indices(k1N)), CartesianRange(indices(k2N))
    for I1 in R1
        for I2 in R2
            out[I1+I2] += k1N[I1]*k2N[I2]
        end
    end
    kernelconv(out, kernels...)
end

default_resource(alg::FIR) = Threads.nthreads() > 1 ? CPUThreads(alg) : CPU1(alg)
default_resource(alg) = CPU1(alg)

## Tiling utilities

function tile_allocate{T}(::Type{T}, kernel::Tuple{Any,Any})
    sz = map(length, calculate_padding(kernel))
    tsz = padded_tilesize(T, sz)
    [Array{T}(tsz) for i = 1:Threads.nthreads()]
end

function tile_allocate{T}(::Type{T}, kernel::Tuple{Any,Any,Vararg{Any}})
    # Allocate a pair of tiles and swap buffers at each stage
    sz = map(length, calculate_padding(kernel))
    tsz = padded_tilesize(T, sz)
    [(Array{T}(tsz), Array{T}(tsz)) for i = 1:Threads.nthreads()]
end

ext: jl


text:
#Include KG_Module
include("../src/KG_Module.jl")

#Get field at certain positions
a , b = get_field_point("./data/binary_run_1/binary_bh_rotating_radius_1.5_alpha_10.0_Rorbit_5.0.h5", 0.0 , 0.0)

#Create plot


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Test

examples = ["ME_simulate", "CS_simulate", "multiple_instances", "parameterize", "modelica_conference_2021"]

@testset "FMI.jl Examples" begin
    for example in examples
        @testset "$(example).jl" begin
            path = joinpath(dirname(@__FILE__), "..", "example", example * ".jl")
            @test include(path) == nothing
        end
    end
end


text:
@testset "transforms" begin

    # Tests for each flow that can't be expressed in terms of the standard interface.
    include("invertibles/affine.jl")
    include("invertibles/planar.jl")
    include("invertibles/radial.jl")

    # Tests which are generic to all of the flows. Note that these tests essentially define
    # this packages' conventions, and ensure that each flow adhere's to them.
    for Transform in __transforms, D in [1, 5], M in [1, 9]

        # Test construction using *_init!, *_init and Foo(::typeof(*_init), ...).
        # Basic tests to make sure that the functions run.
        let rng = MersenneTwister(123456)
            for (init, inplace_init) in [(naive_init, naive_init!),
                                         (identity_init, identity_init!)]
                θ = fill!(Vector{Float64}(nparams(Transform, D)), NaN)
                @test all(.!isnan.(inplace_init(θ, Transform, rng, D)))
                @test length(θ) == length(init(Transform, rng, D))
                @test nparams(Transform(init, rng, D)) == length(θ)
            end
        end

        let rng = MersenneTwister(123456), transform = Transform(naive_init, rng, D)

            # Test naive construction dimensionality.
            @test dim(transform) == D

            # Check that AssertionErrors are thrown for incorrect sized args, and that
            # outputs are not NaN or Inf.
            for foo in [apply, invert, logdetJ]
                if method_exists(foo, Tuple{Transform, AbstractVector{<:Real}})
                    @test_throws AssertionError foo(transform, zeros(D + 1))
                    @test_throws AssertionError foo(transform, zeros(D - 1))
                    @test all(.!isnan.(foo(transform, randn(rng, D))))
                    @test all(.!isinf.(foo(transform, randn(rng, D))))
                end
                if method_exists(foo, Tuple{Transform, AbstractMatrix{<:Real}})
                    @test_throws AssertionError foo(transform, zeros(D + 1, M))
                    @test_throws AssertionError foo(transform, zeros(D - 1, M))
                    @test all(.!isnan.(foo(transform, randn(rng, D, M))))
                    @test all(.!isinf.(foo(transform, randn(rng, D, M))))
                end
            end

            # Check that size conventions are obeyed.
            for foo in [apply, invert]
                if method_exists(foo, Tuple{Transform, AbstractVector{<:Real}})
                    @test size(foo(transform, zeros(D))) == (D,)
                end
                if method_exists(foo, Tuple{Transform, AbstractMatrix{<:Real}})
                    @test size(foo(transform, zeros(D, M))) == (D, M)
                end
            end
            @test size(logdetJ(transform, zeros(D))) == ()
            @test size(logdetJ(transform, zeros(D, M))) == (M,)

            # Check that broadcasting yields the same results as vectorised version.
            for foo in [apply, invert]
                if method_exists(foo, Tuple{Transform, AbstractVector{<:Real}}) &&
                    method_exists(foo, Tuple{Transform, AbstractMatrix{<:Real}})
                    X = randn(rng, D, M)
                    Y_iterated = hcat([foo(transform, X[:, m]) for m in 1:M]...)
                    @test Y_iterated ≈ foo(transform, X)
                end
            end
            if method_exists(logdetJ, Tuple{Transform, AbstractVector{<:Real}}) &&
                method_exists(logdetJ, Tuple{Transform, AbstractMatrix{<:Real}})
                X = randn(rng, D, M)
                Y_iterated = vcat([logdetJ(transform, X[:, m]) for m in 1:M]...)
                @test Y_iterated ≈ logdetJ(transform, X)
            end

            # Check that `invert` is indeed the inverse of `apply`, if both are defined.
            if method_exists(apply, Tuple{Transform, AbstractVector{<:Real}}) &&
                method_exists(invert, Tuple{Transform, AbstractVector{<:Real}})
                let x1 = randn(rng, D)
                    @test apply(transform, invert(transform, x1)) ≈ x1
                    @test invert(transform, apply(transform, x1)) ≈ x1
                end
            end
            if method_exists(apply, Tuple{Transform, AbstractMatrix{<:Real}}) &&
                method_exists(invert, Tuple{Transform, AbstractMatrix{<:Real}})
                let x2 = randn(rng, D, M)
                    @test apply(transform, invert(transform, x2)) ≈ x2
                    @test invert(transform, apply(transform, x2)) ≈ x2
                end
            end

            # Ensure that the naive_init methods don't provide insane initialisations.
            let transform = Transform(naive_init, rng, D)
                for foo in [apply, invert, logdetJ]
                    if method_exists(foo, Tuple{Transform, AbstractVector{<:Real}})
                        @test all(.!isnan.(foo(transform, randn(rng, D))))
                        @test all(.!isinf.(foo(transform, randn(rng, D))))
                    end
                    if method_exists(foo, Tuple{Transform, AbstractMatrix{<:Real}})
                        @test all(.!isnan.(foo(transform, randn(rng, D, M))))
                        @test all(.!isinf.(foo(transform, randn(rng, D, M))))
                    end
                end
            end

            # Check that identity initialisations indeed yield the identity function.
            if method_exists(Transform, Tuple{typeof(identity_init), AbstractRNG, Int})
                identity_transform = Transform(identity_init, rng, D)
                let x1 = randn(rng, D), x2 = randn(rng, D, M)
                    @test apply(identity_transform, x1) ≈ x1
                    @test apply(identity_transform, x2) ≈ x2
                    @test all(abs.(logdetJ(identity_transform, x1)) .< 1e-12)
                    @test all(abs.(logdetJ(identity_transform, x2)) .< 1e-12)
                end
            end

            @test nparams(transform) == nparams(Transform, D)
        end # let
    end # for Transform...
end


text:
#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# # [PSSFSS](https://github.com/simonp0420/PSSFSS) Examples
#


text:
module NanoDates

export NanoDate,
        nnow, ntoday, date_time,
        nanodate2rata, rata2nanodate,
        nanodate2unixnanos, nanodate2unixmicros, nanodate2unixmillis, nanodate2unixseconds,
        unixnanos2nanodate, unixmicros2nanodate, unixmillis2nanodate, unixseconds2nanodate,
        format,
        CapitalT, SmallCapitalT,
        SingleSpace, Underscore,
        SmallWhiteCircle, SmallWhiteStar

using Dates
using Dates: AbstractDateTime, CompoundPeriod,
             format, value, days,
             toms, tons, UTD, UTM

export AbstractDateTime, CompoundPeriod,
    value, days, format

include("constants.jl")
include("nanodate.jl")
include("accessors.jl")
include("conversions.jl")
include("compare.jl")
include("interop.jl")
include("arith.jl")
include("fastintmath.jl")
include("compound.jl")
include("strings.jl")
include("timestamp.jl")
include("io.jl")

@inline nanos_elapsed() = time_ns() - NanosAtStart

function __init__()
    global NanosAtStart = time_ns()
end

end  # NanoDates


text:
module SmartMeterAnalytics

using Base.Dates

export tou

include("tou.jl")

end # module

max_stars_repo_name: GordStephen/SmartMeterAnalytics

max_issues_repo_name: GordStephen/SmartMeterAnalytics

max_forks_repo_name: GordStephen/SmartMeterAnalytics


text:
module BenchDictHistogram

using BangBang.Extras: modify!!
using BenchmarkTools
using ConcurrentCollections

function generate(; datasize = 2^19, nkeys = datasize)
    ints = rand(1:nkeys, datasize)
    lastkey = string(nkeys)
    prefix = suffix = ""
    # prefix = "9" ^ 30  # adding works for hashing and comparison
    # suffix = "0" ^ 30  # adding works for hashing (but not for comparison)
    data = prefix .* string.(ints; pad = length(lastkey)) .* suffix
    return data
end

function hist_seq!(dict::AbstractDict{<:Any,Value}, data) where {Value}
    dict′ = nothing
    for k in data
        dict′, _ = modify!!(dict, k) do v
            Base.@_inline_meta
            Some(something(v, zero(Value)) + oneunit(Value))
        end
    end
    @assert dict′ === dict
    return dict
end

function hist_seq!(dict::ConcurrentDict, data)
    for k in data
        modify!(dict, k) do ref
            Base.@_inline_meta
            Some(ref === nothing ? 1 : ref[] + 1)
        end
    end
    return dict
end

function hist_parallel!(dict::ConcurrentDict, data; ntasks = Threads.nthreads())
    # for k in data
    #     dict[k] = 0
    # end
    @sync for chunk in Iterators.partition(data, cld(length(data), ntasks))
        Threads.@spawn hist_seq!(dict, chunk)
    end
    return dict
end

function hist_dac_impl(data, chunk_starts, basesize)
    if length(chunk_starts) == 0
        return Dict{String,Int}()
    elseif length(chunk_starts) == 1
        i = @inbounds chunk_starts[begin]
        chunk = @inbounds data[i:min(i + basesize - 1, end)]
        return hist_seq!(Dict{String,Int}(), chunk)
    else
        h = length(chunk_starts) ÷ 2
        left_chunk = @view chunk_starts[begin:begin+h-1]
        right_chunk = @view chunk_starts[begin+h:end]
        task = Threads.@spawn hist_dac_impl(data, right_chunk, basesize)
        left = hist_dac_impl(data, left_chunk, basesize)
        right = fetch(task)::typeof(left)
        return mergewith!(+, left, right)
    end
end

function hist_parallel_dac(data; ntasks = Threads.nthreads())
    basesize = cld(length(data), ntasks)
    chunk_starts = firstindex(data):basesize:lastindex(data)
    return hist_dac_impl(data, chunk_starts, basesize)
end

function default_ntasks_list()
    ntasks_list = [Threads.nthreads()]
    if Threads.nthreads() > 2
        pushfirst!(ntasks_list, 2)
    end
    return ntasks_list
end

const CACHE = Ref{Any}()

function setup(;
    ntasks_list = default_ntasks_list(),
    datasize = 2^19,
    nkeys_list = [datasize],
)
    CACHE[] = data = Dict(0 => generate(; nkeys = 1))  # dummy data for allocation
    empty!(data)
    T = typeof(data)

    suite = BenchmarkGroup()
    for nkeys in nkeys_list
        data[nkeys] = generate(; datasize = datasize, nkeys = nkeys)

        s0 = suite["nkeys=$nkeys"] = BenchmarkGroup()

        sbs = s0["alg=:base_seq"] = BenchmarkGroup()
        sbs["ntasks=1"] = @benchmarkable(
            # Base.Dict, sequential
            hist_seq!(dict, (CACHE[]::$T)[$nkeys]),
            setup = (dict = Dict{String,Int}()),
            evals = 1,
        )
        scs = s0["alg=:cdict_seq"] = BenchmarkGroup()
        scs["ntasks=1"] = @benchmarkable(
            # ConcurrentDict, sequential
            hist_seq!(dict, (CACHE[]::$T)[$nkeys]),
            setup = (dict = ConcurrentDict{String,Int}()),
            evals = 1,
        )
        sbp = s0["alg=:base_par"] = BenchmarkGroup()
        scp = s0["alg=:cdict_par"] = BenchmarkGroup()
        for ntasks in ntasks_list
            sbp["ntasks=$ntasks"] = @benchmarkable(
                # Base.Dict, parallel
                hist_parallel_dac((CACHE[]::$T)[$nkeys]; ntasks = $ntasks),
                evals = 1,
            )
            scp["ntasks=$ntasks"] = @benchmarkable(
                # ConcurrentDict, parallel
                hist_parallel!(dict, (CACHE[]::$T)[$nkeys]; ntasks = $ntasks),
                setup = (dict = ConcurrentDict{String,Int}()),
                evals = 1,
            )
        end
    end
    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
@testset "algorithms" begin

    number_of_lines = 10
    H1 = hyperplanes_poisson(1, number_of_lines)
    res1a = find_all_intersections(H1)
    res1b = intersections_naive(H1)

    @test length(res1a)==length(res1b)

end

max_stars_repo_name: SimonCoste/LineIntersections

max_issues_repo_name: SimonCoste/LineIntersections

max_forks_repo_name: SimonCoste/LineIntersections


text:
#=
Probe at randomly selected corpus points.
A fixed count of probes, with short circuit for perfect match.

Should be kept very similar to heuristicPatchMatch
TODO refactor so common code is in one place
=#
function matchPatchesAtRandomCorpusPoints(
    synthPatch,  # const IN
    targetImage, # mutable, IN and OUT
    corpusImage, # const IN
    bestKnownPatchDiff, # IN
    synthResult,    # mutable, OUT
    searchResult    # mutable, IN/OUT to be mutated by the search
    )

    # OLD setStartingSearchResult(searchResult, bestKnownPatchDiff)
    #@debug  "Random search diff to beat %f\n" result.bestProbeResult.patchDifference

    for i = 1:parameters.maxProbeCount
        corpusPatchCenterPoint = generateRandomMaskedPoint(corpusImage)
        probeResult = patchMatch(
            targetImage, corpusImage,
            corpusPatchCenterPoint,
            synthPatch,
            searchResult.bestProbeResult.patchDifference)

        if probeResult.betterment != NotBetter
            # Better or equal to any point probed in this pass

            # Call method of SearchResult to lift the probeResult into the searchResult
            setBetterSearchResult(searchResult, probeResult, corpusPatchCenterPoint)

            @debug  "    Better random probe, at" probeResult corpusPatchCenterPoint

            # Is better, might be perfect
            if probeResult.betterment == PerfectMatch
                @debug "Perfect match"
                break
            end
        end
    end

    @debug  "  >>Random search result " searchResult
    return searchResult
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


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

import DotEnv
DotEnv.config()

using MOSInterface

# Interface
interface = Interface()

# Delete model
delete_model_with_name(interface, "Stochastic Unit Commitment")

# New model
model = new_model(interface, "./examples/jump/unit_commitment/unit_commitment_model.jl")

# Get model
model = get_model_with_name(interface, "Stochastic Unit Commitment")

# Set inputs
set_interface_object(model, "n_g", 10)
set_interface_object(model, "n_gsl", 5)

# Show info
println(get_name(model))
@assert(get_system(model) == "jump")
println(get_status(model))

# Solve
MOSInterface.run(model)
println(get_status(model))

max_stars_repo_name: Fuinn/mos-examples

max_issues_repo_name: Fuinn/mos-examples

max_forks_repo_name: Fuinn/mos-examples


text:
using Oceananigans.Fields: FieldSlicer
using Oceananigans.Diagnostics

using Oceananigans.Models.HydrostaticFreeSurfaceModels: VectorInvariant

function nan_checker_aborts_simulation(arch)
    grid = RegularRectilinearGrid(size=(4, 2, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch)
    simulation = Simulation(model, Δt=1, stop_iteration=1)

    model.velocities.u[1, 1, 1] = NaN

    run!(simulation)

    return nothing
end

TestModel_VerticallyStrectedRectGrid(arch, FT, ν=1.0, Δx=0.5) =
    IncompressibleModel(
          grid = VerticallyStretchedRectilinearGrid(FT, architecture = arch, size=(3, 3, 3), x=(0, 3Δx), y=(0, 3Δx), z_faces=0:Δx:3Δx,),
       closure = IsotropicDiffusivity(FT, ν=ν, κ=ν),
  architecture = arch,
    float_type = FT
)

TestModel_RegularRectGrid(arch, FT, ν=1.0, Δx=0.5) =
    IncompressibleModel(
          grid = RegularRectilinearGrid(FT, topology=(Periodic, Periodic, Periodic), size=(3, 3, 3), extent=(3Δx, 3Δx, 3Δx)),
       closure = IsotropicDiffusivity(FT, ν=ν, κ=ν),
  architecture = arch,
    float_type = FT
)

function diagnostic_windowed_spatial_average(arch, FT)
    model = TestModel_RegularRectGrid(arch, FT)
    set!(model.velocities.u, 7)
    slicer = FieldSlicer(i=model.grid.Nx÷2:model.grid.Nx, k=1)
    u_mean = WindowedSpatialAverage(model.velocities.u; dims=(1, 2), field_slicer=slicer)
    return CUDA.@allowscalar u_mean(model)[1] == 7
end

function diffusive_cfl_diagnostic_is_correct(arch, FT)
    Δt = FT(1.3e-6)
    Δx = FT(0.5)
    ν = FT(1.2)
    CFL_by_hand = Δt * ν / Δx^2

    model = TestModel_RegularRectGrid(arch, FT, ν, Δx)
    cfl = DiffusiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

function advective_cfl_diagnostic_is_correct_on_regular_grid(arch, FT)
    model = TestModel_RegularRectGrid(arch, FT)

    Δt = FT(1.3e-6)
    Δx = FT(model.grid.Δx)
    u₀ = FT(1.2)
    CFL_by_hand = Δt * u₀ / Δx

    set!(model, u=u₀)
    cfl = AdvectiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

function advective_cfl_diagnostic_is_correct_on_vertically_stretched_grid(arch, FT)
    model = TestModel_VerticallyStrectedRectGrid(arch, FT)

    Δt = FT(1.3e-6)
    Δx = FT(model.grid.Δx)
    u₀ = FT(1.2)
    CFL_by_hand = Δt * u₀ / Δx

    set!(model, u=u₀)
    cfl = AdvectiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

function accurate_advective_cfl_on_regular_grid(arch, FT)
    model = TestModel_RegularRectGrid(arch, FT)

    Δt = FT(1.7)

    Δx = model.grid.Δx
    Δy = model.grid.Δy
    Δz = model.grid.Δz

    u₀ = FT(1.2)
    v₀ = FT(-2.5)
    w₀ = FT(3.9)

    set!(model, u=u₀, v=v₀, w=w₀)

    CFL_by_hand = Δt * (abs(u₀) / Δx + abs(v₀) / Δy + abs(w₀) / Δz)

    cfl = CFL(FT(Δt), Oceananigans.Diagnostics.accurate_cell_advection_timescale)

    return cfl(model) ≈ CFL_by_hand
end

function accurate_advective_cfl_on_stretched_grid(arch, FT)
    grid = VerticallyStretchedRectilinearGrid(architecture=arch, size=(4, 4, 8), x=(0, 100), y=(0, 100), z_faces=[k^2 for k in 0:8])
    model = IncompressibleModel(grid=grid, architecture=arch)

    Δt = FT(15.5)

    Δx = model.grid.Δx
    Δy = model.grid.Δy

    # At k = 1, w = 0 so the CFL constraint happens at the second face (k = 2).
    Δz_min = CUDA.@allowscalar Oceananigans.Operators.Δzᵃᵃᶠ(1, 1, 2, grid)

    u₀ = FT(1.2)
    v₀ = FT(-2.5)
    w₀ = FT(3.9)

    set!(model, u=u₀, v=v₀, w=w₀, enforce_incompressibility=false)

    CFL_by_hand = Δt * (abs(u₀) / Δx + abs(v₀) / Δy + abs(w₀) / Δz_min)

    cfl = CFL(FT(Δt), Oceananigans.Diagnostics.accurate_cell_advection_timescale)

    return cfl(model) ≈ CFL_by_hand
end

function accurate_advective_cfl_on_lat_lon_grid(arch, FT)
    grid = RegularLatitudeLongitudeGrid(size=(8, 8, 8), longitude=(-10, 10), latitude=(0, 45), z=(-1000, 0))
    model = HydrostaticFreeSurfaceModel(architecture=arch, grid=grid, momentum_advection=VectorInvariant())

    Δt = FT(1000)

    Nx, Ny, Nz = size(grid)

    # Will be the smallest at higher latitudes.
    Δx_min = CUDA.@allowscalar Oceananigans.Operators.Δxᶠᶜᵃ(1, Ny, 1, grid)

    # Will be the same at every grid point.
    Δy_min = CUDA.@allowscalar Oceananigans.Operators.Δyᶜᶠᵃ(1, 1, 1, grid)

    Δz = model.grid.Δz

    u₀ = FT(1.2)
    v₀ = FT(-2.5)
    w₀ = FT(-0.1)

    # Avoid recomputation of w
    set!(model.velocities.u, u₀)
    set!(model.velocities.v, v₀)
    set!(model.velocities.w, w₀)

    CFL_by_hand = Δt * (abs(u₀) / Δx_min + abs(v₀) / Δy_min + abs(w₀) / Δz)

    cfl = CFL(FT(Δt), Oceananigans.Diagnostics.accurate_cell_advection_timescale)

    return cfl(model) ≈ CFL_by_hand
end

get_iteration(model) = model.clock.iteration
get_time(model) = model.clock.time

function diagnostics_getindex(arch, FT)
    model = TestModel_RegularRectGrid(arch, FT)
    simulation = Simulation(model, Δt=0, stop_iteration=0)
    nc = NaNChecker(model, schedule=IterationInterval(1), fields=model.velocities)
    simulation.diagnostics[:nc] = nc

    # The first diagnostic is the NaN checker.
    return simulation.diagnostics[2] == nc
end

function diagnostics_setindex(arch, FT)
    model = TestModel_RegularRectGrid(arch, FT)
    simulation = Simulation(model, Δt=0, stop_iteration=0)

    nc1 = NaNChecker(model, schedule=IterationInterval(1), fields=model.velocities)
    nc2 = NaNChecker(model, schedule=IterationInterval(2), fields=model.velocities)
    nc3 = NaNChecker(model, schedule=IterationInterval(3), fields=model.velocities)

    push!(simulation.diagnostics, nc1, nc2)
    simulation.diagnostics[2] = nc3

    return simulation.diagnostics[:diag2] == nc3
end

@testset "Diagnostics" begin
    @info "Testing diagnostics..."

    for arch in archs
        @testset "NaN Checker [$(typeof(arch))]" begin
            @info "  Testing NaN Checker [$(typeof(arch))]..."
            @test_throws ErrorException nan_checker_aborts_simulation(arch)
        end
    end

    for arch in archs
        @testset "CFL [$(typeof(arch))]" begin
            @info "  Testing CFL diagnostics [$(typeof(arch))]..."
            for FT in float_types
                @test diffusive_cfl_diagnostic_is_correct(arch, FT)
                @test advective_cfl_diagnostic_is_correct_on_regular_grid(arch, FT)
                @test advective_cfl_diagnostic_is_correct_on_vertically_stretched_grid(arch, FT)
                @test accurate_advective_cfl_on_regular_grid(arch, FT)
                @test accurate_advective_cfl_on_stretched_grid(arch, FT)
                @test accurate_advective_cfl_on_lat_lon_grid(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "Miscellaneous timeseries diagnostics [$(typeof(arch))]" begin
            @info "  Testing miscellaneous timeseries diagnostics [$(typeof(arch))]..."
            for FT in float_types
                @test diagnostic_windowed_spatial_average(arch, FT)
                @test diagnostics_getindex(arch, FT)
                @test diagnostics_setindex(arch, FT)
            end
        end
    end
end


text:
using PisingerKnapsack
using Test

include("Cwrapper.jl")


text:
using SymmetricSquares
using Test

@testset "SymmetricSquares" begin
    @testset "Symmetry group" begin
        include("symmetry.jl")
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

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

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

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using DiffEqCallbacks
using DynamicalSystemsBase
using OrdinaryDiffEq: Vern7, Tsit5
using LinearAlgebra, Statistics
using Test

println("\nTesting callbacks...")

@testset "SavingCallback parallel" begin

diffeq = (abstol=1e-14, reltol=1e-14, alg=Vern7(), maxiters=1e9)
ds = Systems.lorenz()
d0 = 1e-9
T = 100.0

save_func(u, t, integrator) = LinearAlgebra.norm(u[1] - u[2])
saved_values = SavedValues(eltype(ds.t0), eltype(ds.u0[1]))
cb = SavingCallback(save_func, saved_values)

u0 = get_state(ds)

pinteg = parallel_integrator(
    ds, [u0, u0 + rand(SVector{3})*d0];
    diffeq = (diffeq..., callback = cb)
)
step!(pinteg, T)
n = saved_values.saveval
t = saved_values.t
@test length(n) > 10000
# test that norm increases:
@test n[10000] > n[1]
@test length(pinteg.sol.u) == 1
end

@testset "SavingCallback tangent" begin

ds = Systems.lorenz()
d0 = 1e-9
T = 100.0

save_func(u, t, integrator) = LinearAlgebra.norm(get_deviations(integrator)[1, :])
saved_values = SavedValues(eltype(ds.t0), eltype(ds.u0[1]))
cb = SavingCallback(save_func, saved_values)

diffeq = (abstol=1e-14, reltol=1e-14, alg=Tsit5(), callback = cb)

u0 = get_state(ds)
pinteg = tangent_integrator(ds; diffeq)
step!(pinteg, T)
n = saved_values.saveval
t = saved_values.t
@test length(n) > 1000
# test that norm increases:
@test n[end] > n[5]
@test length(pinteg.sol.u) == 1

end

#=
@testset "ManifoldProjection" begin
  ds1 = Systems.henonheiles() #with Jac
  ds2 = CDS(ds1.prob) #without Jac
  @inline Vhh(q1, q2) = 1//2 * (q1^2 + q2^2 + 2q1^2 * q2 - 2//3 * q2^3)
  @inline Thh(p1, p2) = 1//2 * (p1^2 + p2^2)
  @inline Hhh(q1, q2, p1, p2) = Thh(p1, p2) + Vhh(q1, q2)
  @inline Hhh(u::AbstractVector) = Hhh(u...)

  tra1 = trajectory(ds1, 100.0)
  tra2 = trajectory(ds2, 100.0)

  E1 = [Hhh(p) for p in tra1]
  E2 = [Hhh(p) for p in tra2]

  @test std(E1) < 1e-12
  @test std(E2) < 1e-12
end
=#


text:
"""
    unsafe_string_or_null(ptr::Cstring) -> Union{String, Missing}

Convert a `Cstring` to a `Union{String, Missing}`, returning `missing` if the pointer is
`C_NULL`.
"""
function unsafe_string_or_null(ptr::Cstring)::Union{String, Missing}
    ptr == C_NULL ? missing : unsafe_string(ptr)
end


text:
# Define constants for speed:
include("define_constants.jl")
# Include definition of Transit structure type:
include("transit_structure.jl")
# Include code which computes linear limb-darkening term:
include("s2.jl")
# Include code which computes I_v, J_v, and derivatives wrt k:
include("IJv_derivative_struct.jl")

# Computes the coefficient for the uniform disk case, S_0:
function compute_uniform!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; r2=r*r; b2=b*b
# Uniform disk case:
# Compute sn[1] and its derivatives:
if b <= 1-r  # k^2 > 1
  t.sn[1] = pi*(1-r2)
  if t.grad
    t.dsndr[1] = -2*pi*r
    t.dsndb[1] = 0.
  end
  t.kap0 = convert(T,pi); t.kck = zero(T)
else
  # Twice area of kite-shaped region connecting centers of circles & intersection points:
  t.kite_area2 = sqrt(sqarea_triangle(one(T),b,r))
  # Angle of section for occultor:
  t.kap0  = atan(t.kite_area2,(r-1)*(r+1)+b2)
  # Angle of section for source:
  t.pimkap1 = atan(t.kite_area2,(r-1)*(r+1)-b2)
  # Flux of visible uniform disk:
  t.sn[1] = t.pimkap1 - r2*t.kap0 + .5*t.kite_area2
  if t.grad
    t.dsndr[1]= -2*r*t.kap0
    t.dsndb[1]= t.kite_area2/b
  end
  t.kck = t.kite_area2*t.fourbrinv
end
return
end

function sqarea_triangle(a,b::T,c::T) where {T <: Real}
# Function which computes sixteen times the square of the area
# of a triangle with sides a, b and c using Kahan method.
# How to compute (sixteen times the) area squared of triangle with
# high precision (Goldberg 1991).
# First, do a quick sort of three numbers:
if c > b
  tmp = c; c = b; b = tmp
end
if c > a
  tmp = a; a = c; c = tmp
end
if b > a
  tmp = b; b = a; a = tmp
end
area = (a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c))
return area
end

# Computes a limb-darkened transit light curve with the dependence:
# I(\mu) = 1-\sum_{n=1}^N u_n (1-\mu)^n
# where \mu = \cos{\theta} = z is the cosine of the angle
# angle from the sub-stellar point (or equivalently the
# height on the star relative to the sky plane if the radius
# of the star is unity.

function transit_poly_c(t::Transit_Struct{T}) where {T <: Real}
# Number of limb-darkening components to include (beyond 0 and 1):
# We are parameterizing these with the function:
# g_n = c_n [(n+2) z^n - n z^{n-2}] for n >= 2
# while g_{0} = c_0 z^0 (uniform source) and g_{1} = c_1 z^1 (linear limb-darkening)
# which gives a Green's integral of:
# P(G_n) = \int_{\pi-\phi}^{2\pi+\phi} (1-r^2-b^2-2br s_\varphi)^{n/2} (r+b s_\varphi) d\varphi
# for which we have a solution in terms of I_v (for even n) and J_v (for odd n).
# The variable "t" is a structure which contains transit parameters
# and intermediate quantities computed from these:
r=t.r; b=t.b; n = t.n; r2=r*r
# Set up a vector for storing results of P(G_n)-Q(G_n); note that
# this is a different vector than the Starry case:

# Check for different cases:
if b >= 1+r
  # unobscured - return one:
  return one(T)
end
if r >= 1+b
  # full obscuration - return zero:
  return zero(T)
end
if b == 0.0
  # Annular eclipse - integrate around the full boundary of both bodies:
  onemr2 = 1-r2
  flux = zero(T); t.sqrt1mr2 = sqrt(onemr2)
  flux = (t.c_n[1]*onemr2+t.twothird*t.c_n[2]*t.sqrt1mr2^3)
  fac= 2r2*onemr2
  @inbounds for i=2:t.n
    flux += -t.c_n[i+1]*fac
    fac *= t.sqrt1mr2
  end
  return flux*pi*t.den
else
# Next, compute k^2 = m:
  t.onembmr2=(r+1-b)*(1-r+b); t.fourbr = 4b*r; t.fourbrinv = inv(t.fourbr)
#  t.sqbr = sqrt(b*r); t.sqbrinv = inv(t.sqbr)
  t.onembmr2inv=inv(t.onembmr2); t.sqonembmr2 = sqrt(t.onembmr2)
  t.onembpr2 = (1-r-b)*(1+b+r)
#  t.k2 = t.onembmr2*t.fourbrinv
  t.k2 = t.onembpr2*t.fourbrinv+1
  if t.k2 > 0
    t.k = sqrt(t.k2)
  else
    println("negative k2: ",t.k2," r: ",r," b: ",b)
    t.k2 = 0.0; t.k = 0.0
  end
  if t.k2 > 1
    if t.k2 > 2.0
      t.kc2 = 1.0-inv(t.k2)
      t.kc = sqrt(t.kc2)
    else
      t.kc2 = t.onembpr2/((1-b+r)*(1-r+b))
      t.kc = sqrt(t.kc2)
    end
  else
    if t.k2 > 0.5
      t.kc2 = (r-1+b)*(b+r+1)*t.fourbrinv
      t.kc = sqrt(t.kc2)
    else
      t.kc2 = 1.0-t.k2
      t.kc = sqrt(t.kc2)
    end
  end
end

# Compute uniform case, sn[1]:
compute_uniform!(t)
# If uniform model, then return:
if t.n == 0
  return t.c_n[1]*t.sn[1]*t.den
end

# Compute linear case, sn[2]:
#s2!(t)
t.sn[2],t.Eofk,t.Em1mKdm = s2_ell(r,b)
#sn2,Eofk,Em1mKdm = s2_ell(r,b)
# Now, compare results:
#if abs(t.sn[2]-sn2) > 1e-8*abs(sn2)
#  println("r: ",r," b: ",b," sn[2]: ",t.sn[2]," sn2: ",sn2)
#end
# Special case of linear limb-darkening:
if t.n == 1
  flux = t.c_n[1]*t.sn[1]+t.c_n[2]*t.sn[2]
  flux *= t.den  # for c_2 and above, the flux is zero.
  return flux
end

# Special case of quadratic limb-darkening:
if t.n == 2
# Transformed expressions from Mandel & Agol:
  r2=r*r; b2=b*b
  eta2 = 0.5*r2*(r2+2*b2)
  if t.k2 > 1
    four_pi_eta = 4pi*eta2
  else
    four_pi_eta = 2*(pi-t.pimkap1+2*eta2*t.kap0-0.25*t.kite_area2*(1.0+5r2+b2))
  end
  t.sn[3] = 2*(t.sn[1]-pi)+four_pi_eta
  flux = t.c_n[1]*t.sn[1]+t.c_n[2]*t.sn[2]+t.c_n[3]*t.sn[3]
  flux *= t.den  # for c_2 and above, the flux is zero.
  return flux
end

# Compute the J_v and I_v functions:
if t.k2 > 0
  if (t.k2 < 0.5 || t.k2 > 2.0) && t.v_max > 3
# This computes I_v,J_v for the largest v, and then works down to smaller values:
#    IJv_lower!(t.k2,t.kck,t.kc,t.kap,t.Eofk,t.Em1mKdm,t)
    IJv_lower!(t)
  else
# This computes I_0,J_0,J_1, and then works upward to larger v:
#    IJv_raise!(t.k2,t.kck,t.kc,t.kap,t.Eofk,t.Em1mKdm,t)
    IJv_raise!(t)
  end
end

# Add up first two terms in flux numerator:
flux = t.c_n[1]*t.sn[1]+t.c_n[2]*t.sn[2]
# Next, loop over the Green's function components:
@inbounds for n=2:t.n
  pofgn = zero(T)
  if iseven(n)
# For even values of n, sum over I_v:
    n0 = convert(Int64,n/2)
    coeff = (-t.fourbr)^n0
    # Compute i=0 term
    pofgn = coeff*((r-b)*t.Iv[n0+1]+2b*t.Iv[n0+2])
# For even n, compute coefficients for the sum over I_v:
    @inbounds for i=1:n0
      coeff *= -(n0-i+1)/i*t.k2
      pofgn += coeff*((r-b)*t.Iv[n0-i+1]+2b*t.Iv[n0-i+2])
    end
    pofgn *= 2r
  else
# Now do the same for odd N_c in sum over J_v:
    n0 = convert(Int64,(n-3)/2)
    coeff = (-t.fourbr)^n0
    # Compute i=0 term
    pofgn = coeff*((r-b)*t.Jv[n0+1]+2b*t.Jv[n0+2])
# For odd n, compute coefficients for the sum over J_v:
    @inbounds for i=1:n0
      coeff *= -(n0-i+1)/i*t.k2
      pofgn += coeff*((r-b)*t.Jv[n0-i+1]+2b*t.Jv[n0-i+2])
    end
    pofgn *= 2r*t.onembmr2*t.sqonembmr2
  end
# Q(G_n) is zero in this case since on limb of star z^n = 0 at the stellar
# boundary for n > 0.
# Compute sn[n]:
  t.sn[n+1] = -pofgn
  flux += t.c_n[n+1]*t.sn[n+1]
end
# That's it!
# flux = sum(t.c_n.*t.sn)/(pi*(t.c_n[1]+t.twothird*t.c_n[2]))  # for c_2 and above, the flux is zero.
# Divide by denominator:
flux *= t.den  # for c_2 and above, the flux is zero.
return flux
end
# That's it!

function transit_poly!(r::T,b::T,u_n::Array{T,1},dfdrbu::Array{T,1}) where {T <: Real}
t = transit_init(r,b,u_n,true)
# Pass c_n (without last two dummy values):
flux = transit_poly_c!(t)
# Now, transform derivaties from c to u:
fill!(dfdrbu,zero(T))
dfdrbu[1] = t.dfdrbc[1]  # r derivative
dfdrbu[2] = t.dfdrbc[2]  # b derivative
# u_n derivatives:
@inbounds for i=1:t.n, j=0:t.n
  dfdrbu[i+2] += t.dfdrbc[j+3]*t.dcdu[j+1,i]
end
return flux
end

function transit_poly(r::T,b::T,u_n::Array{T,1}) where {T <: Real}
t = transit_init(r,b,u_n,false)
# Pass c_n (without last two dummy values):
return transit_poly_c(t)
end

function transit_poly!(t::Transit_Struct{T}) where {T <: Real}
# This function assumes that c_n has already been computed from u_n
# (this can be used to save compute time when limb-darkening is fixed for
# a range of radii/impact parameters).
# Pass transit structure, and compute flux:
if t.grad
  flux = transit_poly_c!(t)
  # Now, transform derivaties from c to u:
  fill!(t.dfdrbu,zero(T))
  t.dfdrbu[1] = t.dfdrbc[1]  # r derivative
  t.dfdrbu[2] = t.dfdrbc[2]  # b derivative
  @inbounds for i=1:t.n, j=0:t.n
    t.dfdrbu[i+2] += t.dfdrbc[j+3]*t.dcdu[j+1,i]
  end
  return flux
else
  return transit_poly_c(t)
end
return
end

function transit_poly_c!(t::Transit_Struct{T}) where {T <: Real}
r = t.r; b=t.b; n = t.n; r2 =r*r; b2=b*b
@assert((length(t.c_n)+2) == length(t.dfdrbc))
@assert(r > 0)
# Number of limb-darkening components to include (beyond 0 and 1):
# We are parameterizing these with the function:
# g_n = c_n [(n+2) z^n - n z^{n-2}] for n >= 2
# while g_{0} = c_0 z^0 (uniform source) and g_{1} = c_1 z^1 (linear limb-darkening)
# which gives a Green's integral of:
# P(G_n) = \int_{\pi-\phi}^{2\pi+\phi} (1-r^2-b^2-2br s_\varphi)^{n/2} (r+b s_\varphi) d\varphi
# for which we have a solution in terms of I_v (for even n) and J_v (for odd n).
# Compute the derivative of the flux with respect to the different coefficients.

# Set up a vector for storing results of P(G_n)-Q(G_n); note that
# this is a different vector than the Starry case:
# Check for different cases:
if b >= 1+r || r ==  0.0
  # unobscured - return one, and zero derivatives:
  fill!(t.dfdrbc,zero(T))
  return one(T)
end
if r >= 1+b
  # full obscuration - return zero, and zero derivatives:
  fill!(t.dfdrbc,zero(T))
  return zero(T)
end
if b == 0.0
  # Annular eclipse - integrate around the full boundary of both bodies:
  flux = zero(T); onemr2 = 1-r2; t.sqrt1mr2 = sqrt(onemr2)
  fill!(t.dfdrbc,zero(T))
  flux = (t.c_n[1]*onemr2+t.twothird*t.c_n[2]*t.sqrt1mr2^3)*pi*t.den
  fac  = 2r2*onemr2*pi*t.den
  facd = -2r*pi*t.den
  t.dfdrbc[1] = t.c_n[1]*facd + t.c_n[2]*facd*t.sqrt1mr2
  @inbounds for i=2:t.n
    flux -= t.c_n[i+1]*fac
    t.dfdrbc[1] += t.c_n[i+1]*facd*(2*onemr2-i*r2)
    t.dfdrbc[i+3] -= fac
    fac *= t.sqrt1mr2
    facd *= t.sqrt1mr2
  end
  #  dfdrbc[2]=0 since the derivative with respect to b is zero.
  t.dfdrbc[3] = (onemr2-flux)*pi*t.den
  t.dfdrbc[4] = t.twothird*(t.sqrt1mr2^3-flux)*pi*t.den
  # Also need to compute derivatives [ ]
  return flux
else
# Next, compute k^2 = m:
  t.onembmr2=(r-b+1)*(1-r+b); t.fourbr = 4b*r; t.fourbrinv = inv(t.fourbr)
#  t.sqbr = sqrt(b*r); t.sqbrinv = inv(t.sqbr)
  t.onembmr2inv=inv(t.onembmr2); t.sqonembmr2 = sqrt(t.onembmr2)
  t.onembpr2 = (1-r-b)*(1+b+r)
  t.k2 = t.onembmr2*t.fourbrinv;
  if t.k2 > 0
    t.k = sqrt(t.k2)
  else
    println("negative k2: ",t.k2," r: ",r," b: ",b)
    t.k2 = 0.0; t.k = 0.0
  end
  dkdr = (b2-r2-1)/(8*t.k*b*r2)
  dkdb = (r2-b2-1)/(8*t.k*b2*r)
  if t.k2 > 1
    if t.k2 > 2.0
      t.kc2 = 1.0-inv(t.k2)
      t.kc = sqrt(t.kc2)
    else
      t.kc2 = t.onembpr2/((1+r-b)*(1-r+b))
      t.kc = sqrt(t.kc2)
    end
  else
    if t.k2 > 0.5
      t.kc2 = (r-1+b)*(b+r+1)*t.fourbrinv
      t.kc = sqrt(t.kc2)
    else
      t.kc2 = (r-1+b)*(b+r+1)*t.fourbrinv
      t.kc = sqrt(t.kc2)
    end
  end
end

# Compute uniform case:
compute_uniform!(t)

# Compute the highest value of v in J_v or I_v that we need:
# Compute sn[2] and its derivatives:
#s2!(t)
if t.n >= 1
  t.sn[2],t.Eofk,t.Em1mKdm = s2!(r,b,t.s2_grad)
  t.dsndr[2] = t.s2_grad[1]
  t.dsndb[2] = t.s2_grad[2]
end

# Special case of quadratic limb-darkening:
if t.n >= 2
  if t.n == 2
# Transformed expressions from Mandel & Agol:
    r2pb2 = (r2+b2)
    eta2 = 0.5*r2*(r2pb2+b2)
    deta2dr =  2*r*r2pb2
    deta2db = 2*b*r2
    if t.k2 > 1
      four_pi_eta = 4pi*(eta2-0.5)
      detadr = 4pi*deta2dr
      detadb = 4pi*deta2db
    else
      four_pi_eta = 2*(-t.pimkap1+2*eta2*t.kap0-0.25*t.kite_area2*(1.0+5r2+b2))
      detadr = 8r*(r2pb2*t.kap0-t.kite_area2)
      detadb = 2/b*(4*b2*r2*t.kap0-(1+r2pb2)*t.kite_area2)
    end
    t.sn[3] = 2*t.sn[1]+four_pi_eta
    t.dsndr[3] = 2*t.dsndr[1]+detadr
    t.dsndb[3] = 2*t.dsndb[1]+detadb
  else
    # Compute the J_v and I_v functions:
    if t.k2 > 0
      if (t.k2 < 0.5 || t.k2 > 2.0) && t.v_max > 3
    # This computes I_v,J_v for the largest v, and then works down to smaller values:
    #    dIJv_lower_dk!(t.k2,t.kck,t.kc,t.kap,t.Eofk,t.Em1mKdm,t)
        dIJv_lower_dk!(t)
      else
    # This computes I_0,J_0,J_1, and then works upward to larger v:
      #dIJv_raise_dk!(t.k2,t.kck,t.kc,t.kap,t.Eofk,t.Em1mKdm,t)
        dIJv_raise_dk!(t)
      end
    end

    rinv = inv(r); binv = inv(b); rmb_on_onembmr2=(r-b)*t.onembmr2inv
    # Next, loop over the Green's function components:
    Iv1 = zero(T); Iv2=zero(T); Jv1=zero(T); Jv2=zero(T)
    nmi = zero(Int64); fac1 = zero(T)
    @inbounds for n=2:t.n
      pofgn = zero(T)
      dpdr = zero(T)
      dpdb = zero(T)
      dpdk = zero(T)
      if iseven(n)
    # For even values of n, sum over I_v:
        n0 = convert(Int64,n/2)
        coeff = (-t.fourbr)^n0
        # Compute i=0 term
        Iv1 = t.Iv[n0+1]; Iv2 = t.Iv[n0+2]
        pofgn = coeff*((r-b)*Iv1+2b*Iv2)
        dpdr = coeff*Iv1
        dpdb = coeff*(-Iv1+2*Iv2)
        dpdr += (n0+1)*pofgn*rinv
        dpdb += n0*pofgn*binv
        k2n = coeff
    # For even n, compute coefficients for the sum over I_v:
        @inbounds for i=1:n0
          nmi = n0-i+1
          Iv2 = Iv1; Iv1 = t.Iv[nmi]
          k2n *= -t.k2
          coeff = t.bincoeff[i+1,n0]*k2n
          term =  coeff*((r-b)*Iv1+2b*Iv2)
          pofgn += term
          dpdr += coeff*Iv1
          dpdb += coeff*(-Iv1+2.0*Iv2)
          fac1 = i*2.0*rmb_on_onembmr2
          dpdr += term*(-fac1+nmi*rinv)
          dpdb += term*( fac1+(nmi-1.0)*binv)
        end
        pofgn *= 2r
        dpdr *= 2r
        dpdb *= 2r
      else
    # Now do the same for odd N_c in sum over J_v:
        n0 = convert(Int64,(n-3)/2)
        coeff = (-t.fourbr)^n0
        # Compute i=0 term
        Jv1 = t.Jv[n0+1]; Jv2 = t.Jv[n0+2]
        pofgn = coeff*((r-b)*Jv1+2b*Jv2)
        dpdr = coeff*Jv1
        dpdb = coeff*(-Jv1+2*Jv2)
        dpdr  += pofgn*(-3*rmb_on_onembmr2+(n0+1)*rinv)
        dpdb  += pofgn*( 3*rmb_on_onembmr2+n0*binv)
        dpdk = coeff*((r-b)*t.dJvdk[n0+1]+2b*t.dJvdk[n0+2])
    # For even n, compute coefficients for the sum over I_v:
        k2n = coeff
        @inbounds for i=1:n0
          nmi = n0-i+1
          #coeff *= -nmi/i*t.k2
          k2n *= -t.k2
          coeff = t.bincoeff[i+1,n0]*k2n
          Jv2 = Jv1; Jv1 = t.Jv[nmi]
          term = coeff*((r-b)*Jv1+2b*Jv2)
          pofgn += term
          dpdr  +=  coeff*Jv1
          dpdb  +=  coeff*(-Jv1+2.0*Jv2)
          fac1 = (i*2.0+3.0)*rmb_on_onembmr2
          dpdr  += term*(-fac1+nmi*rinv)
          dpdb  += term*( fac1+(nmi-1.0)*binv)
          dpdk  += coeff*((r-b)*t.dJvdk[nmi]+2b*t.dJvdk[nmi+1])
        end
        norm = 2r*t.onembmr2*t.sqonembmr2
        pofgn *= norm
        dpdr  *= norm
        dpdb  *= norm
        dpdk  *= norm
      end
    # Q(G_n) is zero in this case since on limb of star z^n = 0 at the stellar
    # boundary for n > 0.
    # Compute sn[n]:
      t.sn[n+1] = -pofgn
      t.dsndr[n+1] = -(dpdr+dpdk*dkdr)
      t.dsndb[n+1] = -(dpdb+dpdk*dkdb)
    end
  end
end
# That's it!
# Compute derivatives with respect to the coefficients:
flux = zero(T)
t.dfdrbc[1]=zero(T)  # Derivative with respect to r
t.dfdrbc[2]=zero(T)  # Derivative with respect to b
@inbounds for i=0:t.n
  # derivatives with respect to the coefficients:
  t.dfdrbc[i+3]= t.sn[i+1]*t.den
  # total flux:
  flux += t.c_n[i+1]*t.dfdrbc[i+3]
  # derivatives with respect to r and b:
  t.dfdrbc[1] += t.c_n[i+1]*t.dsndr[i+1]*t.den
  t.dfdrbc[2] += t.c_n[i+1]*t.dsndb[i+1]*t.den
end
# Include derivatives with respect to first two c_n parameters:
t.dfdrbc[3] -= flux*t.den*pi
t.dfdrbc[4] -= flux*t.den*pi*t.twothird
return flux
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
@testset "64.minimum-path-sum.jl" begin
    @test min_path_sum([[1, 3, 1], [1, 5, 1], [4, 2, 1]]) == 7
    @test min_path_sum([[1, 2, 3], [4, 5, 6]]) == 12
    @test min_path_sum([[0]]) == 0
    @test min_path_sum([[1, 2, 5], [3, 2, 1]]) == 6
end


text:
# variational inference sample code

module VariationalInference
  using DataFrames, CSV, StatsBase, Plots, Distributions
  using SpecialFunctions
  pyplot()

  function main(;test=false)
    # read original data
    data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
    all_data = CSV.read(data_path, DataFrame,
                        header=["date", "time", "ir", "z"],
                        delim=' ')
    data = sort(sample(all_data.z, 3000))
    data = DataFrame(z=data)

    # initialize responsibility
    K = 3 # num of class
    n = Int(ceil(size(data)[1]/K))
    for k in 1:K
      data_k = []
      for i in 1:size(data)[1]
        if k == Int(ceil(i/n))
          append!(data_k, 1.0)
        else
          append!(data_k, 0.0)
        end
      end

      if k == 1
        data = hcat(data, DataFrame(one=data_k))
      elseif k == 2
        data = hcat(data, DataFrame(two=data_k))
      else
        data = hcat(data, DataFrame(three=data_k))
      end
    end

    # iteration
    params_history = Dict()
    params = DataFrame()
    for t in 1:10000
      data, params = one_step(data, K)
      # record per 10 times
      if t % 10 == 0
        params_history[t] = params
        println("$(t) $(params)")

        if length(params_history) < 2
          continue
        end

        if all([(abs(params_history[t-10][k, "tau"] - params_history[t][k, "tau"]) < 10e-5) for k in 1:K])
          break
        end
      end
    end
    draw_sensor_dist(data, params, K,
                     "src/variational_inference/draw_sensor_dist_K3_t340.png",
                     test)
  end

  function update_parameters(ds, k, mu_avg=600, zeta=1, alpha=1, beta=1, tau=1)
    if k == 1
      key_class = "one"
    elseif k == 2
      key_class = "two"
    else
      key_class = "three"
    end

    # calculate R, S, T
    R = sum([ds[i, key_class] for i in 1:size(ds)[1]])
    S = sum([ds[i, key_class]*ds[i, :z] for i in 1:size(ds)[1]])
    T = sum([ds[i, key_class]*(ds[i, :z]^2) for i in 1:size(ds)[1]])

    # store calculated parameters temporarily
    hat = Dict()
    hat["tau"] = R + tau
    hat["zeta"] = R + zeta
    hat["mu_avg"] = (S + zeta*mu_avg)/hat["zeta"]
    hat["alpha"] = R/2 + alpha
    hat["beta"] = (T + zeta*(mu_avg^2) - hat["zeta"]*(hat["mu_avg"]^2))/2 + beta
    hat["z_std"] = sqrt(hat["beta"]/hat["alpha"])

    return DataFrame(hat)
  end

  function draw_sensor_dist(data, params, K, save_name, test)
    pi_ = rand(Dirichlet([params[k, "tau"] for k in 1:K]))
    pdfs = [Normal(params[k, "mu_avg"], params[k, "z_std"]) for k in 1:K]

    xs = range(600, 650, length=Int64(floor(50/0.5)))

    # histogram of observation
    histogram(data.z, bins=(maximum(data.z)-minimum(data.z)), color=:gray, label="histogram")
    # draw gaussian mixture distribution
    ys = [sum([pdf(pdfs[k], x)*pi_[k] for k in 1:K])*size(data)[1] for x in xs]
    plot!(xs, ys, color=:red, label="Mixture")
    # gaussian distribution of each cluster
    for k in 1:K
      ys = [pdf(pdfs[k], x)*pi_[k]*size(data)[1] for x in xs]
      plot!(xs, ys, color=:blue, label="Cluster$(k)")
    end
    if test == false
      save_path = joinpath(split(@__FILE__, "src")[1], save_name)
      savefig(save_path)
    end
  end

  function responsibility(z, K, params)
    tau_sum = sum([params[k, "tau"] for k in 1:K])

    r = Dict()
    for k in 1:K
      log_rho = (digamma(params[k, "alpha"]) - log(params[k, "beta"]))/2 - (1/params[k, "zeta"] + ((params[k, "mu_avg"]-z)^2)*params[k, "alpha"]/params[k, "beta"])/2 + digamma(params[k, "tau"]) - digamma(tau_sum)
      if k == 1
        r["one"] = exp(log_rho)
      elseif k == 2
        r["two"] = exp(log_rho)
      else
        r["three"] = exp(log_rho)
      end

    end

    # sum of r_k
    rks = []
    for k in 1:K
      if k == 1
        append!(rks, r["one"])
      elseif k == 2
        append!(rks, r["two"])
      else
        append!(rks, r["three"])
      end
    end
    w = sum(rks)

    # normalization
    for k in 1:K
      if k == 1
        r["one"] /= w
      elseif k == 2
        r["two"] /= w
      else
        r["three"] /= w
      end
    end

    return r
  end

  function one_step(data, K)
    # update parameters
    params = DataFrame()
    for k in 1:K
      if k == 1
        params = update_parameters(data, k)
      else
        params = vcat(params, update_parameters(data, k))
      end
    end

    # update responsibility
    rs = [responsibility(data[i, :z], K, params) for i in 1:size(data)[1]]
    for k in 1:K
      if k == 1
        data.one = [rs[i]["one"] for i in 1:size(data)[1]]
      elseif k == 2
        data.two = [rs[i]["two"] for i in 1:size(data)[1]]
      else
        data.three = [rs[i]["three"] for i in 1:size(data)[1]]
      end
    end

    return data, params
  end

end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy