From 9b06b7ee3bbaf929e38e832cda285032e4c8df63 Mon Sep 17 00:00:00 2001
From: Thomas Purcell <purcell@fhi-berlin.mpg.de>
Date: Wed, 10 Jun 2020 07:27:28 +0200
Subject: [PATCH] Bug fixes for generation on the fly behavior

Division/multiplcation nodes cant have division features/inverse

empyt unit check
---
 src/descriptor_identifier/Model/Model.cpp     |   2 +-
 src/descriptor_identifier/Model/Model.hpp     |   6 +-
 src/descriptor_identifier/SISSORegressor.cpp  |   2 +-
 .../feature_space/FeatureSpace.cpp            | 320 +++++++++++-------
 .../feature_space/FeatureSpace.hpp            |  15 +-
 src/feature_creation/node/FeatureNode.cpp     |   8 +-
 src/feature_creation/node/FeatureNode.hpp     |  13 +-
 src/feature_creation/node/Node.hpp            |   6 +
 .../node/operator_nodes/OperatorNode.hpp      |   8 +-
 .../absolute_difference.cpp                   |   8 +-
 .../absolute_difference.hpp                   |   4 +-
 .../allowed_operator_nodes/absolute_value.cpp |   8 +-
 .../allowed_operator_nodes/absolute_value.hpp |   4 +-
 .../allowed_operator_nodes/add.cpp            |   8 +-
 .../allowed_operator_nodes/add.hpp            |   4 +-
 .../allowed_operator_nodes/cos.cpp            |   8 +-
 .../allowed_operator_nodes/cos.hpp            |   4 +-
 .../allowed_operator_nodes/cube.cpp           |   8 +-
 .../allowed_operator_nodes/cube.hpp           |   4 +-
 .../allowed_operator_nodes/cube_root.cpp      |   8 +-
 .../allowed_operator_nodes/cube_root.hpp      |   4 +-
 .../allowed_operator_nodes/divide.cpp         |  12 +-
 .../allowed_operator_nodes/divide.hpp         |   4 +-
 .../allowed_operator_nodes/exponential.cpp    |   8 +-
 .../allowed_operator_nodes/exponential.hpp    |   4 +-
 .../allowed_operator_nodes/inverse.cpp        |   8 +-
 .../allowed_operator_nodes/inverse.hpp        |   4 +-
 .../allowed_operator_nodes/log.cpp            |   8 +-
 .../allowed_operator_nodes/log.hpp            |   4 +-
 .../allowed_operator_nodes/multiply.cpp       |  14 +-
 .../allowed_operator_nodes/multiply.hpp       |   4 +-
 .../negative_exponential.cpp                  |   8 +-
 .../negative_exponential.hpp                  |   4 +-
 .../allowed_operator_nodes/sin.cpp            |   8 +-
 .../allowed_operator_nodes/sin.hpp            |   4 +-
 .../allowed_operator_nodes/sixth_power.cpp    |   8 +-
 .../allowed_operator_nodes/sixth_power.hpp    |   4 +-
 .../allowed_operator_nodes/square.cpp         |   8 +-
 .../allowed_operator_nodes/square.hpp         |   4 +-
 .../allowed_operator_nodes/square_root.cpp    |   8 +-
 .../allowed_operator_nodes/square_root.hpp    |   4 +-
 .../allowed_operator_nodes/subtract.cpp       |   8 +-
 .../allowed_operator_nodes/subtract.hpp       |   4 +-
 .../node/operator_nodes/allowed_ops.cpp       |  34 +-
 .../node/operator_nodes/allowed_ops.hpp       |   4 +-
 src/feature_creation/units/Unit.cpp           |   4 +-
 src/inputs/InputParser.cpp                    |   4 +-
 src/inputs/InputParser.hpp                    |   1 +
 48 files changed, 378 insertions(+), 263 deletions(-)

diff --git a/src/descriptor_identifier/Model/Model.cpp b/src/descriptor_identifier/Model/Model.cpp
index b83ede88..1085dcea 100644
--- a/src/descriptor_identifier/Model/Model.cpp
+++ b/src/descriptor_identifier/Model/Model.cpp
@@ -1,6 +1,6 @@
 #include <descriptor_identifier/Model/Model.hpp>
 
-Model::Model(std::vector<double> prop_train, std::vector<double> prop_test, std::vector<std::shared_ptr<Node>> feats) :
+Model::Model(std::vector<double> prop_train, std::vector<double> prop_test, std::vector<std::shared_ptr<FeatureNode>> feats) :
     _n_samp_train(feats[0]->n_samp()),
     _n_samp_test(feats[0]->n_test_samp()),
     _n_dim(feats.size() + 1),
diff --git a/src/descriptor_identifier/Model/Model.hpp b/src/descriptor_identifier/Model/Model.hpp
index 574a8abb..261995eb 100644
--- a/src/descriptor_identifier/Model/Model.hpp
+++ b/src/descriptor_identifier/Model/Model.hpp
@@ -7,7 +7,7 @@
 #include<fstream>
 #include<iostream>
 
-#include <feature_creation/node/Node.hpp>
+#include <feature_creation/node/FeatureNode.hpp>
 
 /**
  * @brief Class to store the models found from SISSO
@@ -19,7 +19,7 @@ class Model
     int _n_samp_test; //!< The number of test samples per feature
     int _n_dim; //!< Dimension of the model
 
-    std::vector<std::shared_ptr<Node>> _feats; //!< List of features in the model
+    std::vector<std::shared_ptr<FeatureNode>> _feats; //!< List of features in the model
 
     std::vector<double> _coefs; //!< Coefficients for teh features
     std::vector<double> _prop_train; //!< The property to be modeled
@@ -39,7 +39,7 @@ public:
      * @param prop The property
      * @param feats The features for the model
      */
