Package 'matter'

Description

Resample the given data at specified points. Interpolation can be performed within a tolerance using several interpolation methods.

Usage

approx1(x, y, xout, interp = "linear", n = length(x),
	tol = NA_real_, tol.ref = "abs", extrap = NA_real_)

Arguments

Details

The algorithm is implemented in C and provides several fast interpolation methods. Note that interpolation is limited to using data within the given tolerance. This is also used to specify the width for kernel-based interpolation methods such as interp = "gaussian". The use of a tolerance also means that interpolating within the range of the data but where no data points are within the tolerance window is considered extrapolation. This can be useful when resampling sparse signals with large empty regions, by setting extrap = 0, and setting an appropriate tolerance.

Value

A vector of the same length as xout, giving the resampled data.

Author(s)

Kylie A. Bemis

See Also

asearch, approx approx2

Examples

x <- c(1.11, 2.22, 3.33, 5.0, 5.1)
y <- x^1.11

approx1(x, y, 2.22) # 2.42359
approx1(x, y, 3.0) # NA
approx1(x, y, 3.0, tol=0.2, tol.ref="x") # 3.801133

Description

Resample the (potentially scattered) 2D data to a rectilinear grid at the specified coordinates. Interpolation can be performed within a tolerance using several interpolation methods.

Usage

approx2(x, y, z, xout, yout,
	interp = "linear", nx = length(x), ny = length(y),
	tol = NA_real_, tol.ref = "abs", extrap = NA_real_)

Arguments

Details

See approx1 for details of the 1D implementation. The 2D implementation is mostly the same, except it uses a kd-tree to quickly find neighboring points.

Note that interp = "linear" and interp = "cubic" use a kernel-based approximation. Traditionally, bilinear and bicubic interpolation use 4 and 16 neighboring points, respectively. However, to support scattered data, approx2 will use as many points as are found within the given tolerance, and scale the kernels accordingly. If the input data falls on a regular grid already, then the tolerance should be specified accordingly. Set tol equal to the sampling rate for interp = "linear" and twice the sampling rate for interp = "cubic".

Value

A vector of the same length as xout, giving the resampled data.

Author(s)

Kylie A. Bemis

See Also

expand.grid, asearch, approx, approx1

Examples

x <- matrix(1:25, nrow=5, ncol=5)

approx2(x, nx=10, ny=10, interp="cubic") # upsampling

Description

Combine peaks from multiple signals.

Usage

# Bin a list of peaks
binpeaks(peaklist, domain = NULL, xlist = peaklist,
    tol = NA_real_, tol.ref = "abs", merge = FALSE,
    na.drop = TRUE)

# Merge peaks
mergepeaks(peaks, n = nobs(peaks), x = peaks,
    tol = NA_real_, tol.ref = "abs",
    na.drop = TRUE)

Arguments

Details

binpeaks() is used to bin a list of peaks from multiple signals to a set of common peaks. The peaks (or their corresponding values) are binned to the given domain values and are averaged within each bin. If domain is not given, then the bins are created from the range of the peak locations and the specified tol.

mergepeaks() is used to merge any peaks with gaps smaller than the given tolerance and whose counts (n) do not indicate that they should be considered separate peaks. The merged peaks are averaged together.

Value

A numeric stream_stat vector, giving the average locations of each peak.

Author(s)

Kylie A. Bemis

Examples

x <- c(0, 1, 1, 2, 3, 2, 1, 4, 5, 1, 1, 0)
y <- c(0, 1, 1, 3, 2, 2, 1, 5, 4, 1, 1, 0)

p1 <- findpeaks(x)
p2 <- findpeaks(y)
binpeaks(list(p1, p2), merge=FALSE)

Description

Summarize a vector in the bins at the specified indices.

Usage

binvec(x, lower, upper, stat = "sum", prob = 0.5)

Arguments

Value

An numeric vector of the summarized (binned) values.

Author(s)

Kylie A. Bemis

Examples

set.seed(1)

x <- sort(runif(20))

binvec(x, c(1,6,11,16), c(5,10,15,20))

binvec(x, seq(from=1, to=16, by=5), stat="mean")

Description

Use a binary search to find approximate matches for the elements of its first argument among those in its second. This implementation allows for returning the index of the nearest match if there are no exact matches. It also allows specifying a tolerance for the comparison.

Usage

bsearch(x, table, tol = 0, tol.ref = "abs",
		nomatch = NA_integer_, nearest = FALSE)

reldiff(x, y, ref = "y")

Arguments

Details

The algorithm is implemented in C and currently only works for 'integer', 'numeric', and 'character' vectors. If there are multiple matches, then the first match that is found will be returned, with no guarantees. If a nonzero tolerance is provided, the closest match will be returned.

The "nearest" match for strings when there are no exact matches is decided by the match with the most initial matching characters. Tolerance is ignored for strings and integers. Behavior is undefined and results may be unexpected if values includes NAs.

Value

A vector of the same length as x, giving the indexes of the matches in table.

Author(s)

Kylie A. Bemis

See Also

asearch, match, pmatch, findInterval

Examples

a <- c(1.11, 2.22, 3.33, 5.0, 5.1)

bsearch(2.22, a) # 2
bsearch(3.0, a) # NA
bsearch(3.0, a, nearest=TRUE) # 3
bsearch(3.0, a, tol=0.1, tol.ref="values") # 3

b <- c("hello", "world!")
bsearch("world!", b) # 2
bsearch("worl", b) # NA
bsearch("worl", b, nearest=TRUE) # 2

Description

This is a generic function for applying cryptographic hash functions and calculating checksums for externally-stored R objects.

Usage

checksum(x, ...)

## S4 method for signature 'character'
checksum(x, algo = "sha1", ...)

## S4 method for signature 'matter_'
checksum(x, algo = "sha1", ...)

Arguments

Details

The method for matter objects calculates checksums of each of the files in the object's paths.

Value

A character vector giving the hash or hashes of the object.

Author(s)

Kylie A. Bemis

See Also

digest

Examples

x <- matter(1:10)
y <- matter(1:10)

checksum(x)
checksum(y) # should be the same

Description

Perform equivalents of apply, lapply, and mapply, but over parallelized chunks of data. This is most useful if accessing the data is potentially time-consuming, such as for file-based matter objects. Operating on chunks reduces the number of I/O operations.

Usage

## Operate on elements/rows/columns
chunkApply(X, MARGIN, FUN, ...,
    simplify = FALSE, outpath = NULL,
    verbose = NA, BPPARAM = bpparam())

chunkLapply(X, FUN, ...,
    simplify = FALSE, outpath = NULL,
    verbose = NA, BPPARAM = bpparam())

chunkMapply(FUN, ...,
    simplify = FALSE, outpath = NULL,
    verbose = NA, BPPARAM = bpparam())


## Operate on complete chunks
chunk_rowapply(X, FUN, ...,
    simplify = "c", nchunks = NA, depends = NULL,
    verbose = NA, RNG = FALSE, BPPARAM = bpparam())

chunk_colapply(X, FUN, ...,
    simplify = "c", nchunks = NA, depends = NULL,
    verbose = NA, RNG = FALSE, BPPARAM = bpparam())

chunk_lapply(X, FUN, ...,
    simplify = "c", nchunks = NA, depends = NULL,
    verbose = NA, RNG = FALSE, BPPARAM = bpparam())

chunk_mapply(FUN, ..., MoreArgs = NULL,
    simplify = "c", nchunks = NA, depends = NULL,
    verbose = NA, RNG = FALSE, BPPARAM = bpparam())

Arguments

Details

For chunkApply(), chunkLapply(), and chunkMapply():

For vectors and lists, the vector is broken into some number of chunks according to chunks. The individual elements of the chunk are then passed to FUN.

For matrices, the matrix is chunked along rows or columns, based on the number of chunks. The individual rows or columns of the chunk are then passed to FUN.

In this way, the first argument of FUN is analogous to using the base apply, lapply, and mapply functions.

For chunk_rowapply(), chunk_colapply(), chunk_lapply(), and chunk_mapply():

In this situation, the entire chunk is passed to FUN, and FUN is responsible for knowing how to handle a sub-vector or sub-matrix of the original object. This may be useful if FUN is already a function that could be applied to the whole object such as rowSums or colSums.

When this is the case, it may be useful to provide a custom simplify function.

For convenience to the programmer, several attributes are made available when operating on a chunk.

• "chunkid": The index of the chunk currently being processed by FUN.

• "index": The indices of the elements of the chunk, as elements/rows/columns in the original matrix/vector.

• "depends" (optional): If depends is given, then this is a list of indices within the chunk. The length of the list is equal to the number of elements/rows/columns in the chunk. Each list element either NULL or a vector of indices giving the elements/rows/columns of the chunk that should be processed for that index. The indices that should be processed will be non-NULL, and indices that should be ignored will be NULL.

The depends argument can be used to iterate over dependent elements of a vector, or dependent rows/columns of a matrix. This can be useful if the calculation for a particular row/column/element depends on the values of others.

