import matplotlib.pyplot as plt from numpy import * import sys import re #limits of function: #put function here func=lambda x: 1 if(len(sys.argv) == 4): a_s = sys.argv[2] b_s = sys.argv[3] func_string = sys.argv[1] print("a_s: {}".format(a_s)) print("b_s: {}".format(b_s)) a = float(re.sub("a=", "", a_s)) b = float(re.sub("b=", "", b_s)) N= 1000 def plotter(): x = linspace(a,b, N) plt.figure(figsize=(10,10)) plt.plot(x,func(x), label=func_string) #plt.plot(x,g(x), label='g(x)') plt.legend(loc='best') plt.grid() plt.show() plt.savefig("/tmp/plot.png", dpi=150) # 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 if __name__ == '__main__': #plotter() x, y = sample_function(func, [a, b], tol=1e-6) xx = np.linspace(a, b, 300) plt.figure(figsize=((10,10))) plt.plot(x, y[0], '-',linewidth=1.5) plt.grid() plt.savefig("/tmp/plot.png", dpi=400) #weightss = np.empty_like(xx) #for i in range(len(xx)): # if i == 0: # weightss[i] = xx[i+1]-x[i] # if i == len(xx)-1: # weightss[i] = xx[-1]-xx[-2] # else: # weightss[i] = 0.5*(xx[i+1]-xx[i-1]) #avg = average(abs(func(xx))) #plt.ylim(top=avg, bottom=-1.0*avg)