-    Model(std::vector<double> prop_train, std::vector<double> prop_test, std::vector<std::shared_ptr<Node>> feats);
+    Model(std::vector<double> prop_train, std::vector<double> prop_test, std::vector<std::shared_ptr<FeatureNode>> feats);
 
 
     /**
diff --git a/src/descriptor_identifier/SISSORegressor.cpp b/src/descriptor_identifier/SISSORegressor.cpp
index 1c547dd1..63dff34d 100644
--- a/src/descriptor_identifier/SISSORegressor.cpp
+++ b/src/descriptor_identifier/SISSORegressor.cpp
@@ -189,7 +189,7 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
 
     inds = util_funcs::argsort(all_min_error);
 
-    std::vector<node_ptr> min_nodes(n_dim);
+    std::vector<std::shared_ptr<FeatureNode>> min_nodes(n_dim);
     std::vector<Model> models;
 
     for(int rr = 0; rr < _n_residual; ++rr)
diff --git a/src/feature_creation/feature_space/FeatureSpace.cpp b/src/feature_creation/feature_space/FeatureSpace.cpp
index 38e8e96a..eccc5af0 100644
--- a/src/feature_creation/feature_space/FeatureSpace.cpp
+++ b/src/feature_creation/feature_space/FeatureSpace.cpp
@@ -25,6 +25,7 @@ FeatureSpace::FeatureSpace(
     int max_phi,
     int n_sis_select,
     int max_store_rung,
+    int n_rung_generate,
     double min_abs_feat_val,
     double max_abs_feat_val
 ):
@@ -40,8 +41,14 @@ FeatureSpace::FeatureSpace(
     _n_sis_select(n_sis_select),
     _n_samp(phi_0[0]->n_samp()),
     _n_feat(phi_0.size()),
-    _n_rung_store(max_store_rung)
+    _n_rung_store(max_store_rung),
+    _n_rung_generate(n_rung_generate)
 {
+    if(_n_rung_generate > 1)
+        throw std::logic_error("A maximum of one rung can be generated on the fly.");
+    else if(_max_phi - _n_rung_generate < _n_rung_store)
+        throw std::logic_error("Requesting to store more rungs than what can be pre-generated.");
+
     for(auto & op : allowed_ops)
     {
         if((op.compare("add") == 0) || (op.compare("sub") == 0) || (op.compare("mult") == 0) || (op.compare("abs_diff") == 0))
@@ -56,11 +63,71 @@ FeatureSpace::FeatureSpace(
     _scores.resize(_phi.size());
 }
 
+
+void FeatureSpace::generate_new_feats(std::vector<node_ptr>::iterator& feat, std::vector<node_ptr>& feat_set, int& feat_ind, double l_bound, double u_bound)
+{
+    int phi_ind = feat - _phi.begin();
+    feat_set.reserve(feat_set.size() + _un_operators.size() + phi_ind * (_com_bin_operators.size() + 2 * _bin_operators.size()));
+    for(auto& op : _un_operators)
+    {
+        try
+        {
+            feat_set.push_back(op(*feat, feat_ind, l_bound, u_bound));
+            ++feat_ind;
+        }
+        catch(const InvalidFeatureException& e)
+        {
+            // Do Nothing
+        }
+    }
+
+    for(auto& op : _com_bin_operators)
+    {
+        for(auto feat_2 = _phi.begin(); feat_2 != feat; ++feat_2)
+        {
+            try
+            {
+                feat_set.push_back(op(*feat, *feat_2, feat_ind, l_bound, u_bound));
+                ++feat_ind;
+            }
+            catch(const InvalidFeatureException& e)
+            {
+                // Do Nothing
+            }
+        }
+    }
+
+    for(auto& op : _bin_operators)
+    {
+        for(auto feat_2 = _phi.begin(); feat_2 != feat; ++feat_2)
+        {
+            try
+            {
+                feat_set.push_back(op(*feat, *feat_2, feat_ind, l_bound, u_bound));
+                ++feat_ind;
+            }
+            catch(const InvalidFeatureException& e)
+            {
+                // Do Nothing
+            }
+            try
+            {
+                feat_set.push_back(op(*feat_2, *feat, feat_ind, l_bound, u_bound));
+                ++feat_ind;
+            }
+            catch(const InvalidFeatureException& e)
+            {
+                // Do Nothing
+            }
+        }
+    }
+}
+
 void FeatureSpace::generate_feature_space()
 {
     double u_bound = 1e50;
     double l_bound = 1e-50;
-    for(int nn = 1; nn <= _max_phi; ++nn)
+    for(int nn = 1; nn <= _max_phi - _n_rung_generate; ++nn)
     {
         if(nn == _max_phi)
         {
@@ -72,65 +139,9 @@ void FeatureSpace::generate_feature_space()
 
         int feat_ind = _phi.size();
 
-        // for(auto feat_1 = _phi.begin() + start_end[0]; feat_1 != _phi.begin() + start_end[1]; ++feat_1)
         for(auto feat_1 = _phi.begin() + _mpi_comm->rank(); feat_1 < _phi.end(); feat_1 += _mpi_comm->size())
-        {
-            int phi_ind = feat_1 - _phi.begin();
-            next_phi.reserve(_un_operators.size() + phi_ind * (_com_bin_operators.size() + 2 * _bin_operators.size()));
-            for(auto& op : _un_operators)
-            {
-                try
-                {
-                    next_phi.push_back(op(*feat_1, nn, feat_ind, l_bound, u_bound));
-                    ++feat_ind;
-                }
-                catch(const InvalidFeatureException& e)
-                {
-                    // Do Nothing
-                }
-            }
+            generate_new_feats(feat_1, next_phi, feat_ind, l_bound, u_bound);
 
-            for(auto& op : _com_bin_operators)
-            {
-                for(auto feat_2 = _phi.begin(); feat_2 != feat_1; ++feat_2)
-                {
-                    try
-                    {
-                        next_phi.push_back(op(*feat_1, *feat_2, nn, feat_ind, l_bound, u_bound));
-                        ++feat_ind;
-                    }
-                    catch(const InvalidFeatureException& e)
-                    {
-                        // Do Nothing
-                    }
-                }
-            }
-
-            for(auto& op : _bin_operators)
-            {
-                for(auto feat_2 = _phi.begin(); feat_2 != feat_1; ++feat_2)
-                {
-                    try
-                    {
-                        next_phi.push_back(op(*feat_1, *feat_2, nn, feat_ind, l_bound, u_bound));
-                        ++feat_ind;
-                    }
-                    catch(const InvalidFeatureException& e)
-                    {
-                        // Do Nothing
-                    }
-                    try
-                    {
-                        next_phi.push_back(op(*feat_2, *feat_1, nn, feat_ind, l_bound, u_bound));
-                        ++feat_ind;
-                    }
-                    catch(const InvalidFeatureException& e)
-                    {
-                        // Do Nothing
-                    }
-                }
-            }
-        }
         _mpi_comm->barrier();
         _start_gen.push_back(_phi.size());
         if((nn < _max_phi) || (_mpi_comm->size() == 1))
@@ -294,21 +305,111 @@ void FeatureSpace::project_r(double* prop, int size)
     }
 }
 
+std::vector<double> FeatureSpace::project_r(double* prop, int size, std::vector<node_ptr>& phi)
+{
+    std::vector<double> scores(phi.size(), 0.0);
+    std::vector<double> scores_temp(phi.size(), 0.0);
+
+    for(int ff = 0; ff < phi.size(); ++ff)
+        scores[ff] = -1.0 * std::abs(util_funcs::r(&prop[0], phi[ff]->value_ptr(), _n_samp));
+
+    for(int pp = 1; pp < size / _n_samp; ++pp)
+    {
+        for(int ff = 0; ff < phi.size(); ++ff)
+            scores_temp[ff] = -1.0 * std::abs(util_funcs::r(&prop[_n_samp*pp], phi[ff]->value_ptr(), _n_samp));
+        std::transform(scores_temp.begin(), scores_temp.end(), scores.begin(), scores.begin(), [](double s1, double s2){return std::min(s1, s2);});
+    }
+    return scores;
+}
+
+void FeatureSpace::project_generated(double* prop, int size, std::vector<std::shared_ptr<FeatureNode>>& phi_sel, std::vector<double>& scores_sel, std::vector<double>& scores_comp)
+{
+    for(auto feat = _phi.begin() + _mpi_comm->rank(); feat < _phi.end(); feat += _mpi_comm->size())
+    {
+        int feat_ind = _phi.size();
+
+        std::vector<node_ptr> generated_phi;
+        generate_new_feats(feat, generated_phi, feat_ind, _l_bound, _u_bound);
+
+        std::vector<double> scores = project_r(prop, size, generated_phi);
+        std::vector<int> inds = util_funcs::argsort(scores);
+
+        int ii = 0;
+        while((ii < inds.size()) && (scores[inds[ii]] < *std::max_element(scores_sel.begin(), scores_sel.end())))
+        {
+            double cur_score = scores[inds[ii]];
+            int end_check = std::find_if(scores_sel.begin(), scores_sel.end(), [&cur_score](double s1){return cur_score < s1;}) - scores_sel.begin();
+            if(end_check >= _n_sis_select)
+                break;
+
+            bool is_valid = valid_score_against_current(end_check, generated_phi[inds[ii]]->value_ptr(), scores[inds[ii]], scores_sel, scores_comp);
+            // Check the feature against those selected from previous SIS iterations
+            if((node_value_arrs::N_SELECTED > _n_sis_select) && is_valid)
+                is_valid = valid_score_against_past(generated_phi[inds[ii]]->value_ptr(), scores_comp);
+
+            if(is_valid)
+            {
+                std::shared_ptr<FeatureNode> new_feat = std::make_shared<FeatureNode>(node_value_arrs::N_SELECTED - _n_sis_select + end_check, generated_phi[inds[ii]]->expr(), generated_phi[inds[ii]]->value(), generated_phi[inds[ii]]->test_value(), generated_phi[inds[ii]]->unit(), true);
+                phi_sel.insert(phi_sel.begin() + end_check, new_feat);
+                scores_sel.insert(scores_sel.begin() + end_check, cur_score);
+                for(int jj = end_check + 1; jj < _n_sis_select; ++jj)
+                {
+                    phi_sel[jj]->reindex(node_value_arrs::N_SELECTED - _n_sis_select + jj);
+                    phi_sel[jj]->set_value();
+                }
+            }
+            ++ii;
+        }
+        if(scores_sel.size() > _n_sis_select)
+        {
+            scores_sel.erase(scores_sel.begin() + _n_sis_select, scores_sel.end());
+            phi_sel.erase(phi_sel.begin() + _n_sis_select, phi_sel.end());
+        }
+    }
+}
+
+bool FeatureSpace::valid_score_against_past(double* val_ptr, std::vector<double>& scores_comp)
+{
+    double cur_feat_mean = util_funcs::mean(val_ptr, _n_samp);
+    double cur_feat_std = util_funcs::stand_dev(val_ptr, _n_samp);
+
+    std::transform(val_ptr, val_ptr + _n_samp, val_ptr, [&cur_feat_mean, &cur_feat_std](double val){return (val - cur_feat_mean) / cur_feat_std;});
+
+    dgemv_('T', _n_samp, scores_comp.size(), 1.0 / static_cast<double>(_n_samp), node_value_arrs::D_MATRIX.data(), _n_samp, val_ptr, 1, 0.0, scores_comp.data(), 1);
+
+    if(1.0 - util_funcs::max_abs_val<double>(scores_comp.data(), scores_comp.size()) < 1e-13)
+        return false;
+
+    std::transform(val_ptr, val_ptr + _n_samp, val_ptr, [&cur_feat_mean, &cur_feat_std](double val){return val * cur_feat_std + cur_feat_mean;});
+    return true;
+}
+
+bool FeatureSpace::valid_score_against_current(int end_check, double* val_ptr, double cur_score, std::vector<double>& scores_sel, std::vector<double>& scores_comp)
+{
+    std::transform(scores_sel.begin(), scores_sel.begin() + end_check, scores_comp.begin(), [&cur_score](double score){return cur_score - score;});
+
+    // If two scores are the same then they are possibly the same feature, if not then they can't be
+    if(*std::min_element(scores_comp.begin(), scores_comp.begin() + end_check) < 1e-10)
+    {
+        int dd = std::min_element(scores_comp.begin(), scores_comp.begin() + end_check) - scores_comp.begin();
+        if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(node_value_arrs::N_SELECTED - _n_sis_select + dd), val_ptr, _n_samp)) < 1e-13)
+            return false;
+    }
+    return true;
+}
+
 void FeatureSpace::sis(std::vector<double>& prop)
 {
     std::vector<double> means(node_value_arrs::N_SELECTED);
     std::vector<double> stand_devs(node_value_arrs::N_SELECTED);
     std::vector<double> scores_comp(std::max(node_value_arrs::N_SELECTED, _n_sis_select), 1.0);
 
-    std::vector<double> scores_selected(_n_sis_select, 0.0);
-    std::vector<int> inds_selected(_n_sis_select, -1);
-    std::vector<node_ptr> phi_selected;
-    phi_selected.reserve(_n_sis_select);
+    std::vector<double> scores_sel(_n_sis_select, 0.0);
+    std::vector<std::shared_ptr<FeatureNode>> phi_sel;
+    phi_sel.reserve(_n_sis_select);
 
     int cur_feat = node_value_arrs::N_SELECTED;
 
-    double cur_feat_mean = 0.0;
-    double cur_feat_std = 0.0;
 
     // Standardize the description matrix
     if(cur_feat > 0)
@@ -332,51 +433,38 @@ void FeatureSpace::sis(std::vector<double>& prop)
 
     while((cur_feat_local != _n_sis_select) && (ii < _scores.size()))
     {
-        bool is_valid = true;
-
+        bool is_valid = valid_score_against_current(cur_feat_local, _phi[inds[ii]]->value_ptr(), _scores[inds[ii]], scores_sel, scores_comp);
         // Check the feature against those selected from previous SIS iterations
-        if(cur_feat > 0)
-        {
-            cur_feat_mean = util_funcs::mean(_phi[inds[ii]]->value_ptr(), _n_samp);
-            cur_feat_std = util_funcs::stand_dev(_phi[inds[ii]]->value_ptr(), _n_samp);
-            std::transform(_phi[inds[ii]]->value_ptr(), _phi[inds[ii]]->value_ptr() + _n_samp, _phi[inds[ii]]->value_ptr(), [&cur_feat_mean, &cur_feat_std](double val){return (val - cur_feat_mean) / cur_feat_std;});
-
-            dgemv_('T', _n_samp, cur_feat, 1.0 / static_cast<double>(_n_samp), node_value_arrs::D_MATRIX.data(), _n_samp, _phi[inds[ii]]->value_ptr(), 1, 0.0, scores_comp.data(), 1);
-
-            if(1.0 - util_funcs::max_abs_val<double>(scores_comp.data(), scores_comp.size()) < 1e-13)
-                is_valid = false;
+        if(cur_feat > 0 && is_valid)
+            is_valid = valid_score_against_past(_phi[inds[ii]]->value_ptr(), scores_comp);
 
-            std::transform(_phi[inds[ii]]->value_ptr(), _phi[inds[ii]]->value_ptr() + _n_samp, _phi[inds[ii]]->value_ptr(), [&cur_feat_mean, &cur_feat_std](double val){return val * cur_feat_std + cur_feat_mean;});
-        }
-
-        // Check the feature against all those selected previously in this SIS iteration
-        cur_score =  _scores[inds[ii]];
-        std::transform(scores_selected.begin(), scores_selected.begin() + cur_feat_local, scores_comp.begin(), [&cur_score](double score){return cur_score - score;});
-        // If two scores are the same then they are possibly the same feature, if not then they can't be
-        if(*std::min_element(scores_comp.begin(), scores_comp.begin() + cur_feat_local) < 1e-10)
-        {
-            int dd = std::min_element(scores_comp.begin(), scores_comp.begin() + cur_feat_local) - scores_comp.begin();
-            if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(cur_feat + dd), _phi[inds[ii]]->value_ptr(), _n_samp)) < 1e-13)
-                is_valid = false;
-        }
         if(is_valid)
         {
-            inds_selected[cur_feat_local] = _phi[inds[ii]]->feat_ind();
-            scores_selected[cur_feat_local] = _scores[inds[ii]];
-            phi_selected.push_back(std::make_shared<FeatureNode>(cur_feat + cur_feat_local, _phi[inds[ii]]->expr(), _phi[inds[ii]]->value(), _phi[inds[ii]]->test_value(), _phi[inds[ii]]->unit(), true));
+            std::cout << _phi[inds[ii]]->expr() << _scores[inds[ii]] << std::endl;
+            scores_sel[cur_feat_local] = _scores[inds[ii]];
+            phi_sel.push_back(std::make_shared<FeatureNode>(cur_feat + cur_feat_local, _phi[inds[ii]]->expr(), _phi[inds[ii]]->value(), _phi[inds[ii]]->test_value(), _phi[inds[ii]]->unit(), true));
             ++cur_feat_local;
         }
         ++ii;
     }
 
+    if(_n_rung_generate > 0)
+    {
+        for(auto& feat : phi_sel)
+        phi_sel.resize(cur_feat_local);
+        scores_sel.resize(cur_feat_local);
+        project_generated(prop.data(), prop.size(), phi_sel, scores_sel, scores_comp);
+    }
+
     // Unstandardize the description matrix
     if(cur_feat > 0)
         for(int dd = 0; dd < cur_feat; ++dd)
             std::transform(node_value_arrs::get_d_matrix_ptr(dd), node_value_arrs::get_d_matrix_ptr(dd) + _n_samp, node_value_arrs::get_d_matrix_ptr(dd), [&means, &stand_devs, &dd](double val){return val * stand_devs[dd] + means[dd];});
 
-    phi_selected.resize(_n_sis_select);
+    phi_sel.resize(_n_sis_select);
+    scores_sel.resize(_n_sis_select);
 
-    // If we are only on one process then phi_selected are the selected features
+    // If we are only on one process then phi_sel are the selected features
     if(_mpi_comm->size() > 1)
     {
         // Prepare to get all scores, features, and indicies
@@ -386,7 +474,7 @@ void FeatureSpace::sis(std::vector<double>& prop)
         int unique_scores = 1;
         if(_mpi_comm->rank() == 0)
         {
-            std::copy_n(scores_selected.data(), _n_sis_select, sent_scores.data());
+            std::copy_n(scores_sel.data(), _n_sis_select, sent_scores.data());
             for(int rr = 1; rr < _mpi_comm->size(); ++rr)
                 _mpi_comm->recv(rr, _mpi_comm->cantorTagGen(rr, 0, 1, 0), &sent_scores[rr * _n_sis_select], _n_sis_select);
 
@@ -405,7 +493,7 @@ void FeatureSpace::sis(std::vector<double>& prop)
         }
         else
         {
-            _mpi_comm->send(0, _mpi_comm->cantorTagGen(_mpi_comm->rank(), 0, 1, 0), scores_selected.data(), _n_sis_select);
+            _mpi_comm->send(0, _mpi_comm->cantorTagGen(_mpi_comm->rank(), 0, 1, 0), scores_sel.data(), _n_sis_select);
             inds = {};
         }
         mpi::broadcast(*_mpi_comm, inds, 0);
@@ -416,21 +504,20 @@ void FeatureSpace::sis(std::vector<double>& prop)
             int compare_ind = ii + _mpi_comm->rank() * _n_sis_select;
             if(std::none_of(inds.begin(), inds.end(), [&compare_ind](int i1){return i1 == compare_ind;}))
             {
-                scores_selected.erase(scores_selected.begin() + ii);
-                phi_selected.erase(phi_selected.begin() + ii);
-                inds_selected.erase(inds_selected.begin() + ii);
+                scores_sel.erase(scores_sel.begin() + ii);
+                phi_sel.erase(phi_sel.begin() + ii);
             }
         }
-        scores_selected.resize(_n_sis_select, 0.0);
-        phi_selected.resize(_n_sis_select, nullptr);
+        scores_sel.resize(_n_sis_select, 0.0);
+        phi_sel.resize(_n_sis_select, nullptr);
 
         if(_mpi_comm->rank() == 0)
         {
             std::vector<double> sent_scores(_n_sis_select * _mpi_comm->size(), 0.0);
-            std::vector<node_ptr> sent_phi(_n_sis_select * _mpi_comm->size());
+            std::vector<std::shared_ptr<FeatureNode>> sent_phi(_n_sis_select * _mpi_comm->size());
 
-            std::copy_n(scores_selected.begin(), _n_sis_select, sent_scores.begin());
-            std::copy_n(phi_selected.begin(), _n_sis_select, sent_phi.begin());
+            std::copy_n(scores_sel.begin(), _n_sis_select, sent_scores.begin());
+            std::copy_n(phi_sel.begin(), _n_sis_select, sent_phi.begin());
 
             for(int rr = 1; rr < _mpi_comm->size(); ++rr)
             {
@@ -452,31 +539,22 @@ void FeatureSpace::sis(std::vector<double>& prop)
             else
             {
                 inds = util_funcs::argsort(sent_scores);
-                // Clear out the previous D matrix values (from creation of each process' phi_selected)
+                // Clear out the previous D matrix values (from creation of each process' phi_sel)
                 std::fill_n(node_value_arrs::get_d_matrix_ptr(cur_feat), _n_sis_select * node_value_arrs::N_SAMPLES, 0.0);
 
                 cur_feat_local = 0;
                 int ii = 0;
                 std::fill_n(scores_comp.begin(), _n_sis_select, 1.0);
-                scores_selected = std::vector<double>(_n_sis_select, 0.0);
+                scores_sel = std::vector<double>(_n_sis_select, 0.0);
 
                 // Get the n_sis_select best features (compare against features sent from other processes)
                 while((cur_feat != node_value_arrs::N_SELECTED) && (ii < sent_scores.size()))
                 {
-                    bool is_valid = true;
-                    cur_score =  sent_scores[inds[ii]];
-                    std::transform(scores_selected.begin(), scores_selected.begin() + cur_feat_local, scores_comp.begin(), [&cur_score](double score){return cur_score - score;});
-                    if(*std::min_element(scores_comp.begin(), scores_comp.begin() + cur_feat_local) < 1e-10)
-                    {
-                        int dd = std::min_element(scores_comp.begin(), scores_comp.begin() + cur_feat_local) - scores_comp.begin();
-                        if(1.0 - std::abs(util_funcs::r(node_value_arrs::get_d_matrix_ptr(cur_feat + dd), sent_phi[inds[ii]]->value().data(), _n_samp)) < 1e-13)
-                            is_valid = false;
-                    }
-
-                    if(is_valid)
+                    if(valid_score_against_current(cur_feat_local, sent_phi[inds[ii]]->value().data(), sent_scores[inds[ii]], scores_sel, scores_comp))
                     {
                         _phi_selected.push_back(sent_phi[inds[ii]]);
                         _phi_selected.back()->reindex(cur_feat);
+                        scores_sel[cur_feat_local] = sent_scores[inds[ii]];
                         ++cur_feat_local;
                         ++cur_feat;
                     }
@@ -489,8 +567,8 @@ void FeatureSpace::sis(std::vector<double>& prop)
         }
         else
         {
-            _mpi_comm->send(0, _mpi_comm->cantorTagGen(_mpi_comm->rank(), 0, 2, 0), scores_selected.data(), _n_sis_select);
-            _mpi_comm->send(0, _mpi_comm->cantorTagGen(_mpi_comm->rank(), 0, 2, 1), phi_selected.data(), _n_sis_select);
+            _mpi_comm->send(0, _mpi_comm->cantorTagGen(_mpi_comm->rank(), 0, 2, 0), scores_sel.data(), _n_sis_select);
+            _mpi_comm->send(0, _mpi_comm->cantorTagGen(_mpi_comm->rank(), 0, 2, 1), phi_sel.data(), _n_sis_select);
             _phi_selected.resize(node_value_arrs::N_SELECTED);
         }
         mpi::broadcast(*_mpi_comm, &_phi_selected[_phi_selected.size() - _n_sis_select], _n_sis_select, 0);
@@ -504,7 +582,7 @@ void FeatureSpace::sis(std::vector<double>& prop)
     else
     {
         // cur_feat += cur_feat_local;
-        for(auto& feat : phi_selected)
+        for(auto& feat : phi_sel)
         {
             _phi_selected.push_back(feat);
             _phi_selected.back()->reindex(cur_feat);
diff --git a/src/feature_creation/feature_space/FeatureSpace.hpp b/src/feature_creation/feature_space/FeatureSpace.hpp
index de7f41a9..4cf22bc2 100644
--- a/src/feature_creation/feature_space/FeatureSpace.hpp
+++ b/src/feature_creation/feature_space/FeatureSpace.hpp
@@ -19,8 +19,8 @@
  */
 class FeatureSpace
 {
+    std::vector<std::shared_ptr<FeatureNode>> _phi_selected; //!< selected features
     std::vector<node_ptr> _phi; //!< all features
-    std::vector<node_ptr> _phi_selected; //!< selected features
     std::vector<node_ptr> _phi_0; //!< initial feature space
 
     std::vector<std::string> _allowed_ops; //!< list of all allowed operators strings
@@ -43,6 +43,7 @@ class FeatureSpace
     int _n_samp; //!< Number of samples
     int _n_feat; //!< Total number of features
     int _n_rung_store; //!< Total rungs stored
+    int _n_rung_generate; //!< Total number of rungs to generate on the fly
 
 public:
     /**
@@ -62,6 +63,7 @@ public:
         int max_phi=1,
         int n_sis_select=1,
         int max_store_rung=2,
+        int n_rung_generate=0,
         double min_abs_feat_val=1e-50,
         double max_abs_feat_val=1e50
     );
@@ -75,7 +77,7 @@ public:
     /**
      * @brief Accessor function for _phi_selected
      */
-    inline std::vector<node_ptr> phi_selected(){return _phi_selected;};
+    inline std::vector<std::shared_ptr<FeatureNode>> phi_selected(){return _phi_selected;};
 
     /**
      * @brief Accessor function for _phi
@@ -105,6 +107,15 @@ public:
      */
     void project_r(double* prop, int size);
 
+    std::vector<double> project_r(double* prop, int size, std::vector<node_ptr>& phi);
+
+    void generate_new_feats(std::vector<node_ptr>::iterator& feat, std::vector<node_ptr>& feat_set, int& feat_ind, double l_bound=1e-50, double u_bound=1e50);
+
+    void project_generated(double* prop, int size, std::vector<std::shared_ptr<FeatureNode>>& phi_selected, std::vector<double>& scores_selected, std::vector<double>& scores_comp);
+
+    bool valid_score_against_past(double* val_ptr, std::vector<double>& scores_comp);
+
+    bool valid_score_against_current(int end_check, double* val_ptr, double cur_score, std::vector<double>& scores_selected, std::vector<double>& scores_comp);
     /**
      * @brief Perform SIS on a feature set with a specified property
      * @details Perform sure-independence screening with either the correct property
diff --git a/src/feature_creation/node/FeatureNode.cpp b/src/feature_creation/node/FeatureNode.cpp
index 7d8ac1f3..a808977b 100644
--- a/src/feature_creation/node/FeatureNode.cpp
+++ b/src/feature_creation/node/FeatureNode.cpp
@@ -5,11 +5,11 @@ FeatureNode::FeatureNode()
 
 FeatureNode::FeatureNode(int feat_ind, std::string expr, std::vector<double> value, std::vector<double> test_value, Unit unit, bool selected) :
     Node(feat_ind, value.size(), test_value.size()),
-    _selected(selected),
-    _expr(expr),
-    _unit(unit),
     _value(value),
-    _test_value(test_value)
+    _test_value(test_value),
+    _unit(unit),
+    _expr(expr),
+    _selected(selected)
 {
     set_value();
     set_test_value();
diff --git a/src/feature_creation/node/FeatureNode.hpp b/src/feature_creation/node/FeatureNode.hpp
index 284cb926..ff643fc2 100644
--- a/src/feature_creation/node/FeatureNode.hpp
+++ b/src/feature_creation/node/FeatureNode.hpp
@@ -34,11 +34,11 @@ class FeatureNode: public Node
     }
 
 protected:
-    bool _selected; //!< True if the features was selected
-    std::string _expr; //!< Expression of the feature
-    Unit _unit; //!< Unit for the feature
     std::vector<double> _value; //!< values for the feature
     std::vector<double> _test_value; //!< test values for the feature
+    Unit _unit; //!< Unit for the feature
+    std::string _expr; //!< Expression of the feature
+    bool _selected; //!< True if the features was selected
 public:
     /**
      * @brief Base Constructor
@@ -57,12 +57,19 @@ public:
      */
     FeatureNode(int feat_ind, std::string expr, std::vector<double> value, std::vector<double> test_value, Unit unit, bool selected=false);
 
+    FeatureNode(const FeatureNode&) = default;
+    FeatureNode(FeatureNode&&) = default;
+
+    FeatureNode& operator=(const FeatureNode&) = default;
+    FeatureNode& operator=(FeatureNode&&) = default;
+
     ~FeatureNode();
 
     /**
      * @brief Get the expression for the overall descriptor (From head node down)
      */
     inline std::string expr(){return _expr;}
+    inline std::string expr()const{return _expr;}
 
     /**
      * @brief Get the unit for the overall descriptor (From head node down)
diff --git a/src/feature_creation/node/Node.hpp b/src/feature_creation/node/Node.hpp
index a4602172..3ee28fe9 100644
--- a/src/feature_creation/node/Node.hpp
+++ b/src/feature_creation/node/Node.hpp
@@ -59,6 +59,12 @@ public:
      */
     Node(int feat_ind, int n_samp, int n_test_samp);
 
+    Node(const Node& o) = default;
+    Node(Node&& o) = default;
+
+    Node& operator= (const Node& o) = default;
+    Node& operator= (Node&& o) = default;
+
     virtual ~Node();
 
     /**
diff --git a/src/feature_creation/node/operator_nodes/OperatorNode.hpp b/src/feature_creation/node/operator_nodes/OperatorNode.hpp
index 1d1ee423..7e3b252b 100644
--- a/src/feature_creation/node/operator_nodes/OperatorNode.hpp
+++ b/src/feature_creation/node/operator_nodes/OperatorNode.hpp
@@ -48,11 +48,17 @@ public:
      * @param rung run the feature is on (depth of the tree)
      * @param feat_ind index of the feature
      */
-    OperatorNode(std::array<node_ptr, N> feats, int rung, int feat_ind) :
+    OperatorNode(std::array<node_ptr, N> feats, int feat_ind) :
         Node(feat_ind, feats[0]->n_samp(), feats[0]->n_test_samp()),
         _feats(feats)
     {}
 
+    OperatorNode(const OperatorNode& o) = default;
+    OperatorNode(OperatorNode&& o) = default;
+
+    OperatorNode& operator=(const OperatorNode& o) = default;
+    OperatorNode& operator=(OperatorNode&& o) = default;
+
     virtual ~OperatorNode()
     {}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.cpp
index 4feb01f2..bf5aa248 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.cpp
@@ -5,8 +5,8 @@
 AbsDiffNode::AbsDiffNode()
 {}
 
-AbsDiffNode::AbsDiffNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+AbsDiffNode::AbsDiffNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != feats[1]->unit())
         throw InvalidFeatureException();
