FeatureSpace.hpp

// Copyright 2021 Thomas A. R. Purcell
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
/** @file feature_creation/feature_space/FeatureSpace.hpp
 *  @brief Create a feature space from an initial set of features and algebraic operators
 *
 *  Use an initial set of features and combine them to generate more complicated algebraical features. SIS is also performed here
 *
 *  @author Thomas A. R. Purcell (tpurcell)
 *  @bug No known bugs.
 */

#ifndef FEATURE_SPACE
#define FEATURE_SPACE

#include <boost/filesystem.hpp>

#include <utility>

#include "feature_creation/node/utils.hpp"

#include "mpi_interface/MPI_Interface.hpp"
#include "mpi_interface/MPI_ops.hpp"
#include "mpi_interface/serialize_tuple.h"

#include "utils/project.hpp"

#ifdef PY_BINDINGS
    namespace np = boost::python::numpy;
    namespace py = boost::python;
#endif

// DocString: cls_feat_space
/**
 * @brief Feature Space for SISSO calculations
 * @details Stores and performs all feature calculations for SIS
 *
 */
class FeatureSpace
{
    std::vector<node_ptr> _phi_selected; //!< selected features
    std::vector<node_ptr> _phi; //!< all features
    const std::vector<node_ptr> _phi_0; //!< initial feature space

    #ifdef PARAMETERIZE
    std::vector<node_ptr> _phi_reparam; //!< The list of nodes used for reparameterization
    std::vector<int> _end_no_params; //!< The list of indexes of each rung where parameterized nodes start
    std::vector<int> _start_gen_reparam; //!< The list of indexes of each rung where parameterized nodes start

    std::vector<un_param_op_node_gen> _un_param_operators; //!< list of all parameterized unary operators with free parameters
    std::vector<bin_param_op_node_gen> _com_bin_param_operators; //!< list of all parameterized commutable binary operators with free parameters
    std::vector<bin_param_op_node_gen> _bin_param_operators; //!< list of all parameterized binary operators with free parameters
    std::vector<std::string> _allowed_param_ops; //!< Map of parameterization operator set (set of operators and non-linear parameters used for a non-linear least squares fit to property)
    #endif

    std::vector<std::string> _allowed_ops; //!< list of all allowed operators strings
    std::vector<un_op_node_gen> _un_operators; //!< list of all unary operators
    std::vector<bin_op_node_gen> _com_bin_operators; //!< list of all commutable binary operators
    std::vector<bin_op_node_gen> _bin_operators; //!< list of all binary operators

    std::vector<double> _prop; //!< The property to fit
    std::vector<double> _scores; //!< projection scores for each feature

    const std::vector<int> _task_sizes; //!< The number of elements in each task (training data)
    std::vector<int> _start_gen; //!< list of the indexes where each generation starts in _phi
    const std::string _project_type; //!< The type of projection that should be done during SIS
    const std::string _feature_space_file; //!< File to store information about the selected features
    const std::string _feature_space_summary_file; //!< File to store information about the selected features

    std::function<bool(const double*, const int, const double, const std::vector<double>&, const double, const int, const int)> _is_valid; //!< Function used to calculate the scores for SIS
    std::function<bool(const double*, const int, const double, const std::vector<node_ptr>&, const std::vector<double>&, const double)> _is_valid_feat_list; //!< Function used to calculate the scores for SIS without changing omp environment

    std::shared_ptr<MPI_Interface> _mpi_comm; //!< MPI communicator

    const double _cross_cor_max; //!< Maximum cross-correlation used for selecting features
    const double _l_bound; //!< lower bound for absolute value of the features
    const double _u_bound; //!< upper bound for absolute value of the features

    int _n_rung_store; //!< Total rungs stored
    int _n_feat; //!< Total number of features
    int _max_phi; //!< Maximum rung for the feature creation

    const int _n_sis_select; //!< Number of features to select for each dimensions
    const int _n_samp; //!< Number of samples (training data)
    const int _n_rung_generate; //!< Total number of rungs to generate on the fly

    int _max_param_depth; //!< Max depth to parameterize a feature (default=_max_rung)
    const bool _reparam_residual; //!< If True then reparameterize using the residuals of each model

public:

    #ifdef PARAMETERIZE
    /**
     * @brief Constructor for the feature space
     * @details constructs the feature space from an initial set of features and a list of allowed operators
     *
     * @param mpi_comm MPI communicator for the calculations
     * @param phi_0 The initial set of features to combine
     * @param allowed_ops list of allowed operators
     * @param allowed_param_ops dictionary of the parameterizable operators and their associated free parameters
     * @param prop The property to be learned (training data)
     * @param task_sizes The number of samples per task
     * @param project_type The projection operator to use
     * @param max_phi highest rung value for the calculation
     * @param n_sis_select number of features to select during each SIS step
     * @param max_store_rung number of rungs to calculate and store the value of the features for all samples
     * @param n_rung_generate number of rungs to generate on the fly during SIS (this must be 1 or 0 right now, possible to be higher with recursive algorithm)
     * @param cross_corr_max Maximum cross-correlation used for selecting features
     * @param min_abs_feat_val minimum absolute feature value
     * @param max_abs_feat_val maximum absolute feature value
     * @param max_param_depth the maximum paremterization depths for features
     * @param reparam_residual If True then reparameterize using the residuals of each model
     */
    FeatureSpace(
        std::shared_ptr<MPI_Interface> mpi_comm,
        std::vector<node_ptr> phi_0,
        std::vector<std::string> allowed_ops,
        std::vector<std::string> allowed_param_ops,
        std::vector<double> prop,
        std::vector<int> task_sizes,
        std::string project_type="regression",
        int max_phi=1,
        int n_sis_select=1,
        int max_store_rung=-1,
        int n_rung_generate=0,
        double cross_corr_max=1.0,
        double min_abs_feat_val=1e-50,
        double max_abs_feat_val=1e50,
        int max_param_depth=-1,
        bool reparam_residual=false
    );
    #else
    /**
     * @brief Constructor for the feature space
     * @details constructs the feature space from an initial set of features and a list of allowed operators
     *
     * @param mpi_comm MPI communicator for the calculations
     * @param phi_0 The initial set of features to combine
     * @param allowed_ops list of allowed operators
     * @param prop The property to be learned (training data)
     * @param task_sizes The number of samples per task
     * @param project_type The projection operator to use
     * @param max_phi highest rung value for the calculation
     * @param n_sis_select number of features to select during each SIS step
     * @param max_store_rung number of rungs to calculate and store the value of the features for all samples
     * @param n_rung_generate number of rungs to generate on the fly during SIS (this must be 1 or 0 right now, possible to be higher with recursive algorithm)
     * @param cross_corr_max Maximum cross-correlation used for selecting features
     * @param min_abs_feat_val minimum absolute feature value
     * @param max_abs_feat_val maximum absolute feature value
     */
    FeatureSpace(
        std::shared_ptr<MPI_Interface> mpi_comm,
        std::vector<node_ptr> phi_0,
        std::vector<std::string> allowed_ops,
        std::vector<double> prop,
        std::vector<int> task_sizes,
        std::string project_type="regression",
        int max_phi=1,
        int n_sis_select=1,
        int max_store_rung=-1,
        int n_rung_generate=0,
        double cross_corr_max=1.0,
        double min_abs_feat_val=1e-50,
        double max_abs_feat_val=1e50
    );
    #endif
    /**
     * @brief Initialize the feature set given a property vector
     */
    void initialize_fs();

    /**
     * @brief Uses _allowed_ops to set the operator lists
     */
    void set_op_lists();

    /**
     * @brief Initializes the output files for SIS
     */
    void initialize_fs_output_files() const;
    /**
     * @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
     */
    void generate_feature_space();

    /**
     * @brief The selected feature space
     */
    inline std::vector<node_ptr> phi_selected() const {return _phi_selected;};

    /**
     * @brief The full feature space
     */
    inline std::vector<node_ptr> phi() const {return _phi;};

    /**
     * @brief The initial feature space
     */
    inline std::vector<node_ptr> phi0() const {return _phi_0;};

    /**
     * @brief The vector of projection scores for SIS
     */
    inline std::vector<double> scores() const {return _scores;}

