Skip to content

Commit

Permalink
decoupled time penalty from space_time_separable_base (possible to su…
Browse files Browse the repository at this point in the history 
…pply 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
  • Loading branch information
@AlePalu
AlePalu committed Dec 22, 2023
1 parent 3fa5006 commit 3b53a77
Show file tree
Hide file tree
Showing 12 changed files with 168 additions and 178 deletions.
2 changes: 1 addition & 1 deletion fdaPDE/core
Submodule core updated 24 files
+3 −4 fdaPDE/finite_elements/basis/lagrangian_basis.h
+4 −4 fdaPDE/finite_elements/solvers/fem_linear_elliptic_solver.h
+2 −2 fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h
+10 −10 fdaPDE/finite_elements/solvers/fem_solver_base.h
+7 −7 fdaPDE/mesh/element.h
+128 −3 fdaPDE/mesh/mesh.h
+11 −11 fdaPDE/pde/pde.h
+33 −0 fdaPDE/pde/symbols.h
+3 −0 fdaPDE/splines.h
+49 −6 fdaPDE/splines/basis/spline_basis.h
+53 −0 fdaPDE/splines/solvers/spline_linear_elliptic_solver.h
+86 −0 fdaPDE/splines/solvers/spline_solver_base.h
+37 −0 fdaPDE/splines/solvers/spline_solver_selector.h
+34 −41 fdaPDE/splines/spline_assembler.h
+5 −0 fdaPDE/splines/spline_symbols.h
+118 −106 fdaPDE/utils/integration/integrator.h
+6 −5 fdaPDE/utils/integration/integrator_tables.h
+1 −1 test/run_tests.sh
+1 −1 test/src/fem_operators_test.cpp
+126 −129 test/src/fem_pde_test.cpp
+4 −3 test/src/integration_test.cpp
+29 −7 test/src/mesh_test.cpp
+37 −7 test/src/point_location_test.cpp
+22 −39 test/src/spline_test.cpp
7 changes: 5 additions & 2 deletions fdaPDE/models/regression/fpirls.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,12 @@ template <typename Model_> class FPIRLS {
if constexpr (!is_space_time<Model_>::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<Model_>::value) { solver_.set_initial_condition(m_->s(), false); }
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_penalty(), m_->sampling());
solver_.set_temporal_locations(m_->time_locs());
}
}
Expand Down
22 changes: 12 additions & 10 deletions fdaPDE/models/regression/gsrpde.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,18 +44,20 @@ class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, Regularization
// constructor
GSRPDE() = default;
// space-only constructor
template <
typename U = RegularizationType,
typename std::enable_if<std::is_same<U, SpaceOnly>::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>(this, tol_, max_iter_);
};
// space-time constructor
template <
typename U = RegularizationType,
typename std::enable_if<!std::is_same<U, SpaceOnly>::value, int>::type = 0>
GSRPDE(const pde_ptr& pde, const DVector<double>& 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>(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) {
fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
};

Expand Down
17 changes: 8 additions & 9 deletions fdaPDE/models/regression/regression_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::is_same<typename U::RegularizationType, SpaceOnly>::value, int>::type = 0>
RegressionBase(const pde_ptr& pde, Sampling s) : Base(pde), SamplingBase<Model>(s) {};
// space-time constructor
template <
typename U = Model, // fake type to enable substitution in SFINAE
typename std::enable_if<!std::is_same<typename U::RegularizationType, SpaceOnly>::value, int>::type = 0>
RegressionBase(const pde_ptr& pde, Sampling s, const DVector<double>& time) :
fdapde_enable_constructor_if(is_space_only, 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)
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) {};
// 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_; }
Expand Down
36 changes: 18 additions & 18 deletions fdaPDE/models/regression/srpde.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,24 +52,24 @@ class SRPDE : public RegressionBase<SRPDE, SpaceOnly> {
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<double, 2, 2>(
-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<double, 2, 2>(
-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() {
Expand Down
56 changes: 29 additions & 27 deletions fdaPDE/models/regression/strpde.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,31 +69,33 @@ class STRPDE<SpaceTimeSeparable, monolithic> :
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<double>& 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<double, 2, 2>(
-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<double, 2, 2>(
-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<double> sol; // room for problem' solution

DVector<double> 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();
Expand Down Expand Up @@ -155,7 +157,7 @@ class STRPDE<SpaceTimeParabolic, monolithic> :

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<double, 2, 2>(
-PsiTD() * W() * Psi(), lambda_D() * (R1() + lambda_T() * L_).transpose(),
lambda_D() * (R1() + lambda_T() * L_), lambda_D() * R0());
Expand Down Expand Up @@ -265,12 +267,12 @@ class STRPDE<SpaceTimeParabolic, iterative> :
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<double>& R0() const { return pde_.R0(); } // mass matrix in space
const SpMatrix<double>& 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<double>& R0() const { return pde_.mass(); } // mass matrix in space
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
Expand Down
7 changes: 4 additions & 3 deletions fdaPDE/models/space_only_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,15 +43,16 @@ template <typename Model> class SpaceOnlyBase : public ModelBase<Model> {

// 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)); }
// getters
SVector<n_lambda> lambda() const { return lambda_; }
double lambda_D() const { return lambda_[0]; } // smoothing parameter
const SpMatrix<double>& R0() const { return pde_.R0(); } // mass matrix in space
const SpMatrix<double>& R1() const { return pde_.R1(); } // discretization of differential operator L
const SpMatrix<double>& R0() const { return pde_.mass(); } // mass matrix in space
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
Expand All @@ -61,7 +62,7 @@ template <typename Model> class SpaceOnlyBase : public ModelBase<Model> {
auto P() {
if (is_empty(P_)) {
fdapde::SparseLU<SpMatrix<double>> invR0_;
invR0_.compute(pde_.R0());
invR0_.compute(pde_.mass());
P_ = R1().transpose() * invR0_.solve(R1()); // R1^T*R0^{-1}*R1
}
return lambda_D() * P_;
Expand Down
10 changes: 7 additions & 3 deletions fdaPDE/models/space_time_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,13 @@ class SpaceTimeBase : public ModelBase<Model> {
SpaceTimeBase(const pde_ptr& pde, const DVector<double>& time) : ModelBase<Model>(pde), time_(time) {};

// setters
void set_lambda(const SVector<n_lambda>& 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<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, lambda_[1])); }
void set_lambda_T(double lambda_T) { set_lambda(SVector<n_lambda>(lambda_[0], lambda_T)); }
void set_time_domain(const DVector<double>& time) { time_ = time; }
// getters
SVector<n_lambda> lambda() const { return lambda_; }
Expand Down
8 changes: 4 additions & 4 deletions fdaPDE/models/space_time_parabolic_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,11 +83,11 @@ class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
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
Expand All @@ -104,7 +104,7 @@ class SpaceTimeParabolicBase : public SpaceTimeBase<Model, SpaceTimeParabolic> {
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_);
}
Expand Down
Loading

0 comments on commit 3b53a77

Please sign in to comment.