@@ -33,8 +33,8 @@ AbsDiffNode::AbsDiffNode(std::array<node_ptr, 2> feats, int rung, int feat_ind,
     set_test_value();
  }
 
-AbsDiffNode::AbsDiffNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat_1, feat_2}, rung, feat_ind)
+AbsDiffNode::AbsDiffNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat_1, feat_2}, feat_ind)
 {
     if(feat_1->unit() != feat_2->unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.hpp
index 0aca9016..2179a7c4 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_difference.hpp
@@ -16,9 +16,9 @@ class AbsDiffNode: public OperatorNode<2>
 public:
     AbsDiffNode();
 
-    AbsDiffNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    AbsDiffNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound);
 
-    AbsDiffNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound);
+    AbsDiffNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.cpp
index 670f4592..90ee52bf 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.cpp
@@ -5,8 +5,8 @@
 AbsNode::AbsNode()
 {}
 
-AbsNode::AbsNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+AbsNode::AbsNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     set_value();
     if(is_nan() || is_const() || (util_funcs::max_abs_val<double>(value_ptr(), _n_samp) > u_bound) || (util_funcs::max_abs_val<double>(value_ptr(), _n_samp) < l_bound))
@@ -15,8 +15,8 @@ AbsNode::AbsNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l
     set_test_value();
 }
 
