search_space.cpp

Go to the documentation of this file.

#include "search_space.h"
#include "../program/program.h" // TODO: dont import this header here
 
namespace Brush{
 
 
float calc_initial_weight(const ArrayXf& value, const ArrayXf& y)
{
    // OBS: only for terminals!
 
    // weights are initialized as the slope of the z-score of x and y.
 
    // If y has different length from X, we get a core dump in this function.
    // That is why Dataset makes a check for this 
    // TODO: need to make SS (or Datasaet) check for this when loading the data
 
    vector<char> dtypes = {'f', 'f'};
 
    MatrixXf data(value.size(), 2);
 
    data.col(0) << value;
    data.col(1) << y;
 
    Normalizer n(true);
    n.fit_normalize(data, dtypes); // normalize works row-wise
 
    // In slope function, argument order matters (if not normalized with z-score)
    // The feature should be the first value, and the true value the second
    // (it will divide covar(arg1, arg2) by var(arg2)).
    // Since z-score normalizes so mean=0 and std=1, then order doesnt matter
    float prob_change = std::abs(slope(data.col(0).array() ,   // x=variable
                                       data.col(1).array() )); // y=target
 
    // having a minimum feature weight if it was not set to zero
    if (std::abs(prob_change)<1e-5)
        prob_change = 1e-5;
 
    // prob_change will evaluate to nan if variance(x)==0. Features with
    // zero variance should not be used (as they behave just like a constant).
    if (std::isnan(prob_change))
        prob_change = 0.0;
 
    return prob_change;
}
 

vector<Node> generate_terminals(const Dataset& d, const bool weights_init)
{
    vector<Node> terminals;
    int i = 0;
    for ( const auto &[feature_name, value]: d.features ) 
    {
        std::visit(
            [&](auto&& arg) {
                using T = std::decay_t<decltype(arg)>;
                using Scalar = typename T::Scalar;
                constexpr bool weighted = !std::is_same_v<Scalar, bool>;
                // if constexpr (T::Scalar == float)
                auto n = Node(
                    NodeType::Terminal, 
                    Signature<T()>{}, 
                    weighted,
                    feature_name
                );
 
                float prob_change = 1.0; // default value
                
                if (d.y.size()>0 && weights_init) 
                {
                    // if the value can be casted to float array, we can calculate slope
                    if (std::holds_alternative<ArrayXb>(value))
                    {
                        ArrayXf tmp = std::get<ArrayXb>(value).template cast<float>();
                        prob_change = calc_initial_weight(tmp, d.y);
                    }
                    else if (std::holds_alternative<ArrayXi>(value))
                    {
                        // for each variable we create a one-vs-all binary variable, then
                        // calculate slope. Final value will be the average of slopes
 
                        auto tmp = std::get<ArrayXi>(value);
 
                        //get number of unique values
                        std::map<float, bool> uniqueMap;
                        for(int i = 0; i < tmp.size(); i++)
                            uniqueMap[(float)tmp(i)] = true;
 
                        ArrayXf slopes = ArrayXf::Ones(uniqueMap.size());
                        int slopesIterator = 0;
                        for (const auto& pair : uniqueMap)
                        {
                            auto one_vs_all = ArrayXf::Ones(tmp.size()).array() * (tmp.array()==pair.first).cast<float>();
 
                            slopes[slopesIterator++] = calc_initial_weight(one_vs_all, d.y);
                        }
                        
                        prob_change = slopes.mean();
                    }
                    else if (std::holds_alternative<ArrayXf>(value)) 
                    {
                        prob_change = calc_initial_weight(std::get<ArrayXf>(value), d.y);
                    }
                    else
                    {
                        auto msg = fmt::format("Brush coudn't calculate the initial weight of variable {}\n",feature_name);
                        HANDLE_ERROR_THROW(msg);
                    }
                }
                
                n.set_prob_change( prob_change );
 
                terminals.push_back(n);
            },
            value
        );
        ++i;
    };
 
    // iterate through terminals and take the average of values of same signature
    auto signature_avg = [terminals, weights_init](DataType ret_type){
        if (!weights_init)
            return 1.0f;
 
        float sum = 0.0;
        int count = 0;
 
        for (const auto& n : terminals) {
            if (n.ret_type == ret_type && n.get_prob_change()>0) {
                sum += n.get_prob_change();
                count++;
            }
        }
 
        if (count==0) // no terminals of datatype. return 1.0 for the constant
            return 1.0f;
 
        return sum / count;
    };
    
    // constants for each type -- floats, integers, boolean. This is useful
    // to ensure the search space will always have something to fill up open spaces
    // when building/modifying trees
    auto cXf = Node(NodeType::Constant, Signature<ArrayXf()>(), true, "constF");
    cXf.set_prob_change(signature_avg(cXf.ret_type));
    terminals.push_back(cXf);
    
    // Reminder: Integers are used to represent categorical variables
    auto cXi = Node(NodeType::Constant, Signature<ArrayXi()>(), true, "constI");
    cXi.set_prob_change(signature_avg(cXi.ret_type));
    terminals.push_back(cXi);
 
    // Boolean constants are needed to prevent crashes when logical operators (And, Or, Not, Geq, Equals) 
    // appear in programs and need boolean terminals as arguments, even in regression problems
    auto cXb = Node(NodeType::Constant, Signature<ArrayXb()>(), true, "constB");
    cXb.set_prob_change(signature_avg(cXb.ret_type));
    terminals.push_back(cXb);
 
    // mean label node. Does not need to be part of symbols. works only for classification
    if (d.classification)
    {
        auto meanlabel = Node(NodeType::MeanLabel, Signature<ArrayXb()>(), true, "MeanLabel");
        meanlabel.set_prob_change(1.0f); // always present in classification problems
        terminals.push_back(meanlabel);
    }
 
    return terminals;
};
 
std::unordered_map<std::size_t, std::string> ArgsName; 
 
void SearchSpace::print() const { 
    std::cout << fmt::format("{}\n", *this) << std::flush; 
}
 
void SearchSpace::init(const Dataset& d, const unordered_map<string,float>& user_ops,
                       bool weights_init)
{
    // fmt::print("constructing search space...\n");
    this->node_map.clear();
    this->node_map_weights.clear();
    this->terminal_map.clear();
    this->terminal_types.clear();
    this->terminal_weights.clear();
    this->op_names.clear();
 
    bool use_all = user_ops.size() == 0;
    if (use_all)
    {
        for (size_t op_index=0; op_index< NodeTypes::Count; op_index++){
            op_names.push_back(    
                NodeTypeName.at(static_cast<NodeType>(1UL << op_index))
            );
        }
    }
    else
    {
        for (const auto& [op, weight] : user_ops)
            op_names.push_back(op);
    }
 
    // create nodes based on data types 
    terminal_types = d.unique_data_types;
 
    vector<Node> terminals = generate_terminals(d, weights_init);
    
    // If it is a classification problem, we need to add the fixed root nodes 
    // (logistic for binary classification, softmax for multiclassification).
    // Sometimes, the user may not specify these two nodes as candidates when 
    // sampling functions, so we check if they are already in the terminal set, and
    // we add them with zero prob if they are not. They need to be in the func set
    // when calling GenerateNodeMap, so the search_space will contain all the hashes
    // and signatures for them (and they can be used only in program root).
    // TODO: fix softmax and add it here
 
    /* fmt::print("generate nodetype\n"); */
 
    if (d.classification)
    {
        std::unordered_map<std::string, float> extended_user_ops=user_ops;
        
        if (!use_all)
        {
            // Brush can start from a population of decision trees. Thess `if`s
            // below ensures that split nodes will always be in the search space
            if (extended_user_ops.find("SplitOn") == extended_user_ops.end()){
                extended_user_ops.insert({"SplitOn", 0.0f});
                op_names.push_back("SplitOn");
            }
            if (extended_user_ops.find("SplitBest") == extended_user_ops.end()){
                extended_user_ops.insert({"SplitBest", 0.0f});
                op_names.push_back("SplitBest");
            }
 
            // We need some ops in the search space so we can have the logit and offset
            if (extended_user_ops.find("OffsetSum") == extended_user_ops.end()){
                extended_user_ops.insert({"OffsetSum", 0.0f});
                op_names.push_back("OffsetSum");
            }
 
            // Convert ArrayXf to std::vector<float> for compatibility with std::set
            std::vector<float> vec(d.y.data(), d.y.data() + d.y.size());
            std::set<float> unique_classes(vec.begin(), vec.end());
 
            if (unique_classes.size()==2 && (extended_user_ops.find("Logistic") == extended_user_ops.end())) {
                extended_user_ops.insert({"Logistic", 0.0f});
                op_names.push_back("Logistic");
            }
            else if (extended_user_ops.find("Softmax") == extended_user_ops.end()) {
                extended_user_ops.insert({"Softmax", 0.0f});
                op_names.push_back("Softmax");
            }
        }
 
        // fmt::print("generate nodetype\n");
        GenerateNodeMap(extended_user_ops, terminal_types, 
                        std::make_index_sequence<NodeTypes::OpCount>());
    }
    else
    {
        GenerateNodeMap(user_ops, terminal_types, 
                    std::make_index_sequence<NodeTypes::OpCount>());
    }
 
    // map terminals
    /* fmt::print("looping through terminals...\n"); */
    for (const auto& term : terminals)
    {
        /* fmt::print("adding {} to search space...\n", term.get_name()); */
        if (terminal_map.find(term.ret_type) == terminal_map.end())
            terminal_map[term.ret_type] = vector<Node>();
 
        /* fmt::print("terminal ret_type: {}\n", DataTypeName[term.ret_type]); */
        terminal_map[term.ret_type].push_back(term);
        terminal_weights[term.ret_type].push_back(term.get_prob_change());
    }
};
 
std::optional<tree<Node>> SearchSpace::sample_subtree(Node root, int max_d, int max_size) const 
{
    // public interface to use PTC2 algorithm
 
    // PTC is designed to not fail (it will persistently try to find nodes with
    // sampling functions). In pop initialization, this shoudnt be a problem, but
    // during evolution, due to dynamic changes in node weights by the learners, 
    // it now may fail. We need to check, before calling it, that it has elements
    // in search space to sample
    auto ret_match = node_map.at(root.ret_type);
 
    vector<float> args_w = get_weights(root.ret_type);
 
    // at least one operator that matches the weight must have positive probability
    if (!has_solution_space(args_w.begin(), args_w.end()))
        return std::nullopt;
 
    // it will always have a terminal (because we create constants).
    // TODO: I guess I can remove this line below and it will still work
    // if ( (terminal_map.find(root.ret_type) == terminal_map.end())
    // ||   (!has_solution_space(terminal_weights.at(root.ret_type).begin(), 
    //                           terminal_weights.at(root.ret_type).end())) )
    //     return std::nullopt;
 
    auto Tree = tree<Node>();
    auto spot = Tree.insert(Tree.begin(), root);
 
    // we should notice the difference between size of a PROGRAM and a TREE.
    // program count weights in its size, while the TREE structure dont. Wenever
    // using size of a program/tree, make sure you use the function from the correct class
    PTC2(Tree, spot, max_d, max_size, false);
    
    return Tree;
};
 
tree<Node>& SearchSpace::PTC2(tree<Node>& Tree,
    tree<Node>::iterator spot, int max_d, int max_size, bool start_from_decision_trees) const
{
    // PTC2 is agnostic of program type
 
    // A comment about PTC2 method:            
    // PTC2 can work with depth or size restriction, but it does not strictly
    // satisfies these conditions all time. Given a `max_size` and `max_depth`
    // parameters, the real maximum size that can occur is `max_size` plus the
    // highest operator arity, and the real maximum depth is `max_depth` plus one.
 
    // Queue of nodes that need children
    vector<tuple<TreeIter, DataType, int>> queue; 
 
    // node depth
    int d = 1;
    // current tree size
    int s = 1;
 
    Node root = spot.node->data;
 
    // updating size accordingly to root node
    if (Is<NodeType::SplitBest>(root.node_type) || Is<NodeType::SplitOn>(root.node_type)){
        s += 2;
    }        
    else if ( root.get_is_weighted()==true
              && Isnt<NodeType::Constant, NodeType::MeanLabel>(root.node_type) )
        s += 2;
        
    //For each argument position a of n, Enqueue(a; g) 
    for (auto a : root.arg_types)
    { 
        auto child_spot = Tree.append_child(spot);
        queue.push_back(make_tuple(child_spot, a, d));
    }
 
    int max_arity = 4;
 
    Node n;
 
    // Now we actually start the PTC2 procedure to create the program tree.
    // we are considering that all nodes in the queue will be weighted, so we
    // multiply the queue size by 3
    while ( 3*queue.size() + s < max_size && queue.size() > 0) 
    {            
        // including the queue size in the max_size, since each element in queue
        // can grow up exponentially
 
        // by default, terminals are weighted (counts as 3 nodes in program size).
        // since every spot in queue has potential to be a terminal, we multiply
        // its size by 3. Subtracting one due to the fact that this loop will
        // always insert a non terminal (which by default has weights off).
        // this way, we can have PTC2 working properly.
        
        auto [qspot, t, d] = RandomDequeue(queue);
 
        if (d >= max_d || s >= max_size)
        {
            auto opt = sample_terminal(t);
 
            // if it returned optional, then there's nothing to sample based on weights.
            // We'll force sampling again with uniform probs
            if (!opt)
                opt = sample_terminal(t, true);
 
            assert (opt && "PTC2: failed to sample terminal to fill in empty leaves due to max depth or size");
 
            // If we successfully get a terminal, use it
            n = opt.value();
 
            Tree.replace(qspot, n);
        }
        else
        {
            //choose a nonterminal of matching type
            std::optional<Node> opt;
    
            if (start_from_decision_trees) { // 
                // this sample_op is likely to never fail. Only case scenario is
                // a logistic root (arrayF arg type) with no arrayF arguments. 
                // (in this case, sample_terminal will also fail)
                opt = sample_op(NodeType::SplitBest, t, true);
            }
            else {
                opt = sample_op(t);
            }
 
            if (!opt) { // there is no operator for this node. sample a terminal instead
                opt = sample_terminal(t);
 
                // didnt work the easy way, lets try the hard way
                if (!opt)
                    opt = sample_terminal(t, true);
            }
 
            if (!opt) { // no operator or terminal - search space is ill defined
                auto msg = fmt::format(
                    "Failed to sample operator AND terminal of data type {} during PTC2.\n"
                    "start_from_decision_trees: {}\n",
                    DataTypeName[t],
                    start_from_decision_trees ? "true" : "false"
                );
 
                HANDLE_ERROR_THROW(msg);
 
                // queue.push_back(make_tuple(qspot, t, d));
                // continue;
            }
 
            n = opt.value();
            
            auto newspot = Tree.replace(qspot, n);
 
            // For each arg of n, add to queue
            for (auto a : n.arg_types)
            {
                auto child_spot = Tree.append_child(newspot);
 
                queue.push_back(make_tuple(child_spot, a, d+1));
            }
        }
 
        // increment is different based on node weights
        ++s;
        
        if (Is<NodeType::SplitBest>(n.node_type) || Is<NodeType::SplitOn>(n.node_type)){
            // cout << "sampled split on node/n";
            s += 2;
        }
        else if ( n.get_is_weighted()==true
                  && Isnt<NodeType::Constant, NodeType::MeanLabel>(n.node_type) )
            s += 2;
    } 
 
    while (queue.size() > 0)
    {
        if (queue.size() == 0)
            break;
 
        auto [qspot, t, d] = RandomDequeue(queue);
 
        auto opt = sample_terminal(t);
 
        if (!opt)
            opt = sample_terminal(t, true);
            
        assert (opt && "PTC2: failed to sample terminal to fill remaining spots in queue");
 
        n = opt.value();
        
        auto newspot = Tree.replace(qspot, n);
    }
    return Tree;
};
 
// TODO: stop using params as a default argument and actually pass it (also update tests)
RegressorProgram SearchSpace::make_regressor(int max_d, int max_size, const Parameters& params)
{
    return make_program<RegressorProgram>(params, max_d, max_size);
};
 
ClassifierProgram SearchSpace::make_classifier(int max_d, int max_size, const Parameters& params)
{
    return make_program<ClassifierProgram>(params, max_d, max_size);
};
 
MulticlassClassifierProgram SearchSpace::make_multiclass_classifier(
    int max_d, int max_size, const Parameters& params)
{
    return make_program<MulticlassClassifierProgram>(params, max_d, max_size);
};
 
RepresenterProgram SearchSpace::make_representer(int max_d, int max_size, const Parameters& params)
{
    return make_program<RepresenterProgram>(params, max_d, max_size);
};
 
} //Brush