diff --git a/src/inputs/InputParser.cpp b/src/inputs/InputParser.cpp
index eec793b062b9bed78a3be49c54c150d9ec98bc84..c24131c6c2e32cc24ca65cd55c7df14ac64350ca 100644
--- a/src/inputs/InputParser.cpp
+++ b/src/inputs/InputParser.cpp
@@ -22,7 +22,8 @@ InputParser::InputParser(pt::ptree IP, std::string fn, std::shared_ptr<MPI_Inter
     _n_models_store(IP.get<int>("n_models_store", _n_residuals)),
     _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))
+    _fix_intercept(IP.get<bool>("fix_intercept", false)),
+    _global_param_opt(IP.get<bool>("global_param_opt", false))
 {
     // Check if param ops are passed without being build with parameterized features
     #ifndef PARAMETERIZE
@@ -32,7 +33,9 @@ InputParser::InputParser(pt::ptree IP, std::string fn, std::shared_ptr<MPI_Inter
     }
     #else
     nlopt::srand(_nlopt_seed);
-    nlopt_wrapper::NLOPT_SEED = _nlopt_seed; 
+    nlopt_wrapper::NLOPT_SEED = _nlopt_seed;
+    nlopt_wrapper::MAX_PARAM_DEPTH = _max_param_depth;
+    nlopt_wrapper::USE_GLOBAL = _global_param_opt;
     #endif
 
     // Read the data file
diff --git a/src/inputs/InputParser.hpp b/src/inputs/InputParser.hpp
index 60975b7863e7a7c3dc5737cb9088ff530ae2313b..b0d2d9baf14849199dd791e9a3bf203af8bab7ba 100644
--- a/src/inputs/InputParser.hpp
+++ b/src/inputs/InputParser.hpp
@@ -72,6 +72,7 @@ public:
     const int _nlopt_seed; //!< The seed used for the nlOpt library
 
     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)
     /**
      * @brief Constructor
      *
diff --git a/src/nl_opt/NLOptWrapper.cpp b/src/nl_opt/NLOptWrapper.cpp
index b3f59f8759c17219552d80d9b1c7b65bc4e0b3b9..3a6f884b7888003887793d4e836c2677f4a1bdbc 100644
--- a/src/nl_opt/NLOptWrapper.cpp
+++ b/src/nl_opt/NLOptWrapper.cpp
@@ -2,6 +2,7 @@
 
 int nlopt_wrapper::MAX_PARAM_DEPTH = -1;
 int nlopt_wrapper::NLOPT_SEED = -1;
+bool nlopt_wrapper::USE_GLOBAL = false;
 
 NLOptimizer::NLOptimizer(
     const std::vector<int>& task_sizes,
@@ -60,10 +61,10 @@ NLOptimizer::NLOptimizer(
     {
         _max_params += std::pow(2, rr);
     }
-    _lb_global.reserve(_max_params);
-    _ub_global.reserve(_max_params);
-    _lb_local.reserve(_max_params);
-    _ub_local.reserve(_max_params);
+    _lb_init.reserve(_max_params);
+    _ub_init.reserve(_max_params);
+    _lb_final.reserve(_max_params);
+    _ub_final.reserve(_max_params);
     _params.reserve(_max_params);
 }
 
@@ -75,7 +76,15 @@ NLOptimizerClassification::NLOptimizerClassification(
     bool reset_max_param_depth,
     const nlopt::algorithm local_opt_alg
 ) :
-    NLOptimizer(task_sizes, prop, n_rung, nlopt_wrapper::objective_class, max_param_depth, reset_max_param_depth, local_opt_alg)
+    NLOptimizer(
+        task_sizes,
+        prop,
+        n_rung,
+        nlopt_wrapper::objective_class,
+        max_param_depth,
+        reset_max_param_depth,
+        local_opt_alg
+    )
 {
     _convex_hull = std::make_shared<ConvexHull1D>(task_sizes, _prop.data());
 }
@@ -90,7 +99,15 @@ NLOptimizerRegression::NLOptimizerRegression(
     bool reset_max_param_depth,
     const nlopt::algorithm local_opt_alg
 ) :
-    NLOptimizer(task_sizes, prop, n_rung, log_reg ? nlopt_wrapper::objective_log_reg : nlopt_wrapper::objective_reg, max_param_depth, reset_max_param_depth, local_opt_alg),
+    NLOptimizer(
+        task_sizes,
+        prop,
+        n_rung,
+        log_reg ? nlopt_wrapper::objective_log_reg : nlopt_wrapper::objective_reg,
+        max_param_depth,
+        reset_max_param_depth,
+        local_opt_alg
+    ),
     _feature_gradient(_max_params * _n_samp, 0.0),
     _residuals(_n_samp, 0.0),
     _cauchy_scaling(cauchy_scaling * cauchy_scaling)
@@ -105,7 +122,16 @@ NLOptimizerLogRegression::NLOptimizerLogRegression(
     bool reset_max_param_depth,
     const nlopt::algorithm local_opt_alg
 ) :
-    NLOptimizerRegression(task_sizes, prop, n_rung, max_param_depth, cauchy_scaling, true, reset_max_param_depth, local_opt_alg)
+    NLOptimizerRegression(
+        task_sizes,
+        prop,
+        n_rung,
+        max_param_depth,
+        cauchy_scaling,
+        true,
+        reset_max_param_depth,
+        local_opt_alg
+    )
 {}
 
 double NLOptimizer::optimize_feature_params(Node* feat)
@@ -122,10 +148,10 @@ double NLOptimizer::optimize_feature_params(Node* feat)
 
     // resize arrays
     _params.resize(n_params);
-    _lb_global.resize(n_params);
-    _ub_global.resize(n_params);
-    _lb_local.resize(n_params);
-    _ub_local.resize(n_params);
+    _lb_init.resize(n_params);
+    _ub_init.resize(n_params);
+    _lb_final.resize(n_params);
+    _ub_final.resize(n_params);
     _free_param.resize(n_params);
 
     // Set params to be 1.0 for all scale parame0.0 for all shift parameters
@@ -134,15 +160,22 @@ double NLOptimizer::optimize_feature_params(Node* feat)
     feat->initialize_params(_params.data() + 2 * _task_sizes.size());
 
     // Set the bounds of the global search to -100 and 100
-    std::fill_n(_lb_global.data(), n_params, -1e2);
-    std::fill_n(_ub_global.data(), n_params,  1e2);
+    std::fill_n(_lb_init.data(), n_params, -1e2);
+    std::fill_n(_ub_init.data(), n_params,  1e2);
 
     // Update bounds based on the feature limits (Some have to be constants)
     feat->set_bounds(
-        _lb_global.data() + 2 * _task_sizes.size(),
-        _ub_global.data() + 2 * _task_sizes.size()
+        _lb_init.data() + 2 * _task_sizes.size(),
+        _ub_init.data() + 2 * _task_sizes.size()
     );
 
+    // Unrestricted local optimization
+    std::copy_n(_lb_init.begin(), n_params, _lb_final.begin());
+    std::copy_n(_ub_init.begin(), n_params, _ub_final.begin());
+
+    std::replace(_lb_final.begin(), _lb_final.end(), -1e2, -1.0 * HUGE_VAL);
+    std::replace(_ub_final.begin(), _ub_final.end(),  1e2, HUGE_VAL);
+
     // Helper variables
     int start = 0;
     double* val_ptr = feat->value_ptr(_params.data() + _task_sizes.size() * 2);
@@ -169,67 +202,66 @@ double NLOptimizer::optimize_feature_params(Node* feat)
     }
 
     std::transform(
-        _lb_global.begin(),
-        _lb_global.end(),
-        _ub_global.begin(),
+        _lb_init.begin(),
+        _lb_init.end(),
+        _ub_init.begin(),
         _free_param.begin(),
         [](double l, double u){return abs(u - l) > 1e-10;}
     );
 
     // Ensure that all values in the parameter vector are inside the bounds
     std::transform(
-        _lb_global.begin(),
-        _lb_global.end(),
+        _lb_init.begin(),
+        _lb_init.end(),
         _params.begin(),
         _params.begin(),
         [](double lb, double p){return p < lb ? lb : p;}
     );
 
     std::transform(
-        _ub_global.begin(),
-        _ub_global.end(),
+        _ub_init.begin(),
+        _ub_init.end(),
         _params.begin(),
         _params.begin(),
         [](double ub, double p){return p > ub ? ub : p;}
     );
 
-    // Set up global optimizer
-    nlopt::opt opt_global(nlopt::GN_ISRES, n_params);
-    opt_global.set_min_objective(_objective, &data);
-    opt_global.set_maxeval(2500);
-    opt_global.set_xtol_rel(1e-2);
-    opt_global.set_lower_bounds(_lb_global);
-    opt_global.set_upper_bounds(_ub_global);
-
     // Set up initial optimizer
     nlopt::opt opt_init(_local_opt_alg, n_params);
     opt_init.set_min_objective(_objective, &data);
     opt_init.set_maxeval(5000);
     opt_init.set_xtol_rel(1e-3);
-    opt_init.set_lower_bounds(_lb_global);
-    opt_init.set_upper_bounds(_ub_global);
+    opt_init.set_lower_bounds(_lb_init);
+    opt_init.set_upper_bounds(_ub_init);
 
     // Set up final optimizer
-    // Unrestricted local optimization
-    std::copy_n(_lb_global.begin(), n_params, _lb_local.begin());
-    std::copy_n(_ub_global.begin(), n_params, _ub_local.begin());
-
-    std::replace(_lb_local.begin(), _lb_local.end(), -1e2, -1.0 * HUGE_VAL);
-    std::replace(_ub_local.begin(), _ub_local.end(),  1e2, HUGE_VAL);
-
     nlopt::opt opt_final(_local_opt_alg, n_params);
     opt_final.set_min_objective(_objective, &data);
     opt_final.set_maxeval(10000);
     opt_final.set_xtol_rel(1e-6);
-    opt_final.set_lower_bounds(_lb_local);
-    opt_final.set_upper_bounds(_ub_local);
+    opt_final.set_lower_bounds(_lb_final);
+    opt_final.set_upper_bounds(_ub_final);
 
     // Try to optimize the parameters, if it fails set minf to a large value
     try
     {
         //Restricted local then global optimization
         nlopt::result result = opt_init.optimize(_params, minf);
-        result = opt_global.optimize(_params, minf);
+
+        if(nlopt_wrapper::USE_GLOBAL)
+        {
+            // Set up global optimizer
+            nlopt::opt opt_global(nlopt::GN_ISRES, n_params);
+            opt_global.set_min_objective(_objective, &data);
+            opt_global.set_maxeval(5000);
+            opt_global.set_xtol_rel(1e-2);
+            opt_global.set_lower_bounds(_lb_init);
+            opt_global.set_upper_bounds(_ub_init);
+
+            // Run the global optimization
+            result = opt_global.optimize(_params, minf);
+        }
+
         nlopt::result result_final = opt_final.optimize(_params, minf);
         if(result_final == nlopt::MAXEVAL_REACHED)
         {
@@ -352,7 +384,8 @@ double nlopt_wrapper::objective_reg(unsigned int n, const double* p, double* gra
                 std::accumulate(
                     d->_optimizer->feature_gradient(pp * n_samp),
                     d->_optimizer->feature_gradient((pp + 1) * n_samp),
-                    0.0
+                    0.0,
+                    [](double total, double res){return total + util_funcs::round2(res, 33);}
                 ),
                 33
             );
@@ -461,15 +494,27 @@ double nlopt_wrapper::objective_log_reg(unsigned int n, const double* p, double*
         );
         for(int pp = 0; pp < 2 * n_task + d->_feat->n_params(); ++ pp)
         {
-            grad[pp] = 1.0 / n_samp * std::accumulate(
-                d->_optimizer->feature_gradient(pp * n_samp),
-                d->_optimizer->feature_gradient((pp + 1) * n_samp),
-                0.0
+            grad[pp] = util_funcs::round2(
+                std::accumulate(
+                    d->_optimizer->feature_gradient(pp * n_samp),
+                    d->_optimizer->feature_gradient((pp + 1) * n_samp),
+                    0.0,
+                    [](double total, double res){return total + util_funcs::round2(res, 33);}
+                ),
+                33
             );
         }
     }
 
-    return std::accumulate(d->_optimizer->residuals(0), d->_optimizer->residuals(n_samp), 0.0);
+    return util_funcs::round2(
+        std::accumulate(
+            d->_optimizer->residuals(0),
+            d->_optimizer->residuals(n_samp),
+            0.0,
+            [](double total, double res){return total + util_funcs::round2(res, 33);}
+        ),
+        33
+    );
 }
 
 
diff --git a/src/nl_opt/NLOptWrapper.hpp b/src/nl_opt/NLOptWrapper.hpp
index ca1ed27a8ada1f5486df309088db8fdd55b71b24..26052ea69fadd05fea354d07c6dbb3f03cc03241 100644
--- a/src/nl_opt/NLOptWrapper.hpp
+++ b/src/nl_opt/NLOptWrapper.hpp
@@ -18,10 +18,10 @@ protected:
 
     std::vector<double> _a; //!< vector to store the A matrix for dgels
     std::vector<double> _params; //!< The working parameters vector
-    std::vector<double> _lb_global; //!< The lower bound for the global optimization
-    std::vector<double> _ub_global; //!< The upper bound for the global optimization
-    std::vector<double> _lb_local; //!< The lower bound for the local optimization
-    std::vector<double> _ub_local; //!< The upper bound for the local optimization
+    std::vector<double> _lb_init; //!< The lower bound for the global optimization
+    std::vector<double> _ub_init; //!< The upper bound for the global optimization
+    std::vector<double> _lb_final; //!< The lower bound for the local optimization
+    std::vector<double> _ub_final; //!< The upper bound for the local optimization
     const std::vector<double> _prop; //!< The property to fit the functions against
     std::vector<double> _prop_copy; //!< Copy of the property to keep for dgels
     std::vector<double> _work; //!< work array for dgels
@@ -35,6 +35,7 @@ protected:
     int _max_param_depth; //!< parameterize features to all depths of the tree
 
     nlopt::algorithm _local_opt_alg; //!< Algorithm used for local optimization
+
 public:
     /**
      * @brief Constructor
@@ -274,6 +275,7 @@ namespace nlopt_wrapper
 {
 extern int MAX_PARAM_DEPTH; //!< The maximum parameter depth for the problem
 extern int NLOPT_SEED; //!< The seed used to initialize the nlopt psuedorandom number seed
+extern bool USE_GLOBAL; //!< If true then use global optimization when optimizing features
 
 typedef struct
 {
@@ -348,6 +350,7 @@ std::shared_ptr<NLOptimizer> get_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param cauchy_scaling scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -356,9 +359,11 @@ inline NLOptimizerRegression get_reg_optimizer(
     py::list prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_list<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
     return NLOptimizerRegression(
@@ -381,6 +386,7 @@ inline NLOptimizerRegression get_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param cauchy_scaling scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -389,9 +395,11 @@ inline NLOptimizerRegression get_reg_optimizer(
     np::ndarray prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_list<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
     return NLOptimizerRegression(
@@ -414,6 +422,7 @@ inline NLOptimizerRegression get_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param cauchy_scaling scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -422,9 +431,11 @@ inline NLOptimizerRegression get_reg_optimizer(
     py::list prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_ndarray<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
     return NLOptimizerRegression(
@@ -447,6 +458,7 @@ inline NLOptimizerRegression get_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param cauchy_scaling scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -455,9 +467,11 @@ inline NLOptimizerRegression get_reg_optimizer(
     np::ndarray prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_ndarray<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
     return NLOptimizerRegression(
@@ -480,6 +494,7 @@ inline NLOptimizerRegression get_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -488,9 +503,11 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
     py::list prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_list<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
     return NLOptimizerLogRegression(
@@ -512,6 +529,7 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -520,9 +538,11 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
     np::ndarray prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_list<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
     return NLOptimizerLogRegression(
@@ -544,6 +564,7 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -552,9 +573,11 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
     py::list prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_ndarray<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
     return NLOptimizerLogRegression(
@@ -576,6 +599,7 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
  * @param scaling factor used for the Cauchy loss function
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
@@ -584,9 +608,11 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
     np::ndarray prop,
     int n_rung,
     int max_param_depth=-1,
-    double cauchy_scaling=0.5
+    double cauchy_scaling=0.5,
+    bool use_global=false
 )
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_ndarray<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
     return NLOptimizerLogRegression(
@@ -607,11 +633,19 @@ inline NLOptimizerLogRegression get_log_reg_optimizer(
  * @param prop The property to fit the functions against
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
-inline NLOptimizerClassification get_class_optimizer(py::list task_sizes, py::list prop, int n_rung, int max_param_depth=-1)
+inline NLOptimizerClassification get_class_optimizer(
+    py::list task_sizes,
+    py::list prop,
+    int n_rung,
+    int max_param_depth=-1,
+    bool use_global=false
+)
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_list<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
     return NLOptimizerClassification(
@@ -631,11 +665,19 @@ inline NLOptimizerClassification get_class_optimizer(py::list task_sizes, py::li
  * @param prop The property to fit the functions against
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
-inline NLOptimizerClassification get_class_optimizer(py::list task_sizes, np::ndarray prop, int n_rung, int max_param_depth=-1)
+inline NLOptimizerClassification get_class_optimizer(
+    py::list task_sizes,
+    np::ndarray prop,
+    int n_rung,
+    int max_param_depth=-1,
+    bool use_global=false
+)
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_list<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
     return NLOptimizerClassification(
@@ -655,11 +697,19 @@ inline NLOptimizerClassification get_class_optimizer(py::list task_sizes, np::nd
  * @param prop The property to fit the functions against
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
-inline NLOptimizerClassification get_class_optimizer(np::ndarray task_sizes, py::list prop, int n_rung, int max_param_depth=-1)
+inline NLOptimizerClassification get_class_optimizer(
+    np::ndarray task_sizes,
+    py::list prop,
+    int n_rung,
+    int max_param_depth=-1,
+    bool use_global=false
+)
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_ndarray<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
     return NLOptimizerClassification(
@@ -679,11 +729,19 @@ inline NLOptimizerClassification get_class_optimizer(np::ndarray task_sizes, py:
  * @param prop The property to fit the functions against
  * @param n_rung Maximum rung of the features
  * @param max_param_depth maximum depth of the binary expression tress to parameterize from the root
+ * @param use_global If true use global optimization
  *
  * @return The correct optimizer
  */