-AbsNode::AbsNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+AbsNode::AbsNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     set_value();
     if(is_nan() || is_const() || (util_funcs::max_abs_val<double>(value_ptr(), _n_samp) > u_bound) || (util_funcs::max_abs_val<double>(value_ptr(), _n_samp) < l_bound))
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.hpp
index c42265ca..a7a111f9 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/absolute_value.hpp
@@ -15,9 +15,9 @@ class AbsNode: public OperatorNode<1>
 public:
     AbsNode();
 
-    AbsNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    AbsNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    AbsNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    AbsNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.cpp
index b764c054..4786123b 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.cpp
@@ -3,8 +3,8 @@
 AddNode::AddNode()
 {}
 
-AddNode::AddNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+AddNode::AddNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != feats[1]->unit())
         throw InvalidFeatureException();
@@ -31,8 +31,8 @@ AddNode::AddNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l
     set_test_value();
  }
 
-AddNode::AddNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat_1, feat_2}, rung, feat_ind)
+AddNode::AddNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat_1, feat_2}, feat_ind)
 {
     if(feat_1->unit() != feat_2->unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.hpp
index b2579b81..d722fc8e 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/add.hpp
@@ -15,9 +15,9 @@ class AddNode: public OperatorNode<2>
 public:
     AddNode();
 
-    AddNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    AddNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound);
 
-    AddNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound);
+    AddNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.cpp
index 2744e853..9bf1fab6 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.cpp
@@ -3,8 +3,8 @@
 CosNode::CosNode()
 {}
 