    /**
     * @brief The MPI Communicator
     */
    inline std::shared_ptr<MPI_Interface> mpi_comm() const {return _mpi_comm;}

    /**
     * @brief The vector storing the number of samples in each task
     */
    inline std::vector<int> task_sizes() const {return _task_sizes;}

    // DocString: feat_space_feature_space_file
    /**
     * @brief The feature space filename
     */
    inline std::string feature_space_file() const {return _feature_space_file;}

    // DocString: feat_space_l_bound
    /**
     * @brief The minimum absolute value of the feature
     */
    inline double l_bound() const {return _l_bound;}

    // DocString: feat_space_u_bound
    /**
     * @brief The maximum absolute value of the feature
     */
    inline double u_bound() const {return _u_bound;}

    // DocString: feat_space_max_phi
    /**
     * @brief The maximum rung of the feature space
     */
    inline int max_phi() const {return _max_phi;}

    // DocString: feat_space_n_sis_select
    /**
     * @brief The number of features selected in each SIS step
     */
    inline int n_sis_select() const {return _n_sis_select;}

    // DocString: feat_space_n_samp
    /**
     * @brief The number of samples per feature
     */
    inline int n_samp() const {return _n_samp;}

    // DocString: feat_space_n_feat
    /**
     * @brief The number of features in the feature space
     */
    inline int n_feat() const {return _n_feat;}

    // DocString: feat_space_n_rung_store
    /**
     * @brief The number of rungs whose feature's data is always stored in memory
     */
    inline int n_rung_store() const {return _n_rung_store;}

    // DocString: feat_space_n_rung_generate
    /**
     * @brief Either 0 or 1, and is the number of rungs to generate on the fly during SIS
     */
    inline int n_rung_generate() const {return _n_rung_generate;}

    /**
     * @brief Generate a new set of non-parameterized features from a single feature
     * @details Take in the feature and perform all valid algebraic operations on it.
     *
     * @param feat The feature to spawn new features from
     * @param feat_set The feature set to pull features from for combinations
     * @param feat_ind starting index for the next feature generated
     * @param optimizer The object used to optimize the parameterized features
     * @param l_bound lower bound for the absolute value of the feature
     * @param u_bound upper bound for the abosulte value of the feature
     */
    void generate_non_param_feats(
        std::vector<node_ptr>::iterator& feat,
        std::vector<node_ptr>& feat_set,
        const std::vector<node_ptr>::iterator& start,
        unsigned long int& feat_ind,
        const double l_bound=1e-50,
        const double u_bound=1e50
    );

#ifdef PARAMETERIZE
    /**
     * @brief Generate a new set of parameterized features from a single feature
     * @details Take in the feature and perform all valid algebraic operations on it.
     *
     * @param feat The feature to spawn new features from
     * @param feat_set The feature set to pull features from for combinations
     * @param feat_ind starting index for the next feature generated
     * @param optimizer The object used to optimize the parameterized features
     * @param l_bound lower bound for the absolute value of the feature
     * @param u_bound upper bound for the abosulte value of the feature
     */
    void generate_param_feats(
        std::vector<node_ptr>::iterator& feat,
        std::vector<node_ptr>& feat_set,
        const std::vector<node_ptr>::iterator& start,
        unsigned long int& feat_ind,
        std::shared_ptr<NLOptimizer> optimizer,
        const double l_bound=1e-50,
        const double u_bound=1e50
    );

    /**
     * @brief Generate a new set of parameterized features for the residuals
     *
     * @param feat The feature to spawn new features from
     * @param feat_set The feature set to pull features from for combinations
     * @param feat_ind starting index for the next feature generated
     * @param optimizer The object used to optimize the parameterized features
     * @param l_bound lower bound for the absolute value of the feature
     * @param u_bound upper bound for the abosulte value of the feature
     */
    void generate_reparam_feats(
        std::vector<node_ptr>::iterator& feat,
        std::vector<node_ptr>& feat_set,
        unsigned long int& feat_ind,
        std::shared_ptr<NLOptimizer> optimizer,
        const double l_bound=1e-50,
        const double u_bound=1e50
    );

