Package {cooltools}

Type:	Package
Title:	Practical Tools for Scientific Computation and Visualisation
Version:	2.33
Description:	Provides utilities for scientific computation and visualisation, with an emphasis on applications in physics and astrophysics. Functionality includes random sampling from spherical and custom distributions, information and entropy analysis, Fourier transforms, two-point correlation estimation, binning and gridding of point sets, two-dimensional interpolation, Monte Carlo integration, vector operations, coordinate transformations, physical constants, and cosmological conversions. Graphics tools support the creation and export of publication-quality plots, animations, colour scales, map projections, and bitmap images. Several of these tools were used by Obreschkow et al. (2020) <doi:10.1093/mnras/staa445>.
Imports:	plotrix, celestial, data.table, pracma, utils, png, jpeg, MASS, raster, sp, cubature, bit64, randtoolbox, Rcpp, FNN, gitcreds, pak
Suggests:	EBImage, hdf5r
LinkingTo:	Rcpp
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
NeedsCompilation:	yes
URL:	https://github.com/obreschkow/cooltools
Config/roxygen2/version:	8.1.0
Packaged:	2026-09-11 09:44:23 UTC; 00077389
Author:	Danail Obreschkow [aut, cre]
Maintainer:	Danail Obreschkow <danail.obreschkow@gmail.com>
Repository:	CRAN
Date/Publication:	2026-09-11 10:20:08 UTC

Practical Tools for Scientific Computations and Visualizations

Description

Collection of routines for efficient scientific computations in physics and astrophysics. These routines include utility functions, advanced computation tools, as well as visualisation tools. They can be used, for example, for generating random numbers from spherical and custom distributions, information and entropy analysis, special Fourier transforms, two-point correlation estimation (e.g. as in Landy & Szalay (1993) <doi:10.1086/172900>), binning & gridding of point sets, 2D interpolation, Monte Carlo integration, vector arithmetic and coordinate transformations. Also included are a non-exhaustive list of important constants and cosmological conversion functions. The graphics routines can be used to produce and export publication-ready scientific plots and movies, e.g. as used in Obreschkow et al. (2020) <doi:10.1093/mnras/staa445>. These routines include special color scales, projection functions, and bitmap handling routines.

Author(s)

Danail Obreschkow

See Also

Useful links:

• https://github.com/obreschkow/cooltools

Package environment

Description

Environment used to store global variables in the package, e.g. used for subplot routine.

Usage

.cooltools.env

Value

None

Author(s)

Danail Obreschkow

Associated Legendre Polynomials

Description

Compute associated Legendre polynomials P_l^m(x), defined as the canonical solutions of the general Legendre equation. These polynomials are used, for instance, to compute the spherical harmonics.

Usage

alp(x, l = 0, m = 0)

Arguments

Value

Returns a vector with the same number of elements as x

Author(s)

Danail Obreschkow

See Also

sphericalharmonics

Compute the total angular momentum of a particle distribution

Description

Computes the total angular momentum of a set of particles about their centre of mass (or about a user-specified origin).

Usage

angularmomentum(m, x, v, x0 = NULL, v0 = NULL)

Arguments

Value

A numeric vector of length three containing the Cartesian components of the total angular momentum,

\mathbf{J}=\sum_i m_i\,(\mathbf{x}_i-\mathbf{x}_0)\times (\mathbf{v}_i-\mathbf{v}_0).

If m is NULL, a vector of zeros is returned.

Examples

m = c(1, 2)
x = rbind(c(1, 0, 0), c(0, 1, 0))
v = rbind(c(0, 1, 0), c(-1, 0, 0))
angularmomentum(m, x, v)

Bilinear interpolation function of data on a regular grid

Description

Generates a fast function f(x,y) that interpolates gridded data, based on the analogous subroutine approxfun in 1D.

Usage

approxfun2(x, y, z, outside = NA)

Arguments

Value

Returns a fast and vectorized interpolation function f(x,y)

Author(s)

Danail Obreschkow

See Also

approxfun

Examples

x = seq(3)
y = seq(4)
z = array(c(x+1,x+2,x+3,x+4),c(3,4))
f = approxfun2(x,y,z)
print(f(1.7,2.4))

Bin two-dimensional data in one dimension

Description

Divides a vector of values x into finite intervals; returns the counts and other statistics in each interval.

Usage

bindata(x, y = NULL, bins = 20, method = "regular", xlim = NULL)

Arguments

Value

Returns a list of items

Author(s)

Danail Obreschkow

See Also

griddata

Examples

# make and plot 100 random (x,y)-points
set.seed(1)
x = runif(200)
y = x+rnorm(200)
plot(x,y,pch=16,cex=0.5)

# bin the data into 10 bins of 20 points each
bin = bindata(x,y,10,'equal')
segments(bin$xleft,bin$y,bin$xright,bin$y,col='red')
segments(bin$x,bin$y16,bin$x,bin$y84,col='red')
segments(bin$x,bin$y-bin$yerr,bin$x,bin$y+bin$yerr,col='red',lwd=3)
points(bin$x,bin$y,pch=16,col='red')

Gaussian blur for arrays of arbitrary dimension

Description

Smooth a numeric array of arbitrary dimension with a separable Gaussian kernel. This function is conceptually similar to gblur from the EBImage package (which is limited to 2D), but works for arrays of any dimension.

Usage

blur(x, sigma, boundary = "replicate")

Arguments

Details

The Gaussian kernel is separable, so smoothing is performed via successive 1D convolutions along each dimension.

This function is similar in spirit to gblur from EBImage, but generalises naturally to arrays of arbitrary dimension.

Boundary handling can be controlled via 'boundary':

• '"none"' uses centered filtering without padding and therefore returns 'NA' values near array boundaries.

• '"circular"' applies periodic wrap-around boundary conditions.

• '"replicate"' extends the array by repeating the edge values.

For finite arrays, the Gaussian kernel is truncated as needed so that its length does not exceed the array extent along any dimension.

Value

A numeric array with the same dimensions as 'x', containing the blurred values.

Cartesian to polar/cylindrical coordinate conversion

Description

Convert 2D/3D Cartesian to polar/cylindrical coordinates

Usage

car2pol(x)

Arguments

Value

Returns a 2/3-vector or a n-by-2/3 matrix representing the polar/cylindrical coordinates r,phi,(z), where phi=0...2*pi is the azimuth measured positively from the x-axis.

Author(s)

Danail Obreschkow

See Also

pol2car

Cartesian to spherical coordinate conversion

Description

Convert 3D Cartesian to spherical coordinates

Usage

car2sph(x)

Arguments

Value

Returns a 3-vector or a n-by-3 element matrix representing the spherical coordinates (r,theta,phi), where theta=0...pi is the polar angle measured from the north pole and phi=0...2*pi is the azimuth measured positively from the x-axis (ISO 80000-2:2019 physics convention).

Author(s)

Danail Obreschkow

See Also

sph2car

Convert complex numbers to color

Description

Converts a complex number (or a vector/array thereof) into a color-string, such that brightness represents the complex amplitude and hue represents the complex phase.

Usage

cmplx2col(z, max = NULL, hue = 0, saturation = 1, gamma = 1)

Arguments

Value

Returns a single color or vector/array of colors

Author(s)

Danail Obreschkow

Vertical color bar

Description

Adds a vertical color bar to a plot with a custom axis.

Usage

colorbar(
  xleft,
  ybottom,
  xright,
  ytop,
  col = gray.colors(256, 0, 1),
  clim = c(0, 1),
  show.border = TRUE,
  text = "",
  line = 2,
  show.axis = TRUE,
  side = "right",
  lwd = 1,
  nticks = 5,
  at = NULL,
  srt = 0,
  ticklength = 0.1,
  shift = 0,
  ...
)

Arguments

Value

None

Author(s)

Danail Obreschkow

Examples

## Plot a spherical function with color bar
nplot(xlim=c(0,1.2), asp=1)
f = function(theta,phi) cos(10*theta+phi)
sp = sphereplot(f, 200, col=planckcolors(200), phi0=0.1, theta0=pi/3,
add=TRUE, center=c(0.5,0.5), radius=0.4, clim=c(-1,1))
colorbar(1,0.1,1.1,0.9,col=sp$col,clim=sp$clim)

Highest-density contour levels

Description

Computes density thresholds defining highest-density regions.

For each probability p, the function returns a level \ell such that the region where the density is at least \ell contains a fraction p of the total probability mass:

\frac{\int_{f(\mathbf{x}) \geq \ell} f(\mathbf{x})\,d\mathbf{x}} {\int f(\mathbf{x})\,d\mathbf{x}} = p.

The input may be a numeric vector or array containing sampled density values, or a function representing a density over a bounded rectangular domain. The density does not need to be normalised, but it must be non-negative and have a finite, positive integral or sum.

Usage

contourlevel(
  f,
  p = c(0.6826895, 0.9544997),
  xmin = NULL,
  xmax = NULL,
  neval = 10000,
  subdivisions = 1000,
  napprox = 30,
  rel.tol = 1e-05,
  abs.tol = 0,
  ...
)

Arguments

Details

The returned contours define highest-density regions: points are included in decreasing order of density until the requested probability mass is enclosed. For multimodal densities, the resulting region may consist of several disconnected components.

When f is a function, the enclosed mass above a trial level is evaluated by numerically integrating

f(\mathbf{x}) I[f(\mathbf{x})\geq\ell].

