FeatureSpace.cpp

#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")
BOOST_CLASS_EXPORT_GUID(AbsDiffNode, "AbsDiffNode")
BOOST_CLASS_EXPORT_GUID(MultNode, "MultNode")
BOOST_CLASS_EXPORT_GUID(DivNode, "DivNode")
BOOST_CLASS_EXPORT_GUID(SqNode, "SqNode")
BOOST_CLASS_EXPORT_GUID(SqrtNode, "SqrtNode")
BOOST_CLASS_EXPORT_GUID(CbNode, "CbNode")
BOOST_CLASS_EXPORT_GUID(CbrtNode, "CbrtNode")
BOOST_CLASS_EXPORT_GUID(SixPowNode, "SixPowNode")
BOOST_CLASS_EXPORT_GUID(ExpNode, "ExpNode")
BOOST_CLASS_EXPORT_GUID(NegExpNode, "NegExpNode")
BOOST_CLASS_EXPORT_GUID(LogNode, "LogNode")
BOOST_CLASS_EXPORT_GUID(AbsNode, "AbsNode")
BOOST_CLASS_EXPORT_GUID(InvNode, "InvNode")
BOOST_CLASS_EXPORT_GUID(SinNode, "SinNode")
BOOST_CLASS_EXPORT_GUID(CosNode, "CosNode")

#ifdef PARAMETERIZE
BOOST_CLASS_EXPORT_GUID(AddParamNode, "AddParamNode")
BOOST_CLASS_EXPORT_GUID(SubParamNode, "SubParamNode")
BOOST_CLASS_EXPORT_GUID(AbsDiffParamNode, "AbsDiffParamNode")
BOOST_CLASS_EXPORT_GUID(MultParamNode, "MultParamNode")
BOOST_CLASS_EXPORT_GUID(DivParamNode, "DivParamNode")
BOOST_CLASS_EXPORT_GUID(SqParamNode, "SqParamNode")
BOOST_CLASS_EXPORT_GUID(SqrtParamNode, "SqrtParamNode")
BOOST_CLASS_EXPORT_GUID(CbParamNode, "CbParamNode")
BOOST_CLASS_EXPORT_GUID(CbrtParamNode, "CbrtParamNode")
BOOST_CLASS_EXPORT_GUID(SixPowParamNode, "SixPowParamNode")
BOOST_CLASS_EXPORT_GUID(ExpParamNode, "ExpParamNode")
BOOST_CLASS_EXPORT_GUID(NegExpParamNode, "NegExpParamNode")
BOOST_CLASS_EXPORT_GUID(LogParamNode, "LogParamNode")
BOOST_CLASS_EXPORT_GUID(AbsParamNode, "AbsParamNode")
BOOST_CLASS_EXPORT_GUID(InvParamNode, "InvParamNode")
BOOST_CLASS_EXPORT_GUID(SinParamNode, "SinParamNode")
BOOST_CLASS_EXPORT_GUID(CosParamNode, "CosParamNode")