When depends is provided, multiple rows/columns/elements will be passed to FUN. Each element of the depends list should be a vector giving the indices that should be passed to FUN.

For example, this can be used to implement a rolling apply function.

Value

Typically, a list if simplify=FALSE. Otherwise, the results may be coerced to a vector or array.

Author(s)

Kylie A. Bemis

See Also

apply, lapply, mapply, RNGkind, RNGStreams

Examples

register(SerialParam())

set.seed(1)
x <- matrix(rnorm(1000^2), nrow=1000, ncol=1000)

out <- chunkApply(x, 1L, mean, nchunks=10)
head(out)

Description

The chunked class implements chunked wrappers for vectors, arrays, and lists for parallel iteration.

Usage

## Instance creation
chunked_vec(x, nchunks = NA,
    depends = NULL, drop= FALSE)

chunked_mat(x, margin, nchunks = NA,
    depends = NULL, drop= FALSE)

chunked_list(..., nchunks = NA,
    depends = NULL, drop= FALSE)

## Additional methods documented below

Arguments

Value

An object derived from class chunked.

Slots

data:

The data.

index:

The chunk indices.

drop:

The value passed to drop when subsetting the chunks.

margin:

The array margin for the chunks.

Creating Objects

chunked_vec, chunked_mat and chunked_list instances can be created through chunked_vec(), chunked_mat(), and chunked_list(), respectively.

Methods

Standard generic methods:

length(x):

Get number of chunks.

lengths(x):

Get chunk sizes for each chunk.

x[i, ...]:

Get chunks.

x[[i]]:

Get a single chunk.

Author(s)

Kylie A. Bemis

See Also

chunkApply

Examples

x <- matrix(runif(200), nrow=10, ncol=20)

y <- chunked_mat(x, margin=2, nchunks=5)
print(y)

Description

Apply the equivalent of scale to either columns or rows of a matrix, using a grouping variable.

Usage

## S4 method for signature 'ANY'
colscale(x, center=TRUE, scale=TRUE,
    group = NULL, ..., BPPARAM = bpparam())

## S4 method for signature 'ANY'
rowscale(x, center=TRUE, scale=TRUE,
    group = NULL, ..., BPPARAM = bpparam())

Arguments

Details

See scale for details.

Value

A matrix-like object (usually of the same class as x) with either ‘col-scaled:center’ and ‘col-scaled:scale’ attributes or ‘row-scaled:center’ and ‘row-scaled:scale’ attributes.

Author(s)

Kylie A. Bemis

See Also

scale

Examples

x <- matter(1:100, nrow=10, ncol=10)

colscale(x)

Description

These functions perform calculation of summary statistics over matrix rows and columns for each level of a grouping variable.

Usage

## S4 method for signature 'ANY'
rowStats(x, stat, ..., BPPARAM = bpparam())

## S4 method for signature 'ANY'
colStats(x, stat, ..., BPPARAM = bpparam())

## S4 method for signature 'matter_mat'
rowStats(x, stat, ..., BPPARAM = bpparam())

## S4 method for signature 'matter_mat'
colStats(x, stat, ..., BPPARAM = bpparam())

## S4 method for signature 'sparse_mat'
rowStats(x, stat, ..., BPPARAM = bpparam())

## S4 method for signature 'sparse_mat'
colStats(x, stat, ..., BPPARAM = bpparam())

.rowStats(x, stat, group = NULL,
    na.rm = FALSE, simplify = TRUE, drop = TRUE,
    iter.dim = 1L, BPPARAM = bpparam(), ...)

.colStats(x, stat, group = NULL,
    na.rm = FALSE, simplify = TRUE, drop = TRUE,
    iter.dim = 2L, BPPARAM = bpparam(), ...)

Arguments

Details

The summary statistics methods are calculated over chunks of the matrix using s_colstats and s_rowstats. For matter objects, the iteration is performed over the major dimension for IO efficiency.

Value

A list for each stat requested, where each element is either a vector (if no grouping variable is provided) or a matrix where each column corresponds to a different level of group.

If drop=TRUE, and only a single statistic is requested, then the result will be unlisted and returned as a vector or matrix.

Author(s)

Kylie A. Bemis

See Also

colSums

Examples

register(SerialParam())

set.seed(1)

x <- matrix(runif(100^2), nrow=100, ncol=100)

g <- as.factor(rep(letters[1:5], each=20))

colStats(x, "mean", group=g)

Description

Apply the equivalent of sweep to either columns or rows of a matrix, using a grouping variable.

Usage

## S4 method for signature 'ANY'
colsweep(x, STATS, FUN = "-", group = NULL, ...)

## S4 method for signature 'matter_mat'
colsweep(x, STATS, FUN = "-", group = NULL, ...)

## S4 method for signature 'sparse_mat'
colsweep(x, STATS, FUN = "-", group = NULL, ...)

## S4 method for signature 'ANY'
rowsweep(x, STATS, FUN = "-", group = NULL, ...)

## S4 method for signature 'matter_mat'
rowsweep(x, STATS, FUN = "-", group = NULL, ...)

## S4 method for signature 'sparse_mat'
rowsweep(x, STATS, FUN = "-", group = NULL, ...)

Arguments

Details

See sweep for details.

Value

A matrix-like object (usually of the same class as x) with the statistics swept out.

Author(s)

Kylie A. Bemis

See Also

sweep

Examples

set.seed(1)

x <- matrix(1:100, nrow=10, ncol=10)

colsweep(x, colStats(x, "mean"))

Description

Convolve a signal with weights at arbitrary indices.

Usage

convolve_at(x, index, weights, ...)

Arguments

Details

This is essentially just a weighted sum defined by x[i] = sum(weights[[i]] * x[index[[i]]]).

Value

A numeric vector the same length as x with the smoothed result.

Author(s)

Kylie A. Bemis

Examples

set.seed(1)
t <- seq(from=0, to=6 * pi, length.out=5000)
y <- sin(t) + 0.6 * sin(2.6 * t)
x <- y + runif(length(y))

i <- roll(seq_along(x), width=15)
wt <- dnorm((-7):7, sd=7/2)
wt <- wt / sum(wt)

xs <- convolve_at(x, i, wt)

plot(x, type="l")
lines(xs, col="red")

Description

Compute Manders overlap coefficient (MOC), and Manders colocalization coefficients (M1 and M2), and Dice similarity coefficient.

Usage

coscore(x, y, threshold = NA)

Arguments

Details

The Dice coefficient and Manders overlap coefficient are symmetric between images, while M1 and M2 measure the overlap relative to x and y respectively.

Value

A numeric vector with elements named "MOC", "M1", "M2", and "Dice", and an attribute named "threshold" giving the numeric thresholds (if applicable) for converting each image to a logical mask.

Author(s)

Kylie A. Bemis

References

K. W. Dunn, M. M. Kamocka, and J. H. McDonald. “A practical guide to evaluating colocalization in biological microscop.” American Journal of Physiology: Cell Physiology, vol. 300, no. 4, pp. C732-C742, 2011.

S. V. Costes, D. Daelemans, E. H. Cho, Z. Dobbin, G. Pavlakis, and S. Lockett. “Automatic and Quantitative Measurement of Protein-Protein Colocalization in Live Cells.” Biophysical Journal, vol. 86, no. 6, pp. 3993-4003, 2004.

K. H. Zou, S. K. Warfield, A. Bharatha, C. M. C. Tempany, M. R. Kaus, S. J. Haker, W. M. Wells, III, F. A. Jolesz, and R. Kikinis. “Statistical Validation of Image Segmentation Quality Based on a Spatial Overlap Index.” Academic Radiology, vol. 11, issue 2, pp. 178-189, 2004.

Examples

set.seed(1)
y <- x <- matrix(0, nrow=32, ncol=32)
x[5:16,5:16] <- 1
x[17:28,17:28] <- 1
x <- x + runif(length(x))
y[4:15,4:15] <- 1
y[18:29,18:29] <- 1
y <- y + runif(length(y))
xl <- x > median(x)
yl <- y > median(y)

coscore(x, x)
coscore(x, y)
coscore(x, y, threshold=median)
coscore(xl, yl)

Description

These functions provide simple color palettes.

Usage

## Continuous color palettes
cpal(palette = "Viridis")

## Discrete color palettes
dpal(palette = "Tableau 10")

# Add transparency to colors
add_alpha(colors, alpha = 1, pow = 1)

Arguments

Value

A character vector of colors or a function for generating n colors.

Author(s)

Kylie A. Bemis

See Also

vizi, image

Examples

f <- cpal("viridis")
cols <- f(10)
add_alpha(cols, 1:10/10)

Description

Perform k-fold cross-validation with an arbitrary modeling function.

Usage

cv_do(fit., x, y, folds, ...,
	mi = !is.null(bags), bags = NULL, pos = 1L,
	predict. = predict, transpose = FALSE, keep.models = TRUE,
	trainProcess = NULL, trainArgs = list(),
	testProcess = NULL, testArgs = list(),
	verbose = NA, nchunks = NA, BPPARAM = bpparam())