This integrand is discontinuous at the contour boundary, so convergence may be slower than for a smooth integrand, particularly in high dimensions or for complicated contours.

One-dimensional densities are integrated using integrate. Multidimensional densities are integrated using cuhre.

If napprox>0, the enclosed-mass function is evaluated on a grid of density levels and inverted by monotonic linear interpolation. If napprox=0, each requested probability is solved separately using uniroot. The latter is generally slower but avoids the interpolation approximation.

Value

A numeric vector of density levels with the same length and ordering as p.

For a vector or array, all entries are assumed to represent cells of equal volume or equal statistical weight. Because the enclosed mass changes in discrete steps, the returned level is the sampled density threshold whose superlevel set first contains at least the requested probability mass.

Author(s)

Danail Obreschkow

See Also

dpqr

Examples

## f(x) is a one-dimensional PDF
# Compute the one- and two-sigma contour levels of a normal distribution,
# i.e. the values l such that
# integral over dnorm(x) >= l of dnorm(x) dx = p,
# where p = 68.3% and 95.4%.
l = contourlevel(dnorm, xmin = -10, xmax = 10, napprox = 0)
print(l)

# Compare these values with dnorm(1) and dnorm(2)
print(dnorm(c(1, 2)))


## f(x) is a two-dimensional likelihood function
# Produce 20%, 40%, 60%, and 80% highest-density contours.
f = function(x) {
  cos(2*x[1]-x[2]-1)^2*exp(-x[1]^2-x[2]^2-x[1]*x[2])
}

p = c(0.2, 0.4, 0.6, 0.8)

# Values l such that
# integral over f(x) >= l of f(x) dx = p * integral f(x) dx
l = contourlevel(f, p, c(-5, -5), c(5, 5))

# Plot the function and contours at the levels l
x = seq(-3, 3, length.out = 200)
m = pracma::meshgrid(x)
z = array(Vectorize(function(x, y) f(c(x, y)))(m$Y, m$X), dim(m$X))

image(x, x, z, col = terrain.colors(100))
contour(x, x, z, levels = l, add = TRUE,
        labels = sprintf("%.0f%%", p*100), labcex = 0.7)


## f is a 20-by-20 array representing a gridded point set
# Produce 1000 points drawn from a two-dimensional normal distribution.
set.seed(1)
x = MASS::mvrnorm(n = 1000, mu = c(0, 0), Sigma = matrix(c(3, 1, 1, 2), 2, 2))

# Grid these points onto a regular 20-by-20 grid
g = griddata(x, min = -6, max = 6)

# Find one- and two-sigma contour levels and draw the contours
l = contourlevel(g$field)

plot(x, xlim = g$grid[[1]]$lim, ylim = g$grid[[2]]$lim, pch = 20, cex = 0.5)
contour(g$grid[[1]]$mid, g$grid[[2]]$mid, g$field,
        levels = l, add = TRUE, col = "red", lwd = c(2, 1), labels = NA)

Fast cosmology conversion functions

Description

Generates all 20 conversion functions between redshift (z), luminosity distance (dl), comoving distance (dc) and angular diameter distance (da), and lookback time (t = light travel time from specified redshift); based on the *celestial* package.

Usage

cosmofct(zmin = 0, zmax = 1, dz = 0.02, H0 = 70, OmegaM = 0.3, ...)

Arguments

Value

Returns a list of 20 vectorized functions; e.g. dc2z to convert from comoving distance to redshift. Also contains the age of the universe at z=0. All distances are in units of Mpc and times are in units of Gyr.

Author(s)

Danail Obreschkow (based on *celestial* package by Aaron Robotham)

Examples

## uses a flat LCDM cosmology with h=0.68, OmegaM=0.32 and OmegaL=0.68
cosmo = cosmofct(0,1,H0=68,OmegaM=0.32)
curve(cosmo$z2dl(x),0,1,xlab='z',ylab='distance',col='red')
curve(cosmo$z2dc(x),0,1,col='black',add=TRUE)
curve(cosmo$z2da(x),0,1,col='blue',add=TRUE)
d = seq(500,5000,500)
points(cosmo$dl2z(d),d,pch=16,col='red')
points(cosmo$dc2z(d),d,pch=16,col='black')
points(cosmo$da2z(d),d,pch=16,col='blue')

Circularly shift each dimension of an array

Description

Circulates each dimension of an array. This routine is identical to circshift, but works with arrays up to 5 dimensions.

Usage

cshift(x, s)

Arguments

Value

Returns a vector or array of the same shape as x.

Author(s)

Danail Obreschkow

Scientific constants

Description

The list cst contains useful scientific constants in SI units, mainly for astrophysics.

Usage

cst

Value

None

Author(s)

Danail Obreschkow

Examples

# print all the constants to console
for (i in seq(length(cst))) cat(sprintf('%6s = %.12e\n',names(cst)[i],cst[i]))

Cube Helix colour palette

Description

Generate the Cube Helix colour palette, designed to appropriately display of intensity images, as the brightness increases monotonically when displayed in greyscale.

Usage

cubehelix(n, r = 1.5, hue = 1, gamma = 1, rev = FALSE)

Arguments

Details

This scheme was published by Green, D. A., 2011, "A colour scheme for the display of astronomical intensity images." Bulletin of the Astronomical Society of India, 39, 289.

Value

Returns an n-vector of RGB colour strings.

Author(s)

Danail Obreschkow

Discrete Fourier Transform

Description

Discrete Fourier Transform (DFT) with longest modes at the center in Fourier space and normalized such that dft(dft(f),inverse)=f. This is the discretization scheme described in Appendix D of Obreschkow et al. 2013, ApJ 762. Relies on fft.

Usage

dft(f, inverse = FALSE, shift = -floor(dim(as.array(f))/2), simplify = TRUE)

Arguments

Value

Returns an array of the same shape as f, containing the (inverse) Fourier Transform.

Author(s)

Danail Obreschkow

See Also

fft

Examples

## DFT of a 2D normal Gaussian function
n = 30
f = array(0,c(n,n))
for (i in seq(n)) {
  for (j in seq(n)) f[i,j] = exp(-(i-6)^2/4-(j-8)^2/2-(i-6)*(j-8)/2)
}
plot(NA,xlim=c(0,2.1),ylim=c(0,1.1),asp=1,bty='n',xaxt='n',yaxt='n',xlab='',ylab='')
rasterImage(f,0,0,1,1,interpolate=FALSE)
g = dft(f)
img = array(hsv((pracma::angle(g)/2/pi)%%1,1,abs(g)/max(abs(g))),c(n,n))
rasterImage(img,1.1,0,2.1,1,interpolate=FALSE)
text(0.5,1,'Input function f',pos=3)
text(1.6,1,'DFT(f)',pos=3)

Produce coordinates for Discrete Fourier Transform

Description

Produces the direct space coordinates (x) and Fourier space frequencies (f) and angular frequencies (k), corresponding to the Discrete Fourier Transform dft of this package.

Usage

dftgrid(N, L, x0 = 0, k0 = -floor(N/2) * 2 * pi/L)

Arguments

Value

Returns a list with:

Author(s)

Danail Obreschkow

Radius of a point distribution

Description

Computes a characteristic radius of a set of points about a centre.

Usage

distradius(x, m = NULL, x0 = NULL, type = c("max", "mean", "median"))

Arguments

Value

A scalar radius.

See Also

vectornorm

Generate the d/p/q/r family for an arbitrary distribution

Description

Generates the standard d, p, q, and r functions associated with an arbitrary one-dimensional probability density function, analogous to dnorm/pnorm/qnorm/rnorm and related families.

Usage

dpqr(fun, min, max)

Arguments

Value

Returns a list with the following functions:

Author(s)

Danail Obreschkow

See Also

rng, contourlevel

Examples

f = function(x) sin(x)
dist = dpqr(f, 0, pi)

x = dist$r(1000)
hist(x, probability = TRUE)
curve(dist$d(x), add = TRUE)

Information entropy

Description

Computes the information entropy H=sum(p*log_b(p)), also known as Shannon entropy, of a probability vector p.

Usage

entropy(p, b = exp(1), normalize = TRUE)

Arguments

Value

Returns the information entropy in units that depend on b. If b=2, the units are bits; if b=exp(1), the units are nats; if b=10, the units are dits.

Author(s)

Danail Obreschkow

Draw a line with uncertainty regions

Description

Draw a line with uncertainty regions

Usage

errlines(
  x,
  y,
  errp,
  errn = errp,
  col = "black",
  alpha = 0.5,
  smooth = FALSE,
  df = NULL,
  ...
)

Arguments

Value

None

Author(s)

Danail Obreschkow

Stop code and produce an error message

Description

Stops the code and produces an error message, like the stop() function of the base library, but terminating a timer started with tick, if running.

Usage

error(...)

Arguments

Value

None. geterrmessage() returns the last error message, as a character string.

Author(s)

Danail Obreschkow

See Also

tick tock progress

Evenly distributed n points on a sphere

Description

Distributes n points on a sphere in a relatively even fashion following the generalised Fibonacci algorithm, described at http://extremelearning.com.au/evenly-distributing-points-on-a-sphere/

Usage

fibonaccisphere(n = 1000, r = 1, out.xyz = TRUE, out.sph = FALSE)

Arguments

Value

matrix of n points (in n rows), in Cartesian and/or spherical coordinates.

Author(s)

Danail Obreschkow

See Also

runif3

Examples

## plot standard projections of a 1000-point Fibonacci sphere
xyz = fibonaccisphere()
plot(xyz, asp=1, pch=16, cex=0.5)

