From 6353a5ac078ac17e6b244150f7090aa7491962a5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 7 Jul 2024 10:57:28 -0700 Subject: [PATCH] Removed experimental extensions from previous mixed bug-fix and feature merge --- expui/CMakeLists.txt | 2 +- expui/SGSmooth.H | 13 - expui/SGSmooth.cc | 549 ------------------------------------------ expui/expMSSA.cc | 259 -------------------- pyEXP/MSSAWrappers.cc | 43 ---- 5 files changed, 1 insertion(+), 865 deletions(-) delete mode 100644 expui/SGSmooth.H delete mode 100644 expui/SGSmooth.cc diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index d506bea8..8690817f 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -34,7 +34,7 @@ endif() set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc - Koopman.cc BiorthBess.cc SGSmooth.cc) + Koopman.cc BiorthBess.cc) add_library(expui ${expui_SOURCES}) set_target_properties(expui PROPERTIES OUTPUT_NAME expui) target_include_directories(expui PUBLIC ${common_INCLUDE}) diff --git a/expui/SGSmooth.H b/expui/SGSmooth.H deleted file mode 100644 index 6615bffa..00000000 --- a/expui/SGSmooth.H +++ /dev/null @@ -1,13 +0,0 @@ -#ifndef __SGSMOOTH_HPP__ -#define __SGSMOOTH_HPP__ - -#include - -//! savitzky golay smoothing. -std::vector sg_smooth(const std::vector &v, const int w, const int deg); - -//! numerical derivative based on savitzky golay smoothing. -std::vector sg_derivative(const std::vector &v, const int w, - const int deg, const double h=1.0); - -#endif // __SGSMOOTH_HPP__ diff --git a/expui/SGSmooth.cc b/expui/SGSmooth.cc deleted file mode 100644 index 1b5150b5..00000000 --- a/expui/SGSmooth.cc +++ /dev/null @@ -1,549 +0,0 @@ -//! -// Sliding window signal processing (and linear algebra toolkit). -// -// supported operations: -//
    -//
  • Savitzky-Golay smoothing. -//
  • computing a numerical derivative based of Savitzky-Golay smoothing. -//
  • required linear algebra support for SG smoothing using STL based -// vector/matrix classes -//
