Source code for lsst.sims.movingObjects.chebyshevUtils

"""Utilities to estimate and evaluate Chebyshev coefficients of a function.

Implementation of Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310
"""

import numpy as np

__all__ = ['chebeval', 'chebfit', 'makeChebMatrix', 'makeChebMatrixOnlyX']

# Evaluation routine.

[docs]def chebeval(x, p, interval=(-1., 1.), doVelocity=True, mask=False): """Evaluate a Chebyshev series and first derivative at points x. If p is of length n + 1, this function returns: y_hat(x) = p_0 * T_0(x*) + p_1 * T_1(x*) + ... + p_n * T_n(x*) where T_n(x*) are the orthogonal Chebyshev polynomials of the first kind, defined on the interval [-1, 1] and p_n are the coefficients. The scaled variable x* is defined on the [-1, 1] interval such that (x*) = (2*x - a - b)/(b - a), and x is defined on the [a, b] interval. Parameters ---------- x: scalar or numpy.ndarray Points at which to evaluate the polynomial. p: numpy.ndarray Chebyshev polynomial coefficients, as returned by chebfit. interval: 2-element list/tuple Bounds the x-interval on which the Chebyshev coefficients were fit. doVelocity: bool If True, compute the first derivative at points x. mask: bool If True, return Nans when the x goes beyond 'interval'. If False, extrapolate fit beyond 'interval' limits. Returns ------- scalar or numpy.ndarray, scalar or numpy.ndarray Y (position) and velocity values (if computed) """ if len(interval) != 2: raise RuntimeError("interval must have length 2") intervalBegin = np.float(interval[0]) intervalEnd = np.float(interval[-1]) t = 2. * np.array(x, dtype=np.float64) - intervalBegin - intervalEnd t /= intervalEnd - intervalBegin y = 0. v = 0. y0 = np.ones_like(t) y1 = t v0 = np.zeros_like(t) v1 = np.ones_like(t) v2 = 4. * t t = 2. * t N = len(p) if doVelocity: for i in np.arange(0, N, 2): if i == N - 1: y1 = 0. v1 = 0. j = min(i + 1, N - 1) y += p[i] * y0 + p[j] * y1 v += p[i] * v0 + p[j] * v1 y2 = t * y1 - y0 y3 = t * y2 - y1 v2 = t * v1 - v0 + 2 * y1 v3 = t * v2 - v1 + 2 * y2 y0 = y2 y1 = y3 v0 = v2 v1 = v3 if mask: mask = np.where((x < intervalBegin) | (x > intervalEnd), True, False) y = np.where(mask, np.nan, y) v = np.where(mask, np.nan, v) return y, 2 * v / (intervalEnd - intervalBegin) else: for i in np.arange(0, N, 2): if i == N - 1: y1 = 0. j = min((i + 1), (N - 1)) y += p[i] * y0 + p[j] * y1 y0 = t * y1 - y0 y1 = t * y0 - y1 if mask: mask = np.where((x < intervalBegin) | (x > intervalEnd), True, False) y = np.where(mask, np.nan, y) return y, None
# Fitting routines.
[docs]def makeChebMatrix(nPoints, nPoly, weight=0.16): """Compute C1^(-1)C2 using Newhall89 approach. Utility function for fitting chebyshev polynomials to x(t) and dx/dt(t) forcing equality at the end points. This function computes the matrix (C1^(-1)C2). Multiplying this matrix by the x and dx/dt values to be fit produces the chebyshev coefficient. This function need only be called once for a given polynomial degree and number of points. The matrices returned are of shape(nPoints+1)x(nPoly). The coefficients fitting the nPoints+1 points, X, are found by: A = xMultiplier * x + dxMultiplier * dxdt if derivative information is known, or A = xMultiplier * x if no derivative information is known. The xMultiplier matrices are different, depending on whether derivative information is known. Use function makeChebMatrixOnlyX if derviative is not known. See Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310 for details. Parameters ---------- nPoints: int Number of point to be fits. Must be greater than 2. nPoly: int Number of polynomial terms. Polynomial degree + 1 weight: float, optional Weight to allow control of relative effectos of position and velocity values. Newhall80 found best results are obtained with velocity weighted at 0.4 relative to position, giving W the form (1.0, 0.16, 1.0, 0.16,...) Returns ------- numpy.ndarray xMultiplier, C1^(-1)C2 even rows of shape (nPoints+1)x(nPoly) to be multiplied by x values. numpy.ndarray dxMultiplier, C1^(-1)C2 odd rows of shape (nPoints+1)x(nPoly) to be multiplied by dx/dy values """ tmat = np.zeros([nPoints, nPoly]) tdot = np.zeros([nPoints, nPoly]) cj = np.zeros([nPoly]) xj = np.linspace(1, -1, nPoints) for i in np.arange(0, nPoly): cj[:] = 0 cj[i] = 1 y, v = chebeval(xj, cj) tmat[:, i] = y tdot[:, i] = v # make matrix T*W tw = np.zeros([nPoly, nPoints, 2]) tw[:, :, 0] = tmat.transpose() tw[:, :, 1] = tdot.transpose() * weight # make matrix T*WT twt = np.dot(tw[:, :, 0], tmat) + np.dot(tw[:, :, 1], tdot) tw = tw.reshape(nPoly, 2 * nPoints) # insert matrix T*W in matrix C2 c2 = np.zeros([nPoly + 4, 2 * nPoints]) c2[0:nPoly] = tw c2[nPoly, 0] = 1 c2[nPoly + 1, 1] = 1 c2[nPoly + 2, -2] = 1 c2[nPoly + 3, -1] = 1 # insert matrix T*WT in matrix C1 c1 = np.zeros([nPoly + 4, nPoly + 4]) c1[0:nPoly, 0:nPoly] = twt c1[nPoly + 0, 0:nPoly] = tmat[0] c1[nPoly + 1, 0:nPoly] = tdot[0] c1[nPoly + 2, 0:nPoly] = tmat[-1] c1[nPoly + 3, 0:nPoly] = tdot[-1] c1[0:nPoly, nPoly:] = c1[nPoly:, 0:nPoly].transpose() c1inv = np.linalg.inv(c1) c1c2 = np.dot(c1inv, c2) c1c2 = c1c2.reshape(nPoly + 4, nPoints, 2) c1c2 = c1c2[:, ::-1, :] c1c2 = c1c2.reshape(nPoly + 4, 2 * nPoints) # separate even rows for x, and odd rows for dx/dt return c1c2[0:nPoly, 0::2], c1c2[0:nPoly, 1::2]
[docs]def makeChebMatrixOnlyX(nPoints, nPoly): """Compute C1^(-1)C2 using Newhall89 approach without dx/dt Compute xMultiplier using only the equality constraint of the x-values at the endpoints. To be used when first derivatives are not available. If chebyshev approximations are strung together piecewise only the x-values and not the first derivatives will be continuous at the boundaries. Multiplying this matrix by the x-values to be fit produces the chebyshev coefficients. This function need only be called once for a given polynomial degree and number of points. See Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310. Parameters ---------- nPoints : int Number of point to be fits. Must be greater than 2. nPoly : int Number of polynomial terms. Polynomial degree + 1 Returns ------- numpy.ndarray xMultiplier, Even rows of C1^(-1)C2 w/ shape (nPoints+1)x(nPoly) to be multiplied by x values """ tmat = np.zeros([nPoints, nPoly]) cj = np.zeros([nPoly]) xj = np.linspace(1, -1, nPoints) for i in range(0, nPoly): cj[:] = 0 cj[i] = 1 tmat[:, i], v = chebeval(xj, cj) # Augment matrix T to get matrix C2 c2 = np.zeros([nPoly + 2, nPoints]) c2[0:nPoly] = tmat.transpose() c2[nPoly, 0] = 1 c2[nPoly + 1, nPoints - 1] = 1 # Augment matrix T*WT to get the matrix C1 c1 = np.zeros([nPoly + 2, nPoly + 2]) c1[0:nPoly, 0:nPoly] = np.dot(tmat.transpose(), tmat) c1[nPoly + 0, 0:nPoly] = tmat[0] c1[nPoly + 1, 0:nPoly] = tmat[-1] c1[0:nPoly, nPoly:] = c1[nPoly:, 0:nPoly].transpose() c1inv = np.linalg.inv(c1) # C1^(-1) C2 c1c2 = np.dot(c1inv, c2) c1c2 = c1c2.reshape(nPoly + 2, nPoints) c1c2 = c1c2[:, ::-1] return c1c2[0:nPoly]
[docs]def chebfit(t, x, dxdt=None, xMultiplier=None, dxMultiplier=None, nPoly=7): """Fit Chebyshev polynomial constrained at endpoints using Newhall89 approach. Return Chebyshev coefficients and statistics from fit to array of positions (x) and optional velocities (dx/dt). If both the function and its derivative are specified, then the value and derivative of the interpolating polynomial at the endpoints will be exactly equal to the input endpoint values. Many approximations may be piecewise strung together and the function value and its first derivative will be continuous across boundaries. If derivatives are not provided, only the function value will be continuous across boundaries. If xMultiplier and dxMultiplier are not provided or are an inappropriate shape for t and x, they will be recomputed. See Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310 for details. Parameters ---------- t : numpy.ndarray Array of regularly sampled independent variable (e.g. time) x : numpy.ndarray Array of regularly sampled dependent variable (e.g. declination) dxdt : numpy.ndarray, optional Optionally, array of first derivatives of x with respect to t, at the same grid points. (e.g. sky velocity ddecl/dt) xMultiplier : numpy.ndarray, optional Optional 2D Matrix with rows of C1^(-1)C2 corresponding to x. Use makeChebMatrix to compute dxMultiplier : numpy.ndarray, optional Optional 2D Matrix with rows of C1^(-1)C2 corresponding to dx/dt. Use makeChebMatrix to compute nPoly : int, optional Number of polynomial terms. Degree + 1. Must be >=2 and <=2*nPoints, when derivative information is specified, or <=nPoints, when no derivative information is specified. Default = 7. Returns ------- numpy.ndarray Array of chebyshev coefficients with length=nPoly. numpy.ndarray Array of residuals of the tabulated function x minus the approximated function. float The rms of the residuals in the fit. float The maximum of the residals to the fit. """ nPoints = len(t) if len(x) != nPoints: raise ValueError("length of x (%s) != length of t (%s)" % (len(x), nPoints)) if dxdt is None: if nPoly > nPoints: raise RuntimeError('Without velocity constraints, nPoly (%d) must be less than %s' % (nPoly, nPoints)) if nPoly < 2: raise RuntimeError('Without velocity constraints, nPoly (%d) must be greater than 2' % nPoly) else: if nPoly > 2 * nPoints: raise RuntimeError('nPoly (%d) must be less than %s (%d)' % (nPoly, '2 * nPoints', 2 * (nPoints))) if nPoly < 4: raise RuntimeError('nPoly (%d) must be greater than 4' % nPoly) # Recompute C1invX2 if xMultiplier and dxMultiplier are None or # they are not appropriate for sizes of input positions and velocities. if xMultiplier is None: redoX = True else: redoX = (xMultiplier.shape[1] != nPoints) | (xMultiplier.shape[0] != nPoly) if dxMultiplier is None: redoV = True else: redoV = (dxMultiplier.shape[1] != nPoints) | (dxMultiplier.shape[0] != nPoly) if (dxdt is None) & redoX: xMultiplier = makeChebMatrixOnlyX(nPoints, nPoly) if (dxdt is not None) & (redoV | redoX): xMultiplier, dxMultiplier = makeChebMatrix(nPoints, nPoly) if x.size != nPoints: raise RuntimeError("Not enough elements in X") tInterval = np.array([t[0], t[-1]]) - t[0] tScaled = t - t[0] # Compute the X portion of the coefficients a_n = np.dot(xMultiplier, x) # Compute statistics # for x and dxdt if it is available if dxdt is not None: a_n = a_n + np.dot(dxMultiplier, dxdt * (tInterval[1] - tInterval[0]) / 2.) xApprox, dxApprox = chebeval(tScaled, a_n, interval=tInterval) else: # Statistics for x only xApprox, _ = chebeval(tScaled, a_n, interval=tInterval, doVelocity=False) residuals = x - xApprox se = np.sum(residuals**2) rms = np.sqrt(se / (nPoints - 1)) maxresid = np.max(np.abs(residuals)) return a_n, residuals, rms, maxresid