FeatureSpace::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,
    int max_phi,
    int n_sis_select,
    int max_store_rung,
    int n_rung_generate,
    double cross_corr_max,
    double min_abs_feat_val,
    double max_abs_feat_val,
    int max_param_depth,
    bool reparam_residual
):
    _phi(phi_0),
    _phi_0(phi_0),
    _allowed_param_ops(allowed_param_ops),
    _allowed_ops(allowed_ops),
    _prop(prop),
    _scores(phi_0.size(), 0.0),
    _task_sizes(task_sizes),
    _start_gen(1, 0),
    _feature_space_file("feature_space/selected_features.txt"),
    _feature_space_summary_file("feature_space/SIS_summary.txt"),
    _project_type(project_type),
    _mpi_comm(mpi_comm),
    _cross_cor_max(cross_corr_max),
    _l_bound(min_abs_feat_val),
    _u_bound(max_abs_feat_val),
    _max_phi(max_phi),
    _n_sis_select(n_sis_select),
    _n_samp(phi_0[0]->n_samp()),
    _n_feat(phi_0.size()),
    _n_rung_store(max_store_rung),
    _n_rung_generate(n_rung_generate),
    _max_param_depth(max_param_depth),
    _reparam_residual(reparam_residual)
{
    initialize_fs();
}
#else
FeatureSpace::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,
    int max_phi,
    int n_sis_select,
    int max_store_rung,
    int n_rung_generate,
    double cross_corr_max,
    double min_abs_feat_val,
    double max_abs_feat_val
):
    _phi(phi_0),
    _phi_0(phi_0),
    _allowed_ops(allowed_ops),
    _prop(prop),
    _scores(phi_0.size(), 0.0),
    _task_sizes(task_sizes),
    _start_gen(1, 0),
    _feature_space_file("feature_space/selected_features.txt"),
    _feature_space_summary_file("feature_space/SIS_summary.txt"),
    _project_type(project_type),
    _mpi_comm(mpi_comm),
    _cross_cor_max(cross_corr_max),
    _l_bound(min_abs_feat_val),
    _u_bound(max_abs_feat_val),
    _max_phi(max_phi),
    _n_sis_select(n_sis_select),
    _n_samp(phi_0[0]->n_samp()),
    _n_feat(phi_0.size()),
    _n_rung_store(max_store_rung),
    _n_rung_generate(n_rung_generate),
    _max_param_depth(0),
    _reparam_residual(false)
{
    initialize_fs();
}
#endif

void FeatureSpace::initialize_fs()
{
    #ifdef PARAMETERIZE
    if(_max_param_depth == -1)
    {
        _max_param_depth = _max_phi;
    }
    if((_max_param_depth < 0) || (_max_param_depth > _max_phi))
    {
        throw std::logic_error("Invalid parameter depth.");
    }
    nlopt_wrapper::MAX_PARAM_DEPTH = _max_param_depth;
    #endif

    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.");
    }

    node_value_arrs::set_task_sz_train(_task_sizes);
    int n_max_ops = 0;
    for(int rr = 0; rr < _max_phi - _n_rung_store; ++rr)
    {
        n_max_ops += std::pow(2, rr);
    }

    initialize_fs_output_files();
    project_funcs::set_project_fxn(_project_type, _task_sizes.size(), _project, _project_no_omp);
    comp_feats::set_is_valid_fxn(_project_type, _cross_cor_max, _n_samp, _is_valid, _is_valid_feat_list);
    set_op_lists();

    double start = omp_get_wtime();
    generate_feature_space();
    _mpi_comm->barrier();
    double duration = omp_get_wtime() - start;
    if(_mpi_comm->rank() == 0)
    {
        std::cout << "time to generate feat sapce: " << duration << " s" << std::endl;
    }

    mpi_reduce_op::set_op(_project_type, _cross_cor_max, _n_sis_select);

    _scores.reserve(_phi.size());
    _scores.resize(_phi.size());
}

void FeatureSpace::set_op_lists()
{
    for(auto & op : _allowed_ops)
    {
        if((op.compare("add") == 0) || (op.compare("mult") == 0) || (op.compare("abs_diff") == 0) || (op.compare("sub") == 0))
        {
            _com_bin_operators.push_back(allowed_op_maps::binary_operator_map[op]);
        }
        else if((op.compare("div") == 0))
        {
            _bin_operators.push_back(allowed_op_maps::binary_operator_map[op]);
        }
        else
        {
            _un_operators.push_back(allowed_op_maps::unary_operator_map[op]);
        }
    }

    #ifdef PARAMETERIZE
    for(auto& op : _allowed_param_ops)
    {
        if((op.compare("add") == 0) || (op.compare("sub") == 0) || (op.compare("abs_diff") == 0))
        {
            _com_bin_param_operators.push_back(allowed_op_maps::binary_param_operator_map[op]);
        }
        else if((op.compare("div") == 0) || (op.compare("mult") == 0))
        {
            _bin_param_operators.push_back(allowed_op_maps::binary_param_operator_map[op]);
        }
        else
        {
            _un_param_operators.push_back(allowed_op_maps::unary_param_operator_map[op]);
        }
    }
    #endif
}