-CosNode::CosNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+CosNode::CosNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != Unit())
         throw InvalidFeatureException();
@@ -19,8 +19,8 @@ CosNode::CosNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l
        set_test_value();
  }
 
-CosNode::CosNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+CosNode::CosNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->unit() != Unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.hpp
index c006a46b..66799f4e 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cos.hpp
@@ -15,9 +15,9 @@ class CosNode: public OperatorNode<1>
 public:
     CosNode();
 
-    CosNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    CosNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    CosNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    CosNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return Unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.cpp
index f9b47cf7..5ad863e1 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.cpp
@@ -3,8 +3,8 @@
 CbNode::CbNode()
 {}
 
-CbNode::CbNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+CbNode::CbNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->type() == NODE_TYPE::CBRT)
         throw InvalidFeatureException();
@@ -16,8 +16,8 @@ CbNode::CbNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_b
     set_test_value();
 }
 
-CbNode::CbNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+CbNode::CbNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->type() == NODE_TYPE::CBRT)
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.hpp
index 3d7eac11..b110dc91 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube.hpp
@@ -15,9 +15,9 @@ class CbNode: public OperatorNode<1>
 public:
     CbNode();
 
-    CbNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    CbNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    CbNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    CbNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit()^(3.0);}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.cpp
index a679933c..2cc44918 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.cpp
@@ -3,8 +3,8 @@
 CbrtNode::CbrtNode()
 {}
 