## S3 method for class 'cv'
fitted(object, type = c("response", "class"),
	simplify = TRUE, ...)

Arguments

Details

The cross-validation is not performed in parallel, because it is assumed the pre-processing functions, modeling function, and prediction function may make use of parallelization. Therefore, these functions need to be able to handle (or ignore) the arguments nchunks and BPPARAM, which will be passed to them.

If bags is specified, then multiple instance learning is assumed, where observations from the same bag are all assumed to have the same label. The labels for bags are automatically pooled (from y) so that if any observation in a bag is pos, then the entire bag is labeled pos. If mi=TRUE then mi_learn will be called by cv_do; otherwise it is assumed fn will handle the multiple instance learning. The accuracy metrics are calculated with the original y labels.

Value

An object of class cv, with the following components:

• average: The average accuracy metrics.

• scores: The fold-specific accuracy metrics.

• folds: The fold memberships.

• fitted.values: The fold-specific predictions.

• models: (Optional) The fitted models.

Author(s)

Kylie A. Bemis

See Also

predscore

Examples

register(SerialParam())

set.seed(1)
n <- 100
p <- 5
nfolds <- 3
y <- rep(c(rep.int("yes", 60), rep.int("no", 40)), nfolds)
x <- matrix(rnorm(nfolds * n * p), nrow=nfolds * n, ncol=p)
x[,1L] <- x[,1L] + 2 * ifelse(y == "yes", runif(n), -runif(n))
x[,2L] <- x[,2L] + 2 * ifelse(y == "no", runif(n), -runif(n))
folds <- rep(paste0("set", seq_len(nfolds)), each=n)

cv_do(pls_nipals, x, y, k=1:5, folds=folds)

Description

Some arithmetic, comparison, and logical operations are available as delayed operations on matter_arr and sparse_arr objects.

Details

Currently the following delayed operations are supported:

‘Arith’: ‘+’, ‘-’, ‘*’, ‘/’, ‘^’, '

‘Math’: ‘exp’, ‘log’, ‘log2’, ‘log10’

Arithmetic operations are applied in C++ layer immediately after the elements are read from virtual memory. This means that operations that are implemented in C and/or C++ for efficiency (such as summary statistics) will also reflect the execution of the deferred arithmetic operations.

Value

A new matter object with the registered deferred operation. Data in storage is not modified; only object metadata is changed.

Author(s)

Kylie A. Bemis

See Also

Arith, Compare, Logic, Ops, Math

Examples

x <- matter(1:100)
y <- 2 * x + 1

x[1:10]
y[1:10]

mean(x)
mean(y)

Description

These functions are provided for compatibility with older versions of matter, and will be removed in the future.

Description

Downsamples a signal for the purposes of visualization. A subset of the original samples are used to represent the signal. The downsampled signal is intended to resemble the original when visualized and should not typically be used for downstream processing.

Usage

downsample(x, n = length(x) / 10L, domain = NULL,
	method = c("lttb", "ltob", "dynamic"))

Arguments

Details

This function implements the downsampling methods from Sveinn Steinarsson's 2013 MSc thesis Downsampling Time Series for Visual Representation, including largest-triangle-three-buckets (LTTB), largest-triangle-one-bucket (LTOB), and dynamic binning.

Value

A vector of length n, giving the downsampled signal.

Author(s)

Kylie A. Bemis

References

S. Steinarsson. “Downsampling Time Series for Visual Representation.” MSc, School of Engineering and Natural Sciences, University of Iceland, Reykjavik, Iceland, June 2013.

See Also

approx1

Examples

set.seed(1)
t <- seq(from=0, to=6 * pi, length.out=2000)
x <- sin(t) + 0.6 * sin(2.6 * t)
x <- x + runif(length(x))
xs <- downsample(x, n=200)
s <- attr(xs, "sample")

plot(x, type="l")
points(s, xs, col="red", type="b")

Description

The drle class stores delta-run-length-encoded vectors. These differ from other run-length-encoded vectors provided by other packages in that they allow for runs of values that each differ by a common difference (delta).

Usage

## Instance creation
drle(x, type = "drle", cr_threshold = 0)

is.drle(x)
## Additional methods documented below

Arguments

Value

An object of class drle.

Slots

values:

The values that begin each run.

lengths:

The length of each run.

deltas:

The difference between the values of each run.

Creating Objects

drle instances can be created through drle().

Methods

Standard generic methods:

x[i]:

Get the elements of the uncompressed vector.

length(x):

Get the length of the uncompressed vector.

c(x, ...):

Combine vectors.

Author(s)

Kylie A. Bemis

See Also

rle

Examples

## Create a drle vector
x <- c(1,1,1,1,1,6,7,8,9,10,21,32,33,34,15)
y <- drle(x)

# Check that their elements are equal
x == y[]

Description

Enhance the contrast in a 2D signal.

Usage

# Adjust by extreme values
enhance_adj(x, frac = 0.01)

# Histogram equalization
enhance_hist(x, nbins = 256L)

# Contrast-limited adaptive histogram equalization (CLAHE)
enhance_adapt(x, width = sqrt(nrow(x) * ncol(x)) %/% 5L,
    clip = 0.1, nbins = 256L)

Arguments

Details

enhance_adj() performs a simple adjustment of the overall image intensity range by clamping a fraction of the highest and lowest pixel values. This is useful for suppressing very bright hotspots, but may not be sufficient for images with globally poor contrast.

enhance_hist() performs histogram equalization. Histogram equalization transforms the pixel values so that the histogram of the image is approximately flat. This is done by replacing the original pixel values with their associated probability in the image's empirical cumulative distribution.

enhance_adapt() performs contrast-limited adaptive histogram equalization (CLAHE) from Zuiderveld (1994). While ordinary histogram equalization performs a global transformation on the image, adaptive histogram equalization calculates a histogram in a local neighborhood around each pixel to perform the transformation, thereby enhancing the local contrast across the image. However, this can amplify local noise, so to avoid this, the histogram is clipped to a maximum allowed bin value before transforming the pixel values. To speed up the computation, it is implemented here using a sliding window technique as described by Wang and Tao (2006).

These methods rescale the output image so that its median equals the median of the original image and it has equal interquartile range (IQR).

Value

A numeric matrix the same dimensions as x with the smoothed result.

Author(s)

Kylie A. Bemis

References

K. Zuiderveld. “Contrast Limited Adaptive Histogram Equalization.” Graphics Gems IV, Academic Press, pp. 474-485, 1994.

Z. Wang and J. Tao. “A Fast Implementation of Adaptive Histogram Equalization.” IEEE 8th international Conference on Signal Processing, Nov 2006.

Examples

set.seed(1)
x <- matrix(0, nrow=32, ncol=32)
x[9:24,9:24] <- 10
x <- x + runif(length(x))
y <- x + rlnorm(length(x))
z <- enhance_hist(y)

par(mfcol=c(1,3))
image(x, col=hcl.colors(256), main="original")
image(y, col=hcl.colors(256), main="multiplicative noise")
image(z, col=hcl.colors(256), main="histogram equalization")

Description

Estimate the continuum (baseline) of a signal.

Usage

# Continuum based on local extrema
estbase_loc(x,
    smooth = c("none", "loess", "spline"),
    span = 1/10, spar = NULL, upper = FALSE)

# Convex hull
estbase_hull(x, upper = FALSE)

# Sensitive nonlinear iterative peak clipping (SNIP)
estbase_snip(x, width = 100L, decreasing = TRUE)

# Running medians
estbase_med(x, width = 100L)

Arguments

Details

estbase_loc() uses a simple method based on linearly interpolating from local extrema. It typically performs well enough for most situations. Signals with strong noise or wide peaks may require stronger smoothing after the interpolation step.

estbase_hull() estimates the continuum by finding the lower or upper convex hull using the monotonic chain algorithm of A. M. Andrew (1979).

estbase_snip() performs sensitive nonlinear iterative peak (SNIP) clipping using the adaptive clipping window from M. Morhac (2009).

estbase_med() estimates the continuum from running medians.

Value

A numeric vector the same length as x with the estimated continuum.

Author(s)

Kylie A. Bemis

References

A. M. Andrew. “Another efficient algorithm for convex hulls in two dimensions.” Information Processing Letters, vol. 9, issue 5, pp. 216-219, Dec. 1979.

M. Morhac. “An algorithm for determination of peak regions and baseline elimination in spectroscopic data.” Nuclear Instruments and Methods in Physics Research A, vol. 600, issue 2, pp. 478-487, Mar. 2009.

Examples

set.seed(1)
t <- seq(from=0, to=6 * pi, length.out=2000)
x <- sin(t) + 0.6 * sin(2.6 * t)
lo <- estbase_hull(x)
hi <- estbase_hull(x, upper=TRUE)

plot(x, type="l")
lines(lo, col="red")
lines(hi, col="blue")

Description

Estimate the raster dimensions of a scattered 2D signal based on its pixel coordinates.

