diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1547b6baec37e1ae7ca3ef63fd14ad4fdc8071cd..f8f3f85be623a01ed3bd77cd8aa65688c3ca1985 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -16,64 +16,77 @@ set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
-# Check python  settings
-find_package(PythonInterp 3 REQUIRED)
-execute_process(COMMAND ${PYTHON_EXECUTABLE}
-                -c "import distutils.sysconfig as cg; print(cg.get_python_inc())"
-                OUTPUT_VARIABLE PYTHON_INCLUDE_PATH
-                OUTPUT_STRIP_TRAILING_WHITESPACE)
-
-set(PYTHON_INCLUDE_PATH ${PYTHON_INCLUDE_PATH} CACHE PATH "Python Include Directory")
-mark_as_advanced(PYTHON_INCLUDE_PATH)
-
-include_directories(${PYTHON_INCLUDE_PATH})
-
-if(NOT PYTHON_INSTDIR)
-execute_process(COMMAND ${PYTHON_EXECUTABLE}
-                -c "import distutils.sysconfig as cg; print(cg.get_python_lib(1,0))"
-                OUTPUT_VARIABLE PYTHON_INSTDIR OUTPUT_STRIP_TRAILING_WHITESPACE)
+option(USE_PYTHON "Whether to compile with python binding support" OFF)
+
+if(USE_PYTHON)
+    message(STATUS "USE PYTHON True")
+    set(USE_PYTHON TRUE)
+else(USE_PYTHON)
+    message(STATUS "USE PYTHON False")
+    set(USE_PYTHON FALSE)
 endif()
 
