Skip to content

Commit

Permalink
quantile spatial regression models are now named QSRPDE, GCV for QSRP…
Browse files Browse the repository at this point in the history 
…DE restored. FPIRLS is default constructable and is stored as a model data member
  • Loading branch information
@AlePalu
AlePalu committed Dec 10, 2023
1 parent 73a533c commit 3fa5006
Show file tree
Hide file tree
Showing 67 changed files with 334 additions and 326 deletions.
64 changes: 32 additions & 32 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,8 @@
!test/data/models/gsrpde/*
!test/data/models/gcv/
!test/data/models/gcv/*
!test/data/gcv/sqrpde/
!test/data/gcv/sqrpde/*
!test/data/gcv/qsrpde/
!test/data/gcv/qsrpde/*
!test/data/models/fpca/
!test/data/models/fpca/*
!.github/
Expand Down Expand Up @@ -114,20 +114,20 @@
!test/data/models/gsrpde/2D_test5/*
!test/data/models/gsrpde/2D_test6/
!test/data/models/gsrpde/2D_test6/*
!test/data/models/sqrpde/
!test/data/models/sqrpde/*
!test/data/models/sqrpde/2D_test1/
!test/data/models/sqrpde/2D_test1/*
!test/data/models/sqrpde/2D_test2/
!test/data/models/sqrpde/2D_test2/*
!test/data/models/sqrpde/2D_test3/
!test/data/models/sqrpde/2D_test3/*
!test/data/models/sqrpde/2D_test4/
!test/data/models/sqrpde/2D_test4/*
!test/data/models/sqrpde/2D_test5/
!test/data/models/sqrpde/2D_test5/*
!test/data/models/sqrpde/2D_test6/
!test/data/models/sqrpde/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/
Expand All @@ -144,22 +144,22 @@
!test/data/models/gcv/2D_test7/*
!test/data/models/gcv/2D_test8/
!test/data/models/gcv/2D_test8/*
!test/data/gcv/sqrpde/2D_test1/
!test/data/gcv/sqrpde/2D_test1/*
!test/data/gcv/sqrpde/2D_test2/
!test/data/gcv/sqrpde/2D_test2/*
!test/data/gcv/sqrpde/2D_test3/
!test/data/gcv/sqrpde/2D_test3/*
!test/data/gcv/sqrpde/2D_test4/
!test/data/gcv/sqrpde/2D_test4/*
!test/data/gcv/sqrpde/2D_test5/
!test/data/gcv/sqrpde/2D_test5/*
!test/data/gcv/sqrpde/2D_test6/
!test/data/gcv/sqrpde/2D_test6/*
!test/data/gcv/sqrpde/2D_test7/
!test/data/gcv/sqrpde/2D_test7/*
!test/data/gcv/sqrpde/2D_test8/
!test/data/gcv/sqrpde/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/
Expand Down
6 changes: 3 additions & 3 deletions fdaPDE/models/model_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,9 @@ template <typename Model> class ModelBase {
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_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
Expand Down
45 changes: 24 additions & 21 deletions fdaPDE/models/regression/fpirls.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,55 +37,58 @@ template <typename Model_> class FPIRLS {
private:
using Model = typename std::decay<Model_>::type;
using SolverWrapper = RegressionModel<typename std::decay_t<Model_>::RegularizationType>;
Model& m_;
Model* m_;
// algorithm's parameters
double tolerance_; // treshold on objective functional J to convergence
std::size_t max_iter_; // maximum number of iterations before forced stop
std::size_t k_ = 0; // FPIRLS iteration index
SolverWrapper solver_; // internal solver
public:
// constructor
FPIRLS(Model& m, double tolerance, std::size_t max_iter) :
m_(m), tolerance_(tolerance), max_iter_(max_iter) {
if (!solver_) { // default solver initialization
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
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
solver_ = SolverType(m_.pde(), m_.sampling());
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); }
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_separable<Model_>::value) {
solver_.set_temporal_locations(m_.time_locs());
solver_.set_temporal_locations(m_->time_locs());
}
}
// solver initialization
solver_.set_lambda(m_.lambda());
solver_.set_spatial_locations(m_.locs());
solver_.set_data(m_.data()); // possible covariates are passed from here
solver_.set_spatial_locations(m_->locs());
solver_.set_data(m_->data()); // possible covariates are passed from here
}
};

solver_.set_lambda(m_->lambda());
}
// executes the FPIRLS algorithm
void compute() {
m_.fpirls_init();
// update solver
m_->fpirls_init();
// objective functional value at consecutive iterations
double J_old = tolerance_ + 1;
double J_new = 0;
double J_old = tolerance_ + 1, J_new = 0;
k_ = 0;
while (k_ < max_iter_ && std::abs(J_new - J_old) > tolerance_) {
m_.fpirls_compute_step(); // model specific computation of py_ and pW_
m_->fpirls_compute_step(); // model specific computation of py_ and pW_
// solve weighted least square problem
// \argmin_{\beta, f} [ \norm(W^{1/2}(y - X\beta - f_n))^2 + \lambda \int_D (Lf - u)^2 ]
solver_.data().template insert<double>(OBSERVATIONS_BLK, m_.py());
solver_.data().template insert<double>(WEIGHTS_BLK, m_.pW());
solver_.data().template insert<double>(OBSERVATIONS_BLK, m_->py());
solver_.data().template insert<double>(WEIGHTS_BLK, m_->pW());
// update solver and solve
solver_.init();
solver_.solve();
m_.fpirls_update_step(solver_.fitted(), solver_.beta()); // model specific update step
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());
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;
}
Expand Down
3 changes: 1 addition & 2 deletions fdaPDE/models/regression/gcv.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,7 @@ class GCV {
using This = GCV;
using VectorType = DVector<double>;
using MatrixType = DMatrix<double>;
using ModelType = RegressionView<void>;
ModelType model_; // model to calibrate
RegressionView<void> model_; // model to calibrate
EDFStrategy trS_; // strategy used to evaluate the trace of smoothing matrix S
std::vector<double> edfs_; // equivalent degrees of freedom q + Tr[S]
std::vector<double> gcvs_; // computed values of GCV index
Expand Down
59 changes: 32 additions & 27 deletions fdaPDE/models/regression/gsrpde.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,19 +30,9 @@ namespace models {
// base class for GSRPDE model
template <typename RegularizationType_>
class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, RegularizationType_> {
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

// FPIRLS parameters (set to default)
std::size_t max_iter_ = 200;
double tol_ = 1e-4;
public:
using RegularizationType = RegularizationType_;
using RegularizationType = std::decay_t<RegularizationType_>;
using This = GSRPDE<RegularizationType>;
using Base = RegressionBase<GSRPDE<RegularizationType>, RegularizationType>;
// import commonly defined symbols from base
IMPORT_REGRESSION_SYMBOLS;
Expand All @@ -57,37 +47,40 @@ class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, Regularization
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) {};
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) {};
Base(pde, s, time), distr_(distr) {
fpirls_ = FPIRLS<This>(this, tol_, max_iter_);
};

// 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() { return; };
void init_model() { fpirls_.init(); };
void solve() { // finds a solution to the smoothing problem
// execute FPIRLS for minimization of functional \norm{V^{-1/2}(y - \mu)}^2 + \lambda \int_D (Lf - u)^2
FPIRLS<decltype(*this)> fpirls(*this, tol_, max_iter_); // FPIRLS engine
fpirls.compute();
// fpirls converged: extract matrix W and solution estimates
W_ = fpirls.solver().W();
f_ = fpirls.solver().f();
g_ = fpirls.solver().g();
fpirls_.compute();
// fpirls_ converged: extract matrix W and solution estimates
W_ = fpirls_.solver().W();
f_ = fpirls_.solver().f();
g_ = fpirls_.solver().g();
// store parametric part
if (has_covariates()) {
beta_ = fpirls.solver().beta();
XtWX_ = fpirls.solver().XtWX();
invXtWX_ = fpirls.solver().invXtWX();
U_ = fpirls.solver().U();
V_ = fpirls.solver().V();
beta_ = fpirls_.solver().beta();
XtWX_ = fpirls_.solver().XtWX();
invXtWX_ = fpirls_.solver().invXtWX();
U_ = fpirls_.solver().U();
V_ = fpirls_.solver().V();
}
invA_ = fpirls.solver().invA();
invA_ = fpirls_.solver().invA();
return;
}
// required by FPIRLS (see fpirls.h for details)
Expand Down Expand Up @@ -123,6 +116,18 @@ class GSRPDE : public RegressionBase<GSRPDE<RegularizationType_>, Regularization
}

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

// FPIRLS parameters (set to default)
FPIRLS<This> fpirls_;
std::size_t max_iter_ = 200;
double tol_ = 1e-4;
};

} // namespace models
Expand Down
Loading

0 comments on commit 3fa5006

Please sign in to comment.