void FeatureSpace::initialize_fs_output_files() const
{
    if(_mpi_comm->rank() == 0)
    {
        std::ofstream out_file_stream = std::ofstream();
        std::ofstream sum_file_stream = std::ofstream();
        out_file_stream.open(_feature_space_file);
        sum_file_stream.open(_feature_space_summary_file);
        out_file_stream << std::setw(14) <<std::left << "# FEAT_ID" << "Feature Postfix Expression (RPN)" << std::endl;
        sum_file_stream << std::setw(14) <<std::left << "# FEAT_ID" << std::setw(24) << std::left << "Score" << "Feature Expression" << std::endl;
        out_file_stream.close();
        sum_file_stream.close();
    }
}

#ifdef PARAMETERIZE
void FeatureSpace::generate_new_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,
    const double u_bound
)
{

    unsigned long int phi_ind = feat - _phi.begin();
    feat_set.reserve(feat_set.size() + _un_operators.size() + phi_ind * (_com_bin_operators.size() + 2 * _bin_operators.size()));

    for(auto& op : _un_operators)
    {
        op(feat_set, *feat, feat_ind, l_bound, u_bound);
    }

    for(auto& op : _com_bin_operators)
    {
        for(auto feat_2 = _phi.begin(); feat_2 < feat; ++feat_2)
        {
            op(feat_set, *feat, *feat_2, feat_ind, l_bound, u_bound);
        }
    }

    for(auto& op : _bin_operators)
    {
        for(auto feat_2 = _phi.begin(); feat_2 < feat; ++feat_2)
        {
            op(feat_set, *feat, *feat_2, feat_ind, l_bound, u_bound);
            op(feat_set, *feat_2, *feat, feat_ind, l_bound, u_bound);
        }
    }
    for(auto& op : _un_param_operators)
    {
        op(feat_set, *feat, feat_ind, l_bound, u_bound, optimizer);
    }

    for(auto& op : _com_bin_param_operators)
    {
        for(auto feat_2 = _phi.begin(); feat_2 != feat; ++feat_2)
        {
            op(feat_set, *feat, *feat_2, feat_ind, l_bound, u_bound, optimizer);
        }
    }
    for(auto& op : _bin_param_operators)
    {
        for(auto feat_2 = _phi.begin(); feat_2 != feat; ++feat_2)
        {
            op(feat_set, *feat, *feat_2, feat_ind, l_bound, u_bound, optimizer);
            op(feat_set, *feat_2, *feat, feat_ind, l_bound, u_bound, optimizer);
        }
    }
}
#else
void FeatureSpace::generate_new_feats(
    std::vector<node_ptr>::iterator& feat,
    std::vector<node_ptr>& feat_set,
    unsigned long int& feat_ind,
    const double l_bound,
    const double u_bound
)
{
    unsigned long int phi_ind = feat - _phi.begin();
    feat_set.reserve(feat_set.size() + _un_operators.size() + phi_ind * (_com_bin_operators.size() + 2 * _bin_operators.size()));

    for(auto& op : _un_operators)
    {
        op(feat_set, *feat, feat_ind, l_bound, u_bound);
    }

    for(auto& op : _com_bin_operators)
    {
        for(auto feat_2 = _phi.begin(); feat_2 < feat; ++feat_2)
        {
            op(feat_set, *feat, *feat_2, feat_ind, l_bound, u_bound);
        }
    }

    for(auto& op : _bin_operators)
    {
        for(auto feat_2 = _phi.begin(); feat_2 < feat; ++feat_2)
        {
            op(feat_set, *feat, *feat_2, feat_ind, l_bound, u_bound);
            op(feat_set, *feat_2, *feat, feat_ind, l_bound, u_bound);
        }
    }
}
#endif