-execute_process(COMMAND ${PYTHON_EXECUTABLE}
-                -c "import sys; print('{}{}'.format(sys.version_info[0], sys.version_info[1]))"
-                OUTPUT_VARIABLE BOOST_PYTHON_VERSION
-                OUTPUT_STRIP_TRAILING_WHITESPACE)
-
-execute_process(COMMAND ${PYTHON_EXECUTABLE}
-                        -c "import sys; print(sys.version[:3])"
-                        OUTPUT_VARIABLE PYTHON_VERSION
-                        OUTPUT_STRIP_TRAILING_WHITESPACE)
-
-execute_process(COMMAND ${PYTHON_EXECUTABLE}
-                -c "import distutils.sysconfig as cg; print(cg.get_config_var('LIBDIR'))"
-                OUTPUT_VARIABLE PYTHON_LIBDIR
-                OUTPUT_STRIP_TRAILING_WHITESPACE)
-
-message(STATUS "PYTHON_LIBDIR = ${PYTHON_LIBDIR}")
-message(STATUS "PYTHON_INSTDIR = ${PYTHON_INSTDIR}")
-
-find_library(PYTHON_LIBRARIES
-             NAMES python${PYTHON_VERSION}m
-             PATHS ${PYTHON_LIBDIR}  )
-
-if(NOT PYTHON_LIBRARIES)
-  message(FATAL_ERROR "Python libraries not found!")
-endif(NOT PYTHON_LIBRARIES)
-message(STATUS "PYTHON_LIBRARIES = ${PYTHON_LIBRARIES}")
-
-mark_as_advanced(PYTHON_INCLUDE_PATH PYTHON_LIBRARIES)
+if(USE_PYTHON)
+    # Check python  settings
+    find_package(PythonInterp 3 REQUIRED)
+    execute_process(COMMAND ${PYTHON_EXECUTABLE}
+                    -c "import distutils.sysconfig as cg; print(cg.get_python_inc())"
+                    OUTPUT_VARIABLE PYTHON_INCLUDE_PATH
+                    OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+    set(PYTHON_INCLUDE_PATH ${PYTHON_INCLUDE_PATH} CACHE PATH "Python Include Directory")
+    mark_as_advanced(PYTHON_INCLUDE_PATH)
+
+    include_directories(${PYTHON_INCLUDE_PATH})
+
+    if(NOT PYTHON_INSTDIR)
+    execute_process(COMMAND ${PYTHON_EXECUTABLE}
+                    -c "import distutils.sysconfig as cg; print(cg.get_python_lib(1,0))"
+                    OUTPUT_VARIABLE PYTHON_INSTDIR OUTPUT_STRIP_TRAILING_WHITESPACE)
+    endif()
+
+    execute_process(COMMAND ${PYTHON_EXECUTABLE}
+                    -c "import sys; print('{}{}'.format(sys.version_info[0], sys.version_info[1]))"
+                    OUTPUT_VARIABLE BOOST_PYTHON_VERSION
+                    OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+    execute_process(COMMAND ${PYTHON_EXECUTABLE}
+                            -c "import sys; print(sys.version[:3])"
+                            OUTPUT_VARIABLE PYTHON_VERSION
+                            OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+    execute_process(COMMAND ${PYTHON_EXECUTABLE}
+                    -c "import distutils.sysconfig as cg; print(cg.get_config_var('LIBDIR'))"
+                    OUTPUT_VARIABLE PYTHON_LIBDIR
+                    OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+    message(STATUS "PYTHON_LIBDIR = ${PYTHON_LIBDIR}")
+    message(STATUS "PYTHON_INSTDIR = ${PYTHON_INSTDIR}")
+
+    find_library(PYTHON_LIBRARIES
+                 NAMES python${PYTHON_VERSION}m
+                 PATHS ${PYTHON_LIBDIR}  )
+
+    if(NOT PYTHON_LIBRARIES)
+      message(FATAL_ERROR "Python libraries not found!")
+    endif(NOT PYTHON_LIBRARIES)
+    message(STATUS "PYTHON_LIBRARIES = ${PYTHON_LIBRARIES}")
+
+    mark_as_advanced(PYTHON_INCLUDE_PATH PYTHON_LIBRARIES)
+endif()
 
 # Check for Boost
 find_package(Boost REQUIRED COMPONENTS filesystem system mpi serialization)
 set(Boost_LIBS ${Boost_LIBRARIES})
 
-if(${Boost_VERSION} VERSION_LESS 106700)
- find_package(Boost ${NEEDED_Boost_VERSION} REQUIRED COMPONENTS python numpy)
-else()
- find_package(Boost ${NEEDED_Boost_VERSION} REQUIRED COMPONENTS python${BOOST_PYTHON_VERSION} numpy${BOOST_PYTHON_VERSION})
+if(USE_PYTHON)
+    if(${Boost_VERSION} VERSION_LESS 106700)
+     find_package(Boost ${NEEDED_Boost_VERSION} REQUIRED COMPONENTS python numpy)
+    else()
+     find_package(Boost ${NEEDED_Boost_VERSION} REQUIRED COMPONENTS python${BOOST_PYTHON_VERSION} numpy${BOOST_PYTHON_VERSION})
+    endif()
 endif()
 
-
 # Append Python library to the list of Boost libraries.
 list(APPEND Boost_LIBS ${Boost_LIBRARIES})
 set(Boost_LIBRARIES ${Boost_LIBS})
@@ -90,7 +103,7 @@ set(BLA_VENDOR Intel10_64lp_seq)
 find_package(LAPACK)
 if(NOT LAPACK_FOUND)
     set(BLA_VENDOR All)
-    find_package(LAPACK REQUIRED)	
+    find_package(LAPACK REQUIRED)
 endif()
 message(STATUS "LAPACK_LIBRARIES = ${LAPACK_LIBRARIES}")
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index feb824390d88501ce3bedbb35d1518ff1ea01e30..e039b9956d33d258e36f7e18e38137b3204e1a77 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,18 +1,3 @@
-file(GLOB_RECURSE SISSOPP_SOURCES *.cpp)
-file(GLOB_RECURSE SISSOLIB_SOURCES *.cpp)
-
-file(GLOB_RECURSE NOT_SISSOPP_SOURCES python/*.cpp)
-
-list(REMOVE_ITEM SISSOPP_SOURCES ${NOT_SISSOPP_SOURCES})
-list(REMOVE_ITEM SISSOLIB_SOURCES main.cpp)
-
-configure_file(${CMAKE_CURRENT_SOURCE_DIR}/sisso++_config.hpp.in ${CMAKE_CURRENT_BINARY_DIR}/sisso++_config.hpp)
-
-add_executable(sisso++  ${SISSOPP_SOURCES})
-
-add_library(_sisso SHARED ${SISSOLIB_SOURCES})
-
-
 include_directories(${CMAKE_CURRENT_LIST_DIR}/descriptor_identifier/)
 include_directories(${CMAKE_CURRENT_LIST_DIR}/descriptor_identifier/Model)
 include_directories(${CMAKE_CURRENT_LIST_DIR}/feature_creation/node/)
@@ -27,6 +12,14 @@ include_directories(${CMAKE_CURRENT_LIST_DIR}/mpi_interface/)
 include_directories(${CMAKE_CURRENT_LIST_DIR}/python/)
 include_directories(${CMAKE_CURRENT_LIST_DIR}/utils/)
 
+file(GLOB_RECURSE SISSOPP_SOURCES *.cpp)
+file(GLOB_RECURSE NOT_SISSOPP_SOURCES python/*.cpp)
+list(REMOVE_ITEM SISSOPP_SOURCES ${NOT_SISSOPP_SOURCES})
+
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/sisso++_config.hpp.in ${CMAKE_CURRENT_BINARY_DIR}/sisso++_config.hpp)
+
+add_executable(sisso++  ${SISSOPP_SOURCES})
+
 set_target_properties(sisso++
     PROPERTIES
     ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
@@ -34,22 +27,29 @@ set_target_properties(sisso++
     RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
 )
 target_link_libraries(sisso++  ${LAPACK_LIBRARIES} ${MPI_LIBRARIES} ${PYTHON_LIBRARIES} ${Boost_LIBS})
-
-set_target_properties(_sisso
-    PROPERTIES
-    ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
-    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
-    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
-    PREFIX ""
-    SUFFIX ".so"
-)
-target_link_libraries(_sisso  ${LAPACK_LIBRARIES} ${MPI_LIBRARIES} ${PYTHON_LIBRARIES} ${Boost_LIBS})
-
-configure_file(${CMAKE_CURRENT_SOURCE_DIR}/python/__init__.py ${CMAKE_CURRENT_LIST_DIR}/python/__init__.py COPYONLY)
-
-install(TARGETS _sisso DESTINATION "${PYTHON_INSTDIR}/sisso")
-install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/python/ DESTINATION ${PYTHON_INSTDIR}/sisso
-    FILES_MATCHING PATTERN "*.py"
-    PATTERN "CMakeFiles" EXCLUDE)
-
-install(TARGETS sisso++ DESTINATION ${CMAKE_CURRENT_LIST_DIR}/../bin/)
\ No newline at end of file
+install(TARGETS sisso++ DESTINATION ${CMAKE_CURRENT_LIST_DIR}/../bin/)
+
+if(USE_PYTHON)
+    file(GLOB_RECURSE SISSOLIB_SOURCES *.cpp)
+    list(REMOVE_ITEM SISSOLIB_SOURCES main.cpp)
+
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DPY_BINDINGS")
+    add_library(_sisso SHARED ${SISSOLIB_SOURCES})
+
+    configure_file(${CMAKE_CURRENT_SOURCE_DIR}/python/__init__.py ${CMAKE_CURRENT_LIST_DIR}/python/__init__.py COPYONLY)
+
+    set_target_properties(_sisso
+        PROPERTIES
+        ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+        LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+        RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+        PREFIX ""
+        SUFFIX ".so"
+    )
+    target_link_libraries(_sisso  ${LAPACK_LIBRARIES} ${MPI_LIBRARIES} ${PYTHON_LIBRARIES} ${Boost_LIBS})
+
+    install(TARGETS _sisso DESTINATION "${PYTHON_INSTDIR}/sisso")
+    install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/python/ DESTINATION ${PYTHON_INSTDIR}/sisso
+        FILES_MATCHING PATTERN "*.py"
+        PATTERN "CMakeFiles" EXCLUDE)
+endif()
diff --git a/src/descriptor_identifier/Model/Model.cpp b/src/descriptor_identifier/Model/Model.cpp
index 01122e53e8942967e986c38aaa79dbfddb2be120..00c227fbe710874162aaf73e34582d1285932087 100644
--- a/src/descriptor_identifier/Model/Model.cpp
+++ b/src/descriptor_identifier/Model/Model.cpp
@@ -331,30 +331,3 @@ void Model::to_file(std::string filename, bool train, std::vector<int> test_inds
     }
     out_file_stream.close();
 }
-
-void Model::register_python()
-{
-    using namespace boost::python;
-    class_<Model>("Model", init<std::vector<double>, std::vector<double>, std::vector<model_node_ptr>, std::vector<int>, std::vector<int>>())
-        .def(init<std::string>())
-        .def(init<std::string, std::string>())
-        .def("predict", &Model::predict)
-        .def("fit", &Model::predict_train)
-        .def("__str__", &Model::toString)
-        .def("__repr__", &Model::toString)
-        .def_readonly("_n_samp_train", &Model::_n_samp_train)
-        .def_readonly("_n_samp_test", &Model::_n_samp_test)
-        .def_readonly("_n_dim", &Model::_n_dim)
-        .add_property("prop_train_est", &Model::prop_train_est)
-        .add_property("prop_test_est", &Model::prop_test_est)
-        .add_property("prop_train", &Model::prop_train)
-        .add_property("prop_test", &Model::prop_test)
-        .add_property("train_error", &Model::train_error)
-        .add_property("test_error", &Model::test_error)
-        .add_property("feats", &Model::feats)
-        .add_property("coefs", &Model::coefs)
-        .add_property("rmse", &Model::rmse)
-        .add_property("test_rmse", &Model::test_rmse)
-        .add_property("max_ae", &Model::max_ae)
-        .add_property("test_max_ae", &Model::test_max_ae);
-}
diff --git a/src/descriptor_identifier/Model/Model.hpp b/src/descriptor_identifier/Model/Model.hpp
index e0a47000f1809b8e5eb4b0a80c2601bd41e46f3b..99baeb4d63612156f437b326c267846208d01d8f 100644
--- a/src/descriptor_identifier/Model/Model.hpp
+++ b/src/descriptor_identifier/Model/Model.hpp
@@ -4,7 +4,6 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/algorithm/string/trim.hpp>
 #include <boost/filesystem.hpp>
-#include <boost/python.hpp>
 
 #include<iomanip>
 #include<fstream>
@@ -13,8 +12,10 @@
 #include <feature_creation/node/ModelNode.hpp>
 #include <utils/string_utils.hpp>
 
-namespace python = boost::python;
-namespace np = boost::python::numpy;
+#ifdef PY_BINDINGS
+    namespace np = boost::python::numpy;
+    namespace py = boost::python;
+#endif
 
 typedef std::shared_ptr<ModelNode> model_node_ptr;
 /**
@@ -85,6 +86,9 @@ public:
      */
     inline double rmse(){return util_funcs::norm(_train_error.data(), _n_samp_train) / std::sqrt(static_cast<double>(_n_samp_train));}
     inline double test_rmse(){return util_funcs::norm(_test_error.data(), _n_samp_test) / std::sqrt(static_cast<double>(_n_samp_test));}
+    inline int n_samp_train(){return _n_samp_train;}
+    inline int n_samp_test(){return _n_samp_test;}
+    inline int n_dim(){return _n_dim;}
 
     /**
      * @brief The max Absolute error of the array
@@ -99,36 +103,35 @@ public:
         return std::abs(*std::max_element(_test_error.data(), _test_error.data() + _n_samp_test, [](double d1, double d2){return std::abs(d1) < std::abs(d2);}));
     }
 
-    inline python::list coefs()
-    {
-        python::list coef_lst;
-        for(auto& task_coefs : _coefs)
-            coef_lst.append<python::list>(python_conv_utils::to_list<double>(task_coefs));
-        return coef_lst;
-    }
-
-    inline python::list feats()
-    {
-        python::list feat_lst;
-        for(auto& feat : _feats)
-            feat_lst.append<ModelNode>(*feat);
-        return feat_lst;
-    }
-
-    inline np::ndarray prop_train_est(){return python_conv_utils::to_ndarray<double>(_prop_train_est);}
-    inline np::ndarray prop_test_est(){return python_conv_utils::to_ndarray<double>(_prop_test_est);}
-    inline np::ndarray prop_train(){return python_conv_utils::to_ndarray<double>(_prop_train);}
-    inline np::ndarray prop_test(){return python_conv_utils::to_ndarray<double>(_prop_test);}
-    inline np::ndarray train_error(){return python_conv_utils::to_ndarray<double>(_train_error);}
-    inline np::ndarray test_error(){return python_conv_utils::to_ndarray<double>(_test_error);}
-
     /**
      * @brief Print model to a file
      */
     void to_file(std::string filename, bool train = true, std::vector<int> test_inds = {});
 
-
-    static void register_python();
+    #ifdef PY_BINDINGS
+        inline py::list coefs()
+        {
+            py::list coef_lst;
+            for(auto& task_coefs : _coefs)
+                coef_lst.append<py::list>(python_conv_utils::to_list<double>(task_coefs));
+            return coef_lst;
+        }
+
+        inline py::list feats()
+        {
+            py::list feat_lst;
+            for(auto& feat : _feats)
+                feat_lst.append<ModelNode>(*feat);
+            return feat_lst;
+        }
+
+        inline np::ndarray prop_train_est(){return python_conv_utils::to_ndarray<double>(_prop_train_est);}
+        inline np::ndarray prop_test_est(){return python_conv_utils::to_ndarray<double>(_prop_test_est);}
+        inline np::ndarray prop_train(){return python_conv_utils::to_ndarray<double>(_prop_train);}
+        inline np::ndarray prop_test(){return python_conv_utils::to_ndarray<double>(_prop_test);}
+        inline np::ndarray train_error(){return python_conv_utils::to_ndarray<double>(_train_error);}
+        inline np::ndarray test_error(){return python_conv_utils::to_ndarray<double>(_test_error);}
+    #endif
 };
 
 /**
diff --git a/src/descriptor_identifier/SISSORegressor.cpp b/src/descriptor_identifier/SISSORegressor.cpp
index 189ff0ffaf58284db7007ebc9978f88d3bb65e42..298a5267a0e4d25b118974488b27b36c28fbda82 100644
--- a/src/descriptor_identifier/SISSORegressor.cpp
+++ b/src/descriptor_identifier/SISSORegressor.cpp
@@ -1,6 +1,5 @@
 #include <descriptor_identifier/SISSORegressor.hpp>
 
-
 SISSORegressor::SISSORegressor(
     std::shared_ptr<FeatureSpace> feat_space,
     std::vector<double> prop,
@@ -39,82 +38,6 @@ SISSORegressor::SISSORegressor(
     _work = std::vector<double>(_lwork, 0.0);
 }
 
-SISSORegressor::SISSORegressor(
-    std::shared_ptr<FeatureSpace> feat_space,
-    np::ndarray prop,
-    np::ndarray prop_test,
-    python::list task_sizes_train,
-    python::list task_sizes_test,
-    python::list leave_out_inds,
-    int n_dim,
-    int n_residual
-) :
-    _prop(python_conv_utils::from_ndarray<double>(prop)),
-    _prop_test(python_conv_utils::from_ndarray<double>(prop_test)),
-    _a((n_dim + 1) * prop.shape(0)),
-    _b(prop.shape(0)),
-    _error(prop.shape(0), 0.0),
-    _s(n_dim + 1),
-    _task_sizes_train(python_conv_utils::from_list<int>(task_sizes_train)),
-    _task_sizes_test(python_conv_utils::from_list<int>(task_sizes_test)),
-    _leave_out_inds(python_conv_utils::from_list<int>(leave_out_inds)),
-    _feat_space(feat_space),
-    _mpi_comm(feat_space->mpi_comm()),
-    _n_samp(prop.shape(0)),
-    _n_dim(n_dim),
-    _n_residual(n_residual),
-    _lwork(-1),
-    _rank(0)
-{
-    node_value_arrs::initialize_d_matrix_arr();
-    // Initialize a, b, ones, s, and _error arrays
-    std::fill_n(_a.data(), (_n_dim + 1) * _n_samp, 0.0);
-    std::fill_n(_b.data(), _n_samp, 0.0);
-    std::fill_n(_s.data(), _n_dim + 1, 0.0);
-
-    // // Get optimal work size
-    _lwork = get_opt_lwork(_n_dim + 1);
-    _work = std::vector<double>(_lwork, 0.0);
-}
-
-SISSORegressor::SISSORegressor(
-    std::shared_ptr<FeatureSpace> feat_space,
-    python::list prop,
-    python::list prop_test,
-    python::list task_sizes_train,
-    python::list task_sizes_test,
-    python::list leave_out_inds,
-    int n_dim,
-    int n_residual
-) :
-    _prop(python_conv_utils::from_list<double>(prop)),
-    _prop_test(python_conv_utils::from_list<double>(prop_test)),
-    _a((n_dim + 1) * boost::python::len(prop)),
-    _b(boost::python::len(prop)),
-    _error(boost::python::len(prop), 0.0),
-    _s(n_dim + 1),
-    _task_sizes_train(python_conv_utils::from_list<int>(task_sizes_train)),
-    _task_sizes_test(python_conv_utils::from_list<int>(task_sizes_test)),
-    _leave_out_inds(python_conv_utils::from_list<int>(leave_out_inds)),
-    _feat_space(feat_space),
-    _mpi_comm(feat_space->mpi_comm()),
-    _n_samp(boost::python::len(prop)),
-    _n_dim(n_dim),
-    _n_residual(n_residual),
-    _lwork(-1),
-    _rank(0)
-{
-    node_value_arrs::initialize_d_matrix_arr();
-    // Initialize a, b, ones, s, and _error arrays
-    std::fill_n(_a.data(), (_n_dim + 1) * _n_samp, 0.0);
-    std::fill_n(_b.data(), _n_samp, 0.0);
-    std::fill_n(_s.data(), _n_dim + 1, 0.0);
-
-    // // Get optimal work size
-    _lwork = get_opt_lwork(_n_dim + 1);
-    _work = std::vector<double>(_lwork, 0.0);
-}
-
 void  SISSORegressor::set_a(std::vector<int>& inds, int start, int n_samp)
 {
     for(int ii = 0; ii < inds.size(); ++ii)
@@ -305,31 +228,3 @@ void SISSORegressor::l0_norm(std::vector<double>& prop, int n_dim)
 
     _models.push_back(models);
 }
-
-python::list SISSORegressor::models_py()
-{
-    python::list model_list;
-    for(auto& m_list : _models)
-        model_list.append<python::list>(python_conv_utils::to_list<Model>(m_list));
-
-    return model_list;
-}
-
-void SISSORegressor::register_python()
-{
-    using namespace boost::python;
-    class_<SISSORegressor>("SISSORegressor", init<std::shared_ptr<FeatureSpace>, np::ndarray, np::ndarray, python::list, python::list, python::list, int, int>())
-        .def(init<std::shared_ptr<FeatureSpace>, python::list, python::list, python::list, python::list, python::list, int, int>())
-        .def("fit", &SISSORegressor::fit)
-        .add_property("prop", &SISSORegressor::prop_py)
-        .add_property("prop_test", &SISSORegressor::prop_test_py)
-        .add_property("models", &SISSORegressor::models_py)
-        .add_property("n_samp", &SISSORegressor::n_samp)
-        .add_property("n_dim", &SISSORegressor::n_dim)
-        .add_property("n_residual", &SISSORegressor::n_residual)
-        .add_property("feat_space", &SISSORegressor::feat_space)
-        .add_property("error", &SISSORegressor::error_py)
-        .add_property("task_sizes_train", &SISSORegressor::task_sizes_train)
-        .add_property("task_sizes_test", &SISSORegressor::task_sizes_test)
-    ;
-}
\ No newline at end of file
diff --git a/src/descriptor_identifier/SISSORegressor.hpp b/src/descriptor_identifier/SISSORegressor.hpp
index 177588d8fd0527968f0e02bcef1700ddc0da1528..6f07967d33e9d93b27a54502faa4927d674f8ae2 100644
--- a/src/descriptor_identifier/SISSORegressor.hpp
+++ b/src/descriptor_identifier/SISSORegressor.hpp
@@ -5,8 +5,11 @@
 #include <descriptor_identifier/Model/Model.hpp>
 #include <ctime>
 
-namespace python = boost::python;
-namespace np = boost::python::numpy;
+#ifdef PY_BINDINGS
+    namespace np = boost::python::numpy;
+    namespace py = boost::python;
+#endif
+
 /**
  * @brief SISSO Regressor class, to find the best models, and store them
  *
@@ -56,28 +59,6 @@ public:
         int n_dim,
         int n_residual);
 
-    SISSORegressor(
-        std::shared_ptr<FeatureSpace> feat_space,
-        np::ndarray prop,
-        np::ndarray prop_test,
-        python::list task_sizes_train,
-        python::list task_sizes_test,
-        python::list leave_out_inds,
-        int n_dim,
-        int n_residual
-    );
-
-    SISSORegressor(
-        std::shared_ptr<FeatureSpace> feat_space,
-        python::list prop,
-        python::list prop_test,
-        python::list task_sizes_train,
-        python::list task_sizes_test,
-        python::list leave_out_inds,
-        int n_dim,
-        int n_residual
-    );
-
     /**
      * @brief Get the optimal size of the working array
      *
@@ -132,22 +113,21 @@ public:
      */
     inline std::vector<double> prop(){return _prop;}
 
-    inline np::ndarray prop_py(){return python_conv_utils::to_ndarray<double>(_prop);}
-
     /**
      * @brief Acessor function for prop
      */
     inline std::vector<double> prop_test(){return _prop_test;}
 
-    inline np::ndarray prop_test_py(){return python_conv_utils::to_ndarray<double>(_prop_test);}
+    /**
+     * @brief Acessor function for {
+     */
+    inline std::vector<double> error(){return _error;}
 
     /**
      * @brief Acessor function for models
      */
     inline std::vector<std::vector<Model>> models(){return _models;}
 
-    python::list models_py();
-
     /**
      * @brief Acessor function for n_samp
      */
@@ -163,17 +143,37 @@ public:
      */
     inline int n_residual(){return _n_residual;}
 
-    inline python::list task_sizes_train(){python_conv_utils::to_list<int>(_task_sizes_train);}
-    inline python::list task_sizes_test(){python_conv_utils::to_list<int>(_task_sizes_test);}
-
-    /**
-     * @brief Acessor function for {
-     */
-    inline std::vector<double> error(){return _error;}
-
-    inline np::ndarray error_py(){return python_conv_utils::to_ndarray<double>(_error);}
-
-    static void register_python();
+    // Python interface functions
+    #ifdef PY_BINDINGS
+        SISSORegressor(
+            std::shared_ptr<FeatureSpace> feat_space,
+            np::ndarray prop,
+            np::ndarray prop_test,
+            py::list task_sizes_train,
+            py::list task_sizes_test,
+            py::list leave_out_inds,
+            int n_dim,
+            int n_residual
+        );
+
+        SISSORegressor(
+            std::shared_ptr<FeatureSpace> feat_space,
+            py::list prop,
+            py::list prop_test,
+            py::list task_sizes_train,
+            py::list task_sizes_test,
+            py::list leave_out_inds,
+            int n_dim,
+            int n_residual
+        );
+
+        py::list models_py();
+        inline np::ndarray prop_py(){return python_conv_utils::to_ndarray<double>(_prop);}
+        inline np::ndarray prop_test_py(){return python_conv_utils::to_ndarray<double>(_prop_test);}
+        inline py::list task_sizes_train(){python_conv_utils::to_list<int>(_task_sizes_train);}
+        inline py::list task_sizes_test(){python_conv_utils::to_list<int>(_task_sizes_test);}
+        inline np::ndarray error_py(){return python_conv_utils::to_ndarray<double>(_error);}
+    #endif
 };
 
 #endif
diff --git a/src/feature_creation/feature_space/FeatureSpace.cpp b/src/feature_creation/feature_space/FeatureSpace.cpp
index 311dd69f027c6e12c44bc8a19f56c538ce30d29f..7b0929b2984d35466ca9d8eaec916395398df12b 100644
--- a/src/feature_creation/feature_space/FeatureSpace.cpp
+++ b/src/feature_creation/feature_space/FeatureSpace.cpp
@@ -1,4 +1,5 @@
 #include <feature_creation/feature_space/FeatureSpace.hpp>
+
 BOOST_CLASS_EXPORT_GUID(FeatureNode, "FeatureNode")
 BOOST_CLASS_EXPORT_GUID(AddNode, "AddNode")
 BOOST_CLASS_EXPORT_GUID(SubNode, "SubNode")
@@ -52,86 +53,6 @@ FeatureSpace::FeatureSpace(
     initialize_fs(prop);
 }
 
-FeatureSpace::FeatureSpace(
-    python::list phi_0,
-    python::list allowed_ops,
-    python::list prop,
-    python::list task_sizes,
-    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
-):
-    _phi(python_conv_utils::shared_ptr_vec_from_list<Node, FeatureNode>(phi_0)),
-    _phi_0(_phi),
-    _allowed_ops(python_conv_utils::from_list<std::string>(allowed_ops)),
-    _scores(python::len(phi_0), 0.0),
-    _task_sizes(python_conv_utils::from_list<int>(task_sizes)),
-    _start_gen(1, 0),
-    _feature_space_file("feature_space/selected_features.txt"),
-    _mpi_comm(mpi_setup::comm),
-    _l_bound(min_abs_feat_val),
-    _u_bound(max_abs_feat_val),
-    _max_phi(max_phi),
-    _n_sis_select(n_sis_select),
-    _n_feat(python::len(phi_0)),
-    _n_rung_store(max_store_rung),
-    _n_rung_generate(n_rung_generate)
-{
-    _n_samp = _phi_0[0]->n_samp();
-    initialize_fs(python_conv_utils::from_list<double>(prop));
-}
-
-FeatureSpace::FeatureSpace(
-    python::list phi_0,
-    python::list allowed_ops,
-    np::ndarray prop,
-    python::list task_sizes,
-    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
-):
-    _phi(python_conv_utils::shared_ptr_vec_from_list<Node, FeatureNode>(phi_0)),
-    _phi_0(_phi),
-    _allowed_ops(python_conv_utils::from_list<std::string>(allowed_ops)),
-    _scores(python::len(phi_0), 0.0),
-    _task_sizes(python_conv_utils::from_list<int>(task_sizes)),
-    _start_gen(1, 0),
-    _feature_space_file("feature_space/selected_features.txt"),
-    _mpi_comm(mpi_setup::comm),
-    _l_bound(min_abs_feat_val),
-    _u_bound(max_abs_feat_val),
-    _max_phi(max_phi),
-    _n_sis_select(n_sis_select),
-    _n_feat(python::len(phi_0)),
-    _n_rung_store(max_store_rung),
-    _n_rung_generate(n_rung_generate)
-{
-    _n_samp = _phi_0[0]->n_samp();
-    initialize_fs(python_conv_utils::from_ndarray<double>(prop));
-}
-
-boost::python::list FeatureSpace::phi0_py()
-{
-    python::list feat_lst;
-    for(auto& feat : _phi_0)
-        feat_lst.append<FeatureNode>(FeatureNode(feat->feat_ind(), feat->expr(), feat->value(), feat->test_value(), feat->unit()));
-    return feat_lst;
-}
-
-boost::python::list FeatureSpace::phi_selected_py()
-{
-    python::list feat_lst;
-    for(auto& feat : _phi_selected)
-        feat_lst.append<ModelNode>(ModelNode(feat->d_mat_ind(), feat->rung(), feat->expr(), feat->value(), feat->test_value(), feat->unit()));
-    return feat_lst;
-}
-
 void FeatureSpace::initialize_fs(std::vector<double> prop)
 {
     if(_n_rung_store == -1)
@@ -811,32 +732,3 @@ void FeatureSpace::sis(std::vector<double>& prop)
     if(_mpi_comm->rank() == 0)
         out_file_stream.close();
 }
-
-void FeatureSpace::register_python()
-{
-    using namespace boost::python;
-    void (FeatureSpace::*sis_list)(list) = &FeatureSpace::sis;
-    void (FeatureSpace::*sis_ndarray)(np::ndarray) = &FeatureSpace::sis;
-
-    class_<FeatureSpace>("FeatureSpace", init<list, list, np::ndarray, list, optional<int, int, int, int, double, double>>())
-        .def(init<list, list, list, list, optional<int, int, int, int, double, double>>())
-        .def("sis", sis_list)
-        .def("sis", sis_ndarray)
-        .def("feat_in_phi", &FeatureSpace::feat_in_phi)
-        .add_property("phi_selected", &FeatureSpace::phi_selected_py)
-        .add_property("phi0", &FeatureSpace::phi0_py)
-        .add_property("scores", &FeatureSpace::scores_py)
-        .add_property("task_sizes", &FeatureSpace::task_sizes_py)
-        .add_property("allowed_ops", &FeatureSpace::allowed_ops_py)
-        .add_property("start_gen", &FeatureSpace::start_gen_py)
-        .def_readonly("_feature_space_file", &FeatureSpace::_feature_space_file)
-        .def_readonly("_l_bound", &FeatureSpace::_l_bound)
-        .def_readonly("_u_bound", &FeatureSpace::_u_bound)
-        .def_readonly("_max_phi", &FeatureSpace::_max_phi)
-        .def_readonly("_n_sis_select", &FeatureSpace::_n_sis_select)
-        .def_readonly("_n_samp", &FeatureSpace::_n_samp)
-        .def_readonly("_n_feat", &FeatureSpace::_n_feat)
-        .def_readonly("_n_rung_store", &FeatureSpace::_n_rung_store)
-        .def_readonly("_n_rung_generate", &FeatureSpace::_n_rung_generate)
-    ;
-}
diff --git a/src/feature_creation/feature_space/FeatureSpace.hpp b/src/feature_creation/feature_space/FeatureSpace.hpp
index 3b727b07f4878342c490c57221d21c74e855b8b5..f9e568ab209cbe39cd5ea0694921a45c25bdea97 100644
--- a/src/feature_creation/feature_space/FeatureSpace.hpp
+++ b/src/feature_creation/feature_space/FeatureSpace.hpp
@@ -10,13 +10,14 @@
 
 #include <boost/serialization/shared_ptr.hpp>
 #include <boost/filesystem.hpp>
-#include <boost/python.hpp>
 
 #include <iostream>
 #include <iomanip>
 
-namespace python = boost::python;
-namespace np = boost::python::numpy;
+#ifdef PY_BINDINGS
+    namespace np = boost::python::numpy;
+    namespace py = boost::python;
+#endif
 
 /**
  * @brief Feature Space for SISSO calculations
@@ -79,52 +80,6 @@ public:
         double max_abs_feat_val=1e50
     );
 
-    /**
-     * @brief Constructor for the feature space
-     * @details constructs the feature space from an initial set of features and a list of allowed operatiors
-     *
-     * @param mpi_comm MPI communicator for the calculations
-     * @param allowed_ops list of allowed operators
-     * @param max_phi highest rung value for the calculation
-     * @param n_sis_select number of features to select during each SIS step
-     * @param max_abs_feat_val maximum absolute feature value
-     */
-    FeatureSpace(
-        python::list phi_0,
-        python::list allowed_ops,
-        python::list prop,
-        python::list task_sizes,
-        int max_phi=1,
-        int n_sis_select=1,
-        int max_store_rung=-1,
-        int n_rung_generate=0,
-        double min_abs_feat_val=1e-50,
-        double max_abs_feat_val=1e50
-    );
-
-    /**
-     * @brief Constructor for the feature space
-     * @details constructs the feature space from an initial set of features and a list of allowed operatiors
-     *
-     * @param mpi_comm MPI communicator for the calculations
-     * @param allowed_ops list of allowed operators
-     * @param max_phi highest rung value for the calculation
-     * @param n_sis_select number of features to select during each SIS step
-     * @param max_abs_feat_val maximum absolute feature value
-     */
-    FeatureSpace(
-        python::list phi_0,
-        python::list allowed_ops,
-        np::ndarray prop,
-        python::list task_sizes,
-        int max_phi=1,
-        int n_sis_select=1,
-        int max_store_rung=-1,
-        int n_rung_generate=0,
-        double min_abs_feat_val=1e-50,
-        double max_abs_feat_val=1e50
-    );
-
     void initialize_fs(std::vector<double> prop);
 
     /**
@@ -138,8 +93,6 @@ public:
      */
     inline std::vector<node_ptr> phi_selected(){return _phi_selected;};
 
-    boost::python::list phi_selected_py();
-
     /**
      * @brief Accessor function for _phi
      */
@@ -150,14 +103,11 @@ public:
      */
     inline std::vector<node_ptr> phi0(){return _phi_0;};
 
-    boost::python::list phi0_py();
     /**
      * @brief Accessor function for _scores
      */
     inline std::vector<double> scores(){return _scores;}
 
-    inline np::ndarray scores_py(){return python_conv_utils::to_ndarray<double>(_scores);};
-
     /**
      * @brief Accessor function for _mpi_comm
      */
@@ -165,11 +115,15 @@ public:
 
     inline std::vector<int> task_sizes(){return _task_sizes;}
 
-    inline boost::python::list task_sizes_py(){return python_conv_utils::to_list<int>(_task_sizes);};
-
-    inline boost::python::list allowed_ops_py(){return python_conv_utils::to_list<std::string>(_allowed_ops);}
-
-    inline boost::python::list start_gen_py(){return python_conv_utils::to_list<int>(_start_gen);}
+    inline std::string feature_space_file(){return _feature_space_file;}
+    inline double l_bound(){return _l_bound;}
+    inline double u_bound(){return _u_bound;}
+    inline int max_phi(){return _max_phi;}
+    inline int n_sis_select(){return _n_sis_select;}
+    inline int n_samp(){return _n_samp;}
+    inline int n_feat(){return _n_feat;}
+    inline int n_rung_store(){return _n_rung_store;}
+    inline int n_rung_generate(){return _n_rung_generate;}
 
     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);
 
@@ -186,18 +140,6 @@ public:
      */
     void sis(std::vector<double>& prop);
 
-    inline void sis(np::ndarray prop)
-    {
-        std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
-        sis(prop_vec);
-    }
-
-    inline void sis(python::list prop)
-    {
-        std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
-        sis(prop_vec);
-    }
-
     /**
      * @brief Is a feature in this process' _phi?
      *
@@ -206,7 +148,74 @@ public:
      */
     inline bool feat_in_phi(int ind){return (ind >= _phi[0]->feat_ind()) && (ind <= _phi.back()->feat_ind());}
 
-    static void register_python();
+    // Python Interface Functions
+    #ifdef PY_BINDINGS
+        /**
+         * @brief Constructor for the feature space
+         * @details constructs the feature space from an initial set of features and a list of allowed operatiors
+         *
+         * @param mpi_comm MPI communicator for the calculations
+         * @param allowed_ops list of allowed operators
+         * @param max_phi highest rung value for the calculation
+         * @param n_sis_select number of features to select during each SIS step
+         * @param max_abs_feat_val maximum absolute feature value
+         */
+        FeatureSpace(
+            py::list phi_0,
+            py::list allowed_ops,
+            py::list prop,
+            py::list task_sizes,
+            int max_phi=1,
+            int n_sis_select=1,
+            int max_store_rung=-1,
+            int n_rung_generate=0,
+            double min_abs_feat_val=1e-50,
+            double max_abs_feat_val=1e50
+        );
+
+        /**
+         * @brief Constructor for the feature space
+         * @details constructs the feature space from an initial set of features and a list of allowed operatiors
+         *
+         * @param mpi_comm MPI communicator for the calculations
+         * @param allowed_ops list of allowed operators
+         * @param max_phi highest rung value for the calculation
+         * @param n_sis_select number of features to select during each SIS step
+         * @param max_abs_feat_val maximum absolute feature value
+         */
+        FeatureSpace(
+            py::list phi_0,
+            py::list allowed_ops,
+            np::ndarray prop,
+            py::list task_sizes,
+            int max_phi=1,
+            int n_sis_select=1,
+            int max_store_rung=-1,
+            int n_rung_generate=0,
+            double min_abs_feat_val=1e-50,
+            double max_abs_feat_val=1e50
+        );
+
+
+        inline void sis(np::ndarray prop)
+        {
+            std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
+            sis(prop_vec);
+        }
+
+        inline void sis(py::list prop)
+        {
+            std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
+            sis(prop_vec);
+        }
+
+        py::list phi_selected_py();
+        py::list phi0_py();
+        inline np::ndarray scores_py(){return python_conv_utils::to_ndarray<double>(_scores);};
+        inline py::list task_sizes_py(){return python_conv_utils::to_list<int>(_task_sizes);};
+        inline py::list allowed_ops_py(){return python_conv_utils::to_list<std::string>(_allowed_ops);}
+        inline py::list start_gen_py(){return python_conv_utils::to_list<int>(_start_gen);}
+    #endif
 };
 
 #endif
\ No newline at end of file
diff --git a/src/feature_creation/node/FeatureNode.cpp b/src/feature_creation/node/FeatureNode.cpp
index 6184a14100a565d8a189181a1c5bc4ebc5f35f1e..9b0031d4dc42a0609f98521a8e7c2d293d2236ff 100644
--- a/src/feature_creation/node/FeatureNode.cpp
+++ b/src/feature_creation/node/FeatureNode.cpp
@@ -17,45 +17,6 @@ FeatureNode::FeatureNode(int feat_ind, std::string expr, std::vector<double> val
     }
 }
 
-FeatureNode::FeatureNode(int feat_ind, std::string expr, np::ndarray value, np::ndarray test_value, Unit unit) :
-    Node(feat_ind, value.shape(0), test_value.shape(0)),
-    _value(python_conv_utils::from_ndarray<double>(value)),
-    _test_value(python_conv_utils::from_ndarray<double>(test_value)),
-    _unit(unit),
-    _expr(expr)
-{
-    // Automatically resize the storage arrays
-    if(node_value_arrs::N_STORE_FEATURES == 0)
-        node_value_arrs::initialize_values_arr(_n_samp, _n_test_samp, 1);
-    else if((_n_samp != node_value_arrs::N_SAMPLES) || (_n_test_samp != node_value_arrs::N_SAMPLES_TEST))
-        throw std::logic_error("Number of samples in current feature is not the same as the others, (" + std::to_string(_n_samp) + " and " + std::to_string(_n_test_samp) + " vs. "  + std::to_string(node_value_arrs::N_SAMPLES) + " and " + std::to_string(node_value_arrs::N_SAMPLES_TEST) + ")");
-    else if(feat_ind >= node_value_arrs::N_STORE_FEATURES)
-        node_value_arrs::resize_values_arr(0, node_value_arrs::N_STORE_FEATURES + 1, true);
-
-    set_value();
-    set_test_value();
-}
-
-FeatureNode::FeatureNode(int feat_ind, std::string expr, python::list value, python::list test_value, Unit unit) :
-    Node(feat_ind, python::len(value), python::len(test_value)),
-    _value(python_conv_utils::from_list<double>(value)),
-    _test_value(python_conv_utils::from_list<double>(test_value)),
-    _unit(unit),
-    _expr(expr)
-{
-
-    // Automatically resize the storage arrays
-    if(node_value_arrs::N_STORE_FEATURES == 0)
-        node_value_arrs::initialize_values_arr(_n_samp, _n_test_samp, 1);
-    else if((_n_samp != node_value_arrs::N_SAMPLES) || (_n_test_samp != node_value_arrs::N_SAMPLES_TEST))
-        throw std::logic_error("Number of samples in current feature is not the same as the others, (" + std::to_string(_n_samp) + " and " + std::to_string(_n_test_samp) + " vs. "  + std::to_string(node_value_arrs::N_SAMPLES) + " and " + std::to_string(node_value_arrs::N_SAMPLES_TEST) + ")");
-    else if(feat_ind >= node_value_arrs::N_STORE_FEATURES)
-        node_value_arrs::resize_values_arr(0, node_value_arrs::N_STORE_FEATURES + 1, true);
-
-    set_value();
-    set_test_value();
-}
-
 FeatureNode::~FeatureNode()
 {}
 
@@ -79,22 +40,4 @@ void FeatureNode::update_div_mult_leaves(std::map<std::string, double>& div_mult
     expected_abs_tot += std::abs(fact);
 }
 
-void FeatureNode::register_python()
-{
-    std::string (FeatureNode::*expr_1)() = &FeatureNode::expr;
-    std::string (FeatureNode::*expr_const)() const = &FeatureNode::expr;
-
-    using namespace boost::python;
-    class_<FeatureNode, bases<Node>>("FeatureNode", init<int, std::string, np::ndarray, np::ndarray, Unit>())
-        .def(init<int, std::string, python::list, python::list, Unit>())
-        .def("is_nan", &FeatureNode::is_nan)
-        .def("is_const", &FeatureNode::is_const)
-        .def("set_value", &FeatureNode::set_value)
-        .def("set_test_value", &FeatureNode::set_test_value)
-        .add_property("expr", expr_1)
-        .add_property("expr", expr_const)
-        .add_property("unit", &FeatureNode::unit)
-        .add_property("rung", &FeatureNode::rung)
-    ;
-}
 // BOOST_CLASS_EXPORT(FeatureNode)
diff --git a/src/feature_creation/node/FeatureNode.hpp b/src/feature_creation/node/FeatureNode.hpp
index 7979fdee72f513412c4a31ee5756971f3b4a6f75..8c3a14e4af6bfac6eb6f9744882c88aa746b48bb 100644
--- a/src/feature_creation/node/FeatureNode.hpp
+++ b/src/feature_creation/node/FeatureNode.hpp
@@ -1,7 +1,6 @@
 #ifndef FEATURE_NODE
 #define FEATURE_NODE
 
-#include <python/conversion_utils.hpp>
 #include <utils/math_funcs.hpp>
 #include <utils/enum.hpp>
 
@@ -9,12 +8,13 @@
 
 #include <boost/serialization/export.hpp>
 #include <boost/serialization/base_object.hpp>
-#include <boost/python/numpy.hpp>
 
 #include <feature_creation/node/Node.hpp>
 
-namespace np = boost::python::numpy;
-namespace python = boost::python;
+#ifdef PY_BINDINGS
+    namespace np = boost::python::numpy;
+    namespace py = boost::python;
+#endif
 
 /**
  * @brief Node that describe the leaves of the operator graph (Initial features in Phi_0)
@@ -61,9 +61,10 @@ public:
      * @param unit Unit of the feature
      */
     FeatureNode(int feat_ind, std::string expr, std::vector<double> value, std::vector<double> test_value, Unit unit, bool set_val = true);
-    FeatureNode(int feat_ind, std::string expr, np::ndarray value, np::ndarray test_value, Unit unit);
-    FeatureNode(int feat_ind, std::string expr, python::list value, python::list test_value, Unit unit);
-
+    #ifdef PY_BINDINGS
+        FeatureNode(int feat_ind, std::string expr, np::ndarray value, np::ndarray test_value, Unit unit);
+        FeatureNode(int feat_ind, std::string expr, py::list value, py::list test_value, Unit unit);
+    #endif
     FeatureNode(const FeatureNode&) = default;
     FeatureNode(FeatureNode&&) = default;
 
@@ -154,8 +155,6 @@ public:
      * @param add_sub_leaves [description]
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
diff --git a/src/feature_creation/node/ModelNode.cpp b/src/feature_creation/node/ModelNode.cpp
index 246751e88afd4f9b59b65fc27a14103dba4e61b5..13b962916bd0b8fb6636b487f38f9de9f238928d 100644
--- a/src/feature_creation/node/ModelNode.cpp
+++ b/src/feature_creation/node/ModelNode.cpp
@@ -30,18 +30,3 @@ void ModelNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_l
 
     expected_abs_tot += std::abs(fact);
 }
-
-void ModelNode::register_python()
-{
-    std::string (ModelNode::*expr_1)() = &ModelNode::expr;
-    std::string (ModelNode::*expr_const)() const = &ModelNode::expr;
-
-    using namespace boost::python;
-    class_<ModelNode, bases<FeatureNode>>("ModelNode", init<int, int, std::string, std::vector<double>, std::vector<double>, Unit>())
-        .def("is_nan", &ModelNode::is_nan)
-        .def("is_const", &ModelNode::is_const)
-        .def("set_value", &ModelNode::set_value)
-        .def("set_test_value", &ModelNode::set_test_value)
-        .add_property("rung", &ModelNode::rung)
-    ;
-}
diff --git a/src/feature_creation/node/ModelNode.hpp b/src/feature_creation/node/ModelNode.hpp
index c945fa4efade4e3c6e7ef1a9ba1ccd3c3ce74a53..d80a2f8ee09e9d784edb054c30afa26266f8681d 100644
--- a/src/feature_creation/node/ModelNode.hpp
+++ b/src/feature_creation/node/ModelNode.hpp
@@ -108,8 +108,6 @@ public:
      * @param add_sub_leaves [description]
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
diff --git a/src/feature_creation/node/Node.cpp b/src/feature_creation/node/Node.cpp
index 851b6b07e50bc60b88076e8f3eb8ce5d9ab86393..22432f60db9b30290c3e5c6317b2f5b5a0cca688 100644
--- a/src/feature_creation/node/Node.cpp
+++ b/src/feature_creation/node/Node.cpp
@@ -15,107 +15,5 @@ Node::Node(int feat_ind, int n_samp, int n_test_samp) :
 Node::~Node()
 {}
 
-using namespace boost::python;
-struct NodeWrap : Node, wrapper<Node>
-{
-    std::string expr()
-    {
-        return this->get_override("expr")();
-    }
-
-    Unit unit()
-    {
-        return this->get_override("unit")();
-    }
-
-    std::vector<double> value()
-    {
-        return this->get_override("value")();
-    }
-
-    std::vector<double> test_value()
-    {
-        return this->get_override("test_value")();
-    }
-
-    void set_value(int offset = -1)
-    {
-        this->get_override("set_value")();
-    }
-
-    double* value_ptr(int offset = -1)
-    {
-        return this->get_override("value_ptr")();
-    }
-
-    void set_test_value(int offset = -1)
-    {
-        this->get_override("set_test_value")();
-    }
-
-    double* test_value_ptr(int offset = -1)
-    {
-        return this->get_override("test_value_ptr")();
-    }
-
-    bool is_nan()
-    {
-        return this->get_override("is_nan")();
-    }
-
-    bool is_const()
-    {
-        return this->get_override("is_const")();
-    }
-
-    NODE_TYPE type()
-    {
-        return this->get_override("type")();
-    }
-
-    int rung(int cur_rung = 0)
-    {
-        return this->get_override("rung")();
-    }
-
-    void update_add_sub_leaves(std::map<std::string, int>& add_sub_leaves, int pl_mn, int& expected_abs_tot)
-    {
-        this->get_override("update_add_sub_leaves");
-    }
-
-    void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot)
-    {
-        this->get_override("update_div_mult_leaves");
-    }
-
-};
-
-void Node::register_python()
-{
-    void (Node::*reindex_1)(int) = &Node::reindex;
-    void (Node::*reindex_2)(int, int) = &Node::reindex;
-    class_<NodeWrap, boost::noncopyable>("Node", no_init)
-        .def("reindex", reindex_1)
-        .def("reindex", reindex_2)
-        .def("__str__", &Node::expr)
-        .def("__repr__", &Node::expr)
-        .add_property("n_samp", &Node::n_samp)
-        .add_property("n_test_samp", &Node::n_test_samp)
-        .add_property("feat_ind", &Node::feat_ind)
-        .add_property("arr_ind", &Node::arr_ind)
-        .add_property("selected", &Node::selected, &Node::set_selected)
-        .add_property("d_mat_ind", &Node::d_mat_ind, &Node::set_d_mat_ind)
-        .add_property("value", &Node::value_py)
-        .add_property("test_value", &Node::test_value_py)
-        .def("expr", pure_virtual(&Node::expr))
-        .def("unit", pure_virtual(&Node::unit))
-        .def("set_value", pure_virtual(&Node::set_value))
-        .def("set_test_value", pure_virtual(&Node::set_test_value))
-        .def("is_nan", pure_virtual(&Node::is_nan))
-        .def("is_const", pure_virtual(&Node::is_const))
-        .def("rung", pure_virtual(&Node::rung))
-    ;
-}
-
 BOOST_SERIALIZATION_ASSUME_ABSTRACT(Node)
 
diff --git a/src/feature_creation/node/Node.hpp b/src/feature_creation/node/Node.hpp
index 3854df4734327a34b192620cc453070cccd15cbf..ecaf7df29c09bf5e32409accf9f4e3fc3628d317 100644
--- a/src/feature_creation/node/Node.hpp
+++ b/src/feature_creation/node/Node.hpp
@@ -3,20 +3,23 @@
 
 #include <feature_creation/units/Unit.hpp>
 #include <feature_creation/node/value_storage/nodes_value_containers.hpp>
-#include <python/conversion_utils.hpp>
+#ifdef PY_BINDINGS
+    #include <python/conversion_utils.hpp>
+#endif
 #include <utils/math_funcs.hpp>
 #include <utils/enum.hpp>
 
 #include <cmath>
 #include <memory>
 
-#include <boost/python/numpy.hpp>
 #include <boost/serialization/assume_abstract.hpp>
 #include <boost/serialization/serialization.hpp>
 #include <boost/serialization/string.hpp>
 #include <boost/serialization/unique_ptr.hpp>
 
-namespace np = boost::python::numpy;
+#ifdef PY_BINDINGS
+    namespace np = boost::python::numpy;
+#endif
 
 /**
  * @brief Base class for a Node
@@ -134,14 +137,12 @@ public:
      */
     virtual std::vector<double> value() = 0;
 
-    inline np::ndarray value_py(){return python_conv_utils::to_ndarray<double>(value());}
 
     /**
      * @brief Get the test_value of the descriptor
      */
     virtual std::vector<double> test_value() = 0;
 
-    inline np::ndarray test_value_py(){return python_conv_utils::to_ndarray<double>(test_value());}
 
     /**
      * @brief Set the value for the feature
@@ -200,8 +201,12 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     virtual void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot) = 0;
+    #ifdef PY_BINDINGS
 
-    static void register_python();
+
+        inline np::ndarray value_py(){return python_conv_utils::to_ndarray<double>(value());}
+        inline np::ndarray test_value_py(){return python_conv_utils::to_ndarray<double>(test_value());}
+    #endif
 };
 
 typedef std::shared_ptr<Node> node_ptr;
diff --git a/src/feature_creation/node/operator_nodes/OperatorNode.hpp b/src/feature_creation/node/operator_nodes/OperatorNode.hpp
index c73f4ab2e74a1f9a1571223c51c816fc51fa73f3..c5b506adc2771c0185fb7b7f96981e2d54a4680d 100644
--- a/src/feature_creation/node/operator_nodes/OperatorNode.hpp
+++ b/src/feature_creation/node/operator_nodes/OperatorNode.hpp
@@ -12,8 +12,6 @@
 #include <boost/serialization/array.hpp>
 #include <boost/serialization/vector.hpp>
 
-using namespace boost::python;
-
 /**
  * @brief Base class to describe operator nodes
  * @details
@@ -166,64 +164,6 @@ public:
      * @param fact amount to increment the dictionary by
      */
     virtual void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot) = 0;
-
-    struct OperatorNodeWrap : OperatorNode<N>, wrapper<OperatorNode<N>>
-    {
-
-        void set_value(int offset = -1)
-        {
-            this->get_override("set_value")();
-        }
-
-        void set_test_value(int offset = -1)
-        {
-            this->get_override("set_test_value")();
-        }
-
-        NODE_TYPE type()
-        {
-            return this->get_override("type")();
-        }
-
-        int rung(int cur_rung = 0)
-        {
-            return this->get_override("rung")();
-        }
-
-        std::string expr()
-        {
-            return this->get_override("expr")();
-        }
-
-        Unit unit()
-        {
-            return this->get_override("unit")();
-        }
-
-        void update_add_sub_leaves(std::map<std::string, int>& add_sub_leaves, int pl_mn, int& expected_abs_tot)
-        {
-            this->get_override("update_add_sub_leaves")();
-        }
-        void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot)
-        {
-            this->get_override("update_div_mult_leaves")();
-        }
-    };
-
-    static void register_python()
-    {
-        // class_<OperatorNode, bases<Node>>("OperatorNode", init<std::array<node_ptr, N>, int>())
-        // ;
-        class_<OperatorNodeWrap, bases<Node>, boost::noncopyable>("OperatorNode")
-            .def("is_nan", &OperatorNode::is_nan)
-            .def("is_const", &OperatorNode::is_const)
-            .def("set_value", pure_virtual(&OperatorNode::set_value))
-            .def("set_test_value", pure_virtual(&OperatorNode::set_test_value))
-            .def("rung", pure_virtual(&OperatorNode::rung))
-            .def("expr", pure_virtual(&OperatorNode::expr))
-            .def("unit", pure_virtual(&OperatorNode::unit))
-        ;
-    }
 };
 
 #endif
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 78e7078067132c7e4ae7172b8e18bdb63a06537f..708e9c52fb1d5941064e27d84cb4b3ec8e69f685 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
@@ -84,15 +84,3 @@ void AbsDiffNode::update_div_mult_leaves(std::map<std::string, double>& div_mult
 
     expected_abs_tot += std::abs(fact);
 }
-
-void AbsDiffNode::register_python()
-{
-    using namespace boost::python;
-    class_<AbsDiffNode, bases<OperatorNode<2>>>("AbsDiffNode", init<node_ptr, node_ptr, int, double, double>())
-        .def("set_value", &AbsDiffNode::set_value)
-        .def("set_test_value", &AbsDiffNode::set_test_value)
-        .add_property("expr", &AbsDiffNode::expr)
-        .add_property("unit", &AbsDiffNode::unit)
-        .add_property("rung", &AbsDiffNode::rung)
-    ;
-}
\ No newline at end of file
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 4eab084620c715a8437b15611a504f5bd1a88b4f..994e80298d9d422f37c8ff3d53aee2b6d1087f80 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
@@ -67,8 +67,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 145e5ac87d6db0e85068fa9ab1dbe0863070f3b9..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
@@ -46,15 +46,3 @@ void AbsNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void AbsNode::register_python()
-{
-    using namespace boost::python;
-    class_<AbsNode, bases<OperatorNode<1>>>("AbsNode", init<node_ptr, int, double, double>())
-        .def("set_value", &AbsNode::set_value)
-        .def("set_test_value", &AbsNode::set_test_value)
-        .add_property("expr", &AbsNode::expr)
-        .add_property("unit", &AbsNode::unit)
-        .add_property("rung", &AbsNode::rung)
-    ;
-}
\ No newline at end of file
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 80479324fbe354c385747cd9fe6010a24a10724e..44d7fc2a189ddcb4126acb7118e0e9ff8fcf71f2 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
@@ -66,7 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 031a1fbc2dd2dd7300a3c6ed6481002cb15e671c..b5d0ff5995275f91ca4ffca8471a7a04dcc05765 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
@@ -75,15 +75,3 @@ void AddNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void AddNode::register_python()
-{
-    using namespace boost::python;
-    class_<AddNode, bases<OperatorNode<2>>>("AddNode", init<node_ptr, node_ptr, int, double, double>())
-        .def("set_value", &AddNode::set_value)
-        .def("set_test_value", &AddNode::set_test_value)
-        .add_property("expr", &AddNode::expr)
-        .add_property("unit", &AddNode::unit)
-        .add_property("rung", &AddNode::rung)
-    ;
-}
\ No newline at end of file
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 04296e3c7595c9a8690d3f7a28cd1994de56d59d..823f3788f3e2cf274b5d806f9a2ff4fce4b30641 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
@@ -66,8 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 d611f00da473c05dca463720b1309b4283852ccd..f6afc5165a92cb5716c21e2dd0d2cdcdc4123603 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
@@ -56,15 +56,3 @@ void CosNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void CosNode::register_python()
-{
-    using namespace boost::python;
-    class_<CosNode, bases<OperatorNode<1>>>("CosNode", init<node_ptr, int, double, double>())
-        .def("set_value", &CosNode::set_value)
-        .def("set_test_value", &CosNode::set_test_value)
-        .add_property("expr", &CosNode::expr)
-        .add_property("unit", &CosNode::unit)
-        .add_property("rung", &CosNode::rung)
-    ;
-}
\ No newline at end of file
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 c701887866e5b6a94adc3d6c4c0837f84c0733bd..110af9629c6e2ecbccec399321a02659e6720735 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
@@ -66,8 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 5a57120bb79de5748f21986c49cae4612ace50a1..6d53cea31bed572543f1a7d9951616cb9b9eae84 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
@@ -44,15 +44,3 @@ void CbNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_leav
 {
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact * 3.0, expected_abs_tot);
 }
-
-void CbNode::register_python()
-{
-    using namespace boost::python;
-    class_<CbNode, bases<OperatorNode<1>>>("CbNode", init<node_ptr, int, double, double>())
-        .def("set_value", &CbNode::set_value)
-        .def("set_test_value", &CbNode::set_test_value)
-        .add_property("expr", &CbNode::expr)
-        .add_property("unit", &CbNode::unit)
-        .add_property("rung", &CbNode::rung)
-    ;
-}
\ No newline at end of file
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 d1e25bb1bbbcc0ab959256e944fc9c16e73b96c6..5dd0071b6ec93e3135ea89a636023c252d501543 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
@@ -66,7 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 b31ef4de2613227409833ad311bc9191712d2857..47f4daff556b126c7fd163a0ed2b42ca5a1cdcdd 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
@@ -44,15 +44,3 @@ void CbrtNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_le
 {
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact / 3.0, expected_abs_tot);
 }
-
-void CbrtNode::register_python()
-{
-    using namespace boost::python;
-    class_<CbrtNode, bases<OperatorNode<1>>>("CbrtNode", init<node_ptr, int, double, double>())
-        .def("set_value", &CbrtNode::set_value)
-        .def("set_test_value", &CbrtNode::set_test_value)
-        .add_property("expr", &CbrtNode::expr)
-        .add_property("unit", &CbrtNode::unit)
-        .add_property("rung", &CbrtNode::rung)
-    ;
-}
\ No newline at end of file
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 5d6ed1621ebcb24ffbe163859ae2a473232c2d3b..4e2344178a96e24366ed7d50271a4fc5d015b381 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
@@ -66,8 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 5b371c2c87891c14d91948f6cb94c4d636f22d75..098e623cf20a60d54bcbfd992c45fe31e6bc05e0 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
@@ -75,15 +75,3 @@ void DivNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact, expected_abs_tot);
     _feats[1]->update_div_mult_leaves(div_mult_leaves, -1.0*fact, expected_abs_tot);
 }
-
-void DivNode::register_python()
-{
-    using namespace boost::python;
-    class_<DivNode, bases<OperatorNode<2>>>("DivNode", init<node_ptr, node_ptr, int, double, double>())
-        .def("set_value", &DivNode::set_value)
-        .def("set_test_value", &DivNode::set_test_value)
-        .add_property("expr", &DivNode::expr)
-        .add_property("unit", &DivNode::unit)
-        .add_property("rung", &DivNode::rung)
-    ;
-}
\ No newline at end of file
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 9a6cabb95b6682d54b136da8c91526abcb6eadf0..7f0af1bb7c053f442c8b7335ae2ff55f8030161c 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
@@ -66,8 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 08a128280509a628bc7c02909c305ad030066484..c78c0e37448fef4e12919e91036a991b89b57748 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
@@ -56,15 +56,3 @@ void ExpNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void ExpNode::register_python()
-{
-    using namespace boost::python;
-    class_<ExpNode, bases<OperatorNode<1>>>("ExpNode", init<node_ptr, int, double, double>())
-        .def("set_value", &ExpNode::set_value)
-        .def("set_test_value", &ExpNode::set_test_value)
-        .add_property("expr", &ExpNode::expr)
-        .add_property("unit", &ExpNode::unit)
-        .add_property("rung", &ExpNode::rung)
-    ;
-}
\ No newline at end of file
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 49f4bdb9bff0c6bbcf786e51ad4d936bb08b6c29..1b6c9eafa044837d90065ae15791526851ce02f9 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
@@ -66,8 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 371e10c31ef898c2359a5fc1abe6eb8bd2de062a..12d2d7d3df7a2713545097cc5378bb08b439320a 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
@@ -44,15 +44,3 @@ void InvNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 {
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact * -1.0, expected_abs_tot);
 }
-
-void InvNode::register_python()
-{
-    using namespace boost::python;
-    class_<InvNode, bases<OperatorNode<1>>>("InvNode", init<node_ptr, int, double, double>())
-        .def("set_value", &InvNode::set_value)
-        .def("set_test_value", &InvNode::set_test_value)
-        .add_property("expr", &InvNode::expr)
-        .add_property("unit", &InvNode::unit)
-        .add_property("rung", &InvNode::rung)
-    ;
-}
\ No newline at end of file
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 9303c75bf015612afa7d6bed6f1e448d558f6c2e..fcd2565d61b19043cc7d5ae099a5a77e2d338f4d 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
@@ -67,7 +67,5 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 #endif
\ No newline at end of file
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 642f8014107d1aa1112037455c87941477ae83b9..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
@@ -56,15 +56,3 @@ void LogNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void LogNode::register_python()
-{
-    using namespace boost::python;
-    class_<LogNode, bases<OperatorNode<1>>>("LogNode", init<node_ptr, int, double, double>())
-        .def("set_value", &LogNode::set_value)
-        .def("set_test_value", &LogNode::set_test_value)
-        .add_property("expr", &LogNode::expr)
-        .add_property("unit", &LogNode::unit)
-        .add_property("rung", &LogNode::rung)
-    ;
-}
\ No newline at end of file
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 af70a775bab2361174229dfd587fb90e8ad3e9d9..81b97ba7722051fd390d40ee5c26d94aee3590c4 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
@@ -66,7 +66,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 36efccdf4445dea36592a0d37630585f5eebeaef..16aed94d8da1652889283d04e320a3c1b917b62d 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
@@ -75,15 +75,3 @@ void MultNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_le
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact, expected_abs_tot);
     _feats[1]->update_div_mult_leaves(div_mult_leaves, fact, expected_abs_tot);
 }
-
-void MultNode::register_python()
-{
-    using namespace boost::python;
-    class_<MultNode, bases<OperatorNode<2>>>("MultNode", init<node_ptr, node_ptr, int, double, double>())
-        .def("set_value", &MultNode::set_value)
-        .def("set_test_value", &MultNode::set_test_value)
-        .add_property("expr", &MultNode::expr)
-        .add_property("unit", &MultNode::unit)
-        .add_property("rung", &MultNode::rung)
-    ;
-}
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 a7be03d55543b635372fd93a5a7b58b363a14044..ef1bb3ce09024805ea53f5cf30b28a19bceea303 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
@@ -67,8 +67,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 e6eccf80e17a248b689b72769616411c9f809ee2..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
@@ -56,15 +56,3 @@ void NegExpNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_
 
     expected_abs_tot += std::abs(fact);
 }
-
-void NegExpNode::register_python()
-{
-    using namespace boost::python;
-    class_<NegExpNode, bases<OperatorNode<1>>>("NegExpNode", init<node_ptr, int, double, double>())
-        .def("set_value", &NegExpNode::set_value)
-        .def("set_test_value", &NegExpNode::set_test_value)
-        .add_property("expr", &NegExpNode::expr)
-        .add_property("unit", &NegExpNode::unit)
-        .add_property("rung", &NegExpNode::rung)
-    ;
-}
\ No newline at end of file
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 664f581b5c73ef1615affd6a4b6b7d15543da3cb..03d48e85f48114146c728d8236551f5a96adc937 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
@@ -67,8 +67,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 c142c1316c1de53b9f3badff35a0941c76186d0b..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
@@ -56,15 +56,3 @@ void SinNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void SinNode::register_python()
-{
-    using namespace boost::python;
-    class_<SinNode, bases<OperatorNode<1>>>("SinNode", init<node_ptr, int, double, double>())
-        .def("set_value", &SinNode::set_value)
-        .def("set_test_value", &SinNode::set_test_value)
-        .add_property("expr", &SinNode::expr)
-        .add_property("unit", &SinNode::unit)
-        .add_property("rung", &SinNode::rung)
-    ;
-}
\ No newline at end of file
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 275366af44b1b0c03cdf539857692f4d2e14ac95..d4d1a77ed647b08a4e91ffb90417c8e9c98a85ab 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
@@ -67,8 +67,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 5a3a0513667bb7100b9a57c974bb0d681a384fb3..2e67b4b7d9bd7cceea7a7d3dfbc6e2ef9d5f696a 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
@@ -44,15 +44,3 @@ void SixPowNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_
 {
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact * 6.0, expected_abs_tot);
 }
-
-void SixPowNode::register_python()
-{
-    using namespace boost::python;
-    class_<SixPowNode, bases<OperatorNode<1>>>("SixPowNode", init<node_ptr, int, double, double>())
-        .def("set_value", &SixPowNode::set_value)
-        .def("set_test_value", &SixPowNode::set_test_value)
-        .add_property("expr", &SixPowNode::expr)
-        .add_property("unit", &SixPowNode::unit)
-        .add_property("rung", &SixPowNode::rung)
-    ;
-}
\ No newline at end of file
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 bad3ff5f7de7d76d3894f4122c2854da7e900c23..770af9f940e6665cbdf485a442fb7a6fd535bea8 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
@@ -68,7 +68,6 @@ public:
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
 
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 09f811777da2c88dac1c0497ff176754520ef26d..f85e1f0bb547f5290304b0cdfb7db4b4741c5f3e 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
@@ -44,15 +44,3 @@ void SqNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_leav
 {
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact * 2.0, expected_abs_tot);
 }
-
-void SqNode::register_python()
-{
-    using namespace boost::python;
-    class_<SqNode, bases<OperatorNode<1>>>("SqNode", init<node_ptr, int, double, double>())
-        .def("set_value", &SqNode::set_value)
-        .def("set_test_value", &SqNode::set_test_value)
-        .add_property("expr", &SqNode::expr)
-        .add_property("unit", &SqNode::unit)
-        .add_property("rung", &SqNode::rung)
-    ;
-}
\ No newline at end of file
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 850acfd5e6edf9d740c5a6e522506ca5f167d86d..8a84950126c1f625ccf319867af5713f41af685e 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
@@ -65,8 +65,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 60fa44e8d2ea8ad37c92d420231faa96f48aa843..ea54c2cd7578bce3439fa57109a3f239a32fe8fc 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
@@ -44,15 +44,3 @@ void SqrtNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_le
 {
     _feats[0]->update_div_mult_leaves(div_mult_leaves, fact / 2.0, expected_abs_tot);
 }
-
-void SqrtNode::register_python()
-{
-    using namespace boost::python;
-    class_<SqrtNode, bases<OperatorNode<1>>>("SqrtNode", init<node_ptr, int, double, double>())
-        .def("set_value", &SqrtNode::set_value)
-        .def("set_test_value", &SqrtNode::set_test_value)
-        .add_property("expr", &SqrtNode::expr)
-        .add_property("unit", &SqrtNode::unit)
-        .add_property("rung", &SqrtNode::rung)
-    ;
-}
\ No newline at end of file
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 e06b3bf01518b25426e7dc6a8f46a7ef2a2fdf60..3f327423669998ce8b3ead9e4532b33f4924690f 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
@@ -67,7 +67,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
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 6ea320c0e6b149c12cc418cf2c48ebb0dadcdb76..7b454133185999e836df8456375e4893c0dcbc79 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
@@ -75,15 +75,3 @@ void SubNode::update_div_mult_leaves(std::map<std::string, double>& div_mult_lea
 
     expected_abs_tot += std::abs(fact);
 }
-
-void SubNode::register_python()
-{
-    using namespace boost::python;
-    class_<SubNode, bases<OperatorNode<2>>>("SubNode", init<node_ptr, node_ptr, int, double, double>())
-        .def("set_value", &SubNode::set_value)
-        .def("set_test_value", &SubNode::set_test_value)
-        .add_property("expr", &SubNode::expr)
-        .add_property("unit", &SubNode::unit)
-        .add_property("rung", &SubNode::rung)
-    ;
-}
\ No newline at end of file
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 49ee521d310e88a73880ed8d4d3658b2be7a379f..40b87a26f95f034473ab936701d74be108315938 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
@@ -67,7 +67,6 @@ public:
      * @param expected_abs_tot The expected absolute sum of all values in div_mult_leaves
      */
     void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot);
-    static void register_python();
 };
 
 #endif
\ No newline at end of file
diff --git a/src/feature_creation/units/Unit.cpp b/src/feature_creation/units/Unit.cpp
index c95599f437c581c9eba381dfd21519ce0495ad00..ce4520ddda7dbe4e351803fd97b6ddb76c5db7a6 100644
--- a/src/feature_creation/units/Unit.cpp
+++ b/src/feature_creation/units/Unit.cpp
@@ -182,24 +182,3 @@ std::ostream& operator<< (std::ostream& outStream, const Unit& unit)
     outStream << unit.toString();
     return outStream;
 }
-
-
-void Unit::register_python()
-{
-    using namespace boost::python;
-    class_<Unit>("Unit", init<>())
-        .def(init<std::map<std::string, double>>())
-        .def(init<std::string>())
-        .def(init<Unit&>())
-        .def("__str__", &Unit::toString)
-        .def("__repr__", &Unit::toString)
-        .def("inverse", &Unit::inverse)
-        .def(self * self)
-        .def(self / self)
-        .def(self *= self)
-        .def(self /= self)
-        .def(self == self)
-        .def(self != self)
-        .def("__pow__", &Unit::operator^)
-    ;
-}
\ No newline at end of file
diff --git a/src/feature_creation/units/Unit.hpp b/src/feature_creation/units/Unit.hpp
index 3848d3a8eba64da4a161732c653f732e8239505c..8fe630594b8bcb5eefba124792403ae3130ccd29 100644
--- a/src/feature_creation/units/Unit.hpp
+++ b/src/feature_creation/units/Unit.hpp
@@ -12,7 +12,6 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/serialization/map.hpp>
 #include <boost/serialization/string.hpp>
-#include <boost/python.hpp>
 
 using StringRange = boost::iterator_range<std::string::const_iterator>;
 
@@ -142,8 +141,6 @@ public:
     {
         ar & _dct;
     }
-
-    static void register_python();
 };
 
 
diff --git a/src/python/_sisso.cpp b/src/python/_sisso.cpp
index 828ef3ee012d539c358a56634823ab7024107c22..c20670e72624bfc5ea234bab33034c8171a5e965 100644
--- a/src/python/_sisso.cpp
+++ b/src/python/_sisso.cpp
@@ -16,9 +16,7 @@ BOOST_PYTHON_MODULE(_sisso)
     mpi_setup::init_mpi_env();
     allowed_op_maps::set_node_maps();
 
-    sisso::register_python();
-    class_<std::map<std::string, double> >("unit_dict")
-        .def(boost::python::map_indexing_suite<std::map<std::wstring, double> >() );
+    sisso::register_all();
 
     Py_AtExit(&finalize);
 }
diff --git a/src/python/bindings.cpp b/src/python/bindings.cpp
index d79bdf65222fcdc21bbd30acf22e5f229d823930..102f389d0de2c1f4269f43e947aa93860ce60c17 100644
--- a/src/python/bindings.cpp
+++ b/src/python/bindings.cpp
@@ -1,31 +1,372 @@
 #include <python/bindings.hpp>
 