-CbrtNode::CbrtNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+CbrtNode::CbrtNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->type() == NODE_TYPE::CB)
         throw InvalidFeatureException();
@@ -16,8 +16,8 @@ CbrtNode::CbrtNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double
     set_test_value();
 }
 
-CbrtNode::CbrtNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+CbrtNode::CbrtNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->type() == NODE_TYPE::CB)
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.hpp
index 1509726f..faff0108 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/cube_root.hpp
@@ -15,9 +15,9 @@ class CbrtNode: public OperatorNode<1>
 public:
     CbrtNode();
 
-    CbrtNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    CbrtNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    CbrtNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    CbrtNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit()^(1.0 / 3.0);}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.cpp
index b0d4dcfc..f5d65667 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.cpp
@@ -3,10 +3,10 @@
 DivNode::DivNode()
 {}
 
-DivNode::DivNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+DivNode::DivNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
-    if((feats[0]->type() == NODE_TYPE::INV) || (feats[1]->type() == NODE_TYPE::INV))
+    if((feats[0]->type() == NODE_TYPE::INV) || (feats[1]->type() == NODE_TYPE::INV) || (feats[0]->type() == NODE_TYPE::DIV) || (feats[1]->type() == NODE_TYPE::DIV))
         throw InvalidFeatureException();
 
     std::map<std::string, double> div_mult_leaves;
@@ -31,10 +31,10 @@ DivNode::DivNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l
     set_test_value();
 }
 
-DivNode::DivNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat_1, feat_2}, rung, feat_ind)
+DivNode::DivNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat_1, feat_2}, feat_ind)
 {
-    if((feat_1->type() == NODE_TYPE::INV) || (feat_2->type() == NODE_TYPE::INV))
+    if((feat_1->type() == NODE_TYPE::INV) || (feat_2->type() == NODE_TYPE::INV) || (feat_1->type() == NODE_TYPE::DIV) || (feat_2->type() == NODE_TYPE::DIV))
         throw InvalidFeatureException();
 
     std::map<std::string, double> div_mult_leaves;
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.hpp
index 7b6e07d7..bf6a279b 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/divide.hpp
@@ -15,9 +15,9 @@ class DivNode: public OperatorNode<2>
 public:
     DivNode();
 
-    DivNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    DivNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound);
 
-    DivNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound);
+    DivNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit() / _feats[1]->unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.cpp
index 92a529b9..13f4f5b7 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.cpp
@@ -3,8 +3,8 @@
 ExpNode::ExpNode()
 {}
 
-ExpNode::ExpNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+ExpNode::ExpNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != Unit())
         throw InvalidFeatureException();
@@ -19,8 +19,8 @@ ExpNode::ExpNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l
        set_test_value();
  }
 
-ExpNode::ExpNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+ExpNode::ExpNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->unit() != Unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.hpp
index 7870d772..13015a62 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/exponential.hpp
@@ -15,9 +15,9 @@ class ExpNode: public OperatorNode<1>
 public:
     ExpNode();
 
-    ExpNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    ExpNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    ExpNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    ExpNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return Unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.cpp
index b7001b0f..65eea3bd 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.cpp
@@ -3,8 +3,8 @@
 InvNode::InvNode()
 {}
 
-InvNode::InvNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+InvNode::InvNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if((feats[0]->type() == NODE_TYPE::DIV) || (feats[0]->type() == NODE_TYPE::EXP) || (feats[0]->type() == NODE_TYPE::NEG_EXP))
         throw InvalidFeatureException();
@@ -16,8 +16,8 @@ InvNode::InvNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l
     set_test_value();
 }
 
-InvNode::InvNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+InvNode::InvNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if((feat->type() == NODE_TYPE::DIV) || (feat->type() == NODE_TYPE::EXP) || (feat->type() == NODE_TYPE::NEG_EXP))
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.hpp
index 384fff3d..7cc7007b 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/inverse.hpp
@@ -16,9 +16,9 @@ class InvNode: public OperatorNode<1>
 public:
     InvNode();
 
-    InvNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    InvNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    InvNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    InvNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit() ^ (-1.0);}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.cpp
index 7de0fbd6..112709b2 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.cpp
@@ -3,8 +3,8 @@
 LogNode::LogNode()
 {}
 
-LogNode::LogNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+LogNode::LogNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != Unit())
         throw InvalidFeatureException();
@@ -19,8 +19,8 @@ LogNode::LogNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l
        set_test_value();
  }
 
-LogNode::LogNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+LogNode::LogNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->unit() != Unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.hpp
index 11f0e846..e6fe8b2c 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/log.hpp
@@ -15,9 +15,9 @@ class LogNode: public OperatorNode<1>
 public:
     LogNode();
 
