From 435a3683403d743e90f673acbbc552e5c3d49b21 Mon Sep 17 00:00:00 2001
From: Thomas <purcell@fhi-berlin.mpg.de>
Date: Wed, 2 Jun 2021 17:01:00 +0200
Subject: [PATCH] Reparameterizing based on residual now an option

Both the previous and current behavior is now possible
---
 .../feature_space/FeatureSpace.cpp            | 26 ++++++++++++++-----
 .../feature_space/FeatureSpace.hpp            | 13 +++++++---
 src/inputs/InputParser.cpp                    |  7 +++--
 src/inputs/InputParser.hpp                    |  1 +
 src/python/feature_creation/FeatureSpace.cpp  | 26 +++++++++++++------
 5 files changed, 53 insertions(+), 20 deletions(-)

diff --git a/src/feature_creation/feature_space/FeatureSpace.cpp b/src/feature_creation/feature_space/FeatureSpace.cpp
index 66dee6c8..8d3d0cfc 100644
--- a/src/feature_creation/feature_space/FeatureSpace.cpp
+++ b/src/feature_creation/feature_space/FeatureSpace.cpp
@@ -53,7 +53,8 @@ FeatureSpace::FeatureSpace(
     double cross_corr_max,
     double min_abs_feat_val,
     double max_abs_feat_val,
-    int max_param_depth
+    int max_param_depth,
+    bool reparam_residual
 ):
     _phi(phi_0),
     _phi_0(phi_0),
@@ -76,7 +77,8 @@ FeatureSpace::FeatureSpace(
     _n_feat(phi_0.size()),
     _n_rung_store(max_store_rung),
     _n_rung_generate(n_rung_generate),
-    _max_param_depth(max_param_depth)
+    _max_param_depth(max_param_depth),
+    _reparam_residual(reparam_residual)
 {
     initialize_fs();
 }
@@ -115,7 +117,9 @@ FeatureSpace::FeatureSpace(
     _n_samp(phi_0[0]->n_samp()),
     _n_feat(phi_0.size()),
     _n_rung_store(max_store_rung),
-    _n_rung_generate(n_rung_generate)
+    _n_rung_generate(n_rung_generate),
+    _max_param_depth(0),
+    _reparam_residual(false)
 {
     initialize_fs();
 }
@@ -692,9 +696,17 @@ void FeatureSpace::project_generated(const double* prop, const int size, std::ve
         int index_base = _phi.size() + _n_sis_select * (omp_get_thread_num() + _mpi_comm->size());
 
         #ifdef PARAMETERIZE
-        std::vector<double> prop_vec(size, 0.0);
-        std::copy_n(prop, size, prop_vec.data());
-        std::shared_ptr<NLOptimizer> optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, prop_vec, _max_phi, _max_param_depth);
+        std::shared_ptr<NLOptimizer> optimizer;
+        if(_reparam_residual)
+        {
+            std::vector<double> prop_vec(size, 0.0);
+            std::copy_n(prop, size, prop_vec.data());
+            optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, prop_vec, _max_phi, _max_param_depth);
+        }
+        else
+        {
+            optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, _prop, _max_phi, _max_param_depth);
+        }
         #endif
 
         #pragma omp for schedule(dynamic)