    /**
     * @brief Generate reparameterized feature set
     *
     * @param prop The property to optimize against
     */
    void generate_reparam_feature_set(const std::vector<double>& prop);
#endif

    /**
     * @brief Calculate the SIS Scores for feature generated on the fly
     * @details Create the next rung of features and calculate their projection scores. Only keep those that can be selected by SIS.
     *
     * @param prop Pointer to the start of the vector storing the data to project the features onto
     * @param size The size of the data to project over
     * @param phi_selected The features that would be selected from the previous rungs
     * @param scores_selected The projection scores of the features that would be selected from the previous rungs
     */
    void project_generated(std::shared_ptr<LossFunction> loss, std::vector<node_ptr>& phi_selected, std::vector<double>& scores_selected);

    /**
     * @brief Perform Sure-Independence Screening over the FeatureSpace. The features are ranked using a projection operator constructed using _project_type and the Property vector
     *
     * @param prop Vector containing the property vector (training data only)
     */
    void sis(const std::vector<double>& prop);

    /**
     * @brief Perform Sure-Independence Screening over the FeatureSpace. The features are ranked using a projection operator defined in loss
     *
     * @param loss The LossFunction used to project over all of the features
     */
    void sis(std::shared_ptr<LossFunction> loss);

    // DocString: feat_space_feat_in_phi
    /**
     * @brief Is a feature in this process' _phi?
     *
     * @param ind (int) The index of the feature
     *
     * @return True if feature is in this rank's _phi
     */
    inline bool feat_in_phi(int ind) const {return (ind >= _phi[0]->feat_ind()) && (ind <= _phi.back()->feat_ind());}

    // DocString: feat_space_remove_feature
    /**
     * @brief Remove a feature from phi
     *
     * @param ind (int) index of feature to remove
     */
    void remove_feature(const int ind);

    // Python Interface Functions
    #ifdef PY_BINDINGS
    #ifdef PARAMETERIZE

    // DocString: feat_space_init_py_list
    /**
     * @brief FeatureSpace constructor given a set of primary features and operators
     *
     * @param phi_0 (list) The set of primary features
     * @param allowed_ops (list) The list of allowed operators
     * @param allowed_param_ops (list) The list of allowed operators to be used with non-linear optimization
     * @param prop (list) List containing the property vector (training data only)
     * @param project_type (str) The type of loss function/projection operator to use
     * @param max_phi (int) The maximum rung of the feature (Height of the binary expression tree -1)
     * @param n_sis_select (int) The number of features to select during each SIS step
     * @param max_store_rung (int) The number of rungs whose feature's data is always stored in memory
     * @param n_rung_generate (int) Either 0 or 1, and is the number of rungs to generate on the fly during SIS
     * @param cross_corr_max (double) The maximum allowed cross-correlation value between selected features
     * @param min_abs_feat_val (double) The minimum allowed absolute feature value for a feature
     * @param max_abs_feat_val (double) The maximum allowed absolute feature value for a feature
     * @param max_param_depth (int) The maximum depth in the binary expression tree to set non-linear optimization
     * @param reparam_residual (bool) If True then reparameterize features using the residuals of each model
     */
    FeatureSpace(
        py::list phi_0,
        py::list allowed_ops,
        py::list allowed_param_ops,
        py::list prop,
        std::string project_type="regression",
        int max_phi=1,
        int n_sis_select=1,
        int max_store_rung=-1,
        int n_rung_generate=0,
        double cross_corr_max=1.0,
        double min_abs_feat_val=1e-50,
        double max_abs_feat_val=1e50,
        int max_param_depth = -1,
        bool reparam_residual=false
    );