Usage

estdim(x, tol = 1e-6)

Arguments

Value

A numeric vector giving the estimated raster dimensions.

Author(s)

Kylie A. Bemis

Examples

co <- expand.grid(x=1:12, y=1:9)
co$x <- jitter(co$x)
co$y <- jitter(co$y)

estdim(co)

Description

Estimate the noise across a signal.

Usage

# Difference-based noise estimation
estnoise_diff(x, nbins = 1L, overlap = 0.5,
    index = NULL)

# Dynamic noise level filtering
estnoise_filt(x, nbins = 1L, overlap = 0.5,
    msnr = 2, threshold = 0.5, isPeaks = FALSE, index = NULL)

# SD-based noise estimation
estnoise_sd(x, nbins = 1, overlap = 0.5,
    k = 9, index = NULL)

# MAD-based noise estimation
estnoise_mad(x, nbins = 1, overlap = 0.5,
    k = 9, index = NULL)

# Quantile-based noise estimation
estnoise_quant(x, nbins = 1, overlap = 0.5,
    k = 9, prob = 0.95, index = NULL)

Arguments

Details

estnoise_diff() estimates the local noise from the mean absolute differences between the signal and the average of its derivative. For noisy signals, the derivative is dominated by the noise, making it a useful estimator of the noise.

estnoise_sd(), estnoise_mad(), and estnoise_quant() all estimate the local noise by first smoothing the signal with a Gaussian filter, and then subtracting the raw signal from the smoothed signal to isolate the noise component. The noise is then summarized with the corresponding statistic.

estnoise_filt() uses the dynamic noise level filtering algorithm of Xu and Freitas (2010) based on the local signal in an approach similar to Gallia et al. (2013). The peaks in the signal are sorted, and the smallest peak is assumed to be noise and is used to estimate the noise level. Each peak is then compared to the previous peak. A peak is labeled a signal peak only if it exceeds a minimum signal-to-noise ratio. Otherwise, the peak is labeled noise, and the noise level is re-estimated. This process continues until a signal peak is found, and the noise level is estimated from the noise peaks.

Value

A numeric vector the same length as x with the estimated local noise level.

Author(s)

Kylie A. Bemis

References

H. Xu and M. A. Freitas. “A Dynamic Noise Level Algorithm for Spectral Screening of Peptide MS/MS Spectra.” BMC Bioinformatics, vol. 11, no. 436, Aug. 2010.

J. Gallia, K. Lavrich, A. Tan-Wilson, and P. H. Madden. “Filtering of MS/MS data for peptide identification.” BMC Genomics, vol. 14, suppl. 7, Nov. 2013.

Examples

# simple signal
set.seed(1)
n <- 500
x <- rnorm(n)
x <- x + 90 * dnorm(seq_along(x), mean=n/4)
x <- x + 80 * dnorm(seq_along(x), mean=n/2)
x <- x + 70 * dnorm(seq_along(x), mean=3*n/4)

ns <- estnoise_quant(x)
plot(x, type="l")
lines(ns, col="blue")

# simulated spectrum
set.seed(1)
x <- simspec(size=5000)

ns <- estnoise_quant(x, nbins=20)
plot(x, type="l")
lines(ns, col="blue")

Description

Estimate the resolution (approximate sampling rate) of a signal based on its domain values.

Usage

estres(x, tol = NA, ref = NA_character_)

Arguments

Value

A single number named "absolute" or "relative" giving the approximate constant sampling rate matching the given domain values. NA if a sampling rate could not be determined.

Author(s)

Kylie A. Bemis

Examples

x <- seq_rel(501, 600, by=1e-3)

estres(x)

Description

The FastMap algorithm performs approximate multidimensional scaling (MDS) based on any distance function. It is faster and more efficient than traditional MDS algorithms, scaling as O(n) rather than O(n^2). FastMap accomplishes this by finding two distant pivot objects on some hyperplane for each projected dimension, and then projecting all other objects onto the line between these pivots.

Usage

# FastMap projection
fastmap(x, k = 3L, distfun = NULL,
	transpose = FALSE, niter = 3L, verbose = NA, ...)

## S3 method for class 'fastmap'
predict(object, newdata, ...)

# Distance functionals
rowDistFun(x, y, metric = "euclidean", p = 2, weights = NULL,
	verbose = NA, nchunks = NA, BPPARAM = bpparam(), ...)
colDistFun(x, y, metric = "euclidean", p = 2, weights = NULL,
	verbose = NA, nchunks = NA, BPPARAM = bpparam(), ...)

Arguments

Details

The pivots are initialized randomly for each new dimension, so the selection of pivots (and therefore the resulting projection) can be sensitive to the random seed for some datasets.

A custom distance function can be passed via distfun. If not provided, then this defaults to rowDistFun() if transpose=FALSE or colDistFun() if transpose=TRUE.

If a custom function is passed, it should take the form function(x, y, ...), and it must return a function of the form function(i). The returned function should return the distances between the ith object(s) in x and all objects in y. rowDistFun() and colDistFun() are examples of functions that satisfy these properties.

Value

An object of class fastmap, with the following components:

• x: The projected variable matrix.

• sdev: The standard deviations of each column of the projected matrix x.

• pivots: A matrix giving the indices of the pivots and the distances between them.

• pivot.array: A subset of the original data matrix containing only the pivots.

• distfun: The function used to generate the distance function.

Author(s)

Kylie A. Bemis

References

C. Faloutsos, and D. Lin. “FastMap: A Fast Algorithm for Indexing, Data-Mining and Visualization of Traditional and Multimedia Datasets.” Proceedings of the 1995 ACM SIGMOD international conference on Management of data, pp. 163 - 174, June 1995.

See Also

cmdscale, prcomp

Examples

register(SerialParam())
set.seed(1)

a <- matrix(sort(runif(500)), nrow=50, ncol=10)
b <- matrix(rev(sort(runif(500))), nrow=50, ncol=10)
x <- cbind(a, b)

fm <- fastmap(x, k=2)

Description

Smooth a uniformly sampled 1D signal.

Usage

# Moving average filter
filt1_ma(x, width = 5L)

# Linear convolution filter
filt1_conv(x, weights)

# Gaussian filter
filt1_gauss(x, width = 5L, sd = (width %/% 2) / 2)

# Bilateral filter
filt1_bi(x, width = 5L, sddist = (width %/% 2) / 2,
    sdrange = mad(x, na.rm = TRUE))

# Bilateral filter with adaptive parameters
filt1_adapt(x, width = 5L, spar = 1)

# Nonlinear diffusion
filt1_diff(x, niter = 3L, kappa = 50,
    rate = 0.25, method = 1L)

# Guided filter
filt1_guide(x, width = 5L, guide = x,
    sdreg = mad(x, na.rm = TRUE))

# Peak-aware guided filter
filt1_pag(x, width = 5L, guide = NULL,
    sdreg = mad(x, na.rm = TRUE), ftol = 1/10)

# Savitzky-Golay filter
filt1_sg(x, width = 5L, order = min(3L, width - 2L),
    deriv = 0, delta = 1)

Arguments

Details

filt1_ma() performs mean filtering in O(n) time. This is fast and especially useful for calculating other filters that can be constructed as a combination of mean filters.

filt1_gauss() performs Gaussian filtering.

filt1_bi() and filt1_adapt() perform edge-preserving bilateral filtering. The latter calculates the kernel parameters adapatively based on the local signal, using a strategy adapted from Joseph and Periyasamy (2018).

filt1_diff() performs the nonlinear diffusion filtering of Perona and Malik (1990). Rather than relying on a filter width, it progressively diffuses (smooths) the signal over multiple iterations. More iterations will result in a smoother image.

filt1_guide() performs edge-preserving guided filtering. Guided filtering uses a local linear model based on the structure of a so-called "guidance signal". By default, the guidance signal is often the same as the input signal. Guided filtering performs similarly to bilateral filtering, but is often faster (though with more memory use), as it is implemented as a combination of mean filters.

filt1_pag() performs peak-aware guided filtering using a regularization parameter that focuses on preserving peaks rather than edges, using a strategy adapted from Liu and He (2022). By default, the guidance signal is generated by smoothing the input signal with nonlinear diffusion.

filt1_sg() performs traditional Savitzky-Golay filtering, which uses a local least-squares polynomial approximation to perform the smoothing. It reduces noise while attempting to retain the peak shape and height.

Value

A numeric vector the same length as x with the smoothed result.

Author(s)

Kylie A. Bemis

References

J. Joseph and R. Perisamy. “An image driven bilateral filter with adaptive range and spatial parameters for denoising Magnetic Resonance Images.” Computers and Electrical Engineering, vol. 69, pp. 782-795, July 2018.

P. Perona and J. Malik. “Scale-space and edge detection using anisotropic diffusion.” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, issue 7, pp. 629-639, July 1990.

K. He, J. Sun, and X. Tang. “Guided Image Filtering.” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 35, no. 6, pp. 1397-1409, June 2013.

