Resampling.jl

## --- Bootstrap resampling

    # Kernel functions for bsr!
    @inline function uniform(rng::AbstractRNG, mu, halfwidth)
        mu + halfwidth * (2*rand(rng)-1)
    end
    export uniform
    @inline function triangular(rng::AbstractRNG, mu, halfwidth)
        mu + halfwidth * (rand(rng) - rand(rng))
    end
    export triangular
    @inline function gaussian(rng::AbstractRNG, mu, sigma)
        mu + sigma * randn(rng)
    end
    export gaussian

    """
    ```julia
    bsr!([f::Function=gaussian], resampled::Array, index::Vector{Int}, data, sigma, p;
        \trng::AbstractRNG=MersenneTwister()
    )
    ```

    Fill `resampled` with data boostrap resampled from a (sample-per-row / element-per-column)
    dataset `data` with uncertainties `sigma` and resampling probabilities `p`, optionally using
    random numbers generated by `f` where `f` is a function of the form `f(rng, data[i], sigma[i])`
    """
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractVecOrMat, sigma::Number, p::Number; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma)
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractVecOrMat, sigma::Number, p::AbstractVector; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma)
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end

    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractMatrix, sigma::AbstractVector, p::Number; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma[j])
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractMatrix, sigma::AbstractVector, p::AbstractVector; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma[j])
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end

    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractArray, sigma::AbstractArray, p::Number; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma[row,j])
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractArray, sigma::AbstractArray, p::AbstractVector; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma[row,j])
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    # default method if f not specified
    bsr!(resampled::Array, index::Vector{Int}, data::AbstractArray, sigma, p; rng=MersenneTwister()) = bsr!(gaussian, resampled, index, data, sigma, p, rng=rng)
    export bsr!

    """
    ```julia
    resampled = bsresample(data::AbstractArray, sigma, nrows, [p];
        \t kernel = gaussian,
        \t rng = MersenneTwister(),
        \t return_index = false
    )
    ```
    Bootstrap resample a (sample-per-row / element-per-column) array of `data`
    with uncertainties `sigma` and resampling probabilities `p`
    """
    function bsresample(data::AbstractArray, sigma, nrows::Integer, p=min(0.2,nrows/size(data,1));
            kernel = gaussian,
            rng = MersenneTwister(),
            return_index = false,
        )
        index = Array{Int}(undef, nrows)
        resampled = Array{float(eltype(data))}(undef, nrows, size(data,2))
        bsr!(kernel, resampled, index, data, sigma, p, rng=rng)
        return return_index ? (resampled, index) : resampled
    end
    """
    ```julia
    resampled = bsresample(dataset::Dict, nrows, [elements], [p];
        \t kernel = gaussian,
        \t rng = MersenneTwister()
    )
    ```
    Bootstrap resample a dictionary-based `dataset` with uncertainties stored either
    in `dataset["err"]` or `dataset["[variable]_sigma"]`
    """
    function bsresample(dataset::Dict, nrows::Integer, elements=dataset["elements"], p=min(0.2,nrows/length(dataset[first(elements)]));
            kernel = gaussian,
            rng = MersenneTwister(),
            sigma = :auto,
        )
        # 2d array of nominal values
        data = unelementify(dataset, elements, floatout=true)

        # 2d array of absolute 1-sigma uncertainties
        if sigma === :auto
            if haskey(dataset, "err") && isa(dataset["err"], Dict)
                sigma = unelementify(dataset["err"], elements, floatout=true)
            else
                sigma = unelementify(dataset, elements.*"_sigma", floatout=true)
            end
        end

        # Resample
        sdata = bsresample(data, sigma, nrows, p, kernel=kernel, rng=rng)
        return elementify(sdata, elements, skipstart=0, importas=:Dict)
    end
    function bsresample(dataset::NamedTuple, nrows, elements, p=min(0.2,nrows/length(dataset[first(elements)]));
            kernel = gaussian,
            rng = MersenneTwister(),
            sigma = :auto,
        )
        # 2d array of nominal values
        data = unelementify(dataset, elements, floatout=true)

        # 2d array of absolute 1-sigma uncertainties
        if sigma === :auto
            elements_sigma = (String(e)*"_sigma" for e in elements)
            sigma = unelementify(dataset, elements_sigma, floatout=true)
        end

        # Resample
        sdata = bsresample(data, sigma, nrows, p, kernel=kernel, rng=rng)
        return elementify(sdata, elements, skipstart=0, importas=:Tuple)
    end

    export bsresample


    function randsample!(resampled::AbstractArray, data::AbstractArray, nrows::Integer, p::Number, rng::AbstractRNG=MersenneTwister(), buffer::Vector{Int}=Array{Int}(undef,size(data,1)))
        # Prepare
        ndata = size(data,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            nrows_accepted = 0
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    buffer[nrows_accepted] = i
                end
            end
            nrows_new = min(nrows_accepted, nrows - n)

            # Columns go in outer loop because of column major indexing
            @inbounds @simd for j=1:ncolumns
                # Optimized inner loop
                for i = 1:nrows_new
                    resampled[n+i,j] = data[buffer[i],j]
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function randsample!(resampled::AbstractArray, data::AbstractArray, nrows::Integer, p::AbstractVector, rng::AbstractRNG=MersenneTwister(), buffer::Vector{Int}=Array{Int}(undef,size(data,1)))
        # Prepare
        ndata = size(data,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            nrows_accepted = 0
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    buffer[nrows_accepted] = i
                end
            end
            nrows_new = min(nrows_accepted, nrows - n)

            # Columns go in outer loop because of column major indexing
            @inbounds @simd for j=1:ncolumns
                # Optimized inner loop
                for i = 1:nrows_new
                    resampled[n+i,j] = data[buffer[i],j]
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end

    """
    ```julia
    randsample(data::Array, nrows, [p])
    ```
    Bootstrap resample (without uncertainty) a `data` array to length `nrows`.
    Optionally provide weights `p` either as a vector (one-weight-per-sample) or scalar.
    """
    function randsample(data::AbstractArray, nrows::Integer, p=min(0.2,nrows/size(data,1));
            rng::AbstractRNG=MersenneTwister(),
            buffer::Vector{Int}=Array{Int}(undef,size(data,1))
        )
        resampled = Array{eltype(data)}(undef,nrows,size(data,2))
        return randsample!(resampled, data, nrows, p, rng, buffer)
    end
    # Second method for randsample that takes a dictionary as input
    """
    ```julia
    randsample(dataset::Dict, nrows, [elements], [p])
    ```
    Bootstrap resample (without uncertainty) a `dataset` dict to length `nrows`.
    Optionally provide weights `p` either as a vector (one-weight-per-sample) or scalar.
    """
    function randsample(dataset::Dict, nrows::Integer, elements=dataset["elements"], p=min(0.2,nrows/length(dataset[first(elements)])))
        data = unelementify(dataset, elements, floatout=true)
        sdata = randsample(data, nrows, p)
        return elementify(sdata, elements, skipstart=0, importas=:Dict)
    end
    export randsample


## --- Bin a dataset by a given independent variable

    """
    ```julia
    (bincenters, N) = bincounts(x::AbstractArray, xmin::Number, xmax::Number, nbins::Integer)
    ```
    Tally the number of samples that fall into each of `nbins` equally spaced `x`
    bins between `xmin` and `xmax`, aligned with bin edges as
    `xmin:(xmax-xmin)/nbins:xmax`
    """
    function bincounts(x::AbstractArray, xmin::Number, xmax::Number, nbins::Integer)
        # Tally the number of samples (either resampled or corrected/original) that fall into each bin
        binwidth = (xmax-xmin)/nbins
        bincenters = (xmin+binwidth/2):binwidth:(xmax-binwidth/2)

        # Calculate index from x value
        scalefactor = nbins / (xmax - xmin)
        index = (x .- xmin) .* scalefactor

        # Add up the results
        N = fill(0,nbins)
        @inbounds for i in index
            if 0 < i <= nbins
                N[ceil(Int, i)] += 1
            end
        end
        return bincenters, N
    end
    export bincounts

    """
    ```julia
    (c,m,e) = binmeans(x, y, xmin, xmax, nbins, [weight]; resamplingratio::Number=1)
    ```
    The means (ignoring NaNs) of `y` values binned by `x`, into each of `nbins`
    equally spaced `x` bins between `xmin` and `xmax`, returning bincenters,
    means, and standard errors of the mean.

    To calculate binned medians only (without uncertainties), see `nanmean`

    ### Examples
    ```julia
    (c,m,e) = binmeans(x, y, 0, 4000, 40)
    ```
    """
    function binmeans(x::AbstractArray, y::AbstractArray, xmin::Number, xmax::Number, nbins::Integer; resamplingratio::Number=1)
        binwidth = (xmax-xmin)/nbins
        bincenters = (xmin+binwidth/2):binwidth:(xmax-binwidth/2)

        # Calculate index from x value
        scalefactor = nbins / (xmax - xmin)
        index_float = (x .- xmin) .* scalefactor

        # Calculate the nanmeans and nansems
        N = fill(0,nbins)
        mu = fill(0.0,nbins)
        s2 = fill(0.0,nbins)
        for i ∈ eachindex(x)
            if (0 < index_float[i] <= nbins) && !isnan(y[i])
                index = ceil(Int, index_float[i])
                N[index] += 1
                delta = y[i] - mu[index]
                mu[index] += delta / N[index]
                s2[index] += delta * (y[i] - mu[index])
            end
        end

        # Calculate standard errors
        se = Array{Float64}(undef, nbins)
        for i=1:nbins
            if N[i] > 0
                s2[i] /= N[i] - 1 # Variance
                se[i] = sqrt(s2[i] * resamplingratio / N[i])
            else
                mu[i] = NaN
                se[i] = NaN
            end
        end
        return bincenters, mu, se
    end
    function binmeans(x::AbstractArray, y::AbstractArray, min::Number, max::Number, nbins::Integer, weight::AbstractArray; resamplingratio::Number=1)
        binwidth = (max-min)/nbins
        binedges = range(min,max,length=nbins+1)
        bincenters = (min+binwidth/2):binwidth:(max-binwidth/2)

        means = Array{Float64}(undef,nbins)
        errors = Array{Float64}(undef,nbins)
        for i = 1:nbins
            t = (binedges[i] .< x .<= binedges[i+1]) .& (y.==y)
            w = weight[t]
            yₜ = y[t]
            means[i] = nanmean(yₜ, w)
            errors[i] = nanstd(yₜ, w) * sqrt(resamplingratio / count(t))
        end

        return bincenters, means, errors
    end
    export binmeans

    """
    ```julia
    binmedians(x::AbstractArray, y::AbstractArray, min::Number, max::Number, nbins::Integer;
        \tresamplingratio::Number=1
    )
    ```

    The medians (ignoring NaNs) of `y` values binned by `x`, into each of `nbins`
    equally spaced `x` bins between `xmin` and `xmax`, returning bincenters, medians,
    and equivalent standard errors of the mean (1.4828 * median abolute deviation)

    To calculate binned medians only (without uncertainties), see `nanmedian`

    ### Examples
    ```julia
    (c,m,e) = binmedians(x, y, 0, 4000, 40)
    ```
    """
    function binmedians(x::AbstractArray, y::AbstractArray, min::Number, max::Number, nbins::Integer; resamplingratio::Number=1)
        binwidth = (max-min)/nbins
        binedges = range(min,max,length=nbins+1)
        bincenters = (min+binwidth/2):binwidth:(max-binwidth/2)

        medians = Array{Float64}(undef,nbins)
        errors = Array{Float64}(undef,nbins)
        for i = 1:nbins
            t = (binedges[i] .< x .<= binedges[i+1])
            yₜ = y[t]
            medians[i] = nanmedian!(yₜ)
            errors[i] = 1.4826 * nanmad!(yₜ) * sqrt(resamplingratio / countnotnans(yₜ))
        end
        return bincenters, medians, errors
    end
    export binmedians

## --- Bin bootstrap resampled data

    """
    ```julia
    bin_bsr([f!::Function=nanbinmean!], x::Vector, y::VecOrMat, xmin, xmax, nbins, [w];
        \tx_sigma = zeros(size(x)),
        \ty_sigma = zeros(size(y)),
        \tnresamplings = 1000,
        \tsem = :sigma,
        \tp = 0.2
    )
    ```
    Returns the bincenters `c`, means or medians `m`, and uncertainties of the
    mean or median for a variable `y` binned by independent variable `x` into
    `nbins` equal bins between `xmin` and `xmax`, after `nresamplings` boostrap
    resamplings with acceptance probability `p`.

    If a 2-d array (matrix) of `y` values is provided, each column will be treated
    as a separate variable, means and uncertainties will be returned column-wise.

    Optional keyword arguments and defaults:

        x_sigma = zeros(size(x))

    A vector representing the uncertainty (standard deviation) of each x value

        y_sigma = zeros(size(y))

    A vector representing the uncertainty (standard deviation) of each y value

        nresamplings = 1000

    The number of resamplings to conduct

        sem = :sigma

    Format of the uncertainty estimate of the distribution of the mean.
    If `:sigma` is chosen, a tuple of three vectors `(c, m, e)` will be returned,
    where `e` is the standard error of the mean.
    If `:CI` or `:pctile` is chosen, a tuple of four vectors `(c, m, el, eu)`
    will be returned, where `el` and `eu` are the lower and upper bounds of the 95%
    confidence interval.

        p = 0.2

    Resampling probabilities, either as a scalar or a vector of the same length as `x`


    ### Examples:
    ```julia
    (c,m,e) = bin_bsr(nanbinmedian!, x, y, 0, 4000, 40, x_sigma=0.05x, p=probability, sem=:sigma)
    ```
    ```julia
    (c,m,el,eu) = bin_bsr(nanbinmean!, x, y, 0, 4000, 40, x_sigma=0.05x, p=probability, sem=:pctile)
    ```
    """
    function bin_bsr(f!::Function, x::AbstractVector, y::AbstractVector, xmin, xmax, nbins::Integer;
            x_sigma = zeros(size(x)),
            y_sigma = zeros(size(y)),
            nresamplings = 1000,
            sem = :credibleinterval,
            p = 0.2
        )

        data = hcat(x, y)
        sigma = hcat(x_sigma, y_sigma)
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        means = Array{Float64}(undef, nbins, nresamplings)
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins)
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            f!(view(means,:,i), N, view(dbs,:,1), view(dbs,:,2), xmin, xmax, nbins)
        end

        # Return summary of results
        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        if sem === :sigma
            # Standard deviation of means (sem)
            e = nanstd(means,dim=2)
            return c, m, e
        elseif sem === :credibleinterval || sem === :CI || sem === :pctile
            # Lower bound of central 95% CI of means
            el = m .- nanpctile!(means,2.5,dim=2)
            # Upper bound of central 95% CI of means
            eu = nanpctile!(means,97.5,dim=2) .- m
            return c, m, el, eu
        else
            return c, means
        end
    end
    function bin_bsr(f!::Function, x::AbstractVector, y::AbstractMatrix, xmin, xmax, nbins::Integer;
            x_sigma = zeros(size(x)),
            y_sigma = zeros(size(y)),
            nresamplings = 1000,
            sem = :credibleinterval,
            p = 0.2
        )

        data = hcat(x, y)
        sigma = hcat(x_sigma, y_sigma)
        dtype = float(eltype(data))
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{dtype}(undef, nrows, ncols)
        means = Array{dtype}(undef, nbins, nresamplings, size(y,2))
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins, size(y,2))
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            f!(view(means,:,i,:), N, view(dbs,:,1), view(dbs,:,2:1+size(y,2)), xmin, xmax, nbins)
        end

        # Return summary of results
        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        if sem === :sigma
            m = Array{dtype}(undef, nbins, size(y,2))
            e = Array{dtype}(undef, nbins, size(y,2))
            for j = 1:size(y,2)
                m[:,j] .= nanmean(view(means,:,:,j),dim=2) # Mean-of-means
                e[:,j] .= nanstd(view(means,:,:,j),dim=2) # Standard deviation of means (sem)
            end
            return c, m, e
        elseif sem === :credibleinterval || sem === :CI || sem === :pctile
            m = Array{dtype}(undef, nbins, size(y,2))
            el = Array{dtype}(undef, nbins, size(y,2))
            eu = Array{dtype}(undef, nbins, size(y,2))
            for j = 1:size(y,2)
                m[:,j] .= nanmean(view(means,:,:,j),dim=2) # Mean-of-means
                el[:,j] .= m[:,j] .- nanpctile!(view(means,:,:,j), 2.5, dim=2)
                eu[:,j] .= nanpctile!(view(means,:,:,j), 97.5, dim=2) .- m[:,j]
            end
            return c, m, el, eu
        else
            return c, means
        end
    end
    function bin_bsr(f!::Function, x::AbstractVector, y::AbstractVector, xmin, xmax, nbins::Integer, w::AbstractVector;
            x_sigma = zeros(size(x)),
            y_sigma = zeros(size(x)),
            nresamplings = 1000,
            sem = :credibleinterval,
            p = 0.2
        )

        data = hcat(x, y, w)
        sigma = hcat(x_sigma, y_sigma, zeros(size(w)));

        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        means = Array{Float64}(undef, nbins, nresamplings)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins)
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            f!(view(means,:,i), N, view(dbs,:,1), view(dbs,:,2), view(dbs,:,3), xmin, xmax, nbins)
        end

        # Return summary of results
        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        if sem === :sigma
            # Standard deviation of means (sem)
            e = nanstd(means,dim=2)
            return c, m, e
        elseif sem === :credibleinterval || sem === :CI || sem === :pctile
            # Lower bound of central 95% CI of means
            el = m .- nanpctile!(means,2.5,dim=2)
            # Upper bound of central 95% CI of means
            eu = nanpctile!(means,97.5,dim=2) .- m
            return c, m, el, eu
        else
            return c, means
        end
    end
    bin_bsr(x::AbstractVector, y::AbstractVecOrMat, args...; sem=:sigma, kwargs...) = bin_bsr(nanbinmean!, x, y, args...; sem=sem, kwargs...)
    export bin_bsr

    bin_bsr_means(args...; kwargs...) = bin_bsr(nanbinmean!, args...; kwargs...)
    export bin_bsr_means

    bin_bsr_medians(args...; kwargs...) = bin_bsr(nanbinmedian!, args...; kwargs...)
    export bin_bsr_medians

    """
    ```julia
    (c, m, el, eu) = bin_bsr_ratios([f!::Function=nanbinmean!], x::Vector, num::Vector, denom::Vector, xmin, xmax, nbins, [w];
        \tx_sigma = zeros(size(x)),
        \tnum_sigma = zeros(size(num)),
        \tdenom_sigma = zeros(size(denom)),
        \tnresamplings = 1000,
        \tp::Union{Number,Vector} = 0.2
    )
    ```

    Returns the bincenters `c`, means `m`, as well as upper (`el`) and lower (`eu`) 95% CIs of the mean
    for a ratio `num`/`den` binned by `x` into `nbins` equal bins between `xmin` and `xmax`,
    after `nresamplings` boostrap resamplings with acceptance probability `p`.
    """
    function bin_bsr_ratios(f!::Function, x::AbstractVector, num::AbstractVector, denom::AbstractVector, xmin, xmax, nbins::Integer;
            x_sigma::AbstractVector=zeros(size(x)),
            num_sigma::AbstractVector=zeros(size(num)),
            denom_sigma::AbstractVector=zeros(size(denom)),
            nresamplings=1000,
            p::Union{Number,AbstractVector}=0.2
        )

        data = hcat(x, num, denom)
        sigma = hcat(x_sigma, num_sigma, denom_sigma)
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        means = Array{Float64}(undef, nbins, nresamplings)
        fractions = Array{Float64}(undef, nrows)
        fraction_means = Array{Float64}(undef, nbins)
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins) # Array of bin counts -- Not used but preallocated for speed
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            @views @turbo @. fractions = dbs[:,2] / (dbs[:,2] + dbs[:,3])
            f!(fraction_means, N, view(dbs,:,1), fractions, xmin, xmax, nbins)
            @. means[:,i] = fraction_means / (1 - fraction_means)
        end

        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        el = m .- nanpctile!(means,2.5,dim=2) # Lower bound of central 95% CI
        eu = nanpctile!(means,97.5,dim=2) .- m # Upper bound of central 95% CI

        return c, m, el, eu
    end
    function bin_bsr_ratios(f!::Function, x::AbstractVector, num::AbstractVector, denom::AbstractVector, xmin, xmax, nbins::Integer, w::AbstractVector;
            x_sigma::AbstractVector=zeros(size(x)),
            num_sigma::AbstractVector=zeros(size(num)),
            denom_sigma::AbstractVector=zeros(size(denom)),
            nresamplings=1000,
            p::Union{Number,AbstractVector}=0.2
        )

        data = hcat(x, num, denom, w)
        sigma = hcat(x_sigma, num_sigma, denom_sigma, zeros(size(w)))
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        means = Array{Float64}(undef, nbins, nresamplings)
        fractions = Array{Float64}(undef, nrows)
        fraction_means = Array{Float64}(undef, nbins)
        rng = MersenneTwister()
        W = Array{Float64}(undef, nbins) # Array of bin weights -- Not used but preallocated for speed
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            @views @turbo @. fractions = dbs[:,2] / (dbs[:,2] + dbs[:,3])
            f!(fraction_means, W, view(dbs,:,1), fractions, view(dbs,:,4), xmin, xmax, nbins)
            @. means[:,i] = fraction_means / (1 - fraction_means)
        end

        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        el = m .- nanpctile!(means,2.5,dim=2) # Lower bound of central 95% CI
        eu = nanpctile!(means,97.5,dim=2) .- m # Upper bound of central 95% CI

        return c, m, el, eu
    end
    bin_bsr_ratios(x::AbstractVector, args...; kwargs...) = bin_bsr_ratios(nanbinmean!, x, args...; kwargs...)
    export bin_bsr_ratios

    """
    ```julia
    (c, m, el, eu) = bin_bsr_ratio_medians(x::Vector, num::Vector, denom::Vector, xmin, xmax, nbins, [w];
        \tx_sigma = zeros(size(x)),
        \tnum_sigma = zeros(size(num)),
        \tdenom_sigma = zeros(size(denom)),
        \tnresamplings = 1000,
        \tp::Union{Number,Vector} = 0.2
    )
    ```
    Equivalent to `bin_bsr_ratios(nanbinmedian!, ...)`
    """
    bin_bsr_ratio_medians(args...; kwargs...) = bin_bsr_ratios(nanbinmedian!,args...; kwargs...)
    export bin_bsr_ratio_medians