    // DocString: feat_space_init_np_array
    /**
     * @brief FeatureSpace constructor given a set of primary features and operators
     *
     * @param phi_0 (list) The set of primary features
     * @param allowed_ops (list) The list of allowed operators
     * @param allowed_param_ops (list) The list of allowed operators to be used with non-linear optimization
     * @param prop (np.ndarray) List containing the property vector (training data only)
     * @param project_type (str) The type of loss function/projection operator to use
     * @param max_phi (int) The maximum rung of the feature (Height of the binary expression tree -1)
     * @param n_sis_select (int) The number of features to select during each SIS step
     * @param max_store_rung (int) The number of rungs whose feature's data is always stored in memory
     * @param n_rung_generate (int) Either 0 or 1, and is the number of rungs to generate on the fly during SIS
     * @param cross_corr_max (double) The maximum allowed cross-correlation value between selected features
     * @param min_abs_feat_val (double) The minimum allowed absolute feature value for a feature
     * @param max_abs_feat_val (double) The maximum allowed absolute feature value for a feature
     * @param max_param_depth (int) The maximum depth in the binary expression tree to set non-linear optimization
     * @param reparam_residual (bool) If True then reparameterize features using the residuals of each model
     */
    FeatureSpace(
        py::list phi_0,
        py::list allowed_ops,
        py::list allowed_param_ops,
        np::ndarray prop,
        std::string project_type="regression",
        int max_phi=1,
        int n_sis_select=1,
        int max_store_rung=-1,
        int n_rung_generate=0,
        double cross_corr_max=1.0,
        double min_abs_feat_val=1e-50,
        double max_abs_feat_val=1e50,
        int max_param_depth = -1,
        bool reparam_residual=false
    );

    #else

    // DocString: feat_space_ini_no_param_py_list
    /**
     * @brief FeatureSpace constructor given a set of primary features and operators
     *
     * @param phi_0 (list) The set of primary features
     * @param allowed_ops (list) The list of allowed operators
     * @param prop (list) List containing the property vector (training data only)
     * @param project_type (str) The type of loss function/projection operator to use
     * @param max_phi (int) The maximum rung of the feature (Height of the binary expression tree -1)
     * @param n_sis_select (int) The number of features to select during each SIS step
     * @param max_store_rung (int) The number of rungs whose feature's data is always stored in memory
     * @param n_rung_generate (int) Either 0 or 1, and is the number of rungs to generate on the fly during SIS
     * @param cross_corr_max (double) The maximum allowed cross-correlation value between selected features
     * @param min_abs_feat_val (double) The minimum allowed absolute feature value for a feature
     * @param max_abs_feat_val (double) The maximum allowed absolute feature value for a feature
     */
    FeatureSpace(
        py::list phi_0,
        py::list allowed_ops,
        py::list prop,
        std::string project_type="regression",
        int max_phi=1,
        int n_sis_select=1,
        int max_store_rung=-1,
        int n_rung_generate=0,
        double cross_corr_max=1.0,
        double min_abs_feat_val=1e-50,
        double max_abs_feat_val=1e50
    );

    // DocString: feat_space_init_no_param_np_array
    /**
     * @brief FeatureSpace constructor given a set of primary features and operators
     *
     * @param phi_0 (list) The set of primary features
     * @param allowed_ops (list) The list of allowed operators
     * @param prop (np.ndarray) List containing the property vector (training data only)
     * @param project_type (str) The type of loss function/projection operator to use
     * @param max_phi (int) The maximum rung of the feature (Height of the binary expression tree -1)
     * @param n_sis_select (int) The number of features to select during each SIS step
     * @param max_store_rung (int) The number of rungs whose feature's data is always stored in memory
     * @param n_rung_generate (int) Either 0 or 1, and is the number of rungs to generate on the fly during SIS
     * @param cross_corr_max (double) The maximum allowed cross-correlation value between selected features
     * @param min_abs_feat_val (double) The minimum allowed absolute feature value for a feature
     * @param max_abs_feat_val (double) The maximum allowed absolute feature value for a feature
     */
    FeatureSpace(
        py::list phi_0,
        py::list allowed_ops,
        np::ndarray prop,
        std::string project_type="regression",
        int max_phi=1,
        int n_sis_select=1,
        int max_store_rung=-1,
        int n_rung_generate=0,
        double cross_corr_max=1.0,
        double min_abs_feat_val=1e-50,
        double max_abs_feat_val=1e50
    );
    #endif