void FeatureSpace::generate_feature_space()
{
    double u_bound = 1e50;
    double l_bound = 1e-50;
    std::vector<int> inds;

    for(int nn = 1; nn <= _max_phi - _n_rung_generate; ++nn)
    {
        node_value_arrs::clear_temp_reg();
        if(nn == _max_phi)
        {
            u_bound = _u_bound;
            l_bound = _l_bound;
        }
        std::vector<node_ptr> next_phi;
        _n_feat = _phi.size();

        unsigned long int feat_ind = _phi.size();

        node_value_arrs::clear_temp_reg();
        double start = omp_get_wtime();
        #ifdef PARAMETERIZE
        #pragma omp parallel firstprivate(feat_ind, l_bound, u_bound)
        {
            std::vector<node_ptr> next_phi_private;
            std::shared_ptr<NLOptimizer> optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, _prop, _max_phi, _max_param_depth);

            #pragma omp for schedule(dynamic)
            for(auto feat_1 = _phi.begin() + _start_gen.back() + _mpi_comm->rank(); feat_1 < _phi.end(); feat_1 += _mpi_comm->size())
            {
                generate_new_feats(feat_1, next_phi_private, feat_ind, optimizer, l_bound, u_bound);
            }

            #pragma omp critical
            next_phi.insert(next_phi.end(), next_phi_private.begin(), next_phi_private.end());
        }
        #else
        #pragma omp parallel firstprivate(feat_ind, l_bound, u_bound)
        {
            std::vector<node_ptr> next_phi_private;
            #pragma omp for schedule(dynamic)
            for(auto feat_1 = _phi.begin() + _start_gen.back() + _mpi_comm->rank(); feat_1 < _phi.end(); feat_1 += _mpi_comm->size())
            {
                generate_new_feats(feat_1, next_phi_private, feat_ind, l_bound, u_bound);
            }

            #pragma omp critical
            next_phi.insert(next_phi.end(), next_phi_private.begin(), next_phi_private.end());
        }
        #endif
        _start_gen.push_back(_phi.size());
        node_value_arrs::clear_temp_reg();
        if((nn < _max_phi) || (nn <= _n_rung_store) || (_mpi_comm->size() == 1))
        {
            int new_phi_size;
            int phi_size_start = _phi.size();
            if(_mpi_comm->rank() == 0)
            {
                std::vector<std::vector<node_ptr>> next_phi_gathered;
                mpi::gather(*_mpi_comm, next_phi, next_phi_gathered, 0);
                feat_ind = _phi.size();
                for(auto& next_phi_vec : next_phi_gathered)
                {
                    _phi.insert(_phi.end(), next_phi_vec.begin(), next_phi_vec.end());
                }
                new_phi_size = _phi.size();

                // Sort the features to ensure consistent feature spaces for all MPI/OpenMP configurations
                std::sort(
                    _phi.begin() + _start_gen.back(),
                    _phi.end(),
                    [feat_ind](node_ptr n1, node_ptr n2){return n1->sort_score(feat_ind) < n2->sort_score(feat_ind);}
                );

                // Reindex sorted features
                std::for_each(
                    _phi.begin() + _start_gen.back(),
                    _phi.end(),
                    [&feat_ind](node_ptr n){n->reindex(feat_ind); ++feat_ind;}
                );

                mpi::broadcast(*_mpi_comm, new_phi_size, 0);

                for(int bb = 0; bb <= (new_phi_size - phi_size_start) / 10000; ++bb)
                {
                    mpi::broadcast(*_mpi_comm, &_phi[phi_size_start + bb * 10000], std::min(10000, new_phi_size - phi_size_start - bb * 10000), 0);
                }
            }
            else
            {
                mpi::gather(*_mpi_comm, next_phi, 0);
                mpi::broadcast(*_mpi_comm, new_phi_size, 0);
                _phi.resize(new_phi_size);

                for(int bb = 0; bb <= (new_phi_size - phi_size_start) / 10000; ++bb)
                {
                    mpi::broadcast(*_mpi_comm, &_phi[phi_size_start + bb * 10000], std::min(10000, new_phi_size - phi_size_start - bb * 10000), 0);
                }
            }

            if(phi_size_start == new_phi_size)
            {
                throw std::logic_error("No features created during this rung (" + std::to_string(nn) + ")");
            }

            node_value_arrs::clear_temp_reg();
            if(nn < _max_phi)
            {
                // Remove identical features
                _scores.resize(_phi.size());
                _mpi_comm->barrier();
                project_funcs::project_r(_prop.data(), _scores.data(), _phi, _task_sizes, 1);
                _scores.erase(_scores.begin(), _scores.begin() + _start_gen[_start_gen.size() - 1]);
                inds = util_funcs::argsort<double>(_scores);

                std::vector<int> del_inds;

                _mpi_comm->barrier();
                node_value_arrs::clear_temp_reg();
                for(int sc = 0; sc < _scores.size() - 1; ++sc)
                {
                    #ifdef PARAMETERIZE
                    if(_phi[inds[sc] + _start_gen.back()]->n_params() > 0)
                    {
                        continue;
                    }
                    #endif

                    if(_scores[inds[sc]] > -1e-10)
                    {
                        double base_val = std::abs(
                            util_funcs::r(
                                _phi[_start_gen.back() + inds[sc]]->value_ptr(),
                                _phi[_start_gen.back() + inds[sc]]->value_ptr(),
                                _n_samp
                            )
                        );
                        for(int sc2 = sc + 1; sc2 < _scores.size(); ++sc2)
                        {
                            double comp = std::abs(
                                base_val - std::abs(
                                    util_funcs::r(
                                        _phi[_start_gen.back() + inds[sc]]->value_ptr(),
                                        _phi[_start_gen.back() + inds[sc2]]->value_ptr(0, true),
                                        _n_samp
                                    )
                                )
                            );
                            if(comp < 1e-10)
                            {
                                del_inds.push_back(-1 * (inds[sc] + _start_gen.back()));
                                break;
                            }
                        }
                    }
                    else if(_scores[inds[sc + 1]] - _scores[inds[sc]] < 1e-10)
                    {
                        double base_val = std::abs(
                            util_funcs::r(
                                _phi[_start_gen.back() + inds[sc]]->value_ptr(),
                                _phi[_start_gen.back() + inds[sc]]->value_ptr(),
                                _n_samp
                            )
                        );
                        double comp = std::abs(
                            base_val - std::abs(
                                util_funcs::r(
                                    _phi[_start_gen.back() + inds[sc]]->value_ptr(),
                                    _phi[_start_gen.back() + inds[sc + 1]]->value_ptr(0, true),
                                    _n_samp
                                )
                            )
                        );
                        if(comp < 1e-10)
                        {
                            del_inds.push_back(-1 * (inds[sc] + _start_gen.back()));
                        }
                    }
                }

                inds = util_funcs::argsort<int>(del_inds);
                for(int ii = 0; ii < inds.size(); ++ii)
                {
                    _phi.erase(_phi.begin() - del_inds[inds[ii]]);
                }

                // Reindex
                for(int ff = _start_gen.back(); ff < _phi.size(); ++ff)
                {
                    _phi[ff]->reindex(ff);
                }
            }
            node_value_arrs::clear_temp_reg();

            if(nn <= _n_rung_store)
            {
                node_value_arrs::resize_values_arr(nn, _phi.size(), (nn != _max_phi));

                for(int ff = _start_gen.back(); ff < _phi.size(); ++ff)
                {
                    _phi[ff]->set_value();
                    _phi[ff]->set_test_value();
                }
            }
        }
        else
        {
            std::vector<size_t> next_phi_sizes(_mpi_comm->size());
            if(_mpi_comm->rank() == 0)
            {
                mpi::gather(*_mpi_comm, next_phi.size(), next_phi_sizes.data(), 0);
                mpi::broadcast(*_mpi_comm, next_phi_sizes.data(), next_phi_sizes.size(), 0);
            }
            else
            {
                mpi::gather(*_mpi_comm, next_phi.size(), 0);
                mpi::broadcast(*_mpi_comm, next_phi_sizes.data(), next_phi_sizes.size(), 0);
            }

            size_t n_feat = std::accumulate(next_phi_sizes.begin(), next_phi_sizes.end(), 0);
            if(n_feat == 0)
            {
                throw std::logic_error("No features created during this rung (" + std::to_string(nn) + ")");
            }

            size_t n_feat_rank = n_feat / _mpi_comm->size();
            size_t n_feat_below_rank = _mpi_comm->rank() * n_feat_rank;
            size_t n_feat_added = 0;
            if(_mpi_comm->rank() < n_feat % _mpi_comm->size())
            {
                ++n_feat_rank;
                n_feat_below_rank += _mpi_comm->rank();
            }
            else
            {
                n_feat_below_rank += n_feat % _mpi_comm->size();
            }

            while((n_feat_added < n_feat_rank) && (next_phi.size() > 0))
            {
                _phi.push_back(next_phi.back());
                next_phi.pop_back();
                ++n_feat_added;
            }

            // This can be calculated without an all_gather, using it to not introduce too many things at one time
            std::vector<size_t> next_phi_needed(_mpi_comm->size());
            std::vector<size_t> next_phi_excess(_mpi_comm->size());
            if(_mpi_comm->rank() == 0)
            {
                mpi::gather(*_mpi_comm, next_phi.size(), next_phi_excess.data(), 0);
                mpi::gather(*_mpi_comm, n_feat_rank - n_feat_added, next_phi_needed.data(), 0);

                mpi::broadcast(*_mpi_comm, next_phi_excess.data(), next_phi_excess.size(), 0);
                mpi::broadcast(*_mpi_comm, next_phi_needed.data(), next_phi_needed.size(), 0);
            }
            else
            {
                mpi::gather(*_mpi_comm, next_phi.size(), 0);
                mpi::gather(*_mpi_comm, n_feat_rank - n_feat_added, 0);

                mpi::broadcast(*_mpi_comm, next_phi_excess.data(), next_phi_excess.size(), 0);
                mpi::broadcast(*_mpi_comm, next_phi_needed.data(), next_phi_needed.size(), 0);
            }

            std::vector<size_t> send_sizes(next_phi_sizes.size(), 0);
            std::vector<size_t> recv_sizes(next_phi_sizes.size(), 0);
            // Is this rank sending or receiving?
            if(next_phi_excess[_mpi_comm->rank()] > 0)
            {
                size_t total_sent = std::accumulate(next_phi_excess.begin(), next_phi_excess.begin() + _mpi_comm->rank(), 0);
                size_t prev_sent_recv = 0;
                size_t send_size = 0;
                int ind = 0;
                while((prev_sent_recv <= total_sent) && (ind < _mpi_comm->size()))
                {
                    prev_sent_recv += next_phi_needed[ind];
                    ++ind;
                }
                send_size = std::min(next_phi.size(), prev_sent_recv - total_sent);
                send_sizes[ind-1] = send_size;
                total_sent = send_size;
                while((total_sent < next_phi.size()) && (ind < _mpi_comm->size()))
                {
                    send_size = std::min(next_phi.size() - total_sent, next_phi_needed[ind]);
                    send_sizes[ind] = send_size;
                    total_sent += send_size;
                    ++ind;
                }

                total_sent = 0;
                for(int pp = 0; pp < send_sizes.size(); ++pp)
                {
                    if(send_sizes[pp] == 0)
                        continue;

                    std::vector<node_ptr> to_send(send_sizes[pp]);
                    std::copy_n(next_phi.begin() + total_sent, send_sizes[pp], to_send.begin());
                    _mpi_comm->send(pp, _mpi_comm->cantorTagGen(_mpi_comm->rank(), pp, 1, 0), to_send);
                    total_sent += send_sizes[pp];
                }
            }
            else
            {
                size_t total_recv = std::accumulate(next_phi_needed.begin(), next_phi_needed.begin() + _mpi_comm->rank(), 0);
                size_t prev_recv_sent = 0;
                size_t recv_size = 0;
                int ind = 0;
                while((prev_recv_sent <= total_recv) && (ind < _mpi_comm->size()))
                {
                    prev_recv_sent += next_phi_excess[ind];
                    ++ind;
                }
                recv_size = std::min(n_feat_rank - n_feat_added, prev_recv_sent - total_recv);
                recv_sizes[ind-1] = recv_size;
                total_recv = recv_size;
                while((total_recv < n_feat_rank - n_feat_added) && (ind < _mpi_comm->size()))
                {
                    recv_size = std::min(n_feat_rank - n_feat_added - total_recv, next_phi_excess[ind]);
                    recv_sizes[ind] = recv_size;
                    total_recv += recv_size;
                    ++ind;
                }

                total_recv = 0;
                for(int pp = 0; pp < recv_sizes.size(); ++pp)
                {
                    if(recv_sizes[pp] == 0)
                    {
                        continue;
                    }

                    std::vector<node_ptr> to_recv;
                    _mpi_comm->recv(pp, _mpi_comm->cantorTagGen(pp, _mpi_comm->rank(), 1, 0), to_recv);
                    for(auto& feat : to_recv)
                    {
                        _phi.push_back(feat);
                    }
                }
            }

            #pragma omp parallel for
            for(int ff = _start_gen.back(); ff < _phi.size(); ++ff)
            {
                _phi[ff]->reindex(ff + n_feat_below_rank, ff);
                _phi[ff]->set_value();
                _phi[ff]->set_test_value();
            }
        }
    }
    _n_feat = _phi.size();
}

