use parabolic penalty to derive time domain and initial condition, fi…

…xed all models and bases to this change, fixed tests, fPCA returns fitted_loadings, other minor fixes

fdaPDE/models/functional/fpca.h

@@ -46,8 +46,10 @@ class FPCA<RegularizationType_, sequential> :
     using This = FPCA<RegularizationType, sequential>;
     using Base = FunctionalBase<This, RegularizationType>;
     IMPORT_MODEL_SYMBOLS;
-    using Base::lambda;
-    using Base::X;
+    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;
@@ -59,8 +61,10 @@ class FPCA<RegularizationType_, sequential> :
 
     void init_model() { return; };
     void solve() {
-        loadings_.resize(X().cols(), n_pc_);
-        scores_.resize(X().rows(), n_pc_);
+        // preallocate space
+        loadings_.resize(n_basis(), n_pc_);
+        fitted_loadings_.resize(n_locs(), n_pc_);
+        scores_.resize(n_stat_units(), n_pc_);
         DMatrix<double> X_ = X();   // copy original data to avoid side effects
 
         // first guess of PCs set to a multivariate PCA (SVD)
@@ -111,14 +115,16 @@ class FPCA<RegularizationType_, sequential> :
             } break;
             }
             // store results
-            loadings_.col(i) = Psi() * solver.f(); // \frac{\Psi * f}{\norm{f}_{L^2}}
+            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) * loadings_.col(i).transpose();   // deflation step
+            X_ -= scores_.col(i) * fitted_loadings_.col(i).transpose();   // X <- X - s*f_n^\top (deflation step)
         }
         return;
     }
 
     // getters
+    const DMatrix<double>& fitted_loadings() const { return fitted_loadings_; }
     const DMatrix<double>& loadings() const { return loadings_; }
     const DMatrix<double>& scores() const { return scores_; }
     // setters
@@ -145,7 +151,8 @@ class FPCA<RegularizationType_, sequential> :
     int seed_ = fdapde::random_seed;
 
     // problem solution
-    DMatrix<double> loadings_;   // evaluation of the PC functions at data locations
+    DMatrix<double> loadings_;          // PC functions' expansion coefficients
+    DMatrix<double> fitted_loadings_;   // evaluation of the PC functions at data locations
     DMatrix<double> scores_;
 };
 
@@ -173,19 +180,22 @@ class FPCA<RegularizationType_, monolithic> :
         RSVD<This> solver(this);
         solver.compute(X(), lambda(), n_pc_);
         // store results
-        loadings_ = Psi(not_nan()) * solver.W();   // evaluation of L^2 normalized PC functions at data locations
+	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;
     }
 
     // getters
+    const DMatrix<double>& fitted_loadings() const { return fitted_loadings_; }
     const DMatrix<double>& loadings() const { return loadings_; }
     const DMatrix<double>& scores() const { return scores_; }
     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<double> loadings_;   // evaluation of the PC functions at data locations
+    DMatrix<double> loadings_;          // PC functions' expansion coefficients
+    DMatrix<double> fitted_loadings_;   // evaluation of the PC functions at data locations
     DMatrix<double> scores_;
 };
 

fdaPDE/models/functional/functional_base.h

@@ -41,18 +41,15 @@ class FunctionalBase : public select_regularization_base<Model, RegularizationTy
     using SamplingBase<Model>::PsiTD;   // block \Psi^T*D
 
     FunctionalBase() = default;
-    // space-only constructor
-    fdapde_enable_constructor_if(is_space_only, Model) FunctionalBase(const pde_ptr& pde, Sampling s) :
+      // space-only and space-time parabolic constructor (they require only one PDE)
+    fdapde_enable_constructor_if(has_single_penalty, Model) FunctionalBase(const pde_ptr& pde, Sampling s) :
         Base(pde), SamplingBase<Model>(s) {};
