From 3fa50063861e02632adc9430019fff238de75fa7 Mon Sep 17 00:00:00 2001 From: AlePalu Date: Sun, 10 Dec 2023 19:08:21 +0100 Subject: [PATCH] quantile spatial regression models are now named QSRPDE, GCV for QSRPDE restored. FPIRLS is default constructable and is stored as a model data member --- .gitignore | 64 ++++---- fdaPDE/models/model_base.h | 6 +- fdaPDE/models/regression/fpirls.h | 45 +++--- fdaPDE/models/regression/gcv.h | 3 +- fdaPDE/models/regression/gsrpde.h | 59 ++++---- .../models/regression/{sqrpde.h => qsrpde.h} | 86 +++++------ .../gcv/{sqrpde => qsrpde}/2D_test1/edfs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test1/gcvs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test1/y.csv | 0 test/data/gcv/qsrpde/2D_test2/edfs.mtx | 15 ++ test/data/gcv/qsrpde/2D_test2/gcvs.mtx | 15 ++ .../gcv/{sqrpde => qsrpde}/2D_test2/y.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test3/X.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test3/edfs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test3/gcvs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test3/locs.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test3/y.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test4/X.csv | 0 test/data/gcv/qsrpde/2D_test4/edfs.mtx | 11 ++ test/data/gcv/qsrpde/2D_test4/gcvs.mtx | 11 ++ .../gcv/{sqrpde => qsrpde}/2D_test4/locs.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test4/y.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test5/edfs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test5/gcvs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test5/y.csv | 0 test/data/gcv/qsrpde/2D_test6/edfs.mtx | 11 ++ test/data/gcv/qsrpde/2D_test6/gcvs.mtx | 11 ++ .../gcv/{sqrpde => qsrpde}/2D_test6/y.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test7/X.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test7/edfs.mtx | 0 .../gcv/{sqrpde => qsrpde}/2D_test7/gcvs.mtx | 0 .../2D_test7/incidence_matrix.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test7/y.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test8/X.csv | 0 test/data/gcv/qsrpde/2D_test8/edfs.mtx | 15 ++ test/data/gcv/qsrpde/2D_test8/gcvs.mtx | 15 ++ .../2D_test8/incidence_matrix.csv | 0 .../gcv/{sqrpde => qsrpde}/2D_test8/y.csv | 0 test/data/gcv/sqrpde/2D_test2/edfs.mtx | 15 -- test/data/gcv/sqrpde/2D_test2/gcvs.mtx | 15 -- test/data/gcv/sqrpde/2D_test4/edfs.mtx | 11 -- test/data/gcv/sqrpde/2D_test4/gcvs.mtx | 11 -- test/data/gcv/sqrpde/2D_test6/edfs.mtx | 11 -- test/data/gcv/sqrpde/2D_test6/gcvs.mtx | 11 -- test/data/gcv/sqrpde/2D_test8/edfs.mtx | 15 -- test/data/gcv/sqrpde/2D_test8/gcvs.mtx | 15 -- .../{sqrpde => qsrpde}/2D_test1/README.md | 0 .../{sqrpde => qsrpde}/2D_test1/sol.mtx | 0 .../models/{sqrpde => qsrpde}/2D_test1/y.csv | 0 .../{sqrpde => qsrpde}/2D_test2/README.md | 0 .../models/{sqrpde => qsrpde}/2D_test2/X.csv | 0 .../{sqrpde => qsrpde}/2D_test2/beta.mtx | 0 .../{sqrpde => qsrpde}/2D_test2/locs.csv | 0 .../{sqrpde => qsrpde}/2D_test2/sol.mtx | 0 .../models/{sqrpde => qsrpde}/2D_test2/y.csv | 0 .../{sqrpde => qsrpde}/2D_test3/README.md | 0 .../{sqrpde => qsrpde}/2D_test3/sol.mtx | 0 .../models/{sqrpde => qsrpde}/2D_test3/y.csv | 0 .../{sqrpde => qsrpde}/2D_test4/README.md | 0 .../models/{sqrpde => qsrpde}/2D_test4/X.csv | 0 .../{sqrpde => qsrpde}/2D_test4/beta.mtx | 0 .../2D_test4/incidence_matrix.csv | 0 .../{sqrpde => qsrpde}/2D_test4/sol.mtx | 0 .../models/{sqrpde => qsrpde}/2D_test4/y.csv | 0 test/main.cpp | 4 +- ...cv_sqrpde_test.cpp => gcv_qsrpde_test.cpp} | 137 +++++++++--------- test/src/{sqrpde_test.cpp => qsrpde_test.cpp} | 48 +++--- 67 files changed, 334 insertions(+), 326 deletions(-) rename fdaPDE/models/regression/{sqrpde.h => qsrpde.h} (76%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test1/edfs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test1/gcvs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test1/y.csv (100%) create mode 100644 test/data/gcv/qsrpde/2D_test2/edfs.mtx create mode 100644 test/data/gcv/qsrpde/2D_test2/gcvs.mtx rename test/data/gcv/{sqrpde => qsrpde}/2D_test2/y.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test3/X.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test3/edfs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test3/gcvs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test3/locs.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test3/y.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test4/X.csv (100%) create mode 100644 test/data/gcv/qsrpde/2D_test4/edfs.mtx create mode 100644 test/data/gcv/qsrpde/2D_test4/gcvs.mtx rename test/data/gcv/{sqrpde => qsrpde}/2D_test4/locs.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test4/y.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test5/edfs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test5/gcvs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test5/y.csv (100%) create mode 100644 test/data/gcv/qsrpde/2D_test6/edfs.mtx create mode 100644 test/data/gcv/qsrpde/2D_test6/gcvs.mtx rename test/data/gcv/{sqrpde => qsrpde}/2D_test6/y.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test7/X.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test7/edfs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test7/gcvs.mtx (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test7/incidence_matrix.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test7/y.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test8/X.csv (100%) create mode 100644 test/data/gcv/qsrpde/2D_test8/edfs.mtx create mode 100644 test/data/gcv/qsrpde/2D_test8/gcvs.mtx rename test/data/gcv/{sqrpde => qsrpde}/2D_test8/incidence_matrix.csv (100%) rename test/data/gcv/{sqrpde => qsrpde}/2D_test8/y.csv (100%) delete mode 100644 test/data/gcv/sqrpde/2D_test2/edfs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test2/gcvs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test4/edfs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test4/gcvs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test6/edfs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test6/gcvs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test8/edfs.mtx delete mode 100644 test/data/gcv/sqrpde/2D_test8/gcvs.mtx rename test/data/models/{sqrpde => qsrpde}/2D_test1/README.md (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test1/sol.mtx (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test1/y.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test2/README.md (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test2/X.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test2/beta.mtx (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test2/locs.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test2/sol.mtx (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test2/y.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test3/README.md (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test3/sol.mtx (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test3/y.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test4/README.md (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test4/X.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test4/beta.mtx (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test4/incidence_matrix.csv (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test4/sol.mtx (100%) rename test/data/models/{sqrpde => qsrpde}/2D_test4/y.csv (100%) rename test/src/{gcv_sqrpde_test.cpp => gcv_qsrpde_test.cpp} (73%) rename test/src/{sqrpde_test.cpp => qsrpde_test.cpp} (79%) diff --git a/.gitignore b/.gitignore index f93cdcac..d9ef1f6c 100644 --- a/.gitignore +++ b/.gitignore @@ -52,8 +52,8 @@ !test/data/models/gsrpde/* !test/data/models/gcv/ !test/data/models/gcv/* -!test/data/gcv/sqrpde/ -!test/data/gcv/sqrpde/* +!test/data/gcv/qsrpde/ +!test/data/gcv/qsrpde/* !test/data/models/fpca/ !test/data/models/fpca/* !.github/ @@ -114,20 +114,20 @@ !test/data/models/gsrpde/2D_test5/* !test/data/models/gsrpde/2D_test6/ !test/data/models/gsrpde/2D_test6/* -!test/data/models/sqrpde/ -!test/data/models/sqrpde/* -!test/data/models/sqrpde/2D_test1/ -!test/data/models/sqrpde/2D_test1/* -!test/data/models/sqrpde/2D_test2/ -!test/data/models/sqrpde/2D_test2/* -!test/data/models/sqrpde/2D_test3/ -!test/data/models/sqrpde/2D_test3/* -!test/data/models/sqrpde/2D_test4/ -!test/data/models/sqrpde/2D_test4/* -!test/data/models/sqrpde/2D_test5/ -!test/data/models/sqrpde/2D_test5/* -!test/data/models/sqrpde/2D_test6/ -!test/data/models/sqrpde/2D_test6/* +!test/data/models/qsrpde/ +!test/data/models/qsrpde/* +!test/data/models/qsrpde/2D_test1/ +!test/data/models/qsrpde/2D_test1/* +!test/data/models/qsrpde/2D_test2/ +!test/data/models/qsrpde/2D_test2/* +!test/data/models/qsrpde/2D_test3/ +!test/data/models/qsrpde/2D_test3/* +!test/data/models/qsrpde/2D_test4/ +!test/data/models/qsrpde/2D_test4/* +!test/data/models/qsrpde/2D_test5/ +!test/data/models/qsrpde/2D_test5/* +!test/data/models/qsrpde/2D_test6/ +!test/data/models/qsrpde/2D_test6/* !test/data/models/gcv/2D_test1/ !test/data/models/gcv/2D_test1/* !test/data/models/gcv/2D_test2/ @@ -144,22 +144,22 @@ !test/data/models/gcv/2D_test7/* !test/data/models/gcv/2D_test8/ !test/data/models/gcv/2D_test8/* -!test/data/gcv/sqrpde/2D_test1/ -!test/data/gcv/sqrpde/2D_test1/* -!test/data/gcv/sqrpde/2D_test2/ -!test/data/gcv/sqrpde/2D_test2/* -!test/data/gcv/sqrpde/2D_test3/ -!test/data/gcv/sqrpde/2D_test3/* -!test/data/gcv/sqrpde/2D_test4/ -!test/data/gcv/sqrpde/2D_test4/* -!test/data/gcv/sqrpde/2D_test5/ -!test/data/gcv/sqrpde/2D_test5/* -!test/data/gcv/sqrpde/2D_test6/ -!test/data/gcv/sqrpde/2D_test6/* -!test/data/gcv/sqrpde/2D_test7/ -!test/data/gcv/sqrpde/2D_test7/* -!test/data/gcv/sqrpde/2D_test8/ -!test/data/gcv/sqrpde/2D_test8/* +!test/data/gcv/qsrpde/2D_test1/ +!test/data/gcv/qsrpde/2D_test1/* +!test/data/gcv/qsrpde/2D_test2/ +!test/data/gcv/qsrpde/2D_test2/* +!test/data/gcv/qsrpde/2D_test3/ +!test/data/gcv/qsrpde/2D_test3/* +!test/data/gcv/qsrpde/2D_test4/ +!test/data/gcv/qsrpde/2D_test4/* +!test/data/gcv/qsrpde/2D_test5/ +!test/data/gcv/qsrpde/2D_test5/* +!test/data/gcv/qsrpde/2D_test6/ +!test/data/gcv/qsrpde/2D_test6/* +!test/data/gcv/qsrpde/2D_test7/ +!test/data/gcv/qsrpde/2D_test7/* +!test/data/gcv/qsrpde/2D_test8/ +!test/data/gcv/qsrpde/2D_test8/* !test/data/models/fpca/2D_test1/ !test/data/models/fpca/2D_test1/* !test/data/models/fpca/2D_test2/ diff --git a/fdaPDE/models/model_base.h b/fdaPDE/models/model_base.h index d402b93a..edccd01e 100644 --- a/fdaPDE/models/model_base.h +++ b/fdaPDE/models/model_base.h @@ -41,9 +41,9 @@ template class ModelBase { if (model().runtime().query(runtime_status::require_pde_init)) { pde_.init(); } // init penalty if (model().runtime().query(runtime_status::require_penalty_init)) { model().init_regularization(); } if (model().runtime().query(runtime_status::require_functional_basis_evaluation)) { - model().init_sampling(true); // init \Psi matrix, always force recomputation - model().init_nan(); // analyze and set missingness pattern - } + model().init_sampling(true); // init \Psi matrix, always force recomputation + model().init_nan(); // analyze and set missingness pattern + } model().init_data(); // specific data-dependent initialization requested by the model model().init_model(); // model initialization diff --git a/fdaPDE/models/regression/fpirls.h b/fdaPDE/models/regression/fpirls.h index 642b2837..a5880de0 100644 --- a/fdaPDE/models/regression/fpirls.h +++ b/fdaPDE/models/regression/fpirls.h @@ -37,7 +37,7 @@ template class FPIRLS { private: using Model = typename std::decay::type; using SolverWrapper = RegressionModel::RegularizationType>; - Model& m_; + Model* m_; // algorithm's parameters double tolerance_; // treshold on objective functional J to convergence std::size_t max_iter_; // maximum number of iterations before forced stop @@ -45,47 +45,50 @@ template class FPIRLS { SolverWrapper solver_; // internal solver public: // constructor - FPIRLS(Model& m, double tolerance, std::size_t max_iter) : - m_(m), tolerance_(tolerance), max_iter_(max_iter) { - if (!solver_) { // default solver initialization + FPIRLS() = default; + FPIRLS(Model* m, double tolerance, std::size_t max_iter) : m_(m), tolerance_(tolerance), max_iter_(max_iter) {}; + + // initialize internal smoothing solver + void init() { + if (!solver_) { // default solver initialization using SolverType = typename std::conditional< is_space_only::value, SRPDE, STRPDE >::type; // define internal problem solver if constexpr (!is_space_time::value) // space-only - solver_ = SolverType(m_.pde(), m_.sampling()); + solver_ = SolverType(m_->pde(), m_->sampling()); else { // space-time - solver_ = SolverType(m_.pde(), m_.sampling(), m_.time_domain()); - if constexpr (is_space_time_parabolic::value) { solver_.set_initial_condition(m_.s(), false); } + solver_ = SolverType(m_->pde(), m_->sampling(), m_->time_domain()); + if constexpr (is_space_time_parabolic::value) { solver_.set_initial_condition(m_->s(), false); } if constexpr (is_space_time_separable::value) { - solver_.set_temporal_locations(m_.time_locs()); + solver_.set_temporal_locations(m_->time_locs()); } } // solver initialization - solver_.set_lambda(m_.lambda()); - solver_.set_spatial_locations(m_.locs()); - solver_.set_data(m_.data()); // possible covariates are passed from here + solver_.set_spatial_locations(m_->locs()); + solver_.set_data(m_->data()); // possible covariates are passed from here } - }; - + solver_.set_lambda(m_->lambda()); + } // executes the FPIRLS algorithm void compute() { - m_.fpirls_init(); + // update solver + m_->fpirls_init(); // objective functional value at consecutive iterations - double J_old = tolerance_ + 1; - double J_new = 0; + double J_old = tolerance_ + 1, J_new = 0; + k_ = 0; while (k_ < max_iter_ && std::abs(J_new - J_old) > tolerance_) { - m_.fpirls_compute_step(); // model specific computation of py_ and pW_ + m_->fpirls_compute_step(); // model specific computation of py_ and pW_ // solve weighted least square problem // \argmin_{\beta, f} [ \norm(W^{1/2}(y - X\beta - f_n))^2 + \lambda \int_D (Lf - u)^2 ] - solver_.data().template insert(OBSERVATIONS_BLK, m_.py()); - solver_.data().template insert(WEIGHTS_BLK, m_.pW()); + solver_.data().template insert(OBSERVATIONS_BLK, m_->py()); + solver_.data().template insert(WEIGHTS_BLK, m_->pW()); // update solver and solve solver_.init(); solver_.solve(); - m_.fpirls_update_step(solver_.fitted(), solver_.beta()); // model specific update step + m_->fpirls_update_step(solver_.fitted(), solver_.beta()); // model specific update step // update value of the objective functional J = data_loss + \int_D (Lf-u)^2 - double J = m_.data_loss() + m_.lambda_D()*solver_.g().dot(m_.R0() * solver_.g()); + double J = m_->data_loss() + m_->lambda_D()*solver_.g().dot(m_->R0() * solver_.g()); // prepare for next iteration k_++; J_old = J_new; J_new = J; } diff --git a/fdaPDE/models/regression/gcv.h b/fdaPDE/models/regression/gcv.h index 8f9253fd..74d4b389 100644 --- a/fdaPDE/models/regression/gcv.h +++ b/fdaPDE/models/regression/gcv.h @@ -51,8 +51,7 @@ class GCV { using This = GCV; using VectorType = DVector; using MatrixType = DMatrix; - using ModelType = RegressionView; - ModelType model_; // model to calibrate + RegressionView model_; // model to calibrate EDFStrategy trS_; // strategy used to evaluate the trace of smoothing matrix S std::vector edfs_; // equivalent degrees of freedom q + Tr[S] std::vector gcvs_; // computed values of GCV index diff --git a/fdaPDE/models/regression/gsrpde.h b/fdaPDE/models/regression/gsrpde.h index 34c55ff6..2bddf52c 100644 --- a/fdaPDE/models/regression/gsrpde.h +++ b/fdaPDE/models/regression/gsrpde.h @@ -30,19 +30,9 @@ namespace models { // base class for GSRPDE model template class GSRPDE : public RegressionBase, RegularizationType_> { - private: - Distribution distr_ {}; - DVector py_; // \tilde y^k = G^k(y-u^k) + \theta^k - DVector pW_; // diagonal of W^k = ((G^k)^{-2})*((V^k)^{-1}) - DVector mu_; // \mu^k = [ \mu^k_1, ..., \mu^k_n ] : mean vector at step k - DMatrix T_; // T = \Psi^T*Q*\Psi + P - fdapde::SparseLU> invA_; // factorization of non-parametric system matrix A - - // FPIRLS parameters (set to default) - std::size_t max_iter_ = 200; - double tol_ = 1e-4; public: - using RegularizationType = RegularizationType_; + using RegularizationType = std::decay_t; + using This = GSRPDE; using Base = RegressionBase, RegularizationType>; // import commonly defined symbols from base IMPORT_REGRESSION_SYMBOLS; @@ -57,37 +47,40 @@ class GSRPDE : public RegressionBase, Regularization template < typename U = RegularizationType, typename std::enable_if::value, int>::type = 0> - GSRPDE(const pde_ptr& pde, Sampling s, const Distribution& distr) : Base(pde, s), distr_(distr) {}; + GSRPDE(const pde_ptr& pde, Sampling s, const Distribution& distr) : Base(pde, s), distr_(distr) { + fpirls_ = FPIRLS(this, tol_, max_iter_); + }; // space-time constructor template < typename U = RegularizationType, typename std::enable_if::value, int>::type = 0> GSRPDE(const pde_ptr& pde, const DVector& time, Sampling s, const Distribution& distr) : - Base(pde, s, time), distr_(distr) {}; + Base(pde, s, time), distr_(distr) { + fpirls_ = FPIRLS(this, tol_, max_iter_); + }; // setters void set_fpirls_tolerance(double tol) { tol_ = tol; } void set_fpirls_max_iter(std::size_t max_iter) { max_iter_ = max_iter; } void init_data() { return; }; - void init_model() { return; }; + void init_model() { fpirls_.init(); }; void solve() { // finds a solution to the smoothing problem // execute FPIRLS for minimization of functional \norm{V^{-1/2}(y - \mu)}^2 + \lambda \int_D (Lf - u)^2 - FPIRLS fpirls(*this, tol_, max_iter_); // FPIRLS engine - fpirls.compute(); - // fpirls converged: extract matrix W and solution estimates - W_ = fpirls.solver().W(); - f_ = fpirls.solver().f(); - g_ = fpirls.solver().g(); + fpirls_.compute(); + // fpirls_ converged: extract matrix W and solution estimates + W_ = fpirls_.solver().W(); + f_ = fpirls_.solver().f(); + g_ = fpirls_.solver().g(); // store parametric part if (has_covariates()) { - beta_ = fpirls.solver().beta(); - XtWX_ = fpirls.solver().XtWX(); - invXtWX_ = fpirls.solver().invXtWX(); - U_ = fpirls.solver().U(); - V_ = fpirls.solver().V(); + beta_ = fpirls_.solver().beta(); + XtWX_ = fpirls_.solver().XtWX(); + invXtWX_ = fpirls_.solver().invXtWX(); + U_ = fpirls_.solver().U(); + V_ = fpirls_.solver().V(); } - invA_ = fpirls.solver().invA(); + invA_ = fpirls_.solver().invA(); return; } // required by FPIRLS (see fpirls.h for details) @@ -123,6 +116,18 @@ class GSRPDE : public RegressionBase, Regularization } virtual ~GSRPDE() = default; + private: + Distribution distr_ {}; + DVector py_; // \tilde y^k = G^k(y-u^k) + \theta^k + DVector pW_; // diagonal of W^k = ((G^k)^{-2})*((V^k)^{-1}) + DVector mu_; // \mu^k = [ \mu^k_1, ..., \mu^k_n ] : mean vector at step k + DMatrix T_; // T = \Psi^T*Q*\Psi + P + fdapde::SparseLU> invA_; // factorization of non-parametric system matrix A + + // FPIRLS parameters (set to default) + FPIRLS fpirls_; + std::size_t max_iter_ = 200; + double tol_ = 1e-4; }; } // namespace models diff --git a/fdaPDE/models/regression/sqrpde.h b/fdaPDE/models/regression/qsrpde.h similarity index 76% rename from fdaPDE/models/regression/sqrpde.h rename to fdaPDE/models/regression/qsrpde.h index 03f4ba83..e4776947 100644 --- a/fdaPDE/models/regression/sqrpde.h +++ b/fdaPDE/models/regression/qsrpde.h @@ -14,8 +14,8 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#ifndef __SQRPDE_H__ -#define __SQRPDE_H__ +#ifndef __QSRPDE_H__ +#define __QSRPDE_H__ #include #include @@ -28,24 +28,11 @@ namespace fdapde { namespace models { template -class SQRPDE : public RegressionBase, RegularizationType_> { - private: - double alpha_; // quantile order - DVector py_ {}; // y - (1-2*alpha)|y - X*beta - f| - DVector pW_ {}; // diagonal of W^k = 1/(2*n*|y - X*beta - f|) - DVector mu_; // \mu^k = [ \mu^k_1, ..., \mu^k_n ] : mean vector at step k - DMatrix T_; // T = \Psi^T*Q*\Psi + P - fdapde::SparseLU> invA_; // factorization of non-parametric system matrix A - - // FPIRLS parameters (set to default) - std::size_t max_iter_ = 200; - double tol_weights_ = 1e-6; - double tol_ = 1e-6; - - double pinball_loss(double x) const { return 0.5 * std::abs(x) + (alpha_ - 0.5) * x; }; // quantile check function +class QSRPDE : public RegressionBase, RegularizationType_> { public: - using RegularizationType = RegularizationType_; - using Base = RegressionBase, RegularizationType>; + using RegularizationType = std::decay_t; + using This = QSRPDE; + using Base = RegressionBase, RegularizationType>; // import commonly defined symbols from base IMPORT_REGRESSION_SYMBOLS; using Base::invXtWX_; // LU factorization of X^T*W*X @@ -55,18 +42,22 @@ class SQRPDE : public RegressionBase, Regularization using Base::W_; // weight matrix using Base::XtWX_; // q x q matrix X^T*W*X // constructor - SQRPDE() = default; + QSRPDE() = default; // space-only constructor template < typename U = RegularizationType, typename std::enable_if::value, int>::type = 0> - SQRPDE(const pde_ptr& pde, Sampling s, double alpha) : Base(pde, s), alpha_(alpha) {}; + QSRPDE(const pde_ptr& pde, Sampling s, double alpha = 0.5) : Base(pde, s), alpha_(alpha) { + fpirls_ = FPIRLS(this, tol_, max_iter_); + }; // space-time constructor template < typename U = RegularizationType, typename std::enable_if::value, int>::type = 0> - SQRPDE(const pde_ptr& pde, const DVector& time, Sampling s, double alpha) : - Base(pde, s, time), alpha_(alpha) {}; + QSRPDE(const pde_ptr& pde, const DVector& time, Sampling s, double alpha = 0.5) : + Base(pde, s, time), alpha_(alpha) { + fpirls_ = FPIRLS(this, tol_, max_iter_); + }; // setter void set_fpirls_tolerance(double tol) { tol_ = tol; } @@ -74,28 +65,27 @@ class SQRPDE : public RegressionBase, Regularization void set_alpha(double alpha) { alpha_ = alpha; } void init_data() { return; } - void init_model() { return; } + void init_model() { fpirls_.init(); } void solve() { // finds a solution to the smoothing problem - // execute FPIRLS for minimization of functional \norm{V^{-1/2}(y - \mu)}^2 + \lambda \int_D (Lf - u)^2 - FPIRLS fpirls(*this, tol_, max_iter_); // FPIRLS engine - fpirls.compute(); - // fpirls converged: extract matrix W and solution estimates - W_ = fpirls.solver().W(); - f_ = fpirls.solver().f(); - g_ = fpirls.solver().g(); + // execute FPIRLS_ for minimization of functional \norm{V^{-1/2}(y - \mu)}^2 + \lambda \int_D (Lf - u)^2 + fpirls_.compute(); + // fpirls_ converged: extract matrix W and solution estimates + W_ = fpirls_.solver().W(); + f_ = fpirls_.solver().f(); + g_ = fpirls_.solver().g(); // store parametric part if (has_covariates()) { - beta_ = fpirls.solver().beta(); - XtWX_ = fpirls.solver().XtWX(); - invXtWX_ = fpirls.solver().invXtWX(); - U_ = fpirls.solver().U(); - V_ = fpirls.solver().V(); + beta_ = fpirls_.solver().beta(); + XtWX_ = fpirls_.solver().XtWX(); + invXtWX_ = fpirls_.solver().invXtWX(); + U_ = fpirls_.solver().U(); + V_ = fpirls_.solver().V(); } - invA_ = fpirls.solver().invA(); + invA_ = fpirls_.solver().invA(); return; } - // required by FPIRLS (see fpirls.h for details) + // required by FPIRLS_ (see fpirls_.h for details) // initalizes mean vector \mu void fpirls_init() { // non-parametric and semi-parametric cases coincide here, since beta^(0) = 0 @@ -109,7 +99,6 @@ class SQRPDE : public RegressionBase, Regularization DVector b(A.rows()); b.block(n_basis(), 0, n_basis(), 1) = lambda_D() * u(); b.block(0, 0, n_basis(), 1) = PsiTD() * y() / n_obs(); - mu_ = Psi(not_nan()) * (invA.solve(b)).head(n_basis()); } // computes W^k = diag(1/(2*n*|y - X*beta - f|)) and y^k = y - (1-2*alpha)|y - X*beta - f| @@ -137,10 +126,25 @@ class SQRPDE : public RegressionBase, Regularization return result * result / n_obs(); } - virtual ~SQRPDE() = default; + virtual ~QSRPDE() = default; + private: + double alpha_ = 0.5; // quantile order (default to median) + DVector py_ {}; // y - (1-2*alpha)|y - X*beta - f| + DVector pW_ {}; // diagonal of W^k = 1/(2*n*|y - X*beta - f|) + DVector mu_; // \mu^k = [ \mu^k_1, ..., \mu^k_n ] : mean vector at step k + DMatrix T_; // T = \Psi^T*Q*\Psi + P + fdapde::SparseLU> invA_; // factorization of non-parametric system matrix A + + // FPIRLS algorithm + FPIRLS fpirls_; + std::size_t max_iter_ = 200; + double tol_weights_ = 1e-6; + double tol_ = 1e-6; + + double pinball_loss(double x) const { return 0.5 * std::abs(x) + (alpha_ - 0.5) * x; }; // quantile check function }; } // namespace models } // namespace fdapde -#endif // __SQRPDE_H__ +#endif // __QSRPDE_H__ diff --git a/test/data/gcv/sqrpde/2D_test1/edfs.mtx b/test/data/gcv/qsrpde/2D_test1/edfs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test1/edfs.mtx rename to test/data/gcv/qsrpde/2D_test1/edfs.mtx diff --git a/test/data/gcv/sqrpde/2D_test1/gcvs.mtx b/test/data/gcv/qsrpde/2D_test1/gcvs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test1/gcvs.mtx rename to test/data/gcv/qsrpde/2D_test1/gcvs.mtx diff --git a/test/data/gcv/sqrpde/2D_test1/y.csv b/test/data/gcv/qsrpde/2D_test1/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test1/y.csv rename to test/data/gcv/qsrpde/2D_test1/y.csv diff --git a/test/data/gcv/qsrpde/2D_test2/edfs.mtx b/test/data/gcv/qsrpde/2D_test2/edfs.mtx new file mode 100644 index 00000000..c1a5acaf --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test2/edfs.mtx @@ -0,0 +1,15 @@ +%%MatrixMarket matrix coordinate real general +13 1 13 +1 1 4.40181065106688493e+02 +2 1 4.40461549188092931e+02 +3 1 4.37525587313592382e+02 +4 1 4.29211236257885275e+02 +5 1 4.13343713146583411e+02 +6 1 3.86252486482486859e+02 +7 1 3.48096097528348707e+02 +8 1 3.00081914257956271e+02 +9 1 2.40458170291202066e+02 +10 1 1.89289651455957738e+02 +11 1 1.47149667350187968e+02 +12 1 1.08751825437762420e+02 +13 1 8.19879784755787995e+01 diff --git a/test/data/gcv/qsrpde/2D_test2/gcvs.mtx b/test/data/gcv/qsrpde/2D_test2/gcvs.mtx new file mode 100644 index 00000000..502ef0d8 --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test2/gcvs.mtx @@ -0,0 +1,15 @@ +%%MatrixMarket matrix coordinate real general +13 1 13 +1 1 1.85022336471077103e-01 +2 1 8.93634967030510485e-02 +3 1 3.41513787709415706e-02 +4 1 1.43211995579107850e-02 +5 1 6.61563999179170734e-03 +6 1 3.60150915516095019e-03 +7 1 2.32787435564876452e-03 +8 1 1.69191334985671236e-03 +9 1 1.59045816914823437e-03 +10 1 1.85863465303526504e-03 +11 1 2.34318846009063612e-03 +12 1 3.29667509011485720e-03 +13 1 4.36067393625100101e-03 diff --git a/test/data/gcv/sqrpde/2D_test2/y.csv b/test/data/gcv/qsrpde/2D_test2/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test2/y.csv rename to test/data/gcv/qsrpde/2D_test2/y.csv diff --git a/test/data/gcv/sqrpde/2D_test3/X.csv b/test/data/gcv/qsrpde/2D_test3/X.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test3/X.csv rename to test/data/gcv/qsrpde/2D_test3/X.csv diff --git a/test/data/gcv/sqrpde/2D_test3/edfs.mtx b/test/data/gcv/qsrpde/2D_test3/edfs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test3/edfs.mtx rename to test/data/gcv/qsrpde/2D_test3/edfs.mtx diff --git a/test/data/gcv/sqrpde/2D_test3/gcvs.mtx b/test/data/gcv/qsrpde/2D_test3/gcvs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test3/gcvs.mtx rename to test/data/gcv/qsrpde/2D_test3/gcvs.mtx diff --git a/test/data/gcv/sqrpde/2D_test3/locs.csv b/test/data/gcv/qsrpde/2D_test3/locs.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test3/locs.csv rename to test/data/gcv/qsrpde/2D_test3/locs.csv diff --git a/test/data/gcv/sqrpde/2D_test3/y.csv b/test/data/gcv/qsrpde/2D_test3/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test3/y.csv rename to test/data/gcv/qsrpde/2D_test3/y.csv diff --git a/test/data/gcv/sqrpde/2D_test4/X.csv b/test/data/gcv/qsrpde/2D_test4/X.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test4/X.csv rename to test/data/gcv/qsrpde/2D_test4/X.csv diff --git a/test/data/gcv/qsrpde/2D_test4/edfs.mtx b/test/data/gcv/qsrpde/2D_test4/edfs.mtx new file mode 100644 index 00000000..34e19e1b --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test4/edfs.mtx @@ -0,0 +1,11 @@ +%%MatrixMarket matrix coordinate real general +9 1 9 +1 1 1.33208451751103894e+02 +2 1 1.15822163388498851e+02 +3 1 9.71434238109929993e+01 +4 1 8.09580592603756912e+01 +5 1 6.71748110324557643e+01 +6 1 5.57170344101969306e+01 +7 1 4.51604587862486753e+01 +8 1 3.66042264945400362e+01 +9 1 3.05617516363425814e+01 diff --git a/test/data/gcv/qsrpde/2D_test4/gcvs.mtx b/test/data/gcv/qsrpde/2D_test4/gcvs.mtx new file mode 100644 index 00000000..601a35cf --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test4/gcvs.mtx @@ -0,0 +1,11 @@ +%%MatrixMarket matrix coordinate real general +9 1 9 +1 1 1.10713333445202709e-02 +2 1 9.53548301273807186e-03 +3 1 9.01831639778432302e-03 +4 1 8.97841892566965567e-03 +5 1 8.74284423052914574e-03 +6 1 8.50055053315700894e-03 +7 1 8.24363093986560554e-03 +8 1 8.29946782863691217e-03 +9 1 8.47340451288999125e-03 diff --git a/test/data/gcv/sqrpde/2D_test4/locs.csv b/test/data/gcv/qsrpde/2D_test4/locs.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test4/locs.csv rename to test/data/gcv/qsrpde/2D_test4/locs.csv diff --git a/test/data/gcv/sqrpde/2D_test4/y.csv b/test/data/gcv/qsrpde/2D_test4/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test4/y.csv rename to test/data/gcv/qsrpde/2D_test4/y.csv diff --git a/test/data/gcv/sqrpde/2D_test5/edfs.mtx b/test/data/gcv/qsrpde/2D_test5/edfs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test5/edfs.mtx rename to test/data/gcv/qsrpde/2D_test5/edfs.mtx diff --git a/test/data/gcv/sqrpde/2D_test5/gcvs.mtx b/test/data/gcv/qsrpde/2D_test5/gcvs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test5/gcvs.mtx rename to test/data/gcv/qsrpde/2D_test5/gcvs.mtx diff --git a/test/data/gcv/sqrpde/2D_test5/y.csv b/test/data/gcv/qsrpde/2D_test5/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test5/y.csv rename to test/data/gcv/qsrpde/2D_test5/y.csv diff --git a/test/data/gcv/qsrpde/2D_test6/edfs.mtx b/test/data/gcv/qsrpde/2D_test6/edfs.mtx new file mode 100644 index 00000000..413371d4 --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test6/edfs.mtx @@ -0,0 +1,11 @@ +%%MatrixMarket matrix coordinate real general +9 1 9 +1 1 2.76108027251456861e+02 +2 1 2.34846489626921681e+02 +3 1 1.98168466917270962e+02 +4 1 1.64992984514900172e+02 +5 1 1.36254920696387046e+02 +6 1 1.07309453428111595e+02 +7 1 7.98421569650793970e+01 +8 1 4.99539189703289566e+01 +9 1 2.74748237704313887e+01 diff --git a/test/data/gcv/qsrpde/2D_test6/gcvs.mtx b/test/data/gcv/qsrpde/2D_test6/gcvs.mtx new file mode 100644 index 00000000..b594fe31 --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test6/gcvs.mtx @@ -0,0 +1,11 @@ +%%MatrixMarket matrix coordinate real general +9 1 9 +1 1 3.04111220611080410e-03 +2 1 2.58017530185898018e-03 +3 1 2.33049696172554733e-03 +4 1 2.19660369381580439e-03 +5 1 2.12397845206960126e-03 +6 1 2.16475571098232594e-03 +7 1 3.12599572828926690e-03 +8 1 8.22940145188545344e-03 +9 1 3.54877353443419086e-02 diff --git a/test/data/gcv/sqrpde/2D_test6/y.csv b/test/data/gcv/qsrpde/2D_test6/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test6/y.csv rename to test/data/gcv/qsrpde/2D_test6/y.csv diff --git a/test/data/gcv/sqrpde/2D_test7/X.csv b/test/data/gcv/qsrpde/2D_test7/X.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test7/X.csv rename to test/data/gcv/qsrpde/2D_test7/X.csv diff --git a/test/data/gcv/sqrpde/2D_test7/edfs.mtx b/test/data/gcv/qsrpde/2D_test7/edfs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test7/edfs.mtx rename to test/data/gcv/qsrpde/2D_test7/edfs.mtx diff --git a/test/data/gcv/sqrpde/2D_test7/gcvs.mtx b/test/data/gcv/qsrpde/2D_test7/gcvs.mtx similarity index 100% rename from test/data/gcv/sqrpde/2D_test7/gcvs.mtx rename to test/data/gcv/qsrpde/2D_test7/gcvs.mtx diff --git a/test/data/gcv/sqrpde/2D_test7/incidence_matrix.csv b/test/data/gcv/qsrpde/2D_test7/incidence_matrix.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test7/incidence_matrix.csv rename to test/data/gcv/qsrpde/2D_test7/incidence_matrix.csv diff --git a/test/data/gcv/sqrpde/2D_test7/y.csv b/test/data/gcv/qsrpde/2D_test7/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test7/y.csv rename to test/data/gcv/qsrpde/2D_test7/y.csv diff --git a/test/data/gcv/sqrpde/2D_test8/X.csv b/test/data/gcv/qsrpde/2D_test8/X.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test8/X.csv rename to test/data/gcv/qsrpde/2D_test8/X.csv diff --git a/test/data/gcv/qsrpde/2D_test8/edfs.mtx b/test/data/gcv/qsrpde/2D_test8/edfs.mtx new file mode 100644 index 00000000..1c76705e --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test8/edfs.mtx @@ -0,0 +1,15 @@ +%%MatrixMarket matrix coordinate real general +13 1 13 +1 1 1.60332053743644067e+01 +2 1 1.50556739023827202e+01 +3 1 1.24187428657688930e+01 +4 1 1.11837428172403630e+01 +5 1 9.38379037069958422e+00 +6 1 9.14477475144124163e+00 +7 1 7.93557228806301840e+00 +8 1 7.77418029620436002e+00 +9 1 7.33261570671088059e+00 +10 1 5.49635230772290218e+00 +11 1 5.57295577546898713e+00 +12 1 4.60232261412499266e+00 +13 1 3.65935263364707009e+00 diff --git a/test/data/gcv/qsrpde/2D_test8/gcvs.mtx b/test/data/gcv/qsrpde/2D_test8/gcvs.mtx new file mode 100644 index 00000000..549a190d --- /dev/null +++ b/test/data/gcv/qsrpde/2D_test8/gcvs.mtx @@ -0,0 +1,15 @@ +%%MatrixMarket matrix coordinate real general +13 1 13 +1 1 9.10963597892400362e-02 +2 1 7.54397843503786109e-02 +3 1 4.64540641459125006e-02 +4 1 4.16522790107617669e-02 +5 1 3.29708740124532440e-02 +6 1 3.33889836137126125e-02 +7 1 2.93213549610455149e-02 +8 1 3.10768381396734072e-02 +9 1 3.71929918656311170e-02 +10 1 5.34641262706091677e-02 +11 1 1.16737447429273444e-01 +12 1 2.98633657325982516e-01 +13 1 6.22200203224899329e-01 diff --git a/test/data/gcv/sqrpde/2D_test8/incidence_matrix.csv b/test/data/gcv/qsrpde/2D_test8/incidence_matrix.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test8/incidence_matrix.csv rename to test/data/gcv/qsrpde/2D_test8/incidence_matrix.csv diff --git a/test/data/gcv/sqrpde/2D_test8/y.csv b/test/data/gcv/qsrpde/2D_test8/y.csv similarity index 100% rename from test/data/gcv/sqrpde/2D_test8/y.csv rename to test/data/gcv/qsrpde/2D_test8/y.csv diff --git a/test/data/gcv/sqrpde/2D_test2/edfs.mtx b/test/data/gcv/sqrpde/2D_test2/edfs.mtx deleted file mode 100644 index 58c7f78e..00000000 --- a/test/data/gcv/sqrpde/2D_test2/edfs.mtx +++ /dev/null @@ -1,15 +0,0 @@ -%%MatrixMarket matrix coordinate real general -13 1 13 -1 1 4.40183346644801702e+02 -2 1 4.40460206573702180e+02 -3 1 4.37535107647263715e+02 -4 1 4.29261291999342177e+02 -5 1 4.13411875955585913e+02 -6 1 3.86343040753663786e+02 -7 1 3.48260312274547175e+02 -8 1 3.00276203324444225e+02 -9 1 2.40735296585845902e+02 -10 1 1.89706108703869290e+02 -11 1 1.47757844233170005e+02 -12 1 1.09359221752539497e+02 -13 1 8.22482852740177606e+01 diff --git a/test/data/gcv/sqrpde/2D_test2/gcvs.mtx b/test/data/gcv/sqrpde/2D_test2/gcvs.mtx deleted file mode 100644 index 7d970303..00000000 --- a/test/data/gcv/sqrpde/2D_test2/gcvs.mtx +++ /dev/null @@ -1,15 +0,0 @@ -%%MatrixMarket matrix coordinate real general -13 1 13 -1 1 1.86057598699536697e-01 -2 1 8.89195064169029242e-02 -3 1 3.43393091485468885e-02 -4 1 1.44435957709837595e-02 -5 1 6.64837130319134370e-03 -6 1 3.61345281843854097e-03 -7 1 2.33612561794895686e-03 -8 1 1.69658842492572304e-03 -9 1 1.59486296635274810e-03 -10 1 1.86480020931401708e-03 -11 1 2.35291796717036414e-03 -12 1 3.30876179449441418e-03 -13 1 4.36700435832937210e-03 diff --git a/test/data/gcv/sqrpde/2D_test4/edfs.mtx b/test/data/gcv/sqrpde/2D_test4/edfs.mtx deleted file mode 100644 index a7033352..00000000 --- a/test/data/gcv/sqrpde/2D_test4/edfs.mtx +++ /dev/null @@ -1,11 +0,0 @@ -%%MatrixMarket matrix coordinate real general -9 1 9 -1 1 1.32843601726190201e+02 -2 1 1.15499362149987803e+02 -3 1 9.68845493564370202e+01 -4 1 8.06846136603180639e+01 -5 1 6.68050887683269252e+01 -6 1 5.53398533185044670e+01 -7 1 4.47615407563435070e+01 -8 1 3.62331815714342511e+01 -9 1 3.02718640585923318e+01 diff --git a/test/data/gcv/sqrpde/2D_test4/gcvs.mtx b/test/data/gcv/sqrpde/2D_test4/gcvs.mtx deleted file mode 100644 index c70bec5e..00000000 --- a/test/data/gcv/sqrpde/2D_test4/gcvs.mtx +++ /dev/null @@ -1,11 +0,0 @@ -%%MatrixMarket matrix coordinate real general -9 1 9 -1 1 1.09476903560507197e-02 -2 1 9.46100978950794799e-03 -3 1 8.97219835557417841e-03 -4 1 8.93661285456918628e-03 -5 1 8.69363699437235640e-03 -6 1 8.45566051039961346e-03 -7 1 8.20076634170631945e-03 -8 1 8.26143847726866325e-03 -9 1 8.44414034530201370e-03 diff --git a/test/data/gcv/sqrpde/2D_test6/edfs.mtx b/test/data/gcv/sqrpde/2D_test6/edfs.mtx deleted file mode 100644 index 45cd975c..00000000 --- a/test/data/gcv/sqrpde/2D_test6/edfs.mtx +++ /dev/null @@ -1,11 +0,0 @@ -%%MatrixMarket matrix coordinate real general -9 1 9 -1 1 2.76660287538494174e+02 -2 1 2.35446691970586812e+02 -3 1 1.98666282041430208e+02 -4 1 1.65328463942675398e+02 -5 1 1.36269310132865797e+02 -6 1 1.07003073047535693e+02 -7 1 7.94409455692377833e+01 -8 1 5.00465100220374168e+01 -9 1 2.71044020764703610e+01 diff --git a/test/data/gcv/sqrpde/2D_test6/gcvs.mtx b/test/data/gcv/sqrpde/2D_test6/gcvs.mtx deleted file mode 100644 index fd08e6a7..00000000 --- a/test/data/gcv/sqrpde/2D_test6/gcvs.mtx +++ /dev/null @@ -1,11 +0,0 @@ -%%MatrixMarket matrix coordinate real general -9 1 9 -1 1 3.06158602840628192e-03 -2 1 2.59526528889497806e-03 -3 1 2.34008165215113598e-03 -4 1 2.20195327326907402e-03 -5 1 2.12417904098732412e-03 -6 1 2.16078600707469200e-03 -7 1 3.11906192918164115e-03 -8 1 8.23329991297972255e-03 -9 1 3.54242432773061930e-02 diff --git a/test/data/gcv/sqrpde/2D_test8/edfs.mtx b/test/data/gcv/sqrpde/2D_test8/edfs.mtx deleted file mode 100644 index f5186f1d..00000000 --- a/test/data/gcv/sqrpde/2D_test8/edfs.mtx +++ /dev/null @@ -1,15 +0,0 @@ -%%MatrixMarket matrix coordinate real general -13 1 13 -1 1 1.56759797446367308e+01 -2 1 1.47286602265912805e+01 -3 1 1.23351484865424208e+01 -4 1 1.11527593301017696e+01 -5 1 9.44449537476370615e+00 -6 1 9.14945671776664327e+00 -7 1 7.44364745004328388e+00 -8 1 7.09696963194546626e+00 -9 1 6.63828884096076433e+00 -10 1 5.16789224208182940e+00 -11 1 4.02199545049469975e+00 -12 1 4.69818139428146075e+00 -13 1 3.34854806941988414e+00 diff --git a/test/data/gcv/sqrpde/2D_test8/gcvs.mtx b/test/data/gcv/sqrpde/2D_test8/gcvs.mtx deleted file mode 100644 index 3c3e7c5d..00000000 --- a/test/data/gcv/sqrpde/2D_test8/gcvs.mtx +++ /dev/null @@ -1,15 +0,0 @@ -%%MatrixMarket matrix coordinate real general -13 1 13 -1 1 7.66663903215633602e-02 -2 1 6.63701242216443860e-02 -3 1 4.54463152958035008e-02 -4 1 4.13610530476003282e-02 -5 1 3.33511973492155583e-02 -6 1 3.34178042578959431e-02 -7 1 2.70688923193140689e-02 -8 1 2.79003313299180898e-02 -9 1 3.34280350226750303e-02 -10 1 5.11223975545720688e-02 -11 1 9.51743321564772343e-02 -12 1 3.02386979096593178e-01 -13 1 5.99189850752087083e-01 diff --git a/test/data/models/sqrpde/2D_test1/README.md b/test/data/models/qsrpde/2D_test1/README.md similarity index 100% rename from test/data/models/sqrpde/2D_test1/README.md rename to test/data/models/qsrpde/2D_test1/README.md diff --git a/test/data/models/sqrpde/2D_test1/sol.mtx b/test/data/models/qsrpde/2D_test1/sol.mtx similarity index 100% rename from test/data/models/sqrpde/2D_test1/sol.mtx rename to test/data/models/qsrpde/2D_test1/sol.mtx diff --git a/test/data/models/sqrpde/2D_test1/y.csv b/test/data/models/qsrpde/2D_test1/y.csv similarity index 100% rename from test/data/models/sqrpde/2D_test1/y.csv rename to test/data/models/qsrpde/2D_test1/y.csv diff --git a/test/data/models/sqrpde/2D_test2/README.md b/test/data/models/qsrpde/2D_test2/README.md similarity index 100% rename from test/data/models/sqrpde/2D_test2/README.md rename to test/data/models/qsrpde/2D_test2/README.md diff --git a/test/data/models/sqrpde/2D_test2/X.csv b/test/data/models/qsrpde/2D_test2/X.csv similarity index 100% rename from test/data/models/sqrpde/2D_test2/X.csv rename to test/data/models/qsrpde/2D_test2/X.csv diff --git a/test/data/models/sqrpde/2D_test2/beta.mtx b/test/data/models/qsrpde/2D_test2/beta.mtx similarity index 100% rename from test/data/models/sqrpde/2D_test2/beta.mtx rename to test/data/models/qsrpde/2D_test2/beta.mtx diff --git a/test/data/models/sqrpde/2D_test2/locs.csv b/test/data/models/qsrpde/2D_test2/locs.csv similarity index 100% rename from test/data/models/sqrpde/2D_test2/locs.csv rename to test/data/models/qsrpde/2D_test2/locs.csv diff --git a/test/data/models/sqrpde/2D_test2/sol.mtx b/test/data/models/qsrpde/2D_test2/sol.mtx similarity index 100% rename from test/data/models/sqrpde/2D_test2/sol.mtx rename to test/data/models/qsrpde/2D_test2/sol.mtx diff --git a/test/data/models/sqrpde/2D_test2/y.csv b/test/data/models/qsrpde/2D_test2/y.csv similarity index 100% rename from test/data/models/sqrpde/2D_test2/y.csv rename to test/data/models/qsrpde/2D_test2/y.csv diff --git a/test/data/models/sqrpde/2D_test3/README.md b/test/data/models/qsrpde/2D_test3/README.md similarity index 100% rename from test/data/models/sqrpde/2D_test3/README.md rename to test/data/models/qsrpde/2D_test3/README.md diff --git a/test/data/models/sqrpde/2D_test3/sol.mtx b/test/data/models/qsrpde/2D_test3/sol.mtx similarity index 100% rename from test/data/models/sqrpde/2D_test3/sol.mtx rename to test/data/models/qsrpde/2D_test3/sol.mtx diff --git a/test/data/models/sqrpde/2D_test3/y.csv b/test/data/models/qsrpde/2D_test3/y.csv similarity index 100% rename from test/data/models/sqrpde/2D_test3/y.csv rename to test/data/models/qsrpde/2D_test3/y.csv diff --git a/test/data/models/sqrpde/2D_test4/README.md b/test/data/models/qsrpde/2D_test4/README.md similarity index 100% rename from test/data/models/sqrpde/2D_test4/README.md rename to test/data/models/qsrpde/2D_test4/README.md diff --git a/test/data/models/sqrpde/2D_test4/X.csv b/test/data/models/qsrpde/2D_test4/X.csv similarity index 100% rename from test/data/models/sqrpde/2D_test4/X.csv rename to test/data/models/qsrpde/2D_test4/X.csv diff --git a/test/data/models/sqrpde/2D_test4/beta.mtx b/test/data/models/qsrpde/2D_test4/beta.mtx similarity index 100% rename from test/data/models/sqrpde/2D_test4/beta.mtx rename to test/data/models/qsrpde/2D_test4/beta.mtx diff --git a/test/data/models/sqrpde/2D_test4/incidence_matrix.csv b/test/data/models/qsrpde/2D_test4/incidence_matrix.csv similarity index 100% rename from test/data/models/sqrpde/2D_test4/incidence_matrix.csv rename to test/data/models/qsrpde/2D_test4/incidence_matrix.csv diff --git a/test/data/models/sqrpde/2D_test4/sol.mtx b/test/data/models/qsrpde/2D_test4/sol.mtx similarity index 100% rename from test/data/models/sqrpde/2D_test4/sol.mtx rename to test/data/models/qsrpde/2D_test4/sol.mtx diff --git a/test/data/models/sqrpde/2D_test4/y.csv b/test/data/models/qsrpde/2D_test4/y.csv similarity index 100% rename from test/data/models/sqrpde/2D_test4/y.csv rename to test/data/models/qsrpde/2D_test4/y.csv diff --git a/test/main.cpp b/test/main.cpp index 7426c978..c9c1b95c 100644 --- a/test/main.cpp +++ b/test/main.cpp @@ -7,11 +7,11 @@ #include "src/srpde_test.cpp" #include "src/strpde_test.cpp" #include "src/gsrpde_test.cpp" -#include "src/sqrpde_test.cpp" +#include "src/qsrpde_test.cpp" // #include "src/fpca_test.cpp" // GCV test suites #include "src/gcv_srpde_test.cpp" -// #include "src/gcv_sqrpde_test.cpp" +#include "src/gcv_qsrpde_test.cpp" #include "src/gcv_srpde_newton_test.cpp" int main(int argc, char **argv){ diff --git a/test/src/gcv_sqrpde_test.cpp b/test/src/gcv_qsrpde_test.cpp similarity index 73% rename from test/src/gcv_sqrpde_test.cpp rename to test/src/gcv_qsrpde_test.cpp index 912f059a..95e84ce4 100644 --- a/test/src/gcv_sqrpde_test.cpp +++ b/test/src/gcv_qsrpde_test.cpp @@ -28,18 +28,15 @@ using fdapde::core::MatrixDataWrapper; using fdapde::core::PDE; using fdapde::core::VectorDataWrapper; -#include "../../fdaPDE/models/regression/sqrpde.h" +#include "../../fdaPDE/models/regression/qsrpde.h" #include "../../fdaPDE/models/regression/gcv.h" #include "../../fdaPDE/models/sampling_design.h" -using fdapde::models::Areal; -using fdapde::models::GeoStatLocations; -using fdapde::models::GeoStatMeshNodes; -using fdapde::models::SQRPDE; +using fdapde::models::QSRPDE; using fdapde::models::SpaceOnly; -using fdapde::models::MonolithicSolver; using fdapde::models::ExactEDF; using fdapde::models::GCV; using fdapde::models::StochasticEDF; +using fdapde::models::Sampling; #include "utils/constants.h" #include "utils/mesh_loader.h" @@ -57,18 +54,18 @@ using fdapde::testing::read_csv; // BC: no // order FE: 1 // GCV optimization: grid exact -TEST(gcv_sqrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridexact) { +TEST(gcv_qsrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridexact) { // define domain MeshLoader domain("unit_square_coarse"); // import data from files - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test1/y.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test1/y.csv"); // 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 double alpha = 0.1; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::mesh_nodes, alpha); // set model's data BlockFrame df; df.insert(OBSERVATIONS_BLK, y); @@ -76,14 +73,14 @@ TEST(gcv_sqrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridexac model.init(); // define GCV function and grid of \lambda_D values auto GCV = model.gcv(); - std::vector> lambdas; + std::vector> lambdas; for (double x = -8.0; x <= -5.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test1/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test1/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test1/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test1/gcvs.mtx")); } // test 2 @@ -94,18 +91,18 @@ TEST(gcv_sqrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridexac // BC: no // order FE: 1 // GCV optimization: grid stochastic -TEST(gcv_sqrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridstochastic) { +TEST(gcv_qsrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridstochastic) { // define domain MeshLoader domain("unit_square_coarse"); // import data from files - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test2/y.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test2/y.csv"); // 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 double alpha = 0.1; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::mesh_nodes, alpha); // set model's data BlockFrame df; df.insert(OBSERVATIONS_BLK, y); @@ -114,14 +111,14 @@ TEST(gcv_sqrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridstoc // define GCV function and grid of \lambda_D values std::size_t seed = 438172; auto GCV = model.gcv(1000, seed); - std::vector> lambdas; + std::vector> lambdas; for (double x = -8.0; x <= -5.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test2/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test2/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test2/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test2/gcvs.mtx")); } // test 3 @@ -132,20 +129,20 @@ TEST(gcv_sqrpde_test, laplacian_nonparametric_samplingatnodes_spaceonly_gridstoc // BC: no // order FE: 1 // GCV optimization: grid exact -TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingatlocations_gridexact) { +TEST(gcv_qsrpde_test, laplacian_semiparametric_samplingatlocations_gridexact) { // define domain MeshLoader domain("c_shaped"); // import data from files - DMatrix locs = read_csv("../data/gcv/sqrpde/2D_test3/locs.csv"); - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test3/y.csv"); - DMatrix X = read_csv("../data/gcv/sqrpde/2D_test3/X.csv"); + DMatrix locs = read_csv("../data/gcv/qsrpde/2D_test3/locs.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test3/y.csv"); + DMatrix X = read_csv("../data/gcv/qsrpde/2D_test3/X.csv"); // 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 double alpha = 0.9; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::pointwise, alpha); model.set_spatial_locations(locs); // set model's data BlockFrame df; @@ -155,14 +152,14 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingatlocations_gridexact) { model.init(); // define GCV function and grid of \lambda_D values auto GCV = model.gcv(); - std::vector> lambdas; + std::vector> lambdas; for (double x = -5.0; x <= -3.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test3/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test3/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test3/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test3/gcvs.mtx")); } // test 4 @@ -173,20 +170,20 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingatlocations_gridexact) { // BC: no // order FE: 1 // GCV optimization: grid stochastic -TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingatlocations_gridstochastic) { +TEST(gcv_qsrpde_test, laplacian_semiparametric_samplingatlocations_gridstochastic) { // define domain MeshLoader domain("c_shaped"); // import data from files - DMatrix locs = read_csv("../data/gcv/sqrpde/2D_test4/locs.csv"); - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test4/y.csv"); - DMatrix X = read_csv("../data/gcv/sqrpde/2D_test4/X.csv"); + DMatrix locs = read_csv("../data/gcv/qsrpde/2D_test4/locs.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test4/y.csv"); + DMatrix X = read_csv("../data/gcv/qsrpde/2D_test4/X.csv"); // 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 double alpha = 0.9; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::pointwise, alpha); model.set_spatial_locations(locs); // set model's data BlockFrame df; @@ -197,14 +194,14 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingatlocations_gridstochasti // define GCV function and grid of \lambda_D value std::size_t seed = 66546513; auto GCV = model.gcv(1000, seed); - std::vector> lambdas; + std::vector> lambdas; for (double x = -5.0; x <= -3.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test4/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test4/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test4/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test4/gcvs.mtx")); } // test 5 @@ -215,11 +212,11 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingatlocations_gridstochasti // BC: no // order FE: 1 // GCV optimization: grid exact -TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gridexact) { +TEST(gcv_qsrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gridexact) { // define domain MeshLoader domain("unit_square_coarse"); // import data from files - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test5/y.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test5/y.csv"); // define regularizing PDE SMatrix<2> K; K << 1, 0, 0, 4; @@ -228,7 +225,7 @@ TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gride PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); // define model double alpha = 0.1; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::mesh_nodes, alpha); // set model's data BlockFrame df; df.insert(OBSERVATIONS_BLK, y); @@ -236,14 +233,14 @@ TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gride model.init(); // define GCV function and grid of \lambda_D values auto GCV = model.gcv(); - std::vector> lambdas; + std::vector> lambdas; for (double x = -7.0; x <= -5.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // optimize gcv field // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test5/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test5/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test5/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test5/gcvs.mtx")); } // test 6 @@ -254,11 +251,11 @@ TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gride // BC: no // order FE: 1 // GCV optimization: grid stochastic -TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gridstochastic) { +TEST(gcv_qsrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_gridstochastic) { // define domain MeshLoader domain("unit_square_coarse"); // import data from files - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test6/y.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test6/y.csv"); // define regularizing PDE SMatrix<2> K; K << 1, 0, 0, 4; @@ -267,7 +264,7 @@ TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_grids PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); // define model double alpha = 0.1; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::mesh_nodes, alpha); // set model's data BlockFrame df; df.insert(OBSERVATIONS_BLK, y); @@ -276,14 +273,14 @@ TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_grids // define GCV function and grid of \lambda_D values std::size_t seed = 438172; auto GCV = model.gcv(1000, seed); - std::vector> lambdas; + std::vector> lambdas; for (double x = -7.0; x <= -5.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // optimize gcv field // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test6/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test6/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test6/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test6/gcvs.mtx")); } // test 7 @@ -294,20 +291,20 @@ TEST(gcv_sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes_grids // BC: yes // order FE: 1 // GCV optimization: grid exact -TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingareal_gridexact) { +TEST(gcv_qsrpde_test, laplacian_semiparametric_samplingareal_gridexact) { // define domain and regularizing PDE MeshLoader domain("c_shaped_areal"); // import data from files - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test7/y.csv"); - DMatrix X = read_csv("../data/gcv/sqrpde/2D_test7/X.csv"); - DMatrix subdomains = read_csv("../data/gcv/sqrpde/2D_test7/incidence_matrix.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test7/y.csv"); + DMatrix X = read_csv("../data/gcv/qsrpde/2D_test7/X.csv"); + DMatrix subdomains = read_csv("../data/gcv/qsrpde/2D_test7/incidence_matrix.csv"); // 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 double alpha = 0.5; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::areal, alpha); model.set_spatial_locations(subdomains); // set model's data BlockFrame df; @@ -317,14 +314,14 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingareal_gridexact) { model.init(); // define GCV function and grid of \lambda_D values auto GCV = model.gcv(); - std::vector> lambdas; + std::vector> lambdas; for (double x = -4.0; x <= -1.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // optimize gcv field // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test7/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test7/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test7/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test7/gcvs.mtx")); } // test 8 @@ -335,21 +332,21 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingareal_gridexact) { // BC: yes // order FE: 1 // GCV optimization: grid stochastic -TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingareal_gridstochastic) { +TEST(gcv_qsrpde_test, laplacian_semiparametric_samplingareal_gridstochastic) { // define domain and regularizing PDE // define domain and regularizing PDE MeshLoader domain("c_shaped_areal"); // import data from files - DMatrix y = read_csv("../data/gcv/sqrpde/2D_test8/y.csv"); - DMatrix X = read_csv("../data/gcv/sqrpde/2D_test8/X.csv"); - DMatrix subdomains = read_csv("../data/gcv/sqrpde/2D_test8/incidence_matrix.csv"); + DMatrix y = read_csv("../data/gcv/qsrpde/2D_test8/y.csv"); + DMatrix X = read_csv("../data/gcv/qsrpde/2D_test8/X.csv"); + DMatrix subdomains = read_csv("../data/gcv/qsrpde/2D_test8/incidence_matrix.csv"); // 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 double alpha = 0.5; - SQRPDE model(problem, alpha); + QSRPDE model(problem, Sampling::areal, alpha); model.set_spatial_locations(subdomains); // set model's data BlockFrame df; @@ -360,12 +357,12 @@ TEST(gcv_sqrpde_test, laplacian_semiparametric_samplingareal_gridstochastic) { // define GCV function and grid of \lambda_D values std::size_t seed = 438172; auto GCV = model.gcv(100, seed); - std::vector> lambdas; + std::vector> lambdas; for (double x = -4.0; x <= -1.0; x += 0.25) lambdas.push_back(SVector<1>(std::pow(10, x))); // optimize GCV - Grid<1> opt; + Grid opt; opt.optimize(GCV, lambdas); // test correctness - EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/sqrpde/2D_test8/edfs.mtx")); - EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/sqrpde/2D_test8/gcvs.mtx")); + EXPECT_TRUE(almost_equal(GCV.edfs(), "../data/gcv/qsrpde/2D_test8/edfs.mtx")); + EXPECT_TRUE(almost_equal(GCV.gcvs(), "../data/gcv/qsrpde/2D_test8/gcvs.mtx")); } diff --git a/test/src/sqrpde_test.cpp b/test/src/qsrpde_test.cpp similarity index 79% rename from test/src/sqrpde_test.cpp rename to test/src/qsrpde_test.cpp index a6a9cacf..17e49258 100644 --- a/test/src/sqrpde_test.cpp +++ b/test/src/qsrpde_test.cpp @@ -23,9 +23,9 @@ using fdapde::core::fem_order; using fdapde::core::laplacian; using fdapde::core::PDE; -#include "../../fdaPDE/models/regression/sqrpde.h" +#include "../../fdaPDE/models/regression/qsrpde.h" #include "../../fdaPDE/models/sampling_design.h" -using fdapde::models::SQRPDE; +using fdapde::models::QSRPDE; using fdapde::models::SpaceTimeSeparable; using fdapde::models::SpaceTimeParabolic; using fdapde::models::SpaceOnly; @@ -45,11 +45,11 @@ using fdapde::testing::read_csv; // covariates: no // BC: no // order FE: 1 -TEST(sqrpde_test, laplacian_nonparametric_samplingatnodes) { +TEST(qsrpde_test, laplacian_nonparametric_samplingatnodes) { // define domain MeshLoader domain("unit_square_coarse"); // import data from files - DMatrix y = read_csv("../data/models/sqrpde/2D_test1/y.csv"); + DMatrix y = read_csv("../data/models/qsrpde/2D_test1/y.csv"); // define regularizing PDE auto L = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); @@ -57,7 +57,7 @@ TEST(sqrpde_test, laplacian_nonparametric_samplingatnodes) { // define model double lambda = 1.778279 * std::pow(0.1, 4); double alpha = 0.1; - SQRPDE model(problem, Sampling::mesh_nodes, alpha); + QSRPDE model(problem, Sampling::mesh_nodes, alpha); model.set_lambda_D(lambda); // set model's data BlockFrame df; @@ -67,7 +67,7 @@ TEST(sqrpde_test, laplacian_nonparametric_samplingatnodes) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.f(), "../data/models/sqrpde/2D_test1/sol.mtx")); + EXPECT_TRUE(almost_equal(model.f(), "../data/models/qsrpde/2D_test1/sol.mtx")); } // test 2 @@ -77,13 +77,13 @@ TEST(sqrpde_test, laplacian_nonparametric_samplingatnodes) { // covariates: yes // BC: no // order FE: 1 -TEST(sqrpde_test, laplacian_semiparametric_samplingatlocations) { +TEST(qsrpde_test, laplacian_semiparametric_samplingatlocations) { // define domain and regularizing PDE MeshLoader domain("c_shaped"); // import data from files - DMatrix locs = read_csv("../data/models/sqrpde/2D_test2/locs.csv"); - DMatrix y = read_csv("../data/models/sqrpde/2D_test2/y.csv"); - DMatrix X = read_csv("../data/models/sqrpde/2D_test2/X.csv"); + DMatrix locs = read_csv("../data/models/qsrpde/2D_test2/locs.csv"); + DMatrix y = read_csv("../data/models/qsrpde/2D_test2/y.csv"); + DMatrix X = read_csv("../data/models/qsrpde/2D_test2/X.csv"); // define regularizing PDE auto L = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); @@ -91,7 +91,7 @@ TEST(sqrpde_test, laplacian_semiparametric_samplingatlocations) { // define statistical model double alpha = 0.9; double lambda = 3.162277660168379 * std::pow(0.1, 4); // use optimal lambda to avoid possible numerical issues - SQRPDE model(problem, Sampling::pointwise, alpha); + QSRPDE model(problem, Sampling::pointwise, alpha); model.set_lambda_D(lambda); model.set_spatial_locations(locs); // set model's data @@ -103,8 +103,8 @@ TEST(sqrpde_test, laplacian_semiparametric_samplingatlocations) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.f() , "../data/models/sqrpde/2D_test2/sol.mtx" )); - EXPECT_TRUE(almost_equal(model.beta(), "../data/models/sqrpde/2D_test2/beta.mtx")); + EXPECT_TRUE(almost_equal(model.f() , "../data/models/qsrpde/2D_test2/sol.mtx" )); + EXPECT_TRUE(almost_equal(model.beta(), "../data/models/qsrpde/2D_test2/beta.mtx")); } // test 3 @@ -114,11 +114,11 @@ TEST(sqrpde_test, laplacian_semiparametric_samplingatlocations) { // covariates: no // BC: no // order FE: 1 -TEST(sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes) { +TEST(qsrpde_test, costantcoefficientspde_nonparametric_samplingatnodes) { // define domain and regularizing PDE MeshLoader domain("unit_square_coarse"); // import data from files - DMatrix y = read_csv("../data/models/sqrpde/2D_test3/y.csv"); + DMatrix y = read_csv("../data/models/qsrpde/2D_test3/y.csv"); // define regularizing PDE SMatrix<2> K; K << 1, 0, 0, 4; @@ -128,7 +128,7 @@ TEST(sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes) { // define statistical model double alpha = 0.1; double lambda = 5.623413251903491 * pow(0.1, 4); - SQRPDE model(problem, Sampling::mesh_nodes, alpha); + QSRPDE model(problem, Sampling::mesh_nodes, alpha); model.set_lambda_D(lambda); // set model data BlockFrame df; @@ -138,7 +138,7 @@ TEST(sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.f(), "../data/models/sqrpde/2D_test3/sol.mtx")); + EXPECT_TRUE(almost_equal(model.f(), "../data/models/qsrpde/2D_test3/sol.mtx")); } // test 4 @@ -148,13 +148,13 @@ TEST(sqrpde_test, costantcoefficientspde_nonparametric_samplingatnodes) { // covariates: yes // BC: no // order FE: 1 -TEST(sqrpde_test, laplacian_semiparametric_samplingareal) { +TEST(qsrpde_test, laplacian_semiparametric_samplingareal) { // define domain and regularizing PDE MeshLoader domain("c_shaped_areal"); // import data from files - DMatrix y = read_csv("../data/models/sqrpde/2D_test4/y.csv"); - DMatrix X = read_csv("../data/models/sqrpde/2D_test4/X.csv"); - DMatrix subdomains = read_csv("../data/models/sqrpde/2D_test4/incidence_matrix.csv"); + DMatrix y = read_csv("../data/models/qsrpde/2D_test4/y.csv"); + DMatrix X = read_csv("../data/models/qsrpde/2D_test4/X.csv"); + DMatrix subdomains = read_csv("../data/models/qsrpde/2D_test4/incidence_matrix.csv"); // define regularizing PDE auto L = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); @@ -162,7 +162,7 @@ TEST(sqrpde_test, laplacian_semiparametric_samplingareal) { // define statistical model double alpha = 0.5; double lambda = 5.623413251903491 * std::pow(0.1, 3); // use optimal lambda to avoid possible numerical issues - SQRPDE model(problem, Sampling::areal, alpha); + QSRPDE model(problem, Sampling::areal, alpha); model.set_lambda_D(lambda); model.set_spatial_locations(subdomains); // set model data @@ -174,6 +174,6 @@ TEST(sqrpde_test, laplacian_semiparametric_samplingareal) { model.init(); model.solve(); // test correctness - EXPECT_TRUE(almost_equal(model.f() , "../data/models/sqrpde/2D_test4/sol.mtx" )); - EXPECT_TRUE(almost_equal(model.beta(), "../data/models/sqrpde/2D_test4/beta.mtx")); + EXPECT_TRUE(almost_equal(model.f() , "../data/models/qsrpde/2D_test4/sol.mtx" )); + EXPECT_TRUE(almost_equal(model.beta(), "../data/models/qsrpde/2D_test4/beta.mtx")); }