-inline NLOptimizerClassification get_class_optimizer(np::ndarray task_sizes, np::ndarray prop, int n_rung, int max_param_depth=-1)
+inline NLOptimizerClassification get_class_optimizer(
+    np::ndarray task_sizes,
+    np::ndarray prop,
+    int n_rung,
+    int max_param_depth=-1,
+    bool use_global=false
+)
 {
+    USE_GLOBAL = use_global;
     std::vector<int> ts_vec = python_conv_utils::from_ndarray<int>(task_sizes);
     std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
     return NLOptimizerClassification(
diff --git a/src/python/bindings_docstring_keyed.cpp b/src/python/bindings_docstring_keyed.cpp
index 55cf3be330bdce2a2c54b0370db9cb35a39a8299..c127c339d37a928b91a3d2d44d4cb29ab50200c7 100644
--- a/src/python/bindings_docstring_keyed.cpp
+++ b/src/python/bindings_docstring_keyed.cpp
@@ -68,20 +68,20 @@ void sisso::register_all()
         sisso::feature_creation::nloptimizer::registerNLOptimizerRegression();
         sisso::feature_creation::nloptimizer::registerNLOptimizerLogRegression();
 
-        NLOptimizerRegression(*get_reg_optimizer_list_list)(py::list, py::list, int, int, double) = &nlopt_wrapper::get_reg_optimizer;
-        NLOptimizerRegression(*get_reg_optimizer_list_arr)(py::list, np::ndarray, int, int, double) = &nlopt_wrapper::get_reg_optimizer;
-        NLOptimizerRegression(*get_reg_optimizer_arr_list)(np::ndarray, py::list, int, int, double) = &nlopt_wrapper::get_reg_optimizer;
-        NLOptimizerRegression(*get_reg_optimizer_arr_arr)(np::ndarray, np::ndarray, int, int, double) = &nlopt_wrapper::get_reg_optimizer;
-
-        NLOptimizerLogRegression(*get_log_reg_optimizer_list_list)(py::list, py::list, int, int, double) = &nlopt_wrapper::get_log_reg_optimizer;
-        NLOptimizerLogRegression(*get_log_reg_optimizer_list_arr)(py::list, np::ndarray, int, int, double) = &nlopt_wrapper::get_log_reg_optimizer;
-        NLOptimizerLogRegression(*get_log_reg_optimizer_arr_list)(np::ndarray, py::list, int, int, double) = &nlopt_wrapper::get_log_reg_optimizer;
-        NLOptimizerLogRegression(*get_log_reg_optimizer_arr_arr)(np::ndarray, np::ndarray, int, int, double) = &nlopt_wrapper::get_log_reg_optimizer;
-
-        NLOptimizerClassification(*get_class_optimizer_list_list)(py::list, py::list, int, int) = &nlopt_wrapper::get_class_optimizer;
-        NLOptimizerClassification(*get_class_optimizer_list_arr)(py::list, np::ndarray, int, int) = &nlopt_wrapper::get_class_optimizer;
-        NLOptimizerClassification(*get_class_optimizer_arr_list)(np::ndarray, py::list, int, int) = &nlopt_wrapper::get_class_optimizer;
-        NLOptimizerClassification(*get_class_optimizer_arr_arr)(np::ndarray, np::ndarray, int, int) = &nlopt_wrapper::get_class_optimizer;
+        NLOptimizerRegression(*get_reg_optimizer_list_list)(py::list, py::list, int, int, double, bool) = &nlopt_wrapper::get_reg_optimizer;
+        NLOptimizerRegression(*get_reg_optimizer_list_arr)(py::list, np::ndarray, int, int, double, bool) = &nlopt_wrapper::get_reg_optimizer;
+        NLOptimizerRegression(*get_reg_optimizer_arr_list)(np::ndarray, py::list, int, int, double, bool) = &nlopt_wrapper::get_reg_optimizer;
+        NLOptimizerRegression(*get_reg_optimizer_arr_arr)(np::ndarray, np::ndarray, int, int, double, bool) = &nlopt_wrapper::get_reg_optimizer;
+
+        NLOptimizerLogRegression(*get_log_reg_optimizer_list_list)(py::list, py::list, int, int, double, bool) = &nlopt_wrapper::get_log_reg_optimizer;
+        NLOptimizerLogRegression(*get_log_reg_optimizer_list_arr)(py::list, np::ndarray, int, int, double, bool) = &nlopt_wrapper::get_log_reg_optimizer;
+        NLOptimizerLogRegression(*get_log_reg_optimizer_arr_list)(np::ndarray, py::list, int, int, double, bool) = &nlopt_wrapper::get_log_reg_optimizer;
+        NLOptimizerLogRegression(*get_log_reg_optimizer_arr_arr)(np::ndarray, np::ndarray, int, int, double, bool) = &nlopt_wrapper::get_log_reg_optimizer;
+
+        NLOptimizerClassification(*get_class_optimizer_list_list)(py::list, py::list, int, int, bool) = &nlopt_wrapper::get_class_optimizer;
+        NLOptimizerClassification(*get_class_optimizer_list_arr)(py::list, np::ndarray, int, int, bool) = &nlopt_wrapper::get_class_optimizer;
+        NLOptimizerClassification(*get_class_optimizer_arr_list)(np::ndarray, py::list, int, int, bool) = &nlopt_wrapper::get_class_optimizer;
+        NLOptimizerClassification(*get_class_optimizer_arr_arr)(np::ndarray, np::ndarray, int, int, bool) = &nlopt_wrapper::get_class_optimizer;
 
         def("get_reg_optimizer", get_reg_optimizer_list_list, "@DocString_nlopt_wrapper_get_reg_optimizer_list_list");
         def("get_reg_optimizer", get_reg_optimizer_list_arr, "@DocString_nlopt_wrapper_get_reg_optimizer_list_arr");