@@ -814,7 +826,7 @@ void FeatureSpace::sis(const std::vector<double>& prop)
 {
     // Reparameterize for the residuals
 #ifdef PARAMETERIZE
-    if(_phi_selected.size() > 0)
+    if(_reparam_residual && (_phi_selected.size() > 0))
     {
         // Make a hard copy of the previously selected features
         for(int ff = _phi_selected.size() - _n_sis_select; ff < _phi_selected.size(); ++ff)
diff --git a/src/feature_creation/feature_space/FeatureSpace.hpp b/src/feature_creation/feature_space/FeatureSpace.hpp
index 6182e6eb..2c0e9902 100644
--- a/src/feature_creation/feature_space/FeatureSpace.hpp
+++ b/src/feature_creation/feature_space/FeatureSpace.hpp
@@ -85,6 +85,7 @@ class FeatureSpace
     const int _n_rung_generate; //!< Total number of rungs to generate on the fly
 
     int _max_param_depth; //!< Max depth to parameterize a feature (default=_max_rung)
+    const bool _reparam_residual; //!< If True then reparameterize using the residuals of each model
 
 public:
 
@@ -108,6 +109,7 @@ public:
      * @param min_abs_feat_val minimum absolute feature value
      * @param max_abs_feat_val maximum absolute feature value
      * @param max_param_depth the maximum paremterization depths for features
+     * @param reparam_residual If True then reparameterize using the residuals of each model
      */
     FeatureSpace(
         std::shared_ptr<MPI_Interface> mpi_comm,
@@ -124,7 +126,8 @@ public:
         double cross_corr_max=1.0,
         double min_abs_feat_val=1e-50,
         double max_abs_feat_val=1e50,
-        int max_param_depth = -1
+        int max_param_depth=-1,
+        bool reparam_residual=false
     );
     #else
     /**
@@ -362,6 +365,7 @@ public:
      * @param min_abs_feat_val minimum absolute feature value
      * @param max_abs_feat_val maximum absolute feature value
      * @param max_param_depth the maximum paremterization depths for features
+     * @param reparam_residual If True then reparameterize using the residuals of each model
      */
     FeatureSpace(
         py::list phi_0,
@@ -377,7 +381,8 @@ public:
         double cross_corr_max=1.0,
         double min_abs_feat_val=1e-50,
         double max_abs_feat_val=1e50,
-        int max_param_depth = -1
+        int max_param_depth = -1,
+        bool reparam_residual=false
     );
 
     /**
@@ -396,6 +401,7 @@ public:
      * @param cross_corr_max Maximum cross-correlation used for selecting features
      * @param min_abs_feat_val minimum absolute feature value
      * @param max_abs_feat_val maximum absolute feature value
+     * @param reparam_residual If True then reparameterize using the residuals of each model
      */
     FeatureSpace(
         py::list phi_0,
@@ -411,7 +417,8 @@ public:
         double cross_corr_max=1.0,
         double min_abs_feat_val=1e-50,
         double max_abs_feat_val=1e50,
-        int max_param_depth = -1
+        int max_param_depth = -1,
+        bool reparam_residual=false
     );
     #else
     /**
diff --git a/src/inputs/InputParser.cpp b/src/inputs/InputParser.cpp
index 4deb53b9..bfe1e331 100644
--- a/src/inputs/InputParser.cpp
+++ b/src/inputs/InputParser.cpp
@@ -23,7 +23,8 @@ InputParser::InputParser(pt::ptree IP, std::string fn, std::shared_ptr<MPI_Inter
     _max_param_depth(IP.get<int>("max_feat_param_depth", _max_rung)),
     _nlopt_seed(IP.get<int>("nlopt_seed", 42)),
     _fix_intercept(IP.get<bool>("fix_intercept", false)),
-    _global_param_opt(IP.get<bool>("global_param_opt", false))
+    _global_param_opt(IP.get<bool>("global_param_opt", false)),
+    _reparam_residual(IP.get<bool>("reparam_residual", false))
 {
     // Check if param ops are passed without being build with parameterized features
     #ifndef PARAMETERIZE
@@ -405,7 +406,9 @@ void InputParser::generate_feature_space(
         _n_rung_generate,
         _cross_cor_max,
         _l_bound,
-        _u_bound
+        _u_bound,
+        _max_param_depth,
+        _reparam_residual
     );
     #else
     _feat_space = std::make_shared<FeatureSpace>(
diff --git a/src/inputs/InputParser.hpp b/src/inputs/InputParser.hpp
index b0d2d9ba..1e82f017 100644
--- a/src/inputs/InputParser.hpp
+++ b/src/inputs/InputParser.hpp
@@ -73,6 +73,7 @@ public:
 
     const bool _fix_intercept; //!< If true force intercept to be 0.0
     const bool _global_param_opt; //!< True if global optimization is requested for non-linear optimization of parameters (Can break reproducibility)
+    const bool _reparam_residual; //!< True if reparameterize features using the residuals
     /**
      * @brief Constructor
      *
diff --git a/src/python/feature_creation/FeatureSpace.cpp b/src/python/feature_creation/FeatureSpace.cpp
index d47d1c1f..70816ce8 100644
--- a/src/python/feature_creation/FeatureSpace.cpp
+++ b/src/python/feature_creation/FeatureSpace.cpp
@@ -15,7 +15,8 @@ FeatureSpace::FeatureSpace(
     double cross_corr_max,
     double min_abs_feat_val,
     double max_abs_feat_val,
-    int max_param_depth
+    int max_param_depth,
+    bool reparam_residual
 ):
     _phi(python_conv_utils::shared_ptr_vec_from_list<Node, FeatureNode>(phi_0)),
     _phi_0(_phi),
@@ -38,7 +39,8 @@ FeatureSpace::FeatureSpace(
     _n_rung_store(max_store_rung),
     _n_rung_generate(n_rung_generate),
     _n_samp(_phi[0]->n_samp()),
-    _max_param_depth(max_param_depth)
+    _max_param_depth(max_param_depth),
+    _reparam_residual(reparam_residual)
 {
     initialize_fs();
 }
@@ -57,7 +59,8 @@ FeatureSpace::FeatureSpace(
     double cross_corr_max,
     double min_abs_feat_val,
     double max_abs_feat_val,
-    int max_param_depth
+    int max_param_depth,
+    bool reparam_residual
 ):
     _phi(python_conv_utils::shared_ptr_vec_from_list<Node, FeatureNode>(phi_0)),
     _phi_0(_phi),
@@ -80,7 +83,8 @@ FeatureSpace::FeatureSpace(
     _n_rung_store(max_store_rung),
     _n_rung_generate(n_rung_generate),
     _n_samp(_phi[0]->n_samp()),
-    _max_param_depth(max_param_depth)
+    _max_param_depth(max_param_depth),
+    _reparam_residual(reparam_residual)
 {
     initialize_fs();
 }
@@ -118,7 +122,9 @@ FeatureSpace::FeatureSpace(
     _n_feat(py::len(phi_0)),
     _n_rung_store(max_store_rung),
     _n_rung_generate(n_rung_generate),
-    _n_samp(_phi[0]->n_samp())
+    _n_samp(_phi[0]->n_samp()),
+    _max_param_depth(0),
+    _reparam_residual(false)
 {
     initialize_fs();
 }
@@ -156,7 +162,9 @@ FeatureSpace::FeatureSpace(
     _n_feat(py::len(phi_0)),
     _n_rung_store(max_store_rung),
     _n_rung_generate(n_rung_generate),
-    _n_samp(_phi[0]->n_samp())
+    _n_samp(_phi[0]->n_samp()),
+    _max_param_depth(0),
+    _reparam_residual(false)
 {
     initialize_fs();
 }
@@ -186,7 +194,8 @@ FeatureSpace::FeatureSpace(
     _n_rung_store(0),
     _n_rung_generate(0),
     _n_samp(_phi_0[0]->n_samp()),
-    _max_param_depth(-1)
+    _max_param_depth(-1),
+    _reparam_residual(false)
 {
     project_funcs::set_project_fxn(project_type, _task_sizes.size(), _project, _project_no_omp);
     comp_feats::set_is_valid_fxn(project_type, _cross_cor_max, _n_samp, _is_valid, _is_valid_feat_list);
@@ -294,7 +303,8 @@ FeatureSpace::FeatureSpace(
     _n_rung_store(0),
     _n_rung_generate(0),
     _n_samp(_phi_0[0]->n_samp()),
-    _max_param_depth(-1)
+    _max_param_depth(-1),
+    _reparam_residual(false)
 {
     project_funcs::set_project_fxn(project_type, _task_sizes.size(), _project, _project_no_omp);
     comp_feats::set_is_valid_fxn(project_type, _cross_cor_max, _n_samp, _is_valid, _is_valid_feat_list);
-- 
GitLab