-    LogNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    LogNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    LogNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    LogNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return Unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.cpp
index 0d202151..01f44057 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.cpp
@@ -3,9 +3,12 @@
 MultNode::MultNode()
 {}
 
-MultNode::MultNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+MultNode::MultNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
+    if((feats[0]->type() == NODE_TYPE::INV) || (feats[1]->type() == NODE_TYPE::INV) || (feats[0]->type() == NODE_TYPE::DIV) || (feats[1]->type() == NODE_TYPE::DIV))
+        throw InvalidFeatureException();
+
     std::map<std::string, double> div_mult_leaves;
     double expected_abs_tot = 0.0;
     update_div_mult_leaves(div_mult_leaves, 1.0, expected_abs_tot);
@@ -28,9 +31,12 @@ MultNode::MultNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double
     set_test_value();
 }
 
-MultNode::MultNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat_1, feat_2}, rung, feat_ind)
+MultNode::MultNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat_1, feat_2}, feat_ind)
 {
+    if((feat_1->type() == NODE_TYPE::INV) || (feat_2->type() == NODE_TYPE::INV) || (feat_1->type() == NODE_TYPE::DIV) || (feat_2->type() == NODE_TYPE::DIV))
+        throw InvalidFeatureException();
+
     std::map<std::string, double> div_mult_leaves;
     double expected_abs_tot = 0.0;
     update_div_mult_leaves(div_mult_leaves, 1.0, expected_abs_tot);
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.hpp
index 58e8df98..f7187cbd 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/multiply.hpp
@@ -16,9 +16,9 @@ class MultNode: public OperatorNode<2>
 public:
     MultNode();
 
-    MultNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    MultNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound);
 
-    MultNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound);
+    MultNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit() * _feats[1]->unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.cpp
index 5183f5f1..063c5dbd 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.cpp
@@ -3,8 +3,8 @@
 NegExpNode::NegExpNode()
 {}
 
-NegExpNode::NegExpNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+NegExpNode::NegExpNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != Unit())
         throw InvalidFeatureException();
@@ -19,8 +19,8 @@ NegExpNode::NegExpNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, do
     set_test_value();
  }
 
-NegExpNode::NegExpNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+NegExpNode::NegExpNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->unit() != Unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.hpp
index 2b509ade..4e8af7a3 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/negative_exponential.hpp
@@ -16,9 +16,9 @@ class NegExpNode: public OperatorNode<1>
 public:
     NegExpNode();
 
-    NegExpNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    NegExpNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    NegExpNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    NegExpNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return Unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.cpp
index d928cbd6..8bb02076 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.cpp
@@ -3,8 +3,8 @@
 SinNode::SinNode()
 {}
 
-SinNode::SinNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+SinNode::SinNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != Unit())
         throw InvalidFeatureException();
@@ -19,8 +19,8 @@ SinNode::SinNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l
     set_test_value();
  }
 
-SinNode::SinNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+SinNode::SinNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->unit() != Unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.hpp
index 4654cc38..62115fae 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sin.hpp
@@ -16,9 +16,9 @@ class SinNode: public OperatorNode<1>
 public:
     SinNode();
 
-    SinNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    SinNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    SinNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    SinNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return Unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.cpp
index 2f7fe9a9..0fd1a4cc 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.cpp
@@ -3,8 +3,8 @@
 SixPowNode::SixPowNode()
 {}
 
-SixPowNode::SixPowNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+SixPowNode::SixPowNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if((feats[0]->type() == NODE_TYPE::CBRT) || (feats[0]->type() == NODE_TYPE::SQRT))
         throw InvalidFeatureException();
@@ -16,8 +16,8 @@ SixPowNode::SixPowNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, do
     set_test_value();
 }
 
-SixPowNode::SixPowNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+SixPowNode::SixPowNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if((feat->type() == NODE_TYPE::CBRT) || (feat->type() == NODE_TYPE::SQRT))
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.hpp
index 8329bae7..4e2bab71 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/sixth_power.hpp
@@ -16,9 +16,9 @@ class SixPowNode: public OperatorNode<1>
 public:
     SixPowNode();
 
-    SixPowNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    SixPowNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    SixPowNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    SixPowNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit()^(6.0);}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.cpp
index 1bdbd3c0..d286b588 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.cpp
@@ -3,8 +3,8 @@
 SqNode::SqNode()
 {}
 
-SqNode::SqNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+SqNode::SqNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->type() == NODE_TYPE::SQRT)
         throw InvalidFeatureException();
@@ -16,8 +16,8 @@ SqNode::SqNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_b
     set_test_value();
 }
 
-SqNode::SqNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+SqNode::SqNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->type() == NODE_TYPE::SQRT)
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.hpp
index 449c6b48..96c3dde4 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square.hpp
@@ -15,9 +15,9 @@ class SqNode: public OperatorNode<1>
 public:
     SqNode();
 
-    SqNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    SqNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    SqNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    SqNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit()^(2.0);}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.cpp
index 090c37d5..727c58a5 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.cpp
@@ -3,8 +3,8 @@
 SqrtNode::SqrtNode()
 {}
 
-SqrtNode::SqrtNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+SqrtNode::SqrtNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->type() == NODE_TYPE::SQRT)
         throw InvalidFeatureException();
@@ -16,8 +16,8 @@ SqrtNode::SqrtNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double
     set_test_value();
 }
 
