# SciPy Central

## Adaptive sampling of 1D functions

Submitted by: pv, 22 November 2012
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 # License: Creative Commons Zero (almost public domain) http://scpyce.org/cc0 import numpy as np def sample_function(func, points, tol=0.05, min_points=16, max_level=16, sample_transform=None): """ Sample a 1D function to given tolerance by adaptive subdivision. The result of sampling is a set of points that, if plotted, produces a smooth curve with also sharp features of the function resolved. Parameters ---------- func : callable Function func(x) of a single argument. It is assumed to be vectorized. points : array-like, 1D Initial points to sample, sorted in ascending order. These will determine also the bounds of sampling. tol : float, optional Tolerance to sample to. The condition is roughly that the total length of the curve on the (x, y) plane is computed up to this tolerance. min_point : int, optional Minimum number of points to sample. max_level : int, optional Maximum subdivision depth. sample_transform : callable, optional Function w = g(x, y). The x-samples are generated so that w is sampled. Returns ------- x : ndarray X-coordinates y : ndarray Corresponding values of func(x) Notes ----- This routine is useful in computing functions that are expensive to compute, and have sharp features --- it makes more sense to adaptively dedicate more sampling points for the sharp features than the smooth parts. Examples -------- >>> def func(x): ... '''Function with a sharp peak on a smooth background''' ... a = 0.001 ... return x + a**2/(a**2 + x**2) ... >>> x, y = sample_function(func, [-1, 1], tol=1e-3) >>> import matplotlib.pyplot as plt >>> xx = np.linspace(-1, 1, 12000) >>> plt.plot(xx, func(xx), '-', x, y[0], '.') >>> plt.show() """ return _sample_function(func, points, values=None, mask=None, depth=0, tol=tol, min_points=min_points, max_level=max_level, sample_transform=sample_transform) def _sample_function(func, points, values=None, mask=None, tol=0.05, depth=0, min_points=16, max_level=16, sample_transform=None): points = np.asarray(points) if values is None: values = np.atleast_2d(func(points)) if mask is None: mask = Ellipsis if depth > max_level: # recursion limit return points, values x_a = points[...,:-1][...,mask] x_b = points[...,1:][...,mask] x_c = .5*(x_a + x_b) y_c = np.atleast_2d(func(x_c)) x_2 = np.r_[points, x_c] y_2 = np.r_['-1', values, y_c] j = np.argsort(x_2) x_2 = x_2[...,j] y_2 = y_2[...,j] # -- Determine the intervals at which refinement is necessary if len(x_2) < min_points: mask = np.ones([len(x_2)-1], dtype=bool) else: # represent the data as a path in N dimensions (scaled to unit box) if sample_transform is not None: y_2_val = sample_transform(x_2, y_2) else: y_2_val = y_2 p = np.r_['0', x_2[None,:], y_2_val.real.reshape(-1, y_2_val.shape[-1]), y_2_val.imag.reshape(-1, y_2_val.shape[-1]) ] sz = (p.shape[0]-1)//2 xscale = x_2.ptp(axis=-1) yscale = abs(y_2_val.ptp(axis=-1)).ravel() p[0] /= xscale p[1:sz+1] /= yscale[:,None] p[sz+1:] /= yscale[:,None] # compute the length of each line segment in the path dp = np.diff(p, axis=-1) s = np.sqrt((dp**2).sum(axis=0)) s_tot = s.sum() # compute the angle between consecutive line segments dp /= s dcos = np.arccos(np.clip((dp[:,1:] * dp[:,:-1]).sum(axis=0), -1, 1)) # determine where to subdivide: the condition is roughly that # the total length of the path (in the scaled data) is computed # to accuracy tol dp_piece = dcos * .5*(s[1:] + s[:-1]) mask = (dp_piece > tol * s_tot) mask = np.r_[mask, False] mask[1:] |= mask[:-1].copy() # -- Refine, if necessary if mask.any(): return _sample_function(func, x_2, y_2, mask, tol=tol, depth=depth+1, min_points=min_points, max_level=max_level, sample_transform=sample_transform) else: return x_2, y_2 
Software license: Creative Commons Zero. No rights reserved.
Users have permission to do anything with the code and other material on this page. (More details)
More information:

Sample a 1D function to given tolerance by adaptive subdivision.

The result of sampling is a set of points that, if plotted, produces a smooth curve with also sharp features of the function resolved.

This routine is useful in computing functions that are expensive to compute, and have sharp features — it makes more sense to adaptively dedicate more sampling points for the sharp features than the smooth parts.

# Examples

>>> def func(x):
...     '''Function with a sharp peak on a smooth background'''
...     a = 0.001
...     return x + a**2/(a**2 + x**2)
...
>>> x, y = sample_function(func, [-1, 1], tol=1e-3)

>>> import matplotlib.pyplot as plt
>>> xx = np.linspace(-1, 1, 12000)
>>> plt.plot(xx, func(xx), '-', x, y[0], '.')
>>> plt.show()


Tags: sampling