core submodule update, fixed code to new geometry module. templated c…

…alibrator::GCV, GCV is no more responsible to resize the gcv scalar field, GCV can accept a RegressionView<void> type as model

fdaPDE/calibration/gcv.h

@@ -25,28 +25,34 @@ namespace fdapde {
 namespace calibration {
 
 // GCV calibrator (fulfills the calibration strategy concept)
-class GCV : public CalibratorBase<GCV> {
+template <typename RegularizationType> class GCV : public CalibratorBase<GCV<RegularizationType>> {
    private:
     models::GCV gcv_ {};
     core::Optimizer<models::GCV> opt_ {};
-    DVector<double> optimum_; 
    public:
     // constructor
     template <typename Optimizer_, typename EDFStrategy_> GCV(Optimizer_&& opt, EDFStrategy_&& edf) : opt_(opt) {
+        if constexpr (std::is_same_v<RegularizationType, models::SpaceOnly>) gcv_.resize(1);
+        else gcv_.resize(2);
         gcv_.set_edf_strategy(edf);
     }
     // selects best smoothing parameter of regression model by minimization of GCV index
     template <typename ModelType_, typename... Args>
     DVector<double> fit(ModelType_& model, Args&&... args) {
         gcv_.set_model(model);
-        optimum_ = opt_.optimize(gcv_, std::forward<Args>(args)...);
-        return optimum_;
+        return opt_.optimize(gcv_, std::forward<Args>(args)...);
     }
-    const DVector<double>& optimum() const { return optimum_; }       // optimal tuning parameter
+    const DVector<double>& optimum() const { return opt_.optimum(); }
     const std::vector<double>& edfs() const { return gcv_.edfs(); }   // equivalent degrees of freedom q + Tr[S]
     const std::vector<double>& gcvs() const { return gcv_.gcvs(); }   // computed values of GCV index
     void set_step(double step) { gcv_.set_step(step); }
+    void resize(int gcv_dynamic_inner_size) {   // set GCV's domain dimension
+        static_assert(std::is_same_v<RegularizationType, void>, "THIS_METHOD_IS_FOR_VOID_REGULARIZATION_ONLY");
+        gcv_.resize(gcv_dynamic_inner_size);
+    }
 };
+// template argument deduction rule
+template <typename Optimizer_, typename EDFStrategy_> GCV(Optimizer_&& opt, EDFStrategy_&& edf) -> GCV<void>;
 
 }   // namespace calibration
 }   // namespace fdapde

fdaPDE/models/model_type_erasure.h

@@ -54,6 +54,10 @@ template <typename RegularizationType> struct StatisticalModel__ { };
 
 // basic model interface
 template <> struct StatisticalModel__<void> {
+    // compile time constants
+    using RegularizationType = void;
+    static constexpr int n_lambda = fdapde::Dynamic;
+    // interface implementation
     template <typename M>
     using fn_ptrs = fdapde::mem_fn_ptrs<
       BASE_MODEL_FN_PTRS, static_cast<void (ModelBase<M>::*)(const DVector<double>&)>(&M::set_lambda),

fdaPDE/models/regression/gcv.h

@@ -49,7 +49,7 @@ class GCV {
     using VectorType = DVector<double>;
     using MatrixType = DMatrix<double>;
     static constexpr int DomainDimension = fdapde::Dynamic;
-    RegressionView<void> model_;    // model to calibrate
+    RegressionView<void> model_;
     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
@@ -81,14 +81,15 @@ class GCV {
    public:
     template <typename ModelType_, typename EDFStrategy_>
     GCV(const ModelType_& model, EDFStrategy_&& trS) : model_(model), trS_(trS), gcv_(this, &This::gcv_impl) {
-        // resize gcv input space dimension
-        if constexpr (is_space_only<std::decay_t<ModelType_>>::value) gcv_.resize(1);
-        else gcv_.resize(2);
         // set model pointer in edf computation strategy
         trS_.set_model(model_);
     };
     template <typename ModelType_> GCV(const ModelType_& model) : GCV(model, StochasticEDF()) { }
-    GCV(const GCV& other) : model_(other.model_), trS_(other.trS_), gcv_(this, &This::gcv_impl) {};
+    GCV(const GCV& other) : model_(other.model_), trS_(other.trS_), gcv_(this, &This::gcv_impl) {
+        // copy other GCV functor configuration
+        gcv_.resize(other.gcv_.inner_size());
+	gcv_.set_step(other.gcv_.step());
+    };
     GCV() : gcv_(this, &This::gcv_impl) { }
 
     // call operator and numerical derivative approximations
@@ -112,15 +113,15 @@ class GCV {
 	edfs_.clear(); gcvs_.clear(); cache_.clear();
     }
     template <typename ModelType_> void set_model(ModelType_&& model) {
-      	// resize gcv input space dimension
-        if constexpr (is_space_only<std::decay_t<ModelType_>>::value) gcv_.resize(1);
-        else gcv_.resize(2);
         model_ = model;
 	if(trS_) trS_.set_model(model_);
         edfs_.clear(); gcvs_.clear(); cache_.clear();
     }
     void set_step(double step) { gcv_.set_step(step); }
-
+    void resize(int gcv_dynamic_inner_size) {
+        fdapde_assert(gcv_dynamic_inner_size == 1 || gcv_dynamic_inner_size == 2);
+        gcv_.resize(gcv_dynamic_inner_size);
+    }
     // getters
     const std::vector<double>& edfs() const { return edfs_; }   // equivalent degrees of freedom q + Tr[S]
     const std::vector<double>& gcvs() const { return gcvs_; }   // computed values of GCV index

fdaPDE/models/regression/regression_base.h

@@ -141,7 +141,9 @@ class RegressionBase :
     double ftPf() const { return ftPf(Base::lambda()); }
     // GCV support
     template <typename EDFStrategy_, typename... Args> GCV gcv(Args&&... args) {
-        return GCV(Base::model(), EDFStrategy_(std::forward<Args>(args)...));
+        GCV gcv(Base::model(), EDFStrategy_(std::forward<Args>(args)...));
+        gcv.resize(Base::n_lambda);
+        return gcv;
     }
     const DMatrix<double>& T() {   // T = \Psi^\top*Q*\Psi + P
         T_ = PsiTD() * lmbQ(Psi()) + P();

fdaPDE/models/sampling_design.h

@@ -18,7 +18,7 @@
 #define __SAMPLING_DESIGN_H__
 
 #include <fdaPDE/linear_algebra.h>
-#include <fdaPDE/mesh.h>
+#include <fdaPDE/geometry.h>
 #include <fdaPDE/utils.h>
 #include <fdaPDE/pde.h>
 using fdapde::core::Kronecker;

test/CMakeLists.txt

@@ -1,9 +1,9 @@
 cmake_minimum_required(VERSION 3.14)
 
 project(TEST)
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
-set(CMAKE_CXX_FLAGS "-g -O2 -march=native -std=c++17 -I${CMAKE_CURRENT_SOURCE_DIR}/../ -I${CMAKE_CURRENT_SOURCE_DIR}/../fdaPDE/core/")
+set(CMAKE_CXX_FLAGS "-g -O2 -march=x86-64 -std=c++20 -I${CMAKE_CURRENT_SOURCE_DIR}/../ -I${CMAKE_CURRENT_SOURCE_DIR}/../fdaPDE/core/")
 
 find_package (Eigen3 3.3 REQUIRED NO_MODULE)
 

test/src/gcv_srpde_newton_test.cpp

@@ -32,6 +32,7 @@ using fdapde::models::ExactEDF;
 using fdapde::models::GCV;
 using fdapde::models::StochasticEDF;
 using fdapde::models::Sampling;
+#include "../../fdaPDE/calibration/gcv.h"
 
 #include "utils/constants.h"
 #include "utils/mesh_loader.h"
@@ -175,5 +176,11 @@ TEST(gcv_srpde_newton_test, laplacian_nonparametric_samplingatnodes_newton_fd_st
   auto best_lambda = opt.optimum();
   DVector<double> expected_lambda = SVector<1>(0.0000075627208132);
   // check optimal lambda
-  EXPECT_TRUE( almost_equal(best_lambda[0], expected_lambda[0]) );
+  EXPECT_TRUE(almost_equal(best_lambda[0], expected_lambda[0]));
+
+  // check consistency with GCV calibrator
+  auto GCV_ = fdapde::calibration::GCV<SpaceOnly> {Newton<fdapde::Dynamic>(10, 0.05, 1), StochasticEDF(100, seed)};
+  GCV_.set_step(4e-8);
+  fdapde::models::RegressionView<void> model_view = model;   // get view to model
+  EXPECT_TRUE(GCV_(pt).fit(model_view) == opt.optimum());    // check if it works with a view
 }

test/src/utils/mesh_loader.h

@@ -17,7 +17,7 @@
 #ifndef __MESH_LOADER_H__
 #define __MESH_LOADER_H__
 
-#include <fdaPDE/mesh.h>
+#include <fdaPDE/geometry.h>
 #include <fdaPDE/utils.h>
 #include <gtest/gtest.h>   // testing framework