unstable support for QSRPDE (quantile regression) for space-time data…

…, use smoothed pinball loss function for norm computation in GCV, fix space_time_separable penalty matrices (inverted kronecker tensor products), KFold CV support (unstable) for FPIRLS-based models (GSR/QSR/etc.)

.gitignore

@@ -1,193 +1,12 @@
-#ignore all
-*
-
-#but not...
-!*.h
-!*.cpp
-!*.tpp
-!*.hpp
-!*.md
-!*.yml
-!fdaPDE/
-|fdaPDE/*
-!fdaPDE/models/
-!fdaPDE/models/*
-!fdaPDE/models/regression/
-!fdaPDE/models/regression/*
-!fdaPDE/models/functional/
-!fdaPDE/models/functional/*
-!fdaPDE/calibration/
-!fdaPDE/calibration/*
-!test/
-!test/*
-!test/src/
-!test/src/*
-!test/src/utils/
-!test/src/utils/*
-!test/data/
-!test/data/*
-!test/data/mesh/
-!test/data/mesh/*
-!test/data/mesh/quasi_circle/
-!test/data/mesh/quasi_circle/*
-!test/data/mesh/unit_sphere/
-!test/data/mesh/unit_sphere/*
-!test/data/mesh/surface/
-!test/data/mesh/surface/*
-!test/data/mesh/network/
-!test/data/mesh/network/*
-!test/data/mesh/c_shaped/
-!test/data/mesh/c_shaped/*
-!test/data/mesh/c_shaped_areal/
-!test/data/mesh/c_shaped_areal/*
-!test/data/models/
-!test/data/models/*
-!test/data/models/srpde/
-!test/data/models/srpde/*
-!test/data/models/strpde/
-!test/data/models/strpde/*
-!test/data/models/gsrpde/
-!test/data/models/gsrpde/*
-!test/data/models/gcv/
-!test/data/models/gcv/*
-!test/data/gcv/qsrpde/
-!test/data/gcv/qsrpde/*
-!test/data/models/fpca/
-!test/data/models/fpca/*
-!test/data/models/centering/
-!test/data/models/centering/*
-!test/data/models/fpls/
-!test/data/models/fpls/*
-!.github/
-!.github/*
-!.github/workflows/
-!.github/workflows/*
-!.gitattributes
-!.clang-format
-
-*.csv
-!test/data/mesh/quasi_circle/*.csv
-!test/data/mesh/unit_sphere/*.csv
-!test/data/mesh/surface/*.csv
-!test/data/mesh/network/*.csv
-!test/data/mesh/c_shaped/*.csv
-!test/data/mesh/unit_square/*.csv
-!test/data/mesh/unit_interval/*.csv
-!test/data/mesh/unit_interval_o2/*.csv
-!test/data/mesh/unit_square_coarse/*.csv
-!test/data/mesh/unit_square_medium/*.csv
-!test/data/mesh/c_shaped_areal/*.csv
-!test/data/*.mtx
-!test/data/models/srpde/
-!test/data/models/srpde/*
-!test/data/models/srpde/2D_test1/
-!test/data/models/srpde/2D_test1/*
-!test/data/models/srpde/2D_test2/
-!test/data/models/srpde/2D_test2/*
-!test/data/models/srpde/2D_test3/
-!test/data/models/srpde/2D_test3/*
-!test/data/models/srpde/2D_test4/
-!test/data/models/srpde/2D_test4/*
-!test/data/models/strpde/
-!test/data/models/strpde/*
-!test/data/models/strpde/2D_test1/
-!test/data/models/strpde/2D_test1/*
-!test/data/models/strpde/2D_test2/
-!test/data/models/strpde/2D_test2/*
-!test/data/models/strpde/2D_test3/
-!test/data/models/strpde/2D_test3/*
-!test/data/models/strpde/2D_test4/
-!test/data/models/strpde/2D_test4/*
-!test/data/models/strpde/2D_test5/
-!test/data/models/strpde/2D_test5/*
-!test/data/models/strpde/2D_test6/
-!test/data/models/strpde/2D_test6/*
-!test/data/models/gsrpde/
-!test/data/models/gsrpde/*
-!test/data/models/gsrpde/2D_test1/
-!test/data/models/gsrpde/2D_test1/*
-!test/data/models/gsrpde/2D_test2/
-!test/data/models/gsrpde/2D_test2/*
-!test/data/models/gsrpde/2D_test3/
-!test/data/models/gsrpde/2D_test3/*
-!test/data/models/gsrpde/2D_test4/
-!test/data/models/gsrpde/2D_test4/*
-!test/data/models/gsrpde/2D_test5/
-!test/data/models/gsrpde/2D_test5/*
-!test/data/models/gsrpde/2D_test6/
-!test/data/models/gsrpde/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/
-!test/data/models/gcv/2D_test2/*
-!test/data/models/gcv/2D_test3/
-!test/data/models/gcv/2D_test3/*
-!test/data/models/gcv/2D_test4/
-!test/data/models/gcv/2D_test4/*
-!test/data/models/gcv/2D_test5/
-!test/data/models/gcv/2D_test5/*
-!test/data/models/gcv/2D_test6/
-!test/data/models/gcv/2D_test6/*
-!test/data/models/gcv/2D_test7/
-!test/data/models/gcv/2D_test7/*
-!test/data/models/gcv/2D_test8/
-!test/data/models/gcv/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/
-!test/data/models/fpca/2D_test2/*
-!test/data/models/fpca/2D_test3/
-!test/data/models/fpca/2D_test3/*
-!test/data/models/fpca/2D_test4/
-!test/data/models/fpca/2D_test4/*
-!test/data/models/centering/2D_test1/
-!test/data/models/centering/2D_test1/*
-!test/data/models/fpls/2D_test1/
-!test/data/models/fpls/2D_test1/*
-!test/data/models/fpls/2D_test2/
-!test/data/models/fpls/2D_test2/*
 *.txt
 *.json
 *.cmake
 *Makefile
 *.Rhistory
 *.tar.gz
-*Main.cpp
-*Main
 *.RData
 *.o
 *.so
 test/fdapde_test
-
-!test/CMakeLists.txt
-tmp/
+test/.cache
+test/build/

fdaPDE/calibration/rmse.h

@@ -42,7 +42,7 @@ class RMSE {
       const BinaryVector<fdapde::Dynamic>& test_mask) {
         model_.set_lambda(lambda);
         // fit model on train set
-        model_.mask_obs(test_mask);   // discard test set from training phase
+        model_.set_mask(test_mask);   // discard test set from training phase
         model_.init();
         model_.solve();
 

fdaPDE/models/model_type_erasure.h

@@ -20,6 +20,7 @@
 #include <fdaPDE/utils.h>
 using fdapde::core::BlockFrame;
 #include "model_traits.h"
+#include "sampling_design.h"
 
 namespace fdapde {
 namespace models {
@@ -68,7 +69,6 @@ template <> struct StatisticalModel__<SpaceOnly> {
     // compile time constants
     using RegularizationType = SpaceOnly;
     static constexpr int n_lambda = 1;
-
     // interface implementation
     template <typename M>
     using fn_ptrs = fdapde::mem_fn_ptrs<
@@ -86,7 +86,6 @@ template <> struct StatisticalModel__<SpaceTimeSeparable> {
     // compile time constants
     using RegularizationType = SpaceTimeSeparable;
     static constexpr int n_lambda = 2;
-
     // interface implementation
     template <typename M>
     using fn_ptrs = fdapde::mem_fn_ptrs<
@@ -108,7 +107,6 @@ template <> struct StatisticalModel__<SpaceTimeParabolic> {
     // compile time constants
     using RegularizationType = SpaceTimeParabolic;
     static constexpr int n_lambda = 2;
-
     // interface implementation
     template <typename M>
     using fn_ptrs = fdapde::mem_fn_ptrs<

fdaPDE/models/regression/fpirls.h

@@ -36,7 +36,8 @@ namespace models {
 template <typename Model_> class FPIRLS {
    private:
     using Model = typename std::decay<Model_>::type;
-    using SolverWrapper = RegressionModel<typename std::decay_t<Model_>::RegularizationType>;
+    using RegularizationType = typename std::decay_t<Model_>::RegularizationType;
+    using SolverWrapper = RegressionModel<RegularizationType>;
     Model* m_;
     // algorithm's parameters
     double tolerance_;       // treshold on objective functional J to convergence
@@ -47,7 +48,7 @@ template <typename Model_> class FPIRLS {
     // constructor
     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
@@ -62,9 +63,10 @@ template <typename Model_> class FPIRLS {
             }
             // solver initialization
             solver_.set_spatial_locations(m_->locs());
-            solver_.set_data(m_->data());   // possible covariates are passed from here
+            solver_.set_data(m_->data());
         }
-	solver_.set_lambda(m_->lambda());
+        solver_.set_lambda(m_->lambda());     // derive smoothing level
+        solver_.set_mask(m_->masked_obs());   // derive missing and masking data pattern
     }
     // executes the FPIRLS algorithm
     void compute() {
@@ -83,10 +85,9 @@ template <typename Model_> class FPIRLS {
 	    solver_.init();
             solver_.solve();
             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());
-            // prepare for next iteration
-            k_++; J_old = J_new; J_new = J;
+            // update objective functional J = data_loss + f^\top * P_{\lambda}(f) * f 
+            k_++; J_old = J_new;
+	    J_new = m_->data_loss() + m_->ftPf(m_->lambda(), solver_.f(), solver_.g());;
         }
         return;
     }

fdaPDE/models/regression/gsrpde.h

@@ -59,10 +59,9 @@ class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, Regularization
     // 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() { fpirls_.init(); };
-    void solve() {   // finds a solution to the smoothing problem
+    void solve() {
+        fdapde_assert(y().rows() != 0);
         // 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
@@ -103,23 +102,21 @@ class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, Regularization
     const DVector<double>& py() const { return py_; }
     const DVector<double>& pW() const { return pW_; }
     const fdapde::SparseLU<SpMatrix<double>>& invA() const { return invA_; }
-
     // GCV support
     double norm(const DMatrix<double>& op1, const DMatrix<double>& op2) const {   // total deviance \sum dev(\hat y - y)
         DMatrix<double> mu = distr_.inv_link(op1);
         double result = 0;
-        for (std::size_t i = 0; i < op2.rows(); ++i) { result += distr_.deviance(mu.coeff(i, 0), op2.coeff(i, 0)); }
+        for (std::size_t i = 0; i < n_locs(); ++i) {
+            if (!Base::masked_obs()[i]) result += distr_.deviance(mu.coeff(i, 0), op2.coeff(i, 0));
+        }
         return result;
     }
-
-    virtual ~GSRPDE() = default;
    private:
     Distribution distr_ {};
     DVector<double> py_;   // \tilde y^k = G^k(y-u^k) + \theta^k
     DVector<double> pW_;   // diagonal of W^k = ((G^k)^{-2})*((V^k)^{-1})
     DVector<double> mu_;   // \mu^k = [ \mu^k_1, ..., \mu^k_n ] : mean vector at step k
-    DMatrix<double> T_;    // T = \Psi^T*Q*\Psi + P
-    fdapde::SparseLU<SpMatrix<double>> invA_;   // factorization of non-parametric system matrix A
+    fdapde::SparseLU<SpMatrix<double>> invA_;
 
     // FPIRLS parameters (set to default)
     FPIRLS<This> fpirls_;