D. Liu and C. He. “Peak-aware guided filtering for spectrum signal denoising.” Chemometrics and Intelligent Laboratory Systems, vol. 222, March 2022.

A. Savitzky and M. J. E. Golay. “Smoothing and Differentiation of Data by Simplified Least Squares Procedures.” Analytical Chemistry, vol. 36, no. 8, pp. 1627-1639, July 1964.

Examples

set.seed(1)
t <- seq(from=0, to=6 * pi, length.out=5000)
y <- sin(t) + 0.6 * sin(2.6 * t)
x <- y + runif(length(y))
xs <- filt1_gauss(x, width=25)

plot(x, type="l")
lines(xs, col="red")

Description

Smooth a uniformly sampled 2D signal.

Usage

# Moving average filter
filt2_ma(x, width = 5L)

# Linear convolution filter
filt2_conv(x, weights)

# Gaussian filter
filt2_gauss(x, width = 5L, sd = (width %/% 2) / 2)

# Bilateral filter
filt2_bi(x, width = 5L, sddist = (width %/% 2) / 2,
    sdrange = mad(x, na.rm = TRUE))

# Bilateral filter with adaptive parameters
filt2_adapt(x, width = 5L, spar = 1)

# Nonlinear diffusion
filt2_diff(x, niter = 3L, kappa = 50,
    rate = 0.25, method = 1L)

# Guided filter
filt2_guide(x, width = 5L, guide = x,
    sdreg = mad(x, na.rm = TRUE))

Arguments

Details

filt2_ma() performs mean filtering in O(n) time. This is fast and especially useful for calculating other filters that can be constructed as a combination of mean filters.

filt2_gauss() performs Gaussian filtering.

filt2_bi() and filt2_adapt() perform edge-preserving bilateral filtering. The latter calculates the kernel parameters adapatively based on the local signal, using a strategy adapted from Joseph and Periyasamy (2018).

filt2_diff() performs the nonlinear diffusion filtering of Perona and Malik (1990). Rather than relying on a filter width, it progressively diffuses (smooths) the signal over multiple iterations. More iterations will result in a smoother image.

filt2_guide() performs edge-preserving guided filtering. Guided filtering uses a local linear model based on the structure of a so-called "guidance signal". By default, the guidance signal is often the same as the input signal. Guided filtering performs similarly to bilateral filtering, but is often faster (though with more memory use), as it is implemented as a combination of mean filters.

Value

A numeric matrix the same dimensions as x with the smoothed result.

Author(s)

Kylie A. Bemis

References

J. Joseph and R. Perisamy. “An image driven bilateral filter with adaptive range and spatial parameters for denoising Magnetic Resonance Images.” Computers and Electrical Engineering, vol. 69, pp. 782-795, July 2018.

P. Perona and J. Malik. “Scale-space and edge detection using anisotropic diffusion.” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, issue 7, pp. 629-639, July 1990.

K. He, J. Sun, and X. Tang. “Guided Image Filtering.” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 35, no. 6, pp. 1397-1409, June 2013.

Examples

set.seed(1)
i <- seq(-4, 4, length.out=12)
j <- seq(1, 3, length.out=9)
co <- expand.grid(i=i, j=j)
x <- matrix(atan(co$i / co$j), nrow=12, ncol=9)
x <- 10 * (x - min(x)) / diff(range(x))
x <- x + 2.5 * runif(length(x))
xs <- filt2_gauss(x)

par(mfcol=c(1,2))
image(x, col=hcl.colors(256), main="original")
image(xs, col=hcl.colors(256), main="smoothed")

Description

Smooth a nonuniformly sampled signal in K-nearest neighbors.

Usage

# Moving average filter
filtn_ma(x, index, k = 5L, metric = "euclidean", p = 2)

# Linear convolution filter
filtn_conv(x, index, weights, metric = "euclidean", p = 2)

# Gaussian filter
filtn_gauss(x, index, k = 5L, sd = median(r) / 2,
    metric = "euclidean", p = 2)

# Bilateral filter
filtn_bi(x, index, k = 5L, sddist = median(r) / 2,
    sdrange = mad(x, na.rm=TRUE), metric = "euclidean", p = 2)

# Bilateral filter with adaptive parameters
filtn_adapt(x, index, k = 5L, spar = 1,
    metric = "euclidean", p = 2)

Arguments

Details

These functions must first perform a K-nearest neighbors search on the sample locations (index) before smoothing the signal (x), but they can be applied to nonuniformly sampled signals in an arbitrary number of dimensions. The nearest neighbor search is performed using a kd-tree. It is efficient for low dimensional signals, but performance will degrade in high dimensions. In general, the complexity of these filters is O(k * n log(n)).

Value

A numeric vector the same dimensions as x with the smoothed result.

Author(s)

Kylie A. Bemis

References

J. Joseph and R. Perisamy. “An image driven bilateral filter with adaptive range and spatial parameters for denoising Magnetic Resonance Images.” Computers and Electrical Engineering, vol. 69, pp. 782-795, July 2018.

Examples

set.seed(1)

# signal intensities
x <- rlnorm(100)

# 3D sampling locations
index <- replicate(3, runif(100))

filtn_ma(x, index, k=5)

Description

Find peaks in a signal based on its local maxima, as determined by a sliding window.

Usage

# Find peaks
findpeaks(x, width = 5L, prominence = NULL,
    snr = NULL, noise = "quant", bounds = TRUE,
    relheight = 0.005, ...)

# Local maxima
locmax(x, width = 5L)

# Local minima
locmin(x, width = 5L)

Arguments

Details

For locmax() and locmin(), a local extremum is defined as an element greater (or less) than all of the elements within width / 2 elements to the left of it, and greater (or less) than or equal to all of the elements within width / 2 elements to the right of it. That is, ties are resolved such that only the leftmost sample is considered the local extremum.

For findpeaks(), the peaks are simply the local maxima of the signal. The peak boundaries are found by descending a local maximum until a local minimum is found on either side, using the same criteria as above. The peaks are optionally filtered based on their prominences.

Optionally, the signal-to-noise ratio (SNR) can be estimated and used for filtering the peaks. These use the functions estnoise_quant, estnoise_diff, estnoise_filt, etc., to estimate the noise in the signal.

Value

For locmax() and locmin(), an logical vector indicating whether each element is a local extremum.

For findpeaks(), an integer vector giving the indices of the peaks, with attributes 'left_bounds' and 'right_bounds' giving the left and right boundaries of the peak as determined using the rule above.

Author(s)

Kylie A. Bemis

See Also

findpeaks_cwt, findpeaks_knn, estnoise_quant, estnoise_sd, estnoise_mad, estnoise_diff, estnoise_filt, peakwidths, peakareas, peakheights, binpeaks, mergepeaks

Examples

# simple signal
x <- c(0, 1, 1, 2, 3, 2, 1, 4, 5, 1, 1, 0)
locmax(x)
findpeaks(x)

# simulated spectrum
set.seed(1)
x <- simspec(size=5000)

# find peaks with snr >= 3
p <- findpeaks(x, snr=3, noise="quant")
plot(x, type="l")
points(p, x[p], col="red")

# find peaks with derivative-based noise
p <- findpeaks(x, snr=3, noise="diff")
plot(x, type="l")
points(p, x[p], col="red")

Description

Find peaks in a signal using continuous wavelet transform (CWT).

Usage

# Find peaks with CWT
findpeaks_cwt(x, snr = 2, wavelet = ricker, scales = NULL,
    maxdists = scales, ngaps = 3L, ridgelen = length(scales) %/% 4L,
    qnoise = 0.95, width = length(x) %/% 20L, bounds = TRUE)

# Find ridges lines in a matrix
findridges(x, maxdists, ngaps)

# Continuous Wavelet Transform
cwt(x, wavelet = ricker, scales = NULL)

Arguments

Details

findpeaks_cwt() uses the peak detection method based on continuous wavelet transform (CWT) proposed by Du, Kibbe, and Lin (2006).

The raw signal is convolved with a wavelet (by default, a Ricker wavelet is used) at a range of different scales. This produces a matrix of CWT coefficients with a number of rows equal to the length of the original signal and each column representing a different scale of convolution.

The convolution at the smallest scales represent a good estimate of noise and peak location. The larger scales represent a smoother signal where larger peaks are prominent and smaller peaks are removed.

The method proceeds by identifying ridge lines in the CWT coefficient matrix using findridges(). Local maxima are identified at each scale and connected across each scale, forming the ridge lines.

Finally, the local noise is estimated from the CWT coefficients at the smallest scale. The peaks are filtered based on signal-to-noise ratio and the length of their ridge lines.

Value

For findpeaks_cwt(), an integer vector giving the indices of the peaks, with attributes 'left_bounds' and 'right_bounds' giving the left and right boundaries of the peak as determined using the rule above.

For findridges(), a list of matrices giving the row and column indices of the entries of each detected ridge line.

Author(s)

Kylie A. Bemis

See Also

findpeaks, peakwidths, peakareas, peakheights, binpeaks, mergepeaks