    // DocString: feat_space_init_file_np_array
    /**
     * @brief FeatureSpace constructor that uses a file containing postfix feature expressions to describe all features in Phi, and a primary feature setn <python/feature_creation/FeatureSpace.cpp>)
     *
     * @param feature_file (str) The file containing the postfix expressions of all features in the FeatureSpace
     * @param phi_0 (list) The set of primary features
     * @param prop (np.ndarray) List containing the property vector (training data only)
     * @param task_sizes (list) The number of samples in the training data per task
     * @param project_type (str) The type of loss function/projection operator to use
     * @param n_sis_select (int) The number of features to select during each SIS step
     * @param cross_corr_max (double) The maximum allowed cross-correlation value between selected features
     */
    FeatureSpace(
        std::string feature_file,
        py::list phi_0,
        np::ndarray prop,
        py::list task_sizes,
        std::string project_type="regression",
        int n_sis_select=1,
        double cross_corr_max=1.0
    );

    // DocString: feat_space_init_file_py_list
    /**
     * @brief FeatureSpace constructor that uses a file containing postfix feature expressions to describe all features in Phi, and a primary feature setn <python/feature_creation/FeatureSpace.cpp>)
     *
     * @param feature_file (str) The file containing the postfix expressions of all features in the FeatureSpace
     * @param phi_0 (list) The set of primary features
     * @param prop (list) List containing the property vector (training data only)
     * @param task_sizes (list) The number of samples in the training data per task
     * @param project_type (str) The type of loss function/projection operator to use
     * @param n_sis_select (int) The number of features to select during each SIS step
     * @param cross_corr_max (double) The maximum allowed cross-correlation value between selected features
     */
    FeatureSpace(
        std::string feature_file,
        py::list phi_0,
        py::list prop,
        py::list task_sizes,
        std::string project_type="regression",
        int n_sis_select=1,
        double cross_corr_max=1.0
    );

    // DocString: feat_space_sis_arr
    /**
     * @brief Perform Sure-Independence Screening over the FeatureSpace. The features are ranked using a projection operator constructed using _project_type and the Property vector
     *
     * @param prop (np.ndarray) Array containing the property vector (training data only)
     */
    inline void sis(np::ndarray prop)
    {
        std::vector<double> prop_vec = python_conv_utils::from_ndarray<double>(prop);
        sis(prop_vec);
    }

    // DocString: feat_space_sis_list
    /**
     * @brief Perform Sure-Independence Screening over the FeatureSpace. The features are ranked using a projection operator constructed using _project_type and the Property vector
     *
     * @param prop (list) List containing the property vector (training data only)
     */
    inline void sis(py::list prop)
    {
        std::vector<double> prop_vec = python_conv_utils::from_list<double>(prop);
        sis(prop_vec);
    }

    // DocString: feat_space_phi_selected_py
    /**
     * @brief A list of all of the selected features
     */
    py::list phi_selected_py();

    // DocString: feat_space_phi0_py
    /**
     * @brief A list containing all of the primary features
     */
    py::list phi0_py();

    // DocString: feat_space_phi_py
    /**
     * @brief A list of all features in the FeatureSpace
     */
    py::list phi_py();

    // DocString: feat_space_scores_py
    /**
     * @brief An array of all stored projection scores from SIS
     */
    inline np::ndarray scores_py(){return python_conv_utils::to_ndarray<double>(_scores);};

    // DocString: feat_space_task_sizes_py
    /**
     * @brief A list of the number of samples in each task for the training data
     */
    inline py::list task_sizes_py(){return python_conv_utils::to_list<int>(_task_sizes);};

    // DocString: feat_space_allowed_ops_py
    /**
     * @brief The list of allowed operators
     */
    inline py::list allowed_ops_py(){return python_conv_utils::to_list<std::string>(_allowed_ops);}

    // DocString: feat_space_start_gen_py
    /**
     * @brief A list containing the index of the first feature of each rung in the feature space.
     */
    inline py::list start_gen_py(){return python_conv_utils::to_list<int>(_start_gen);}

    // DocString: feat_space_get_feature
    /**
     * @brief Access the feature in _phi with an index ind
     *
     * @param ind (int) The index of the feature to get
     * @return A ModelNode of the feature at index ind
     */
    inline ModelNode get_feature(const int ind) const {return ModelNode(_phi[ind]);}
    #endif
};

#endif