Install a GitHub Package

Description

Installs an R package from GitHub using pak::pak(). If installation fails with GitHub HTTP error 401 / bad credentials, the stored GitHub credentials are deleted using gitcreds::gitcreds_delete(), and the installation is attempted once more.

Usage

getgit(reponame, ...)

Arguments

Value

Invisibly returns the result of pak::pak().

Examples

## Not run: 
getgit("obreschkow/cooltools")

## End(Not run)

Compute gradient

Description

Evaluates the gradient of a vector of array of rank 2 or 3

Usage

gradient(f, circular = FALSE)

Arguments

Value

returns a list with the different components of the gradient; each component has the same dimension as f.

Author(s)

Danail Obreschkow

Gaussian Random Field generator

Description

Generates a scalar Gaussian Random Field (GRF) in 1, 2 or 3 dimensions, with a power-law power spectrum of custom slope alpha. The field is normalized such that its mean is zero (up to floating point errors) and its expected standard deviation is one, for alpha=0. The Fourier phases are sampled in order of increasing frequency, such that the random structure is preserved when changing the output size.

Usage

grf(nside = 100, dim = 2, alpha = 0, seed = NULL)

Arguments

Value

Returns a vector, matrix or 3D-array with nside elements per dimension

Author(s)

Danail Obreschkow

Examples

# Generate the same 2D GRF in two different resolutions
nplot(c(0,2.1), asp=1)
lowres = grf(nside=30, alpha=-2, seed=1)
graphics::rasterImage(stretch(lowres),0,0,1,1)
highres = grf(nside=120, alpha=-2, seed=1)
graphics::rasterImage(stretch(highres),1.1,0,2.1,1)

# Check the power spectrum of a general GRF
nside = 50 # change this to any integer >1
alpha = -1.7 # change this to any power
dim = 3 # change this to any positive integer (but keep nside^dim reasonable)
x = grf(nside=nside, dim=dim, alpha=alpha)
fourier.grid = dftgrid(N=nside,L=1)
knorm = vectornorm(expand.grid(rep(list(fourier.grid$k),dim)))
power = abs(dft(x))^2
b = bindata(knorm,power,method='equal')
plot(b$xmedian,b$ymedian,log='xy',pch=16,
     xlab='Fourier mode |k|',ylab='Power p(k)',main='Power spectrum')
graphics::curve(x^alpha/nside^dim,col='blue',add=TRUE) # expected power-law

Distribute a point set onto a regular grid

Description

Distributes a set of points in D dimensions onto a regular, D-dimensional grid, using a fast nearest-neighbor algorithm. Weights can be used optionally.

Usage

griddata(
  x,
  w = NULL,
  n = 10,
  min = NULL,
  max = NULL,
  type = c("counts", "density", "probability", "mean")
)

Arguments

Value

Returns a list of items

Author(s)

Danail Obreschkow

Examples

# On some machines data.table does not perform well in multi-threading mode.
# The following command limits this to one thread as a precaution, but
# may not be needed on most machines.
data.table::setDTthreads(1L)

# Distribute 1-dimensional data onto a regular grid
npoints = 1e4
x = rnorm(npoints)
g = griddata(x,min=-3,max=3,n=100,type='probability')
curve(dnorm(x),-3,3)
points(g$grid$mid,g$field,pch=16)

# Distribute 2-dimensional data onto a regular grid
x = runif(100,max=2)
y = runif(100)
g = griddata(cbind(x,y),min=c(0,0),max=c(2,1),n=c(20,10))
image(g$grid[[1]]$breaks,g$grid[[2]]$breaks,g$field,
      asp=1,col=grey.colors(100,0,1),xlab='x',ylab='y')
points(x,y,col='red',pch=16)

# ... same with weights
w = runif(100)
g = griddata(cbind(x,y),w,min=c(0,0),max=c(2,1),n=c(20,10))
image(g$grid[[1]]$breaks,g$grid[[2]]$breaks,g$field,
      asp=1,col=grey.colors(100,0,1),xlab='x',ylab='y')
points(x,y,col='red',pch=16,cex=w)

# ... mean weight in each cell (dark red, where no data)
par(bg = "darkred")
g = griddata(cbind(x,y),w,min=c(0,0),max=c(2,1),n=c(20,10),type='mean')
image(g$grid[[1]]$breaks,g$grid[[2]]$breaks,g$field,
      asp=1,col=grey.colors(100,0,1),xlab='x',ylab='y')
points(x,y,col='red',pch=16,cex=w)

Generate histogram coordinates from mid points

Description

Converts the mid-point x-values and mean densities of binned data into (x,y)-coordinates of a histogram.

Usage

histcoord(x, y, yleft = 0, yright = 0)

Arguments

Value

(2n+2)-by-2 matrix of (x,y)-coordinates to draw histogram as a connected line

Author(s)

Danail Obreschkow

Examples

x = seq(5)
y = sin(x)
plot(x,y,xlim=c(0,6))
lines(histcoord(x,y))

Moment of inertia tensor

Description

Computes the moment of inertia tensor of a set of point masses in arbitrary dimensions,

\mathbf{I} = \sum_i m_i \left( |\mathbf{x}_i|^2\mathbf{1} - \mathbf{x}_i\mathbf{x}_i^\mathrm{T} \right),

where \mathbf{1} is the identity matrix.

Usage

inertia(x, m = 1)

Arguments

Value

A symmetric d\times d matrix containing the moment of inertia tensor, where d is the dimension of the points.

See Also

moments, quadrupole

Invert and shift colors of an image

Description

Invert the brightness of each color channel in an image and/or circularly shifts the hue value. Optionally, a Gaussian blur can be applied.

Usage

invert(
  img = NULL,
  invert = TRUE,
  colshift = 0,
  blur = 0,
  file.in = "",
  file.out = "",
  format = "png",
  show.image = TRUE
)

Arguments

Value

Returns an n-by-m-by-3 array or n-by-m-by-4 array of the processed image.

Author(s)

Danail Obreschkow

Examples

img = yinyangyong # this is an example image included in the package

# invert brightness of all channels
invert(img)

# invert brightness, but preserve hue
invert(img, colshift=0.5)

Numerical equality check

Description

Checks if two or more numerical values are identical within a relative numerical tolerance

Usage

is.equal(x, ..., stoptext = NULL, eps = sqrt(.Machine$double.eps))

Arguments

Value

Returns a logical value. TRUE means that all values of x are equal within the specified relative tolerance; also returns TRUE if all values are Inf or NA or NaN.

Author(s)

Danail Obreschkow

Examples

# almost identical values
x = c(acos(1/sqrt(2)),pi/4)
print(x)
print(x[1]==x[2])
print(is.equal(x))

# various other examples
print(is.equal(1,2,3))
print(is.equal(1,NA,3))
print(is.equal(Inf,NA))
print(is.equal(NaN,NA))
print(is.equal(NaN,Inf))
print(is.equal(Inf,Inf,Inf))
print(is.equal(NA,NA))
print(is.equal(NaN,NaN))
print(is.equal(1.4,1.4))
print(is.equal(1.4,1.400000001))
print(is.equal(1.4,1.400000001,1.41))
print(is.equal(0,0,0,0))

Jackknife Estimation

Description

Computes the "leave-one-out" Jackknife bias and standard error of an estimator f(x) of a data-vector x, or an estimator f(x,y) of vectors x and y. See Efron and Tibshirani (1993) for details

Usage

jackknife(x, f, y = NULL, ...)

Arguments

Value

Returns a list with the following components:

Author(s)

Danail Obreschkow

Compute the co-rotational kinetic energy fraction

Description

Computes the co-rotational kinetic energy fraction, \kappa_\mathrm{co}, of a particle distribution. This quantity is defined as the fraction of the total kinetic energy invested in ordered co-rotation about the total angular momentum axis. It is commonly used as a measure of the rotational support of galaxies in particle-based simulations.

Usage

kappacorotation(m, x, v, x0 = NULL, v0 = NULL)

Arguments

Details

The rotation axis is taken to be the direction of the total angular momentum of the particle distribution after subtracting the centre-of-mass position and velocity. Only particles with positive angular momentum about this axis contribute to the co-rotational kinetic energy.

Value

A scalar in the range [0,1] giving the fraction of the total kinetic energy contained in co-rotating azimuthal motion,

\kappa_\mathrm{co} = \frac{\sum_{j_z>0} m\,v_\phi^2} {\sum m\,|\mathbf{v}|^2},

where j_z is the specific angular momentum about the total angular momentum axis and v_\phi=j_z/R is the azimuthal velocity. The common factor of 1/2 in the kinetic energies cancels.

If m is NULL, NA is returned.

See Also

angularmomentum #'

Examples

m <- rep(1, 100)
x <- matrix(rnorm(300), ncol = 3)
v <- matrix(rnorm(300), ncol = 3)
kappacorotation(m, x, v)

Multi-dimensional adaptive kernel density estimation

Description

Produces a 2D kernel density estimation on a 2D grid from a D-dimensional (D>=2) point set

Usage

kde2(
  x,
  w = NULL,
  nx = 300,
  xlim = NULL,
  ylim = NULL,
  smoothing = 1,
  sigma = NULL,
  sigma.min = 0,
  sigma.max = Inf,
  reflect = "",
  algorithm = "kdenn",
  probability = FALSE
)

Arguments

Value

Returns a list of items

Author(s)

Danail Obreschkow

See Also

