diff --git a/fdaPDE/core b/fdaPDE/core index bb6063ad..a52aad88 160000 --- a/fdaPDE/core +++ b/fdaPDE/core @@ -1 +1 @@ -Subproject commit bb6063ad50b60472f79d6c0f892a652039922d5c +Subproject commit a52aad88444eabf31b2d30623a14bca057e11a48 diff --git a/fdaPDE/models/functional/center.h b/fdaPDE/models/functional/center.h index 17394327..bd4afdf0 100644 --- a/fdaPDE/models/functional/center.h +++ b/fdaPDE/models/functional/center.h @@ -27,6 +27,10 @@ struct CenterReturnType { DMatrix mean; // the mean field \mu }; + // need to handle case with missing data + // introudce a mean function which just produces the (smooth) mean field (no centering) + // then we can speak about weighted average field + // functional centering of the data matrix X template CenterReturnType center(const DMatrix& X, SmootherType_&& smoother, CalibrationType_&& calibration) { @@ -34,7 +38,12 @@ CenterReturnType center(const DMatrix& X, SmootherType_&& smoother, Cali df.insert(OBSERVATIONS_BLK, X.colwise().sum().transpose() / X.rows()); smoother.set_data(df); // set optimal lambda according to choosen calibration strategy - smoother.set_lambda(calibration.fit(smoother)); + if constexpr (std::is_base_of, CalibrationType_>::value) { // fixed \lambda + fdapde_assert(calibration.cols() == 1 && calibration.rows() == SmootherType_::n_lambda); + smoother.set_lambda(calibration); // center with fixed \lambda + } else { + smoother.set_lambda(calibration.fit(smoother)); + } smoother.init(); smoother.solve(); // compute mean matrix and return diff --git a/fdaPDE/models/functional/fpca.h b/fdaPDE/models/functional/fpca.h index 3e536eb3..b55a0228 100644 --- a/fdaPDE/models/functional/fpca.h +++ b/fdaPDE/models/functional/fpca.h @@ -27,172 +27,53 @@ using fdapde::calibration::Calibration; #include "functional_base.h" #include "power_iteration.h" -#include "rsvd.h" +#include "regularized_svd.h" namespace fdapde { namespace models { -// FPCA (Functional Principal Components Analysis) model signature -template class FPCA; - -// implementation of FPCA for sequential approach, see e.g. -// Lila, E., Aston, J.A.D., Sangalli, L.M. (2016), Smooth Principal Component Analysis over two-dimensional manifolds -// with an application to Neuroimaging, Annals of Applied Statistics, 10 (4), 1854-1879. +// Functional Principal Components Analysis template -class FPCA : - public FunctionalBase, RegularizationType_> { - public: - using RegularizationType = std::decay_t; - using This = FPCA; - using Base = FunctionalBase; - IMPORT_MODEL_SYMBOLS; - using Base::lambda; // vector of smoothing parameters - using Base::n_basis; // number of basis functions - using Base::n_stat_units; // number of statistical units - using Base::X; // n_stat_units \times n_locs data matrix - - // constructors - FPCA() = default; - fdapde_enable_constructor_if(is_space_only, This) FPCA(const pde_ptr& pde, Sampling s, Calibration c) : - Base(pde, s), calibration_(c) {}; - fdapde_enable_constructor_if(is_space_time_separable, This) - FPCA(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s, Calibration c) : - Base(space_penalty, time_penalty, s), calibration_(c) {}; - - void init_model() { return; }; - void solve() { - // preallocate space - loadings_.resize(n_basis(), n_pc_); - fitted_loadings_.resize(n_locs(), n_pc_); - scores_.resize(n_stat_units(), n_pc_); - DMatrix X_ = X(); // copy original data to avoid side effects - - // first guess of PCs set to a multivariate PCA (SVD) - Eigen::JacobiSVD> svd(X_, Eigen::ComputeThinU | Eigen::ComputeThinV); - PowerIteration solver(this, tolerance_, max_iter_); // power iteration solver - solver.set_seed(seed_); - solver.init(); - - // sequential extraction of principal components - for (std::size_t i = 0; i < n_pc_; i++) { - DVector f0 = svd.matrixV().col(i); - switch (calibration_) { - case Calibration::off: { - // find vectors s,f minimizing \norm_F{Y - s^T*f}^2 + (s^T*s)*P(f) fixed \lambda - solver.compute(X_, lambda(), f0); - } break; - case Calibration::gcv: { - // select \lambda minimizing the GCV index - ScalarField gcv([&solver, &X_, &f0](const DVector& lambda) -> double { - solver.compute(X_, lambda, f0); - return solver.gcv(); // return GCV index at convergence - }); - solver.compute(X_, core::Grid {}.optimize(gcv, lambda_grid_), f0); - } break; - case Calibration::kcv: { - auto cv_score = [&solver, &X_, &f0]( - const DVector& lambda, const core::BinaryVector& train_set, - const core::BinaryVector& test_set) -> double { - solver.compute(train_set.blk_repeat(1, X_.cols()).select(X_), lambda, f0); // fit on train set - // reconstruction error on test set: \norm{X_test * (I - fn*fn^\top/J)}_F/n_test, with - // J = \norm{f_n}_2^2 + f^\top*P(\lambda)*f (PS: the division of f^\top*P(\lambda)*f by - // \norm{f}_{L^2} is necessary to obtain J as expected) - return (test_set.blk_repeat(1, X_.cols()).select(X_) * - (DMatrix::Identity(X_.cols(), X_.cols()) - - solver.fn() * solver.fn().transpose() / - (solver.fn().squaredNorm() + solver.ftPf(lambda) / solver.f_squaredNorm()))) - .squaredNorm() / - test_set.count() * X_.cols(); - }; - solver.compute(X_, calibration::KCV {n_folds_, seed_}.fit(*this, lambda_grid_, cv_score), f0); - } break; - } - // store results - loadings_.col(i) = solver.f(); - fitted_loadings_.col(i) = solver.fn(); // \frac{\Psi * f}{\norm{f}_{L^2}} - scores_.col(i) = solver.s() * solver.f_norm(); - X_ -= scores_.col(i) * fitted_loadings_.col(i).transpose(); // X <- X - s*f_n^\top (deflation step) - } - return; - } - - // getters - const DMatrix& fitted_loadings() const { return fitted_loadings_; } - const DMatrix& loadings() const { return loadings_; } - const DMatrix& scores() const { return scores_; } - // setters - void set_npc(std::size_t n_pc) { n_pc_ = n_pc; } - void set_tolerance(double tolerance) { tolerance_ = tolerance; } - void set_max_iter(std::size_t max_iter) { max_iter_ = max_iter; } - void set_seed(std::size_t seed) { seed_ = seed; } - void set_lambda(const std::vector>& lambda_grid) { - fdapde_assert(calibration_ != Calibration::off); - lambda_grid_ = lambda_grid; - } - void set_nfolds(std::size_t n_folds) { - fdapde_assert(calibration_ == Calibration::kcv); - n_folds_ = n_folds; - } +class FPCA : public FunctionalBase, RegularizationType_> { private: - std::size_t n_pc_ = 3; // number of principal components - Calibration calibration_; // PC function's smoothing parameter selection strategy - std::size_t n_folds_ = 10; // for a kcv calibration strategy, the number of folds - std::vector> lambda_grid_; - // power iteration parameters - double tolerance_ = 1e-6; // relative tolerance between Jnew and Jold, used as stopping criterion - std::size_t max_iter_ = 20; // maximum number of allowed iterations - int seed_ = fdapde::random_seed; - - // problem solution - DMatrix loadings_; // PC functions' expansion coefficients - DMatrix fitted_loadings_; // evaluation of the PC functions at data locations - DMatrix scores_; -}; + std::size_t n_pc_ = 3; // number of principal components -// implementation of FPCA for monolithic approach, see e.g. -template -class FPCA : - public FunctionalBase, RegularizationType_> { + struct SolverType__ { // type erased solver strategy + using This_ = FPCA; + template using fn_ptrs = mem_fn_ptrs<&T::template compute, &T::loadings, &T::scores>; + void compute(const DMatrix& X, std::size_t rank, This_& model) { + invoke(*this, X, rank, model); + } + decltype(auto) loadings() const { return invoke&, 1>(*this); } + decltype(auto) scores() const { return invoke&, 2>(*this); } + }; + using SolverType = fdapde::erase; + SolverType solver_; public: using RegularizationType = std::decay_t; - using This = FPCA; + using This = FPCA; using Base = FunctionalBase; IMPORT_MODEL_SYMBOLS; - using Base::lambda; - using Base::X; + using Base::X; // n_stat_units \times n_locs data matrix // constructors FPCA() = default; - fdapde_enable_constructor_if(is_space_only, This) FPCA(const pde_ptr& pde, Sampling s) : Base(pde, s) {}; + fdapde_enable_constructor_if(is_space_only, This) FPCA(const pde_ptr& pde, Sampling s, SolverType solver) : + Base(pde, s), solver_(solver) {}; fdapde_enable_constructor_if(is_space_time_separable, This) - FPCA(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s) : - Base(space_penalty, time_penalty, s) {}; - - void init_model() { return; }; - void solve() { // monolithic approach via Regularized SVD - RSVD solver(this); - solver.compute(X(), lambda(), n_pc_); - // store results - loadings_ = solver.W(); - fitted_loadings_ = Psi(not_nan()) * solver.W(); // evaluation of L^2 normalized PC functions at data locations - scores_ = solver.H().array().rowwise() * solver.w_norm().transpose().array(); - return; - } + FPCA(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s, SolverType solver) : + Base(space_penalty, time_penalty, s), solver_(solver) {}; + void init_model() { fdapde_assert(bool(solver_) == true); }; // check solver is assigned + void solve() { solver_.compute(X(), n_pc_, *this); } // getters - const DMatrix& fitted_loadings() const { return fitted_loadings_; } - const DMatrix& loadings() const { return loadings_; } - const DMatrix& scores() const { return scores_; } + const DMatrix& loadings() const { return solver_.loadings(); } + const DMatrix& scores() const { return solver_.scores(); } + // setters void set_npc(std::size_t n_pc) { n_pc_ = n_pc; } - private: - std::size_t n_pc_ = 3; // number of principal components - // problem solution - DMatrix loadings_; // PC functions' expansion coefficients - DMatrix fitted_loadings_; // evaluation of the PC functions at data locations - DMatrix scores_; + void set_solver(SolverType solver) { solver_ = solver; } }; - } // namespace models } // namespace fdapde diff --git a/fdaPDE/models/functional/functional_base.h b/fdaPDE/models/functional/functional_base.h index 6aeda3cc..5271283d 100644 --- a/fdaPDE/models/functional/functional_base.h +++ b/fdaPDE/models/functional/functional_base.h @@ -49,7 +49,7 @@ class FunctionalBase : public select_regularization_base(s) {}; - // getters + // getters const DMatrix& X() const { return df_.template get(OBSERVATIONS_BLK); } // observation matrix X std::size_t n_stat_units() const { return X().rows(); } std::size_t n_obs() const { return X().size(); }; diff --git a/fdaPDE/models/functional/power_iteration.h b/fdaPDE/models/functional/power_iteration.h index a9953043..32fd066a 100644 --- a/fdaPDE/models/functional/power_iteration.h +++ b/fdaPDE/models/functional/power_iteration.h @@ -36,42 +36,40 @@ template class PowerIteration { using SolverType = typename std::conditional< is_space_only::value, SRPDE, STRPDE>::type; static constexpr int n_lambda = Model::n_lambda; - Model* m_; GCV gcv_; // GCV functor // algorithm's parameters double tolerance_ = 1e-6; // treshold on |Jnew - Jold| used as stopping criterion std::size_t max_iter_ = 20; // maximum number of iterations before forced stop std::size_t k_ = 0; // iteration index SolverType solver_; // internal solver - int seed_; + int seed_ = fdapde::random_seed; + int n_mc_samples_ = 100; DVector s_; // estimated score vector DVector fn_; // field evaluation at data location \frac{\Psi*f}{\norm{f}_{L^2}} DVector f_; // field basis expansion at convergence DVector g_; // PDE misfit at convergence double f_norm_; // L^2 norm of estimated field at converegence + public: // constructors PowerIteration() = default; - PowerIteration(Model* m, double tolerance, std::size_t max_iter, int seed) : - m_(m), tolerance_(tolerance), max_iter_(max_iter), - seed_((seed == fdapde::random_seed) ? std::random_device()() : seed) {}; - PowerIteration(Model* m, double tolerance, std::size_t max_iter) : - PowerIteration(m, tolerance, max_iter, fdapde::random_seed) {}; - - // initializes internal smoothing solver - void init() { - if constexpr (is_space_only::value) { // space-only solver - solver_ = SolverType(m_->pde(), m_->sampling()); - } else { // space-time separable solver - solver_ = SolverType(m_->pde(), m_->time_pde(), m_->sampling()); - solver_.set_temporal_locations(m_->time_locs()); + PowerIteration(const Model& m, double tolerance, std::size_t max_iter, int seed) : + tolerance_(tolerance), max_iter_(max_iter), + seed_((seed == fdapde::random_seed) ? std::random_device()() : seed) { + // initialize internal smoothing solver + if constexpr (is_space_only::value) { solver_ = SolverType(m.pde(), m.sampling()); } + else { + solver_ = SolverType(m.pde(), m.time_pde(), m.sampling()); + solver_.set_temporal_locations(m.time_locs()); } - solver_.set_spatial_locations(m_->locs()); - // initialize GCV functor - gcv_ = solver_.template gcv(100, seed_); - } + solver_.set_spatial_locations(m.locs()); + return; + }; + template PowerIteration(const ModelType& m, double tolerance, std::size_t max_iter) : + PowerIteration(m, tolerance, max_iter, fdapde::random_seed) {}; + void init() { gcv_ = solver_.template gcv(n_mc_samples_, seed_); } // executes the power itration algorithm on data X and smoothing parameter \lambda, starting from f0 void compute(const DMatrix& X, const SVector& lambda, const DVector& f0) { // initialization @@ -83,16 +81,16 @@ template class PowerIteration { double Jnew = 1; fn_ = f0; // set starting point while (!almost_equal(Jnew, Jold, tolerance_) && k_ < max_iter_) { - // compute score vector s as \frac{X*fn}{\norm(X*fn)} + // s = \frac{X*fn}{\norm(X*fn)} s_ = X * fn_; s_ = s_ / s_.norm(); - // compute loadings by solving a proper smoothing problem + // compute loadings (solve smoothing problem) solver_.data().template insert(OBSERVATIONS_BLK, X.transpose() * s_); // X^\top*s solver_.solve(); // prepare for next iteration k_++; fn_ = solver_.fitted(); // \Psi*f - // update value of discretized functional \norm{X - s*f_n^\top}_F + f^\top*P(\lambda)*f + // update value of discretized functional: \norm{X - s*f_n^\top}_F + f^\top*P(\lambda)*f Jold = Jnew; Jnew = (X - s_ * fn_.transpose()).squaredNorm() + solver_.ftPf(lambda); } diff --git a/fdaPDE/models/functional/regularized_svd.h b/fdaPDE/models/functional/regularized_svd.h new file mode 100644 index 00000000..83dde73e --- /dev/null +++ b/fdaPDE/models/functional/regularized_svd.h @@ -0,0 +1,171 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __REGULARIZED_SVD_H__ +#define __REGULARIZED_SVD_H__ + +#include +#include + +#include "../model_traits.h" + +namespace fdapde { +namespace models { + +// Let X be a data matrix made of noisy and discrete measurements of smooth functions sampled from a random field +// \mathcal{X}. RegularizedSVD implements the computation of a low-rank approximation of X using some regularizing term +template class RegularizedSVD; + +// Finds a low-rank approximation of X while penalizing for the eigenfunctions of \mathcal{X} by sequentially solving +// \argmin_{s,f} \norm_F{X - s^\top*f}^2 + (s^\top*s)*P_{\lambda}(f), up to a desired rank +template <> class RegularizedSVD { + public: + // constructors + RegularizedSVD(Calibration c) : calibration_(c) {} + RegularizedSVD() : RegularizedSVD(Calibration::off) {}; + + // sequentially solves \argmin_{s,f} \norm_F{X - s^\top*f}^2 + (s^\top*s)*P_{\lambda}(f), up to the specified rank, + // selecting the level of smoothing of the component according to the desired strategy + template void compute(const DMatrix& X, std::size_t rank, ModelType& model) { + // preallocate space + loadings_.resize(model.n_basis(), rank); + scores_.resize(X.rows(), rank); + loadings_norm_.resize(rank); + DMatrix X_ = X; // copy data + + // first guess of PCs set to a multivariate PCA (SVD) + Eigen::JacobiSVD> svd(X_, Eigen::ComputeThinU | Eigen::ComputeThinV); + // initialize power iteration solver + PowerIteration solver(model, tolerance_, max_iter_, seed_); + solver.init(); + // sequential extraction of principal components + for (std::size_t i = 0; i < rank; i++) { + DVector f0 = svd.matrixV().col(i); + switch (calibration_) { + case Calibration::off: { + // find vectors s,f minimizing \norm_F{Y - s^T*f}^2 + (s^T*s)*P(f) fixed \lambda + solver.compute(X_, model.lambda(), f0); + } break; + case Calibration::gcv: { + // select \lambda minimizing the GCV index + ScalarField gcv([&solver, &X_, &f0](const DVector& lambda) -> double { + solver.compute(X_, lambda, f0); + return solver.gcv(); // return GCV index at convergence + }); + solver.compute(X_, core::Grid {}.optimize(gcv, lambda_grid_), f0); + } break; + case Calibration::kcv: { + // select \lambda minimizing the reconstruction error in cross-validation + auto cv_score = [&solver, &X_, &f0]( + const DVector& lambda, const core::BinaryVector& train_set, + const core::BinaryVector& test_set) -> double { + solver.compute(train_set.blk_repeat(1, X_.cols()).select(X_), lambda, f0); // fit on train set + // reconstruction error on test set: \norm{X_test * (I - fn*fn^\top/J)}_F/n_test, with + // J = \norm{f_n}_2^2 + f^\top*P(\lambda)*f (PS: the division of f^\top*P(\lambda)*f by + // \norm{f}_{L^2} is necessary to obtain J as expected) + return (test_set.blk_repeat(1, X_.cols()).select(X_) * + (DMatrix::Identity(X_.cols(), X_.cols()) - + solver.fn() * solver.fn().transpose() / + (solver.fn().squaredNorm() + solver.ftPf(lambda) / solver.f_squaredNorm()))) + .squaredNorm() / + test_set.count() * X_.cols(); + }; + solver.compute(X_, calibration::KCV {n_folds_, seed_}.fit(model, lambda_grid_, cv_score), f0); + } break; + } + // store results + loadings_.col(i) = solver.f(); + scores_.col(i) = solver.s() * solver.f_norm(); + loadings_norm_[i] = solver.f_norm(); + X_ -= scores_.col(i) * solver.fn().transpose(); // X <- X - s*f_n^\top (deflation step) + } + return; + } + // getters + const DMatrix& scores() const { return scores_; } + const DMatrix& loadings() const { return loadings_; } + const DVector& loadings_norm() const { return loadings_norm_; } + // setters + void set_tolerance(double tolerance) { tolerance_ = tolerance; } + void set_max_iter(std::size_t max_iter) { max_iter_ = max_iter; } + void set_seed(std::size_t seed) { seed_ = seed; } + RegularizedSVD& set_lambda(const std::vector>& lambda_grid) { + fdapde_assert(calibration_ != Calibration::off); + lambda_grid_ = lambda_grid; + return *this; + } + RegularizedSVD& set_nfolds(std::size_t n_folds) { + fdapde_assert(calibration_ == Calibration::kcv); + n_folds_ = n_folds; + return *this; + } + private: + Calibration calibration_; // PC function's smoothing parameter selection strategy + std::size_t n_folds_ = 10; // for a kcv calibration strategy, the number of folds + std::vector> lambda_grid_; + // power iteration parameters + double tolerance_ = 1e-6; // relative tolerance between Jnew and Jold, used as stopping criterion + std::size_t max_iter_ = 20; // maximum number of allowed iterations + int seed_ = fdapde::random_seed; + + // problem solution + DMatrix scores_; + DMatrix loadings_; // PC functions' expansion coefficients + DVector loadings_norm_; // L^2 norm of estimated fields +}; + +// finds a rank r matrix U minimizing \norm{X - U*\Psi^\top}_F^2 + Tr[U*P_{\lambda}(f)*U^\top] +template <> class RegularizedSVD { + public: + // solves \norm{X - U*\Psi^\top}_F^2 + Tr[U*P_{\lambda}(f)*U^\top] retaining the first rank components + template void compute(const DMatrix& X, std::size_t rank, ModelType& model) { + // compute matrix C = \Psi^\top*\Psi + P(\lambda) + DMatrix C = model.Psi().transpose() * model.Psi() + model.P(); + // compute the inverse of the cholesky factor of C, D^{-1} + DMatrix D = C.llt().matrixL(); + DMatrix invD = D.inverse(); + + // compute SVD of X*\Psi*(D^{-1})^\top + Eigen::JacobiSVD> svd( + X * model.Psi() * invD.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV); + // store results + scores_ = svd.matrixU().leftCols(rank); + loadings_ = + (svd.singularValues().head(rank).asDiagonal() * svd.matrixV().leftCols(rank).transpose() * invD).transpose(); + loadings_norm_.resize(rank); + for (std::size_t i = 0; i < rank; ++i) { + loadings_norm_[i] = std::sqrt(loadings_.col(i).dot(model.R0() * loadings_.col(i))); // L^2 norm + loadings_.col(i) = loadings_.col(i) / loadings_norm_[i]; + } + scores_ = scores_.array().rowwise() * loadings_norm_.transpose().array(); + return; + } + // getters + const DMatrix& scores() const { return scores_; } + const DMatrix& loadings() const { return loadings_; } + const DVector& loadings_norm() const { return loadings_norm_; } + private: + // let E*\Sigma*F^\top the reduced (rank r) SVD of X*\Psi*(D^{1})^\top, with D^{-1} the inverse of the cholesky + // factor of \Psi^\top * \Psi + P(\lambda), then + DMatrix scores_; // matrix E in the reduced SVD of X*\Psi*(D^{-1})^\top + DMatrix loadings_; // \Sigma*F^\top*D^{-1} (PC functions expansion coefficients, L^2 normalized) + DVector loadings_norm_; // L^2 norm of estimated fields +}; + +} // namespace models +} // namespace fdapde + +#endif // __REGULARIZED_SVD_H__ diff --git a/fdaPDE/models/functional/rsvd.h b/fdaPDE/models/functional/rsvd.h deleted file mode 100644 index c60d9681..00000000 --- a/fdaPDE/models/functional/rsvd.h +++ /dev/null @@ -1,78 +0,0 @@ -// This file is part of fdaPDE, a C++ library for physics-informed -// spatial and functional data analysis. -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see . - -#ifndef __RSVD_H__ -#define __RSVD_H__ - -#include - -#include - -#include "../model_traits.h" - -namespace fdapde { -namespace models { - -// finds a rank r matrix U minimizing \norm{X - U*\Psi^\top}_F^2 + Tr[U*P_{\lambda}(f)*U^\top] -// with P_{\lambda}(f) the penalty term and \norm{}_F the Frobenius norm -template class RSVD { - private: - using Model = typename std::remove_reference::type; - static constexpr int n_lambda = Model::n_lambda; - const Model* m_; - - // let E*\Sigma*F^\top the reduced (rank r) SVD of X*\Psi*(D^{1})^\top, with D^{-1} the inverse of the cholesky - // factor of \Psi^\top * \Psi + P(\lambda), then - DMatrix H_; // matrix E in the reduced SVD of X*\Psi*(D^{-1})^\top (scores) - DMatrix W_; // \Sigma*F^\top*D^{-1} (field expansion coefficients, L^2 normalized) - DVector w_norm_; // L^2 norm of estimated fields - public: - // constructors - RSVD() = default; - RSVD(const Model* m) : m_(m) { } - - // executes the Regularized Singular Value Decomposition on data X for a given rank r and smoothing level \lambda - void compute(const DMatrix& X, const SVector& lambda, std::size_t r) { - // compute matrix C = \Psi^\top*\Psi + P(\lambda) - DMatrix C = m_->Psi().transpose() * m_->Psi() + m_->P(lambda); - // compute the inverse of the cholesky factor of C, D^{-1} - DMatrix D = C.llt().matrixL(); - DMatrix invD = D.inverse(); - - // compute SVD of X*\Psi*(D^{-1})^\top - Eigen::JacobiSVD> svd( - X * m_->Psi() * invD.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV); - // store results - H_ = svd.matrixU().leftCols(r); - W_ = (svd.singularValues().head(r).asDiagonal() * svd.matrixV().leftCols(r).transpose() * invD).transpose(); - w_norm_.resize(r); - for (std::size_t i = 0; i < r; ++i) { - w_norm_[i] = std::sqrt(W_.col(i).dot(m_->R0() * W_.col(i))); // L2 norm of estimated field - W_.col(i) = W_.col(i) / w_norm_[i]; - } - return; - } - - // getters - const DMatrix& H() const { return H_; } - const DMatrix& W() const { return W_; } - const DVector& w_norm() const { return w_norm_; } -}; - -} // namespace models -} // namespace fdapde - -#endif // __RSVD_H__ diff --git a/fdaPDE/models/model_type_erasure.h b/fdaPDE/models/model_type_erasure.h index 4d4d6f86..1f986b1a 100644 --- a/fdaPDE/models/model_type_erasure.h +++ b/fdaPDE/models/model_type_erasure.h @@ -25,89 +25,103 @@ namespace fdapde { namespace models { // minimal type erased wrappers for models -template struct IStatModel { }; +template struct StatisticalModel__ { }; -// base statistical models interface (do not permute fn_ptrs members!) -#define DEFINE_BASE_STAT_MODEL_INTERFACE \ - void init() { fdapde::invoke(*this); } \ - void solve() { fdapde::invoke(*this); } \ - void set_lambda_D(double lambda) { fdapde::invoke(*this, lambda); } \ +// base statistical models interface (do not permute BASE_MODEL_FN_PTRS members!) +#define BASE_MODEL_INTERFACE \ + void init() { invoke(*this); } \ + void solve() { invoke(*this); } \ void set_data(const BlockFrame& data, bool reindex = false) { \ - fdapde::invoke(*this, data, reindex); \ + invoke(*this, data, reindex); \ } \ - void set_spatial_locations(const DMatrix& locs) { fdapde::invoke(*this, locs); } \ - BlockFrame& data() { return fdapde::invoke&, 5>(*this); } \ - decltype(auto) R0() const { return fdapde::invoke&, 6>(*this); } \ - std::size_t n_locs() const { return fdapde::invoke(*this); } \ - std::size_t n_basis() const { return fdapde::invoke(*this); } + decltype(auto) n_locs() const { return invoke(*this); } \ + decltype(auto) n_basis() const { return invoke(*this); } \ + decltype(auto) data() { return invoke&, 5>(*this); } \ + decltype(auto) data() const { return invoke&, 6>(*this); } \ + decltype(auto) R0() const { return invoke&, 7>(*this); } \ + decltype(auto) R1() const { return invoke&, 8>(*this); } \ + decltype(auto) u() const { return invoke&, 9>(*this); } \ + decltype(auto) Psi() const { return invoke&, 10>(*this, not_nan()); } \ + decltype(auto) PsiTD() const { return invoke&, 11>(*this, not_nan()); } -// space-only statistical model interface -template <> struct IStatModel { - using RegularizationType = SpaceOnly; +#define BASE_MODEL_FN_PTRS \ + &M::init, &M::solve, &M::set_data, &M::n_locs, &M::n_basis, \ + static_cast& (ModelBase::*)()>(&M::data), \ + static_cast& (ModelBase::*)() const>(&M::data), &M::R0, &M::R1, &M::u, \ + static_cast& (SamplingBase::*)(not_nan) const>(&M::Psi), \ + static_cast& (SamplingBase::*)(not_nan) const>(&M::PsiTD) + +// basic model interface +template <> struct StatisticalModel__ { template using fn_ptrs = fdapde::mem_fn_ptrs< - &M::init, &M::solve, // initializtion and solution of the regression problem - &M::set_lambda_D, &M::set_data, &M::set_spatial_locations, - static_cast& (ModelBase::*)()>(&M::data), &M::R0, &M::n_locs, &M::n_basis, - static_cast::*)(const SVector<1>&)>(&M::set_lambda)>; + BASE_MODEL_FN_PTRS, static_cast::*)(const DVector&)>(&M::set_lambda), + static_cast (ModelBase::*)(int) const>(&M::lambda)>; // interface implementation - DEFINE_BASE_STAT_MODEL_INTERFACE; - void set_lambda(const SVector<1>& lambda) { fdapde::invoke(*this, lambda); } + BASE_MODEL_INTERFACE + void set_lambda(const DVector& lambda) { invoke(*this, lambda); } + decltype(auto) lambda() const { return invoke, 13>(*this, fdapde::Dynamic); } }; -// space-time separable statistical model interface -template <> struct IStatModel { - using RegularizationType = SpaceTimeSeparable; +// space-only statistical model interface +template <> struct StatisticalModel__ { + // compile time constants + using RegularizationType = SpaceOnly; + static constexpr int n_lambda = 1; + + // interface implementation template using fn_ptrs = fdapde::mem_fn_ptrs< - &M::init, &M::solve, // initializtion and solution of the regression problem - &M::set_lambda_D, &M::set_data, &M::set_spatial_locations, - static_cast& (ModelBase::*)()>(&M::data), &M::R0, &M::n_locs, &M::n_basis, - static_cast::*)(const SVector<2>&)>(&M::set_lambda), &M::set_lambda_T, - &M::set_temporal_locations>; - // interface implementation - DEFINE_BASE_STAT_MODEL_INTERFACE; - void set_lambda(const SVector<2>& lambda) { fdapde::invoke(*this, lambda); } - void set_lambda_T(double lambda_T) { fdapde::invoke(*this, lambda_T); } - void set_temporal_locations(const DMatrix& locs) { fdapde::invoke(*this, locs); } + BASE_MODEL_FN_PTRS, &M::set_lambda_D, static_cast (SpaceOnlyBase::*)() const>(&M::lambda), + static_cast::*)(const SVector<1>&)>(&M::set_lambda), &M::set_spatial_locations>; + BASE_MODEL_INTERFACE + void set_lambda_D(double lambda_D) { invoke(*this, lambda_D); } + SVector<1> lambda() const { return invoke, 13>(*this); } + void set_lambda(const SVector<1>& lambda) { invoke(*this, lambda); } + void set_spatial_locations(const DMatrix& locs) { invoke(*this, locs); } }; -// space-time parabolic statistical model interface -template <> struct IStatModel { - using RegularizationType = SpaceTimeParabolic; +// space-time separable interface +template <> struct StatisticalModel__ { + // compile time constants + using RegularizationType = SpaceTimeSeparable; + static constexpr int n_lambda = 2; + + // interface implementation template using fn_ptrs = fdapde::mem_fn_ptrs< - &M::init, &M::solve, // initializtion and solution of the regression problem - &M::set_lambda_D, &M::set_data, &M::set_spatial_locations, - static_cast& (ModelBase::*)()>(&M::data), &M::R0, &M::n_locs, &M::n_basis, - static_cast::*)(const SVector<2>&)>(&M::set_lambda), &M::set_lambda_T>; - // interface implementation - DEFINE_BASE_STAT_MODEL_INTERFACE; - void set_lambda(const SVector<2>& lambda) { fdapde::invoke(*this, lambda); } - void set_lambda_T(double lambda_T) { fdapde::invoke(*this, lambda_T); } + BASE_MODEL_FN_PTRS, &M::set_lambda_D, &M::set_lambda_T, + static_cast (SpaceTimeBase::*)() const>(&M::lambda), + static_cast::*)(const SVector<2>&)>(&M::set_lambda), + &M::set_spatial_locations, &M::set_temporal_locations>; + BASE_MODEL_INTERFACE + void set_lambda_D(double lambda_D) { invoke(*this, lambda_D); } + void set_lambda_T(double lambda_T) { invoke(*this, lambda_T); } + SVector<2> lambda() const { return invoke, 14>(*this); } + void set_lambda(const SVector<2>& lambda) { invoke(*this, lambda); } + void set_spatial_locations(const DMatrix& locs) { invoke(*this, locs); } + void set_temporal_locations(const DMatrix& locs) { invoke(*this, locs); } }; -// regularization-independent model interface -template <> struct IStatModel { +// space-time parabolic interface +template <> struct StatisticalModel__ { + // compile time constants + using RegularizationType = SpaceTimeParabolic; + static constexpr int n_lambda = 2; + + // interface implementation template using fn_ptrs = fdapde::mem_fn_ptrs< - &M::init, &M::solve, static_cast::*)(const DVector&)>(&M::set_lambda), - static_cast (ModelBase::*)(int) const>(&M::lambda), - static_cast& (SamplingBase::*)(not_nan) const>(&M::Psi), - static_cast& (SamplingBase::*)(not_nan) const>(&M::PsiTD), &M::set_data, - &M::n_locs, &M::n_basis>; - // interface implementation - void init() { fdapde::invoke(*this); } - void solve() { fdapde::invoke(*this); } - void set_lambda(const DVector& lambda) { fdapde::invoke(*this, lambda); } - decltype(auto) lambda() const { return fdapde::invoke, 3>(*this, fdapde::Dynamic); } - decltype(auto) Psi() const { return fdapde::invoke&, 4>(*this, not_nan()); } - decltype(auto) PsiTD() const { return fdapde::invoke&, 5>(*this, not_nan()); } - void set_data(const BlockFrame& data, bool reindex = false) { - fdapde::invoke(*this, data, reindex); - } - std::size_t n_locs() const { return fdapde::invoke(*this); } - std::size_t n_basis() const { return fdapde::invoke(*this); } + BASE_MODEL_FN_PTRS, &M::set_lambda_D, &M::set_lambda_T, + static_cast (SpaceTimeBase::*)() const>(&M::lambda), + static_cast::*)(const SVector<2>&)>(&M::set_lambda), + &M::set_spatial_locations>; + BASE_MODEL_INTERFACE + void set_lambda_D(double lambda_D) { invoke(*this, lambda_D); } + void set_lambda_T(double lambda_T) { invoke(*this, lambda_T); } + SVector<2> lambda() const { return invoke, 14>(*this); } + void set_lambda(const SVector<2>& lambda) { invoke(*this, lambda); } + void set_spatial_locations(const DMatrix& locs) { invoke(*this, locs); } }; } // namespace models diff --git a/fdaPDE/models/regression/exact_edf.h b/fdaPDE/models/regression/exact_edf.h index 37492d8d..3b9c4f0c 100644 --- a/fdaPDE/models/regression/exact_edf.h +++ b/fdaPDE/models/regression/exact_edf.h @@ -27,8 +27,7 @@ namespace models { // operator: Tr[S] = Tr[\Psi*T^{-1}*\Psi^T*Q] = Tr[Q*\Psi*T^{-1}*\Psi^T] class ExactEDF { private: - using ModelType = fdapde::erase, IRegression>; - ModelType model_; + RegressionView model_; // computes smoothing matrix S = Q*\Psi*T^{-1}*\Psi^T const DMatrix& S() { // factorize matrix T @@ -43,7 +42,7 @@ class ExactEDF { ExactEDF() = default; double compute() { return S().trace(); } // computes Tr[S] - void set_model(const ModelType& model) { model_ = model; } + void set_model(RegressionView model) { model_ = model; } }; } // namespace models diff --git a/fdaPDE/models/regression/gcv.h b/fdaPDE/models/regression/gcv.h index d88dd542..b3e5ab44 100644 --- a/fdaPDE/models/regression/gcv.h +++ b/fdaPDE/models/regression/gcv.h @@ -38,9 +38,7 @@ struct EDFStrategy__ { template using fn_ptrs = fdapde::mem_fn_ptrs<&M::compute, &M::set_model>; // forwardings decltype(auto) compute() { return fdapde::invoke(*this); } - decltype(auto) set_model(const fdapde::erase, IRegression>& model) { - fdapde::invoke(*this, model); - } + void set_model(const RegressionView& model) { fdapde::invoke(*this, model); } }; using EDFStrategy = fdapde::erase; diff --git a/fdaPDE/models/regression/regression_type_erasure.h b/fdaPDE/models/regression/regression_type_erasure.h index 8076ffd4..f49c7e9f 100644 --- a/fdaPDE/models/regression/regression_type_erasure.h +++ b/fdaPDE/models/regression/regression_type_erasure.h @@ -32,54 +32,43 @@ namespace fdapde { namespace models { // type erased wrapper for regression models -struct IRegression { +struct RegressionModel__ { template using fn_ptrs = fdapde::mem_fn_ptrs< &M::f, &M::beta, &M::g, &M::fitted, &M::W, &M::XtWX, &M::U, &M::V, &M::invXtWX, &M::invA, &M::q, &M::n_obs, &M::norm, &M::y, &M::T, &M::lmbQ, &M::has_covariates, &M::nan_mask, &M::mask_obs, &M::X>; // interface implementation - decltype(auto) f() const { return fdapde::invoke& , 0>(*this); } - decltype(auto) beta() const { return fdapde::invoke& , 1>(*this); } - decltype(auto) g() const { return fdapde::invoke& , 2>(*this); } - decltype(auto) fitted() const { return fdapde::invoke , 3>(*this); } - decltype(auto) W() const { return fdapde::invoke&, 4>(*this); } - decltype(auto) XtWX() const { return fdapde::invoke& , 5>(*this); } - decltype(auto) U() const { return fdapde::invoke& , 6>(*this); } - decltype(auto) V() const { return fdapde::invoke& , 7>(*this); } - decltype(auto) invXtWX() const { return fdapde::invoke>&, 8>(*this); } - decltype(auto) invA() const { return fdapde::invoke>& , 9>(*this); } - decltype(auto) q() const { return fdapde::invoke(*this); } - decltype(auto) n_obs() const { return fdapde::invoke(*this); } + decltype(auto) f() const { return invoke& , 0>(*this); } + decltype(auto) beta() const { return invoke& , 1>(*this); } + decltype(auto) g() const { return invoke& , 2>(*this); } + decltype(auto) fitted() const { return invoke , 3>(*this); } + decltype(auto) W() const { return invoke&, 4>(*this); } + decltype(auto) XtWX() const { return invoke& , 5>(*this); } + decltype(auto) U() const { return invoke& , 6>(*this); } + decltype(auto) V() const { return invoke& , 7>(*this); } + decltype(auto) invXtWX() const { return invoke>&, 8>(*this); } + decltype(auto) invA() const { return invoke>& , 9>(*this); } + decltype(auto) q() const { return invoke(*this); } + decltype(auto) n_obs() const { return invoke(*this); } decltype(auto) norm(const DMatrix& op1, const DMatrix& op2) const { - return fdapde::invoke (*this, op1, op2); + return invoke (*this, op1, op2); } - decltype(auto) y() const { return fdapde::invoke&, 13>(*this); } - decltype(auto) T() { return fdapde::invoke&, 14>(*this); } - decltype(auto) lmbQ(const DMatrix& x) const { return fdapde::invoke, 15>(*this, x); } - decltype(auto) has_covariates() const { return fdapde::invoke(*this); } - decltype(auto) nan_mask() const { return fdapde::invoke&, 17>(*this); } - decltype(auto) mask_obs(const BinaryVector& mask) { return fdapde::invoke(*this, mask); } - decltype(auto) X() const { return fdapde::invoke&, 19>(*this); } -}; + decltype(auto) y() const { return invoke&, 13>(*this); } + decltype(auto) T() { return invoke&, 14>(*this); } + decltype(auto) lmbQ(const DMatrix& x) const { return invoke, 15>(*this, x); } + decltype(auto) has_covariates() const { return invoke(*this); } + decltype(auto) nan_mask() const { return invoke&, 17>(*this); } + decltype(auto) mask_obs(const BinaryVector& mask) { return invoke(*this, mask); } + decltype(auto) X() const { return invoke&, 19>(*this); } +}; template -using RegressionModel = fdapde::erase, IRegression>; +using RegressionModel = fdapde::erase, RegressionModel__>; template -using RegressionView = fdapde::erase, IRegression>; +using RegressionView = + fdapde::erase, RegressionModel__>; } // namespace models - -template -typename std::enable_if< is_space_only::value, models::RegressionModel>::type -make_model(PDE&& args, Sampling s) { - return models::RegressionModel(Model(std::forward(args), s)); -} -template -typename std::enable_if::value, models::RegressionModel>::type -make_model(PDE&& args, Sampling s, const DVector& time) { - return models::RegressionModel(Model(std::forward(args), s, time)); -} - } // namespace fdapde #endif // __REGRESSION_TYPE_ERASURE_H__ diff --git a/fdaPDE/models/regression/stochastic_edf.h b/fdaPDE/models/regression/stochastic_edf.h index 4c4e24af..aba148a6 100644 --- a/fdaPDE/models/regression/stochastic_edf.h +++ b/fdaPDE/models/regression/stochastic_edf.h @@ -30,8 +30,7 @@ namespace models { // computes an approximation of the trace of S = \Psi*T^{-1}*\Psi^T*Q using a monte carlo approximation. class StochasticEDF { private: - using ModelType = fdapde::erase, IRegression>; - ModelType model_; + RegressionView model_; std::size_t r_ = 100; // number of monte carlo realizations DMatrix Us_; // sample from Rademacher distribution DMatrix Bs_; // \Psi^T*Q*Us_ @@ -86,7 +85,7 @@ class StochasticEDF { return MCmean / r_; } // setter - void set_model(ModelType& model) { model_ = model; } + void set_model(RegressionView model) { model_ = model; } void set_seed(int seed) { seed_ = seed; } void set_n_mc_samples(int r) { r_ = r; } }; diff --git a/test/src/fpca_test.cpp b/test/src/fpca_test.cpp index 26997a9b..aad71ee1 100644 --- a/test/src/fpca_test.cpp +++ b/test/src/fpca_test.cpp @@ -28,6 +28,7 @@ using fdapde::core::PDE; #include "../../fdaPDE/models/functional/fpca.h" #include "../../fdaPDE/models/sampling_design.h" using fdapde::models::FPCA; +using fdapde::models::RegularizedSVD; using fdapde::models::Sampling; #include "../../fdaPDE/calibration/symbols.h" using fdapde::calibration::Calibration; @@ -56,10 +57,10 @@ TEST(fpca_test, laplacian_samplingatnodes_sequential) { // define regularizing PDE auto L = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + PDE, FEM, fem_order<1>> pde(domain.mesh, L, u); // define model double lambda_D = 1e-2; - FPCA model(problem, Sampling::mesh_nodes, Calibration::off); + FPCA model(pde, Sampling::mesh_nodes, RegularizedSVD{Calibration::off}); model.set_lambda_D(lambda_D); // set model's data BlockFrame df; @@ -69,8 +70,8 @@ TEST(fpca_test, laplacian_samplingatnodes_sequential) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.fitted_loadings(), "../data/models/fpca/2D_test1/loadings_seq.mtx")); - EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test1/scores_seq.mtx" )); + EXPECT_TRUE(almost_equal(model.Psi() * model.loadings(), "../data/models/fpca/2D_test1/loadings_seq.mtx")); + EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test1/scores_seq.mtx" )); } // test 2 @@ -92,7 +93,7 @@ TEST(fpca_test, laplacian_samplingatnodes_monolithic) { PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); // define model double lambda_D = 1e-2; - FPCA model(problem, Sampling::mesh_nodes); + FPCA model(problem, Sampling::mesh_nodes, RegularizedSVD()); model.set_lambda_D(lambda_D); // set model's data BlockFrame df; @@ -102,8 +103,8 @@ TEST(fpca_test, laplacian_samplingatnodes_monolithic) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.fitted_loadings(), "../data/models/fpca/2D_test1/loadings_mon.mtx")); - EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test1/scores_mon.mtx" )); + EXPECT_TRUE(almost_equal(model.Psi() * model.loadings(), "../data/models/fpca/2D_test1/loadings_mon.mtx")); + EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test1/scores_mon.mtx" )); } // test 3 @@ -123,15 +124,16 @@ TEST(fpca_test, laplacian_samplingatlocations_sequential_gcv) { // define regularizing PDE auto L = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); - // define model - FPCA model(problem, Sampling::pointwise, Calibration::gcv); - model.set_spatial_locations(locs); + PDE, FEM, fem_order<1>> pde(domain.mesh, L, u); // grid of smoothing parameters std::vector> lambda_grid; for (double x = -4; x <= -2; x += 0.1) { lambda_grid.push_back(SVector<1>(std::pow(10, x))); } - model.set_lambda(lambda_grid); - model.set_seed(78965); // for reproducibility purposes in testing + // define model + RegularizedSVD rsvd(Calibration::gcv); + rsvd.set_lambda(lambda_grid); + rsvd.set_seed(78965); // for reproducibility purposes in testing + FPCA model(pde, Sampling::pointwise, rsvd); + model.set_spatial_locations(locs); // set model's data BlockFrame df; df.insert(OBSERVATIONS_BLK, y); @@ -140,8 +142,8 @@ TEST(fpca_test, laplacian_samplingatlocations_sequential_gcv) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.fitted_loadings(), "../data/models/fpca/2D_test2/loadings.mtx")); - EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test2/scores.mtx" )); + EXPECT_TRUE(almost_equal(model.Psi() * model.loadings(), "../data/models/fpca/2D_test2/loadings.mtx")); + EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test2/scores.mtx" )); } // test 4 @@ -162,15 +164,15 @@ TEST(fpca_test, laplacian_samplingatlocations_sequential_kcv) { auto L = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + // grid of smoothing parameters + std::vector> lambda_grid; + for (double x = -4; x <= -2; x += 0.1) lambda_grid.push_back(SVector<1>(std::pow(10, x))); // define model - FPCA model(problem, Sampling::pointwise, Calibration::kcv); + RegularizedSVD rsvd(Calibration::kcv); + rsvd.set_lambda(lambda_grid); + rsvd.set_seed(12654); // for reproducibility purposes in testing + FPCA model(problem, Sampling::pointwise, rsvd); model.set_spatial_locations(locs); - // grid of smoothing parameters - std::vector> lambdas; - for (double x = -4; x <= -2; x += 0.1) lambdas.push_back(SVector<1>(std::pow(10, x))); - model.set_lambda(lambdas); - model.set_seed(12654); // for reproducibility purposes in testing - model.set_nfolds(10); // perform a 10 folds cross-validation // set model's data BlockFrame df; df.insert(OBSERVATIONS_BLK, y); @@ -179,8 +181,8 @@ TEST(fpca_test, laplacian_samplingatlocations_sequential_kcv) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.fitted_loadings(), "../data/models/fpca/2D_test3/loadings.mtx")); - EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test3/scores.mtx" )); + EXPECT_TRUE(almost_equal(model.Psi() * model.loadings(), "../data/models/fpca/2D_test3/loadings.mtx")); + EXPECT_TRUE(almost_equal(model.scores(), "../data/models/fpca/2D_test3/scores.mtx" )); } /*