diff --git a/src/descriptor_identifier/Model/Model.cpp b/src/descriptor_identifier/Model/Model.cpp
index b83ede88dcb51c1a8343db064b6775a362abab68..1085dcea0b6cec22dd0a7abb45181eee38c8b0f5 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 574a8abbb16adeeca39203367c6fd86be11c973f..261995eb95ab09f2425073a927b3389468172180 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 1c547dd1750ee2d72186b6aec63b19e843ec46ca..63dff34d862a482d024bcd38de27020ff951e24b 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 38e8e96a6c6b20040c16beadc63d20167db01b9b..eccc5af03b263215bc969d159871864f1fd550e7 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 de7f41a98ff9db7a4a0aebc5273d7b43486a215d..4cf22bc2e938c426250caf8f82e1d4587575bdfa 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 7d8ac1f38402eec8be77f1a0ec6a34ec07f69457..a808977ba028920c6ae967670d6bae4ce0c931d7 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 284cb926a8d00f1f249cf6aac973f4b57e0f49e2..ff643fc29275ff5d6b477caca112a2aa991da7a8 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 a4602172f0ef079f85a3c79d7379d507f568500c..3ee28fe9710413cf68ae75b3d1ad1933303d3538 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 1d1ee42363d1acbe2f46076f57d66b369e06d6e5..7e3b252b8d1019626875a5db3a15c5071653a4ca 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 4feb01f2f9f481e69235597105ce8a93b4b1841e..bf5aa248acc4eff96db28ca79226c24f83f6790d 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 0aca901696c85e45f9eeaf7517b29359cc2f2e80..2179a7c4578b27f7076e8994c07e8cdbbf11ba4f 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 670f4592bf3f15171eeeeaa6f5dfb6dc2aea147d..90ee52bf8f6c2bec4591159db37807a1fa4ebd1d 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 c42265ca26a775239b5e739956f918f504514304..a7a111f9703960660ea5436248b0cc80ec989e8f 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 b764c054221a14c28feabd3420329557018f3041..4786123b978114a1143d1edc17a6683759be4743 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 b2579b811cc2eb09e53f2cd259c4bd53f6a31788..d722fc8ed497e45d2e840926243af6b9a59ff58e 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 2744e85365d73b76ed4557f547ff7f6ab96b8782..9bf1fab672a151297a12d09282f9d7bb0bd01d7e 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 c006a46bf64ae01c83d912ca6350b54cf5407f42..66799f4e61f13b46f3dfd2ff0d3b8e9e4df23a9f 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 f9b47cf72ac18e2653c568661926ec2ca2d8e188..5ad863e100faaecc4a0c0dd7084a4d1f37df2554 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 3d7eac11b6eedf921111ae9cc1cb7a726928aa5d..b110dc91685a2c2a60f6acdda60c8804e8943224 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 a679933caea178b62ca8cdcbced5f195dd76c523..2cc4491830ca98c6a6e893f79718663fb6d4fdb0 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 1509726f0154c30c4ed8fe0912370a8e040f381c..faff010878e5367c8b505c2f9a6f364fe0859150 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 b0d4dcfc136928391bb41e025724db9549587733..f5d656670953325d9871321a286d7407f912db16 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 7b6e07d79131e43f13b1b0ad1bfa0a8802b1847e..bf6a279b5f55c3709d45ced53ea96c650e0a2359 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 92a529b9dd1bb5ebc706e4522648b61765624c93..13f4f5b7f0bbf598e752f2b7f774628233c1ebea 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 7870d772afdd52afba999ad72869e264258a01f3..13015a62dafe94bf6bd68a4758b25375f9be4ec7 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 b7001b0ff3ba80919f4beed68ad179fdaa58e334..65eea3bd56696422f75abc18df3fb00a1504c027 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 384fff3d18995d74e0e834b55ef95dff0613364a..7cc7007bb4bbc32e22118b6bad53367436deed25 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 7de0fbd69edf94c83b1ded58a4bd15b38ece9f42..112709b2fbd91ee01378c1d2bfc3daa6284b9dba 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 11f0e8460e89164646f05b11928d5dfeef40bd65..e6fe8b2ca77f7969bc8455baa9984debca8a58c9 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 0d202151176e450eb2f86b0467bd1aa1a5ba63ae..01f44057e0884a1106cd766ea880ebd24c6a2bfe 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 58e8df98d9e3ec517a3dd9ea0d378d4568acdbfc..f7187cbdb50cdbdf90eca427c55bd31136a9e34c 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 5183f5f131e52495152823be49d29cd1cca7b9b5..063c5dbdd73d66f0b96e162201f39577cfc85d23 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 2b509ade484936ee5b4bc5d02914c0429cab886c..4e8af7a3807f99a23cbbb7a70a93989def66198a 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 d928cbd650d82b879c26e56fa8783d230707abd7..8bb02076c009281fba23b8e99991b5f4239d5641 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 4654cc380d51e382f9a324781a05047d87b16a54..62115fae725087829ef0be4a67ae207657663c17 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 2f7fe9a9f0f85c597d037430a337153de81e509e..0fd1a4cc99c112459edd1d694b30a3407bb678e1 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 8329bae7290c5381ab1a62ac8963747e4a6e572a..4e2bab71a4d5eb672476c434a416f697d6c964ae 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 1bdbd3c0eda1d9f4c03ef272e2c8a29440c04001..d286b588a87588845522c1d9a03699eaccaa7013 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 449c6b48984c22aa5e5251497bdeb0fa35cde007..96c3dde4c3514cfb22f8425a2a585dc3e6454626 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 090c37d5ff603ebd3ae9781088a953ea7f08e454..727c58a5b66fb05f9e88345eaa168d4531167b2e 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 1cfbe8ce1de74e0140f7645ad1b9ab53375777cf..53520d78a289bf764e9b41a25b8cde107e6a84bf 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 9af1f605aa7047d3d85c81bfcf04a5d4f0c609b9..f3a654766d2d022cf2c711fdf39fb429b57e2ef9 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 0ea056cab5ee220c77a872c9a85dada086494c2b..3962cf7309f61e268b6336fd680eeda83f3d00d5 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 f93bd046e01a6059e1e11e4fe26a53d500c8a5c1..68ff28ef91e5e377210a04184409badd68f7d0ae 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 00d1e526da480861413d7e590360837a50c7df7e..045c3792a859974c000db7a90485fe7db05f8b62 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 0347bfdcd2c7c38d9a818c165ad052184dc45025..c384ea8583f3591f784833413a286f3643f356e6 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 58063aadffccc28f009d52068e578ef50d41c200..c20264c13d0107ef10c54de977a7d82885716e6a 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 a3801fac23842fe0a9107a8d36bf9f7104a4c413..8cab6a31110ac5da267d3198157d003b2aa7b95a 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;