griddata

Examples

# make a mock sample of n d-dimensional points from
# three different components (1D line, 2D square, d-D normal distr)
d = 3 # number of dimensions of mock point set; try to choose different values 2, 3, 4, ...
n = 1e4 # number of particles per component
set.seed(1)
x = rbind(cbind(array(rep(runif(n,-1,1),2),c(n,2)),array(0,c(n,d-2))),
          cbind(array(runif(2*n),c(n,2)),array(0,c(n,d-2))),
          array(rnorm(d*n),c(n,d)))

# grid total projected probability density
npixels = 500 # number of pixels along a grid side
q = midseq(-3,3,npixels)
f1 = outer(dnorm(q),dnorm(q),'*')/3+outer(dunif(q),dunif(q),'*')/3
q = seq(round(npixels/3),round(npixels*2/3))
f1[q+npixels*(q-1)] = f1[q+npixels*(q-1)]+(npixels/6)^2/length(q)/3

# recover 2D projected pdf from 3D point sample using different methods
f2 = kde2(x, n=npixels, xlim=c(-3,3), ylim=c(-3,3), algorithm='basic', probability=TRUE)$field
f3 = kde2(x, n=npixels, xlim=c(-3,3), ylim=c(-3,3), algorithm='kdefast', probability=TRUE)$field
f4 = kde2(x, n=npixels, xlim=c(-3,3), ylim=c(-3,3), algorithm='kdenn', probability=TRUE)$field

# plot the 2D fields
img = function(f,x,y,title) {
  graphics::rasterImage(rasterflip(lim(f)^0.3),x,y,x+0.99,y+0.99)
  graphics::text(x+0.05,y+0.9,title,col='orange',pos=4)
}
oldpar = graphics::par(mar=rep(0.1,4))
nplot(c(0,2),c(0,2),asp=1)
img(f1,0,1,'Input pdf')
img(f2,1,1,'Random sample ("basic")')
img(f3,0,0,'Recovered pdf ("kdefast")')
img(f4,1,0,'Recovered pdf ("kdenn")')
graphics::par(oldpar)

Two-point correlation estimation

Description

Evaluates the Landy-Szalay (1993) estimator of the two-point correlation function of a point set D given a random comparison set R. The two point sets D and R can be made of different numbers of points, as the pair-counts are automatically normalized according to the number of points. In fact, it is often preferable to make the R set larger to reduce the R-related shot noise in the two-point estimator.

Usage

landyszalay(D, R, dr = 0.1, cpp = TRUE)

Arguments

Value

Returns a data frame with the two-point statistics of the data points:

Author(s)

Danail Obreschkow

See Also

paircount

Last element of a vector

Description

Returns the last element of a vector or the n-th element counting from the end of a vector.

Usage

last(x, n = 1)

Arguments

Value

scalar of the same type as x

Author(s)

Danail Obreschkow

Change lightness of a color

Description

Change lightness of a color

Usage

lightness(col, light = 0.5)

Arguments

Value

Returns a 9-character color or vector/array of 9-character colors.

Author(s)

Danail Obreschkow

See Also

transparent

Examples

# Generate different lightnesses of the same color
plot(runif(50),runif(50),pch=20,cex=10,col=lightness('purple',runif(50)))

Crop values to a custom range

Description

Limits the values of a vector or array to a desired interval, while keeping the shape of the input argument

Usage

lim(x, min = 0, max = 1, clip = NULL, na = NULL)

Arguments

Value

vector/array of the same shape as x

Author(s)

Danail Obreschkow

See Also

stretch

Linear Function Fitter

Description

Fits a global linear model y = a + b x and returns a function that predicts y for arbitrary x, similar in spirit to approxfun, but using a single least-squares line instead of piecewise interpolation.

Usage

linfun(x, y, na.rm = FALSE, ...)

Arguments

Details

This is a convenience wrapper around lm. It returns a callable function analogous to approxfun, but with a single global linear fit:

f(x) = a + b x,

where a and b are the intercept and slope from a least-squares regression of y on x.

Value

A function f(xnew) that evaluates the fitted linear regression at numeric values xnew.

See Also

[stats::lm()], [stats::approxfun()]

Examples

set.seed(1)
x = 1:10
y = 2 + 3 * x + rnorm(10)
f = linfun(x, y)
plot(x,y)
curve(f, col='red', add = TRUE) # show linear fit
points(6.6,f(6.6),col='red') # show predicted y-value at x = 6.6

Handle spaces in Linux filenames

Description

Convert spaces in filenames (" ") to linux-type spaces "\ ", needed when calling system() on macOS.

Usage

linuxspaces(txt)

Arguments

Value

filename with modified spaces, e.g. "my\ file\ 1.txt"

Author(s)

Danail Obreschkow

Examples

filename = '~/Desktop/my file 1.txt'
command = sprintf('ls -l %s',linuxspaces(filename))
## Not run: 
system(command)

## End(Not run)

Read binary data into array

Description

Reads binary data using the base function readBin and recasts it into an array of custom dimensions.

Usage

loadbin(
  filename,
  dim,
  bytes = 4,
  type = "numeric",
  signed = FALSE,
  endian = "little"
)

Arguments

Value

Returns an array of dimension dim.

Author(s)

Danail Obreschkow

Display a single movie frame

Description

Displays a single movie-frame in the R-console, exactly as used in a movie generated with makemovie.

Usage

makeframe(
  frame.draw,
  frame.index,
  width = 1080,
  height = 720,
  cex = 1,
  oversampling = 1,
  pngfile = NULL
)

Arguments

Value

Returns the displayed image as n-by-my-by 4 array, representing the 4 RGBA channels of height n and width m.

Author(s)

Danail Obreschkow

See Also

makemovie

Examples

## Example: Movie of a manual clock

# Function to draw a single clock face with two hands
frame = function(time) {
  oldpar = graphics::par(mar=c(0,0,0,0))
  nplot(xlim=c(-1.1,1.1),ylim=c(-1.1,1.1),pty='s')
  plotrix::draw.circle(0,0,1,col='#aaaaff')
  radius = c(0.5,0.9)
  speed = 2*pi/c(720,60)
  lwd = c(4,2)
  graphics::arrows(0,0,radius*sin(speed*time),radius*cos(speed*time),lwd=lwd)
  graphics::par(oldpar)
}

# Produce movie
## Not run: 
makeframe(frame,15,200,200)

## End(Not run)

Produce a movie from frame-drawing function

Description

Generates an MP4-movie provided a custom function that plots individual frames. The routine has been developed and tested for MacOS and it requires on a working installation of ffmpeg.

Usage

makemovie(
  frame.draw,
  frame.index,
  output.path,
  output.filename,
  width = 1080,
  height = 720,
  fps = 60,
  keep.frames = FALSE,
  quiet = FALSE,
  separator = "/",
  ffmpeg.cmd = "ffmpeg",
  ffmpeg.opt = "-vcodec libx264 -crf 18 -pix_fmt yuv420p",
  manual = FALSE,
  cex = 1,
  oversampling = 1,
  first.index = 1,
  last.index = length(frame.index)
)

Arguments

Value

Linux command line to convert frames into movie using ffmpeg.

Author(s)

Danail Obreschkow

See Also

makeframe

Examples

## Example: Movie of a manual clock

# Function to draw a single clock face with two hands
frame = function(time) {
  oldpar = graphics::par(mar=c(0,0,0,0))
  nplot(xlim=c(-1.1,1.1),ylim=c(-1.1,1.1),pty='s')
  plotrix::draw.circle(0,0,1,col='#aaaaff')
  radius = c(0.5,0.9)
  speed = 2*pi/c(720,60)
  lwd = c(4,2)
  graphics::arrows(0,0,radius*sin(speed*time),radius*cos(speed*time),lwd=lwd)
  graphics::par(oldpar)
}

# Produce movie
## Not run: 
makemovie(frame,seq(0,60,0.5),'~/testmovie','movie.mp4',200,200)

## End(Not run)

Monte Carlo and Quasi-Monte Carlo integration in any dimension

Description

Numerical integration using a Monte Carlo (MC) or Quasi-Monte Carlo (QMC) algorithm, based on a Halton sequence. These algorithms are of low order (1/sqrt(n) for MC, log(n)/n for QMC in one dimension) compared to the typical orders of 1D deterministic integrators, such as those available in the integrate function. The MC and QMC integrators are suitable to compute D-dimensional integrals with D>>1, since the order of most deterministic methods deteriorates exponentially with D, whereas the order of MC remains 1/sqrt(n), irrespective of D, and the order of QMC only deteriorates slowly with D as log(n)^D/n.

Usage

mcintegral(f, a, b, n = 1e+05, qmc = FALSE, seed = NULL, warn = TRUE)

Arguments

Value

Returns a list of items:

Author(s)

Danail Obreschkow

Examples

## Numerically integrate sin(x)
f = function(x) sin(x)
m = mcintegral(f,0,pi)
cat(sprintf('Integral = %.3f\u00B1%.3f (true value = 2)\n',m$value,m$sigma))

## Numerically compute the volume of a unit sphere
sphere = function(x) as.numeric(rowSums(x^2)<=1) # this is vectorized
vmc = mcintegral(sphere,rep(-1,3),rep(1,3),seed=1)
vqmc = mcintegral(sphere,rep(-1,3),rep(1,3),qmc=TRUE)
cat(sprintf('Volume of unit sphere = %.3f\u00B1%.3f (MC)\n',vmc$value,vmc$error))
cat(sprintf('Volume of unit sphere = %.3f\u00B1%.3f (QMC)\n',vqmc$value,vqmc$error))
cat(sprintf('Volume of unit sphere = %.3f (exact)\n',4*pi/3))