void FeatureSpace::project_generated(const double* prop, const int size, std::vector<node_ptr>& phi_sel, std::vector<double>& scores_sel)
{
    std::vector<double> scores_sel_all(node_value_arrs::N_SELECTED);
    if(node_value_arrs::N_SELECTED > _n_sis_select)
    {
        scores_sel_all.resize(_phi_selected.size());
        _project(prop, scores_sel_all.data(), _phi_selected, _task_sizes, size / _n_samp);
    }
    std::copy_n(scores_sel.data(), scores_sel.size(), &scores_sel_all[node_value_arrs::N_SELECTED - _n_sis_select]);

    int worst_score_ind = std::max_element(scores_sel.begin(), scores_sel.end()) - scores_sel.begin();
    double worst_score = scores_sel[worst_score_ind];

    #pragma omp parallel firstprivate(worst_score, worst_score_ind, prop, size, scores_sel_all)
    {
        std::vector<node_ptr> phi_sel_private(phi_sel);
        std::vector<double> scores_sel_private(scores_sel);
        int index_base = _phi.size() + _n_sis_select * (omp_get_thread_num() + _mpi_comm->size());

        #ifdef PARAMETERIZE
        std::shared_ptr<NLOptimizer> optimizer;
        if(_reparam_residual)
        {
            std::vector<double> prop_vec(size, 0.0);
            std::copy_n(prop, size, prop_vec.data());
            optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, prop_vec, _max_phi, _max_param_depth);
        }
        else
        {
            optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, _prop, _max_phi, _max_param_depth);
        }
        #endif

        #pragma omp for schedule(dynamic)
        for(auto feat = _phi.begin() + _start_gen.back() + _mpi_comm->rank(); feat < _phi.end(); feat += _mpi_comm->size())
        {
            unsigned long int feat_ind = _phi.size() + _n_sis_select * (omp_get_num_threads() + _mpi_comm->size());

            node_value_arrs::clear_temp_reg_thread();
            std::vector<node_ptr> generated_phi;

            bool is_sel = (*feat)->selected();
            (*feat)->set_selected(false);
            #ifdef PARAMETERIZE
            generate_new_feats(feat, generated_phi, feat_ind, optimizer, _l_bound, _u_bound);
            #else
            generate_new_feats(feat, generated_phi, feat_ind, _l_bound, _u_bound);
            #endif
            (*feat)->set_selected(is_sel);

            if(generated_phi.size() == 0)
            {
                continue;
            }

            node_value_arrs::clear_temp_reg_thread();
            std::vector<double> scores(generated_phi.size());
            _project_no_omp(prop, scores.data(), generated_phi, _task_sizes, size / _n_samp);

            std::vector<int> inds = util_funcs::argsort<double>(scores);

            int ii = 0;
            while((ii < inds.size()) && (scores[inds[ii]] < -1.0))
            {
                ++ii;
            }

            while((ii < inds.size()) && ((scores[inds[ii]] < worst_score) || (phi_sel_private.size() < _n_sis_select)))
            {
                double cur_score = scores[inds[ii]];
                bool valid_feat = _is_valid(
                    generated_phi[inds[ii]]->value_ptr(0),
                    _n_samp,
                    _cross_cor_max,
                    scores_sel_all,
                    cur_score,
                    node_value_arrs::N_SELECTED - _n_sis_select,
                    0
                );
                valid_feat = valid_feat && _is_valid_feat_list(
                    generated_phi[inds[ii]]->value_ptr(0),
                    _n_samp,
                    _cross_cor_max,
                    phi_sel_private,
                    scores_sel_private,
                    cur_score
                );
                if(valid_feat)
                {
                    if(scores_sel_private.size() == _n_sis_select)
                    {
                        generated_phi[inds[ii]]->reindex(index_base + worst_score_ind);
                        generated_phi[inds[ii]]->set_value();
                        phi_sel_private[worst_score_ind] = generated_phi[inds[ii]];
                        scores_sel_private[worst_score_ind] = cur_score;
                    }
                    else
                    {
                        generated_phi[inds[ii]]->reindex(index_base + scores_sel_private.size());
                        generated_phi[inds[ii]]->set_value();
                        phi_sel_private.push_back(generated_phi[inds[ii]]);
                        scores_sel_private.push_back(cur_score);
                    }
                    worst_score_ind = std::max_element(scores_sel_private.begin(), scores_sel_private.end()) - scores_sel_private.begin();
                    worst_score = scores_sel_private[worst_score_ind];
                }
                ++ii;
            }
        }

        #pragma omp critical
        {
            index_base = _phi.size() + _n_sis_select * _mpi_comm->rank();
            worst_score_ind = std::max_element(scores_sel.begin(), scores_sel.end()) - scores_sel.begin();
            for(int sc = 0; sc < scores_sel_private.size(); ++sc)
            {
                if(
                    ((phi_sel.size() < _n_sis_select) || (scores_sel_private[sc] < scores_sel[worst_score_ind])) &&
                    _is_valid_feat_list(phi_sel_private[sc]->value_ptr(), _n_samp, _cross_cor_max, phi_sel, scores_sel, scores_sel_private[sc])
                )
                {
                    if(phi_sel.size() == _n_sis_select)
                    {
                        scores_sel[worst_score_ind] = scores_sel_private[sc];
                        phi_sel[worst_score_ind] = phi_sel_private[sc];
                        if(phi_sel[worst_score_ind]->rung() == _max_phi)
                        {
                            phi_sel[worst_score_ind]->reindex(index_base + worst_score_ind);
                        }
                    }
                    else
                    {
                        scores_sel.push_back(scores_sel_private[sc]);
                        phi_sel.push_back(phi_sel_private[sc]);
                        if(phi_sel.back()->rung() == _max_phi)
                        {
                            phi_sel.back()->reindex(index_base + phi_sel.size());
                        }
                    }
                    worst_score_ind = std::max_element(scores_sel.begin(), scores_sel.end()) - scores_sel.begin();
                }
            }
        }
    }
}

