refactor

Author: mathias-von-ottenbreit

cpp/APLRRegressor.h

@@ -483,7 +483,8 @@ void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors
     intercept_steps=VectorXd::Constant(m,0);
 
     terms_eligible_current.reserve(X_train.cols()*reserved_terms_times_num_x);
-    for (size_t i = 0; i < static_cast<size_t>(X_train.cols()); ++i)
+    size_t X_train_cols{static_cast<size_t>(X_train.cols())};
+    for (size_t i = 0; i < X_train_cols; ++i)
     {
         bool term_has_one_unique_value{check_if_base_term_has_only_one_unique_value(i)};
         Term copy_of_base_term{Term(i)};
@@ -495,7 +496,7 @@ void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors
     }
 
     predictor_indexes.resize(X_train.cols());
-    for (size_t i = 0; i < static_cast<size_t>(X_train.cols()); ++i)
+    for (size_t i = 0; i < X_train_cols; ++i)
     {
         predictor_indexes[i]=i;
     }
@@ -857,7 +858,7 @@ void APLRRegressor::add_necessary_given_terms_to_interaction(Term &interaction,
             }
 
             bool given_term_provides_an_unique_zero{false};
-            for (size_t row = 0; row < static_cast<size_t>(X_train.rows()); ++row)
+            for (Eigen::Index row = 0; row < X_train.rows(); ++row)
             {
                 given_term_provides_an_unique_zero = combined_value_indicator_for_the_other_given_terms[row]>0 && value_indicator_for_each_given_term.col(col)[row]==0;
                 if(given_term_provides_an_unique_zero)
@@ -892,7 +893,7 @@ void APLRRegressor::add_promising_interactions_and_select_the_best_one()
     size_t best_term_before_interactions{best_term_index};
     bool best_term_before_interactions_was_not_selected{best_term_before_interactions==std::numeric_limits<size_t>::max()};
     bool error_is_less_than_for_best_term_before_interactions;
-    for (size_t j = 0; j < static_cast<size_t>(sorted_indexes_of_errors_for_interactions_to_consider.size()); ++j) //for each interaction to consider starting from lowest to highest error
+    for (Eigen::Index j = 0; j < sorted_indexes_of_errors_for_interactions_to_consider.size(); ++j) //for each interaction to consider starting from lowest to highest error
     {
         bool allowed_to_add_one_interaction{interactions_eligible<max_interactions};
         if(allowed_to_add_one_interaction)
@@ -1121,7 +1122,7 @@ void APLRRegressor::revert_scaling_if_using_log_link_function()
     {
         y_train/=scaling_factor_for_log_link_function;
         intercept+=std::log(1/scaling_factor_for_log_link_function);
-        for (size_t i = 0; i < static_cast<size_t>(intercept_steps.size()); ++i)
+        for (Eigen::Index i = 0; i < intercept_steps.size(); ++i)
         {
             intercept_steps[i]+=std::log(1/scaling_factor_for_log_link_function);
         }
@@ -1134,7 +1135,8 @@ void APLRRegressor::name_terms(const MatrixXd &X, const std::vector<std::string>
     if(x_names_not_provided)
     {
         std::vector<std::string> temp(X.cols());
-        for (size_t i = 0; i < static_cast<size_t>(X.cols()); ++i)
+        size_t X_cols{static_cast<size_t>(X.cols())};
+        for (size_t i = 0; i < X_cols; ++i)
         {
             temp[i]="X"+std::to_string(i+1);
         }
@@ -1207,7 +1209,7 @@ void APLRRegressor::calculate_feature_importance_on_validation_set()
 {
     feature_importance=VectorXd::Constant(number_of_base_terms,0);
     MatrixXd li{calculate_local_feature_importance(X_validation)};
-    for (size_t i = 0; i < static_cast<size_t>(li.cols()); ++i) //for each column calculate mean abs values
+    for (Eigen::Index i = 0; i < li.cols(); ++i) //for each column calculate mean abs values
     {
         feature_importance[i]=li.col(i).cwiseAbs().mean();
     }
@@ -1302,7 +1304,7 @@ VectorXd APLRRegressor::calculate_linear_predictor(const MatrixXd &X)
 
 void APLRRegressor::cap_predictions_to_minmax_in_training(VectorXd &predictions)
 {
-    for (size_t i = 0; i < static_cast<size_t>(predictions.rows()); ++i)
+    for (Eigen::Index i = 0; i < predictions.rows(); ++i)
     {
         if(std::isgreater(predictions[i],max_training_prediction_or_response))
             predictions[i]=max_training_prediction_or_response;

cpp/functions.h

@@ -150,7 +150,7 @@ VectorXd calculate_exp_of_linear_predictor_adjusted_for_numerical_problems(const
     VectorXd exp_of_linear_predictor{linear_predictor.array().exp()};
     double min_exp_of_linear_predictor{std::exp(min_exponent)};
     double max_exp_of_linear_predictor{std::exp(max_exponent)};
-    for (size_t i = 0; i < static_cast<size_t>(linear_predictor.rows()); ++i)
+    for (Eigen::Index i = 0; i < linear_predictor.rows(); ++i)
     {            
         bool linear_predictor_is_too_small{std::isless(linear_predictor[i], min_exponent)};
         if(linear_predictor_is_too_small)
@@ -270,7 +270,7 @@ void throw_error_if_matrix_has_nan_or_infinite_elements(const T &x, const std::s
 VectorXi calculate_indicator(const VectorXd &v)
 {
     VectorXi indicator{VectorXi::Constant(v.rows(),1)};
-    for (size_t i = 0; i < static_cast<size_t>(v.size()); ++i)
+    for (Eigen::Index i = 0; i < v.size(); ++i)
     {
         if(is_approximately_zero(v[i]))
             indicator[i]=0;
@@ -281,7 +281,7 @@ VectorXi calculate_indicator(const VectorXd &v)
 VectorXi calculate_indicator(const VectorXi &v)
 {
     VectorXi indicator{VectorXi::Constant(v.rows(),1)};
-    for (size_t i = 0; i < static_cast<size_t>(v.size()); ++i)
+    for (Eigen::Index i = 0; i < v.size(); ++i)
     {
         if(v[i]==0)
             indicator[i]=0;
@@ -369,7 +369,7 @@ double trapezoidal_integration(const VectorXd &y, const VectorXd &x)
     if(y_is_large_enough && x_and_y_have_the_same_size)
     {
         output=0;
-        for (size_t i = 1; i < static_cast<size_t>(y.size()); ++i)
+        for (Eigen::Index i = 1; i < y.size(); ++i)
         {
             double delta_y{(y[i]+y[i-1])/2};
             double delta_x{x[i]-x[i-1]};

cpp/term.h

@@ -36,8 +36,6 @@ class Term
     double error_where_given_terms_are_zero;
     SortedData sorted_vectors;
     VectorXd negative_gradient_discretized;
-    VectorXd errors_initial;
-    double error_initial;
     std::vector<size_t> observations_in_bins;
     int monotonic_constraint;
 
@@ -49,8 +47,8 @@ class Term
     void setup_bins();
     void discretize_data_by_bin();
     void estimate_split_point_on_discretized_data();
-    void calculate_coefficient_and_error_on_discretized_data(bool direction_right, double split_point);
-    void estimate_coefficient_and_error_on_all_data();
+    void estimate_coefficient_and_error(const VectorXd &x, const VectorXd &y, const VectorXd &sample_weight,double error_added=0.0);
+    double estimate_coefficient(const VectorXd &x, const VectorXd &y, const VectorXd &sample_weight=VectorXd(0));
     void cleanup_after_estimate_split_point();
     void cleanup_after_fit();
     void cleanup_when_this_term_was_added_as_a_given_term();
@@ -174,7 +172,8 @@ void Term::estimate_split_point(const MatrixXd &X,const VectorXd &negative_gradi
     }
     discretize_data_by_bin();
     estimate_split_point_on_discretized_data();
-    estimate_coefficient_and_error_on_all_data();
+    estimate_coefficient_and_error(calculate_without_interactions(sorted_vectors.values_sorted),sorted_vectors.negative_gradient_sorted,
+        sorted_vectors.sample_weight_sorted,error_where_given_terms_are_zero);
     cleanup_after_estimate_split_point();
     determine_if_can_be_used_as_a_given_term(X.col(base_term));
 }
@@ -188,7 +187,7 @@ void Term::calculate_rows_to_zero_out_and_not_due_to_given_terms(const MatrixXd
         for (auto &given_term:given_terms)
         {
             VectorXd values_given_term{given_term.calculate(X)};
-            for (size_t i = 0; i < static_cast<size_t>(X.rows()); ++i)
+            for (Eigen::Index i = 0; i < X.rows(); ++i)
             {
                 if(is_approximately_zero(values_given_term[i]))
                 {
@@ -200,7 +199,7 @@ void Term::calculate_rows_to_zero_out_and_not_due_to_given_terms(const MatrixXd
         rows_to_zero_out_and_not_due_to_given_terms.zeroed.resize(X.rows()-rows_to_zero_out_and_not_due_to_given_terms.not_zeroed.rows());
         size_t count_zeroed{0};
         size_t count_not_zeroed{0};
-        for (size_t i = 0; i < static_cast<size_t>(X.rows()); ++i)
+        for (Eigen::Index i = 0; i < X.rows(); ++i)
         {
             bool value_is_non_zero{non_zero_values[i]==1};
             if(value_is_non_zero)
@@ -233,7 +232,7 @@ VectorXd Term::calculate(const MatrixXd &X)
         for(auto &given_term:given_terms)
         {
             VectorXd values_given_term{given_term.calculate(X)};
-            for (size_t i = 0; i < static_cast<size_t>(values.size()); ++i)
+            for (Eigen::Index i = 0; i < values.size(); ++i)
             {
                 if(is_approximately_zero(values_given_term[i]))
                     values[i]=0;
@@ -276,14 +275,14 @@ void Term::calculate_error_where_given_terms_are_zero(const VectorXd &negative_g
     {
         if(sample_weight.size()==0)
         {
-            for (size_t i = 0; i < static_cast<size_t>(rows_to_zero_out_and_not_due_to_given_terms.zeroed.size()); ++i)
+            for (Eigen::Index i = 0; i < rows_to_zero_out_and_not_due_to_given_terms.zeroed.size(); ++i)
             {
                 error_where_given_terms_are_zero+=calculate_error_one_observation(negative_gradient[rows_to_zero_out_and_not_due_to_given_terms.zeroed[i]],0.0,NAN_DOUBLE);
             }
         }
         else
         {
-            for (size_t i = 0; i < static_cast<size_t>(rows_to_zero_out_and_not_due_to_given_terms.zeroed.size()); ++i)
+            for (Eigen::Index i = 0; i < rows_to_zero_out_and_not_due_to_given_terms.zeroed.size(); ++i)
             {
                 error_where_given_terms_are_zero+=calculate_error_one_observation(negative_gradient[rows_to_zero_out_and_not_due_to_given_terms.zeroed[i]],0.0,sample_weight[rows_to_zero_out_and_not_due_to_given_terms.zeroed[i]]);
             }
@@ -504,40 +503,32 @@ void Term::discretize_data_by_bin()
 }
 
 void Term::estimate_split_point_on_discretized_data()
-{
-    errors_initial=calculate_errors(negative_gradient_discretized,VectorXd::Constant(negative_gradient_discretized.size(),0.0),
-        sample_weight_discretized,FAMILY_GAUSSIAN);
-    error_initial=calculate_sum_error(errors_initial);
-
-    double split_point_temp;
-    
-    bool SPLIT_POINT_NAN{false};
-    calculate_coefficient_and_error_on_discretized_data(SPLIT_POINT_NAN, NAN_DOUBLE);
-    double error_cp_nan{split_point_search_errors_sum};
+{    
+    split_point=NAN_DOUBLE;
+    estimate_coefficient_and_error(calculate_without_interactions(values_discretized),negative_gradient_discretized,sample_weight_discretized);
+    double error_split_point_nan{split_point_search_errors_sum};
 
-    bool DIRECTION_LEFT{false};
     double split_point_left{NAN_DOUBLE};
-    double error_min_left{error_cp_nan};
-    for (size_t i = 0; i < bins_split_points_left.size(); ++i)
+    double error_min_left{error_split_point_nan};
+    for(auto &bin:bins_split_points_left)
     {
-        split_point_temp=bins_split_points_left[i];
-     
-        calculate_coefficient_and_error_on_discretized_data(DIRECTION_LEFT, split_point_temp);
+        split_point=bin;
+        direction_right=false;
+        estimate_coefficient_and_error(calculate_without_interactions(values_discretized),negative_gradient_discretized,sample_weight_discretized);
         if(std::islessequal(split_point_search_errors_sum,error_min_left))
         {
             error_min_left=split_point_search_errors_sum;
             split_point_left=split_point;
         }
     }
 
-    bool DIRECTION_RIGHT{true};
     double split_point_right{NAN_DOUBLE};
-    double error_min_right{error_cp_nan};
-    for (size_t i = 0; i < bins_split_points_right.size(); ++i)
+    double error_min_right{error_split_point_nan};
+    for(auto &bin:bins_split_points_right)
     {
-        split_point_temp=bins_split_points_right[i];
-     
-        calculate_coefficient_and_error_on_discretized_data(DIRECTION_RIGHT, split_point_temp);
+        split_point=bin;
+        direction_right=true;
+        estimate_coefficient_and_error(calculate_without_interactions(values_discretized),negative_gradient_discretized,sample_weight_discretized);
         if(std::islessequal(split_point_search_errors_sum,error_min_right))
         {
             error_min_right=split_point_search_errors_sum;
@@ -560,52 +551,52 @@ void Term::estimate_split_point_on_discretized_data()
     }
 }
 
-void Term::calculate_coefficient_and_error_on_discretized_data(bool direction_right, double split_point)
+void Term::estimate_coefficient_and_error(const VectorXd &x, const VectorXd &y, const VectorXd &sample_weight, double error_added)
 {
-    this->direction_right=direction_right;
-    this->split_point=split_point;
-    
-    VectorXd values_sorted{calculate_without_interactions(values_discretized)};
-
-    size_t index_start{0};
-    size_t index_end{max_index_discretized};
-
-    double xwx{0};
-    double xwy{0};
-    for (size_t i = index_start; i <= index_end; ++i)
+    coefficient = estimate_coefficient(x,y,sample_weight);
+    if(std::isfinite(coefficient))
     {
-        xwx+=values_sorted[i]*values_sorted[i]*sample_weight_discretized[i];
-        xwy+=values_sorted[i]*negative_gradient_discretized[i]*sample_weight_discretized[i];
-    }
-    if(xwx!=0)
-    {
-        double error_reduction{0};
-        coefficient=xwy/xwx*v;
+        coefficient*=v;
         if(coefficient_adheres_to_monotonic_constraint())
         {
-            double predicted;
-            double sample_weight_one_obs{NAN_DOUBLE};
-            for (size_t i = index_start; i <= index_end; ++i)
-            {
-                predicted=values_sorted[i]*coefficient;
-                if(sample_weight_discretized.size()>0)
-                    sample_weight_one_obs=sample_weight_discretized[i];
-
-                error_reduction+=errors_initial[i]-calculate_error_one_observation(negative_gradient_discretized[i],predicted,sample_weight_one_obs);
-            }
-            split_point_search_errors_sum=error_initial-error_reduction;
+            VectorXd predictions{x*coefficient};
+            split_point_search_errors_sum=calculate_sum_error(calculate_errors(y,predictions,sample_weight,FAMILY_GAUSSIAN))+error_added;
         }
         else
         {
             coefficient=0;
-            split_point_search_errors_sum=error_initial;
+            split_point_search_errors_sum=std::numeric_limits<double>::infinity();
         }
     }
     else
     {
         coefficient=0;
-        split_point_search_errors_sum=error_initial;
+        split_point_search_errors_sum=std::numeric_limits<double>::infinity();
+    }    
+}
+
+double Term::estimate_coefficient(const VectorXd &x, const VectorXd &y, const VectorXd &sample_weight)
+{
+    double numerator{0};
+    double denominator{0};
+    bool sample_weight_is_provided{sample_weight.size()>0};
+    if(sample_weight_is_provided)
+    {
+        for (Eigen::Index i = 0; i < y.size(); ++i)
+        {
+            numerator+=x[i]*y[i]*sample_weight[i];
+            denominator+=x[i]*x[i]*sample_weight[i];
+        }
+    }
+    else
+    {
+        for (Eigen::Index i = 0; i < y.size(); ++i)
+        {
+            numerator+=x[i]*y[i];
+            denominator+=x[i]*x[i];
+        }
     }
+    return numerator/denominator;
 }
 
 bool Term::coefficient_adheres_to_monotonic_constraint()
@@ -620,43 +611,6 @@ bool Term::coefficient_adheres_to_monotonic_constraint()
     return coefficient_adheres;
 }
 
-void Term::estimate_coefficient_and_error_on_all_data()
-{
-    sorted_vectors.values_sorted=calculate_without_interactions(sorted_vectors.values_sorted);
-    double xwx{0};
-    double xwy{0};
-    if(sorted_vectors.sample_weight_sorted.size()>0)
-    {
-        xwx=(sorted_vectors.values_sorted.array()*sorted_vectors.values_sorted.array()*sorted_vectors.sample_weight_sorted.array()).sum();
-        xwy=(sorted_vectors.values_sorted.array()*sorted_vectors.negative_gradient_sorted.array()*sorted_vectors.sample_weight_sorted.array()).sum();
-    }
-    else
-    {
-        xwx=(sorted_vectors.values_sorted.array()*sorted_vectors.values_sorted.array()).sum();
-        xwy=(sorted_vectors.values_sorted.array()*sorted_vectors.negative_gradient_sorted.array()).sum();
-    }
-    if(xwx!=0)
-    {
-        coefficient=xwy/xwx*v;
-        if(coefficient_adheres_to_monotonic_constraint())
-        {
-            VectorXd predictions{sorted_vectors.values_sorted*coefficient};
-            split_point_search_errors_sum=calculate_sum_error(calculate_errors(sorted_vectors.negative_gradient_sorted,predictions,
-                sorted_vectors.sample_weight_sorted,FAMILY_GAUSSIAN))+error_where_given_terms_are_zero;
-        }
-        else
-        {
-            coefficient=0;
-            split_point_search_errors_sum=std::numeric_limits<double>::infinity();
-        }
-    }
-    else
-    {
-        coefficient=0;
-        split_point_search_errors_sum=std::numeric_limits<double>::infinity();
-    }
-}
-
 void Term::cleanup_after_estimate_split_point()
 {
     rows_to_zero_out_and_not_due_to_given_terms.not_zeroed.resize(0);
@@ -665,7 +619,6 @@ void Term::cleanup_after_estimate_split_point()
     sorted_vectors.negative_gradient_sorted.resize(0);
     sorted_vectors.sample_weight_sorted.resize(0);
     negative_gradient_discretized.resize(0);
-    errors_initial.resize(0);
 }
 
 void Term::determine_if_can_be_used_as_a_given_term(const VectorXd &x)

setup.py

@@ -15,7 +15,7 @@
 
 setuptools.setup(
     name='aplr',
-    version='2.1.0',
+    version='2.2.0',
     description='Automatic Piecewise Linear Regression',
     ext_modules=[sfc_module],
     author="Mathias von Ottenbreit",