-// -// \brief Linear Algebra "Toolkit". -// -// modified by Rob Patro, 2016 -// modified by MDW, 2024 - -// system headers -#include -#include -#include // for size_t -#include // for fabs -#include - -//! default convergence -static const double TINY_FLOAT = 1.0e-300; - -//! comfortable array of doubles -using float_vect = std::vector; -//! comfortable array of ints; -using int_vect = std::vector; - -/*! matrix class. - * - * This is a matrix class derived from a vector of float_vects. Note that - * the matrix elements indexed [row][column] with indices starting at 0 (c - * style). Also note that because of its design looping through rows should - * be faster than looping through columns. - * - * \brief two dimensional floating point array - */ -class float_mat : public std::vector { -private: - //! disable the default constructor - explicit float_mat() {}; - //! disable assignment operator until it is implemented. - float_mat &operator =(const float_mat &) { return *this; }; -public: - //! constructor with sizes - float_mat(const size_t rows, const size_t cols, const double def=0.0); - //! copy constructor for matrix - float_mat(const float_mat &m); - //! copy constructor for vector - float_mat(const float_vect &v); - - //! use default destructor - // ~float_mat() {}; - - //! get size - size_t nr_rows(void) const { return size(); }; - //! get size - size_t nr_cols(void) const { return front().size(); }; -}; - - - -// constructor with sizes -float_mat::float_mat(const size_t rows,const size_t cols,const double defval) - : std::vector(rows) { - int i; - for (i = 0; i < rows; ++i) { - (*this)[i].resize(cols, defval); - } - if ((rows < 1) || (cols < 1)) { - std::ostringstream msg; - msg << "cannot build matrix with " << rows << " rows and " << cols - << "columns"; - throw std::runtime_error(msg.str()); - } -} - -// copy constructor for matrix -float_mat::float_mat(const float_mat &m) : std::vector(m.size()) { - - float_mat::iterator inew = begin(); - float_mat::const_iterator iold = m.begin(); - for (/* empty */; iold < m.end(); ++inew, ++iold) { - const size_t oldsz = iold->size(); - inew->resize(oldsz); - const float_vect oldvec(*iold); - *inew = oldvec; - } -} - -// copy constructor for vector -float_mat::float_mat(const float_vect &v) - : std::vector(1) { - - const size_t oldsz = v.size(); - front().resize(oldsz); - front() = v; -} - -////////////////////// -// Helper functions // -////////////////////// - -//! permute() orders the rows of A to match the integers in the index array. -void permute(float_mat &A, int_vect &idx) -{ - int_vect i(idx.size()); - int j,k; - - for (j = 0; j < A.nr_rows(); ++j) { - i[j] = j; - } - - // loop over permuted indices - for (j = 0; j < A.nr_rows(); ++j) { - if (i[j] != idx[j]) { - - // search only the remaining indices - for (k = j+1; k < A.nr_rows(); ++k) { - if (i[k] ==idx[j]) { - std::swap(A[j],A[k]); // swap the rows and - i[k] = i[j]; // the elements of - i[j] = idx[j]; // the ordered index. - break; // next j - } - } - } - } -} - -/*! \brief Implicit partial pivoting. - * - * The function looks for pivot element only in rows below the current - * element, A[idx[row]][column], then swaps that row with the current one in - * the index map. The algorithm is for implicit pivoting (i.e., the pivot is - * chosen as if the max coefficient in each row is set to 1) based on the - * scaling information in the vector scale. The map of swapped indices is - * recorded in swp. The return value is +1 or -1 depending on whether the - * number of row swaps was even or odd respectively. */ -static int partial_pivot(float_mat &A, const size_t row, const size_t col, - float_vect &scale, int_vect &idx, double tol) -{ - if (tol <= 0.0) - tol = TINY_FLOAT; - - int swapNum = 1; - - // default pivot is the current position, [row,col] - int pivot = row; - double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]]; - - // loop over possible pivots below current - int j; - for (j = row + 1; j < A.nr_rows(); ++j) { - - const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]]; - - // if this elem is larger, then it becomes the pivot - if (tmp > piv_elem) { - pivot = j; - piv_elem = tmp; - } - } - -#if 0 - if (piv_elem < tol) { - throw std::runtime_error("partial_pivot(): Zero pivot encountered."); -#endif - - if(pivot > row) { // bring the pivot to the diagonal - j = idx[row]; // reorder swap array - idx[row] = idx[pivot]; - idx[pivot] = j; - swapNum = -swapNum; // keeping track of odd or even swap - } - return swapNum; -} - -/*! \brief Perform backward substitution. - * - * Solves the system of equations A*b=a, ASSUMING that A is upper - * triangular. If diag==1, then the diagonal elements are additionally - * assumed to be 1. Note that the lower triangular elements are never - * checked, so this function is valid to use after a LU-decomposition in - * place. A is not modified, and the solution, b, is returned in a. */ -static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false) -{ - int r,c,k; - - for (r = (A.nr_rows() - 1); r >= 0; --r) { - for (c = (A.nr_cols() - 1); c > r; --c) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] -= A[r][c] * a[c][k]; - } - } - if(!diag) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] /= A[r][r]; - } - } - } -} - -/*! \brief Perform forward substitution. - * - * Solves the system of equations A*b=a, ASSUMING that A is lower - * triangular. If diag==1, then the diagonal elements are additionally - * assumed to be 1. Note that the upper triangular elements are never - * checked, so this function is valid to use after a LU-decomposition in - * place. A is not modified, and the solution, b, is returned in a. */ -static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true) -{ - int r,k,c; - for (r = 0;r < A.nr_rows(); ++r) { - for(c = 0; c < r; ++c) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] -= A[r][c] * a[c][k]; - } - } - if(!diag) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] /= A[r][r]; - } - } - } -} - -/*! \brief Performs LU factorization in place. - * - * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of - * swapped indeces is recorded in idx. The return value is +1 or -1 - * depending on whether the number of row swaps was even or odd - * respectively. idx must be preinitialized to a valid set of indices - * (e.g., {1,2, ... ,A.nr_rows()}). */ -static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT) -{ - if ( tol <= 0.0) - tol = TINY_FLOAT; - - if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) { - throw std::runtime_error("lu_factorize(): cannot handle empty " - "or nonsquare matrices."); - } - - float_vect scale(A.nr_rows()); // implicit pivot scaling - int i,j; - for (i = 0; i < A.nr_rows(); ++i) { - double maxval = 0.0; - for (j = 0; j < A.nr_cols(); ++j) { - if (fabs(A[i][j]) > maxval) - maxval = fabs(A[i][j]); - } - if (maxval == 0.0) { - throw std::runtime_error("lu_factorize(): zero pivot found."); - } - scale[i] = 1.0 / maxval; - } - - int swapNum = 1; - int c,r; - for (c = 0; c < A.nr_cols() ; ++c) { // loop over columns - swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal - for(r = 0; r < A.nr_rows(); ++r) { // loop over rows - int lim = (r < c) ? r : c; - for (j = 0; j < lim; ++j) { - A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c]; - } - if (r > c) - A[idx[r]][c] /= A[idx[c]][c]; - } - } - permute(A,idx); - return swapNum; -} - -/*! \brief Solve a system of linear equations. - * Solves the inhomogeneous matrix problem with lu-decomposition. Note that - * inversion may be accomplished by setting a to the identity_matrix. */ -static float_mat lin_solve(const float_mat &A, const float_mat &a, - double tol=TINY_FLOAT) -{ - float_mat B(A); - float_mat b(a); - int_vect idx(B.nr_rows()); - int j; - - for (j = 0; j < B.nr_rows(); ++j) { - idx[j] = j; // init row swap label array - } - lu_factorize(B,idx,tol); // get the lu-decomp. - permute(b,idx); // sort the inhomogeneity to match the lu-decomp - lu_forwsubst(B,b); // solve the forward problem - lu_backsubst(B,b); // solve the backward problem - return b; -} - -/////////////////////// -// related functions // -/////////////////////// - -//! Returns the inverse of a matrix using LU-decomposition. -static float_mat invert(const float_mat &A) -{ - const int n = A.size(); - float_mat E(n, n, 0.0); - float_mat B(A); - int i; - - for (i = 0; i < n; ++i) { - E[i][i] = 1.0; - } - - return lin_solve(B, E); -} - -//! returns the transposed matrix. -static float_mat transpose(const float_mat &a) -{ - float_mat res(a.nr_cols(), a.nr_rows()); - int i,j; - - for (i = 0; i < a.nr_rows(); ++i) { - for (j = 0; j < a.nr_cols(); ++j) { - res[j][i] = a[i][j]; - } - } - return res; -} - -//! matrix multiplication. -float_mat operator *(const float_mat &a, const float_mat &b) -{ - float_mat res(a.nr_rows(), b.nr_cols()); - if (a.nr_cols() != b.nr_rows()) { - throw std::runtime_error("incompatible matrices in multiplication"); - } - - int i,j,k; - - for (i = 0; i < a.nr_rows(); ++i) { - for (j = 0; j < b.nr_cols(); ++j) { - double sum(0.0); - for (k = 0; k < a.nr_cols(); ++k) { - sum += a[i][k] * b[k][j]; - } - res[i][j] = sum; - } - } - return res; -} - - -//! calculate savitzky golay coefficients. -static float_vect sg_coeff(const float_vect &b, const size_t deg) -{ - const size_t rows(b.size()); - const size_t cols(deg + 1); - float_mat A(rows, cols); - float_vect res(rows); - - // generate input matrix for least squares fit - int i,j; - for (i = 0; i < rows; ++i) { - for (j = 0; j < cols; ++j) { - A[i][j] = pow(double(i), double(j)); - } - } - - float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); - - for (i = 0; i < b.size(); ++i) { - res[i] = c[0][0]; - for (j = 1; j <= deg; ++j) { - res[i] += c[j][0] * pow(double(i), double(j)); - } - } - return res; -} - -/*! \brief savitzky golay smoothing. - * - * This method means fitting a polynome of degree 'deg' to a sliding window - * of width 2w+1 throughout the data. The needed coefficients are - * generated dynamically by doing a least squares fit on a "symmetric" unit - * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome - * yields the sg-coefficients. at the border non symmectric vectors b are - * used. */ -float_vect sg_smooth(const float_vect &v, const int width, const int deg) -{ - float_vect res(v.size(), 0.0); - if ((width < 1) || (deg < 0) || (v.size() < (2 * width + 2))) { - throw std::runtime_error("sgsmooth: parameter error."); - } - - const int window = 2 * width + 1; - const int endidx = v.size() - 1; - - // do a regular sliding window average - int i,j; - if (deg == 0) { - // handle border cases first because we need different coefficients -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i < width; ++i) { - const double scale = 1.0/double(i+1); - const float_vect c1(width, scale); - for (j = 0; j <= i; ++j) { - res[i] += c1[j] * v[j]; - res[endidx - i] += c1[j] * v[endidx - j]; - } - } - - // now loop over rest of data. reusing the "symmetric" coefficients. - const double scale = 1.0/double(window); - const float_vect c2(window, scale); -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i <= (v.size() - window); ++i) { - for (j = 0; j < window; ++j) { - res[i + width] += c2[j] * v[i + j]; - } - } - return res; - } - - // handle border cases first because we need different coefficients -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i < width; ++i) { - float_vect b1(window, 0.0); - b1[i] = 1.0; - - const float_vect c1(sg_coeff(b1, deg)); - for (j = 0; j < window; ++j) { - res[i] += c1[j] * v[j]; - res[endidx - i] += c1[j] * v[endidx - j]; - } - } - - // now loop over rest of data. reusing the "symmetric" coefficients. - float_vect b2(window, 0.0); - b2[width] = 1.0; - const float_vect c2(sg_coeff(b2, deg)); - -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i <= (v.size() - window); ++i) { - for (j = 0; j < window; ++j) { - res[i + width] += c2[j] * v[i + j]; - } - } - return res; -} - -/*! least squares fit a polynome of degree 'deg' to data in 'b'. - * then calculate the first derivative and return it. */ -static float_vect lsqr_fprime(const float_vect &b, const int deg) -{ - const int rows(b.size()); - const int cols(deg + 1); - float_mat A(rows, cols); - float_vect res(rows); - - // generate input matrix for least squares fit - int i,j; - for (i = 0; i < rows; ++i) { - for (j = 0; j < cols; ++j) { - A[i][j] = pow(double(i), double(j)); - } - } - - float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); - - for (i = 0; i < b.size(); ++i) { - res[i] = c[1][0]; - for (j = 1; j < deg; ++j) { - res[i] += c[j + 1][0] * double(j+1) - * pow(double(i), double(j)); - } - } - return res; -} - -/*! \brief savitzky golay smoothed numerical derivative. - * - * This method means fitting a polynome of degree 'deg' to a sliding window - * of width 2w+1 throughout the data. - * - * In contrast to the sg_smooth function we do a brute force attempt by - * always fitting the data to a polynome of degree 'deg' and using the - * result. */ -float_vect sg_derivative(const float_vect &v, const int width, - const int deg, const double h) -{ - float_vect res(v.size(), 0.0); - if ((width < 1) || (deg < 1) || (v.size() < (2 * width + 2))) { - throw std::runtime_error("sgsderiv: parameter error"); - } - - const int window = 2 * width + 1; - - // handle border cases first because we do not repeat the fit - // lower part - float_vect b(window, 0.0); - int i,j; - - for (i = 0; i < window; ++i) { - b[i] = v[i] / h; - } - const float_vect c(lsqr_fprime(b, deg)); - for (j = 0; j <= width; ++j) { - res[j] = c[j]; - } - // upper part. direction of fit is reversed - for (i = 0; i < window; ++i) { - b[i] = v[v.size() - 1 - i] / h; - } - const float_vect d(lsqr_fprime(b, deg)); - for (i = 0; i <= width; ++i) { - res[v.size() - 1 - i] = -d[i]; - } - - // now loop over rest of data. wasting a lot of least squares calcs - // since we only use the middle value. -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 1; i < (v.size() - window); ++i) { - for (j = 0; j < window; ++j) { - b[j] = v[i + j] / h; - } - res[i + width] = lsqr_fprime(b, deg)[width]; - } - return res; -} - -// Local Variables: -// mode: c++ -// c-basic-offset: 4 -// fill-column: 76 -// indent-tabs-mode: nil -// End: diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index 7a19c6f6..b136d2e7 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -1554,265 +1554,6 @@ namespace MSSA { } - std::tuple - expMSSA::getKoopmanModes(double tol, int D, bool debug) - { - bool use_fullKh = true; // Use the non-reduced computation of - // Koopman/eDMD - // Number of channels - // - nkeys = mean.size(); - - // Make sure parameters are sane - // - if (numW<=0) numW = numT/2; - if (numW > numT/2) numW = numT/2; - - numK = numT - numW + 1; - - Eigen::VectorXd S1; - Eigen::MatrixXd Y1; - Eigen::MatrixXd V1; - Eigen::MatrixXd VT1; - Eigen::MatrixXd VT2; - - // Make a new trajetory matrix with smoothing - // - Y1.resize(numK, numW*nkeys + D*(nkeys-1)); - - size_t n=0, offset=0; - for (auto k : mean) { - for (int i=0; i 0) { - // Back blending - for (int j=0; j(D-j)/D; - } - } - // Main series - for (int j=0; j(D-j)/D; - } - } - } - offset += numW + D; - n++; - } - - double Scale = Y1.norm(); - - // auto YY = Y1/Scale; - auto YY = Y1; - - // Use one of the built-in Eigen3 algorithms - // - /* - if (params["Jacobi"]) { - // -->Using Jacobi - Eigen::JacobiSVD - svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); - S1 = svd.singularValues(); - V1 = svd.matrixV(); - } else if (params["BDCSVD"]) { - */ - // -->Using BDC - Eigen::BDCSVD - svd(YY, Eigen::ComputeFullU | Eigen::ComputeFullV); - // svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); - S1 = svd.singularValues(); - V1 = svd.matrixV(); - /* - } else { - // -->Use Random approximation algorithm from Halko, Martinsson, - // and Tropp - int srank = std::min(YY.cols(), YY.rows()); - RedSVD::RedSVD svd(YY, srank); - S1 = svd.singularValues(); - V1 = svd.matrixV(); - } - */ - - std::cout << "shape V1 = " << V1.rows() << " x " - << V1.cols() << std::endl; - - std::cout << "shape Y1 = " << Y1.rows() << " x " - << Y1.cols() << std::endl; - - int lags = V1.rows(); - int rank = V1.cols(); - - std::ofstream out; - if (debug) out.open("debug.txt"); - - if (out) out << "rank=" << rank << " lags=" << lags << std::endl; - - VT1.resize(rank, lags-1); - VT2.resize(rank, lags-1); - - for (int j=0; j uu; - for (int i=0; iUsing Jacobi - Eigen::JacobiSVD - // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); - svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); - SS = svd.singularValues(); - UU = svd.matrixU(); - VV = svd.matrixV(); - } else if (params["BDCSVD"]) { - */ - { - // -->Using BDC - Eigen::BDCSVD - // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); - svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); - SS = svd.singularValues(); - UU = svd.matrixU(); - VV = svd.matrixV(); - } - /* - } else { - // -->Use Random approximation algorithm from Halko, Martinsson, - // and Tropp - // RedSVD::RedSVD svd(VT1, std::min(rank, numK-1)); - RedSVD::RedSVD svd(VT1, std::max(VT1.rows(), VT2.cols())); - SS = svd.singularValues(); - UU = svd.matrixU(); - VV = svd.matrixV(); - } - */ - - if (out) out << "Singular values" << std::endl << SS << std::endl; - - // Compute inverse - for (int i=0; i tol) SS(i) = 1.0/SS(i); - else SS(i) = 0.0; - } - - // Compute full Koopman operator - if (use_fullKh) { - - Eigen::MatrixXd DD(VV.cols(), UU.cols()); - DD.setZero(); - for (int i=0; i es(AT); - - L = es.eigenvalues(); - Phi = es.eigenvectors(); - - if (out) { - out << std::endl << "Eigenvalues" << std::endl << L << std::endl - << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; - } - - } - // Compute the reduced Koopman operator - else { - - Eigen::MatrixXd AT = UU.transpose() * (VT2 * VV) * SS.asDiagonal(); - - // Compute spectrum - Eigen::EigenSolver es(AT, true); - - L = es.eigenvalues(); - auto W = es.eigenvectors(); - - // Compute the EDMD modes - // - Eigen::VectorXcd Linv(L); - for (int i=0; i tol) Linv(i) = 1.0/Linv(i); - else Linv(i) = 0.0; - } - - Phi = VT2 * VV * SS.asDiagonal() * W * Linv.asDiagonal(); - - if (out) { - out << std::endl << "Eigenvalues" << std::endl << L << std::endl - << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; - } - } - - // Cache window size - // - window = D; - - return {L, Phi}; - } - - std::map - expMSSA::getReconstructedKoopman(int mode) - { - // Copy the original map for return - // - auto newdata = data; - - size_t n=0, offset=0; - - for (auto u : mean) { - - double disp = totVar; - if (type == TrendType::totPow) disp = totPow; - if (disp==0.0) disp = var[u.first]; - - std::complex phase = 1.0; - for (int i=0; i