Mid-points of regular grid

Description

Compute the mid-point positions of a one-dimensional regular grid of n equal intervals.

Usage

midseq(min, max, n = 1)

Arguments

Value

vector of mid points

Author(s)

Danail Obreschkow

Mollweide projection

Description

Performs a Mollweide projection (also known as Babinet projection, homalographic projection, homolographic projection, and elliptical projection) of longitude and latitude coordinates. The most important feature of the Mollweide projection is that it preserves surface areas, which makes it a commonly used projection in geography, astronomy and cosmology. The total surface area of the standard projection is equal to the surface area of the unit sphere (4pi); and the shape of the fully projected sphere is an ellipse (with axes lengths 2*sqrt(2) and sqrt(2)).

Usage

mollweide(lon, lat, lon0 = 0, radius = 1, deg = FALSE)

Arguments

Value

Returns an n-by-2 matrix of 2D Cartesian coordinates x and y.

Author(s)

Danail Obreschkow

Examples

lon = runif(1e4,0,2*pi)
lat = asin(runif(1e4,-1,1)) # = uniform sampling of the sphere
plot(mollweide(lon,lat),xlim=c(-3,3),ylim=c(-1.5,1.5),pch=16,cex=0.5)
plotrix::draw.ellipse(0,0,2*sqrt(2),sqrt(2),border='orange',lwd=2)

Second moment tensor

Description

Computes the second moment tensor of a set of point masses in arbitrary dimensions,

\mathbf{M} = \sum_i m_i\,\mathbf{x}_i\mathbf{x}_i^\mathrm{T}.

Usage

moments(x, m = 1)

Arguments

Value

A symmetric d\times d matrix containing the second moment tensor, where d is the dimension of the points.

See Also

quadrupole, inertia

Mutual information of two random variables

Description

Computes the mutual information of two random variables X and Y, given their 2D density represented in a matrix.

Usage

mutual(x, y = NULL, b = exp(1), n = NULL, xlim = NULL, ylim = NULL)

Arguments

Value

Returns a list of items:

Author(s)

Danail Obreschkow

Non-uniform Discrete Fourier Transform

Description

Compute the one-dimensional Non-uniform Discrete Fourier Transform (NDFT). This is needed if the data are sampled at irregularly spaced times.

Usage

ndft(
  f,
  x = seq(0, length(f) - 1)/length(f),
  nu = seq(0, length(f) - 1),
  inverse = FALSE,
  weighing = TRUE,
  simplify = TRUE
)

Arguments

Details

The one-dimensional NDFT of a vector f=(f_1,...,f_N) is defined as

F_j=\sum_i w_i f_i exp(-2\pi i*x_i*\nu_j)

where w_i are optional weights, proportional to the interval around x_i, only used if weighing=TRUE. Likewise, the inverse NDFT is defined as

f_i=\sum_j w_j F_j exp(+2\pi i*x_i*nu_j)

where w_j are optional weights, proportional to the interval around \nu_j. In this implementation NDFTs are computed using a brute force algorithm, scaling as O(N*N), which is considerably worse than the O(N)*log(N) scaling of FFT algorithms. It is therefore important to pick the required frequencies wisely to minimise computing times.

Value

Returns a vector of the same length as x (if inverse=FALSE) or nu (if inverse=TRUE).

Author(s)

Danail Obreschkow

See Also

fft

Examples

# Define an example signal
nu1 = 1 # [Hz] first frequency
nu2 = 8 # [Hz] second frequency in the signal
s = function(t) sin(2*pi*nu1*t)+0.7*cos(2*pi*nu2*t+5)

# Discretize signal
N = 50 # number of samples
t.uniform = seq(0,N-1)/N
t.nonuniform = t.uniform^1.3
s.uniform = s(t.uniform)
s.nonuniform = s(t.nonuniform)

# Plot signal
oldpar = par(mfrow = c(1, 2))
curve(s,0,1,500,xaxs='i',main='Time signal',xlab='Time t',ylab='s(t)',col='grey')
points(t.uniform,s.uniform,pch=16,cex=0.8)
points(t.nonuniform,s.nonuniform,pch=4,col='blue')
legend('topright',c('Continuous signal','Uniform sample','Non-uniform sample'),
       lwd=c(1,NA,NA),pch=c(NA,16,4),col=c('grey','black','blue'),pt.cex=c(1,0.8,1))

# Uniform and non-uniform DFT
nu = seq(0,N-1) # discrete frequencies
spectrum.uniform = stats::fft(s.uniform)
spectrum.nonuniform = ndft(s.nonuniform,t.nonuniform,nu)
spectrum.wrong = stats::fft(s.nonuniform)

# Evaluate power
power.uniform = Mod(spectrum.uniform)^2
power.nonuniform = Mod(spectrum.nonuniform)^2
power.wrong = Mod(spectrum.wrong)^2

# Plot DFT and NDFT up to Nyquist frequency
plot(nu,power.uniform,pch=16,cex=0.8,xlim=c(0,N/2),xaxs='i',
               main='Power spectrum',xlab=expression('Frequency'~nu~'[Hz]'),ylab='Power')
points(nu,power.nonuniform,pch=4,col='blue')
points(nu,power.wrong,pch=1,col='red')
abline(v=c(nu1,nu2),col='grey',lty=2)
legend('topright',c('DFT of uniform sample','NDFT of non-uniform sample',
'DFT of non-uniform sample (wrong)','Input frequencies'),
       lwd=c(NA,NA,NA,1),lty=c(NA,NA,NA,2),pch=c(16,4,1,NA),
       col=c('black','blue','red','grey'),pt.cex=c(0.8,1,1,NA))
par(oldpar)

Draw a regular n-gon

Description

Draws a regular polygon with n sides, such as a triangle (n=3) or hexagon (n=6).

Usage

ngon(x = 0, y = 0, s = 1, n = 6, angle = 0, fix.aspect = TRUE, ...)

Arguments

Value

None.

Author(s)

Danail Obreschkow

Examples

## Plot random points on the unit sphere in Mollweide projection
# hexagon at the center of a plot
nplot(bty='o', asp=0.5)
ngon(x=0.5, y=0.5, s=0.1, n=6, fix.aspect=FALSE)

ngon(x=0.5, y=0.5, s=0.1, n=6, border='red')

plot(NA,xlim=c(1,1e3),ylim=c(1,1e5),log='xy')
ngon(x=10^runif(10,0,3), y=10^runif(10,0,5), s=1, n=6,  border='red',lwd=3)

Make empty plot area

Description

Open an empty plot

Usage

nplot(
  xlim = c(0, 1),
  ylim = c(0, 1),
  xlab = "",
  ylab = "",
  xaxs = "i",
  yaxs = "i",
  xaxt = "n",
  yaxt = "n",
  bty = "n",
  ...
)

Arguments

Value

None

Author(s)

Danail Obreschkow

Count the number of point-pairs in distance bins

Description

Count the number of point-pairs in equally spaced distances bins. Code works in any dimension. If only one point set is provided, the distances of this point set with itself are used (counting each pairs only once, i.e. only ij, not ji).

Usage

paircount(x, y = NULL, dr, rmax, cpp = TRUE)

Arguments

Value

Returns a data frame of two column vectors:

Author(s)

Danail Obreschkow

See Also

landyszalay

Convert pdf to jpg

Description

Calls the console "convert" function to convert a pdf-file into a jpeg-image. Requires "convert" to be installed already.

Usage

pdf2jpg(
  pdf.filename,
  jpg.filename,
  quality = 100,
  background = "white",
  dim = c(1600, 1200),
  remove.pdf = FALSE,
  verbose = TRUE
)

Arguments

Value

None

Author(s)

Danail Obreschkow

Planck CMB colour palette

Description

Generates a colour scale that closely mimics the official colours used to display the Planck CMB temperature map from -300uK to +300uK.

Usage

planckcolors(n)

Arguments

Value

Returns an n-vector of RGB colour strings.

Author(s)

Danail Obreschkow

Examples

nplot()
rasterImage(rbind(planckcolors(1e3)),0,0,1,1)

Polar/cylindrical to Cartesian coordinate conversion

Description

Convert polar/cylindrical coordinates to 2D/3D Cartesian coordinates

Usage

pol2car(x)

Arguments

Value

Returns a 2/3-element vector or a n-by-2/3 element matrix representing the Cartesian coordinates (x,y)/(x,y,z).

Author(s)

Danail Obreschkow

See Also

car2pol

Show progress while timer in running

Description

Show progress while a timer is running following a call of tick(). To suppress all progress updates, set the environment variable .cooltools.env$suppressProgress to TRUE.

Usage

progress(txt = NULL)

Arguments

Value

None

Author(s)

Danail Obreschkow

See Also

tick tock error

Examples

tick('Test')
Sys.sleep(.1)
for (i in seq(3)) {
  progress(i)
  Sys.sleep(.1)
}
tock()

Quadrupole tensor

Description

Computes the trace-free mass quadrupole tensor of a set of point masses in arbitrary dimensions,

\mathbf{Q} = \sum_i m_i \left( d\,\mathbf{x}_i\mathbf{x}_i^\mathrm{T} - |\mathbf{x}_i|^2\mathbf{I} \right),

where d is the dimensionality of the points and \mathbf{I} is the d\times d identity matrix.

Usage

quadrupole(x, m = 1)

Arguments

Value

A symmetric d\times d matrix containing the trace-free quadrupole tensor.

See Also

moments, inertia

Suppress console output

Description

Executes routines and command lines and load packages while suppressing the console output.

Usage

quiet(x)

Arguments

Value

Returns whatever the called routine returns in invisible form.

Author(s)

Danail Obreschkow

Examples

# Test function
test = function(x) {
  cat('This routine is likes to talk a lot!\n')
  return(x^2)
}

# Standard call call:
y = test(5)
print(y)

# Quiet call:
y = quiet(test(6))
print(y)

Flip array to be displayed with rasterImage()

Description

Flips the array A to be displayed with rasterImage, such that the first index runs from left to right and second index runs from bottom to top, like in standard Cartesian coordinates. In this way rasterImage(rasterflip(A)) has the same orientation as image(A).

Usage

rasterflip(A)

Arguments

Value

A m-by-n array of a monochromatic image or m-by-n-by-k array of a color image.

Author(s)

Danail Obreschkow

Read data from an HDF5 file

Description

Reads data from an HDF5 file, including attributes, into a nested list that preserves the hierarchical structure of the HDF5 data. The routine can read 64-bit integers and represent them as 64-bit integers in R via the bit64 package. It can also read 32-bit floating-points values, not part of native R. These are automatically converted into 64-bit doubles.

Usage

readhdf5(file, subtree = "*", group.attr.as.data = FALSE, empty = FALSE)

Arguments

Details

This function, based on the hdf5r package, recursively parses and reads HDF5 files into nested lists that preserve the original hierarchy. Attributes in groups and datasets are included in the output.

Value

Nested list representing the contents of the HDF5 file. Groups are nested sublists, datasets are represented by their data. Attributes of groups and datasets are attached as attributes to the corresponding sublists and data elements.

See Also

writehdf5 for examples.

Re-bin density histograms

Description

Transform density histogram data into histogram data with different bins

Usage

rebindensity(x, y, xout)

Arguments

Value

m-vector of y-values associated with the bins specified by xout.

Author(s)

Danail Obreschkow

Examples

# original binning
x = seq(0.5,4.5)
y = seq(5)
plot(x,y,xlim=c(-1,6),ylim=c(0,6),pch=16)
lines(histcoord(x,y),lwd=3)

# rebinning
xout = seq(-0.9,6,0.3)
yout = rebindensity(x,y,xout)
points(xout,yout,col='red',pch=16)
lines(histcoord(xout,yout),col='red')

Random number generator for a custom d-dimensional distribution

Description

Brute-force algorithm for drawing random numbers from a d-dimensional distribution.

Usage

rng(f, n, min, max, fmax = NULL, quasi = FALSE, start = 1, warn = TRUE)

Arguments

Value

Returns list of items:

Author(s)

Danail Obreschkow

See Also

dpqr

Examples

## 1D random number generation from a sine-function
f = function(x) sin(x)
out.pseudo = rng(f,1e3,0,pi)
out.quasi = rng(f,1e3,0,pi,quasi=TRUE)
hist(out.pseudo$x,100,freq=FALSE,border=NA,xlab='x',main='sine-distribution')
hist(out.quasi$x,100,freq=FALSE,border=NA,col='#ff000066',add=TRUE)
curve(sin(x)/2,0,pi,add=TRUE)

## 2D quasi-random sampling of a disk with exponentially declining surface density
f = function(x) exp(-sqrt(x[,1]^2+x[,2]^2))
out = rng(f,1e4,min=c(-5,-5),max=c(5,5),quasi=TRUE)
plot(out$x,cex=0.3,pch=16,asp=1,main='Quasi-random exponential disk')

## 5D random number generation (5-dimensional sphere)
f = function(x) as.numeric(sum(x^2)<=1)
out = rng(f,1e4,rep(-1,5),rep(1,5))
cat(sprintf('Number of successes over number of trials : %.4f\n',out$n/out$ntrials))
cat(sprintf('Expected ratio for n=\u221E : %.4f\n',pi^(5/2)/gamma(1+5/2)/2^5))

2D rotation matrix

Description

Compute a 2D rotation matrix given a rotation angle

Usage

rotation2(angle)

Arguments

Value

Returns a 2-by-2 anti-symmetric rotation matrix

Author(s)

Danail Obreschkow

See Also

rotation3

3D rotation matrix

Description

Compute a 3D rotation matrix given an axis and an angle

Usage

rotation3(u, angle = NULL)

Arguments

Value

Returns a 3-by-3 rotation matrix

Author(s)

Danail Obreschkow

See Also

rotation2

Generate randomly oriented vectors in 2D

Description

Generate randomly oriented vectors in 2D, following an isotropic distribution (optionally truncated to a region).

Usage

runif2(n = 1, r = c(0, 1), azimuth = c(0, 2 * pi), quasi = FALSE, start = 1)

Arguments

Value

Returns an n-by-2 array of n vectors.

Author(s)

Danail Obreschkow

See Also

runif3

Examples

## generate 500 unit vectors with radii between 0.5 and 1
x = runif2(500,r=c(0.5,1))
oldpar = par(pty='s')
plot(x,pch=20)
par(oldpar)

Generate randomly oriented vectors in 3D

Description

Generate randomly oriented vectors in 3D, following an isotropic distribution (optionally truncated to a region).

Usage

runif3(
  n = 1,
  r = c(0, 1),
  azimuth = c(0, 2 * pi),
  polarangle = c(0, pi),
  quasi = FALSE,
  start = 1
)

Arguments

Value

Returns an n-by-3 array of n vectors.

Author(s)

Danail Obreschkow

See Also

runif2

Examples

## draw 20 unit vectors on a sphere
x = runif3(20,r=c(1,1))
print(rowSums(x^2))

Scalar product

Description

Computes the scalar (dot) product of vectors. The inputs may be individual vectors or matrices whose rows represent vectors.

Usage

scalarproduct(x, y)

Arguments

Value

Returns either a scalar (vector-vector case) or a vector containing the row-wise scalar products.

See Also

scalarproduct, vectornorm, unitvector

Sigmoid function

Description

Returns S-shaped functions, monotonically increasing between two constants with a point symmetry at the half-way point.

Usage

sigmoid(x, type = 1, fmin = 0, fmax = 1, xmid = 0, slope = 1)

Arguments

Value

Returns a vector with the same number of elements as x

Author(s)

Danail Obreschkow

Examples

col = c('red','orange','#cccc00','#00bb00','blue','purple')
for (i in seq(6)) {
  graphics::curve(sigmoid(x,i),-3,3,add=i>1,col=col[i],ylab='sigmoid(x)')
}
graphics::lines(c(-1.5,1.5),c(-1,2),lty=2)

Round a vector of floating-point values while preserving their sum

Description

Rounds the values in a vector up and down, preserving the sum of the vector and minimizing the total rounding error under this condition. An example where this is useful is when rounding a vector of percentages, where the total should add up to 100 percent.

Usage

smartround(x, digits = 0)

Arguments

Value

Returns a vector of rounded values of the same length as x.

Author(s)

Danail Obreschkow

Examples

x = runif(5)
x = x/sum(x)*100
print(x)
print(sum(x))
y = smartround(x)
print(y)
print(sum(y))
y2 = smartround(x,2)
print(y2)
print(sum(y2))

Draw smoothed contours

Description

Draw smoothed iso-countours for a density field. The contours are computed using the contourLines routine and smoothed using the smooth.spline function. Both open and closed contour lines are handled correctly.

Usage

smoothcontour(
  x = seq(0, 1, length.out = nrow(z)),
  y = seq(0, 1, length.out = ncol(z)),
  z,
  levels,
  smoothing = 0.5,
  min.radius = 1,
  lwd = 1,
  lty = 1,
  col = "black",
  ...
)

Arguments

Value

None

Author(s)

Danail Obreschkow

See Also

contourLines, smooth.spline

Examples

set.seed(1)
f = function(x) cos(2*x[1]-x[2]-1)^2*exp(-x[1]^2-x[2]^2-x[1]*x[2])
x = seq(-3,3,length=100)
m = pracma::meshgrid(x)
z = array(Vectorize(function(x,y) f(c(x,y)))(m$Y,m$X)+rnorm(4e4,sd=0.1),dim(m$X))
image(x,x,z,col=terrain.colors(100))
contour(x,x,z,levels=c(0.2,0.5),add=TRUE)
smoothcontour(x,x,z,levels=c(0.2,0.5),lwd=3,smoothing=0.8,min.radius=2)

Smoothed Function

Description

Generates a cubic smoothed spline function y=f(x) approximating supplied (x,y)-data with a custom number of degrees of freedom. The routines builds on smooth.spline, but directly returns the smoothed function rather than the fitted model.

Usage

smoothfun(x, y = NULL, w = NULL, df = NULL, ...)

Arguments

Value

Returns a fast and vectorized smoothed function f(x).

Author(s)

Danail Obreschkow

See Also

smooth.spline

Examples

# make random data set
set.seed(1)
x = runif(100)
y = sin(2*pi*x)+rnorm(100, sd=0.5)
plot(x,y,pch=16)

# smoothed spline
f = smoothfun(x, y)
curve(f, add=TRUE, col='red')

# smoothed spline with custom degree of freedom
g = smoothfun(x, y, df=5)
curve(g, add=TRUE, col='blue')

Sort List Entries Alphabetically

