better penalties management across space-only and space-time models, …

…fixed some minor issues

fdaPDE/models/model_base.h

@@ -17,11 +17,8 @@
 #ifndef __MODEL_BASE_H__
 #define __MODEL_BASE_H__
 
-#include <fdaPDE/mesh.h>
 #include <fdaPDE/utils.h>
-#include <fdaPDE/pde.h>
 using fdapde::core::BlockFrame;
-using fdapde::core::pde_ptr;
 
 #include "model_macros.h"
 #include "model_traits.h"
@@ -33,22 +30,18 @@ namespace models {
 // abstract base interface for any fdaPDE statistical model.
 template <typename Model> class ModelBase {
    public:
-    // constructors
     ModelBase() = default;
-    ModelBase(const pde_ptr& pde) : pde_(pde) { pde_.init(); };
     // full model stack initialization
     void init() {
-        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_data();    // specific data-dependent initialization requested by the model
         model().init_model();   // model initialization
     }
-
+  
     // setters
     void set_data(const BlockFrame<double, int>& df, bool reindex = false) {
         df_ = df;
@@ -60,20 +53,14 @@ template <typename Model> class ModelBase {
             df_.insert(INDEXES_BLK, idx);
         }
     }
-    void set_pde(const pde_ptr& pde) {
-        pde_ = pde;
-        model().runtime().set(runtime_status::require_pde_init);
-    }
     void set_lambda(const DVector<double>& lambda) {   // dynamic sized version of set_lambda provided by upper layers
 	model().set_lambda_D(lambda[0]);
 	if constexpr(is_space_time<Model>::value) model().set_lambda_T(lambda[1]);
     }
-
     // getters
     const BlockFrame<double, int>& data() const { return df_; }
     BlockFrame<double, int>& data() { return df_; }   // direct write-access to model's internal data storage
     const DMatrix<int>& idx() const { return df_.get<int>(INDEXES_BLK); }   // data indices
-    const pde_ptr& pde() const { return pde_; }   // regularizing term Lf - u (defined on some domain D)
     std::size_t n_locs() const { return model().n_spatial_locs() * model().n_temporal_locs(); }
     DVector<double> lambda(int) const { return model().lambda(); }   // int supposed to be fdapde::Dynamic
     // access to model runtime status
@@ -82,7 +69,6 @@ template <typename Model> class ModelBase {
 
     virtual ~ModelBase() = default;
    protected:
-    pde_ptr pde_ {};                     // regularizing term Lf - u
     BlockFrame<double, int> df_ {};      // blockframe for data storage
     model_runtime_handler runtime_ {};   // model's runtime status
 

fdaPDE/models/model_macros.h

@@ -28,7 +28,7 @@
                                  /* space-time problems) */                                                            \
     using Base::R0;              /* mass matrix in space (tensorized for space-time problems) */                       \
     using Base::u;               /* discretization of forcing term */                                                  \
-    using Base::pde;             /* differential operator L (regularizing term) */                                     \
+    using Base::pde;             /* differential operator Lf- u in space */                                            \
     using Base::data;            /* BlockFrame object containing data */
 
 // imports all basic symbols a model can expect from a RegressionBase

fdaPDE/models/regression/fpirls.h

@@ -63,7 +63,7 @@ template <typename Model_> class FPIRLS {
                     solver_.set_initial_condition(m_->s(), false);
                 }
                 if constexpr (is_space_time_separable<Model_>::value) {
-                    solver_ = SolverType(m_->pde(), m_->time_penalty(), m_->sampling());
+                    solver_ = SolverType(m_->pde(), m_->time_pde(), m_->sampling());
                     solver_.set_temporal_locations(m_->time_locs());
                 }
             }

fdaPDE/models/regression/regression_base.h

@@ -57,7 +57,6 @@ class RegressionBase :
     using Base::df_;           // BlockFrame for problem's data storage
     using Base::idx;           // indices of observations
     using Base::n_basis;       // number of basis function over domain D
-    using Base::pde_;          // differential operator L
     using Base::P;             // discretized penalty matrix
     using SamplingBase<Model>::D;     // for areal sampling, matrix of subdomains measures, identity matrix otherwise
     using SamplingBase<Model>::Psi;   // matrix of spatial basis evaluation at locations p_1 ... p_n
@@ -71,12 +70,10 @@ class RegressionBase :
     // space-time constructors
     fdapde_enable_constructor_if(is_space_time_separable, Model)
       RegressionBase(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s) :
-    Base(space_penalty, time_penalty), SamplingBase<Model>(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) {};
-    // copy constructor, copy only pde object (as a consequence also the problem domain)
-    RegressionBase(const RegressionBase& rhs) : Base(rhs), SamplingBase<Model>(rhs) { pde_ = rhs.pde_; }
 
     // getters
     const DMatrix<double>& y() const { return df_.template get<double>(OBSERVATIONS_BLK); }   // observation vector y

fdaPDE/models/regression/strpde.h

@@ -256,13 +256,15 @@ class STRPDE<SpaceTimeParabolic, iterative> :
     using Base::lambda_D;          // smoothing parameter in space
     using Base::lambda_T;          // smoothing parameter in time
     using Base::n_temporal_locs;   // number of time instants m defined over [0,T]
+    using Base::pde_;              // parabolic differential operator df/dt + Lf - u
     // constructor
     STRPDE() = default;
-    STRPDE(const pde_ptr& pde, Sampling s, const DMatrix<double>& time) : Base(pde, s, time) {};
+    STRPDE(const pde_ptr& pde, Sampling s, const DMatrix<double>& time) : Base(pde, s, time) { pde_.init(); };
 
     // redefine SpaceTimeParabolicBase properties affected by iterative approach
     void tensorize_psi() { return; } // avoid tensorization of \Psi matrix
     void init_regularization() {
+        pde_.init();
         // compute time step (assuming equidistant points)
         DeltaT_ = time_[1] - time_[0];
         u_ = pde_.force();   // compute forcing term
@@ -274,9 +276,8 @@ class STRPDE<SpaceTimeParabolic, iterative> :
     const SpMatrix<double>& R1() const { return pde_.stiff(); }   // discretization of differential operator L
     std::size_t n_basis() const { return pde_.n_dofs(); }         // number of basis functions
 
-    void init_model() { return; };          // update model object in case of **structural** changes in its definition
-    void update_to_weights() { return; };   // update model object in case of changes in the weights matrix
-    void solve() {                          // finds a solution to the smoothing problem
+    void init_model() { return; };
+    void solve() {
         // compute starting point (f^(k,0), g^(k,0)) k = 1 ... m for iterative minimization of functional J(f,g)
         A_ = SparseBlockMatrix<double, 2, 2>(
           PsiTD() * Psi(),   lambda_D() * R1().transpose(),

fdaPDE/models/space_only_base.h

@@ -18,7 +18,9 @@
 #define __SPACE_ONLY_BASE_H__
 
 #include <fdaPDE/utils.h>
+#include <fdaPDE/pde.h>
 #include "model_base.h"
+using fdapde::core::pde_ptr;
 
 namespace fdapde {
 namespace models {
@@ -27,36 +29,40 @@ namespace models {
 template <typename Model> class SpaceOnlyBase : public ModelBase<Model> {
    protected:
     typedef ModelBase<Model> Base;
-    using Base::pde_;      // regularizing PDE
-    using Base::model;     // underlying model object
+    using Base::model;   // underlying model object
     static constexpr int n_lambda = n_smoothing_parameters<SpaceOnly>::value;
-
+
+    pde_ptr pde_ {};       // differential penalty in space Lf - u
     SpMatrix<double> P_;   // discretization of penalty term: R1^T*R0^{-1}*R1
     SVector<n_lambda> lambda_ = SVector<n_lambda>::Zero();
    public:
     using Base::lambda;       // dynamic sized smoothing parameter vector
     using Base::set_lambda;   // dynamic sized setter for \lambda
     // constructor
     SpaceOnlyBase() = default;
-    SpaceOnlyBase(const pde_ptr& pde) : ModelBase<Model>(pde) {};
-    void init_regularization() { return; }   // do nothing
-
+    SpaceOnlyBase(const pde_ptr& space_penalty) : pde_(space_penalty) {};
+    void init_regularization() { pde_.init(); } 
     // setters
     void set_lambda(const SVector<n_lambda>& lambda) {
         if(lambda_ == lambda) return;
         model().runtime().set(runtime_status::is_lambda_changed);
         lambda_ = lambda;
     }
     void set_lambda_D(double lambda_D) { set_lambda(SVector<n_lambda>(lambda_D)); }
+    void set_penalty(const pde_ptr& pde) {
+        pde_ = pde;
+        model().runtime().set(runtime_status::require_penalty_init);
+    }
     // getters
     SVector<n_lambda> lambda() const { return lambda_; }
-    double lambda_D() const { return lambda_[0]; }                 // smoothing parameter
-    const SpMatrix<double>& R0() const { return pde_.mass(); }     // mass matrix in space
+    double lambda_D() const { return lambda_[0]; }
+    const SpMatrix<double>& R0() const { return pde_.mass(); }     // mass matrix
     const SpMatrix<double>& R1() const { return pde_.stiff(); }    // discretization of differential operator L
     const DMatrix<double>& u() const { return pde_.force(); }      // discretization of forcing term u
     inline std::size_t n_temporal_locs() const { return 1; }       // number of time instants
     std::size_t n_basis() const { return pde_.n_dofs(); };         // number of basis functions
     std::size_t n_spatial_basis() const { return n_basis(); }      // number of basis functions in space
+    const pde_ptr& pde() const { return pde_; }                    // regularizing term Lf - u
 
     // computes and cache R1^T*R0^{-1}*R1. Returns the discretized penalty P = \lambda_D*(R1^T*R0^{-1}*R1)
     auto P() {
@@ -66,8 +72,7 @@ template <typename Model> class SpaceOnlyBase : public ModelBase<Model> {
             P_ = R1().transpose() * invR0_.solve(R1());   // R1^T*R0^{-1}*R1
         }
         return lambda_D() * P_;
-    }
-
+    }  
     // destructor
     virtual ~SpaceOnlyBase() = default;
 };

fdaPDE/models/space_time_base.h

@@ -30,7 +30,6 @@ class SpaceTimeBase : public ModelBase<Model> {
    protected:
     typedef ModelBase<Model> Base;
     using Base::model;     // underlying model object
-    using Base::pde_;      // regularizing PDE
     static constexpr int n_lambda = n_smoothing_parameters<RegularizationType>::value;
 
     DVector<double> time_;   // time domain [0, T]
@@ -40,8 +39,7 @@ class SpaceTimeBase : public ModelBase<Model> {
     using Base::set_lambda;   // dynamic sized setter for \lambda
     // constructor
     SpaceTimeBase() = default;
-    SpaceTimeBase(const pde_ptr& pde, const DVector<double>& time) : ModelBase<Model>(pde), time_(time) {};
-
+    SpaceTimeBase(const DVector<double>& time) : Base(), time_(time) {};
     // setters
     void set_lambda(const SVector<n_lambda>& lambda) {
         if(lambda_ == lambda) return;
@@ -57,9 +55,7 @@ class SpaceTimeBase : public ModelBase<Model> {
     inline double lambda_T() const { return lambda_[1]; }
     const DVector<double>& time_domain() const { return time_; }           // number of nodes in time
     const DVector<double>& time_locs() const { return time_; }             // time locations where we have observations
-    inline std::size_t n_temporal_locs() const { return time_.rows(); }    // number of time instants
-    std::size_t n_spatial_basis() const { return pde_.n_dofs(); }          // number of basis functions in space
-
+    inline std::size_t n_temporal_locs() const { return time_.rows(); }    // number of time instants  
     // destructor
     virtual ~SpaceTimeBase() = default;
 };

fdaPDE/models/space_time_parabolic_base.h

@@ -17,12 +17,11 @@
 #ifndef __SPACE_TIME_PARABOLIC_BASE_H__
 #define __SPACE_TIME_PARABOLIC_BASE_H__
 
-#include <fdaPDE/finite_elements.h>
+#include <fdaPDE/pde.h>
 #include <fdaPDE/linear_algebra.h>
 #include <fdaPDE/utils.h>
-using fdapde::core::is_parabolic;
+using fdapde::core::pde_ptr;
 using fdapde::core::Kronecker;
-using fdapde::core::SparseKroneckerTensorProduct;
 
 #include "space_time_base.h"
 
@@ -33,6 +32,7 @@ namespace models {
 template <typename Model>
 class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
    protected:
+    pde_ptr pde_ {};   // parabolic differential penalty df/dt + Lf - u
     // let m the number of time points
     DMatrix<double> s_;     // N x 1 initial condition vector
     DMatrix<double> u_;     // discretized forcing [1/DeltaT * (u_1 + R_0*s) \ldots u_n]
@@ -51,22 +51,21 @@ class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
     using Base::lambda_D;   // smoothing parameter in space
     using Base::lambda_T;   // smoothing parameter in time
     using Base::model;      // underlying model object
-    using Base::pde_;       // regularizing term in space
     using Base::time_;      // time interval [0,T]
     using Base::df_;        // model's data
-
     // constructor
     SpaceTimeParabolicBase() = default;
-    SpaceTimeParabolicBase(const pde_ptr& pde, const DVector<double>& time) : Base(pde, time) { }
+    SpaceTimeParabolicBase(const pde_ptr& parabolic_penalty, const DVector<double>& time) :
+      pde_(parabolic_penalty), Base(time) { }
     // init data structure related to parabolic regularization
     void init_regularization() {
+        pde_.init(); 
         std::size_t m_ = time_.rows();   // number of time points
         DeltaT_ = time_[1] - time_[0];   // time step (assuming equidistant points)
 
         // assemble once the m x m identity matrix and cache for fast access
         Im_.resize(m_, m_);
         Im_.setIdentity();
-
         // assemble matrix associated with derivation in time L_
         // [L_]_{ii} = 1/DeltaT for i \in {1 ... m} and [L_]_{i,i-1} = -1/DeltaT for i \in {1 ... m-1}
         std::vector<fdapde::Triplet<double>> triplet_list;
@@ -78,7 +77,6 @@ class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
             triplet_list.emplace_back(i, i, invDeltaT);
             triplet_list.emplace_back(i, i - 1, -invDeltaT);
         }
-        // finalize construction
         L_.resize(m_, m_);
         L_.setFromTriplets(triplet_list.begin(), triplet_list.end());
         L_.makeCompressed();
@@ -89,11 +87,29 @@ class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
 	u_ = pde_.force();
 	u_.block(0, 0, model().n_basis(), 1) += (1.0 / DeltaT_) * (pde_.mass() * s_);
     }
-
+    // setters
+    void set_penalty(const pde_ptr& pde) {
+        pde_ = pde;
+        model().runtime().set(runtime_status::require_penalty_init);
+    }
+    // shift = true, cause the removal of the first time instant of data, in case it has been used to estimate the IC
+    void set_initial_condition(const DMatrix<double>& s, bool shift = true) {
+        s_ = s;
+        if (shift) { // left shrink time domain by one step
+            std::size_t m = time_.rows();       // number of time instants
+            time_ = time_.tail(m - 1).eval();   // correct time interval [0,T] (eval() to avoid aliasing)
+            pde_.set_forcing(pde_.forcing_data().rightCols(m - 1));
+            model().runtime().set(runtime_status::require_penalty_init);   // force pde (re-)initialization
+            // remove from data the first time instant, reindex points
+            model().set_data(df_.tail(model().n_spatial_locs()).extract(), true);
+        }
+    }
     // getters
+    const pde_ptr& pde() const { return pde_; }   // regularizing term df/dt + Lf - u
     const SpMatrix<double>& R0() const { return R0_; }
     const SpMatrix<double>& R1() const { return R1_; }
     std::size_t n_basis() const { return pde_.n_dofs(); }   // number of basis functions
+    std::size_t n_spatial_basis() const { return pde_.n_dofs(); }
     const SpMatrix<double>& L() const { return L_; }
     const DMatrix<double>& u() const { return u_; }   // discretized force corrected by initial conditions
     const DMatrix<double>& s() { return s_; }         // initial condition
@@ -108,19 +124,6 @@ class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
         }
         return lambda_D() * (R1() + lambda_T() * penT_).transpose() * invR0_.solve(R1() + lambda_T() * penT_);
     }
-    // setters
-    // shift = true, cause the removal of the first time instant of data, in case it has been used to estimate the IC
-    void set_initial_condition(const DMatrix<double>& s, bool shift = true) {
-        s_ = s;
-        if (shift) { // left shrink time domain by one step
-            std::size_t m = time_.rows();       // number of time instants
-            time_ = time_.tail(m - 1).eval();   // correct time interval [0,T] (eval() to avoid aliasing)
-            pde_.set_forcing(pde_.forcing_data().rightCols(m - 1));
-            model().runtime().set(runtime_status::require_pde_init);   // force pde (re-)initialization
-            // remove from data the first time instant, reindex points
-            model().set_data(df_.tail(model().n_spatial_locs()).extract(), true);
-        }
-    }
 
     // destructor
     virtual ~SpaceTimeParabolicBase() = default;