## --- Quick Monte Carlo binning/interpolation functions

    function mcfit(x::AbstractVector, σx::AbstractVector, y::AbstractVector, σy::AbstractVector,
            xmin::Number, xmax::Number, nbins::Integer=10;
            binwidth::Number=(xmax-xmin)/(nbins-1),
            minrows::Number=100000
        )
        # Run a simplified Monte Carlo fit with nbins of witdth binwidth between xmin and xmax

        # Remove missing data
        hasdata = .!(isnan.(x) .| isnan.(y))
        x′ = x[hasdata]
        y′ = y[hasdata]
        σx′ = σx[hasdata]
        σy′ = σy[hasdata]

        # Fill in variances where not provided explicitly
        σx′[isnan.(σx′)] .= nanstd(x′)
        σy′[isnan.(σy′)] .= nanstd(y′)

        # Increase x uncertainty if x sampling is sparse
        xsorted = [xmin; sort(x′); xmax]
        minerr = maximum(xsorted[2:end] - xsorted[1:end-1]) / 2
        σx′[σx′ .< minerr] .= minerr

        # Bin centers
        c = xmin:(xmax-xmin)/(nbins-1):xmax
        halfwidth = binwidth / 2

        # Run the Monte Carlo
        N = fill(0, nbins)
        m = fill(zero(float(eltype(y′))), nbins)
        xresampled = similar(x′, float(eltype(x′)))
        yresampled = similar(y′, float(eltype(y′)))
        @inbounds for n = 1:ceil(Int, minrows/length(x′))
            randn!(xresampled)
            xresampled .= x′ .+ σx′ .* xresampled
            randn!(yresampled)
            yresampled .= y′ .+ σy′ .* yresampled
            for j = 1:nbins
                l = (c[j] - halfwidth)
                u = (c[j] + halfwidth)
                for i ∈ eachindex(xresampled)
                    if l < xresampled[i] <= u
                        m[j] += yresampled[i]
                        N[j] += 1
                    end
                end
            end
        end
        m ./= N

        return c, m
    end
    export mcfit