Examples

# simple signal
x <- c(0, 1, 1, 2, 3, 2, 1, 4, 5, 1, 1, 0)
locmax(x)
findpeaks(x)

# simulated spectrum
set.seed(1)
x <- simspec(size=5000)

# find peaks with snr >= 3
p <- findpeaks_cwt(x, snr=3)
plot(x, type="l")
points(p, x[p], col="red")

# plot ridges
ridges <- attr(p, "ridges")
plot(c(0, length(x)), c(0, 25), type="n")
for ( ri in ridges )
    lines(ri, type="o", pch=20, cex=0.5)

Description

Find peaks in a signal based on its local maxima, as determined by K-nearest neighbors.

Usage

# Find peaks with KNN
findpeaks_knn(x, index, k = 5L,
    snr = NULL, noise = "quant", relheight = 0.005,
    arr.ind = is.array(x), useNames = TRUE, ...)

# Local maxima
knnmax(x, index, k = 5L)

# Local minima
knnmin(x, index, k = 5L)

Arguments

Details

For knnmax() and knnmin(), a local extremum is defined as an element greater (or less) than all of its nearest neighbors with a lesser index, and greater (or less) than or equal to all of its nearest neighbors with a greater or equal index. That is, ties are resolved such that the sample with the lowest index is considered the local extremum.

For findpeaks_knn(), the peaks are simply the local maxima of the signal. The peaks are optionally filtered based on their relative heights.

Optionally, the signal-to-noise ratio (SNR) can be estimated and used for filtering the peaks. These use the functions estnoise_quant, estnoise_diff, estnoise_filt, etc., to estimate the noise in the signal.

Value

For knnmax() and knnmin(), an logical vector or array indicating whether each element is a local extremum.

For findpeaks_knn(), an integer vector (or matrix if arr.ind=TRUE) giving the indices of the peaks.

Author(s)

Kylie A. Bemis

See Also

findpeaks, estnoise_quant, estnoise_sd, estnoise_mad, estnoise_diff, estnoise_filt

Examples

# simple 2D signal
x <- c(
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0,
    0, 0, 0, 2, 0, 0, 1, 4, 2, 0, 1, 1, 0, 0, 0,
    0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 3, 2, 1, 0, 0,
    0, 0, 0, 0, 1, 3, 3, 0, 0, 1, 4, 4, 3, 1, 0,
    0, 0, 0, 0, 0, 3, 2, 0, 1, 0, 3, 2, 3, 0, 0,
    0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 2, 2, 3, 0, 0)
x <- matrix(x, nrow=7, ncol=15, byrow=TRUE)

# find peaks in 2D using KNN
findpeaks_knn(x)

Description

Check if a series of x-y points are contained in a closed 2D polygon.

Usage

inpoly(points, poly)

Arguments

Details

This function works by extending a horizontal ray from each point and counting the number of times it crosses an edge of the polygon.

Value

A logical vector that is TRUE for points that are fully inside the polygon, a vertex, or on an edge, and FALSE for points fully outside the polygon.

Note

There are various public implementations of this function with no clear original source. The version implemented here is loosely based on code by W. Randolph Franklin with modifications so that vertices and points on edges are considered inside the polygon.

Author(s)

W. R. Franklin and Kylie A. Bemis

References

W. R. Franklin. “PNPOLY - Point Inclusion in Polygon Test.” https://wrfranklin.org/Research/Short_Notes/pnpoly.html, 1970.

M. Shimrat, "Algorithm 112, Position of Point Relative to Polygon", Comm. ACM 5(8), pp. 434, Aug 1962.

E. Haines. “Point in Polygon Strategies.” http://www.acm.org/pubs/tog/editors/erich/ptinpoly/, 1994.

See Also

kdsearch

Examples

poly <- data.frame(
        x=c(3,5,5,3),
        y=c(3,3,5,5))
xy <- data.frame(
    x=c(4,6,4,2,3,5,5,3,4,5,4,4,
        3.5,4.5,4.0,3.5),
    y=c(2,4,6,4,3,3,5,5,3,4,5,3,
        4.0,4.0,4.5,4.0),
    ref=c(
        rep("out", 4),
        rep("vertex", 4),
        rep("edge", 4),
        rep("in", 4)))

xy$test <- inpoly(xy[,1:2], poly)
xy

Description

Search a matrix of K-dimensional data points and return the indices of the nearest neighbors or of all data points that are within a specified tolerance in each dimension.

Usage

# Nearest neighbor search
knnsearch(x, data, k = 1L, metric = "euclidean", p = 2)

# Range search
kdsearch(x, data, tol = 0, tol.ref = "abs")

# K-D tree
kdtree(data)

Arguments

Details

knnsearch() performs k-nearest neighbor searches. kdsearch() performs range searches for points within a given tolerance of the query points.

The algorithms are implemented in C and work by building a kd-tree to perform the search. If multiple calls to kdsearch() or knnsearch() are expected on the same data, it can be much faster to build the tree once with kdtree().

A kd-tree is essentially a multidimensional generalization of a binary search tree. Building the search tree is O(n * log n) and searching for a single data point is O(log n).

For knnsearch(), ties are broken based on the original ordering of the rows in data.

Value

For knnsearch(), a matrix with rows equal to the number of rows of x and columns equal to k giving the indices of the k-nearest neighbors.

For kdsearch(), a list with length equal to the number of rows of x, where each list element is a vector of indexes of the matches in data.

Author(s)

Kylie A. Bemis

See Also

bsearch, approx2,

Examples

d <- expand.grid(x=1:10, y=1:10)
x <- rbind(c(1.11, 2.22), c(3.33, 4.44))

knnsearch(x, d, k=3)

Description

The matter_arr class implements out-of-memory arrays.

Usage

## Instance creation
matter_arr(data, type = "double", path = NULL,
    dim = NA_integer_, dimnames = NULL, offset = 0, extent = NA_real_,
    readonly = NA, append = FALSE, rowMaj = FALSE, ...)

matter_mat(data, type = "double", path = NULL,
    nrow = NA_integer_, ncol = NA_integer_, dimnames = NULL,
    offset = 0, extent = NA_real_, readonly = NA,
    append = FALSE, rowMaj = FALSE, ...)

matter_vec(data, type = "double", path = NULL,
    length = NA_integer_, names = NULL, offset = 0, extent = NA_real_,
    readonly = NA, append = FALSE, rowMaj = FALSE, ...)

## Additional methods documented below

Arguments

Value

An object of class matter_arr.

Slots

data:

This slot stores any information necessary to access the data for the object (which may include the data itself and/or paths to file locations, etc.).

type:

The storage mode of the accessed data when read into R. This is a 'factor' with levels 'raw', 'logical', 'integer', 'numeric', or 'character'.

dim:

Either 'NULL' for vectors, or an integer vector of length one of more giving the maximal indices in each dimension for matrices and arrays.

names:

The names of the data elements for vectors.

dimnames:

Either 'NULL' or the names for the dimensions. If not 'NULL', then this should be a list of character vectors of the length given by 'dim' for each dimension. This is always 'NULL' for vectors.

ops:

Deferred arithmetic operations.

transpose:

Indicates whether the data is stored in row-major order (TRUE) or column-major order (FALSE). For a matrix, switching the order that the data is read is equivalent to transposing the matrix (without changing any data).

indexed:

For matter_mat only. Indicates whether the pointers to rows or columns are indexed for quick access or not.

Extends

matter

Creating Objects

matter_arr instances can be created through matter_arr() or matter(). Matrices and vectors can also be created through matter_mat() and matter_vec().

Methods

Standard generic methods:

length(x), length(x) <- value:

Get or set length.

dim(x), dim(x) <- value:

Get or set 'dim'.

names(x), names(x) <- value:

Get or set 'names'.

dimnames(x), dimnames(x) <- value:

Get or set 'dimnames'.

x[...], x[...] <- value:

Get or set the elements of the array.

cbind(x, ...), rbind(x, ...):

Combine matrices by row or column.

t(x):

Transpose a matrix. This is a quick operation which only changes metadata and does not touch the data representation.

rowMaj(x):

Check the data orientation.

Author(s)

Kylie A. Bemis

See Also

matter

Examples

x <- matter_arr(1:1000, dim=c(10,10,10))
x

Description

The matter_fct class implements out-of-memory factors.

Usage

## Instance creation
matter_fct(data, levels, path = NULL,
    length = NA_integer_, names = NULL, offset = 0, extent = NA_real_,
    readonly = NA, append = FALSE, labels = as.character(levels), ...)

## Additional methods documented below

Arguments

Value

An object of class matter_fct.

Slots

data:

This slot stores any information necessary to access the data for the object (which may include the data itself and/or paths to file locations, etc.).

type:

The storage mode of the accessed data when read into R. This is a 'factor' with levels 'raw', 'logical', 'integer', 'numeric', or 'character'.

dim:

Either 'NULL' for vectors, or an integer vector of length one of more giving the maximal indices in each dimension for matrices and arrays.