-    // space-time constructors
-    fdapde_enable_constructor_if(is_space_time_separable, Model)
+    // space-time separable constructor
+    fdapde_enable_constructor_if(has_double_penalty, Model)
       FunctionalBase(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s) :
         Base(space_penalty, time_penalty), SamplingBase<Model>(s) {};
-    fdapde_enable_constructor_if(is_space_time_parabolic, Model)
-      FunctionalBase(const pde_ptr& pde, Sampling s, const DVector<double>& time) :
-        Base(pde, time), SamplingBase<Model>(s) {};
 
-    // getters
+  // getters
     const DMatrix<double>& X() const { return df_.template get<double>(OBSERVATIONS_BLK); }   // observation matrix X
     std::size_t n_stat_units() const { return X().rows(); }
     std::size_t n_obs() const { return X().size(); };

fdaPDE/models/functional/power_iteration.h

@@ -61,18 +61,11 @@ template <typename Model_> class PowerIteration {
 
     // initializes internal smoothing solver
     void init() {
-        // define internal problem solver
-        if constexpr (!is_space_time<Model>::value)   // space-only
+        if constexpr (is_space_only<Model>::value) {   // space-only solver
             solver_ = SolverType(m_->pde(), m_->sampling());
-        else {   // space-time
-            if constexpr (is_space_time_parabolic<Model>::value) {
-                solver_ = SolverType(m_->pde(), m_->sampling(), m_->time_domain());
-                solver_.set_initial_condition(m_->s(), false);
-            }
-            if constexpr (is_space_time_separable<Model>::value) {
-                solver_ = SolverType(m_->pde(), m_->time_pde(), m_->sampling());
-                solver_.set_temporal_locations(m_->time_locs());
-            }
+        } else {   // space-time separable solver
+            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

fdaPDE/models/model_traits.h

@@ -41,24 +41,29 @@ template <typename T> struct is_stat_model {
     static constexpr bool value = fdapde::is_base_of_template<ModelBase, T>::value;
 };
 
+// space-only regularization trait
 template <typename Model> struct is_space_only {
-    using ModelType = typename std::decay<Model>::type;
-    static constexpr bool value =
-      std::is_same<typename ModelType::RegularizationType, SpaceOnly>::value;
+    static constexpr bool value = std::is_same<typename std::decay<Model>::type::RegularizationType, SpaceOnly>::value;
 };
+// traits for space-time regularizations
 template <typename Model> struct is_space_time {
     static constexpr bool value = !is_space_only<Model>::value;
 };
-// specific traits for space-time regularizations
-template <typename Model> struct is_space_time_separable {   // separable regularization
-    using ModelType = typename std::decay<Model>::type;
+template <typename Model> struct is_space_time_separable {
     static constexpr bool value =
-      std::is_same<typename ModelType::RegularizationType, SpaceTimeSeparable>::value;
+      std::is_same<typename std::decay<Model>::type::RegularizationType, SpaceTimeSeparable>::value;
 };
-template <typename Model> struct is_space_time_parabolic {   // parabolic regularization
-    using ModelType = typename std::decay<Model>::type;
+template <typename Model> struct is_space_time_parabolic {
     static constexpr bool value =
-      std::is_same<typename ModelType::RegularizationType, SpaceTimeParabolic>::value;
+      std::is_same<typename std::decay<Model>::type::RegularizationType, SpaceTimeParabolic>::value;
+};
+
+template <typename Model_> struct has_single_penalty {
+    using Model = typename std::decay<Model_>::type;
+    static constexpr bool value = (is_space_only<Model>::value || is_space_time_parabolic<Model>::value);
+};
+template <typename Model_> struct has_double_penalty {
+    static constexpr bool value = is_space_time_separable<typename std::decay<Model_>::type>::value;
 };
 
 // selects the regularization base depending on the Model regularization type

fdaPDE/models/model_type_erasure.h

@@ -37,57 +37,54 @@ template <typename RegularizationType> struct IStatModel { };
     }                                                                                                                  \
     void set_spatial_locations(const DMatrix<double>& locs) { fdapde::invoke<void, 4>(*this, locs); }                  \
     BlockFrame<double, int>& data() { return fdapde::invoke<BlockFrame<double, int>&, 5>(*this); }                     \
-    decltype(auto) R0() const { return fdapde::invoke<const SpMatrix<double>&, 6>(*this); }
+    decltype(auto) R0() const { return fdapde::invoke<const SpMatrix<double>&, 6>(*this); }                            \
+    std::size_t n_locs() const { return fdapde::invoke<std::size_t, 7>(*this); }                                       \
+    std::size_t n_basis() const { return fdapde::invoke<std::size_t, 8>(*this); }
 
 // space-only statistical model interface
 template <> struct IStatModel<SpaceOnly> {
+    using RegularizationType = SpaceOnly;
     template <typename M>
     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<BlockFrame<double, int>& (ModelBase<M>::*)()>(&M::data), &M::R0,
+      static_cast<BlockFrame<double, int>& (ModelBase<M>::*)()>(&M::data), &M::R0, &M::n_locs, &M::n_basis,
       static_cast<void (SpaceOnlyBase<M>::*)(const SVector<1>&)>(&M::set_lambda)>;
     // interface implementation
     DEFINE_BASE_STAT_MODEL_INTERFACE;
-    void set_lambda(const SVector<1>& lambda) { fdapde::invoke<void, 7>(*this, lambda); }
+    void set_lambda(const SVector<1>& lambda) { fdapde::invoke<void, 9>(*this, lambda); }
 };
 
 // space-time separable statistical model interface
 template <> struct IStatModel<SpaceTimeSeparable> {
+    using RegularizationType = SpaceTimeSeparable;
     template <typename M>
     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<BlockFrame<double, int>& (ModelBase<M>::*)()>(&M::data), &M::R0,
+      static_cast<BlockFrame<double, int>& (ModelBase<M>::*)()>(&M::data), &M::R0, &M::n_locs, &M::n_basis,
       static_cast<void (SpaceTimeBase<M, SpaceTimeSeparable>::*)(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<void, 7>(*this, lambda); }
-    void set_lambda_T(double lambda_T) { fdapde::invoke<void, 8>(*this, lambda_T); }
-    void set_temporal_locations(const DMatrix<double>& locs) { fdapde::invoke<void, 9>(*this, locs); }
-
-    using RegularizationType = SpaceTimeSeparable;
+    void set_lambda(const SVector<2>& lambda) { fdapde::invoke<void, 9>(*this, lambda); }
+    void set_lambda_T(double lambda_T) { fdapde::invoke<void, 10>(*this, lambda_T); }
+    void set_temporal_locations(const DMatrix<double>& locs) { fdapde::invoke<void, 11>(*this, locs); }
 };
 
 // space-time parabolic statistical model interface
 template <> struct IStatModel<SpaceTimeParabolic> {
+    using RegularizationType = SpaceTimeParabolic;
     template <typename M>
     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<BlockFrame<double, int>& (ModelBase<M>::*)()>(&M::data), &M::R0,
-      static_cast<void (SpaceTimeBase<M, SpaceTimeParabolic>::*)(const SVector<2>&)>(&M::set_lambda), &M::set_lambda_T,
-      &M::set_initial_condition>;
+      static_cast<BlockFrame<double, int>& (ModelBase<M>::*)()>(&M::data), &M::R0, &M::n_locs, &M::n_basis,
+      static_cast<void (SpaceTimeBase<M, SpaceTimeParabolic>::*)(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<void, 7>(*this, lambda); }
-    void set_lambda_T(double lambda_T) { fdapde::invoke<void, 8>(*this, lambda_T); }
-    void set_initial_condition(const DMatrix<double>& s, bool shift = true) {
-        fdapde::invoke<void, 9>(*this, s, shift);
-    }
-
-    using RegularizationType = SpaceTimeParabolic;
+    void set_lambda(const SVector<2>& lambda) { fdapde::invoke<void, 9>(*this, lambda); }
+    void set_lambda_T(double lambda_T) { fdapde::invoke<void, 10>(*this, lambda_T); }
 };
 
 // regularization-independent model interface
@@ -97,7 +94,8 @@ template <> struct IStatModel<void> {
       &M::init, &M::solve, static_cast<void (ModelBase<M>::*)(const DVector<double>&)>(&M::set_lambda),
       static_cast<DVector<double> (ModelBase<M>::*)(int) const>(&M::lambda),
       static_cast<const SpMatrix<double>& (SamplingBase<M>::*)(not_nan) const>(&M::Psi),
-      static_cast<const SpMatrix<double>& (SamplingBase<M>::*)(not_nan) const>(&M::PsiTD), &M::set_data>;
+      static_cast<const SpMatrix<double>& (SamplingBase<M>::*)(not_nan) const>(&M::PsiTD), &M::set_data,
+      &M::n_locs, &M::n_basis>;
     // interface implementation
     void init()  { fdapde::invoke<void, 0>(*this); }
     void solve() { fdapde::invoke<void, 1>(*this); }
@@ -108,6 +106,8 @@ template <> struct IStatModel<void> {
     void set_data(const BlockFrame<double, int>& data, bool reindex = false) {
         fdapde::invoke<void, 6>(*this, data, reindex);
     }
+    std::size_t n_locs() const { return fdapde::invoke<std::size_t, 7>(*this); }
+    std::size_t n_basis() const { return fdapde::invoke<std::size_t, 8>(*this); }
 };
 
 }   // namespace models

fdaPDE/models/regression/fpirls.h

@@ -54,18 +54,11 @@ template <typename Model_> class FPIRLS {
             using SolverType = typename std::conditional<
               is_space_only<Model>::value, SRPDE,
               STRPDE<typename Model::RegularizationType, fdapde::monolithic> >::type;
-            // define internal problem solver
-            if constexpr (!is_space_time<Model_>::value)   // space-only
+            if constexpr (is_space_only<Model_>::value || is_space_time_parabolic<Model_>::value) {
                 solver_ = SolverType(m_->pde(), m_->sampling());
-            else {   // space-time
-                if constexpr (is_space_time_parabolic<Model_>::value) {
-                    solver_ = SolverType(m_->pde(), m_->sampling(), m_->time_domain());
-                    solver_.set_initial_condition(m_->s(), false);
-                }
-                if constexpr (is_space_time_separable<Model_>::value) {
-                    solver_ = SolverType(m_->pde(), m_->time_pde(), m_->sampling());
-                    solver_.set_temporal_locations(m_->time_locs());
-                }
+            } else {   // space-time separable
+                solver_ = SolverType(m_->pde(), m_->time_pde(), m_->sampling());
+                solver_.set_temporal_locations(m_->time_locs());
             }
             // solver initialization
             solver_.set_spatial_locations(m_->locs());

fdaPDE/models/regression/gsrpde.h

@@ -43,21 +43,16 @@ class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, Regularization
     using Base::XtWX_;      // q x q matrix X^T*W*X
     // constructor
     GSRPDE() = default;
-    // space-only constructor
-    fdapde_enable_constructor_if(is_space_only, This)
+    // space-only and space-time parabolic constructor
+    fdapde_enable_constructor_if(has_single_penalty, This)
       GSRPDE(const pde_ptr& pde, Sampling s, const Distribution& distr) :
         Base(pde, s), distr_(distr) {
         fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
     };
-    // space-time constructors
-    fdapde_enable_constructor_if(is_space_time_separable, This)
-    GSRPDE(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s, const Distribution& distr) :
-      Base(space_penalty, time_penalty, s), distr_(distr) {
-        fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
-    };
-    fdapde_enable_constructor_if(is_space_time_parabolic, This)
-    GSRPDE(const pde_ptr& space_penalty, const DVector<double>& time, Sampling s, const Distribution& distr) :
-      Base(space_penalty, s, time), distr_(distr) {
+    // space-time separable constructor
+    fdapde_enable_constructor_if(has_double_penalty, This)
+      GSRPDE(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s, const Distribution& distr) :
+        Base(space_penalty, time_penalty, s), distr_(distr) {
         fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
     };
 

fdaPDE/models/regression/qsrpde.h

@@ -43,19 +43,15 @@ class QSRPDE : public RegressionBase<QSRPDE<RegularizationType_>, Regularization
     using Base::XtWX_;      // q x q matrix X^T*W*X
     // constructor
     QSRPDE() = default;
-    // space-only constructor
-    template <
-      typename U = RegularizationType,
-      typename std::enable_if<std::is_same<U, SpaceOnly>::value, int>::type = 0>
+    // space-only and space-time parabolic constructor
+    fdapde_enable_constructor_if(has_single_penalty, This)
     QSRPDE(const pde_ptr& pde, Sampling s, double alpha = 0.5) : Base(pde, s), alpha_(alpha) {
         fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
     };
-    // space-time constructor
-    template <
-      typename U = RegularizationType,
-      typename std::enable_if<!std::is_same<U, SpaceOnly>::value, int>::type = 0>
-    QSRPDE(const pde_ptr& pde, const DVector<double>& time, Sampling s, double alpha = 0.5) :
-        Base(pde, s, time), alpha_(alpha) {
+    // space-time separable constructor
+    fdapde_enable_constructor_if(has_double_penalty, This)
+    QSRPDE(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s, double alpha = 0.5) :
+        Base(space_penalty, time_penalty, s), alpha_(alpha) {
         fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
     };
 

fdaPDE/models/regression/regression_base.h

@@ -70,16 +70,13 @@ class RegressionBase :
     using Base::model;
 
     RegressionBase() = default;
-    // space-only constructor
-    fdapde_enable_constructor_if(is_space_only, Model) RegressionBase(const pde_ptr& pde, Sampling s) :
+    // space-only and space-time parabolic constructor (they require only one PDE)
+    fdapde_enable_constructor_if(has_single_penalty, Model) RegressionBase(const pde_ptr& pde, Sampling s) :
         Base(pde), SamplingBase<Model>(s) {};
-    // space-time constructors
-    fdapde_enable_constructor_if(is_space_time_separable, Model)
+    // space-time separable constructor
+    fdapde_enable_constructor_if(has_double_penalty, Model)
       RegressionBase(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s) :
         Base(space_penalty, time_penalty), SamplingBase<Model>(s) {};
-    fdapde_enable_constructor_if(is_space_time_parabolic, Model)
-      RegressionBase(const pde_ptr& pde, Sampling s, const DVector<double>& time) :
-        Base(pde, time), SamplingBase<Model>(s) {};
 
     // getters
     const DMatrix<double>& y() const { return df_.template get<double>(OBSERVATIONS_BLK); }   // observation vector y

fdaPDE/models/regression/stochastic_edf.h

@@ -51,7 +51,7 @@ class StochasticEDF {
             // compute sample from Rademacher distribution
             std::mt19937 rng(seed_);
             std::bernoulli_distribution Be(0.5);   // bernulli distribution with parameter p = 0.5
-            Us_.resize(model_.n_locs(), r_);        // preallocate memory for matrix Us
+            Us_.resize(model_.n_locs(), r_);       // preallocate memory for matrix Us
             // fill matrix
             for (std::size_t i = 0; i < model_.n_locs(); ++i) {
                 for (std::size_t j = 0; j < r_; ++j) {