-SqrtNode::SqrtNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat}, rung, feat_ind)
+SqrtNode::SqrtNode(node_ptr feat, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat}, feat_ind)
 {
     if(feat->type() == NODE_TYPE::SQRT)
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.hpp
index 1cfbe8ce..53520d78 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/square_root.hpp
@@ -16,9 +16,9 @@ class SqrtNode: public OperatorNode<1>
 public:
     SqrtNode();
 
-    SqrtNode(std::array<node_ptr, 1> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    SqrtNode(std::array<node_ptr, 1> feats, int feat_ind, double l_bound, double u_bound);
 
-    SqrtNode(node_ptr feat, int rung, int feat_ind, double l_bound, double u_bound);
+    SqrtNode(node_ptr feat, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit()^(0.5);}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.cpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.cpp
index 9af1f605..f3a65476 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.cpp
@@ -3,8 +3,8 @@
 SubNode::SubNode()
 {}
 
-SubNode::SubNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode(feats, rung, feat_ind)
+SubNode::SubNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound):
+    OperatorNode(feats, feat_ind)
 {
     if(feats[0]->unit() != feats[1]->unit())
         throw InvalidFeatureException();
@@ -31,8 +31,8 @@ SubNode::SubNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l
        set_test_value();
  }
 
-SubNode::SubNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound):
-    OperatorNode({feat_1, feat_2}, rung, feat_ind)
+SubNode::SubNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound):
+    OperatorNode({feat_1, feat_2}, feat_ind)
 {
     if(feat_1->unit() != feat_2->unit())
         throw InvalidFeatureException();
diff --git a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.hpp b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.hpp
index 0ea056ca..3962cf73 100644
--- a/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_operator_nodes/subtract.hpp
@@ -16,9 +16,9 @@ class SubNode: public OperatorNode<2>
 public:
     SubNode();
 
-    SubNode(std::array<node_ptr, 2> feats, int rung, int feat_ind, double l_bound, double u_bound);
+    SubNode(std::array<node_ptr, 2> feats, int feat_ind, double l_bound, double u_bound);
 
-    SubNode(node_ptr feat_1, node_ptr feat_2, int rung, int feat_ind, double l_bound, double u_bound);
+    SubNode(node_ptr feat_1, node_ptr feat_2, int feat_ind, double l_bound, double u_bound);
 
     inline Unit unit(){return _feats[0]->unit();}
 
diff --git a/src/feature_creation/node/operator_nodes/allowed_ops.cpp b/src/feature_creation/node/operator_nodes/allowed_ops.cpp
index f93bd046..68ff28ef 100644
--- a/src/feature_creation/node/operator_nodes/allowed_ops.cpp
+++ b/src/feature_creation/node/operator_nodes/allowed_ops.cpp
@@ -5,22 +5,22 @@ std::map<std::string, bin_op_node_gen> allowed_op_maps::binary_operator_map;
 
 void allowed_op_maps::set_node_maps()
 {
-    allowed_op_maps::binary_operator_map["add"] = [](node_ptr f1, node_ptr f2, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<AddNode>(f1, f2, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::binary_operator_map["sub"] = [](node_ptr f1, node_ptr f2, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<SubNode>(f1, f2, rung, feat_ind, l_bound, u_bound);};;
-    allowed_op_maps::binary_operator_map["abs_diff"] = [](node_ptr f1, node_ptr f2, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<AbsDiffNode>(f1, f2, rung, feat_ind, l_bound, u_bound);};;
-    allowed_op_maps::binary_operator_map["mult"] = [](node_ptr f1, node_ptr f2, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<MultNode>(f1, f2, rung, feat_ind, l_bound, u_bound);};;
-    allowed_op_maps::binary_operator_map["div"] = [](node_ptr f1, node_ptr f2, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<DivNode>(f1, f2, rung, feat_ind, l_bound, u_bound);};;
+    allowed_op_maps::binary_operator_map["add"] = [](node_ptr f1, node_ptr f2, int feat_ind, double l_bound, double u_bound){return std::make_shared<AddNode>(f1, f2, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::binary_operator_map["sub"] = [](node_ptr f1, node_ptr f2, int feat_ind, double l_bound, double u_bound){return std::make_shared<SubNode>(f1, f2, feat_ind, l_bound, u_bound);};;
+    allowed_op_maps::binary_operator_map["abs_diff"] = [](node_ptr f1, node_ptr f2, int feat_ind, double l_bound, double u_bound){return std::make_shared<AbsDiffNode>(f1, f2, feat_ind, l_bound, u_bound);};;
+    allowed_op_maps::binary_operator_map["mult"] = [](node_ptr f1, node_ptr f2, int feat_ind, double l_bound, double u_bound){return std::make_shared<MultNode>(f1, f2, feat_ind, l_bound, u_bound);};;
+    allowed_op_maps::binary_operator_map["div"] = [](node_ptr f1, node_ptr f2, int feat_ind, double l_bound, double u_bound){return std::make_shared<DivNode>(f1, f2, feat_ind, l_bound, u_bound);};;
 
-    allowed_op_maps::unary_operator_map["exp"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<ExpNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["neg_exp"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<NegExpNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["inv"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<InvNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["sq"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<SqNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["cb"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<CbNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["six_pow"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<SixPowNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["sqrt"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<SqrtNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["cbrt"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<CbrtNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["log"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<LogNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["abs"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<AbsNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["sin"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<SinNode>(f1, rung, feat_ind, l_bound, u_bound);};
-    allowed_op_maps::unary_operator_map["cos"] = [](node_ptr f1, int rung, int feat_ind, double l_bound, double u_bound){return std::make_shared<CosNode>(f1, rung, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["exp"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<ExpNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["neg_exp"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<NegExpNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["inv"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<InvNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["sq"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<SqNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["cb"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<CbNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["six_pow"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<SixPowNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["sqrt"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<SqrtNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["cbrt"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<CbrtNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["log"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<LogNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["abs"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<AbsNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["sin"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<SinNode>(f1, feat_ind, l_bound, u_bound);};
+    allowed_op_maps::unary_operator_map["cos"] = [](node_ptr f1, int feat_ind, double l_bound, double u_bound){return std::make_shared<CosNode>(f1, feat_ind, l_bound, u_bound);};
 }
diff --git a/src/feature_creation/node/operator_nodes/allowed_ops.hpp b/src/feature_creation/node/operator_nodes/allowed_ops.hpp
index 00d1e526..045c3792 100644
--- a/src/feature_creation/node/operator_nodes/allowed_ops.hpp
+++ b/src/feature_creation/node/operator_nodes/allowed_ops.hpp
@@ -22,8 +22,8 @@
 #include <map>
 #include <iostream>
 
-typedef std::function<node_ptr(node_ptr, int, int, double, double)> un_op_node_gen;
-typedef std::function<node_ptr(node_ptr, node_ptr, int, int, double, double)> bin_op_node_gen;
+typedef std::function<node_ptr(node_ptr, int, double, double)> un_op_node_gen;
+typedef std::function<node_ptr(node_ptr, node_ptr, int, double, double)> bin_op_node_gen;
 
 namespace allowed_op_maps
 {
diff --git a/src/feature_creation/units/Unit.cpp b/src/feature_creation/units/Unit.cpp
index 0347bfdc..c384ea85 100644
--- a/src/feature_creation/units/Unit.cpp
+++ b/src/feature_creation/units/Unit.cpp
@@ -129,14 +129,14 @@ bool Unit::equal(Unit unit_2)
 {
     for(auto& el : unit_2.dct())
     {
-        if(_dct.count(el.first) == 0)
+        if((_dct.count(el.first) == 0) && (el.second != 0))
             return false;
         else if(_dct[el.first] != el.second)
             return false;
     }
     for(auto& el : _dct)
     {
-        if(unit_2.dct().count(el.first) == 0)
+        if((unit_2.dct().count(el.first) == 0) && (el.second != 0))
             return false;
         else if(unit_2.dct()[el.first] != el.second)
             return false;
diff --git a/src/inputs/InputParser.cpp b/src/inputs/InputParser.cpp
index 58063aad..c20264c1 100644
--- a/src/inputs/InputParser.cpp
+++ b/src/inputs/InputParser.cpp
@@ -12,6 +12,7 @@ InputParser::InputParser(boost::property_tree::ptree IP, std::string fn, std::sh
     _n_sis_select(IP.get<int>("_n_sis_select")),
     _max_rung(IP.get<int>("max_rung")),
     _max_store_rung(IP.get<int>("n_rung_store", _max_rung - 1)),
+    _n_rung_generate(IP.get<int>("n_rung_generate", 0)),
     _n_samp(-1),
     _n_residuals(IP.get<int>("n_residual", 1))
 {
@@ -174,12 +175,11 @@ void InputParser::generate_feature_space(std::shared_ptr<MPI_Interface> comm)
     }
 
     node_value_arrs::initialize_values_arr(_prop_train.size(), _prop_test.size(), headers.size());
-
     std::vector<node_ptr> phi_0;
     for(int ff = 0; ff < headers.size(); ++ff)
         phi_0.push_back(std::make_shared<FeatureNode>(ff, headers[ff], data[ff], test_data[ff], units[ff]));
 
-    _feat_space = std::make_shared<FeatureSpace>(comm, phi_0, _opset, _max_rung, _n_sis_select, _max_store_rung, _l_bound, _u_bound);
+    _feat_space = std::make_shared<FeatureSpace>(comm, phi_0, _opset, _max_rung, _n_sis_select, _max_store_rung, _n_rung_generate, _l_bound, _u_bound);
 }
 
 void stripComments(std::string& filename)
diff --git a/src/inputs/InputParser.hpp b/src/inputs/InputParser.hpp
index a3801fac..8cab6a31 100644
--- a/src/inputs/InputParser.hpp
+++ b/src/inputs/InputParser.hpp
@@ -42,6 +42,7 @@ public:
     int _n_dim;
     int _max_rung;
     int _max_store_rung;
+    int _n_rung_generate;
     int _n_sis_select;
     int _n_samp;
     int _n_residuals;
-- 
GitLab