names:

The names of the data elements for vectors.

dimnames:

Either 'NULL' or the names for the dimensions. If not 'NULL', then this should be a list of character vectors of the length given by 'dim' for each dimension. This is always 'NULL' for vectors.

levels:

The levels of the factor.

labels:

The labels for the levels.

Extends

matter, matter_vec

Creating Objects

matter_fct instances can be created through matter_fct() or matter().

Methods

Standard generic methods:

length(x), length(x) <- value:

Get or set length.

names(x), names(x) <- value:

Get or set 'names'.

x[i], x[i] <- value:

Get or set the elements of the factor.

levels(x), levels(x) <- value:

Get or set the levels of the factor.

Author(s)

Kylie A. Bemis

See Also

matter, matter_vec

Examples

x <- matter_fct(rep(c("a", "a", "b"), 5), levels=c("a", "b", "c"))
x

Description

The matter_list class implements out-of-memory lists.

Usage

## Instance creation
matter_list(data, type = "double", path = NULL,
    lengths = NA_integer_, names = NULL, offset = 0, extent = NA_real_,
    readonly = NA, append = FALSE, ...)

## Additional methods documented below

Arguments

Value

An object of class matter_list.

Slots

data:

This slot stores any information necessary to access the data for the object (which may include the data itself and/or paths to file locations, etc.).

type:

The storage mode of the accessed data when read into R. This is a 'factor' with levels 'raw', 'logical', 'integer', 'numeric', or 'character'.

dim:

Either 'NULL' for vectors, or an integer vector of length one of more giving the maximal indices in each dimension for matrices and arrays.

names:

The names of the data elements for vectors.

dimnames:

Either 'NULL' or the names for the dimensions. If not 'NULL', then this should be a list of character vectors of the length given by 'dim' for each dimension. This is always 'NULL' for vectors.

Extends

matter

Creating Objects

matter_list instances can be created through matter_list() or matter().

Methods

Standard generic methods:

x[[i]], x[[i]] <- value:

Get or set a single element of the list.

x[[i, j]]:

Get the jth sub-elements of the ith element of the list.

x[i], x[i] <- value:

Get or set the ith elements of the list.

lengths(x):

Get the lengths of all elements in the list.

Author(s)

Kylie A. Bemis

See Also

matter

Examples

x <- matter_list(list(c(TRUE,FALSE), 1:5, c(1.11, 2.22, 3.33)), lengths=c(2,5,3))
x[]
x[1]
x[[1]]

x[[3,1]]
x[[2,1:3]]

Description

The matter_str class implements out-of-memory strings.

Usage

## Instance creation
matter_str(data, encoding, path = NULL,
    nchar = NA_integer_, names = NULL, offset = 0, extent = NA_real_,
    readonly = NA, append = FALSE, ...)

## Additional methods documented below

Arguments

Value

An object of class matter_str.

Slots

data:

This slot stores any information necessary to access the data for the object (which may include the data itself and/or paths to file locations, etc.).

type:

The storage mode of the accessed data when read into R. This is a 'factor' with levels 'raw', 'logical', 'integer', 'numeric', or 'character'.

dim:

Either 'NULL' for vectors, or an integer vector of length one of more giving the maximal indices in each dimension for matrices and arrays.

names:

The names of the data elements for vectors.

dimnames:

Either 'NULL' or the names for the dimensions. If not 'NULL', then this should be a list of character vectors of the length given by 'dim' for each dimension. This is always 'NULL' for vectors.

encoding:

The string encoding used.

Extends

matter_list

Creating Objects

matter_str instances can be created through matter_str() or matter().

Methods

Standard generic methods:

x[i], x[i] <- value:

Get or set the string elements of the vector.

lengths(x):

Get the number of characters (in bytes) of all string elements in the vector.

Author(s)

Kylie A. Bemis

See Also

matter

Examples

x <- matter_str(rep(c("hello", "world!"), 50))
x

Description

The matter class and its subclasses are designed for easy on-demand read/write access to binary virtual memory data structures, and working with them as vectors, matrices, arrays, lists, and data frames.

Usage

## Instance creation
matter(...)

# Check if an object is a matter object
is.matter(x)

# Coerce an object to a matter object
as.matter(x)

## Additional methods documented below

Arguments

Value

An object of class matter.

Slots

data:

This slot stores any information necessary to access the data for the object (which may include the data itself and/or paths to file locations, etc.).

type:

The storage mode of the accessed data when read into R. This is a 'factor' with levels 'raw', 'logical', 'integer', 'numeric', or 'character'.

dim:

Either 'NULL' for vectors, or an integer vector of length one of more giving the maximal indices in each dimension for matrices and arrays.

names:

The names of the data elements for vectors.

dimnames:

Either 'NULL' or the names for the dimensions. If not 'NULL', then this should be a list of character vectors of the length given by 'dim' for each dimension. This is always 'NULL' for vectors.

Creating Objects

matter is a virtual class and cannot be instantiated directly, but instances of its subclasses can be created through matter().

Methods

Class-specific methods:

atomdata(x):

Access the 'data' slot.

adata(x):

An alias for atomdata(x).

type(x), type(x) <- value:

Get or set data 'type'.

Standard generic methods:

length(x), length(x) <- value:

Get or set length.

dim(x), dim(x) <- value:

Get or set 'dim'.

names(x), names(x) <- value:

Get or set 'names'.

dimnames(x), dimnames(x) <- value:

Get or set 'dimnames'.

Author(s)

Kylie A. Bemis

See Also

matter_arr, matter_mat, matter_vec, matter_fct, matter_list, matter_str

Examples

## Create a matter_vec vector
x <- matter(1:100, length=100)
x

## Create a matter_mat matrix
y <- matter(1:100, nrow=10, ncol=10)
y

Description

The matter package provides the following options:

• options(matter.compress.atoms=3): The compression ratio threshold to be used to determine when to compress atoms in a matter object. Setting to 0 or FALSE means that atoms are never compressed.

• options(matter.default.nchunks=20L): The default number of chunks to use when iterating over matter objects.

• options(matter.default.verbose=FALSE): The default verbosity for printing progress messages.

• options(matter.matmul.bpparam=NULL): An optional BiocParallelParam passed to bplapply used when performing matrix multiplication with matter_mat and sparse_mat objects.

• options(matter.show.head=TRUE): Should a preview of the beginning of the data be displayed when the object is printed?

• options(matter.show.head.n=6): The number of elements, rows, and/or columns to be displayed by the object preview.

• options(matter.coerce.altrep=FALSE): When coercing matter objects to native R objects (such as matrix), should a matter-backed ALTREP object be returned instead? The initial coercion will be cheap, and the result will look like a native R object. This does not guarantee that the full data is never read into memory. Not all functions are ALTREP-aware at the C-level, so some operations may still trigger the full data to be read into memory. This should only ever happen once, as long as the object is not duplicated, though.

• options(matter.wrap.altrep=FALSE): When coercing to a matter-backed ALTREP object, should the object be wrapped in an ALTREP wrapper? (This is always done in cases where the coercion preserves existing attributes.) This allows setting of attributes without triggering a (potentially expensive) duplication of the object when safe to do so.

• options(matter.dump.dir=tempdir()): Temporary directory where matter object files should be dumped when created without user-specified file paths.

Description

The matter package defines a number of data types for translating between data elements stored in virtual memory and data elements loaded into R. These are typically set and stored via the datamode argument and slot.

At the R level, matter objects may be any of the following data modes:

• raw: matter objects of this mode are typically vectors of raw bytes.

• logical: Any matter object that represents a logical vector or has had any Compare or Logic delayed operations applied to it will be of this type.

• integer: matter objects represented as integers in R.

• numeric: matter objects represented as doubles in R.

• character: matter objects representated as character vectors in R.

In virtual memory, matter objects may be composed of atomic units of the following data types:

• char: 8-bit signed integer; defined as char.

• uchar: 8-bit unsigned integer; used for ‘Rbyte’ or ‘raw’; defined as unsigned char.

• int16: 16-bit signed integer; defined as int16_t. May be aliased as ‘short’ and ‘16-bit integer’.

• uint16: 16-bit unsigned integer; defined as uint16_t. May be aliased as ‘ushort’ and ‘16-bit unsigned integer’.

• int32: 32-bit signed integer; defined as int32_t. May be aliased as ‘int’ and ‘32-bit integer’.

• uint32: 32-bit unsigned integer; defined as uint32_t. May be aliased as ‘uint’ and ‘32-bit unsigned integer’.

• int64: 64-bit signed integer; defined as int64_t. May be aliased as ‘long’ and ‘64-bit integer’.

• uint64: 64-bit unsigned integer; defined as uint64_t. May be aliased as ‘ulong’ and ‘64-bit unsigned integer’.

• float32: 32-bit float; defined as float. May be aliased as ‘float’ and ‘32-bit float’.

• float64: 64-bit float; defined as double. May be aliased as ‘double’ and ‘64-bit float’.