## --- Downsample an image / array

    function downsample(A::Matrix, factor::Integer, jfactor::Integer=factor)
        rows = size(A,1) ÷ factor
        cols = size(A,2) ÷ jfactor

        result = similar(A, rows,cols)
        @turbo for i=1:rows
            for j=1:cols
                iₛ = i*factor
                jₛ = j*factor
                result[i,j] = A[iₛ, jₛ]
            end
        end
        return result
    end
    function downsample(A::AbstractMatrix, factor::Integer, jfactor::Integer=factor)
        rows = size(A,1) ÷ factor
        cols = size(A,2) ÷ jfactor

        result = similar(A, rows,cols)
        @inbounds for i=1:rows
            for j=1:cols
                result[i,j] = A[i*factor, j*factor]
            end
        end
        return result
    end
    downsample(A::AbstractArray, factor::Integer) = A[factor:factor:end]
    export downsample

## --- Spatiotemporal sample weighting

    const PI_180 = pi/180

    """
    ```julia
    k = invweight(lat::AbstractArray, lon::AbstractArray, age::AbstractArray;
        \tlp::Number=2,
        \tspatialscale=1.8,
        \tagescale=38.0
    )
    ```

    Find the inverse weights `k` (proportional to spatiotemporal sample density) for
    a set of geological samples with specified latitude (`lat`), logitude (`lon`),
    and `age` (of crystallization, deposition, etc.).

    The default `spatialscale` and `agescale` are taken from Keller and Schoene 2012.
    However, alternative scalings can be supplied. If an array is supplied for either
    `spatialscale`, `agescale`, or both, a 3-d matrix of `k` values will be returned,
    with dimensions length(`spatialscale`)*length(`agescale`)*nrows.
    """
    function invweight(lat::AbstractArray, lon::AbstractArray, age::AbstractArray;
        lp::Number=2, spatialscale=1.8, agescale=38.0)

        # Convert lat and lon to radians
        latᵣ = materialize(vec(lat*PI_180))
        lonᵣ = materialize(vec(lon*PI_180))
        spatialscaleᵣ = materialize(spatialscale*PI_180)
        ageᵥ = materialize(age)

        if any(isnan, latᵣ) || any(isnan, lonᵣ) || any(isnan, ageᵥ)
            k = fill(Inf, length(latᵣ))
            t = @turbo @. !(isnan(latᵣ) | isnan(lonᵣ) | isnan(ageᵥ))
            k[t] .= invweight_nonans(latᵣ[t], lonᵣ[t], ageᵥ[t], lp, spatialscaleᵣ, materialize(agescale))
            return k
        else
            return invweight_nonans(latᵣ, lonᵣ, ageᵥ, lp, spatialscaleᵣ, materialize(agescale))
        end
    end
    function invweight_nonans(latᵣ::AbstractArray, lonᵣ::AbstractArray, age::AbstractArray, lp::Number, spatialscaleᵣ::Number, agescale::Number)

        # Precalculate some sines and cosines
        latsin = @turbo sin.(latᵣ)
        latcos = @turbo cos.(latᵣ)

        # Allocate and fill ks
        N = length(latᵣ)
        k = Array{Float64}(undef, N)
        p = Progress(N÷10, desc="Calculating weights: ")
        @inbounds @batch for i ∈ 1:N
            # Calculate weight
            kᵢ = 0.0
            @turbo for j ∈ 1:N
                Δdᵣ = acos(min( latsin[i] * latsin[j] + latcos[i] * latcos[j] * cos(lonᵣ[i] - lonᵣ[j]), 1.0 ))
                Δa = abs(age[i] - age[j])
                kᵢⱼ = 1.0 / ((Δdᵣ/spatialscaleᵣ)^lp + 1.0) + 1.0 / ((Δa/agescale)^lp + 1.0)
                kᵢ += kᵢⱼ
            end
            k[i] = kᵢ
            (i % 10 == 0) && next!(p)
        end
        return k
    end
    function invweight_nonans(latᵣ::AbstractArray, lonᵣ::AbstractArray, age::AbstractArray, lp::Number, spatialscaleᵣ, agescale)

        # Precalculate some sines and cosines
        latsin = @turbo sin.(latᵣ)
        latcos = @turbo cos.(latᵣ)

        # Allocate and fill ks
        N = length(latᵣ)
        spatialdistᵣ = similar(latᵣ)
        k = Array{Float64}(undef, length(spatialscaleᵣ), length(agescale), N)
        p = Progress(N÷10, desc="Calculating weights: ")
        @inbounds for i ∈ 1:N
            # Calculate weight
            @turbo @. spatialdistᵣ = acos(min( latsin[i] * latsin + latcos[i] * latcos * cos(lonᵣ[i] - lonᵣ), 1.0 ))
            @batch for g ∈ eachindex(spatialscaleᵣ)
                for h ∈ eachindex(agescale)
                    kᵢ = 0.0
                    @inbounds for j ∈ 1:N
                        kᵢⱼ = 1.0 / ((spatialdistᵣ[j]/spatialscaleᵣ[g])^lp + 1.0) + 1.0 / ((abs(age[i] - age[j])/agescale[h])^lp + 1.0)
                        kᵢ += kᵢⱼ
                    end
                    k[g,h,i] = kᵢ
                end
            end
            (i % 10 == 0) && next!(p)
        end
        return k
    end
    export invweight


    """
    ```julia
    k = invweight_location(lat::AbstractArray, lon::AbstractArray;
        \tlp::Number=2,
        \tspatialscale::Number=1.8
    )
    ```

    Find the inverse weights `k` (proportional to spatial sample density) for
    a set of geological samples with specified latitude (`lat`), and logitude (`lon`).
    """
    function invweight_location(lat::AbstractArray, lon::AbstractArray;
        lp::Number=2, spatialscale::Number=1.8)

        # Convert lat and lon to radians
        latᵣ = materialize(vec(lat*PI_180))
        lonᵣ = materialize(vec(lon*PI_180))
        spatialscaleᵣ = materialize(spatialscale*PI_180)

        if any(isnan, latᵣ) || any(isnan, lonᵣ)
            k = fill(Inf, length(latᵣ))
            t = @turbo @. !(isnan(latᵣ) | isnan(lonᵣ))
            k[t] .= invweight_location_nonans(latᵣ[t], lonᵣ[t], lp, spatialscaleᵣ)
            return k
        else
            return invweight_location_nonans(latᵣ, lonᵣ, lp, spatialscaleᵣ)
        end
    end
    # Inner function that can't handle NaNs
    function invweight_location_nonans(latᵣ::AbstractArray, lonᵣ::AbstractArray, lp::Number, spatialscaleᵣ::Number)

        # Precalculate some sines and cosines
        latsin = @turbo sin.(latᵣ)
        latcos = @turbo cos.(latᵣ)

        # Allocate and fill ks
        N = length(latᵣ)
        k = Array{Float64}(undef,N)
        p = Progress(N÷10, desc="Calculating weights: ")
        @inbounds @batch for i ∈ 1:N
            # Otherwise, calculate weight
            kᵢ = 0.0
            @turbo for j ∈ 1:N
                lc = latsin[i] * latsin[j] + latcos[i] * latcos[j] * cos(lonᵣ[i] - lonᵣ[j])
                Δdᵣ = acos(min(lc , 1.0))
                kᵢⱼ = 1.0 / ( (Δdᵣ/spatialscaleᵣ)^lp + 1.0)
                kᵢ += kᵢⱼ
            end
            k[i] = kᵢ
            (i % 10 == 0) && next!(p)
        end
        return k
    end
    export invweight_location


    """
    ```julia
    k = invweight(nums::AbstractArray, scale::Number; lp=2)
    ```

    Find the inverse weights for a single array `nums` for a given `scale`, and
    exponent `lp` (default lp = 2).

    Returns an array k where k[i] is the "inverse weight" for element i of the
    input array.
    """
    function invweight(nums::AbstractArray, scale::Number; lp=2)
        numsₘ = materialize(nums)
        if any(isnan, numsₘ)
            k = fill(Inf, length(numsₘ))
            t = @turbo @. !isnan(numsₘ)
            k[t] .= invweight_nonans(numsₘ[t], scale, lp,)
            return k
        else
            return invweight_nonans(numsₘ, scale, lp,)
        end
    end
    function invweight_nonans(nums::AbstractArray, scale::Number, lp::Number)

        N = length(nums)
        k = Array{Float64}(undef, N)
        p = Progress(N÷10, desc="Calculating weights: ")
        @inbounds @batch for i ∈ 1:N
            if isnan(nums[i])
                # If there is missing data, set k=inf for weight=0
                k[i] = Inf
            else
                # Otherwise, calculate weight
                kᵢ = 0.0
                @turbo for j ∈ 1:N
                    kᵢⱼ = 1.0 / ( (abs(nums[i] - nums[j])/scale)^lp + 1.0)
                    kᵢ += kᵢⱼ
                end
                k[i] = kᵢ
            end
            (i % 10 == 0) && next!(p)
        end
        return k
    end
    export invweight


    """
    ```julia
    k = invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0)
    ```

    Find the inverse weights `k` (proportional to temporal sample density) for
    a set of geological samples with specified `age` (of crystallization, deposition, etc.).
    """
    function invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0)
        return invweight(age, agescale, lp=lp)
    end
    export invweight_age

## --- End of file