Description

This function sorts the entries of a structured list in alphabetical order, including all nested sublists recursively, making it easier to navigate nested list structures and compare different lists. Unnammed items are listed last. Attributes of the original list are preserved.

Usage

sortlist(lst, recursive = TRUE, convert.data.frames = FALSE)

Arguments

Value

Structured list with alphabetically sorted entries. If recursive = TRUE, all sublists are also sorted alphabetically. If the list contains unnamed elements, they are left unsorted.

Examples

# Create a nested list
my_list = list(
b = list(d = 4, a = 2),
empty = list(),
'unnamed',
a = 1,
f = list(b = 3, a = list(d = seq(5), 'unnamed', c = 2))
)

# Sort the list recursively
sorted_list_recursive = sortlist(my_list)
print(sorted_list_recursive)

Spectrum colour palette

Description

Generates smooth rainbow color scale from red to purple, similar to rainbow, but with improved smoothness

Usage

spectrumcolors(n, alpha = 1, rev = FALSE)

Arguments

Value

Returns an n-vector of RGB colour strings.

Author(s)

Danail Obreschkow

#' @examples nplot() rasterImage(rbind(spectrumcolors(1e3)),0,0,1,0.5) rasterImage(rbind(rainbow(1e3,end=5/6)),0,0.5,1,1) text(0.5,0.25,'spectrum') text(0.5,0.75,'rainbow') abline(h=0.5)

Spherical to Cartesian coordinate conversion

Description

Convert 3D spherical to Cartesian coordinates

Usage

sph2car(x)

Arguments

Value

Returns a 3-element vector or a n-by-3 element matrix representing the Cartesian coordinates (x,y,z)

Author(s)

Danail Obreschkow

See Also

car2sph

Plot a spherical function or point set

Description

Plots a spherical function or a point set in a 2D projection using only standard R graphics. This avoids compatibility issues of rgl, e.g. knitting markdown documents.

Usage

sphereplot(
  f = NULL,
  n = 100,
  theta0 = pi/2,
  phi0 = 0,
  angle = 0,
  projection = "globe",
  col = gray.colors(256, 0, 1),
  clim = NULL,
  add = FALSE,
  center = c(0, 0),
  radius = 1,
  nv = 500,
  show.border = TRUE,
  show.grid = TRUE,
  grid.phi = seq(0, 330, 30)/180 * pi,
  grid.theta = seq(30, 150, 30)/180 * pi,
  pch = 16,
  pt.col = "black",
  pt.cex = 0.5,
  lwd = 0.5,
  lty = 1,
  line.col = "black",
  background = "white",
  ...
)

Arguments

Value

Returns a list containing the vector col of colors and 2-vector clim of values corresponding to the first and last color.

Author(s)

Danail Obreschkow

Examples

## Plot random points on the unit sphere in Mollweide projection
set.seed(1)
f = cbind(acos(runif(5000,-1,1)),runif(5000,0,2*pi))
sphereplot(f,theta0=pi/3,projection='mollweide',pt.col='red')

## Plot real spherical harmonics up to third degree
oldpar = par(mar=c(0,0,0,0))
nplot(xlim=c(-4,3.5),ylim=c(0,4),asp=1)
for (l in seq(0,3)) { # degree of spherical harmonic
  for (m in seq(-l,l)) { # order of spherical harmonic

    # make spherical harmonic function in real-valued convention
    f = function(theta,phi) sphericalharmonics(l,m,cbind(theta,phi))

    # plot spherical harmonic
    sphereplot(f, 50, col=planckcolors(100), phi0=0.1, theta0=pi/3, add=TRUE, clim=c(-0.7,0.7),
               center=c(m,3.5-l), radius=0.47)

    if (l==3) text(m,-0.15,sprintf('m=%+d',m))
  }
  text(-l-0.55,3.5-l,sprintf('l=%d',l),pos=2)
}
par(oldpar)

Spherical Harmonics

Description

Evaluates spherical harmonics Y, either in the real-valued or complex-valued basis.

Usage

sphericalharmonics(l, m, x, basis = "real")

Arguments

Value

Returns an n-vector of the spherical harmonics; for points x=c(0,0,0), a value of 0 is returned

Author(s)

Danail Obreschkow

Examples

## Check orthonormalization of all spherical harmonics up to 3rd degree

# make indices l and m up to 3rd degree
l = c(0,rep(1,3),rep(2,5),rep(3,7))
m = c(0,seq(-1,1),seq(-2,2),seq(-3,3))

# check orthonormalization for all pairs
for (i in seq(16)) {
  for (j in seq(16)) {

    # compute scalar product
    f = function(theta,phi) {
      Yi = sphericalharmonics(l[i],m[i],cbind(theta,phi))
      Yj = sphericalharmonics(l[j],m[j],cbind(theta,phi))
      return(Re(Yi*Conj(Yj))*sin(theta))
    }
    g = Vectorize(function(phi) integrate(f,0,pi,phi)$value)
    scalar.product = integrate(g,0,2*pi)$value

    # compare scalar product to expected value
    ok = abs(scalar.product-(i==j))<1e-6
    cat(sprintf('(l=%1d,m=%+1d|l=%1d,m=%+1d)=%5.3f  %s\n',l[i],m[i],l[j],m[j],
                scalar.product+1e-10,ifelse(ok,'ok','wrong')))
  }
}

Draw streamlines of a two-dimensional vector field

Description

Draws streamlines of a planar vector field f(x,y) by integrating forward and backward from a set of starting points. Starting points can be provided as a single integer, a length-2 integer vector, or an n \times 2 matrix of coordinates.

Usage

streamlines(
  f,
  xlim,
  ylim,
  points = 20,
  nsteps = 100,
  arrows = FALSE,
  add = FALSE,
  ...
)

Arguments

Details

The integration uses a fixed step length based on the plot diagonal,

d = \sqrt{dx^2 + dy^2} / nsteps,

where 'dx = diff(range(xlim))' and 'dy = diff(range(ylim))'. At each step, only the direction of the vector field is used, via [unitvector()], so the result shows streamlines of the field rather than trajectories with speed information.

If 'add = FALSE', a new plot is initialized using [plot()].

Arrows are drawn along the local field direction. If 'arrows = TRUE', each arrow is centered on the streamline starting point. If 'arrows' is a matrix of coordinates, the arrows are placed at those coordinates, regardless of whether they lie exactly on a streamline.

Value

Invisibly returns a list with components:

• 'forward': array of forward-integrated coordinates;

• 'backward': array of backward-integrated coordinates;

• 'start': matrix of starting coordinates.

Examples

f = function(x,y) {
  vx = -x-y
  vy = -y+(x-0.5)
  cbind(vx,vy)
}

streamlines(f, c(-1,1), c(-1,1), 50, 500, col = "#5555ff")
streamlines(f, c(-1,1), c(-1,1), 20, 200, arrows = TRUE, add = FALSE,
            col = "#5555ff")

Stretch values to a custom range

Description

Shifts and stretches the values of a vector or array to a desired interval, while maintaining the shape of the input argument

Usage

stretch(x, min = 0, max = 1, invert = FALSE, gamma = NULL, na = NULL)

Arguments

Value

vector/array of the same shape as x

Author(s)

Danail Obreschkow

See Also

lim

Insert a sub-panel into plot

Description

Insert a sub-panel into an existing plotting area. To open a subplot, call par(subplot(...)), to close it, you must call par(subplot('off')).

Usage

subplot(xleft = NA, ybottom, xright, ytop)

Arguments

Value

graphical parameters to user with par.

Author(s)

Danail Obreschkow

Examples

# main plot
f = function(x) x*sin(1/(x+.Machine$double.eps))
curve(f,0,1,n=1000,ylim=c(-1,1),xlab='',ylab='',xaxs='i',yaxs='i')
rect(0.02,-0.1,0.1,0.1,border='blue',col=transparent('blue',0.2))

# subplot
par(subplot(0.55,-0.8,0.93,0))
curve(f,0.02,0.1,n=500,ylim=c(-0.1,0.1),xlab='',ylab='',xaxs='i',yaxs='i')
rect(0.02,-0.1,0.1,0.1,border=NA,col= transparent('blue',0.2))
par(subplot('off'))

Start timer

Description

Start timer and write a custom text into the console.

Usage

tick(txt = "Start")

Arguments

Value

None

Author(s)

Danail Obreschkow

See Also

tock progress error

Examples

tick('Sum 10 million random numbers')
x = sum(runif(1e7))
tock()

Stop timer

Description

Stop timer and write the computation in seconds since the last call of tick().

Usage

tock(txt = "", fmt = " (%.2fs). %s\n")

Arguments

Value

Returns the elapsed time in seconds since calling tick().

Author(s)

Danail Obreschkow

See Also

tick progress error

Examples

tick('Sum 10 million random numbers')
x = sum(runif(1e7))
tock()

Add transparency to a color

Description

Add transparency to a color

Usage

transparent(col, alpha = 0.5)

Arguments

Value

Returns a 9-character color or vector/array of 9-character colors.

Author(s)

Danail Obreschkow

See Also

lightness

Examples

# Add different transparencies of the same color
plot(runif(50),runif(50),pch=20,cex=10,col=transparent('purple',runif(50)))

# Add the same transparency to different colors
plot(runif(50),runif(50),pch=20,cex=10,col=transparent(rainbow(50)))

Zoom, translate and rotate array image

Description

Zoom/rotate/translate an image relative to its center for images represented as simple arrays.