-void sisso::register_python()
-{
-    Model::register_python();
-    SISSORegressor::register_python();
-    FeatureSpace::register_python();
-    Unit::register_python();
-    Node::register_python();
-    FeatureNode::register_python();
-    ModelNode::register_python();
-    OperatorNode<1>::register_python();
-    OperatorNode<2>::register_python();
-    AddNode::register_python();
-    SubNode::register_python();
-    DivNode::register_python();
-    MultNode::register_python();
-    AbsDiffNode::register_python();
-    AbsNode::register_python();
-    InvNode::register_python();
-    LogNode::register_python();
-    ExpNode::register_python();
-    NegExpNode::register_python();
-    SinNode::register_python();
-    CosNode::register_python();
-    CbNode::register_python();
-    CbrtNode::register_python();
-    SqNode::register_python();
-    SqrtNode::register_python();
-    SixPowNode::register_python();
-}
\ No newline at end of file
+using namespace boost::python;
+
+void sisso::register_all()
+{
+    sisso::descriptor_identifier::registerModel();
+    sisso::descriptor_identifier::registerSISSORegressor();
+    sisso::feature_creation::registerFeatureSpace();
+    sisso::feature_creation::registerUnit();
+    sisso::feature_creation::node::registerNode();
+    sisso::feature_creation::node::registerFeatureNode();
+    sisso::feature_creation::node::registerModelNode();
+    sisso::feature_creation::node::registerOperatorNode<1>();
+    sisso::feature_creation::node::registerOperatorNode<2>();
+    sisso::feature_creation::node::registerAddNode();
+    sisso::feature_creation::node::registerSubNode();
+    sisso::feature_creation::node::registerDivNode();
+    sisso::feature_creation::node::registerMultNode();
+    sisso::feature_creation::node::registerAbsDiffNode();
+    sisso::feature_creation::node::registerAbsNode();
+    sisso::feature_creation::node::registerInvNode();
+    sisso::feature_creation::node::registerLogNode();
+    sisso::feature_creation::node::registerExpNode();
+    sisso::feature_creation::node::registerNegExpNode();
+    sisso::feature_creation::node::registerSinNode();
+    sisso::feature_creation::node::registerCosNode();
+    sisso::feature_creation::node::registerCbNode();
+    sisso::feature_creation::node::registerCbrtNode();
+    sisso::feature_creation::node::registerSqNode();
+    sisso::feature_creation::node::registerSqrtNode();
+    sisso::feature_creation::node::registerSixPowNode();
+}
+
+void sisso::feature_creation::registerFeatureSpace()
+{
+    void (FeatureSpace::*sis_list)(list) = &FeatureSpace::sis;
+    void (FeatureSpace::*sis_ndarray)(np::ndarray) = &FeatureSpace::sis;
+
+    class_<FeatureSpace>("FeatureSpace", init<list, list, np::ndarray, list, optional<int, int, int, int, double, double>>())
+        .def(init<list, list, list, list, optional<int, int, int, int, double, double>>())
+        .def("sis", sis_list)
+        .def("sis", sis_ndarray)
+        .def("feat_in_phi", &FeatureSpace::feat_in_phi)
+        .add_property("phi_selected", &FeatureSpace::phi_selected_py)
+        .add_property("phi0", &FeatureSpace::phi0_py)
+        .add_property("scores", &FeatureSpace::scores_py)
+        .add_property("task_sizes", &FeatureSpace::task_sizes_py)
+        .add_property("allowed_ops", &FeatureSpace::allowed_ops_py)
+        .add_property("start_gen", &FeatureSpace::start_gen_py)
+        .add_property("feature_space_file", &FeatureSpace::feature_space_file)
+        .add_property("l_bound", &FeatureSpace::l_bound)
+        .add_property("u_bound", &FeatureSpace::u_bound)
+        .add_property("max_phi", &FeatureSpace::max_phi)
+        .add_property("n_sis_select", &FeatureSpace::n_sis_select)
+        .add_property("n_samp", &FeatureSpace::n_samp)
+        .add_property("n_feat", &FeatureSpace::n_feat)
+        .add_property("n_rung_store", &FeatureSpace::n_rung_store)
+        .add_property("n_rung_generate", &FeatureSpace::n_rung_generate)
+    ;
+}
+
+void sisso::feature_creation::registerUnit()
+{
+    class_<Unit>("Unit", init<>())
+        .def(init<std::map<std::string, double>>())
+        .def(init<std::string>())
+        .def(init<Unit&>())
+        .def("__str__", &Unit::toString)
+        .def("__repr__", &Unit::toString)
+        .def("inverse", &Unit::inverse)
+        .def(self * self)
+        .def(self / self)
+        .def(self *= self)
+        .def(self /= self)
+        .def(self == self)
+        .def(self != self)
+        .def("__pow__", &Unit::operator^)
+    ;
+}
+
+void sisso::feature_creation::node::registerNode()
+{
+    void (Node::*reindex_1)(int) = &Node::reindex;
+    void (Node::*reindex_2)(int, int) = &Node::reindex;
+    class_<sisso::feature_creation::node::NodeWrap, boost::noncopyable>("Node", no_init)
+        .def("reindex", reindex_1)
+        .def("reindex", reindex_2)
+        .def("__str__", &Node::expr)
+        .def("__repr__", &Node::expr)
+        .add_property("n_samp", &Node::n_samp)
+        .add_property("n_test_samp", &Node::n_test_samp)
+        .add_property("feat_ind", &Node::feat_ind)
+        .add_property("arr_ind", &Node::arr_ind)
+        .add_property("selected", &Node::selected, &Node::set_selected)
+        .add_property("d_mat_ind", &Node::d_mat_ind, &Node::set_d_mat_ind)
+        .add_property("value", &Node::value_py)
+        .add_property("test_value", &Node::test_value_py)
+        .def("expr", pure_virtual(&Node::expr))
+        .def("unit", pure_virtual(&Node::unit))
+        .def("set_value", pure_virtual(&Node::set_value))
+        .def("set_test_value", pure_virtual(&Node::set_test_value))
+        .def("is_nan", pure_virtual(&Node::is_nan))
+        .def("is_const", pure_virtual(&Node::is_const))
+        .def("rung", pure_virtual(&Node::rung))
+    ;
+}
+
+void sisso::feature_creation::node::registerFeatureNode()
+{
+    std::string (FeatureNode::*expr_1)() = &FeatureNode::expr;
+    std::string (FeatureNode::*expr_const)() const = &FeatureNode::expr;
+
+    using namespace boost::python;
+    class_<FeatureNode, bases<Node>>("FeatureNode", init<int, std::string, np::ndarray, np::ndarray, Unit>())
+        .def(init<int, std::string, py::list, py::list, Unit>())
+        .def("is_nan", &FeatureNode::is_nan)
+        .def("is_const", &FeatureNode::is_const)
+        .def("set_value", &FeatureNode::set_value)
+        .def("set_test_value", &FeatureNode::set_test_value)
+        .add_property("expr", expr_1)
+        .add_property("expr", expr_const)
+        .add_property("unit", &FeatureNode::unit)
+        .add_property("rung", &FeatureNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerModelNode()
+{
+    std::string (ModelNode::*expr_1)() = &ModelNode::expr;
+    std::string (ModelNode::*expr_const)() const = &ModelNode::expr;
+
+    using namespace boost::python;
+    class_<ModelNode, bases<FeatureNode>>("ModelNode", init<int, int, std::string, std::vector<double>, std::vector<double>, Unit>())
+        .def("is_nan", &ModelNode::is_nan)
+        .def("is_const", &ModelNode::is_const)
+        .def("set_value", &ModelNode::set_value)
+        .def("set_test_value", &ModelNode::set_test_value)
+        .add_property("rung", &ModelNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerAddNode()
+{
+    class_<AddNode, bases<OperatorNode<2>>>("AddNode", init<node_ptr, node_ptr, int, double, double>())
+        .def("set_value", &AddNode::set_value)
+        .def("set_test_value", &AddNode::set_test_value)
+        .add_property("expr", &AddNode::expr)
+        .add_property("unit", &AddNode::unit)
+        .add_property("rung", &AddNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerSubNode()
+{
+    class_<SubNode, bases<OperatorNode<2>>>("SubNode", init<node_ptr, node_ptr, int, double, double>())
+        .def("set_value", &SubNode::set_value)
+        .def("set_test_value", &SubNode::set_test_value)
+        .add_property("expr", &SubNode::expr)
+        .add_property("unit", &SubNode::unit)
+        .add_property("rung", &SubNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerDivNode()
+{
+    class_<DivNode, bases<OperatorNode<2>>>("DivNode", init<node_ptr, node_ptr, int, double, double>())
+        .def("set_value", &DivNode::set_value)
+        .def("set_test_value", &DivNode::set_test_value)
+        .add_property("expr", &DivNode::expr)
+        .add_property("unit", &DivNode::unit)
+        .add_property("rung", &DivNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerMultNode()
+{
+    class_<MultNode, bases<OperatorNode<2>>>("MultNode", init<node_ptr, node_ptr, int, double, double>())
+        .def("set_value", &MultNode::set_value)
+        .def("set_test_value", &MultNode::set_test_value)
+        .add_property("expr", &MultNode::expr)
+        .add_property("unit", &MultNode::unit)
+        .add_property("rung", &MultNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerAbsDiffNode()
+{
+    class_<AbsDiffNode, bases<OperatorNode<2>>>("AbsDiffNode", init<node_ptr, node_ptr, int, double, double>())
+        .def("set_value", &AbsDiffNode::set_value)
+        .def("set_test_value", &AbsDiffNode::set_test_value)
+        .add_property("expr", &AbsDiffNode::expr)
+        .add_property("unit", &AbsDiffNode::unit)
+        .add_property("rung", &AbsDiffNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerAbsNode()
+{
+    class_<AbsNode, bases<OperatorNode<1>>>("AbsNode", init<node_ptr, int, double, double>())
+        .def("set_value", &AbsNode::set_value)
+        .def("set_test_value", &AbsNode::set_test_value)
+        .add_property("expr", &AbsNode::expr)
+        .add_property("unit", &AbsNode::unit)
+        .add_property("rung", &AbsNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerInvNode()
+{
+    class_<InvNode, bases<OperatorNode<1>>>("InvNode", init<node_ptr, int, double, double>())
+        .def("set_value", &InvNode::set_value)
+        .def("set_test_value", &InvNode::set_test_value)
+        .add_property("expr", &InvNode::expr)
+        .add_property("unit", &InvNode::unit)
+        .add_property("rung", &InvNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerLogNode()
+{
+    class_<LogNode, bases<OperatorNode<1>>>("LogNode", init<node_ptr, int, double, double>())
+        .def("set_value", &LogNode::set_value)
+        .def("set_test_value", &LogNode::set_test_value)
+        .add_property("expr", &LogNode::expr)
+        .add_property("unit", &LogNode::unit)
+        .add_property("rung", &LogNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerExpNode()
+{
+    class_<ExpNode, bases<OperatorNode<1>>>("ExpNode", init<node_ptr, int, double, double>())
+        .def("set_value", &ExpNode::set_value)
+        .def("set_test_value", &ExpNode::set_test_value)
+        .add_property("expr", &ExpNode::expr)
+        .add_property("unit", &ExpNode::unit)
+        .add_property("rung", &ExpNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerNegExpNode()
+{
+    class_<NegExpNode, bases<OperatorNode<1>>>("NegExpNode", init<node_ptr, int, double, double>())
+        .def("set_value", &NegExpNode::set_value)
+        .def("set_test_value", &NegExpNode::set_test_value)
+        .add_property("expr", &NegExpNode::expr)
+        .add_property("unit", &NegExpNode::unit)
+        .add_property("rung", &NegExpNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerSinNode()
+{
+    class_<SinNode, bases<OperatorNode<1>>>("SinNode", init<node_ptr, int, double, double>())
+        .def("set_value", &SinNode::set_value)
+        .def("set_test_value", &SinNode::set_test_value)
+        .add_property("expr", &SinNode::expr)
+        .add_property("unit", &SinNode::unit)
+        .add_property("rung", &SinNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerCosNode()
+{
+    class_<CosNode, bases<OperatorNode<1>>>("CosNode", init<node_ptr, int, double, double>())
+        .def("set_value", &CosNode::set_value)
+        .def("set_test_value", &CosNode::set_test_value)
+        .add_property("expr", &CosNode::expr)
+        .add_property("unit", &CosNode::unit)
+        .add_property("rung", &CosNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerCbNode()
+{
+    class_<CbNode, bases<OperatorNode<1>>>("CbNode", init<node_ptr, int, double, double>())
+        .def("set_value", &CbNode::set_value)
+        .def("set_test_value", &CbNode::set_test_value)
+        .add_property("expr", &CbNode::expr)
+        .add_property("unit", &CbNode::unit)
+        .add_property("rung", &CbNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerCbrtNode()
+{
+    class_<CbrtNode, bases<OperatorNode<1>>>("CbrtNode", init<node_ptr, int, double, double>())
+        .def("set_value", &CbrtNode::set_value)
+        .def("set_test_value", &CbrtNode::set_test_value)
+        .add_property("expr", &CbrtNode::expr)
+        .add_property("unit", &CbrtNode::unit)
+        .add_property("rung", &CbrtNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerSqNode()
+{
+    class_<SqNode, bases<OperatorNode<1>>>("SqNode", init<node_ptr, int, double, double>())
+        .def("set_value", &SqNode::set_value)
+        .def("set_test_value", &SqNode::set_test_value)
+        .add_property("expr", &SqNode::expr)
+        .add_property("unit", &SqNode::unit)
+        .add_property("rung", &SqNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerSqrtNode()
+{
+    class_<SqrtNode, bases<OperatorNode<1>>>("SqrtNode", init<node_ptr, int, double, double>())
+        .def("set_value", &SqrtNode::set_value)
+        .def("set_test_value", &SqrtNode::set_test_value)
+        .add_property("expr", &SqrtNode::expr)
+        .add_property("unit", &SqrtNode::unit)
+        .add_property("rung", &SqrtNode::rung)
+    ;
+}
+
+void sisso::feature_creation::node::registerSixPowNode()
+{
+    class_<SixPowNode, bases<OperatorNode<1>>>("SixPowNode", init<node_ptr, int, double, double>())
+        .def("set_value", &SixPowNode::set_value)
+        .def("set_test_value", &SixPowNode::set_test_value)
+        .add_property("expr", &SixPowNode::expr)
+        .add_property("unit", &SixPowNode::unit)
+        .add_property("rung", &SixPowNode::rung)
+    ;
+}
+
+void sisso::descriptor_identifier::registerModel()
+{
+    class_<Model>("Model", init<std::vector<double>, std::vector<double>, std::vector<model_node_ptr>, std::vector<int>, std::vector<int>>())
+        .def(init<std::string>())
+        .def(init<std::string, std::string>())
+        .def("predict", &Model::predict)
+        .def("fit", &Model::predict_train)
+        .def("__str__", &Model::toString)
+        .def("__repr__", &Model::toString)
+        .add_property("n_samp_train", &Model::n_samp_train)
+        .add_property("n_samp_test", &Model::n_samp_test)
+        .add_property("n_dim", &Model::n_dim)
+        .add_property("prop_train_est", &Model::prop_train_est)
+        .add_property("prop_test_est", &Model::prop_test_est)
+        .add_property("prop_train", &Model::prop_train)
+        .add_property("prop_test", &Model::prop_test)
+        .add_property("train_error", &Model::train_error)
+        .add_property("test_error", &Model::test_error)
+        .add_property("feats", &Model::feats)
+        .add_property("coefs", &Model::coefs)
+        .add_property("rmse", &Model::rmse)
+        .add_property("test_rmse", &Model::test_rmse)
+        .add_property("max_ae", &Model::max_ae)
+        .add_property("test_max_ae", &Model::test_max_ae);
+}
+
+void sisso::descriptor_identifier::registerSISSORegressor()
+{
+    class_<SISSORegressor>("SISSORegressor", init<std::shared_ptr<FeatureSpace>, np::ndarray, np::ndarray, py::list, py::list, py::list, int, int>())
+        .def(init<std::shared_ptr<FeatureSpace>, py::list, py::list, py::list, py::list, py::list, int, int>())
+        .def("fit", &SISSORegressor::fit)
+        .add_property("prop", &SISSORegressor::prop_py)
+        .add_property("prop_test", &SISSORegressor::prop_test_py)
+        .add_property("models", &SISSORegressor::models_py)
+        .add_property("n_samp", &SISSORegressor::n_samp)
+        .add_property("n_dim", &SISSORegressor::n_dim)
+        .add_property("n_residual", &SISSORegressor::n_residual)
+        .add_property("feat_space", &SISSORegressor::feat_space)
+        .add_property("error", &SISSORegressor::error_py)
+        .add_property("task_sizes_train", &SISSORegressor::task_sizes_train)
+        .add_property("task_sizes_test", &SISSORegressor::task_sizes_test)
+    ;
+}
diff --git a/src/python/bindings.hpp b/src/python/bindings.hpp
index ee80b8c6553c182cdab158cdcfc16d21339e9bd8..08940f5dcc16b6f1046aa03a8ac6675f14a045d5 100644
--- a/src/python/bindings.hpp
+++ b/src/python/bindings.hpp
@@ -1,13 +1,96 @@
-
-#ifndef PYTHON_BINDINGS
+#ifndef PYTHON_BINDING
 #define PYTHON_BINDINGS
 
 #include <descriptor_identifier/SISSORegressor.hpp>
 #include <feature_creation/feature_space/FeatureSpace.hpp>
 
+namespace py = boost::python;
+namespace np = boost::python::numpy;
+
 namespace sisso
 {
-    void register_python();
+    void register_all();
+
+    namespace feature_creation
+    {
+        static void registerFeatureSpace();
+        static void registerUnit();
+        namespace node
+        {
+            struct NodeWrap :  Node, py::wrapper<Node>
+            {
+            public:
+                inline std::string expr(){return this->get_override("expr")();}
+                inline Unit unit(){return this->get_override("unit")();}
+                inline std::vector<double> value(){return this->get_override("value")();}
+                inline std::vector<double> test_value(){return this->get_override("test_value")();}
+                inline void set_value(int offset = -1){this->get_override("set_value")();}
+                inline double* value_ptr(int offset = -1){return this->get_override("value_ptr")();}
+                inline void set_test_value(int offset = -1){this->get_override("set_test_value")();}
+                inline double* test_value_ptr(int offset = -1){return this->get_override("test_value_ptr")();}
+                inline bool is_nan(){return this->get_override("is_nan")();}
+                inline bool is_const(){return this->get_override("is_const")();}
+                inline NODE_TYPE type(){return this->get_override("type")();}
+                inline int rung(int cur_rung = 0){return this->get_override("rung")();}
+                inline void update_add_sub_leaves(std::map<std::string, int>& add_sub_leaves, int pl_mn, int& expected_abs_tot){this->get_override("update_add_sub_leaves");}
+                inline void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot){this->get_override("update_div_mult_leaves");}
+            };
+            template<int N>
+            struct OperatorNodeWrap : OperatorNode<N>, py::wrapper<OperatorNode<N>>
+            {
+                inline void set_value(int offset = -1){this->get_override("set_value")();}
+                inline void set_test_value(int offset = -1){this->get_override("set_test_value")();}
+                inline NODE_TYPE type(){return this->get_override("type")();}
+                inline int rung(int cur_rung = 0){return this->get_override("rung")();}
+                inline std::string expr(){return this->get_override("expr")();}
+                inline Unit unit(){return this->get_override("unit")();}
+                inline void update_add_sub_leaves(std::map<std::string, int>& add_sub_leaves, int pl_mn, int& expected_abs_tot){this->get_override("update_add_sub_leaves")();}
+                inline void update_div_mult_leaves(std::map<std::string, double>& div_mult_leaves, double fact, double& expected_abs_tot){this->get_override("update_div_mult_leaves")();}
+            };
+
+
+            static void registerNode();
+            static void registerFeatureNode();
+            static void registerModelNode();
+            template<int N>
+            static void registerOperatorNode()
+            {
+                py::class_<OperatorNodeWrap<N>, py::bases<Node>, boost::noncopyable>("OperatorNode")
+                    .def("is_nan", &OperatorNode<N>::is_nan)
+                    .def("is_const", &OperatorNode<N>::is_const)
+                    .def("set_value", py::pure_virtual(&OperatorNode<N>::set_value))
+                    .def("set_test_value", py::pure_virtual(&OperatorNode<N>::set_test_value))
+                    .def("rung", py::pure_virtual(&OperatorNode<N>::rung))
+                    .def("expr", py::pure_virtual(&OperatorNode<N>::expr))
+                    .def("unit", py::pure_virtual(&OperatorNode<N>::unit))
+                ;
+            }
+            static void registerAddNode();
+            static void registerSubNode();
+            static void registerDivNode();
+            static void registerMultNode();
+            static void registerAbsDiffNode();
+            static void registerAbsNode();
+            static void registerInvNode();
+            static void registerLogNode();
+            static void registerExpNode();
+            static void registerNegExpNode();
+            static void registerSinNode();
+            static void registerCosNode();
+            static void registerCbNode();
+            static void registerCbrtNode();
+            static void registerSqNode();
+            static void registerSqrtNode();
+            static void registerSixPowNode();
+
+        }
+    }
+
+    namespace descriptor_identifier
+    {
+        static void registerModel();
+        static void registerSISSORegressor();
+    }
 }
 
 #endif
\ No newline at end of file
diff --git a/src/python/descriptor_identifier/SISSORegressor.cpp b/src/python/descriptor_identifier/SISSORegressor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..99d5dbc23ef78bc3b518253025c024ae029fcc85
--- /dev/null
+++ b/src/python/descriptor_identifier/SISSORegressor.cpp
@@ -0,0 +1,90 @@
+#include <descriptor_identifier/SISSORegressor.hpp>
+
+SISSORegressor::SISSORegressor(
+    std::shared_ptr<FeatureSpace> feat_space,
+    np::ndarray prop,
+    np::ndarray prop_test,
+    py::list task_sizes_train,
+    py::list task_sizes_test,
+    py::list leave_out_inds,
+    int n_dim,
+    int n_residual
+) :
+    _prop(python_conv_utils::from_ndarray<double>(prop)),
+    _prop_test(python_conv_utils::from_ndarray<double>(prop_test)),
+    _a((n_dim + 1) * prop.shape(0)),
+    _b(prop.shape(0)),
+    _error(prop.shape(0), 0.0),
+    _s(n_dim + 1),
+    _task_sizes_train(python_conv_utils::from_list<int>(task_sizes_train)),
+    _task_sizes_test(python_conv_utils::from_list<int>(task_sizes_test)),
+    _leave_out_inds(python_conv_utils::from_list<int>(leave_out_inds)),
+    _feat_space(feat_space),
+    _mpi_comm(feat_space->mpi_comm()),
+    _n_samp(prop.shape(0)),
+    _n_dim(n_dim),
+    _n_residual(n_residual),
+    _lwork(-1),
+    _rank(0)
+{
+    node_value_arrs::initialize_d_matrix_arr();
+    // Initialize a, b, ones, s, and _error arrays
+    std::fill_n(_a.data(), (_n_dim + 1) * _n_samp, 0.0);
+    std::fill_n(_b.data(), _n_samp, 0.0);
+    std::fill_n(_s.data(), _n_dim + 1, 0.0);
+
+    // // Get optimal work size
+    _lwork = get_opt_lwork(_n_dim + 1);
+    _work = std::vector<double>(_lwork, 0.0);
+}
+
+SISSORegressor::SISSORegressor(
+    std::shared_ptr<FeatureSpace> feat_space,
+    py::list prop,
+    py::list prop_test,
+    py::list task_sizes_train,
+    py::list task_sizes_test,
+    py::list leave_out_inds,
+    int n_dim,
+    int n_residual
+) :
+    _prop(python_conv_utils::from_list<double>(prop)),
+    _prop_test(python_conv_utils::from_list<double>(prop_test)),
+    _a((n_dim + 1) * py::len(prop)),
+    _b(py::len(prop)),
+    _error(py::len(prop), 0.0),
+    _s(n_dim + 1),
+    _task_sizes_train(python_conv_utils::from_list<int>(task_sizes_train)),
+    _task_sizes_test(python_conv_utils::from_list<int>(task_sizes_test)),
+    _leave_out_inds(python_conv_utils::from_list<int>(leave_out_inds)),
+    _feat_space(feat_space),
+    _mpi_comm(feat_space->mpi_comm()),
+    _n_samp(py::len(prop)),
+    _n_dim(n_dim),
+    _n_residual(n_residual),
+    _lwork(-1),
+    _rank(0)
+{
+    node_value_arrs::initialize_d_matrix_arr();
+    // Initialize a, b, ones, s, and _error arrays
+    std::fill_n(_a.data(), (_n_dim + 1) * _n_samp, 0.0);
+    std::fill_n(_b.data(), _n_samp, 0.0);
+    std::fill_n(_s.data(), _n_dim + 1, 0.0);
+
+    // // Get optimal work size
+    _lwork = get_opt_lwork(_n_dim + 1);
+    _work = std::vector<double>(_lwork, 0.0);
+}
+
+py::list SISSORegressor::models_py()
+{
+    py::list model_list;
+    for(auto& mod_list : _models)
+    {
+        py::list lst;
+        for(auto& mod : mod_list)
+            lst.append<Model>(mod);
+        model_list.append<py::list>(lst);
+    }
+    return model_list;
+}
diff --git a/src/python/feature_creation/FeatureNode.cpp b/src/python/feature_creation/FeatureNode.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b9330df9f09c4ffd29ae1b020bc97cb7df5cce67
--- /dev/null
+++ b/src/python/feature_creation/FeatureNode.cpp
@@ -0,0 +1,40 @@
+#include <feature_creation/node/FeatureNode.hpp>
+
+FeatureNode::FeatureNode(int feat_ind, std::string expr, np::ndarray value, np::ndarray test_value, Unit unit) :
+    Node(feat_ind, value.shape(0), test_value.shape(0)),
+    _value(python_conv_utils::from_ndarray<double>(value)),
+    _test_value(python_conv_utils::from_ndarray<double>(test_value)),
+    _unit(unit),
+    _expr(expr)
+{
+    // Automatically resize the storage arrays
+    if(node_value_arrs::N_STORE_FEATURES == 0)
+        node_value_arrs::initialize_values_arr(_n_samp, _n_test_samp, 1);
+    else if((_n_samp != node_value_arrs::N_SAMPLES) || (_n_test_samp != node_value_arrs::N_SAMPLES_TEST))
+        throw std::logic_error("Number of samples in current feature is not the same as the others, (" + std::to_string(_n_samp) + " and " + std::to_string(_n_test_samp) + " vs. "  + std::to_string(node_value_arrs::N_SAMPLES) + " and " + std::to_string(node_value_arrs::N_SAMPLES_TEST) + ")");
+    else if(feat_ind >= node_value_arrs::N_STORE_FEATURES)
+        node_value_arrs::resize_values_arr(0, node_value_arrs::N_STORE_FEATURES + 1, true);
+
+    set_value();
+    set_test_value();
+}
+
+FeatureNode::FeatureNode(int feat_ind, std::string expr, py::list value, py::list test_value, Unit unit) :
+    Node(feat_ind, py::len(value), py::len(test_value)),
+    _value(python_conv_utils::from_list<double>(value)),
+    _test_value(python_conv_utils::from_list<double>(test_value)),
+    _unit(unit),
+    _expr(expr)
+{
+
+    // Automatically resize the storage arrays
+    if(node_value_arrs::N_STORE_FEATURES == 0)
+        node_value_arrs::initialize_values_arr(_n_samp, _n_test_samp, 1);
+    else if((_n_samp != node_value_arrs::N_SAMPLES) || (_n_test_samp != node_value_arrs::N_SAMPLES_TEST))
+        throw std::logic_error("Number of samples in current feature is not the same as the others, (" + std::to_string(_n_samp) + " and " + std::to_string(_n_test_samp) + " vs. "  + std::to_string(node_value_arrs::N_SAMPLES) + " and " + std::to_string(node_value_arrs::N_SAMPLES_TEST) + ")");
+    else if(feat_ind >= node_value_arrs::N_STORE_FEATURES)
+        node_value_arrs::resize_values_arr(0, node_value_arrs::N_STORE_FEATURES + 1, true);
+
+    set_value();
+    set_test_value();
+}
diff --git a/src/python/feature_creation/FeatureSpace.cpp b/src/python/feature_creation/FeatureSpace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..589bb14ab2360a60dd932228408abd2df007d459
--- /dev/null
+++ b/src/python/feature_creation/FeatureSpace.cpp
@@ -0,0 +1,81 @@
+#include <feature_creation/feature_space/FeatureSpace.hpp>
+
+FeatureSpace::FeatureSpace(
+    py::list phi_0,
+    py::list allowed_ops,
+    py::list prop,
+    py::list task_sizes,
+    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
+):
+    _phi(python_conv_utils::shared_ptr_vec_from_list<Node, FeatureNode>(phi_0)),
+    _phi_0(_phi),
+    _allowed_ops(python_conv_utils::from_list<std::string>(allowed_ops)),
+    _scores(py::len(phi_0), 0.0),
+    _task_sizes(python_conv_utils::from_list<int>(task_sizes)),
+    _start_gen(1, 0),
+    _feature_space_file("feature_space/selected_features.txt"),
+    _mpi_comm(mpi_setup::comm),
+    _l_bound(min_abs_feat_val),
+    _u_bound(max_abs_feat_val),
+    _max_phi(max_phi),
+    _n_sis_select(n_sis_select),
+    _n_feat(py::len(phi_0)),
+    _n_rung_store(max_store_rung),
+    _n_rung_generate(n_rung_generate)
+{
+    _n_samp = _phi_0[0]->n_samp();
+    initialize_fs(python_conv_utils::from_list<double>(prop));
+}
+
+FeatureSpace::FeatureSpace(
+    py::list phi_0,
+    py::list allowed_ops,
+    np::ndarray prop,
+    py::list task_sizes,
+    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
+):
+    _phi(python_conv_utils::shared_ptr_vec_from_list<Node, FeatureNode>(phi_0)),
+    _phi_0(_phi),
+    _allowed_ops(python_conv_utils::from_list<std::string>(allowed_ops)),
+    _scores(py::len(phi_0), 0.0),
+    _task_sizes(python_conv_utils::from_list<int>(task_sizes)),
+    _start_gen(1, 0),
+    _feature_space_file("feature_space/selected_features.txt"),
+    _mpi_comm(mpi_setup::comm),
+    _l_bound(min_abs_feat_val),
+    _u_bound(max_abs_feat_val),
+    _max_phi(max_phi),
+    _n_sis_select(n_sis_select),
+    _n_feat(py::len(phi_0)),
+    _n_rung_store(max_store_rung),
+    _n_rung_generate(n_rung_generate)
+{
+    _n_samp = _phi_0[0]->n_samp();
+    initialize_fs(python_conv_utils::from_ndarray<double>(prop));
+}
+
+py::list FeatureSpace::phi0_py()
+{
+    py::list feat_lst;
+    for(auto& feat : _phi_0)
+        feat_lst.append<FeatureNode>(FeatureNode(feat->feat_ind(), feat->expr(), feat->value(), feat->test_value(), feat->unit()));
+    return feat_lst;
+}
+
+py::list FeatureSpace::phi_selected_py()
+{
+    py::list feat_lst;
+    for(auto& feat : _phi_selected)
+        feat_lst.append<ModelNode>(ModelNode(feat->d_mat_ind(), feat->rung(), feat->expr(), feat->value(), feat->test_value(), feat->unit()));
+    return feat_lst;
+}
diff --git a/src/utils/mkl_interface.hpp b/src/utils/mkl_interface.hpp
index edb228dec62520ba943b846334b148ad4186dc15..6382c7eff7c7662d28d613faa65c531b98f43ae4 100644
--- a/src/utils/mkl_interface.hpp
+++ b/src/utils/mkl_interface.hpp
@@ -2,7 +2,6 @@
  *  @brief wrapper functions to a mkl blas functionality
  *
  *  @author Thomas A. Purcell (tpurcell90)
- *  @author Joshua E. Szekely (jeszekely)
  *  @bug No known bugs
  */