calibrator call operator returns a ConfiguredCalibrator holding an lv…

…alue reference to itself if called on an lvalue ref-qualified calibrator, holds a copy if the caller is an rvalue. KFoldCV recives the grid of lambda values as DMatrix<double> (as core::Grid), code adapted to this change. calibrator::GCV exposes setters for optimizer and edf computation strategy.

fdaPDE/calibration/calibration_base.h

@@ -24,24 +24,23 @@ namespace calibration {
 
 template <typename Calibrator, typename... Args_> class ConfiguredCalibrator {
    private:
-    using CalibratorType = std::decay_t<Calibrator>;
+    using CalibratorType = std::remove_const_t<Calibrator>;
     using ArgsPackType = std::tuple<std::decay_t<Args_>...>;
     CalibratorType c_;
     ArgsPackType args_;   // parameter pack forwarded to calibration strategy at fit time
     template <typename ModelType_, std::size_t... I>
     auto call_(ModelType_& model, [[maybe_unused]] std::index_sequence<I...> seq) {
         return c_.fit(model, std::get<I>(args_)...);   // forward stored parameters to calibrator
     }
-
     // templated member detection trait for set_model()
     template <typename T, typename M, typename = void> struct callback_require_model : std::false_type { };
     template <typename T, typename M>
     struct callback_require_model<
-      T, M, std::void_t<decltype(std::declval<T>().template set_model<M>(std::declval<M>()))>> : std::true_type { };
+      T, M, std::void_t<decltype(std::declval<T>().template set_model<M>(std::declval<M>()))>
+      > : std::true_type { };
    public:
-    ConfiguredCalibrator() = default;
-    ConfiguredCalibrator(const CalibratorType& c, Args_&&... args) :
-        c_(c), args_(std::make_tuple(std::forward<Args_>(args)...)) {};
+    ConfiguredCalibrator(CalibratorType c, Args_&&... args) :
+        c_(c), args_(std::make_tuple(std::forward<Args_>(args)...)) { }
     template <typename ModelType_> DVector<double> fit(ModelType_&& model) {
         // forward model instance to callbacks which require it
         auto set_model = [&model](auto&& arg) {
@@ -56,18 +55,24 @@ template <typename Calibrator, typename... Args_> class ConfiguredCalibrator {
 };
 
 template <typename T> struct CalibratorBase {
-    template <typename... Args> ConfiguredCalibrator<T, Args...> operator()(Args&&... args) const {
-        return ConfiguredCalibrator<T, Args...>(static_cast<const T&>(*this), std::forward<Args>(args)...);
+    template <typename... Args>
+    ConfiguredCalibrator<std::add_lvalue_reference_t<T>, Args...> operator()(Args&&... args) & {
+        return ConfiguredCalibrator<std::add_lvalue_reference_t<T>, Args...>(
+          static_cast<std::add_lvalue_reference_t<T>>(*this), std::forward<Args>(args)...);
+    }
+    // rvalue ref-qualified overload pushes calibrator by copy
+    template <typename... Args> ConfiguredCalibrator<T, Args...> operator()(Args&&... args) && {
+        return {static_cast<T&>(*this), std::forward<Args>(args)...};
     }
 };
 
 // a type-erasure wrapper for a (configured) calibration algorithm for models of type ModelType
 template <typename ModelType_> struct Calibrator__ {
-    template <typename T> using fn_ptrs = fdapde::mem_fn_ptrs<&T::template fit<ModelType_>, &T::optimum>;
-    decltype(auto) fit(ModelType_& model) { return fdapde::invoke<DVector<double>, 0>(*this, model); }
-    decltype(auto) optimum(ModelType_& model) { return fdapde::invoke<DVector<double>, 1>(*this, model); }
+    template <typename T> using fn_ptrs = mem_fn_ptrs<&T::template fit<ModelType_>, &T::optimum>;
+    decltype(auto) fit(ModelType_& model) { return invoke<DVector<double>, 0>(*this, model); }
+    decltype(auto) optimum() { return invoke<DVector<double>, 1>(*this); }
 };
-template <typename ModelType_> using Calibrator = fdapde::erase<fdapde::heap_storage, Calibrator__<ModelType_>>;
+template <typename ModelType_> using Calibrator = erase<heap_storage, Calibrator__<ModelType_>>;
 
 }   // namespace calibration
 }   // namespace fdapde

fdaPDE/calibration/gcv.h

@@ -34,31 +34,33 @@ template <typename RegularizationType> class GCV : public CalibratorBase<GCV<Reg
     };
    public:
     // constructor
-    GCV() = default;
-    template <typename Optimizer_, typename EDFStrategy_> GCV(Optimizer_&& opt, EDFStrategy_&& edf) : opt_(opt) {
+    GCV() {
         if constexpr (!is_void<RegularizationType>::value) {
-            if constexpr (std::is_same_v<RegularizationType, models::SpaceOnly>) {
-                gcv_.resize(1);
-            } else {   // space-time models
-                gcv_.resize(2);
-            }
+            constexpr int n_lambda = std::is_same_v<RegularizationType, models::SpaceOnly> ? 1 : 2;
+            gcv_.resize(n_lambda);
         }
+    }
+    template <typename Optimizer_, typename EDFStrategy_> GCV(Optimizer_&& opt, EDFStrategy_&& edf) : GCV() {
+        opt_ = opt;
         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) {
-        fdapde_assert(gcv_.inner_size() != 0);
+        fdapde_assert(gcv_.inner_size() != 0 && bool(opt_) == true);
         gcv_.set_model(model);
         return opt_.optimize(gcv_, std::forward<Args>(args)...);
     }
     DVector<double> optimum() { return opt_.optimum(); }              // optimal \lambda found
     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
+    // setters
     void set_step(double step) { gcv_.set_step(step); }
+    template <typename Optimizer_> void set_optimizer(Optimizer_&& opt) { opt_ = opt; }
     template <typename RegularizationType_> void set() {   // set GCV's domain dimension
         fdapde_static_assert(is_void<RegularizationType>::value, THIS_METHOD_IS_FOR_VOID_REGULARIZATION_ONLY);
         gcv_.resize(models::n_smoothing_parameters<RegularizationType_>::value);
     }
+    template <typename EDFStrategy_> void set_edf_strategy(EDFStrategy_&& edf) { gcv_.set_edf_strategy(edf); }
 };
 // template argument deduction rule
 template <typename Optimizer_, typename EDFStrategy_> GCV(Optimizer_&& opt, EDFStrategy_&& edf) -> GCV<void>;

fdaPDE/calibration/kfold_cv.h

@@ -71,24 +71,25 @@ class KCV : public CalibratorBase<KCV> {
 
     // selects best smoothing parameter of model according to a K-fold cross validation strategy
     template <typename ModelType, typename ScoreType>
-    DVector<double> fit(ModelType& model, const std::vector<DVector<double>>& lambdas, ScoreType cv_score) {
+    DVector<double> fit(ModelType& model, const DMatrix<double>& lambdas, ScoreType cv_score) {
+        fdapde_assert(lambdas.cols() == 1 || lambdas.cols() == 2);
         // reserve space for CV scores
-        scores_.resize(K_, lambdas.size());
+        scores_.resize(K_, lambdas.rows());
         if (shuffle_) {   // perform a first shuffling of the data if required
             model.set_data(model.data().shuffle(seed_));
 	}
         // cycle over all tuning parameters
-        for (std::size_t j = 0; j < lambdas.size(); ++j) {
+        for (int j = 0; j < lambdas.rows(); ++j) {
             for (int fold = 0; fold < K_; ++fold) {   // fixed a tuning parameter, cycle over all data splits
                 // compute train-test partition and evaluate CV score
                 TrainTestPartition partition_mask = split(model.data(), fold);
-                scores_.coeffRef(fold, j) = cv_score(lambdas[j], partition_mask.first, partition_mask.second);
+                scores_.coeffRef(fold, j) = cv_score(lambdas.row(j), partition_mask.first, partition_mask.second);
             }
         }
         // reserve space for storing results
-        avg_scores_ = DVector<double>::Zero(lambdas.size());
-        std_scores_ = DVector<double>::Zero(lambdas.size());
-        for (std::size_t j = 0; j < lambdas.size(); ++j) {
+        avg_scores_ = DVector<double>::Zero(lambdas.rows());
+        std_scores_ = DVector<double>::Zero(lambdas.rows());
+        for (int j = 0; j < lambdas.rows(); ++j) {
             // record the average score and its standard deviation computed among the K_ runs
             double avg_score = 0;
             for (int i = 0; i < K_; ++i) avg_score += scores_(i, j);
@@ -104,7 +105,7 @@ class KCV : public CalibratorBase<KCV> {
         // store optimal lambda according to given metric F
         Eigen::Index opt_score;
         avg_scores_.minCoeff(&opt_score);
-        optimum_ = lambdas[opt_score];
+        optimum_ = lambdas.row(opt_score);
 	return optimum_;
     }
 

fdaPDE/models/functional/fpca.h

@@ -47,7 +47,7 @@ class FPCA : public FunctionalBase<FPCA<RegularizationType_>, RegularizationType
         decltype(auto) scores()   const { return invoke<const DMatrix<double>&, 2>(*this); }
     };
     using SolverType = fdapde::erase<heap_storage, SolverType__>;
-    SolverType solver_;
+    SolverType solver_;   // RegularizedSVD solver
    public:
     using RegularizationType = std::decay_t<RegularizationType_>;
     using This = FPCA<RegularizationType>;
@@ -58,13 +58,13 @@ class FPCA : public FunctionalBase<FPCA<RegularizationType_>, RegularizationType
     // constructors
     FPCA() = default;
     // space-only constructor
-    FPCA(const Base::PDE& pde, Sampling s, SolverType solver)
-    requires(is_space_only<This>::value) : Base(pde, s), solver_(solver) { }
+    FPCA(const Base::PDE& pde, Sampling s, SolverType solver) requires(is_space_only<This>::value) :
+      Base(pde, s), solver_(solver) { }
     // space-time separable constructor
     FPCA(const Base::PDE& space_penalty, const Base::PDE& time_penalty, Sampling s, SolverType solver)
-    requires(is_space_time_separable<This>::value) : Base(space_penalty, time_penalty, s), solver_(solver) { }
+      requires(is_space_time_separable<This>::value) : Base(space_penalty, time_penalty, s), solver_(solver) { }
 
-    void init_model() { fdapde_assert(bool(solver_) == true); };   // check solver is assigned
+    void init_model() { fdapde_assert(bool(solver_) == true); };
     void solve() { solver_.compute(X(), *this, n_pc_); }
     // getters
     const DMatrix<double>& loadings() const { return solver_.loadings(); }

fdaPDE/models/functional/fpls.h

@@ -49,11 +49,11 @@ class FPLS : public FunctionalBase<FPLS<RegularizationType_>, RegularizationType
     // constructors
     FPLS() = default;
     // space-only constructor
-    FPLS(const Base::PDE& pde, Sampling s, RegularizedSVD<sequential> rsvd)
-    requires(is_space_only<This>::value) : Base(pde, s), rsvd_(rsvd) { }
+    FPLS(const Base::PDE& pde, Sampling s, RegularizedSVD<sequential> rsvd) requires(is_space_only<This>::value) :
+      Base(pde, s), rsvd_(rsvd) { }
     // space-time separable constructor
     FPLS(const Base::PDE& space_penalty, const Base::PDE& time_penalty, Sampling s, RegularizedSVD<sequential> rsvd)
-    requires(is_space_time_separable<This>::value) : Base(space_penalty, time_penalty, s), rsvd_(rsvd) { }
+      requires(is_space_time_separable<This>::value) : Base(space_penalty, time_penalty, s), rsvd_(rsvd) { }
 
     void init_model() {
         // initialize smoothing solver for regression step

fdaPDE/models/functional/regularized_svd.h

@@ -41,7 +41,7 @@ template <> class RegularizedSVD<sequential> {
    private:
     Calibration calibration_;    // PC function's smoothing parameter selection strategy
     int n_folds_ = 10;   // for a kcv calibration strategy, the number of folds
-    std::vector<DVector<double>> lambda_grid_;
+    DMatrix<double> lambda_grid_;
     // power iteration parameters
     double tolerance_ = 1e-6;   // relative tolerance between Jnew and Jold, used as stopping criterion
     int max_iter_ = 20;         // maximum number of allowed iterations
@@ -163,13 +163,13 @@ template <> class RegularizedSVD<sequential> {
     const DMatrix<double>& loadings() const { return loadings_; }
     const DVector<double>& loadings_norm() const { return loadings_norm_; }
     const std::vector<DVector<double>>& selected_lambdas() const { return selected_lambdas_; }
-    const std::vector<DVector<double>>& lambda_grid() const { return lambda_grid_; }
+    const DMatrix<double>& lambda_grid() const { return lambda_grid_; }
     Calibration calibration() const { return calibration_; }
     // setters
     void set_tolerance(double tolerance) { tolerance_ = tolerance; }
     void set_max_iter(int max_iter) { max_iter_ = max_iter; }
     void set_seed(int seed) { seed_ = seed; }
-    RegularizedSVD& set_lambda(const std::vector<DVector<double>>& lambda_grid) {
+    RegularizedSVD& set_lambda(const DMatrix<double>& lambda_grid) {
         fdapde_assert(calibration_ != Calibration::off);
         lambda_grid_ = lambda_grid;
         return *this;

fdaPDE/models/regression/gcv.h

@@ -39,19 +39,18 @@ class GCV {
     using This = GCV;
     using VectorType = DVector<double>;
     using MatrixType = DMatrix<double>;
-    // type-erased wrapper for Tr[S] computation strategies
+    // type-erased wrapper for Tr[S] computation strategy
     struct EDFStrategy__ {
         template <typename M> using fn_ptrs = fdapde::mem_fn_ptrs<&M::compute, &M::set_model>;
         // forwardings
         decltype(auto) compute() { return fdapde::invoke<double, 0>(*this); }
         void set_model(const RegressionView<void>& model) { fdapde::invoke<void, 1>(*this, model); }
     };
-    using EDFStrategy = fdapde::erase<fdapde::heap_storage, EDFStrategy__>;
 
     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
+    erase<heap_storage, 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
     // cache pairs (lambda, Tr[S]) for fast access if GCV is queried at an already computed point
     std::map<VectorType, double, fdapde::d_vector_compare<double>> cache_;
 
@@ -79,6 +78,7 @@ class GCV {
     }
    public:
     static constexpr int DomainDimension = fdapde::Dynamic;
+    using EDFStrategy = erase<heap_storage, EDFStrategy__>;
     // constructors
     template <typename ModelType_, typename EDFStrategy_>
     GCV(const ModelType_& model, EDFStrategy_&& trS) : model_(model), trS_(trS), gcv_(this, &This::gcv_impl) {

test/src/centering_test.cpp

@@ -62,8 +62,8 @@ TEST(centering_test, srpde_gcv_stochastic_grid) {
     DMatrix<double> u = DMatrix<double>::Zero(domain.mesh.n_elements() * 3, 1);
     PDE<decltype(domain.mesh), decltype(L), DMatrix<double>, FEM, fem_order<1>> pde(domain.mesh, L, u);
     // perform centering
-    std::vector<DVector<double>> lambda_grid;
-    for (double x = -6.0; x <= 0.0; x += 0.5) lambda_grid.push_back(SVector<1>(std::pow(10, x)));
+    DMatrix<double> lambda_grid(12, 1);
+    for (int i = 0; i < 12; ++i) lambda_grid(i, 0) = std::pow(10, -6 + 0.5 * i);
     auto centered_data = center(
       X, SRPDE {pde, Sampling::mesh_nodes},
       fdapde::calibration::GCV<SpaceOnly> {Grid<fdapde::Dynamic> {}, StochasticEDF(100)}(lambda_grid));

test/src/fpca_test.cpp

@@ -126,8 +126,8 @@ TEST(fpca_test, laplacian_samplingatlocations_sequential_gcv) {
     DMatrix<double> u = DMatrix<double>::Zero(domain.mesh.n_elements() * 3, 1);
     PDE<decltype(domain.mesh), decltype(L), DMatrix<double>, FEM, fem_order<1>> pde(domain.mesh, L, u);
     // grid of smoothing parameters
-    std::vector<DVector<double>> lambda_grid;
-    for (double x = -4; x <= -2; x += 0.1) { lambda_grid.push_back(SVector<1>(std::pow(10, x))); }
+    DMatrix<double> lambda_grid(20, 1);
+    for (int i = 0; i < 20; ++i) lambda_grid(i, 0) = std::pow(10, -4 + 0.1 * i);
     // define model
     RegularizedSVD<fdapde::sequential> rsvd(Calibration::gcv);
     rsvd.set_lambda(lambda_grid);
@@ -165,8 +165,8 @@ TEST(fpca_test, laplacian_samplingatlocations_sequential_kcv) {
     DMatrix<double> u = DMatrix<double>::Zero(domain.mesh.n_elements() * 3, 1);
     PDE<decltype(domain.mesh), decltype(L), DMatrix<double>, FEM, fem_order<1>> problem(domain.mesh, L, u);
     // grid of smoothing parameters
-    std::vector<DVector<double>> lambda_grid;
-    for (double x = -4; x <= -2; x += 0.1) lambda_grid.push_back(SVector<1>(std::pow(10, x)));
+    DMatrix<double> lambda_grid(20, 1);
+    for (int i = 0; i < 20; ++i) lambda_grid(i, 0) = std::pow(10, -4 + 0.1 * i);
     // define model
     RegularizedSVD<fdapde::sequential> rsvd(Calibration::kcv);
     rsvd.set_lambda(lambda_grid);

test/src/fpls_test.cpp

@@ -115,8 +115,8 @@ TEST(fpls_test, laplacian_samplingatnodes_sequential_gcv) {
     // define model
     std::size_t seed = 476813;
     // grid for smoothing parameter selection
-    std::vector<DVector<double>> lambda_grid;
-    for (double x = -4; x <= 0; x += 1) lambda_grid.push_back(SVector<1>(std::pow(10, x)));
+    DMatrix<double> lambda_grid(5, 1);
+    for (int i = 0; i < 5; ++i) lambda_grid(i, 0) = std::pow(10, -4 + 1.0 * i);
     RegularizedSVD<fdapde::sequential> rsvd {Calibration::gcv};
     rsvd.set_tolerance(1e-2);
     rsvd.set_max_iter(20);