Usage

transzoom(
  img = NULL,
  zoom = 1,
  shift = c(0, 0),
  angle = 0,
  size = NULL,
  col = "black",
  filter = "bilinear",
  file.in = "",
  file.out = "",
  format = "png",
  show.image = TRUE
)

Arguments

Value

Returns an n-by-m-by-3 array or n-by-m-by-4 array of the processed image.

Author(s)

Danail Obreschkow

Examples

img = yinyangyong # this is an example image included in the package
transzoom(img, zoom=2) # zoom by a factor 2

Map positive 64-bit integers onto unique doubles

Description

Uses a reversible one-to-one mapping to converts non-negative 64-bit integers, produced via the bit64 package, to double-precision floating-point values for faster handling of 64-bit integer indices, e.g. comparison of indices, count number of unique indices, etc. The double values are completely different from the corresponding integer values.

Usage

uniquedouble(x, inverse = FALSE)

Arguments

Value

Scalar/vector/array of double values, if inverse=FALSE, or integer64 values, if inverse=TRUE.

Author(s)

Danail Obreschkow

Examples

# produce example vector of non-negative 64-bit integers
n = 1e3
x = c(bit64::as.integer64(0:n),
      bit64::as.integer64('9218868437227405311')+c(-n:n),
      bit64::lim.integer64()[2]-c(0:n))

# map this vector onto unique doubles
y = uniquedouble(x)

# check if the double-vector can be inverted back to the integer-vector
cat(sprintf("Check inversion = %s.\n", all(x==uniquedouble(y, inverse=TRUE))))

# measure the time taken to compare all the int64 values
tick("Compare original integer64 values. ")
comparison = TRUE
for(i in seq_along(x)) comparison = comparison & x[i]==x[i]
tock(comparison, fmt="Time taken = %.3fs. Test result = %s.\n")

# measure the time taken to compare all the corresponding double values
tick("Compare corresponding double values. ")
comparison = TRUE
for(i in seq_along(y)) comparison = comparison & y[i]==y[i]
tock(comparison, fmt="Time taken = %.3fs. Test result = %s.\n")

Normalize vectors to unit length

Description

Computes unit vectors from one or more vectors.

Usage

unitvector(x)

Arguments

Value

An object of the same type and dimensions as x, with each vector normalized to unit length. Zero vectors are returned unchanged.

See Also

vectornorm

Retrieve Custom Object Attributes

Description

This function takes a generic R object and returns the custom attributes added by the user, distinguishing them from the object's inherent attributes.

Usage

userattributes(obj)

Arguments

Value

Named list of custom attributes, mimicking the output of the base function attributes, but only including custom attributes

See Also

attributes

Examples

# Create a list containing objects of different classes
example = list(1,
               as.integer(1:3),
               bit64::as.integer64(1),
               diag(5),
               array(0,c(2,3,5)),
               list(a=1,b='c'),
               data.frame(x=1, y=1:10, char=LETTERS[1:10])
)

# Loop over all the objects to check if the custom attribute is correctly extracted.
for (i in seq_along(example)) {
  x = example[[i]]
  attr(x,'MyAttribute') = 'Test'
    cat(sprintf('Class = %s, Custom attribute(s) = %s\n',class(x)[1],
                paste(names(userattributes(x)), collapse = ', ')))
}

Vector norm

Description

Compute the norm of a vector or a matrix of vectors

Usage

vectornorm(x)

Arguments

Value

Returns a scalar or an n-element vector with the vector norms

Author(s)

Danail Obreschkow

See Also

vectornorm

Vector product

Description

Computes the vector (cross) product of two three-dimensional vectors. The inputs may be individual vectors or matrices whose rows represent vectors.

Usage

vectorproduct(x, y, normalize = FALSE)

Arguments

Value

Returns either a 3-element vector or an n\times3 matrix containing the row-wise vector products.

See Also

scalarproduct, vectornorm, unitvector

Convert wavelength to RGB

Description

Converts a given wavelength of light to an approximate RGB color value, using black in the invisible range.

Usage

wavelength2col(wavelength)

Arguments

Value

Returns a color string or vector of color strings with the same number of elements as wavelength.

Author(s)

Danail Obreschkow

Source

Smoothed implementation of the original Fortran version by Dan Bruton (http://www.physics.sfasu.edu/astro/color/spectra.html) and the R-function by Michael Friendly (https://gist.github.com/friendly).

Examples

lambda = seq(300,800)
col = matrix(wavelength2col(lambda),nrow=1)
plot(NA,xlim=range(lambda),ylim=c(0,1),xaxs='i',xlab='wavelength [nm]',yaxs='i',yaxt='n',ylab='')
rasterImage(col,min(lambda),0,max(lambda),1)

Write structured list to an HDF5 file

Description

Writes a named list, including its hierarchy of nested sublists, to an HDF5 file with HDF5 groups, subgroups and datasets preserving the hierarchical structure. The routine supports standard HDF5 data types, including double-precision floating point numbers (64 bit), integers (32 bit), characters (8 bit), and booleans, as well as vectors and arrays of these types. It also supports 64-bit integers (H5T_NATIVE_INT64), available via the bit64 package, as well as 32-bit floating-point objects of class float32. Custom attributes of sublists and data, assigned in R via attr, are automatically transcribed to group and dataset attributes in the HDF5 file. Such attributes can also be provided for empty groups (produced by list()) and datasets (produced by numeric(0)).

Usage

writehdf5(obj, file, inherent.attributes = FALSE, level = 6, overwrite = TRUE)

Arguments

Details

The function relies on the hdf5r package and on the bit64 package.

The nested list obj should only contain data types available to HDF5. The only exception are data frames, which, if included, are automatically converted to lists. Normally, the list obj and its nested sublists, should all be named lists, i.e. list(a=1, b=2) rather than list(1,2). Unnamed elements are automatically assigned a name 'unnamed_#' in the HDF5 file.

Some data types in R have inherent attributes, such as 'names' for data frames and 'dim' for arrays. By default, these inherent attributes are not written to the HDF5 file. They are, however, automatically regenerated when the HDF5 file is loaded back via readhdf5. The argument inherent.attributes can be used to force writing all attributes, including the inherent ones, to the HDF5 file.

If a structured R list is saved with writehdf5 and then reloaded using readhdf5, the recovered list is identical to the input list up to the ordering of list elements, which is alphabetic by default in HDF5. The only other difference between the input and recovered data occurs when the input data contain data frames, as these are automatically converted do lists. A workaround is to set inherent.attributes=TRUE when writing the HDF5 file. In this case, the reloaded HDF5 data becomes a data frame again.

Value

None

See Also

readhdf5

Examples

# Create example data
input = list(
  group_empty = list(),
  dataset_empty = numeric(0),
  group_dataframe = data.frame(
    ID = as.integer(1:3),
    Name = c("Alice", "Bob", "Charlie"),
    Salary = c(2341.2, 3534.2, 4541.9),
    Employed = c(TRUE, TRUE, FALSE)
  ),
  dataset_parameters = 1.234,
  group_integers = list(int32 = as.integer(123),
                        int64 = bit64::as.integer64(123),
                        vector = bit64::as.integer64((1:3)+12345678912345)),
  dataset_nonnumeric = c(NA, NaN, Inf, -Inf),
  group_mixed = list(
    header = 'test header',
    subgroup1 = list(
      dataset1 = c("A", "%}~&^", "x1y2z3", "", " "),
      dataset2 = matrix(as.integer(1:10), nrow = 2),
      dataset3 = array(runif(30), dim = c(2, 5, 3))
    ),
    subgroup2 = list(date = as.character(as.Date("2025-01-01")),
                     location = 'Perth')
  )
)

# Add attributes to some datasets
attr(input$dataset_empty,'Comment') = 'This is a test file.'
attr(input$group_mixed$subgroup1$dataset3,'Type') = '3D array'
attr(input$group_integers$vector,'Comment') = 'Vector of long integers'
attr(input$group_integers$vector,'Package') = 'bit64'

# Add attributes to some groups
attr(input$group_dataframe,'Company branch') = 'Sales'
attr(input$group_integers,'Comment') = 'Testing different integers'
attr(input$group_empty,'Timestamp') = date()
attr(input$group_empty,'Working directory') = getwd()

if (requireNamespace("hdf5r", quietly = TRUE)) {

  # Write list to HDF5 file
  filename = tempfile()
  writehdf5(input, filename)

  # Read HDF5 file into a new list
  output = readhdf5(filename)

  # Check if input and output lists are identical
  # (up to alphabetic ordering and data frame-to-list conversion)
  print(all.equal(sortlist(input, convert.data.frames = TRUE), output))

  # Write list to HDF5 file again, this time with inherent attributes
  filename = tempfile()
  writehdf5(input, filename, inherent.attributes = TRUE)

  # Read HDF5 file into a new list
  output = readhdf5(filename)

  # Check if input and output lists are identical
  # (up to alphabetic ordering only)
  print(all.equal(sortlist(input), output))
}

Yin-Yang-Yong image

Description

A data set containing a 500-by-500-by-4 array, representing an image of a Yin-Yang-Yong symbol with transparency. This symbol was 'reinvented' many times. The particular version in this package is the one drawn by the package author Danail Obreschkow in 1999.

Usage

yinyangyong

Format

An object of class array of dimension 500 x 500 x 4.

Value

None

Author(s)

Danail Obreschkow

Examples

nplot(asp=1)
rasterImage(yinyangyong,0,0,1,1)