While a substantial effort is made to coerce data elements properly between data types, sometimes this cannot be done losslessly. Loss of precision is silent, while values outside of the representable range will generate a warning (sometimes many such warnings) and will be set to NA if available or 0 otherwise.

Note that the unsigned data types do not support NA; coercion between signed integer types attempts to preserve missingness. The special values NaN, Inf, and -Inf are only supported by the floating-point types, and will be set to NA for signed integral types, and to 0 for unsigned integral types.

Description

These are utility functions for checking memory used by objects and by R during the execution of an expression.

Usage

mem(x, reset = FALSE)

mtime(expr)

Arguments

Details

These are wrappers around the built-in gc and link{proc.time} functions. Note that they only count memory managed by R.

Value

For mtime, a vector giving [1] the amount of memory used at the start of execution, [2] the amount of memory used at the end of execution, [3] the maximum amount of memory used during execution, [4] the memory overhead as defined by the maximum memory used minus the starting memory use, and [5] the execution time in seconds.

For mem, either a single numeric value giving the memory used by an object, or a vector providing a more readable version of the information returned by gc (see its help page for details).

Author(s)

Kylie A. Bemis

See Also

gc,

Examples

x <- 1:100

mem(x)

mtime(mean(x + 1))

Description

Multiple instance learning is a strategy for training classifiers when the class labels are observed at a coarser level than the individual data points. For example, if an entire image is classified as "positive" or "negative" but the classifier is trained and predicts at the pixel level.

Usage

mi_learn(fn, x, y, bags, pos = 1L, ...,
	score = fitted, threshold = 0.01, verbose = NA)

Arguments

Details

This is a generic wrapper for applying a multiple instance learning strategy for any classifier that satisfies certain criteria. The labels must be binary (positive and negative).

The multiple instance learning algorithm here assumes that if a single data point is positive, then the entire bag to which it belongs is labeled as positive. For example, if a single pixel in an image indicates the presence of disease, then the entire image is labeled as disease.

The model returned by fn must support returning either a vector of probabilities or a 2-column score matrix when passed to score.

Value

A model object returned by fn.

Author(s)

Kylie A. Bemis

References

D. Guo, M. C. Foell, V. Volkmann, K. Enderle-Ammour, P. Bronsert, O. Shilling, and O. Vitek. “Deep multiple instance learning classifies subtissue locations in mass spectrometry images from tissue-level annotations.” Bioinformatics, vol. 36, issue Supplement_1, pp. i300-i308, 2020.

See Also

nscentroids

Examples

register(SerialParam())
set.seed(1)
n <- 100
p <- 5
g <- 5
bags <- rep(paste0("s", seq_len(g)), each=n %/% g)
bags <- factor(rep_len(bags, n))
x <- matrix(rnorm(n * p), nrow=n, ncol=p)
colnames(x) <- paste0("x", seq_len(p))

# create bagged labels
y <- ifelse(bags %in% c("s1", "s2", "s3"), "pos", "neg")
y <- factor(y, levels=c("pos", "neg"))
ipos <- which(y == "pos")
ineg <- which(y == "neg")
z <- y

# create "true" labels (with some within-bag noise)
z[ipos] <- sample(c("pos", "neg"), length(ipos), replace=TRUE)
jpos <- which(z == "pos")
jneg <- which(z == "neg")

# create data
x[jpos,] <- x[jpos,] + rnorm(p * length(jpos), mean=1)
x[jneg,] <- x[jneg,] - rnorm(p * length(jneg), mean=1)

# fit ordinary NSC and mi-NSC
fit0 <- nscentroids(x=x, y=y)
fit1 <- mi_learn(nscentroids, x=x, y=y, bags=bags, priors=1)

# improved performance on "true" labels
mean(fitted(fit0, "class") == z)
mean(fitted(fit1, "class") == z)

Description

Nonnegative matrix factorization (NMF) decomposes a nonnegative data matrix into a matrix of basis variables and a matrix of activations (or coefficients). The factorization is approximate and may be less accurate than alternative methods such as PCA, but can greatly improve the interpretability of the reduced dimensions.

Usage

# Alternating least squares
nnmf_als(x, k = 3L, s = 1e-9, transpose = FALSE,
	niter = 100L, tol = 1e-5, verbose = NA, ...)

# Multiplicative updates
nnmf_mult(x, k = 3L, s = 1e-9, cost = c("euclidean", "KL", "IS"),
	transpose = FALSE, niter = 100L,
	tol = 1e-5, verbose = NA, ...)

## S3 method for class 'nnmf'
predict(object, newdata, ...)

# Nonnegative double SVD
nndsvd(x, k = 3L, ...)

Arguments

Details

These functions implement nonnegative matrix factorization (NMF) using either alternating least squares as described by Berry at al. (2007) or multiplicative updates from Lee and Seung (2000) and further described by Burred (2014). The algorithms are initialized using nonnegative double singular value decomposition (NNDSVD) from Boutsidis and Gallopoulos (2008).

The algorithm using multiplicative updates (nnmf_mult()) tends to be more stable but converges more slowly. The alternative least squares algorithm (nnmf_als()) tends to converge faster to more accurate results, but can be less numerically stable than the multiplicative updates algorithm.

Note for nnmf_mult() that method = "euclidean" is the only method that can handle out-of-memory matter_mat and sparse_mat matrices. x will be coerced to an in-memory matrix for other methods.

Value

An object of class nnmf, with the following components:

• activation: The (transposed) coefficient matrix (H).

• x: The basis variable matrix (W).

• iter: The number of iterations performed.

Author(s)

Kylie A. Bemis

References

M. W. Berry, M. Browne, A. N. Langville, V. P. Pauca, R. J. Plemmons. “Algorithms and applications for approximate nonnegative matrix factorization.” Computational Statistics and Data Analysis, vol. 52, issue 1, pp. 155-173, Sept. 2007.

D. D. Lee and H. S. Seung. “Algorithms for non-negative matrix factorization.” Proceedings of the 13th International Conference on Neural Information Processing Systems (NIPS), pp. 535-541, Jan. 2000.

C. Boutsidis and E. Gallopoulos. “SVD based initialization: A head start for nonnegative matrix factorization.” Pattern Recognition, vol. 41, issue 4, pp. 1350-1362, Apr. 2008.

J. J. Burred. “Detailed derivation of multiplicative update rules for NMF.” Techical report, Paris, March 2014.

See Also

svd, prcomp

Examples

set.seed(1)

a <- matrix(sort(runif(500)), nrow=50, ncol=10)
b <- matrix(rev(sort(runif(500))), nrow=50, ncol=10)
x <- cbind(a, b)

mf <- nnmf_als(x, k=3)

Description

Nearest shrunken centroids performs regularized classification of high-dimensional data. Originally developed for classification of microarrays, it calculates test statistics for each feature/dimension based on the deviation between the class centroids and the global centroid. It applies regularization (via soft thresholding) to these test statistics to produce shrunken centroids for each class.

Usage

# Nearest shrunken centroids
nscentroids(x, y, s = 0, distfun = NULL,
	priors = table(y), center = NULL, transpose = FALSE,
	verbose = NA, nchunks = NA, BPPARAM = bpparam(), ...)

## S3 method for class 'nscentroids'
fitted(object, type = c("response", "class"), ...)

## S3 method for class 'nscentroids'
predict(object, newdata,
	type = c("response", "class"), ...)

## S3 method for class 'nscentroids'
logLik(object, ...)

Arguments

Details

This functions implements nearest shrunken centroids based on the original algorithm by Tibshirani et al. (2002). It provides a sparse strategy for classification based on regularized class centroids. The class centroids are shrunken toward the global centroid. The shrunken test statistics used to perform the regularization can then be interpreted to determine which features are relevant to the classification. (Important features will have nonzero test statitistics after soft thresholding.)

Unlike the original algorithm, this implementation allows specifying a custom dissimilarity function. If not provided, then this defaults to rowDistFun() if transpose=FALSE or colDistFun() if transpose=TRUE.

If a custom function is passed, it should take the form function(x, y, ...), and it must return a function of the form function(i). The returned function should return the distances between the ith object(s) in x and all objects in y. In addition, it must support an argument called weights that takes a vector of feature weights used to scale the features during the distance calculation. rowDistFun() and colDistFun() are examples of functions that satisfy these properties.

Value

An object of class nscentroids, with the following components:

• class: The predicted classes.

• probability: A matrix of posterior class probabilities.

• centers: The shrunken class centroids used for classification.

• statistic: The shrunken test statistics.

• sd: The pooled within-class standard deviations for each feature.

• priors: The prior class probabilities.

• s: The regularization (soft thresholding) parameter.

• distfun: The function used to generate the dissimilarity function.

Author(s)

Kylie A. Bemis

References

R. Tibshirani, T. Hastie, B. Narasimhan, and G. Chu. “Diagnosis of multiple cancer types by shrunken centroids of gene expression.” Proceedings of the National Academy of Sciences of the USA, vol. 99, no. 10, pp. 6567-6572, 2002.