From 3b53a77757e346aaf7a55dc5b07be09c2ed0f8c5 Mon Sep 17 00:00:00 2001 From: AlePalu Date: Fri, 22 Dec 2023 13:04:14 +0100 Subject: [PATCH] decoupled time penalty from space_time_separable_base (possible to supply any PDE operator, built on a B-Spline basis, as time-penalty in a spatio-temporal separable problem), exploit introduced Mesh<1,1> as time domain, models adapted to the new interface --- fdaPDE/core | 2 +- fdaPDE/models/regression/fpirls.h | 7 +- fdaPDE/models/regression/gsrpde.h | 22 +++--- fdaPDE/models/regression/regression_base.h | 17 ++-- fdaPDE/models/regression/srpde.h | 36 ++++----- fdaPDE/models/regression/strpde.h | 56 ++++++------- fdaPDE/models/space_only_base.h | 7 +- fdaPDE/models/space_time_base.h | 10 ++- fdaPDE/models/space_time_parabolic_base.h | 8 +- fdaPDE/models/space_time_separable_base.h | 71 ++++++----------- test/src/gsrpde_test.cpp | 18 ++--- test/src/strpde_test.cpp | 92 +++++++++++----------- 12 files changed, 168 insertions(+), 178 deletions(-) diff --git a/fdaPDE/core b/fdaPDE/core index bde70420..b2c6b62d 160000 --- a/fdaPDE/core +++ b/fdaPDE/core @@ -1 +1 @@ -Subproject commit bde70420d2452e6f69cedbcbde11870281a5c27a +Subproject commit b2c6b62d4b7066989d2f015f121598d23083f868 diff --git a/fdaPDE/models/regression/fpirls.h b/fdaPDE/models/regression/fpirls.h index a5880de0..8f74f4d9 100644 --- a/fdaPDE/models/regression/fpirls.h +++ b/fdaPDE/models/regression/fpirls.h @@ -58,9 +58,12 @@ template class FPIRLS { if constexpr (!is_space_time::value) // space-only 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); } + if constexpr (is_space_time_parabolic::value) { + solver_ = SolverType(m_->pde(), m_->sampling(), m_->time_domain()); + solver_.set_initial_condition(m_->s(), false); + } if constexpr (is_space_time_separable::value) { + solver_ = SolverType(m_->pde(), m_->time_penalty(), m_->sampling()); solver_.set_temporal_locations(m_->time_locs()); } } diff --git a/fdaPDE/models/regression/gsrpde.h b/fdaPDE/models/regression/gsrpde.h index 2bddf52c..a0211100 100644 --- a/fdaPDE/models/regression/gsrpde.h +++ b/fdaPDE/models/regression/gsrpde.h @@ -44,18 +44,20 @@ class GSRPDE : public RegressionBase, Regularization // constructor GSRPDE() = default; // space-only constructor - 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) { + fdapde_enable_constructor_if(is_space_only, This) + 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) { + // 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, tol_, max_iter_); + }; + fdapde_enable_constructor_if(is_space_time_parabolic, This) + GSRPDE(const pde_ptr& space_penalty, const DVector& time, Sampling s, const Distribution& distr) : + Base(space_penalty, s, time), distr_(distr) { fpirls_ = FPIRLS(this, tol_, max_iter_); }; diff --git a/fdaPDE/models/regression/regression_base.h b/fdaPDE/models/regression/regression_base.h index 16f40305..0174c360 100644 --- a/fdaPDE/models/regression/regression_base.h +++ b/fdaPDE/models/regression/regression_base.h @@ -66,15 +66,14 @@ class RegressionBase : RegressionBase() = default; // space-only constructor - template < - typename U = Model, // fake type to enable substitution in SFINAE - typename std::enable_if::value, int>::type = 0> - RegressionBase(const pde_ptr& pde, Sampling s) : Base(pde), SamplingBase(s) {}; - // space-time constructor - template < - typename U = Model, // fake type to enable substitution in SFINAE - typename std::enable_if::value, int>::type = 0> - RegressionBase(const pde_ptr& pde, Sampling s, const DVector& time) : + fdapde_enable_constructor_if(is_space_only, Model) RegressionBase(const pde_ptr& pde, Sampling s) : + Base(pde), SamplingBase(s) {}; + // 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(s) { }; + fdapde_enable_constructor_if(is_space_time_parabolic, Model) + RegressionBase(const pde_ptr& pde, Sampling s, const DVector& time) : Base(pde, time), SamplingBase(s) {}; // copy constructor, copy only pde object (as a consequence also the problem domain) RegressionBase(const RegressionBase& rhs) : Base(rhs), SamplingBase(rhs) { pde_ = rhs.pde_; } diff --git a/fdaPDE/models/regression/srpde.h b/fdaPDE/models/regression/srpde.h index 964862bd..979c179f 100644 --- a/fdaPDE/models/regression/srpde.h +++ b/fdaPDE/models/regression/srpde.h @@ -52,24 +52,24 @@ class SRPDE : public RegressionBase { SRPDE(const pde_ptr& pde, Sampling s) : Base(pde, s) {}; void init_model() { - // a change in the smoothing parameter must reset the whole linear system - if(runtime().query(runtime_status::is_lambda_changed)) { - // assemble system matrix for nonparameteric part - A_ = SparseBlockMatrix( - -PsiTD() * W() * Psi(), lambda_D() * R1().transpose(), - lambda_D() * R1(), lambda_D() * R0() ); - invA_.compute(A_); - // prepare rhs of linear system - b_.resize(A_.rows()); - b_.block(n_basis(), 0, n_basis(), 1) = lambda_D() * u(); - return; - } - if(runtime().query(runtime_status::require_W_update)) { - // adjust north-west block of matrix A_ only - A_.block(0, 0) = -PsiTD() * W() * Psi(); - invA_.compute(A_); - return; - } + // a change in the smoothing parameter must reset the whole linear system + if (runtime().query(runtime_status::is_lambda_changed)) { + // assemble system matrix for nonparameteric part + A_ = SparseBlockMatrix( + -PsiTD() * W() * Psi(), lambda_D() * R1().transpose(), + lambda_D() * R1(), lambda_D() * R0() ); + invA_.compute(A_); + // prepare rhs of linear system + b_.resize(A_.rows()); + b_.block(n_basis(), 0, n_basis(), 1) = lambda_D() * u(); + return; + } + if (runtime().query(runtime_status::require_W_update)) { + // adjust north-west block of matrix A_ only + A_.block(0, 0) = -PsiTD() * W() * Psi(); + invA_.compute(A_); + return; + } } // finds a solution to the smoothing problem void solve() { diff --git a/fdaPDE/models/regression/strpde.h b/fdaPDE/models/regression/strpde.h index 44ae7e46..c527e87e 100644 --- a/fdaPDE/models/regression/strpde.h +++ b/fdaPDE/models/regression/strpde.h @@ -69,31 +69,33 @@ class STRPDE : using Base::P1; // time penalization matrix: [P1_]_{ij} = \int_{[0,T]} (\phi_i)_tt*(\phi_j)_tt // constructor STRPDE() = default; - STRPDE(const pde_ptr& pde, Sampling s, const DMatrix& time) : Base(pde, s, time) {}; - - void init_model() { // update model object in case of **structural** changes in its definition - // assemble system matrix for the nonparameteric part of the model - if (is_empty(K_)) K_ = Kronecker(P1(), pde().R0()); - A_ = SparseBlockMatrix( - -PsiTD() * W() * Psi() - lambda_T() * K_, lambda_D() * R1().transpose(), lambda_D() * R1(), - lambda_D() * R0()); - // cache system matrix for reuse - invA_.compute(A_); - // prepare rhs of linear system - b_.resize(A_.rows()); - b_.block(A_.rows() / 2, 0, A_.rows() / 2, 1) = lambda_D() * u(); - return; - } - void update_to_weights() { // update model object in case of changes in the weights matrix - // adjust north-west block of matrix A_ and factorize - A_.block(0, 0) = -PsiTD() * W() * Psi() - lambda_T() * K_; - invA_.compute(A_); - return; + STRPDE(const pde_ptr& space_penalty, const pde_ptr& time_penalty, Sampling s) : + Base(space_penalty, time_penalty, s) {}; + + void init_model() { + // a change in the smoothing parameter must reset the whole linear system + if (runtime().query(runtime_status::is_lambda_changed)) { + // assemble system matrix for the nonparameteric part + if (is_empty(K_)) K_ = Kronecker(P1(), pde().mass()); + A_ = SparseBlockMatrix( + -PsiTD() * W() * Psi() - lambda_T() * K_, lambda_D() * R1().transpose(), + lambda_D() * R1(), lambda_D() * R0() ); + invA_.compute(A_); + // prepare rhs of linear system + b_.resize(A_.rows()); + b_.block(A_.rows() / 2, 0, A_.rows() / 2, 1) = lambda_D() * u(); + return; + } + if (runtime().query(runtime_status::require_W_update)) { + // adjust north-west block of matrix A_ and factorize + A_.block(0, 0) = -PsiTD() * W() * Psi() - lambda_T() * K_; + invA_.compute(A_); + return; + } } void solve() { // finds a solution to the smoothing problem BLOCK_FRAME_SANITY_CHECKS; - DVector sol; // room for problem' solution - + DVector sol; // room for problem' solution if (!Base::has_covariates()) { // nonparametric case // update rhs of STR-PDE linear system b_.block(0, 0, A_.rows() / 2, 1) = -PsiTD() * W() * y(); @@ -155,7 +157,7 @@ class STRPDE : void init_model() { // update model object in case of **structural** changes in its definition // assemble system matrix for the nonparameteric part of the model - if (is_empty(L_)) L_ = Kronecker(L(), pde().R0()); + if (is_empty(L_)) L_ = Kronecker(L(), pde().mass()); A_ = SparseBlockMatrix( -PsiTD() * W() * Psi(), lambda_D() * (R1() + lambda_T() * L_).transpose(), lambda_D() * (R1() + lambda_T() * L_), lambda_D() * R0()); @@ -265,12 +267,12 @@ class STRPDE : DeltaT_ = time_[1] - time_[0]; u_ = pde_.force(); // compute forcing term // correct first n rows of discretized force as (u_1 + R0*s/DeltaT) - u_.block(0, 0, n_basis(), 1) += (1.0 / DeltaT_) * (pde_.R0() * s_); + u_.block(0, 0, n_basis(), 1) += (1.0 / DeltaT_) * (pde_.mass() * s_); } // getters - const SpMatrix& R0() const { return pde_.R0(); } // mass matrix in space - const SpMatrix& R1() const { return pde_.R1(); } // discretization of differential operator L - std::size_t n_basis() const { return pde_.n_dofs(); } // number of basis functions + const SpMatrix& R0() const { return pde_.mass(); } // mass matrix in space + const SpMatrix& 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 diff --git a/fdaPDE/models/space_only_base.h b/fdaPDE/models/space_only_base.h index d6543a09..32bed0cb 100644 --- a/fdaPDE/models/space_only_base.h +++ b/fdaPDE/models/space_only_base.h @@ -43,6 +43,7 @@ template class SpaceOnlyBase : public ModelBase { // setters void set_lambda(const SVector& lambda) { + if(lambda_ == lambda) return; model().runtime().set(runtime_status::is_lambda_changed); lambda_ = lambda; } @@ -50,8 +51,8 @@ template class SpaceOnlyBase : public ModelBase { // getters SVector lambda() const { return lambda_; } double lambda_D() const { return lambda_[0]; } // smoothing parameter - const SpMatrix& R0() const { return pde_.R0(); } // mass matrix in space - const SpMatrix& R1() const { return pde_.R1(); } // discretization of differential operator L + const SpMatrix& R0() const { return pde_.mass(); } // mass matrix in space + const SpMatrix& R1() const { return pde_.stiff(); } // discretization of differential operator L const DMatrix& 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 @@ -61,7 +62,7 @@ template class SpaceOnlyBase : public ModelBase { auto P() { if (is_empty(P_)) { fdapde::SparseLU> invR0_; - invR0_.compute(pde_.R0()); + invR0_.compute(pde_.mass()); P_ = R1().transpose() * invR0_.solve(R1()); // R1^T*R0^{-1}*R1 } return lambda_D() * P_; diff --git a/fdaPDE/models/space_time_base.h b/fdaPDE/models/space_time_base.h index 1fcde680..64c9af51 100644 --- a/fdaPDE/models/space_time_base.h +++ b/fdaPDE/models/space_time_base.h @@ -43,9 +43,13 @@ class SpaceTimeBase : public ModelBase { SpaceTimeBase(const pde_ptr& pde, const DVector& time) : ModelBase(pde), time_(time) {}; // setters - void set_lambda(const SVector& lambda) { lambda_ = lambda; } - void set_lambda_D(double lambda_D) { lambda_[0] = lambda_D; } - void set_lambda_T(double lambda_T) { lambda_[1] = lambda_T; } + void set_lambda(const SVector& 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(lambda_D, lambda_[1])); } + void set_lambda_T(double lambda_T) { set_lambda(SVector(lambda_[0], lambda_T)); } void set_time_domain(const DVector& time) { time_ = time; } // getters SVector lambda() const { return lambda_; } diff --git a/fdaPDE/models/space_time_parabolic_base.h b/fdaPDE/models/space_time_parabolic_base.h index a47f0544..02586eeb 100644 --- a/fdaPDE/models/space_time_parabolic_base.h +++ b/fdaPDE/models/space_time_parabolic_base.h @@ -83,11 +83,11 @@ class SpaceTimeParabolicBase : public SpaceTimeBase { L_.setFromTriplets(triplet_list.begin(), triplet_list.end()); L_.makeCompressed(); // compute tensorized matrices - R0_ = Kronecker(Im_, pde_.R0()); - R1_ = Kronecker(Im_, pde_.R1()); + R0_ = Kronecker(Im_, pde_.mass()); + R1_ = Kronecker(Im_, pde_.stiff()); // correct first n rows of discretized force as (u_1 + R0*s/DeltaT) u_ = pde_.force(); - u_.block(0, 0, model().n_basis(), 1) += (1.0 / DeltaT_) * (pde_.R0() * s_); + u_.block(0, 0, model().n_basis(), 1) += (1.0 / DeltaT_) * (pde_.mass() * s_); } // getters @@ -104,7 +104,7 @@ class SpaceTimeParabolicBase : public SpaceTimeBase { auto P() { if (is_empty(pen_)) { // compute once and cache result invR0_.compute(R0()); - penT_ = Kronecker(L_, pde_.R0()); + penT_ = Kronecker(L_, pde_.mass()); } return lambda_D() * (R1() + lambda_T() * penT_).transpose() * invR0_.solve(R1() + lambda_T() * penT_); } diff --git a/fdaPDE/models/space_time_separable_base.h b/fdaPDE/models/space_time_separable_base.h index 5a678172..f1e1c038 100644 --- a/fdaPDE/models/space_time_separable_base.h +++ b/fdaPDE/models/space_time_separable_base.h @@ -36,18 +36,15 @@ template class SpaceTimeSeparableBase : public SpaceTimeBase { protected: // let \phi_i the i-th basis function in time - SpMatrix P0_; // mass matrix in time: [P0_]_{ij} = \int_{[0,T]} \phi_i*\phi_j - SpMatrix P1_; // penalty matrix in time: [P1_]_{ij} = \int_{[0,T]} (\phi_i)_tt*(\phi_j)_tt SpMatrix Phi_; // [Phi_]_{ij} = \phi_i(t_j) DVector u_; // stacked discretized forcing [u_1 \ldots u_n, \ldots, u_1 \ldots u_n] - SpMatrix R0_; // P0_ \kron R0 (R0: discretization of the identity operator) - SpMatrix R1_; // P0_ \kron R1 (R1: discretization of the differential operator L in the regularizing PDE) + SpMatrix R0_; // P0 \kron R0 (R0: mass matrix in space) + SpMatrix R1_; // P0 \kron R1 (R1: stiff matrix in space) - using TimeBasis = SplineBasis<3>; // basis used for discretization in time - TimeBasis basis_; + pde_ptr time_penalty_ {}; // time regularizing term DVector time_locs_; // time instants t_1, ..., t_m - SpMatrix P_D_; // discretization of space regularization: (R1^T*R0^{-1}*R1) \kron Rt - SpMatrix P_T_; // discretization of time regularization: (R0 \kron Pt) + SpMatrix P_D_; // discretization of space regularization: (R1^T*R0^{-1}*R1) \kron P0 + SpMatrix P_T_; // discretization of time regularization: (R0 \kron P1) public: typedef SpaceTimeBase Base; using Base::lambda_D; // smoothing parameter in space @@ -58,62 +55,45 @@ class SpaceTimeSeparableBase : public SpaceTimeBase { // constructor SpaceTimeSeparableBase() = default; - SpaceTimeSeparableBase(const pde_ptr& pde, const DVector& time) : Base(pde, time) { } + SpaceTimeSeparableBase(const pde_ptr& space_penalty, const pde_ptr& time_penalty) : + Base(space_penalty, time_penalty.dof_coords()), time_penalty_(time_penalty) { } // init data structure related to separable regularization void init_regularization() { - basis_ = TimeBasis(time_); - + time_penalty_.init(); // initialize time-penalty operator // compute \Phi = [\Phi]_{ij} = \phi_i(t_j) // assume time instants equal to time nodes if not provided DVector time_locs = is_empty(time_locs_) ? time_ : time_locs_; - int m = time_locs.rows(); - int M = basis_.size(); - Phi_.resize(m, M); - - // start matrix assembly - std::vector> triplet_list; - triplet_list.reserve(m * M); - for (int i = 0; i < M; ++i) { - for (int j = 0; j < m; ++j) { // evaluate basis function at given m time locations - triplet_list.emplace_back(j, i, basis_[i](SVector<1>(time_locs[j]))); - } - } - // finalize construction - Phi_.setFromTriplets(triplet_list.begin(), triplet_list.end()); - Phi_.prune(0.0); // remove zeros - Phi_.makeCompressed(); - - // discretize operators in time using spline-approximation - core::IntegratorTable<1, 3, core::GaussLegendre> quadrature; - Assembler, TimeBasis, decltype(quadrature)> assembler(time_, quadrature); - P0_ = assembler.discretize_operator(reaction(1.0)); // mass matrix - P1_ = assembler.discretize_operator(-bilaplacian()); - + Phi_ = time_penalty_.eval_basis(core::eval::pointwise, time_locs)->Psi; // compute tensorized matrices - R0_ = Kronecker(P0_, pde_.R0()); - R1_ = Kronecker(P0_, pde_.R1()); + R0_ = Kronecker(time_penalty_.mass(), pde_.mass()); + R1_ = Kronecker(time_penalty_.mass(), pde_.stiff()); } // setters void set_temporal_locations(const DVector& time_locations) { time_locs_ = time_locations; } + void set_time_penalty(const pde_ptr& time_penalty) { + time_penalty_ = time_penalty; + model().runtime().set(runtime_status::require_penalty_init); + } // getters const SpMatrix& R0() const { return R0_; } const SpMatrix& R1() const { return R1_; } - std::size_t n_basis() const { return pde_.n_dofs() * basis_.size(); } // number of basis functions - std::size_t n_temporal_basis() const { return basis_.size(); } // number of time basis functions + std::size_t n_basis() const { return pde_.n_dofs() * time_penalty_.n_dofs(); } // number of basis functions + std::size_t n_temporal_basis() const { return time_penalty_.n_dofs(); } // number of time basis functions // matrices proper of separable regularization - const SpMatrix& P0() const { return P0_; } - const SpMatrix& P1() const { return P1_; } + const SpMatrix& P0() const { return time_penalty_.mass(); } + const SpMatrix& P1() const { return time_penalty_.stiff(); } const SpMatrix& Phi() const { return Phi_; } inline std::size_t n_temporal_locs() const { return is_empty(time_locs_) ? time_.rows() : time_locs_.rows(); } const DVector& time_locs() const { return is_empty(time_locs_) ? time_ : time_locs_; } - + const pde_ptr& time_penalty() const { return time_penalty_; } + // return stacked version of discretized forcing field const DVector& u() { if (is_empty(u_)) { // compute once and cache std::size_t N = Base::n_spatial_basis(); u_.resize(n_basis()); // in separable regularization PDE doesn't depend on time. stack forcing term m times - for (std::size_t i = 0; i < basis_.size(); ++i) { u_.segment(i * N, N) = pde_.force(); } + for (std::size_t i = 0; i < n_temporal_basis(); ++i) { u_.segment(i * N, N) = pde_.force(); } } return u_; } @@ -123,9 +103,10 @@ class SpaceTimeSeparableBase : public SpaceTimeBase { auto P() { if (is_empty(P_D_)) { // compute once and cache result fdapde::SparseLU> invR0_; - invR0_.compute(pde_.R0()); - P_D_ = Kronecker(pde_.R1().transpose() * invR0_.solve(pde_.R1()), P0_); // (R1^T*R0^{-1}*R1) \kron P0 - P_T_ = Kronecker(pde_.R0(), P1_); // (R0 \kron P1) + invR0_.compute(pde_.mass()); + P_D_ = + Kronecker(pde_.stiff().transpose() * invR0_.solve(pde_.stiff()), P0()); // (R1^T*R0^{-1}*R1) \kron P0 + P_T_ = Kronecker(pde_.mass(), P1()); // (R0 \kron P1) } return lambda_D() * P_D_ + lambda_T() * P_T_; } diff --git a/test/src/gsrpde_test.cpp b/test/src/gsrpde_test.cpp index 362b5675..30fd34c0 100644 --- a/test/src/gsrpde_test.cpp +++ b/test/src/gsrpde_test.cpp @@ -191,24 +191,24 @@ TEST(gsrpde_test, laplacian_nonparametric_samplingatlocations_gamma) { // time penalization: separable (mass penalization) // distribution: gamma TEST(gsrpde_test, laplacian_semiparametric_samplingatlocations_separable_monolithic_gamma) { - // define temporal domain - DVector time_mesh; - time_mesh.resize(4); - for (std::size_t i = 0; i < 4; ++i) time_mesh[i] = (1. / 3) * i; - // define spatial domain + // define temporal and spatial domain + Mesh<1, 1> time_mesh(0, 1, 3); MeshLoader domain("c_shaped"); // import data from files DMatrix locs = read_csv("../data/models/gsrpde/2D_test5/locs.csv"); DMatrix y = read_csv("../data/models/gsrpde/2D_test5/y.csv" ); DMatrix X = read_csv("../data/models/gsrpde/2D_test5/X.csv" ); - // define regularizing PDE - auto L = -laplacian(); + // define regularizing PDE in space + auto Ld = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + PDE, decltype(Ld), DMatrix, FEM, fem_order<1>> space_penalty(domain.mesh, Ld, u); + // define regularizing PDE in time + auto Lt = -bilaplacian(); + PDE, decltype(Lt), DMatrix, SPLINE, spline_order<3>> time_penalty(time_mesh, Lt); // define model double lambda_D = std::pow(0.1, 2.5); double lambda_T = std::pow(0.1, 2.5); - GSRPDE model(problem, time_mesh, Sampling::pointwise, Gamma()); + GSRPDE model(space_penalty, time_penalty, Sampling::pointwise, Gamma()); model.set_lambda_D(lambda_D); model.set_lambda_T(lambda_T); model.set_spatial_locations(locs); diff --git a/test/src/strpde_test.cpp b/test/src/strpde_test.cpp index 5dd61729..1e4222fe 100644 --- a/test/src/strpde_test.cpp +++ b/test/src/strpde_test.cpp @@ -26,6 +26,8 @@ using fdapde::core::laplacian; using fdapde::core::MatrixDataWrapper; using fdapde::core::PDE; using fdapde::core::VectorDataWrapper; +using fdapde::core::Mesh; +using fdapde::core::spline_order; #include "../../fdaPDE/models/regression/strpde.h" #include "../../fdaPDE/models/sampling_design.h" @@ -51,23 +53,21 @@ using fdapde::testing::read_csv; // order FE: 1 // time penalization: separable (mass penalization) TEST(strpde_test, laplacian_nonparametric_samplingatnodes_separable_monolithic) { - // define temporal domain - DVector time_mesh; - time_mesh.resize(11); - std::size_t i = 0; - for (double x = 0; x <= 2; x += 0.2, ++i) time_mesh[i] = x; - // define spatial domain + // define temporal and spatial domain + Mesh<1, 1> time_mesh(0, 2, 10); MeshLoader domain("unit_square_coarse"); // import data from files DMatrix y = read_csv("../data/models/strpde/2D_test1/y.csv"); - // define regularizing PDE - auto L = -laplacian(); - DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3 * time_mesh.rows(), 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + // define regularizing PDE in space + auto Ld = -laplacian(); + DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3 * time_mesh.n_nodes(), 1); + PDE, decltype(Ld), DMatrix, FEM, fem_order<1>> space_penalty(domain.mesh, Ld, u); + // define regularizing PDE in time + auto Lt = -bilaplacian(); + PDE, decltype(Lt), DMatrix, SPLINE, spline_order<3>> time_penalty(time_mesh, Lt); // define model - double lambda_D = 0.01; - double lambda_T = 0.01; - STRPDE model(problem, Sampling::mesh_nodes, time_mesh); + double lambda_D = 0.01, lambda_T = 0.01; + STRPDE model(space_penalty, time_penalty, Sampling::mesh_nodes); model.set_lambda_D(lambda_D); model.set_lambda_T(lambda_T); // set model's data @@ -90,24 +90,24 @@ TEST(strpde_test, laplacian_nonparametric_samplingatnodes_separable_monolithic) // order FE: 1 // time penalization: separable (mass penalization) TEST(strpde_test, laplacian_semiparametric_samplingatlocations_separable_monolithic) { - // define temporal domain - DVector time_mesh; - time_mesh.resize(5); - for (std::size_t i = 0; i < 5; ++i) time_mesh[i] = (fdapde::testing::pi / 4) * i; - // define spatial domain + // define temporal and spatial domain + Mesh<1, 1> time_mesh(0, fdapde::testing::pi, 4); MeshLoader domain("c_shaped"); // import data from files DMatrix locs = read_csv("../data/models/strpde/2D_test2/locs.csv"); DMatrix y = read_csv("../data/models/strpde/2D_test2/y.csv"); DMatrix X = read_csv("../data/models/strpde/2D_test2/X.csv"); - // define regularizing PDE - auto L = -laplacian(); + // define regularizing PDE in space + auto Ld = -laplacian(); DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3, 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + PDE, decltype(Ld), DMatrix, FEM, fem_order<1>> space_penalty(domain.mesh, Ld, u); + // define regularizing PDE in time + auto Lt = -bilaplacian(); + PDE, decltype(Lt), DMatrix, SPLINE, spline_order<3>> time_penalty(time_mesh, Lt); // define model double lambda_D = 0.01; double lambda_T = 0.01; - STRPDE model(problem, Sampling::pointwise, time_mesh); + STRPDE model(space_penalty, time_penalty, Sampling::pointwise); model.set_lambda_D(lambda_D); model.set_lambda_T(lambda_T); model.set_spatial_locations(locs); @@ -223,24 +223,24 @@ TEST(strpde_test, laplacian_nonparametric_samplingatnodes_parabolic_iterative) { // order FE: 1 // time penalization: separable (mass penalization) TEST(strpde_test, laplacian_nonparametric_samplingatnodes_timelocations_separable_monolithic) { - // define temporal domain - DVector time_mesh; - time_mesh.resize(11); - std::size_t i = 0; - for (double x = 0; x <= 2; x += 0.2, ++i) time_mesh[i] = x; - // define spatial domain + // define temporal and spatial domain + Mesh<1, 1> time_mesh(0, 2, 10); MeshLoader domain("unit_square_coarse"); // import data from files DMatrix time_locs = read_csv("../data/models/strpde/2D_test5/time_locations.csv"); DMatrix y = read_csv("../data/models/strpde/2D_test5/y.csv"); - // define regularizing PDE - auto L = -laplacian(); - DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3 * time_mesh.rows(), 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + // define regularizing PDE in space + auto Ld = -laplacian(); + DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3 * time_mesh.n_nodes(), 1); + PDE, decltype(Ld), DMatrix, FEM, fem_order<1>> space_penalty(domain.mesh, Ld, u); + // define regularizing PDE in time + auto Lt = -bilaplacian(); + PDE, decltype(Lt), DMatrix, SPLINE, spline_order<3>> time_penalty(time_mesh, Lt); + // define model double lambda_D = 0.01; double lambda_T = 0.01; - STRPDE model(problem, Sampling::mesh_nodes, time_mesh); + STRPDE model(space_penalty, time_penalty, Sampling::mesh_nodes); model.set_lambda_D(lambda_D); model.set_lambda_T(lambda_T); model.set_temporal_locations(time_locs); @@ -266,26 +266,24 @@ TEST(strpde_test, laplacian_nonparametric_samplingatnodes_timelocations_separabl // order FE: 1 // time penalization: separable (mass penalization) TEST(strpde_test, laplacian_nonparametric_samplingatlocations_timelocations_separable_monolithic_missingdata) { - // define temporal domain - DVector time_mesh; - time_mesh.resize(21); - for (std::size_t i = 0; i < 21; ++i) time_mesh[i] = 1.0 / 20 * i; - // define spatial domain and regularizing PDE + // define temporal and spatial domain + Mesh<1, 1> time_mesh(0, 1, 20); MeshLoader domain("c_shaped"); // import data from files DMatrix time_locs = read_csv("../data/models/strpde/2D_test6/time_locations.csv"); DMatrix space_locs = read_csv("../data/models/strpde/2D_test6/locs.csv"); DMatrix y = read_csv("../data/models/strpde/2D_test6/y.csv" ); - // define regularizing PDE - auto L = -laplacian(); - DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3 * time_mesh.rows(), 1); - PDE, FEM, fem_order<1>> problem(domain.mesh, L, u); + // define regularizing PDE in space + auto Ld = -laplacian(); + DMatrix u = DMatrix::Zero(domain.mesh.n_elements() * 3 * time_mesh.n_nodes(), 1); + PDE, decltype(Ld), DMatrix, FEM, fem_order<1>> space_penalty(domain.mesh, Ld, u); + // define regularizing PDE in time + auto Lt = -bilaplacian(); + PDE, decltype(Lt), DMatrix, SPLINE, spline_order<3>> time_penalty(time_mesh, Lt); // define model - double lambda_D = 1e-3; - double lambda_T = 1e-3; - STRPDE model(problem, Sampling::pointwise, time_mesh); - model.set_lambda_D(lambda_D); - model.set_lambda_T(lambda_T); + STRPDE model(space_penalty, time_penalty, Sampling::pointwise); + model.set_lambda_D(1e-3); + model.set_lambda_T(1e-3); model.set_spatial_locations(space_locs); model.set_temporal_locations(time_locs); // set model's data