void FeatureSpace::sis(const std::vector<double>& prop)
{
    // Reparameterize for the residuals
#ifdef PARAMETERIZE
    if(_reparam_residual && (_phi_selected.size() > 0))
    {
        // Make a hard copy of the previously selected features
        for(int ff = _phi_selected.size() - _n_sis_select; ff < _phi_selected.size(); ++ff)
        {
            _phi_selected[ff] = _phi_selected[ff]->hard_copy();
        }

        // Reparameterize based on residuals
        #pragma omp parallel
        {
            std::shared_ptr<NLOptimizer> optimizer = nlopt_wrapper::get_optimizer(_project_type, _task_sizes, prop, _max_phi, _max_param_depth);
            #pragma omp for schedule(dynamic)
            for(int ff = _start_gen[1]; ff < _phi.size(); ++ff)
            {
                if(_phi[ff]->n_params() > 0)
                {
                    _phi[ff]->get_parameters(optimizer);
                }
            }
        }
    }
#endif
    // Create output directories if needed
    boost::filesystem::path p(_feature_space_file.c_str());
    boost::filesystem::create_directories(p.remove_filename());

    // Set up output file generation
    std::ofstream out_file_stream = std::ofstream();
    out_file_stream.open(_feature_space_file, std::ios::app);

    std::ofstream sum_file_stream = std::ofstream();
    sum_file_stream.open(_feature_space_summary_file, std::ios::app);

    // Set up selection data structures
    std::vector<node_ptr> phi_sel;
    std::vector<double> scores_sel(_n_sis_select, std::numeric_limits<double>::infinity());
    phi_sel.reserve(_n_sis_select);

    int cur_feat = node_value_arrs::N_SELECTED;
    node_value_arrs::resize_d_matrix_arr(_n_sis_select);
    _phi_selected.reserve(_phi_selected.size() + _n_sis_select);