diff --git a/xs/src/BSpline/BSpline.cpp b/xs/src/BSpline/BSpline.cpp new file mode 100644 index 000000000..6cce57319 --- /dev/null +++ b/xs/src/BSpline/BSpline.cpp @@ -0,0 +1,879 @@ +// -*- mode: c++; c-basic-offset: 4; -*- +/************************************************************************ + * Copyright 2009 University Corporation for Atmospheric Research. + * All rights reserved. + * + * Use of this code is subject to the standard BSD license: + * + * http://www.opensource.org/licenses/bsd-license.html + * + * See the COPYRIGHT file in the source distribution for the license text, + * or see this web page: + * + * http://www.eol.ucar.edu/homes/granger/bspline/doc/ + * + *************************************************************************/ + +/** + * @file + * + * This file defines the implementation for the BSpline and BSplineBase + * templates. + **/ +#include "BSpline.h" +#include "BandedMatrix.h" + +#include +#include +#include +#include +#include +#include +#include + + + + + +/* + * Original BSplineBase.cpp start here + */ + +/* + * This class simulates a namespace for private symbols used by this template + * implementation which should not pollute the global namespace. + */ +class my +{ + public: + template static inline + T abs(const T t) + { + return (t < 0) ? -t : t; + } + + template static inline + const T& min(const T& a, + const T& b) + { + return (a < b) ? a : b; + } + + template static inline + const T& max(const T& a, + const T& b) + { + return (a > b) ? a : b; + } +}; + +////////////////////////////////////////////////////////////////////// +template class Matrix : public BandedMatrix +{ + public: + Matrix &operator +=(const Matrix &B) + { + Matrix &A = *this; + typename Matrix::size_type M = A.num_rows(); + typename Matrix::size_type N = A.num_cols(); + + assert(M==B.num_rows()); + assert(N==B.num_cols()); + + typename Matrix::size_type i, j; + for (i=0; i::operator= (e); + return *this; + } + + }; + ////////////////////////////////////////////////////////////////////// + // Our private state structure, which hides our use of some matrix + // template classes. + +template struct BSplineBaseP +{ + typedef Matrix MatrixT; + + MatrixT Q; // Holds P+Q and its factorization + std::vector X; + std::vector Nodes; +}; + +////////////////////////////////////////////////////////////////////// + +// This array contains the beta parameter for the boundary condition +// constraints. The boundary condition type--0, 1, or 2--is the first +// index into the array, followed by the index of the endpoints. See the +// Beta() method. +template const double BSplineBase::BoundaryConditions[3][4] = + { + // 0 1 M-1 M + { + -4, + -1, + -1, + -4 }, + { + 0, + 1, + 1, + 0 }, + { + 2, + -1, + -1, + 2 } }; +////////////////////////////////////////////////////////////////////// +template inline bool BSplineBase::Debug(int on) +{ + static bool debug = false; + if (on >= 0) + debug = (on > 0); + return debug; +} + +////////////////////////////////////////////////////////////////////// +template const double BSplineBase::BS_PI = 3.1415927; +////////////////////////////////////////////////////////////////////// +template const char * BSplineBase::ImplVersion() +{ + return ("$Id: BSpline.cpp 6352 2008-05-05 04:40:39Z martinc $"); +} + +////////////////////////////////////////////////////////////////////// +template const char * BSplineBase::IfaceVersion() +{ + return (_BSPLINEBASE_IFACE_ID); +} + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +template BSplineBase::~BSplineBase() +{ + delete base; +} + +// This is a member-wise copy except for replacing our +// private base structure with the source's, rather than just copying +// the pointer. But we use the compiler's default copy constructor for +// constructing our BSplineBaseP. +template BSplineBase::BSplineBase(const BSplineBase &bb) : + K(bb.K), BC(bb.BC), OK(bb.OK), base(new BSplineBaseP(*bb.base)) +{ + xmin = bb.xmin; + xmax = bb.xmax; + alpha = bb.alpha; + waveLength = bb.waveLength; + DX = bb.DX; + M = bb.M; + NX = base->X.size(); +} +////////////////////////////////////////////////////////////////////// +template BSplineBase::BSplineBase(const T *x, + int nx, + double wl, + int bc, + int num_nodes) : + NX(0), K(2), OK(false), base(new BSplineBaseP) +{ + setDomain(x, nx, wl, bc, num_nodes); +} + +////////////////////////////////////////////////////////////////////// +// Methods +template bool BSplineBase::setDomain(const T *x, + int nx, + double wl, + int bc, + int num_nodes) +{ + if ((nx <= 0) || (x == 0) || (wl< 0) || (bc< 0) || (bc> 2)) { + return false; + } + OK = false; + waveLength = wl; + BC = bc; + // Copy the x array into our storage. + base->X.resize(nx); + std::copy(x, x+nx, base->X.begin()); + NX = base->X.size(); + + // The Setup() method determines the number and size of node intervals. + if (Setup(num_nodes)) { + if (Debug()) { + std::cerr << "Using M node intervals: " << M << " of length DX: " + << DX << std::endl; + std::cerr << "X min: " << xmin << " ; X max: " << xmax << std::endl; + std::cerr << "Data points per interval: " << (float)NX/(float)M + << std::endl; + std::cerr << "Nodes per wavelength: " << (float)waveLength + /(float)DX << std::endl; + std::cerr << "Derivative constraint degree: " << K << std::endl; + } + + // Now we can calculate alpha and our Q matrix + alpha = Alpha(waveLength); + if (Debug()) { + std::cerr << "Cutoff wavelength: " << waveLength << " ; " + << "Alpha: " << alpha << std::endl; + std::cerr << "Calculating Q..." << std::endl; + } + calculateQ(); + if (Debug() && M < 30) { + std::cerr.fill(' '); + std::cerr.precision(2); + std::cerr.width(5); + std::cerr << base->Q << std::endl; + } + + if (Debug()) + std::cerr << "Calculating P..." << std::endl; + addP(); + if (Debug()) { + std::cerr << "Done." << std::endl; + if (M < 30) { + std::cerr << "Array Q after addition of P." << std::endl; + std::cerr << base->Q; + } + } + + // Now perform the LU factorization on Q + if (Debug()) + std::cerr << "Beginning LU factoring of P+Q..." << std::endl; + if (!factor()) { + if (Debug()) + std::cerr << "Factoring failed." << std::endl; + } else { + if (Debug()) + std::cerr << "Done." << std::endl; + OK = true; + } + } + return OK; +} +////////////////////////////////////////////////////////////////////// +/* + * Calculate the alpha parameter given a wavelength. + */ +template double BSplineBase::Alpha(double wl) +{ + // K is the degree of the derivative constraint: 1, 2, or 3 + double a = (double) (wl / (2 * BS_PI * DX)); + a *= a; // a^2 + if (K == 2) + a = a * a; // a^4 + else if (K == 3) + a = a * a * a; // a^6 + return a; +} +////////////////////////////////////////////////////////////////////// +/* + * Return the correct beta value given the node index. The value depends + * on the node index and the current boundary condition type. + */ +template inline double BSplineBase::Beta(int m) +{ + if (m > 1 && m < M-1) + return 0.0; + if (m >= M-1) + m -= M-3; + assert(0 <= BC && BC <= 2); + assert(0 <= m && m <= 3); + return BoundaryConditions[BC][m]; +} +////////////////////////////////////////////////////////////////////// +/* + * Given an array of y data points defined over the domain + * of x data points in this BSplineBase, create a BSpline + * object which contains the smoothed curve for the y array. + */ +template BSpline * BSplineBase::apply(const T *y) +{ + return new BSpline (*this, y); +} +////////////////////////////////////////////////////////////////////// +/* + * Evaluate the closed basis function at node m for value x, + * using the parameters for the current boundary conditions. + */ +template double BSplineBase::Basis(int m, + T x) +{ + double y = 0; + double xm = xmin + (m * DX); + double z = my::abs((double)(x - xm) / (double)DX); + if (z < 2.0) { + z = 2 - z; + y = 0.25 * (z*z*z); + z -= 1.0; + if (z > 0) + y -= (z*z*z); + } + + // Boundary conditions, if any, are an additional addend. + if (m == 0 || m == 1) + y += Beta(m) * Basis(-1, x); + else if (m == M-1 || m == M) + y += Beta(m) * Basis(M+1, x); + + return y; +} +////////////////////////////////////////////////////////////////////// +/* + * Evaluate the deriviative of the closed basis function at node m for + * value x, using the parameters for the current boundary conditions. + */ +template double BSplineBase::DBasis(int m, + T x) +{ + double dy = 0; + double xm = xmin + (m * DX); + double delta = (double)(x - xm) / (double)DX; + double z = my::abs(delta); + if (z < 2.0) { + z = 2.0 - z; + dy = 0.25 * z * z; + z -= 1.0; + + if (z > 0) { + dy -= z * z; + } + dy *= ((delta > 0) ? -1.0 : 1.0) * 3.0 / DX; + } + + // Boundary conditions, if any, are an additional addend. + if (m == 0 || m == 1) + dy += Beta(m) * DBasis(-1, x); + else if (m == M-1 || m == M) + dy += Beta(m) * DBasis(M+1, x); + + return dy; +} +////////////////////////////////////////////////////////////////////// +template double BSplineBase::qDelta(int m1, + int m2) +/* + * Return the integral of the product of the basis function derivative + * restricted to the node domain, 0 to M. + */ +{ + // These are the products of the Kth derivative of the + // normalized basis functions + // given a distance m nodes apart, qparts[K-1][m], 0 <= m <= 3 + // Each column is the integral over each unit domain, -2 to 2 + static const double qparts[3][4][4] = + { + { + { + 0.11250, + 0.63750, + 0.63750, + 0.11250 }, + { + 0.00000, + 0.13125, + -0.54375, + 0.13125 }, + { + 0.00000, + 0.00000, + -0.22500, + -0.22500 }, + { + 0.00000, + 0.00000, + 0.00000, + -0.01875 } }, + { + { + 0.75000, + 2.25000, + 2.25000, + 0.75000 }, + { + 0.00000, + -1.12500, + -1.12500, + -1.12500 }, + { + 0.00000, + 0.00000, + 0.00000, + 0.00000 }, + { + 0.00000, + 0.00000, + 0.00000, + 0.37500 } }, + { + { + 2.25000, + 20.25000, + 20.25000, + 2.25000 }, + { + 0.00000, + -6.75000, + -20.25000, + -6.75000 }, + { + 0.00000, + 0.00000, + 6.75000, + 6.75000 }, + { + 0.00000, + 0.00000, + 0.00000, + -2.25000 } } }; + + if (m1 > m2) + std::swap(m1, m2); + + if (m2 - m1 > 3) + return 0.0; + + double q = 0; + for (int m = my::max(m1-2, 0); m < my::min(m1+2, M); ++m) + q += qparts[K-1][m2-m1][m-m1+2]; + return q * alpha; +} +////////////////////////////////////////////////////////////////////// +template void BSplineBase::calculateQ() +{ + Matrix &Q = base->Q; + Q.setup(M+1, 3); + Q = 0; + if (alpha == 0) + return; + + // First fill in the q values without the boundary constraints. + int i; + for (i = 0; i <= M; ++i) { + Q[i][i] = qDelta(i, i); + for (int j = 1; j < 4 && i+j <= M; ++j) { + Q[i][i+j] = Q[i+j][i] = qDelta(i, i+j); + } + } + + // Now add the boundary constraints: + // First the upper left corner. + float b1, b2, q; + for (i = 0; i <= 1; ++i) { + b1 = Beta(i); + for (int j = i; j < i+4; ++j) { + b2 = Beta(j); + assert(j-i >= 0 && j - i < 4); + q = 0.0; + if (i+1 < 4) + q += b2*qDelta(-1, i); + if (j+1 < 4) + q += b1*qDelta(-1, j); + q += b1*b2*qDelta(-1, -1); + Q[j][i] = (Q[i][j] += q); + } + } + + // Then the lower right + for (i = M-1; i <= M; ++i) { + b1 = Beta(i); + for (int j = i - 3; j <= i; ++j) { + b2 = Beta(j); + q = 0.0; + if (M+1-i < 4) + q += b2*qDelta(i, M+1); + if (M+1-j < 4) + q += b1*qDelta(j, M+1); + q += b1*b2*qDelta(M+1, M+1); + Q[j][i] = (Q[i][j] += q); + } + } +} +////////////////////////////////////////////////////////////////////// +template void BSplineBase::addP() +{ + // Add directly to Q's elements + Matrix &P = base->Q; + std::vector &X = base->X; + + // For each data point, sum the product of the nearest, non-zero Basis + // nodes + int mx, m, n, i; + for (i = 0; i < NX; ++i) { + // Which node does this put us in? + T &x = X[i]; + mx = (int)((x - xmin) / DX); + + // Loop over the upper triangle of nonzero basis functions, + // and add in the products on each side of the diagonal. + for (m = my::max(0, mx-1); m <= my::min(M, mx+2); ++m) { + float pn; + float pm = Basis(m, x); + float sum = pm * pm; + P[m][m] += sum; + for (n = m+1; n <= my::min(M, mx+2); ++n) { + pn = Basis(n, x); + sum = pm * pn; + P[m][n] += sum; + P[n][m] += sum; + } + } + } +} +////////////////////////////////////////////////////////////////////// +template bool BSplineBase::factor() +{ + Matrix &LU = base->Q; + + if (LU_factor_banded(LU, 3) != 0) { + if (Debug()) + std::cerr << "LU_factor_banded() failed." << std::endl; + return false; + } + if (Debug() && M < 30) + std::cerr << "LU decomposition: " << std::endl << LU << std::endl; + return true; +} +////////////////////////////////////////////////////////////////////// +template inline double BSplineBase::Ratiod(int &ni, + double &deltax, + double &ratiof) +{ + deltax = (xmax - xmin) / ni; + ratiof = waveLength / deltax; + double ratiod = (double) NX / (double) (ni + 1); + return ratiod; +} +////////////////////////////////////////////////////////////////////// +// Setup the number of nodes (and hence deltax) for the given domain and +// cutoff wavelength. According to Ooyama, the derivative constraint +// approximates a lo-pass filter if the cutoff wavelength is about 4*deltax +// or more, but it should at least be 2*deltax. We can increase the number +// of nodes to increase the number of nodes per cutoff wavelength. +// However, to get a reasonable representation of the data, the setup +// enforces at least as many nodes as data points in the domain. (This +// constraint assumes reasonably even distribution of data points, since +// its really the density of data points which matters.) +// +// Return zero if the setup fails, non-zero otherwise. +// +// The algorithm in this routine is mostly taken from the FORTRAN +// implementation by James Franklin, NOAA/HRD. +// +template bool BSplineBase::Setup(int num_nodes) +{ + std::vector &X = base->X; + + // Find the min and max of the x domain + xmin = X[0]; + xmax = X[0]; + + int i; + for (i = 1; i < NX; ++i) { + if (X[i] < xmin) + xmin = X[i]; + else if (X[i] > xmax) + xmax = X[i]; + } + if (Debug()) + std::cerr << "Xmax=" << xmax << ", Xmin=" << xmin << std::endl; + + // Number of node intervals (number of spline nodes - 1). + int ni; + double deltax; + + if (num_nodes >= 2) { + // We've been told explicitly the number of nodes to use. + ni = num_nodes - 1; + if (waveLength == 0) { + waveLength = 1.0; + } + if (Debug()) + { + std::cerr << "Num nodes explicitly given as " << num_nodes + << ", wavelength set to " << waveLength << std::endl; + } + } else if (waveLength == 0) { + // Turn off frequency constraint and just set two node intervals per + // data point. + ni = NX * 2; + waveLength = 1; + if (Debug()) + { + std::cerr << "Frequency constraint disabled, using 2 intervals " + << "per node: " << ni << " intervals, wavelength=" + << waveLength << std::endl; + } + } else if (waveLength > xmax - xmin) { + if (Debug()) + { + std::cerr << "Wavelength " << waveLength << " exceeds X span: " + << xmin << " - " << xmax << std::endl; + } + return (false); + } else { + if (Debug()) + { + std::cerr << "Searching for a reasonable number of " + << "intervals for wavelength " << waveLength + << " while keeping at least 2 intervals per " + << "wavelength ..." << std::endl; + } + // Minimum acceptable number of node intervals per cutoff wavelength. + static const double fmin = 2.0; + + // Start out at a minimum number of intervals, meaning the maximum + // number of points per interval, then work up to the maximum + // number of intervals for which the intervals per wavelength is + // still adequate. I think the minimum must be more than 2 since + // the basis function is evaluated on multiple nodes. + ni = 5; + + double ratiof; // Nodes per wavelength for current deltax + double ratiod; // Points per node interval + + // Increase the number of node intervals until we reach the minimum + // number of intervals per cutoff wavelength, but only as long as + // we can maintain at least one point per interval. + do { + if (Ratiod(++ni, deltax, ratiof) < 1.0) + { + if (Debug()) + { + std::cerr << "At " << ni << " intervals, fewer than " + << "one point per interval, and " + << "intervals per wavelength is " + << ratiof << "." << std::endl; + } + return false; + } + } while (ratiof < fmin); + + // Now increase the number of intervals until we have at least 4 + // intervals per cutoff wavelength, but only as long as we can + // maintain at least 2 points per node interval. There's also no + // point to increasing the number of intervals if we already have + // 15 or more nodes per cutoff wavelength. + // + do { + if ((ratiod = Ratiod(++ni, deltax, ratiof)) < 1.0 || ratiof > 15.0) { + --ni; + break; + } + } while (ratiof < 4 || ratiod > 2.0); + + if (Debug()) + { + std::cerr << "Found " << ni << " intervals, " + << "length " << deltax << ", " + << ratiof << " nodes per wavelength " << waveLength + << ", " + << ratiod << " data points per interval." << std::endl; + } + } + + // Store the calculations in our state + M = ni; + DX = (xmax - xmin) / ni; + + return (true); +} +////////////////////////////////////////////////////////////////////// +template const T * BSplineBase::nodes(int *nn) +{ + if (base->Nodes.size() == 0) { + base->Nodes.reserve(M+1); + for (int i = 0; i <= M; ++i) { + base->Nodes.push_back(xmin + (i * DX)); + } + } + + if (nn) + *nn = base->Nodes.size(); + + assert(base->Nodes.size() == (unsigned)(M+1)); + return &base->Nodes[0]; +} +////////////////////////////////////////////////////////////////////// +template std::ostream &operator<<(std::ostream &out, + const std::vector &c) +{ + for (typename std::vector::const_iterator it = c.begin(); it < c.end(); ++it) + out << *it << ", "; + out << std::endl; + return out; +} + + + + + +/* + * Original BSpline.cpp start here + */ + + +////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////// +// BSpline Class +////////////////////////////////////////////////////////////////////// + +template struct BSplineP { + std::vector spline; + std::vector A; +}; + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +/* + * This BSpline constructor constructs and sets up a new base and + * solves for the spline curve coeffiecients all at once. + */ +template BSpline::BSpline(const T *x, + int nx, + const T *y, + double wl, + int bc_type, + int num_nodes) : + BSplineBase(x, nx, wl, bc_type, num_nodes), s(new BSplineP) { + solve(y); +} +////////////////////////////////////////////////////////////////////// +/* + * Create a new spline given a BSplineBase. + */ +template BSpline::BSpline(BSplineBase &bb, + const T *y) : + BSplineBase(bb), s(new BSplineP) { + solve(y); +} +////////////////////////////////////////////////////////////////////// +/* + * (Re)calculate the spline for the given set of y values. + */ +template bool BSpline::solve(const T *y) { + if (!OK) + return false; + + // Any previously calculated curve is now invalid. + s->spline.clear(); + OK = false; + + // Given an array of data points over x and its precalculated + // P+Q matrix, calculate the b vector and solve for the coefficients. + std::vector &B = s->A; + std::vector &A = s->A; + A.clear(); + A.resize(M+1); + + if (Debug()) + std::cerr << "Solving for B..." << std::endl; + + // Find the mean of these data + mean = 0.0; + int i; + for (i = 0; i < NX; ++i) { + mean += y[i]; + } + mean = mean / (double)NX; + if (Debug()) + std::cerr << "Mean for y: " << mean << std::endl; + + int mx, m, j; + for (j = 0; j < NX; ++j) { + // Which node does this put us in? + T &xj = base->X[j]; + T yj = y[j] - mean; + mx = (int)((xj - xmin) / DX); + + for (m = my::max(0, mx-1); m <= my::min(mx+2, M); ++m) { + B[m] += yj * Basis(m, xj); + } + } + + if (Debug() && M < 30) { + std::cerr << "Solution a for (P+Q)a = b" << std::endl; + std::cerr << " b: " << B << std::endl; + } + + // Now solve for the A vector in place. + if (LU_solve_banded(base->Q, A, 3) != 0) { + if (Debug()) + std::cerr << "LU_solve_banded() failed." << std::endl; + } else { + OK = true; + if (Debug()) + std::cerr << "Done." << std::endl; + if (Debug() && M < 30) { + std::cerr << " a: " << A << std::endl; + std::cerr << "LU factor of (P+Q) = " << std::endl << base->Q + << std::endl; + } + } + return (OK); +} +////////////////////////////////////////////////////////////////////// +template BSpline::~BSpline() { + delete s; +} +////////////////////////////////////////////////////////////////////// +template T BSpline::coefficient(int n) { + if (OK) + if (0 <= n && n <= M) + return s->A[n]; + return 0; +} +////////////////////////////////////////////////////////////////////// +template T BSpline::evaluate(T x) { + T y = 0; + if (OK) { + int n = (int)((x - xmin)/DX); + for (int i = my::max(0, n-1); i <= my::min(M, n+2); ++i) { + y += s->A[i] * Basis(i, x); + } + y += mean; + } + return y; +} +////////////////////////////////////////////////////////////////////// +template T BSpline::slope(T x) { + T dy = 0; + if (OK) { + int n = (int)((x - xmin)/DX); + for (int i = my::max(0, n-1); i <= my::min(M, n+2); ++i) { + dy += s->A[i] * DBasis(i, x); + } + } + return dy; +} + + + + + + + +/// Instantiate BSplineBase +template class BSplineBase; +//template class BSplineBase; + +/// Instantiate BSpline +template class BSpline; +//template class BSpline; diff --git a/xs/src/BSpline/BSpline.h b/xs/src/BSpline/BSpline.h new file mode 100644 index 000000000..9114385c1 --- /dev/null +++ b/xs/src/BSpline/BSpline.h @@ -0,0 +1,452 @@ +/* -*- mode: c++; c-basic-offset: 4; -*- */ +/************************************************************************ + * Copyright 2009 University Corporation for Atmospheric Research. + * All rights reserved. + * + * Use of this code is subject to the standard BSD license: + * + * http://www.opensource.org/licenses/bsd-license.html + * + * See the COPYRIGHT file in the source distribution for the license text, + * or see this web page: + * + * http://www.eol.ucar.edu/homes/granger/bspline/doc/ + * + *************************************************************************/ + +#ifndef BSPLINE_H +#define BSPLINE_H + +#include + + +/* + * Warning: original BSplineBase.h/cpp merged into BSpline.h/cpp to avoid dependency issues caused by Build::WithXSpp which tries to compile all .cpp files in /src + * Original BSplineBase.h starts here!!! + */ + + + +#ifndef _BSPLINEBASE_IFACE_ID +#define _BSPLINEBASE_IFACE_ID "$Id: BSpline.h 6353 2008-05-05 19:30:48Z martinc $" +#endif + +/** + * @file + * + * This file defines the BSpline library interface. + * + */ +template class BSpline; + +/* + * Opaque member structure to hide the matrix implementation. + */ +template struct BSplineBaseP; + +/** + * @class BSplineBase + * + * The base class for a spline object containing the nodes for a given + * domain, cutoff wavelength, and boundary condition. + + * To smooth a single curve, the BSpline interface contains a constructor + * which both sets up the domain and solves for the spline. Subsequent + * curves over the same domain can be created by apply()ing them to the + * BSpline object, where apply() is a BSplineBase method. [See apply().] + * New curves can also be smoothed within the same BSpline object by + * calling solve() with the new set of y values. [See BSpline.] A + * BSplineBase can be created on its own, in which case all of the + * computations dependent on the x values, boundary conditions, and cutoff + * wavelength have already been completed. + * + * The solution of the cubic b-spline is divided into two parts. The first + * is the setup of the domain given the x values, boundary conditions, and + * wavelength. The second is the solution of the spline for a set of y + * values corresponding to the x values in the domain. The first part is + * done in the creation of the BSplineBase object (or when calling the + * setDomain method). The second part is done when creating a BSpline + * object (or calling solve() on a BSpline object). + * + * A BSpline object can be created with either one of its constructors, or + * by calling apply() on an existing BSplineBase object. Once a spline has + * been solved, it can be evaluated at any x value. The following example + * creates a spline curve and evaluates it over the domain: + * +@verbatim + + vector x; + vector y; + { ... } + int bc = BSplineBase::BC_ZERO_SECOND; + BSpline::Debug = true; + BSpline spline (x.begin(), x.size(), y.begin(), wl, bc); + if (spline.ok()) + { + ostream_iterator of(cout, "\t "); + float xi = spline.Xmin(); + float xs = (spline.Xmax() - xi) / 2000.0; + for (; xi <= spline.Xmax(); xi += xs) + { + *of++ = spline.evaluate (xi); + } + } + +@endverbatim + * + * In the usual usage, the BSplineBase can compute a reasonable number of + * nodes for the spline, balancing between a few desirable factors. There + * needs to be at least 2 nodes per cutoff wavelength (preferably 4 or + * more) for the derivative constraint to reliably approximate a lo-pass + * filter. There should be at least 1 and preferably about 2 data points + * per node (measured just by their number and not by any check of the + * density of points across the domain). Lastly, of course, the fewer the + * nodes then the faster the computation of the spline. The computation of + * the number of nodes happens in the Setup() method during BSplineBase + * construction and when setDomain() is called. If the setup fails to find + * a desirable number of nodes, then the BSplineBase object will return + * false from ok(). + * + * The ok() method returns false when a BSplineBase or BSpline could not + * complete any operation successfully. In particular, as mentioned above, + * ok() will return false if some problem was detected with the domain + * values or if no reasonable number of nodes could be found for the given + * cutoff wavelength. Also, ok() on a BSpline object will return false if + * the matrix equation could not be solved, such as after BSpline + * construction or after a call to apply(). + * + * If letting Setup() determine the number of nodes is not acceptable, the + * constructors and setDomain() accept the parameter num_nodes. By + * default, num_nodes is passed as zero, forcing Setup() to calculate the + * number of nodes. However, if num_nodes is passed as 2 or greater, then + * Setup() will bypass its own algorithm and accept the given number of + * nodes instead. Obviously, it's up to the programmer to understand the + * affects of the number of nodes on the representation of the data and on + * the solution (or non-solution) of the spline. Remember to check the + * ok() method to detect when the spline solution has failed. + * + * The interface for the BSplineBase and BSpline templates is defined in + * the header file BSpline.h. The implementation is defined in BSpline.cpp. + * Source files which will instantiate the template should include the + * implementation file and @em not the interface. If the implementation + * for a specific type will be linked from elsewhere, such as a + * static library or Windows DLL, source files should only include the + * interface file. On Windows, applications should link with the import + * library BSpline.lib and make sure BSpline.dll is on the path. The DLL + * contains an implementation for BSpline and BSpline. + * For debugging, an application can include the implementation to get its + * own instantiation. + * + * The algorithm is based on the cubic spline described by Katsuyuki Ooyama + * in Montly Weather Review, Vol 115, October 1987. This implementation + * has benefited from comparisons with a previous FORTRAN implementation by + * James L. Franklin, NOAA/Hurricane Research Division. In particular, the + * algorithm in the Setup() method is based mostly on his implementation + * (VICSETUP). The Setup() method finds a suitable default for the number + * of nodes given a domain and cutoff frequency. This implementation + * adopts most of the same constraints, including a constraint that the + * cutoff wavelength not be greater than the span of the domain values: wl + * < max(x) - min(x). If this is not an acceptable constraint, then use the + * num_nodes parameter to specify the number of nodes explicitly. + * + * The cubic b-spline is formulated as the sum of some multiple of the + * basis function centered at each node in the domain. The number of nodes + * is determined by the desired cutoff wavelength and a desirable number of + * x values per node. The basis function is continuous and differentiable + * up to the second degree. A derivative constraint is included in the + * solution to achieve the effect of a low-pass frequency filter with the + * given cutoff wavelength. The derivative constraint can be disabled by + * specifying a wavelength value of zero, which reduces the analysis to a + * least squares fit to a cubic b-spline. The domain nodes, boundary + * constraints, and wavelength determine a linear system of equations, + * Qa=b, where a is the vector of basis function coefficients at each node. + * The coefficient vector is solved by first LU factoring along the + * diagonally banded matrix Q in BSplineBase. The BSpline object then + * computes the B vector for a set of y values and solves for the + * coefficient vector with the LU matrix. Only the diagonal bands are + * stored in memory and calculated during LU factoring and back + * substitution, and the basis function is evaluated as few times as + * possible in computing the diagonal matrix and B vector. + * + * @author Gary Granger (http://www.eol.ucar.edu/homes/granger) + * +@verbatim +Copyright (c) 1998-2009 +University Corporation for Atmospheric Research, UCAR +All rights reserved. +@endverbatim + **/ + +template +class BSplineBase +{ +public: + // Datum type + typedef T datum_type; + + /// Return a string describing the implementation version. + static const char *ImplVersion(); + + /// Return a string describing the interface version. + static const char *IfaceVersion(); + + /** + * Call this class method with a value greater than zero to enable + * debug messages, or with zero to disable messages. Calling with + * no arguments returns true if debugging enabled, else false. + */ + static bool Debug (int on = -1); + + /** + * Boundary condition types. + */ + enum BoundaryConditionTypes + { + /// Set the endpoints of the spline to zero. + BC_ZERO_ENDPOINTS = 0, + /// Set the first derivative of the spline to zero at the endpoints. + BC_ZERO_FIRST = 1, + /// Set the second derivative to zero. + BC_ZERO_SECOND = 2 + }; + +public: + + /** + * Construct a spline domain for the given set of x values, cutoff + * wavelength, and boundary condition type. The parameters are the + * same as for setDomain(). Call ok() to check whether domain + * setup succeeded after construction. + */ + BSplineBase (const T *x, int nx, + double wl, int bc_type = BC_ZERO_SECOND, + int num_nodes = 0); + + /// Copy constructor + BSplineBase (const BSplineBase &); + + /** + * Change the domain of this base. [If this is part of a BSpline + * object, this method {\em will not} change the existing curve or + * re-apply the smoothing to any set of y values.] + * + * The x values can be in any order, but they must be of sufficient + * density to support the requested cutoff wavelength. The setup of + * the domain may fail because of either inconsistency between the x + * density and the cutoff wavelength, or because the resulting matrix + * could not be factored. If setup fails, the method returns false. + * + * @param x The array of x values in the domain. + * @param nx The number of values in the @p x array. + * @param wl The cutoff wavelength, in the same units as the + * @p x values. A wavelength of zero disables + * the derivative constraint. + * @param bc_type The enumerated boundary condition type. If + * omitted it defaults to BC_ZERO_SECOND. + * @param num_nodes The number of nodes to use for the cubic b-spline. + * If less than 2 a reasonable number will be + * calculated automatically, if possible, taking + * into account the given cutoff wavelength. + * + * @see ok(). + */ + bool setDomain (const T *x, int nx, double wl, + int bc_type = BC_ZERO_SECOND, + int num_nodes = 0); + + /** + * Create a BSpline smoothed curve for the given set of NX y values. + * The returned object will need to be deleted by the caller. + * @param y The array of y values corresponding to each of the nX() + * x values in the domain. + * @see ok() + */ + BSpline *apply (const T *y); + + /** + * Return array of the node coordinates. Returns 0 if not ok(). The + * array of nodes returned by nodes() belongs to the object and should + * not be deleted; it will also be invalid if the object is destroyed. + */ + const T *nodes (int *nnodes); + + /** + * Return the number of nodes (one more than the number of intervals). + */ + int nNodes () { return M+1; } + + /** + * Number of original x values. + */ + int nX () { return NX; } + + /// Minimum x value found. + T Xmin () { return xmin; } + + /// Maximum x value found. + T Xmax () { return xmin + (M * DX); } + + /** + * Return the Alpha value for a given wavelength. Note that this + * depends on the current node interval length (DX). + */ + double Alpha (double wavelength); + + /** + * Return alpha currently in use by this domain. + */ + double Alpha () { return alpha; } + + /** + * Return the current state of the object, either ok or not ok. + * Use this method to test for valid state after construction or after + * a call to setDomain(). ok() will return false if either fail, such + * as when an appropriate number of nodes and node interval cannot be + * found for a given wavelength, or when the linear equation for the + * coefficients cannot be solved. + */ + bool ok () { return OK; } + + virtual ~BSplineBase(); + +protected: + + typedef BSplineBaseP Base; + + // Provided + double waveLength; // Cutoff wavelength (l sub c) + int NX; + int K; // Degree of derivative constraint (currently fixed at 2) + int BC; // Boundary conditions type (0,1,2) + + // Derived + T xmax; + T xmin; + int M; // Number of intervals (M+1 nodes) + double DX; // Interval length in same units as X + double alpha; + bool OK; + Base *base; // Hide more complicated state members + // from the public interface. + + bool Setup (int num_nodes = 0); + void calculateQ (); + double qDelta (int m1, int m2); + double Beta (int m); + void addP (); + bool factor (); + double Basis (int m, T x); + double DBasis (int m, T x); + + static const double BoundaryConditions[3][4]; + static const double BS_PI; + + double Ratiod (int&, double &, double &); +}; + + + +/* + * Original BSpline.h start here + */ + + + +template struct BSplineP; + + +/** + * Used to evaluate a BSpline. + * Inherits the BSplineBase domain information and interface and adds + * smoothing. See the BSplineBase documentation for a summary of the + * BSpline interface. + */ +template +class BSpline : public BSplineBase +{ +public: + /** + * Create a single spline with the parameters required to set up + * the domain and subsequently smooth the given set of y values. + * The y values must correspond to each of the values in the x array. + * If either the domain setup fails or the spline cannot be solved, + * the state will be set to not ok. + * + * @see ok(). + * + * @param x The array of x values in the domain. + * @param nx The number of values in the @p x array. + * @param y The array of y values corresponding to each of the + * nX() x values in the domain. + * @param wl The cutoff wavelength, in the same units as the + * @p x values. A wavelength of zero disables + * the derivative constraint. + * @param bc_type The enumerated boundary condition type. If + * omitted it defaults to BC_ZERO_SECOND. + * @param num_nodes The number of nodes to use for the cubic b-spline. + * If less than 2 a "reasonable" number will be + * calculated automatically, taking into account + * the given cutoff wavelength. + */ + BSpline (const T *x, int nx, /* independent variable */ + const T *y, /* dependent values @ ea X */ + double wl, /* cutoff wavelength */ + int bc_type = BSplineBase::BC_ZERO_SECOND, + int num_nodes = 0); + + /** + * A BSpline curve can be derived from a separate @p base and a set + * of data points @p y over that base. + */ + BSpline (BSplineBase &base, const T *y); + + /** + * Solve the spline curve for a new set of y values. Returns false + * if the solution fails. + * + * @param y The array of y values corresponding to each of the nX() + * x values in the domain. + */ + bool solve (const T *y); + + /** + * Return the evaluation of the smoothed curve + * at a particular @p x value. If current state is not ok(), returns 0. + */ + T evaluate (T x); + + /** + * Return the first derivative of the spline curve at the given @p x. + * Returns zero if the current state is not ok(). + */ + T slope (T x); + + /** + * Return the @p n-th basis coefficient, from 0 to M. If the current + * state is not ok(), or @p n is out of range, the method returns zero. + */ + T coefficient (int n); + + virtual ~BSpline(); + + using BSplineBase::Debug; + using BSplineBase::Basis; + using BSplineBase::DBasis; + +protected: + + using BSplineBase::OK; + using BSplineBase::M; + using BSplineBase::NX; + using BSplineBase::DX; + using BSplineBase::base; + using BSplineBase::xmin; + using BSplineBase::xmax; + + // Our hidden state structure + BSplineP *s; + T mean; // Fit without mean and add it in later + +}; + +#endif diff --git a/xs/src/BSpline/BandedMatrix.h b/xs/src/BSpline/BandedMatrix.h new file mode 100644 index 000000000..bbe526a8b --- /dev/null +++ b/xs/src/BSpline/BandedMatrix.h @@ -0,0 +1,375 @@ +/* -*- mode: c++; c-basic-offset: 4; -*- */ +/************************************************************************ + * Copyright 2009 University Corporation for Atmospheric Research. + * All rights reserved. + * + * Use of this code is subject to the standard BSD license: + * + * http://www.opensource.org/licenses/bsd-license.html + * + * See the COPYRIGHT file in the source distribution for the license text, + * or see this web page: + * + * http://www.eol.ucar.edu/homes/granger/bspline/doc/ + * + *************************************************************************/ + +/** + * Template for a diagonally banded matrix. + **/ +#ifndef _BANDEDMATRIX_ID +#define _BANDEDMATRIX_ID "$Id$" + +#include +#include +#include + +template class BandedMatrixRow; + + +template class BandedMatrix +{ +public: + typedef unsigned int size_type; + typedef T element_type; + + // Create a banded matrix with the same number of bands above and below + // the diagonal. + BandedMatrix (int N_ = 1, int nbands_off_diagonal = 0) : bands(0) + { + if (! setup (N_, nbands_off_diagonal)) + setup (); + } + + // Create a banded matrix by naming the first and last non-zero bands, + // where the diagonal is at zero, and bands below the diagonal are + // negative, bands above the diagonal are positive. + BandedMatrix (int N_, int first, int last) : bands(0) + { + if (! setup (N_, first, last)) + setup (); + } + + // Copy constructor + BandedMatrix (const BandedMatrix &b) : bands(0) + { + Copy (*this, b); + } + + inline bool setup (int N_ = 1, int noff = 0) + { + return setup (N_, -noff, noff); + } + + bool setup (int N_, int first, int last) + { + // Check our limits first and make sure they make sense. + // Don't change anything until we know it will work. + if (first > last || N_ <= 0) + return false; + + // Need at least as many N_ as columns and as rows in the bands. + if (N_ < abs(first) || N_ < abs(last)) + return false; + + top = last; + bot = first; + N = N_; + out_of_bounds = T(); + + // Finally setup the diagonal vectors + nbands = last - first + 1; + if (bands) delete[] bands; + bands = new std::vector[nbands]; + int i; + for (i = 0; i < nbands; ++i) + { + // The length of each array varies with its distance from the + // diagonal + int len = N - (abs(bot + i)); + bands[i].clear (); + bands[i].resize (len); + } + return true; + } + + BandedMatrix & operator= (const BandedMatrix &b) + { + return Copy (*this, b); + } + + BandedMatrix & operator= (const T &e) + { + int i; + for (i = 0; i < nbands; ++i) + { + std::fill_n (bands[i].begin(), bands[i].size(), e); + } + out_of_bounds = e; + return (*this); + } + + ~BandedMatrix () + { + if (bands) + delete[] bands; + } + +private: + // Return false if coordinates are out of bounds + inline bool check_bounds (int i, int j, int &v, int &e) const + { + v = (j - i) - bot; + e = (i >= j) ? j : i; + return !(v < 0 || v >= nbands || + e < 0 || (unsigned int)e >= bands[v].size()); + } + + static BandedMatrix & Copy (BandedMatrix &a, const BandedMatrix &b) + { + if (a.bands) delete[] a.bands; + a.top = b.top; + a.bot = b.bot; + a.N = b.N; + a.out_of_bounds = b.out_of_bounds; + a.nbands = a.top - a.bot + 1; + a.bands = new std::vector[a.nbands]; + int i; + for (i = 0; i < a.nbands; ++i) + { + a.bands[i] = b.bands[i]; + } + return a; + } + +public: + T &element (int i, int j) + { + int v, e; + if (check_bounds(i, j, v, e)) + return (bands[v][e]); + else + return out_of_bounds; + } + + const T &element (int i, int j) const + { + int v, e; + if (check_bounds(i, j, v, e)) + return (bands[v][e]); + else + return out_of_bounds; + } + + inline T & operator() (int i, int j) + { + return element (i-1,j-1); + } + + inline const T & operator() (int i, int j) const + { + return element (i-1,j-1); + } + + size_type num_rows() const { return N; } + + size_type num_cols() const { return N; } + + const BandedMatrixRow operator[] (int row) const + { + return BandedMatrixRow(*this, row); + } + + BandedMatrixRow operator[] (int row) + { + return BandedMatrixRow(*this, row); + } + + +private: + + int top; + int bot; + int nbands; + std::vector *bands; + int N; + T out_of_bounds; + +}; + + +template +std::ostream &operator<< (std::ostream &out, const BandedMatrix &m) +{ + unsigned int i, j; + for (i = 0; i < m.num_rows(); ++i) + { + for (j = 0; j < m.num_cols(); ++j) + { + out << m.element (i, j) << " "; + } + out << std::endl; + } + return out; +} + + + +/* + * Helper class for the intermediate in the [][] operation. + */ +template class BandedMatrixRow +{ +public: + BandedMatrixRow (const BandedMatrix &_m, int _row) : bm(_m), i(_row) + { } + + BandedMatrixRow (BandedMatrix &_m, int _row) : bm(_m), i(_row) + { } + + ~BandedMatrixRow () {} + + typename BandedMatrix::element_type & operator[] (int j) + { + return const_cast &>(bm).element (i, j); + } + + const typename BandedMatrix::element_type & operator[] (int j) const + { + return bm.element (i, j); + } + +private: + const BandedMatrix &bm; + int i; +}; + + +/* + * Vector multiplication + */ + +template +Vector operator* (const Matrix &m, const Vector &v) +{ + typename Matrix::size_type M = m.num_rows(); + typename Matrix::size_type N = m.num_cols(); + + assert (N <= v.size()); + //if (M > v.size()) + // return Vector(); + + Vector r(N); + for (unsigned int i = 0; i < M; ++i) + { + typename Matrix::element_type sum = 0; + for (unsigned int j = 0; j < N; ++j) + { + sum += m[i][j] * v[j]; + } + r[i] = sum; + } + return r; +} + + + +/* + * LU factor a diagonally banded matrix using Crout's algorithm, but + * limiting the trailing sub-matrix multiplication to the non-zero + * elements in the diagonal bands. Return nonzero if a problem occurs. + */ +template +int LU_factor_banded (MT &A, unsigned int bands) +{ + typename MT::size_type M = A.num_rows(); + typename MT::size_type N = A.num_cols(); + if (M != N) + return 1; + + typename MT::size_type i,j,k; + typename MT::element_type sum; + + for (j = 1; j <= N; ++j) + { + // Check for zero pivot + if ( A(j,j) == 0 ) + return 1; + + // Calculate rows above and on diagonal. A(1,j) remains as A(1,j). + for (i = (j > bands) ? j-bands : 1; i <= j; ++i) + { + sum = 0; + for (k = (j > bands) ? j-bands : 1; k < i; ++k) + { + sum += A(i,k)*A(k,j); + } + A(i,j) -= sum; + } + + // Calculate rows below the diagonal. + for (i = j+1; (i <= M) && (i <= j+bands); ++i) + { + sum = 0; + for (k = (i > bands) ? i-bands : 1; k < j; ++k) + { + sum += A(i,k)*A(k,j); + } + A(i,j) = (A(i,j) - sum) / A(j,j); + } + } + + return 0; +} + + + +/* + * Solving (LU)x = B. First forward substitute to solve for y, Ly = B. + * Then backwards substitute to find x, Ux = y. Return nonzero if a + * problem occurs. Limit the substitution sums to the elements on the + * bands above and below the diagonal. + */ +template +int LU_solve_banded(const MT &A, Vector &b, unsigned int bands) +{ + typename MT::size_type i,j; + typename MT::size_type M = A.num_rows(); + typename MT::size_type N = A.num_cols(); + typename MT::element_type sum; + + if (M != N || M == 0) + return 1; + + // Forward substitution to find y. The diagonals of the lower + // triangular matrix are taken to be 1. + for (i = 2; i <= M; ++i) + { + sum = b[i-1]; + for (j = (i > bands) ? i-bands : 1; j < i; ++j) + { + sum -= A(i,j)*b[j-1]; + } + b[i-1] = sum; + } + + // Now for the backward substitution + b[M-1] /= A(M,M); + for (i = M-1; i >= 1; --i) + { + if (A(i,i) == 0) // oops! + return 1; + sum = b[i-1]; + for (j = i+1; (j <= N) && (j <= i+bands); ++j) + { + sum -= A(i,j)*b[j-1]; + } + b[i-1] = sum / A(i,i); + } + + return 0; +} + + +#endif /* _BANDEDMATRIX_ID */ + diff --git a/xs/src/BSpline/COPYRIGHT b/xs/src/BSpline/COPYRIGHT new file mode 100644 index 000000000..50c77bb13 --- /dev/null +++ b/xs/src/BSpline/COPYRIGHT @@ -0,0 +1,44 @@ +Copyright (c) 1998-2009,2015 +University Corporation for Atmospheric Research, UCAR +All rights reserved. + +This software is licensed with the standard BSD license: + + http://www.opensource.org/licenses/bsd-license.html + +When citing this software, here is a suggested reference: + + This software is written by Gary Granger of the National Center for + Atmospheric Research (NCAR), sponsored by the National Science Foundation + (NSF). The algorithm is based on the cubic spline described by Katsuyuki + Ooyama in Montly Weather Review, Vol 115, October 1987. This + implementation has benefited from comparisons with a previous FORTRAN + implementation by James L. Franklin, NOAA/Hurricane Research Division. + +The text of the license is reproduced below: + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + * Neither the name of the UCAR nor the names of its contributors may be + used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE.