From f3bfec003d0031d84c94923f1b8211b4b192259f Mon Sep 17 00:00:00 2001
From: Thomas Purcell <purcell@fhi-berlin.mpg.de>
Date: Sun, 5 Jul 2020 14:04:46 +0200
Subject: [PATCH] interface is working, but with some issues

---
 src/__init__.py                               |   0
 src/descriptor_identifier/SISSORegressor.cpp  |  82 ++++++++++--
 src/descriptor_identifier/SISSORegressor.hpp  |  17 +++
 .../feature_space/FeatureSpace.cpp            | 120 +++++++++++++++---
 .../feature_space/FeatureSpace.hpp            |  84 +++++++++++-
 src/feature_creation/node/FeatureNode.cpp     |  45 ++++++-
 src/feature_creation/node/FeatureNode.hpp     |   8 ++
 src/feature_creation/node/ModelNode.cpp       |   1 -
 src/feature_creation/node/Node.cpp            |   5 +-
 src/feature_creation/node/Node.hpp            |   8 ++
 .../node/operator_nodes/OperatorNode.hpp      |   3 -
 .../absolute_difference.cpp                   |   2 -
 .../allowed_operator_nodes/absolute_value.cpp |   2 -
 .../allowed_operator_nodes/add.cpp            |   2 -
 .../allowed_operator_nodes/cos.cpp            |   2 -
 .../allowed_operator_nodes/cube.cpp           |   2 -
 .../allowed_operator_nodes/cube_root.cpp      |   2 -
 .../allowed_operator_nodes/divide.cpp         |   2 -
 .../allowed_operator_nodes/exponential.cpp    |   2 -
 .../allowed_operator_nodes/inverse.cpp        |   2 -
 .../allowed_operator_nodes/log.cpp            |   2 -
 .../allowed_operator_nodes/multiply.cpp       |   2 -
 .../negative_exponential.cpp                  |   2 -
 .../allowed_operator_nodes/sin.cpp            |   2 -
 .../allowed_operator_nodes/sixth_power.cpp    |   2 -
 .../allowed_operator_nodes/square.cpp         |   2 -
 .../allowed_operator_nodes/square_root.cpp    |   2 -
 .../allowed_operator_nodes/subtract.cpp       |   2 -
 .../value_storage/nodes_value_containers.cpp  |  12 +-
 src/feature_creation/units/Unit.cpp           |   3 +-
 src/main.cpp                                  |  19 +--
 src/mpi_interface/MPI_Interface.cpp           |   2 +
 src/python/__init__.py                        |   0
 src/python/_sisso.cpp                         |  12 ++
 src/python/conversion_utils.hpp               |  64 ++++++++++
 src/sisso++.cmakein                           |   0
 src/src/python/__init__.py                    |   0
 src/utils/project.cpp                         |   1 +
 38 files changed, 428 insertions(+), 92 deletions(-)
 create mode 100644 src/__init__.py
 create mode 100644 src/python/__init__.py
 create mode 100644 src/python/conversion_utils.hpp
 create mode 100644 src/sisso++.cmakein
 create mode 100644 src/src/python/__init__.py

diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/descriptor_identifier/SISSORegressor.cpp b/src/descriptor_identifier/SISSORegressor.cpp
index 03d4c3b6..55fa4a03 100644
--- a/src/descriptor_identifier/SISSORegressor.cpp
+++ b/src/descriptor_identifier/SISSORegressor.cpp
@@ -28,6 +28,62 @@ SISSORegressor::SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::ve
     _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, 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)),
+    _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)
+{
+
+    // 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, 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)),
+    _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)
+{
+
+    // 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)
@@ -194,20 +250,30 @@ 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>, std::vector<double>, std::vector<double>, std::vector<int>, std::vector<int>, int, int>())
+    class_<SISSORegressor>("SISSORegressor", init<std::shared_ptr<FeatureSpace>, np::ndarray, np::ndarray, python::list, python::list, int, int>())
+        .def(init<std::shared_ptr<FeatureSpace>, python::list, python::list, python::list, python::list, int, int>())
         .def("fit", &SISSORegressor::fit)
-        .add_property("prop", &SISSORegressor::prop)
-        .add_property("prop_test", &SISSORegressor::prop_test)
-        .add_property("models", &SISSORegressor::models)
+        .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)
-        .def_readonly("_task_sizes_train", &SISSORegressor::_task_sizes_train)
-        .def_readonly("_task_sizes_test", &SISSORegressor::_task_sizes_test)
-        .def_readonly("_mpi_comm", &SISSORegressor::_mpi_comm);
+        .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 f47bdf1f..da4b400f 100644
--- a/src/descriptor_identifier/SISSORegressor.hpp
+++ b/src/descriptor_identifier/SISSORegressor.hpp
@@ -5,6 +5,8 @@
 #include <descriptor_identifier/Model/Model.hpp>
 #include <ctime>
 
+namespace python = boost::python;
+namespace np = boost::python::numpy;
 /**
  * @brief SISSO Regressor class, to find the best models, and store them
  *
@@ -45,6 +47,10 @@ public:
      */
     SISSORegressor(std::shared_ptr<FeatureSpace> feat_space, std::vector<double> prop, std::vector<double> prop_test, std::vector<int> task_sizes_train, std::vector<int> task_sizes_test, 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, 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, int n_dim, int n_residual);
+
     /**
      * @brief Get the optimal size of the working array
      *
@@ -99,16 +105,22 @@ 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 models
      */
     inline std::vector<std::vector<Model>> models(){return _models;}
 
+    python::list models_py();
+
     /**
      * @brief Acessor function for n_samp
      */
@@ -124,11 +136,16 @@ 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();
 };
 
diff --git a/src/feature_creation/feature_space/FeatureSpace.cpp b/src/feature_creation/feature_space/FeatureSpace.cpp
index 5cc4487e..7e9d4168 100644
--- a/src/feature_creation/feature_space/FeatureSpace.cpp
+++ b/src/feature_creation/feature_space/FeatureSpace.cpp
@@ -51,6 +51,86 @@ FeatureSpace::FeatureSpace(
     _max_temp_store(max_temp_store)
 
 {
+    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,
+    int max_temp_store,
+    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>(allowed_ops)),
+    _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),
+    _max_temp_store(max_temp_store)
+{
+    _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,
+    int max_temp_store,
+    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),
+    _max_temp_store(max_temp_store)
+{
+    _n_samp = _phi_0[0]->n_samp();
+    initialize_fs(python_conv_utils::from_ndarray<double>(prop));
+}
+
+void FeatureSpace::initialize_fs(std::vector<double> prop)
+{
+    if(_n_rung_store == -1)
+        _n_rung_store = _max_phi - 1;
+    if(_n_rung_generate > 1)
+        throw std::logic_error("A maximum of one rung can be generated on the fly.");
+    else if(_max_phi - _n_rung_generate < _n_rung_store)
+        throw std::logic_error("Requesting to store more rungs than what can be pre-generated.");
+
     if(_mpi_comm->rank() == 0)
     {
         std::ofstream out_file_stream = std::ofstream();
@@ -58,17 +138,13 @@ FeatureSpace::FeatureSpace(
         out_file_stream << std::setw(14) <<std::left << "# FEAT_ID" << std::setw(24) << std::left << "Score" << "Feature Expression" << std::endl;
         out_file_stream.close();
     }
+
     if(_max_temp_store != -1)
         _max_temp_store /= 3;
 
     _project = project_funcs::project_r;
 
-    if(_n_rung_generate > 1)
-        throw std::logic_error("A maximum of one rung can be generated on the fly.");
-    else if(_max_phi - _n_rung_generate < _n_rung_store)
-        throw std::logic_error("Requesting to store more rungs than what can be pre-generated.");
-
-    for(auto & op : allowed_ops)
+    for(auto & op : _allowed_ops)
     {
         if((op.compare("add") == 0) || (op.compare("sub") == 0) || (op.compare("mult") == 0) || (op.compare("abs_diff") == 0))
             _com_bin_operators.push_back(allowed_op_maps::binary_operator_map[op]);
@@ -77,11 +153,12 @@ FeatureSpace::FeatureSpace(
         else
             _un_operators.push_back(allowed_op_maps::unary_operator_map[op]);
     }
+
     generate_feature_space(prop);
     _scores.reserve(_phi.size());
     _scores.resize(_phi.size());
-}
 
+}
 
 void FeatureSpace::generate_new_feats(std::vector<node_ptr>::iterator& feat, std::vector<node_ptr>& feat_set, int& feat_ind, double l_bound, double u_bound)
 {
@@ -735,24 +812,27 @@ void FeatureSpace::sis(std::vector<double>& prop)
     }
     if(_mpi_comm->rank() == 0)
         out_file_stream.close();
+
 }
 
 void FeatureSpace::register_python()
 {
     using namespace boost::python;
-    class_<FeatureSpace>("FeatureSpace", init<std::shared_ptr<MPI_Interface>, std::vector<node_ptr>, std::vector<std::string>, std::vector<double>, std::vector<int>, optional<int, int, int, int, int, double, double>>())
-        .def("sis", &FeatureSpace::sis)
+    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, int, double, double>>())
+        .def(init<list, list, list, list, optional<int, 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)
-        .add_property("phi", &FeatureSpace::phi)
-        .add_property("phi0", &FeatureSpace::phi0)
-        .add_property("scores", &FeatureSpace::scores)
-        .add_property("task_sizes", &FeatureSpace::task_sizes)
-        .def_readonly("_allowed_ops", &FeatureSpace::_allowed_ops)
-        .def_readonly("_un_operators", &FeatureSpace::_un_operators)
-        .def_readonly("_com_bin_operators", &FeatureSpace::_com_bin_operators)
-        .def_readonly("_bin_operators", &FeatureSpace::_bin_operators)
-        .def_readonly("_start_gen", &FeatureSpace::_start_gen)
+        .add_property("phi_selected", &FeatureSpace::phi_selected_py)
+        .add_property("phi", &FeatureSpace::phi_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)
@@ -762,7 +842,7 @@ void FeatureSpace::register_python()
         .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)
-        .def_readonly("_mpi_comm", &FeatureSpace::_mpi_comm);
+    ;
 }
 
 
diff --git a/src/feature_creation/feature_space/FeatureSpace.hpp b/src/feature_creation/feature_space/FeatureSpace.hpp
index ffd9e2d6..8caf3cef 100644
--- a/src/feature_creation/feature_space/FeatureSpace.hpp
+++ b/src/feature_creation/feature_space/FeatureSpace.hpp
@@ -14,7 +14,9 @@
 #include <iostream>
 #include <iomanip>
 
-// namespace mpi = boost::mpi;
+namespace python = boost::python;
+namespace np = boost::python::numpy;
+
 /**
  * @brief Feature Space for SISSO calculations
  * @details Stores and performs all feature calculations for SIS
@@ -70,13 +72,63 @@ public:
         std::vector<int> task_sizes,
         int max_phi=1,
         int n_sis_select=1,
-        int max_store_rung=2,
+        int max_store_rung=-1,
+        int n_rung_generate=0,
+        int max_temp_store=-1,
+        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,
+        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,
+        int max_temp_store=-1,
+        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,
         int max_temp_store=-1,
         double min_abs_feat_val=1e-50,
         double max_abs_feat_val=1e50
     );
 
+    void initialize_fs(std::vector<double> prop);
+
     /**
      * @brief Generate the full feature set from the allowed operators and initial feature set
      * @details populates phi with all features from an initial set and the allowed operators
@@ -88,20 +140,28 @@ public:
      */
     inline std::vector<node_ptr> phi_selected(){return _phi_selected;};
 
+    inline boost::python::list phi_selected_py(){return python_conv_utils::to_list<node_ptr>(_phi_selected);};
+
     /**
      * @brief Accessor function for _phi
      */
     inline std::vector<node_ptr> phi(){return _phi;};
 
+    inline boost::python::list phi_py(){return python_conv_utils::to_list<node_ptr>(_phi);};
+
     /**
      * @brief Accessor function for _phi_0
      */
     inline std::vector<node_ptr> phi0(){return _phi_0;};
 
+    inline boost::python::list phi0_py(){return python_conv_utils::to_list<node_ptr>(_phi_0);};
+
     /**
      * @brief Accessor function for _scores
      */
-    inline std::vector<double> scores(){return _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
@@ -110,6 +170,12 @@ 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);}
+
     void generate_new_feats(std::vector<node_ptr>::iterator& feat, std::vector<node_ptr>& feat_set, int& feat_ind, double l_bound=1e-50, double u_bound=1e50);
 
     void project_generated(double* prop, int size, std::vector<node_ptr>& phi_selected, std::vector<double>& scores_selected, std::vector<double>& scores_comp);
@@ -125,6 +191,18 @@ 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?
      *
diff --git a/src/feature_creation/node/FeatureNode.cpp b/src/feature_creation/node/FeatureNode.cpp
index 457ee7f5..abdadcee 100644
--- a/src/feature_creation/node/FeatureNode.cpp
+++ b/src/feature_creation/node/FeatureNode.cpp
@@ -14,6 +14,45 @@ FeatureNode::FeatureNode(int feat_ind, std::string expr, std::vector<double> val
     set_test_value();
 }
 
+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, 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, node_value_arrs::N_STORE_FEATURES + 1, true);
+
+    set_value();
+    set_test_value();
+}
+
 FeatureNode::~FeatureNode()
 {}
 
@@ -43,17 +82,15 @@ void FeatureNode::register_python()
     std::string (FeatureNode::*expr_const)() const = &FeatureNode::expr;
 
     using namespace boost::python;
-    class_<FeatureNode, bases<Node>>("FeatureNode", init<int, std::string, std::vector<double>, std::vector<double>, Unit>())
+    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("value", &FeatureNode::value)
-        .add_property("test_value", &FeatureNode::test_value)
         .add_property("expr", expr_1)
         .add_property("expr", expr_const)
         .add_property("unit", &FeatureNode::unit)
-        .add_property("type", &FeatureNode::type)
         .add_property("rung", &FeatureNode::rung)
     ;
 }
diff --git a/src/feature_creation/node/FeatureNode.hpp b/src/feature_creation/node/FeatureNode.hpp
index 0eaaa89d..1ef04770 100644
--- a/src/feature_creation/node/FeatureNode.hpp
+++ b/src/feature_creation/node/FeatureNode.hpp
@@ -1,15 +1,21 @@
 #ifndef FEATURE_NODE
 #define FEATURE_NODE
 
+#include <python/conversion_utils.hpp>
 #include <utils/math_funcs.hpp>
 #include <utils/enum.hpp>
+
 #include <memory>
 
 #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;
+
 /**
  * @brief Node that describe the leaves of the operator graph (Initial features in Phi_0)
  */
@@ -55,6 +61,8 @@ 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);
+    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);
 
     FeatureNode(const FeatureNode&) = default;
     FeatureNode(FeatureNode&&) = default;
diff --git a/src/feature_creation/node/ModelNode.cpp b/src/feature_creation/node/ModelNode.cpp
index cc48ac80..75d25160 100644
--- a/src/feature_creation/node/ModelNode.cpp
+++ b/src/feature_creation/node/ModelNode.cpp
@@ -42,7 +42,6 @@ void ModelNode::register_python()
         .def("is_const", &ModelNode::is_const)
         .def("set_value", &ModelNode::set_value)
         .def("set_test_value", &ModelNode::set_test_value)
-        .add_property("type", &ModelNode::type)
         .add_property("rung", &ModelNode::rung)
     ;
 }
diff --git a/src/feature_creation/node/Node.cpp b/src/feature_creation/node/Node.cpp
index c11fabee..3bde8ded 100644
--- a/src/feature_creation/node/Node.cpp
+++ b/src/feature_creation/node/Node.cpp
@@ -103,15 +103,14 @@ void Node::register_python()
         .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("value", pure_virtual(&Node::value))
-        .def("test_value", pure_virtual(&Node::test_value))
         .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("type", pure_virtual(&Node::type))
         .def("rung", pure_virtual(&Node::rung))
         ;
 }
diff --git a/src/feature_creation/node/Node.hpp b/src/feature_creation/node/Node.hpp
index 131dd159..3854df47 100644
--- a/src/feature_creation/node/Node.hpp
+++ b/src/feature_creation/node/Node.hpp
@@ -3,17 +3,21 @@
 
 #include <feature_creation/units/Unit.hpp>
 #include <feature_creation/node/value_storage/nodes_value_containers.hpp>
+#include <python/conversion_utils.hpp>
 #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;
+
 /**
  * @brief Base class for a Node
  * @details Class used to describe a Node on the descriptor graph. Features are treated as an operation graph, these are the nodes on that graph.
@@ -130,11 +134,15 @@ 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
      */
diff --git a/src/feature_creation/node/operator_nodes/OperatorNode.hpp b/src/feature_creation/node/operator_nodes/OperatorNode.hpp
index e8f4a578..c73f4ab2 100644
--- a/src/feature_creation/node/operator_nodes/OperatorNode.hpp
+++ b/src/feature_creation/node/operator_nodes/OperatorNode.hpp
@@ -217,11 +217,8 @@ public:
         class_<OperatorNodeWrap, bases<Node>, boost::noncopyable>("OperatorNode")
             .def("is_nan", &OperatorNode::is_nan)
             .def("is_const", &OperatorNode::is_const)
-            .add_property("value", &OperatorNode::value)
-            .add_property("test_value", &OperatorNode::test_value)
             .def("set_value", pure_virtual(&OperatorNode::set_value))
             .def("set_test_value", pure_virtual(&OperatorNode::set_test_value))
-            .def("type", pure_virtual(&OperatorNode::type))
             .def("rung", pure_virtual(&OperatorNode::rung))
             .def("expr", pure_virtual(&OperatorNode::expr))
             .def("unit", pure_virtual(&OperatorNode::unit))
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 59a62cfd..78e70780 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
@@ -89,12 +89,10 @@ void AbsDiffNode::register_python()
 {
     using namespace boost::python;
     class_<AbsDiffNode, bases<OperatorNode<2>>>("AbsDiffNode", init<node_ptr, node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 2>, 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("type", &AbsDiffNode::type)
         .add_property("rung", &AbsDiffNode::rung)
     ;
 }
\ 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 fe476da4..145e5ac8 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
@@ -51,12 +51,10 @@ void AbsNode::register_python()
 {
     using namespace boost::python;
     class_<AbsNode, bases<OperatorNode<1>>>("AbsNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &AbsNode::type)
         .add_property("rung", &AbsNode::rung)
     ;
 }
\ 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 d6216084..031a1fbc 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
@@ -80,12 +80,10 @@ void AddNode::register_python()
 {
     using namespace boost::python;
     class_<AddNode, bases<OperatorNode<2>>>("AddNode", init<node_ptr, node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 2>, 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("type", &AddNode::type)
         .add_property("rung", &AddNode::rung)
     ;
 }
\ 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 2e14154a..d611f00d 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
@@ -61,12 +61,10 @@ void CosNode::register_python()
 {
     using namespace boost::python;
     class_<CosNode, bases<OperatorNode<1>>>("CosNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &CosNode::type)
         .add_property("rung", &CosNode::rung)
     ;
 }
\ 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 055b3361..5a57120b 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
@@ -49,12 +49,10 @@ void CbNode::register_python()
 {
     using namespace boost::python;
     class_<CbNode, bases<OperatorNode<1>>>("CbNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &CbNode::type)
         .add_property("rung", &CbNode::rung)
     ;
 }
\ 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 6cf69aee..b31ef4de 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
@@ -49,12 +49,10 @@ void CbrtNode::register_python()
 {
     using namespace boost::python;
     class_<CbrtNode, bases<OperatorNode<1>>>("CbrtNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &CbrtNode::type)
         .add_property("rung", &CbrtNode::rung)
     ;
 }
\ 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 d07c20c2..5b371c2c 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
@@ -80,12 +80,10 @@ void DivNode::register_python()
 {
     using namespace boost::python;
     class_<DivNode, bases<OperatorNode<2>>>("DivNode", init<node_ptr, node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 2>, 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("type", &DivNode::type)
         .add_property("rung", &DivNode::rung)
     ;
 }
\ 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 5d97e99b..08a12828 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
@@ -61,12 +61,10 @@ void ExpNode::register_python()
 {
     using namespace boost::python;
     class_<ExpNode, bases<OperatorNode<1>>>("ExpNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &ExpNode::type)
         .add_property("rung", &ExpNode::rung)
     ;
 }
\ 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 5e35fdb6..371e10c3 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
@@ -49,12 +49,10 @@ void InvNode::register_python()
 {
     using namespace boost::python;
     class_<InvNode, bases<OperatorNode<1>>>("InvNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &InvNode::type)
         .add_property("rung", &InvNode::rung)
     ;
 }
\ 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 7b6eab0c..642f8014 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
@@ -61,12 +61,10 @@ void LogNode::register_python()
 {
     using namespace boost::python;
     class_<LogNode, bases<OperatorNode<1>>>("LogNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &LogNode::type)
         .add_property("rung", &LogNode::rung)
     ;
 }
\ 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 66a8bbd1..36efccdf 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
@@ -80,12 +80,10 @@ void MultNode::register_python()
 {
     using namespace boost::python;
     class_<MultNode, bases<OperatorNode<2>>>("MultNode", init<node_ptr, node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 2>, 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("type", &MultNode::type)
         .add_property("rung", &MultNode::rung)
     ;
 }
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 c92c426a..e6eccf80 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
@@ -61,12 +61,10 @@ void NegExpNode::register_python()
 {
     using namespace boost::python;
     class_<NegExpNode, bases<OperatorNode<1>>>("NegExpNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &NegExpNode::type)
         .add_property("rung", &NegExpNode::rung)
     ;
 }
\ 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 02e1f415..c142c131 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
@@ -61,12 +61,10 @@ void SinNode::register_python()
 {
     using namespace boost::python;
     class_<SinNode, bases<OperatorNode<1>>>("SinNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &SinNode::type)
         .add_property("rung", &SinNode::rung)
     ;
 }
\ 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 c66bd122..5a3a0513 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
@@ -49,12 +49,10 @@ void SixPowNode::register_python()
 {
     using namespace boost::python;
     class_<SixPowNode, bases<OperatorNode<1>>>("SixPowNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &SixPowNode::type)
         .add_property("rung", &SixPowNode::rung)
     ;
 }
\ 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 49a83ddc..09f81177 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
@@ -49,12 +49,10 @@ void SqNode::register_python()
 {
     using namespace boost::python;
     class_<SqNode, bases<OperatorNode<1>>>("SqNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &SqNode::type)
         .add_property("rung", &SqNode::rung)
     ;
 }
\ 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 2a5a606c..60fa44e8 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
@@ -49,12 +49,10 @@ void SqrtNode::register_python()
 {
     using namespace boost::python;
     class_<SqrtNode, bases<OperatorNode<1>>>("SqrtNode", init<node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 1>, 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("type", &SqrtNode::type)
         .add_property("rung", &SqrtNode::rung)
     ;
 }
\ 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 2ed49e6c..6ea320c0 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
@@ -80,12 +80,10 @@ void SubNode::register_python()
 {
     using namespace boost::python;
     class_<SubNode, bases<OperatorNode<2>>>("SubNode", init<node_ptr, node_ptr, int, double, double>())
-        .def(init<std::array<node_ptr, 2>, 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("type", &SubNode::type)
         .add_property("rung", &SubNode::rung)
     ;
 }
\ No newline at end of file
diff --git a/src/feature_creation/node/value_storage/nodes_value_containers.cpp b/src/feature_creation/node/value_storage/nodes_value_containers.cpp
index 07706742..654a6c20 100644
--- a/src/feature_creation/node/value_storage/nodes_value_containers.cpp
+++ b/src/feature_creation/node/value_storage/nodes_value_containers.cpp
@@ -1,11 +1,11 @@
 #include <feature_creation/node/value_storage/nodes_value_containers.hpp>
 
-int node_value_arrs::N_SELECTED;
-int node_value_arrs::N_SAMPLES;
-int node_value_arrs::N_STORE_FEATURES;
-int node_value_arrs::N_TEMP_STORE_FEATURES;
-int node_value_arrs::N_RUNGS_STORED;
-int node_value_arrs::N_SAMPLES_TEST;
+int node_value_arrs::N_SELECTED = 0;
+int node_value_arrs::N_SAMPLES = 0;
+int node_value_arrs::N_STORE_FEATURES = 0;
+int node_value_arrs::N_TEMP_STORE_FEATURES = 0;
+int node_value_arrs::N_RUNGS_STORED = 0;
+int node_value_arrs::N_SAMPLES_TEST = 0;
 
 std::vector<int> node_value_arrs::TEMP_STORAGE_REG;
 std::vector<int> node_value_arrs::TEMP_STORAGE_TEST_REG;
diff --git a/src/feature_creation/units/Unit.cpp b/src/feature_creation/units/Unit.cpp
index 3fbcda80..c95599f4 100644
--- a/src/feature_creation/units/Unit.cpp
+++ b/src/feature_creation/units/Unit.cpp
@@ -183,6 +183,7 @@ std::ostream& operator<< (std::ostream& outStream, const Unit& unit)
     return outStream;
 }
 
+
 void Unit::register_python()
 {
     using namespace boost::python;
@@ -200,5 +201,5 @@ void Unit::register_python()
         .def(self == self)
         .def(self != self)
         .def("__pow__", &Unit::operator^)
-        .add_property("dct", &Unit::dct);
+    ;
 }
\ No newline at end of file
diff --git a/src/main.cpp b/src/main.cpp
index 8170f789..b1535fa8 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,8 +6,7 @@ int main(int argc, char const *argv[])
 {
     allowed_op_maps::set_node_maps();
 
-    mpi::environment env;
-    std::shared_ptr<MPI_Interface> comm = std::make_shared<MPI_Interface>();
+    mpi_setup::init_mpi_env();
 
     std::string filename;
     double duration= 0.0;
@@ -16,11 +15,11 @@ int main(int argc, char const *argv[])
     else
         filename = argv[1];
 
-    if(comm->rank() == 0)
+    if(mpi_setup::comm->rank() == 0)
         stripComments(filename);
     else
         filename = "stripped_" + filename;
-    comm->barrier();
+    mpi_setup::comm->barrier();
 
     //construct the parser and pass it to the inputs
     boost::property_tree::ptree propTree;
@@ -28,19 +27,19 @@ int main(int argc, char const *argv[])
     std::clock_t start;
     start = std::clock();
 
-    InputParser IP(propTree, filename, comm);
-    if(comm->rank() == 0)
+    InputParser IP(propTree, filename, mpi_setup::comm);
+    if(mpi_setup::comm->rank() == 0)
         boost::filesystem::remove(filename);
-    comm->barrier();
+    mpi_setup::comm->barrier();
     duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
-    if(comm->rank() == 0)
+    if(mpi_setup::comm->rank() == 0)
         std::cout<< "time input_parsing/Feature space generation: "<< duration << std::endl;
 
     node_value_arrs::initialize_d_matrix_arr();
     SISSORegressor sisso(IP._feat_space, IP._prop_train, IP._prop_test, IP._task_sizes_train, IP._task_sizes_test, IP._n_dim, IP._n_residuals);
     sisso.fit();
 
-    if(comm->rank() == 0)
+    if(mpi_setup::comm->rank() == 0)
     {
         for(int ii = 0; ii < sisso.models().size(); ++ii)
         {
@@ -58,5 +57,7 @@ int main(int argc, char const *argv[])
             }
         }
     }
+
+    mpi_setup::finalize_mpi_env();
     return 0;
 }
diff --git a/src/mpi_interface/MPI_Interface.cpp b/src/mpi_interface/MPI_Interface.cpp
index ddb8b988..6f9ce6b7 100644
--- a/src/mpi_interface/MPI_Interface.cpp
+++ b/src/mpi_interface/MPI_Interface.cpp
@@ -21,6 +21,7 @@ void mpi_setup::init_mpi_env(int& argc, char **&argv)
 {
     if(env == 0)
         env = new boost::mpi::environment(argc, argv);
+    comm = std::make_shared<MPI_Interface>();
 }
 
 void mpi_setup::init_mpi_env()
@@ -29,6 +30,7 @@ void mpi_setup::init_mpi_env()
     {
         #ifdef BOOST_MPI_HAS_NOARG_INITIALIZATION
             env = new boost::mpi::environment();
+            comm = std::make_shared<MPI_Interface>();
         #else
             throw std::runtime_error("MPI cannot be initialized without arguments");
         #endif
diff --git a/src/python/__init__.py b/src/python/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/python/_sisso.cpp b/src/python/_sisso.cpp
index 0209be86..828ef3ee 100644
--- a/src/python/_sisso.cpp
+++ b/src/python/_sisso.cpp
@@ -1,13 +1,25 @@
 #include <python/bindings.hpp>
+
+#include <boost/python/suite/indexing/map_indexing_suite.hpp>
+#include <boost/python/numpy.hpp>
+
 #include <mpi_interface/MPI_Interface.hpp>
 #include <Python.h>
 
+
 static void finalize();
 
 BOOST_PYTHON_MODULE(_sisso)
 {
+    Py_Initialize();
+    np::initialize();
     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> >() );
+
     Py_AtExit(&finalize);
 }
 
diff --git a/src/python/conversion_utils.hpp b/src/python/conversion_utils.hpp
new file mode 100644
index 00000000..182db98d
--- /dev/null
+++ b/src/python/conversion_utils.hpp
@@ -0,0 +1,64 @@
+#ifndef UTILS_CONVERSION
+#define UTILS_CONVERSION
+
+#include <algorithm>
+#include <vector>
+#include <boost/python.hpp>
+#include <boost/python/numpy.hpp>
+
+#include <iostream>
+
+namespace python_conv_utils
+{
+    using namespace boost::python;
+    namespace np = boost::python::numpy;
+
+    template<typename T>
+    std::vector<T> from_list(list lst)
+    {
+        std::vector<T> vec(len(lst));
+        for(int ll = vec.size() - 1; ll >= 0; --ll)
+            vec[ll] = extract<T>(lst.pop());
+        return vec;
+    }
+
+    template<typename T>
+    std::vector<T> from_ndarray(np::ndarray arr)
+    {
+        if(arr.get_dtype() != np::dtype::get_builtin<T>())
+            throw std::logic_error("arr is not of the correct type.");
+
+        std::vector<T> vec(arr.shape(0));
+        std::copy_n(reinterpret_cast<T*>(arr.get_data()), vec.size(), vec.data());
+
+        return vec;
+    }
+
+    template<typename T_ptr, typename T_base>
+    std::vector<std::shared_ptr<T_ptr>> shared_ptr_vec_from_list(list lst)
+    {
+        std::vector<std::shared_ptr<T_ptr>> vec(len(lst));
+        for(int ll = vec.size() - 1; ll >= 0; --ll)
+            vec[ll] = std::make_shared<T_base>(extract<T_base>(lst.pop()));
+        return vec;
+    }
+
+    template<typename T>
+    list to_list(std::vector<T> vec)
+    {
+        list lst;
+        for(auto& item : vec)
+            lst.append<T>(item);
+        return lst;
+    }
+
+    template<typename T>
+    np::ndarray to_ndarray(std::vector<T> vec)
+    {
+        np::ndarray arr = np::zeros(make_tuple(vec.size()), np::dtype::get_builtin<T>());
+        std::copy_n(vec.data(), vec.size(), reinterpret_cast<T*>(arr.get_data()));
+        return arr;
+    }
+}
+
+#endif
\ No newline at end of file
diff --git a/src/sisso++.cmakein b/src/sisso++.cmakein
new file mode 100644
index 00000000..e69de29b
diff --git a/src/src/python/__init__.py b/src/src/python/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/utils/project.cpp b/src/utils/project.cpp
index 9893bac5..4ca64b45 100644
--- a/src/utils/project.cpp
+++ b/src/utils/project.cpp
@@ -3,6 +3,7 @@
 void project_funcs::project_r(double* prop, double* scores, std::vector<node_ptr>& phi, std::vector<int>& sizes, int n_prop)
 {
     int n_samp = std::accumulate(sizes.begin(), sizes.end(), 0);
+
     std::transform(phi.begin(), phi.end(), scores, [&prop, &sizes](node_ptr feat){return -1 * util_funcs::r(prop, feat->value_ptr(), sizes.data(), sizes.size());});
 
     for(int pp = 1; pp < n_prop; ++pp)
-- 
GitLab