From b548cc60426a7acda4b2cfa40d5715bcb41452e0 Mon Sep 17 00:00:00 2001 From: Guilherme Seidyo Imai Aldeia Date: Tue, 22 Sep 2026 16:29:08 -0400 Subject: [PATCH 1/3] precision, recall, AUPRC, AUROC --- pybrush/BrushEstimator.py | 40 +++++- pybrush/EstimatorInterface.py | 21 ++- src/eval/metrics.cpp | 260 +++++++++++++++++++++++++++++++--- src/eval/metrics.h | 87 +++++++++++- src/eval/scorer.h | 7 + src/ind/individual.h | 5 +- src/selection/lexicase.cpp | 11 +- 7 files changed, 397 insertions(+), 34 deletions(-) diff --git a/pybrush/BrushEstimator.py b/pybrush/BrushEstimator.py index b3c93a1e6..1521b68fb 100644 --- a/pybrush/BrushEstimator.py +++ b/pybrush/BrushEstimator.py @@ -16,6 +16,7 @@ from sklearn.metrics import accuracy_score, balanced_accuracy_score, log_loss from sklearn.metrics import average_precision_score, mean_squared_error +from sklearn.metrics import precision_score, recall_score, roc_auc_score from pybrush import Parameters, Dataset, SearchSpace, brush_rng, individual from pybrush._brush import set_random_state as set_brush_random_state @@ -24,6 +25,28 @@ from pandas.api.types import is_float_dtype, is_bool_dtype, is_integer_dtype +def _ovr_macro(binary_metric): + """Macro one-vs-rest average of a binary ranking metric, skipping classes + that are absent (or the only class present) in `y_true`. Matches brush's + multiclass `roc_auc` and `average_precision_score`.""" + + # Auxiliary function to use in final model selection. + # The goal is to have in python an equivalent implementation of the metrics + # in c++, so we can do any kind of final model selection after running the + # evolution. + + def metric(y_true, y_score, sample_weight=None): + y_true, y_score = np.asarray(y_true), np.asarray(y_score) + scores = [] + for label in range(y_score.shape[1]): + y_bin = (y_true == label).astype(float) + if 0 < y_bin.sum() < len(y_bin): + scores.append(binary_metric(y_bin, y_score[:, label], + sample_weight=sample_weight)) + return np.mean(scores) if scores else 0.5 + return metric + + class BrushEstimator(EstimatorInterface, BaseEstimator): """ This is the base class for Brush estimators using the c++ engine. @@ -206,7 +229,8 @@ def partial_fit(self, X, y, *, # The logistic root is not affected by locking or unlocking. # It is fixed due to prob_change==0.0. - # This updates the parameters (such as class weights) + # This updates the parameters (such as class weights and the scorer) + self.parameters_ = new_parameters self.engine_.params = new_parameters # replicating the best individual @@ -291,15 +315,22 @@ def _update_final_model(self, data=None): "multi_log": log_loss, "accuracy": accuracy_score, "balanced_accuracy": balanced_accuracy_score, - "average_precision_score": average_precision_score + "average_precision_score": average_precision_score, + "precision": precision_score, + "recall": recall_score, + "roc_auc": roc_auc_score, } loss_f = loss_f_dict[self.parameters_.scorer] + multiclass = self.mode == 'classification' and self.parameters_.n_classes > 2 + if multiclass and self.parameters_.scorer in ["roc_auc", "average_precision_score"]: + loss_f = _ovr_macro(loss_f) + def eval(ind, sample=None): if sample is None: sample = np.arange(len(y)) - if self.parameters_.scorer in ["log", "multi_log", "average_precision_score"]: + if self.parameters_.scorer in ["log", "multi_log", "average_precision_score", "roc_auc"]: y_pred = np.array(ind.predict_proba(data)) else: # accuracy, balanced accuracy, or regression metrics y_pred = np.array(ind.predict(data)) @@ -307,6 +338,9 @@ def eval(ind, sample=None): metric_kwargs = {} if self.parameters_.scorer == "multi_log": metric_kwargs["labels"] = np.arange(self.parameters_.n_classes) + elif self.parameters_.scorer in ["precision", "recall"]: + metric_kwargs["zero_division"] = 0 + metric_kwargs["average"] = "macro" if multiclass else "binary" # y_pred = np.nan_to_num(y_pred) # Protecting the evaluation diff --git a/pybrush/EstimatorInterface.py b/pybrush/EstimatorInterface.py index d5e67ebe2..675584f0a 100644 --- a/pybrush/EstimatorInterface.py +++ b/pybrush/EstimatorInterface.py @@ -76,8 +76,18 @@ class EstimatorInterface(): used as secondary objective. scorer : str, default None The metric to use for the "scorer" objective. If None, it will be set to - "mse" for regression and "log" for binary classification. - Available options are `["mse", "log", "multi_log", "accuracy", "balanced_accuracy", "average_precision_score"]` + "mse" for regression, "log" for binary classification, and "multi_log" + for multiclass classification. + Available options are `["mse"]` for regression, and + `["log", "accuracy", "balanced_accuracy", "average_precision_score", + "precision", "recall", "roc_auc"]` for binary classification (with + `"multi_log"` replacing `"log"` for multiclass classification). For + multiclass problems, `"precision"` and `"recall"` are macro averages, and + `"roc_auc"` and `"average_precision_score"` are macro one-vs-rest averages. + The scorer drives selection, survival, and the archive, but it is not + used to fit parameters: weights are always optimized with the log loss, + and split thresholds with the gini impurity. You can change the scorer + between calls to `partial_fit`. algorithm : {"nsga2island", "nsga2", "gaisland", "ga"}, default "nsga2" Which Evolutionary Algorithm framework to use to evolve the population. This is used only in DeapEstimators. @@ -349,10 +359,13 @@ def _wrap_parameters(self, y, **extra_kwargs): "Invalid scorer for the regression mode" elif params.n_classes == 2: assert self.scorer in ['log', 'balanced_accuracy', 'accuracy', - 'average_precision_score'], \ + 'average_precision_score', 'precision', + 'recall', 'roc_auc'], \ "Invalid scorer for binary classification" else: - assert self.scorer in ['multi_log', 'balanced_accuracy', 'accuracy'], \ + assert self.scorer in ['multi_log', 'balanced_accuracy', 'accuracy', + 'average_precision_score', 'precision', + 'recall', 'roc_auc'], \ "Invalid scorer for multiclass classification" params.scorer = self.scorer diff --git a/src/eval/metrics.cpp b/src/eval/metrics.cpp index 79cb4c630..7864c2d99 100644 --- a/src/eval/metrics.cpp +++ b/src/eval/metrics.cpp @@ -131,29 +131,25 @@ float bal_zero_one_loss(const VectorXf& y, return (TPR + TNR) / 2.0; } -float average_precision_score(const VectorXf& y, const VectorXf& predict_proba, - VectorXf& loss, - const vector& class_weights) { - - // AP is implemented as AUC PR in sklearn. - // AP summarizes a precision-recall curve as the weighted mean of precisions - // achieved at each threshold, with the increase in recall from the previous threshold used as the weight +// anonymous namespace. make the headers private to metrics.cpp. only affects linkage +namespace { - // Assuming y contains binary labels (0 or 1) - int num_instances = y.size(); - - float eps = 1e-6f; // first we set the loss vector values - loss.resize(num_instances); - for (int i = 0; i < num_instances; ++i) { - float p = predict_proba(i); +// per-sample weights from class weights (all ones if no class weights) +vector sample_weights(const VectorXf& y, const vector& class_weights) +{ + vector w(y.size(), 1.0f); + if (!class_weights.empty()) + for (int i = 0; i < y.size(); ++i) + w[i] = class_weights.at(static_cast(y(i))); + return w; +} - // The loss vector is used in lexicase selection. we need to set something useful here - // that does make sense on individual level. Using log loss here. - if (p < eps || 1 - p < eps) - loss(i) = -(y(i)*log(eps) + (1-y(i))*log(1-eps)); - else - loss(i) = -(y(i)*log(p) + (1-y(i))*log(1-p)); - } +// Binary average precision. `y` holds 0/1 labels and `w` per-sample weights. +float binary_average_precision(const VectorXf& y, const VectorXf& predict_proba, + const vector& w) +{ + int num_instances = y.size(); + float eps = 1e-6f; // get argsort of predict proba (descending) vector order(num_instances); @@ -171,7 +167,7 @@ float average_precision_score(const VectorXf& y, const VectorXf& predict_proba, y_sorted[i] = y(idx); p_sorted[i] = predict_proba(idx); - w_sorted[i] = class_weights.empty() ? 1.0f : class_weights.at(y(idx)); + w_sorted[i] = w[idx]; ysum += y_sorted[i] * w_sorted[i]; } @@ -232,6 +228,170 @@ float average_precision_score(const VectorXf& y, const VectorXf& predict_proba, return average_precision; } +// Binary AUROC (trapezoidal rule over the ROC curve, treating tied scores as +// a single threshold, like sklearn). `y` holds 0/1 labels and `w` per-sample +// weights. Returns 0.5 if only one class is present (AUROC is undefined). +float binary_roc_auc(const VectorXf& y, const VectorXf& predict_proba, + const vector& w) +{ + int num_instances = y.size(); + + vector order(num_instances); + iota(order.begin(), order.end(), 0); + stable_sort(order.begin(), order.end(), [&](int i, int j) { + return predict_proba(i) > predict_proba(j); // descending + }); + + float pos = 0.0f; + float neg = 0.0f; + for (int i = 0; i < num_instances; ++i) { + // remember: this is for the binary case! + pos += y(i) * w[i]; + neg += (1.0f - y(i)) * w[i]; + } + + // default case, copying sklearn, returns 0.5 if only one class exists in the y + if (pos == 0.0f || neg == 0.0f) + return 0.5f; + + float tp = 0.0f, fp = 0.0f; + float tp_prev = 0.0f, fp_prev = 0.0f; + float area = 0.0f; + for (int i = 0; i < num_instances; ++i) { + int idx = order[i]; + tp += y(idx) * w[idx]; + fp += (1.0f - y(idx)) * w[idx]; + + // only add a point to the curve at the end of a block of tied scores + bool last_of_block = (i == num_instances - 1) + || (predict_proba(order[i+1]) != predict_proba(idx)); + + if (last_of_block) { + area += (fp - fp_prev) * (tp + tp_prev) / 2.0f; + tp_prev = tp; + fp_prev = fp; + } + } + + return area / (pos * neg); +} + +// Weighted confusion matrix entries for class `label` (one-vs-rest). +void confusion(const VectorXf& y, const ArrayXi& yhat, int label, + const vector& w, float& TP, float& FP, float& FN) +{ + // Used to calculate precision and recall for multiclass settings. + // TP, FP, FN, passed as reference + + TP = FP = FN = 0.0f; + for (int i = 0; i < y.size(); ++i) { + bool is_true = static_cast(y(i)) == label; + bool is_pred = yhat(i) == label; + + if ( is_true && is_pred) TP += w[i]; + else if (!is_true && is_pred) FP += w[i]; + else if ( is_true && !is_pred) FN += w[i]; + } +} + +ArrayXi argmax_rows(const ArrayXXf& predict_proba) +{ + // converting the pred proba matrix to predictions + + ArrayXi yhat(predict_proba.rows()); + for (int i = 0; i < predict_proba.rows(); ++i) + predict_proba.row(i).maxCoeff(&yhat(i)); + + return yhat; +} + +// Macro average of precision or recall over the classes present in either +// the true or predicted labels (sklearn's default label set). +float multi_macro_precision_recall(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, const vector& class_weights, + bool precision) +{ + if (predict_proba.rows() != y.rows()) + HANDLE_ERROR_THROW("Multiclass probabilities and labels have different numbers of rows"); + + ArrayXi yhat = argmax_rows(predict_proba); + + // again setting the loss here as hit or miss, a.k.a. accuracy + loss = (yhat != y.cast().array()).cast(); + + vector w = sample_weights(y, class_weights); + + float sum = 0.0f; + int n_labels = 0; + for (int label = 0; label < predict_proba.cols(); ++label) { + bool present = (y.cast().array() == label).any() || (yhat == label).any(); + if (!present) + continue; + + float TP, FP, FN; + confusion(y, yhat, label, w, TP, FP, FN); + + float denom = precision ? TP + FP : TP + FN; + sum += denom == 0.0f ? 0.0f : TP / denom; + ++n_labels; + } + return n_labels == 0 ? 0.0f : sum / n_labels; +} + +} // anonymous namespace + +float average_precision_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, + const vector& class_weights) { + + // AP is implemented as AUC PR in sklearn. + // AP summarizes a precision-recall curve as the weighted mean of precisions + // achieved at each threshold, with the increase in recall from the previous threshold used as the weight + + // The loss vector is used in lexicase selection. we need to set something useful here + // that does make sense on individual level. Using log loss here. + loss = log_loss(y, predict_proba, class_weights); + + return binary_average_precision(y, predict_proba, sample_weights(y, class_weights)); +} + +// implementing precision_score and recall_score for the binary case. +// it will be used per-class in the multiclass case below. +float precision_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + ArrayXi yhat = (predict_proba.array() > 0.5).cast(); + + // Again updating the loss vector. Doing the same way as binary accuracy (zero_one_loss) here + loss = (yhat != y.cast().array()).cast(); + + float TP, FP, FN; + confusion(y, yhat, 1, sample_weights(y, class_weights), TP, FP, FN); + + return (TP + FP) == 0.0f ? 0.0f : TP / (TP + FP); +} + +float recall_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + ArrayXi yhat = (predict_proba.array() > 0.5).cast(); + loss = (yhat != y.cast().array()).cast(); + + float TP, FP, FN; + confusion(y, yhat, 1, sample_weights(y, class_weights), TP, FP, FN); + + return (TP + FN) == 0.0f ? 0.0f : TP / (TP + FN); +} + +float roc_auc_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + // AUROC is not decomposable per sample; log loss is used for lexicase + loss = log_loss(y, predict_proba, class_weights); + + return binary_roc_auc(y, predict_proba, sample_weights(y, class_weights)); +} + // multinomial log loss VectorXf multi_log_loss(const VectorXf& y, const ArrayXXf& predict_proba, const vector& class_weights) @@ -340,5 +500,61 @@ float multi_bal_zero_one_loss(const VectorXf& y, return present_classes == 0 ? 0.0f : recall_sum / present_classes; } +float multi_precision_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + return multi_macro_precision_recall(y, predict_proba, loss, class_weights, true); +} + +float multi_recall_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + return multi_macro_precision_recall(y, predict_proba, loss, class_weights, false); +} + +float multi_roc_auc_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + loss = multi_log_loss(y, predict_proba, class_weights); + + vector w = sample_weights(y, class_weights); + + float sum = 0.0f; + int n_labels = 0; + for (int label = 0; label < predict_proba.cols(); ++label) { + VectorXf y_bin = (y.cast().array() == label).cast(); + + // one-vs-rest AUROC is undefined if the class is absent (or is the only one) + if (y_bin.sum() == 0.0f || y_bin.sum() == y_bin.size()) + continue; + + sum += binary_roc_auc(y_bin, predict_proba.col(label).matrix(), w); + ++n_labels; + } + return n_labels == 0 ? 0.5f : sum / n_labels; +} + +float multi_average_precision_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, const vector& class_weights) +{ + loss = multi_log_loss(y, predict_proba, class_weights); + + vector w = sample_weights(y, class_weights); + + float sum = 0.0f; + int n_labels = 0; + for (int label = 0; label < predict_proba.cols(); ++label) { + VectorXf y_bin = (y.cast().array() == label).cast(); + + // recall is undefined if the class is absent + if (y_bin.sum() == 0.0f) + continue; + + sum += binary_average_precision(y_bin, predict_proba.col(label).matrix(), w); + ++n_labels; + } + return n_labels == 0 ? 0.0f : sum / n_labels; +} + } // metrics } // Brush diff --git a/src/eval/metrics.h b/src/eval/metrics.h index 1bafe7923..1c96dee35 100644 --- a/src/eval/metrics.h +++ b/src/eval/metrics.h @@ -82,9 +82,53 @@ float zero_one_loss(const VectorXf& y, const VectorXf& predict_proba, * @return The final accuracy. */ float bal_zero_one_loss(const VectorXf& y, const VectorXf& predict_proba, - VectorXf& loss, + VectorXf& loss, const vector& class_weights=vector() ); - + +/** + * @brief Precision for binary classification (threshold 0.5, positive label 1). + * @details Equivalent to sklearn's `precision_score(zero_division=0)`. Class + * weights are used as sample weights. The loss vector holds the per-sample + * misclassification indicator (used in lexicase selection). + * @param y The true labels. + * @param predict_proba The predicted probabilities. + * @param loss Reference to store the calculated losses for each sample. + * @param class_weights The optional class weights. + * @return The precision. + */ +float precision_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + +/** + * @brief Recall for binary classification (threshold 0.5, positive label 1). + * @details Equivalent to sklearn's `recall_score(zero_division=0)`. The loss + * vector holds the per-sample misclassification indicator. + * @param y The true labels. + * @param predict_proba The predicted probabilities. + * @param loss Reference to store the calculated losses for each sample. + * @param class_weights The optional class weights. + * @return The recall. + */ +float recall_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + +/** + * @brief Area under the ROC curve for binary classification. + * @details Equivalent to sklearn's `roc_auc_score`. Returns 0.5 when only one + * class is present (where the metric is undefined). The loss vector holds the + * per-sample log loss (used in lexicase selection). + * @param y The true labels. + * @param predict_proba The predicted probabilities. + * @param loss Reference to store the calculated losses for each sample. + * @param class_weights The optional class weights. + * @return The AUROC. + */ +float roc_auc_score(const VectorXf& y, const VectorXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + // multiclass classification --------------------------------------------------- /** @@ -126,6 +170,45 @@ float multi_bal_zero_one_loss(const VectorXf& y, const ArrayXXf& predict_proba, VectorXf& loss, const vector& class_weights=vector() ); +/** + * @brief Macro-averaged precision for multi-classification. + * @details Equivalent to sklearn's `precision_score(average='macro', + * zero_division=0)`: averages over classes present in either the true or the + * predicted labels. The loss vector holds the misclassification indicator. + */ +float multi_precision_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + +/** + * @brief Macro-averaged recall for multi-classification. + * @details Equivalent to sklearn's `recall_score(average='macro', + * zero_division=0)`. The loss vector holds the misclassification indicator. + */ +float multi_recall_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + +/** + * @brief Macro-averaged one-vs-rest AUROC for multi-classification. + * @details Mean of the binary AUROC of each class against the rest, skipping + * classes that are absent from `y`. The loss vector holds the per-sample + * multinomial log loss. + */ +float multi_roc_auc_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + +/** + * @brief Macro-averaged one-vs-rest average precision for multi-classification. + * @details Mean of the binary average precision of each class against the + * rest, skipping classes that are absent from `y`. The loss vector holds the + * per-sample multinomial log loss. + */ +float multi_average_precision_score(const VectorXf& y, const ArrayXXf& predict_proba, + VectorXf& loss, + const vector& class_weights=vector() ); + } // metrics } // Brush diff --git a/src/eval/scorer.h b/src/eval/scorer.h index 19a4548c6..fac32ea68 100644 --- a/src/eval/scorer.h +++ b/src/eval/scorer.h @@ -91,6 +91,9 @@ typedef float (*funcPointer)(const VectorXf&, score_hash["average_precision_score"] = &average_precision_score; score_hash["accuracy"] = &zero_one_loss; score_hash["balanced_accuracy"] = &bal_zero_one_loss; + score_hash["precision"] = &precision_score; + score_hash["recall"] = &recall_score; + score_hash["roc_auc"] = &roc_auc_score; this->set_scorer(scorer); }; @@ -171,6 +174,10 @@ typedef float (*funcPointer)(const VectorXf&, score_hash["multi_log"] = &mean_multi_log_loss; score_hash["accuracy"] = &multi_zero_one_loss; score_hash["balanced_accuracy"] = &multi_bal_zero_one_loss; + score_hash["precision"] = &multi_precision_score; + score_hash["recall"] = &multi_recall_score; + score_hash["roc_auc"] = &multi_roc_auc_score; + score_hash["average_precision_score"] = &multi_average_precision_score; this->set_scorer(scorer); }; diff --git a/src/ind/individual.h b/src/ind/individual.h index 734892ea6..78a35fbda 100644 --- a/src/ind/individual.h +++ b/src/ind/individual.h @@ -169,7 +169,10 @@ class Individual{ {"multi_log", -1.0}, {"average_precision_score", +1.0}, {"balanced_accuracy", +1.0}, - {"accuracy", +1.0} + {"accuracy", +1.0}, + {"precision", +1.0}, + {"recall", +1.0}, + {"roc_auc", +1.0} // {"scorer", -1.0} }; diff --git a/src/selection/lexicase.cpp b/src/selection/lexicase.cpp index 37aa7e062..64fae3c79 100644 --- a/src/selection/lexicase.cpp +++ b/src/selection/lexicase.cpp @@ -48,10 +48,17 @@ vector Lexicase::select(Population& pop, int island, // define epsilon ArrayXf epsilon = ArrayXf::Zero(N); - // if output is continuous, use epsilon lexicase + // if output is continuous (per sample!!), use epsilon lexicase. + // basically, every classification metric that updates the loss as hit/miss should + // not be considered here. If the clf metric updates the reference loss vector + // by assigning it float predict probas, then it will need the epsilon lexicase + // to work. + // The classification scorer names [average_precision_score, roc_auc] are the same + // for the binary and multiclassifier if (!params.classification || params.scorer.compare("log")==0 || params.scorer.compare("multi_log")==0 - || params.scorer.compare("average_precision_score")==0 ) + || params.scorer.compare("average_precision_score")==0 + || params.scorer.compare("roc_auc")==0 ) { // for each sample, calculate epsilon for (int i = 0; i Date: Tue, 22 Sep 2026 17:08:08 -0400 Subject: [PATCH 2/3] Test cases for new metrics --- src/program/operator.h | 6 +- tests/cpp/test_brush.cpp | 9 ++ tests/cpp/test_optimization.cpp | 44 ++++++++- tests/python/test_classification_metrics.py | 100 ++++++++++++++++++++ tests/python/test_final_model_selection.py | 14 ++- tests/python/test_params.py | 3 + tests/python/test_sklearn_interface.py | 4 + 7 files changed, 175 insertions(+), 5 deletions(-) create mode 100644 tests/python/test_classification_metrics.py diff --git a/src/program/operator.h b/src/program/operator.h index e84d5843f..2d9e2a3fc 100644 --- a/src/program/operator.h +++ b/src/program/operator.h @@ -235,7 +235,11 @@ struct Operator auto inputs = get_kids(d, tn, weights); if constexpr (is_one_of_v) { - if (tn.data.get_is_weighted()) + // Nodes like Floor and Ceil can have is_weighted set (it only + // checks the return type), but get_weights/set_weights skip them. + // Consuming a weight here would shift every following weight and + // read past the end of the optimizer's parameter array. + if (tn.data.get_is_weighted() && IsWeighable(NT)) { auto w = util::get_weight(tn, weights); return this->apply(inputs)*w; diff --git a/tests/cpp/test_brush.cpp b/tests/cpp/test_brush.cpp index 49b7ab05a..7e9af49a4 100644 --- a/tests/cpp/test_brush.cpp +++ b/tests/cpp/test_brush.cpp @@ -216,6 +216,15 @@ TEST_P(EngineTest, ClassificationEngineWorks) Brush::ClassifierEngine est4(params, ss); est4.run(data); + // Test the scorers that are only used for selection (never for fitting weights) + for (const auto& scorer : {"precision", "recall", "roc_auc"}) { + params.set_scorer(scorer); + std::cout << "Bandit type: " << bandit_type << std::endl; + std::cout << "Metric: " << scorer << std::endl; + Brush::ClassifierEngine est5(params, ss); + est5.run(data); + } + std::cout << "Parameters probs:" << std::endl; std::cout << "cx: " << est.params.get_cx_prob() << std::endl; for (const auto& [name, prob] : est.params.get_mutation_probs()) diff --git a/tests/cpp/test_optimization.cpp b/tests/cpp/test_optimization.cpp index d07c5fee6..0eff51099 100644 --- a/tests/cpp/test_optimization.cpp +++ b/tests/cpp/test_optimization.cpp @@ -371,4 +371,46 @@ INSTANTIATE_TEST_SUITE_P(OptimizerTestParameters, OptimizerTest, } ) ) -); \ No newline at end of file +); +TEST(Optimizer, NonWeighableNodesDoNotConsumeWeights) +{ + // Floor and Ceil can be flagged `is_weighted` (it depends only on the return + // type), but get_weights/set_weights skip them. If evaluation consumed a + // weight for them, every following weight would shift and the last one + // would be read past the end of the optimizer's parameter array, making the + // fit depend on uninitialized memory. + ArrayXXf X(10, 1); + X << 0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f; + + ArrayXf y = 2.0f * X.col(0) + 1.0f; + + Dataset data(X, y, {}, {}, {"ArrayF"}); + + // y = floor(1.0 * 1.0) + w * x_0, so the fit should find w = 2 + json PRGjson = {{"Tree", { + {{"node_type", "Add"}, {"is_weighted", false}}, + {{"node_type", "Floor"}, {"is_weighted", true}}, + {{"node_type", "Constant"}, {"is_weighted", true}, {"W", 1.0}}, + {{"node_type", "Terminal"}, {"feature", "x_0"}, {"is_weighted", true}, {"W", 1.0}} + }}, {"is_fitted_", false}}; + + RegressorProgram first = PRGjson; + ASSERT_EQ(first.get_n_weights(), 2); // Constant and Terminal, not Floor + + first.fit(data); + ArrayXf weights = first.get_weights(); + + ASSERT_TRUE(first.predict(data).isApprox(y, 1e-3)) + << "prediction=" << first.predict(data).transpose() + << ", weights=" << weights.transpose(); + ASSERT_NEAR(weights(1), 2.0f, 1e-3f); + + // fitting the same program again must give the same weights + for (int i = 0; i < 20; ++i) { + RegressorProgram again = PRGjson; + again.fit(data); + ASSERT_TRUE((again.get_weights() == weights).all()) + << "run " << i << ": " << again.get_weights().transpose() + << " vs " << weights.transpose(); + } +} diff --git a/tests/python/test_classification_metrics.py b/tests/python/test_classification_metrics.py new file mode 100644 index 000000000..e6efd92d5 --- /dev/null +++ b/tests/python/test_classification_metrics.py @@ -0,0 +1,100 @@ +import pytest +import numpy as np + +from pybrush import BrushClassifier + +from sklearn.datasets import make_classification +from sklearn.metrics import (average_precision_score, precision_score, + recall_score, roc_auc_score) + + +def _ovr_macro(binary_metric, y, proba, sample_weight): + scores = [] + for label in range(proba.shape[1]): + y_bin = (y == label).astype(float) + if 0 < y_bin.sum() < len(y_bin): + scores.append(binary_metric(y_bin, proba[:, label], sample_weight=sample_weight)) + return np.mean(scores) + + +def _sklearn_score(scorer, ind, data, n_classes, sample_weight): + y = np.array(data.y) + multiclass = n_classes > 2 + + if scorer in ["precision", "recall"]: + metric = precision_score if scorer == "precision" else recall_score + return metric(y, np.array(ind.predict(data)).astype(float), + average="macro" if multiclass else "binary", + zero_division=0, sample_weight=sample_weight) + + metric = roc_auc_score if scorer == "roc_auc" else average_precision_score + proba = np.array(ind.predict_proba(data)).astype(float) + if multiclass: + return _ovr_macro(metric, y, proba, sample_weight) + return metric(y, proba, sample_weight=sample_weight) + + +@pytest.mark.parametrize("n_classes", [2, 3]) +@pytest.mark.parametrize("scorer", ["precision", "recall", "roc_auc", "average_precision_score"]) +@pytest.mark.parametrize("class_weights", ["unbalanced", "support"]) +def test_metrics_match_sklearn(n_classes, scorer, class_weights): + X, y = make_classification(n_samples=150, n_features=6, n_informative=4, + n_classes=n_classes, weights=None, + random_state=42) + + est = BrushClassifier( + max_gens=3, + pop_size=30, + scorer=scorer, + class_weights=class_weights, + functions=['Add', 'Sub', 'Mul', 'SplitBest'], + validation_size=0.3, + final_model_selection="", + random_state=42, + verbosity=0, + ).fit(X, y) + + data = est.validation_ + y_val = np.array(data.y) + + sample_weight = None + if class_weights == "support": + classes, counts = np.unique(y_val, return_counts=True) + support = {int(c): len(y_val) / (n_classes * n) for c, n in zip(classes, counts)} + sample_weight = np.array([support[int(label)] for label in y_val]) + + brush_scores = [ind.fitness.loss_v for ind in est.archive_] + sklearn_scores = [_sklearn_score(scorer, ind, data, n_classes, sample_weight) + for ind in est.archive_] + + assert np.allclose(brush_scores, sklearn_scores, atol=1e-4), \ + f"brush={brush_scores}, sklearn={sklearn_scores}" + + +@pytest.mark.parametrize("n_classes", [2, 3]) +def test_partial_fit_changes_scorer(n_classes): + X, y = make_classification(n_samples=150, n_features=6, n_informative=4, + n_classes=n_classes, random_state=0) + + # no validation split and no class weights, so the stored fitness is the + # plain metric on the full data + est = BrushClassifier(max_gens=3, pop_size=20, scorer="roc_auc", + validation_size=0.0, class_weights="unbalanced", + random_state=0, verbosity=0).fit(X, y) + assert est.engine_.params.scorer == "roc_auc" + + est.scorer = "average_precision_score" + est.partial_fit(X, y, lock_nodes_depth=est.max_depth + 1, + keep_leaves_unlocked=False, keep_current_weights=False) + + assert est.engine_.params.scorer == "average_precision_score" + assert est.parameters_.scorer == "average_precision_score" + assert est.best_estimator_.fitness.weights[0] == +1.0 + + proba = est.predict_proba(X) + if n_classes == 2: + expected = average_precision_score(y, proba[:, 1]) + else: + expected = _ovr_macro(average_precision_score, y, proba, None) + + assert np.isclose(est.best_estimator_.fitness.loss, expected, atol=1e-4) diff --git a/tests/python/test_final_model_selection.py b/tests/python/test_final_model_selection.py index 963a272de..2c9c74cac 100644 --- a/tests/python/test_final_model_selection.py +++ b/tests/python/test_final_model_selection.py @@ -8,6 +8,7 @@ from sklearn.model_selection import GridSearchCV from sklearn.datasets import make_classification, make_regression from sklearn.metrics import mean_squared_error, log_loss, accuracy_score, balanced_accuracy_score, average_precision_score +from sklearn.metrics import precision_score, recall_score, roc_auc_score def test_smallest_complexity_selection_regression(): @@ -84,7 +85,8 @@ def test_classification_selection(): idx = np.argmin([p.fitness.linear_complexity for p in model.archive_]) -@pytest.mark.parametrize("scorer", ['log', 'accuracy', 'balanced_accuracy', 'average_precision_score']) +@pytest.mark.parametrize("scorer", ['log', 'accuracy', 'balanced_accuracy', 'average_precision_score', + 'precision', 'recall', 'roc_auc']) @pytest.mark.parametrize("class_weights", ['unbalanced', 'support', [1.0, 1.0], [1.0, 1.3]]) def test_final_model_selection_best_validation_ci_replicated(scorer, class_weights): # Small dataset for testing @@ -122,14 +124,19 @@ def test_final_model_selection_best_validation_ci_replicated(scorer, class_weigh print("Unique values in validation data", np.unique(y, return_counts=True)) + # Defining the sklearn equivalent for each metric -- we will check if they + # match later loss_f_dict = { "mse": mean_squared_error, "log": log_loss, "accuracy": accuracy_score, "balanced_accuracy": balanced_accuracy_score, "average_precision_score": average_precision_score, + "precision": lambda *a, **kw: precision_score(*a, zero_division=0, **kw), + "recall": lambda *a, **kw: recall_score(*a, zero_division=0, **kw), + "roc_auc": roc_auc_score, } - loss_f = loss_f_dict[est.parameters_.scorer] + loss_f = loss_f_dict[est.parameters_.scorer] # The reference -- should reproduce sklearn implementation values def eval_with_sklearn(individual, sample=None, log=False): @@ -141,7 +148,8 @@ def eval_with_sklearn(individual, sample=None, log=False): y_pred = None - if est.parameters_.scorer in ["log", "average_precision_score"]: + # Those scorers are also the ones that we can use in e-lexicase + if est.parameters_.scorer in ["log", "average_precision_score", "roc_auc"]: y_pred = np.array(individual.predict_proba(data)).astype(float) if est.parameters_.scorer == "log": diff --git a/tests/python/test_params.py b/tests/python/test_params.py index 1f79ddfa1..4a4a9b1fc 100644 --- a/tests/python/test_params.py +++ b/tests/python/test_params.py @@ -260,6 +260,9 @@ def test_population_split_nodes_with_and_without_SplitOn_function(): ("balanced_accuracy", [+1.0, -1.0]), # higher is better ("log", [-1.0, -1.0]), # lower is better ("average_precision_score", [+1.0, -1.0]), # higher is better + ("precision", [+1.0, -1.0]), # higher is better + ("recall", [+1.0, -1.0]), # higher is better + ("roc_auc", [+1.0, -1.0]), # higher is better ]) def test_fitness_weights_match_scorer_sign(scorer, expected_weights): """Ensure fitness.weights has correct sign according to the scorer function, diff --git a/tests/python/test_sklearn_interface.py b/tests/python/test_sklearn_interface.py index 419248fed..7abf0aa0a 100644 --- a/tests/python/test_sklearn_interface.py +++ b/tests/python/test_sklearn_interface.py @@ -67,6 +67,10 @@ def test_fixed_seed_produces_identical_brush_runs(): n_jobs=2, random_state=42, verbosity=0, + constants_simplification=False, + inexact_simplification=False, + bandit='dummy', + shuffle_split=False, ) first = BrushClassifier(**config).fit(X, y) From cb8ff1a95e99cb50d5de9346436ef64b6e07427c Mon Sep 17 00:00:00 2001 From: Guilherme Seidyo Imai Aldeia Date: Tue, 22 Sep 2026 17:26:40 -0400 Subject: [PATCH 3/3] more tests and example notebook --- docs/guide/index.md | 1 + docs/guide/switching_metrics.ipynb | 805 +++++++++++++++++++++++++++++ tests/cpp/test_evaluation.cpp | 206 ++++++++ 3 files changed, 1012 insertions(+) create mode 100644 docs/guide/switching_metrics.ipynb diff --git a/docs/guide/index.md b/docs/guide/index.md index eb09ae2aa..bbc5e43e7 100644 --- a/docs/guide/index.md +++ b/docs/guide/index.md @@ -16,6 +16,7 @@ multiclassification json saving_loading_populations locking_mechanism +switching_metrics archive deap ``` diff --git a/docs/guide/switching_metrics.ipynb b/docs/guide/switching_metrics.ipynb new file mode 100644 index 000000000..6811a337e --- /dev/null +++ b/docs/guide/switching_metrics.ipynb @@ -0,0 +1,805 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "50f0c7bc", + "metadata": {}, + "source": [ + "# Switching the classification metric with `partial_fit`\n", + "\n", + "Brush uses the `scorer` to evaluate programs: it drives selection, survival,\n", + "the archive, and the final model choice. For classification, the available\n", + "scorers are:\n", + "\n", + "| scorer | binary | multiclass |\n", + "|---|---|---|\n", + "| `\"log\"` / `\"multi_log\"` | log loss | multinomial log loss |\n", + "| `\"accuracy\"` | accuracy | accuracy |\n", + "| `\"balanced_accuracy\"` | balanced accuracy | balanced accuracy |\n", + "| `\"precision\"` | precision (threshold 0.5) | macro precision |\n", + "| `\"recall\"` | recall (threshold 0.5) | macro recall |\n", + "| `\"roc_auc\"` | AUROC | macro one-vs-rest AUROC |\n", + "| `\"average_precision_score\"` | average precision (AUPRC) | macro one-vs-rest average precision |\n", + "\n", + "The scorer is **not** used to fit parameters. Weights are always optimized with\n", + "the log loss, and split thresholds with the gini impurity. This keeps parameter\n", + "fitting smooth and well-behaved while letting you pick whichever metric you care\n", + "about for model selection.\n", + "\n", + "Because the scorer is just an estimator attribute, you can change it between\n", + "calls to `partial_fit`. This notebook:\n", + "\n", + "1. Fits a model on an imbalanced binary problem using AUROC.\n", + "2. Locks every internal node of the best program, leaving only the leaves and\n", + " the weights free to change.\n", + "3. Switches the scorer to average precision and calls `partial_fit`.\n", + "4. Compares all metrics before and after the switch." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "37c250d0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:04.638378Z", + "iopub.status.busy": "2026-09-22T19:12:04.638227Z", + "iopub.status.idle": "2026-09-22T19:12:05.855283Z", + "shell.execute_reply": "2026-09-22T19:12:05.854853Z" + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import graphviz\n", + "\n", + "from sklearn.datasets import make_classification\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.metrics import (log_loss, accuracy_score, balanced_accuracy_score,\n", + " precision_score, recall_score, roc_auc_score,\n", + " average_precision_score)\n", + "\n", + "from pybrush import BrushClassifier" + ] + }, + { + "cell_type": "markdown", + "id": "db851330", + "metadata": {}, + "source": [ + "## 1. An imbalanced binary problem\n", + "\n", + "AUROC and average precision disagree the most when positives are rare: AUROC\n", + "rewards ranking negatives correctly, while average precision focuses on how\n", + "clean the top of the ranking is." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c6262431", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:05.856401Z", + "iopub.status.busy": "2026-09-22T19:12:05.856312Z", + "iopub.status.idle": "2026-09-22T19:12:05.865082Z", + "shell.execute_reply": "2026-09-22T19:12:05.864748Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "train prevalence: 0.109\n", + "test prevalence: 0.107\n" + ] + } + ], + "source": [ + "X, y = make_classification(\n", + " n_samples=1000, n_features=6, n_informative=4, n_redundant=1,\n", + " weights=[0.9, 0.1], class_sep=0.8, flip_y=0.02, random_state=42,\n", + ")\n", + "\n", + "X_train, X_test, y_train, y_test = train_test_split(\n", + " X, y, test_size=0.3, stratify=y, random_state=42)\n", + "\n", + "print('train prevalence:', y_train.mean().round(3))\n", + "print('test prevalence: ', y_test.mean().round(3))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ed5afaf0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:05.865913Z", + "iopub.status.busy": "2026-09-22T19:12:05.865867Z", + "iopub.status.idle": "2026-09-22T19:12:05.867551Z", + "shell.execute_reply": "2026-09-22T19:12:05.867314Z" + } + }, + "outputs": [], + "source": [ + "def all_metrics(est, X, y):\n", + " proba = est.predict_proba(X)[:, 1]\n", + " pred = est.predict(X)\n", + " return {\n", + " 'log_loss': log_loss(y, proba),\n", + " 'accuracy': accuracy_score(y, pred),\n", + " 'balanced_accuracy': balanced_accuracy_score(y, pred),\n", + " 'precision': precision_score(y, pred, zero_division=0),\n", + " 'recall': recall_score(y, pred, zero_division=0),\n", + " 'roc_auc': roc_auc_score(y, proba),\n", + " 'average_precision': average_precision_score(y, proba),\n", + " }\n", + "\n", + "def report(est, label):\n", + " return pd.DataFrame({\n", + " (label, 'train'): all_metrics(est, X_train, y_train),\n", + " (label, 'test'): all_metrics(est, X_test, y_test),\n", + " })" + ] + }, + { + "cell_type": "markdown", + "id": "c072f652", + "metadata": {}, + "source": [ + "## 2. Fit using AUROC as the scorer" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "bfa036f8", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:05.868283Z", + "iopub.status.busy": "2026-09-22T19:12:05.868241Z", + "iopub.status.idle": "2026-09-22T19:12:14.024733Z", + "shell.execute_reply": "2026-09-22T19:12:14.024331Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "scorer: roc_auc\n", + "internal scorer value (train): 0.9027\n", + "model: Logistic(Add(0.98,Add(-1.05*Mul(-1.32*x_4,-1.01*x_2),x_0)))\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
roc_auc
traintest
log_loss0.3982100.387116
accuracy0.8100000.843333
balanced_accuracy0.8009870.871035
precision0.3389830.397260
recall0.7894740.906250
roc_auc0.8869560.938666
average_precision0.6782970.789859
\n", + "
" + ], + "text/plain": [ + " roc_auc \n", + " train test\n", + "log_loss 0.398210 0.387116\n", + "accuracy 0.810000 0.843333\n", + "balanced_accuracy 0.800987 0.871035\n", + "precision 0.338983 0.397260\n", + "recall 0.789474 0.906250\n", + "roc_auc 0.886956 0.938666\n", + "average_precision 0.678297 0.789859" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "est = BrushClassifier(\n", + " functions=['SplitBest', 'Add', 'Sub', 'Mul', 'Div', 'Logabs', 'Exp'],\n", + " scorer='roc_auc',\n", + " pop_size=200,\n", + " max_gens=30,\n", + " max_depth=5,\n", + " max_size=20,\n", + " random_state=42,\n", + " verbosity=0,\n", + ")\n", + "est.fit(X_train, y_train)\n", + "\n", + "print('scorer:', est.parameters_.scorer)\n", + "print('internal scorer value (train):', round(est.best_estimator_.fitness.loss, 4))\n", + "print('model:', est.best_estimator_.get_model())\n", + "\n", + "before = report(est, 'roc_auc')\n", + "before" + ] + }, + { + "cell_type": "markdown", + "id": "3aa6132b", + "metadata": {}, + "source": [ + "## 3. Lock the internal nodes, switch to average precision, and refit\n", + "\n", + "A `lock_nodes_depth` larger than the tree depth locks every operator and split\n", + "in the program. With `keep_leaves_unlocked=True`, the leaves (features and\n", + "constants) stay free, so the search can still swap a leaf or grow it into a\n", + "small subtree. With `keep_current_weights=False`, the weights are re-optimized.\n", + "\n", + "We then set `est.scorer = 'average_precision_score'` before calling\n", + "`partial_fit`. The new scorer is picked up by the engine, so every program is\n", + "re-evaluated, selected, and archived by average precision from here on." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b33832db", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:14.025731Z", + "iopub.status.busy": "2026-09-22T19:12:14.025682Z", + "iopub.status.idle": "2026-09-22T19:12:20.981227Z", + "shell.execute_reply": "2026-09-22T19:12:20.980941Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "scorer: average_precision_score\n", + "internal scorer value (train): 0.917\n", + "model before: Logistic(Add(0.98,Add(-1.05*Mul(-1.32*x_4,-1.01*x_2),x_0)))\n", + "model after: Logistic(Add(-2.77,Add(-1.53*Mul(x_4,x_2),3.87*Exp(0.42*x_0))))\n" + ] + } + ], + "source": [ + "structure_before = est.best_estimator_.get_model()\n", + "dot_before = est.best_estimator_.get_model('dot')\n", + "\n", + "est.scorer = 'average_precision_score'\n", + "est.partial_fit(\n", + " X_train, y_train,\n", + " lock_nodes_depth=est.max_depth + 1, # deeper than any tree: lock every internal node\n", + " keep_leaves_unlocked=True, # ...but leave the leaves free to change\n", + " keep_current_weights=False, # weights can still be optimized\n", + ")\n", + "\n", + "print('scorer:', est.parameters_.scorer)\n", + "print('internal scorer value (train):', round(est.best_estimator_.fitness.loss, 4))\n", + "print('model before:', structure_before)\n", + "print('model after: ', est.best_estimator_.get_model())\n", + "\n", + "after = report(est, 'average_precision_score')" + ] + }, + { + "cell_type": "markdown", + "id": "ed1c72d3", + "metadata": {}, + "source": [ + "The two programs side by side. The locked internal nodes are the same in\n", + "both, and only the leaves (and the weights) differ." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "f7a39b67", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:20.982237Z", + "iopub.status.busy": "2026-09-22T19:12:20.982178Z", + "iopub.status.idle": "2026-09-22T19:12:21.166282Z", + "shell.execute_reply": "2026-09-22T19:12:21.165925Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + "

Before (AUROC)

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "G\n", + "\n", + "^ split feature fixed, * split threshold fixed\n", + "\n", + "\n", + "1473c7000\n", + "\n", + "Logistic\n", + "\n", + "\n", + "\n", + "105024bc0\n", + "\n", + "Add\n", + "\n", + "\n", + "\n", + "1473c7000->105024bc0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "1473c70b0\n", + "\n", + "Add\n", + "\n", + "\n", + "\n", + "105024bc0->1473c70b0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "105024bc0Offset\n", + "\n", + "0.98\n", + "\n", + "\n", + "\n", + "105024bc0->105024bc0Offset\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "1473884b0\n", + "\n", + "Mul\n", + "\n", + "\n", + "\n", + "1473c70b0->1473884b0\n", + "\n", + "\n", + "-1.05\n", + "\n", + "\n", + "\n", + "x_0\n", + "\n", + "x_0\n", + "\n", + "\n", + "\n", + "1473c70b0->x_0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "x_4\n", + "\n", + "x_4\n", + "\n", + "\n", + "\n", + "1473884b0->x_4\n", + "\n", + "\n", + "-1.32\n", + "\n", + "\n", + "\n", + "x_2\n", + "\n", + "x_2\n", + "\n", + "\n", + "\n", + "1473884b0->x_2\n", + "\n", + "\n", + "-1.01\n", + "\n", + "\n", + "\n", + "
\n", + "

After (average precision)

\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "G\n", + "\n", + "^ split feature fixed, * split threshold fixed\n", + "\n", + "\n", + "15f54b570\n", + "\n", + "Logistic\n", + "\n", + "\n", + "\n", + "15f54ba00\n", + "\n", + "Add\n", + "\n", + "\n", + "\n", + "15f54b570->15f54ba00\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "15f5c09c0\n", + "\n", + "Add\n", + "\n", + "\n", + "\n", + "15f54ba00->15f5c09c0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "15f54ba00Offset\n", + "\n", + "-2.77\n", + "\n", + "\n", + "\n", + "15f54ba00->15f54ba00Offset\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "15f5c0a70\n", + "\n", + "Mul\n", + "\n", + "\n", + "\n", + "15f5c09c0->15f5c0a70\n", + "\n", + "\n", + "-1.53\n", + "\n", + "\n", + "\n", + "15f5f0460\n", + "\n", + "Exp\n", + "\n", + "\n", + "\n", + "15f5c09c0->15f5f0460\n", + "\n", + "\n", + "3.87\n", + "\n", + "\n", + "\n", + "x_4\n", + "\n", + "x_4\n", + "\n", + "\n", + "\n", + "15f5c0a70->x_4\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "x_2\n", + "\n", + "x_2\n", + "\n", + "\n", + "\n", + "15f5c0a70->x_2\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "x_0\n", + "\n", + "x_0\n", + "\n", + "\n", + "\n", + "15f5f0460->x_0\n", + "\n", + "\n", + "0.42\n", + "\n", + "\n", + "\n", + "
\n", + "
\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from IPython.display import HTML\n", + "\n", + "def svg(dot):\n", + " return graphviz.Source(dot).pipe(format='svg').decode()\n", + "\n", + "HTML(f\"\"\"\n", + "
\n", + "

Before (AUROC)

{svg(dot_before)}
\n", + "

After (average precision)

{svg(est.best_estimator_.get_model('dot'))}
\n", + "
\n", + "\"\"\")" + ] + }, + { + "cell_type": "markdown", + "id": "26aa005f", + "metadata": {}, + "source": [ + "## 4. Compare every metric before and after the switch" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "32d527ab", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-22T19:12:21.167352Z", + "iopub.status.busy": "2026-09-22T19:12:21.167284Z", + "iopub.status.idle": "2026-09-22T19:12:21.173975Z", + "shell.execute_reply": "2026-09-22T19:12:21.173673Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
roc_aucaverage_precision_scoredelta
traintesttraintesttraintest
log_loss0.39820.38710.38190.3787-0.0163-0.0084
accuracy0.81000.84330.84000.86330.03000.0200
balanced_accuracy0.80100.87100.81780.88220.01680.0112
precision0.33900.39730.38460.43280.04560.0356
recall0.78950.90620.78950.90620.00000.0000
roc_auc0.88700.93870.88920.93970.00220.0010
average_precision0.67830.78990.70410.79830.02580.0085
\n", + "
" + ], + "text/plain": [ + " roc_auc average_precision_score delta \\\n", + " train test train test train \n", + "log_loss 0.3982 0.3871 0.3819 0.3787 -0.0163 \n", + "accuracy 0.8100 0.8433 0.8400 0.8633 0.0300 \n", + "balanced_accuracy 0.8010 0.8710 0.8178 0.8822 0.0168 \n", + "precision 0.3390 0.3973 0.3846 0.4328 0.0456 \n", + "recall 0.7895 0.9062 0.7895 0.9062 0.0000 \n", + "roc_auc 0.8870 0.9387 0.8892 0.9397 0.0022 \n", + "average_precision 0.6783 0.7899 0.7041 0.7983 0.0258 \n", + "\n", + " \n", + " test \n", + "log_loss -0.0084 \n", + "accuracy 0.0200 \n", + "balanced_accuracy 0.0112 \n", + "precision 0.0356 \n", + "recall 0.0000 \n", + "roc_auc 0.0010 \n", + "average_precision 0.0085 " + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "comparison = pd.concat([before, after], axis=1)\n", + "comparison[('delta', 'train')] = comparison[('average_precision_score', 'train')] - comparison[('roc_auc', 'train')]\n", + "comparison[('delta', 'test')] = comparison[('average_precision_score', 'test')] - comparison[('roc_auc', 'test')]\n", + "comparison.round(4)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "brush", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/cpp/test_evaluation.cpp b/tests/cpp/test_evaluation.cpp index 7d748613f..50c92bcf1 100644 --- a/tests/cpp/test_evaluation.cpp +++ b/tests/cpp/test_evaluation.cpp @@ -116,6 +116,212 @@ TEST(Evaluation, MulticlassSoftmaxAndMetrics) EXPECT_NEAR(multi_bal_zero_one_loss(y, probabilities, loss), 1.0f, 1e-6f); } +// Expected values in the tests below were computed with sklearn.metrics +// (precision_score, recall_score, roc_auc_score, average_precision_score). +class BinaryMetrics : public ::testing::Test { +protected: + VectorXf y, proba, misclassified; + + void SetUp() override { + y.resize(10); proba.resize(10); misclassified.resize(10); + + y << 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0; + proba << 0.1, 0.9, 0.4, 0.6, 0.8, 0.3, 0.7, 0.2, 0.05, 0.55; + + // threshold 0.5 -> TP = 3 (idx 1, 4, 9), FP = 2 (idx 3, 6), FN = 2 (idx 2, 7) + misclassified << 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0; + } +}; + +TEST_F(BinaryMetrics, Precision) +{ + VectorXf loss; + EXPECT_NEAR(precision_score(y, proba, loss), 0.6f, 1e-6f); // 3 / (3 + 2) + ASSERT_TRUE(loss == misclassified); + + // class weights act as sample weights: 3*2 / (3*2 + 2*1) + EXPECT_NEAR(precision_score(y, proba, loss, {1.0f, 2.0f}), 0.75f, 1e-6f); + // the loss vector is not weighted (it is used by lexicase) + ASSERT_TRUE(loss == misclassified); +} + +TEST_F(BinaryMetrics, Recall) +{ + VectorXf loss; + EXPECT_NEAR(recall_score(y, proba, loss), 0.6f, 1e-6f); // 3 / (3 + 2) + ASSERT_TRUE(loss == misclassified); + + // recall only looks at positives, so class weights cancel out + EXPECT_NEAR(recall_score(y, proba, loss, {1.0f, 2.0f}), 0.6f, 1e-6f); +} + +TEST_F(BinaryMetrics, PrecisionRecallWithoutPredictedPositives) +{ + // nothing is predicted as positive: zero_division=0, like sklearn + VectorXf loss; + VectorXf low = VectorXf::Constant(10, 0.1f); + EXPECT_NEAR(precision_score(y, low, loss), 0.0f, 1e-6f); + EXPECT_NEAR(recall_score(y, low, loss), 0.0f, 1e-6f); + ASSERT_TRUE(loss == y); // every positive is a miss +} + +TEST_F(BinaryMetrics, RocAuc) +{ + VectorXf loss; + EXPECT_NEAR(roc_auc_score(y, proba, loss), 0.72f, 1e-6f); + + // per-sample loss is the log loss + VectorXf expected_loss = log_loss(y, proba); + ASSERT_TRUE(loss.isApprox(expected_loss, 1e-6f)); + + // AUROC is invariant to scaling all weights of one class + EXPECT_NEAR(roc_auc_score(y, proba, loss, {1.0f, 2.0f}), 0.72f, 1e-6f); +} + +TEST_F(BinaryMetrics, RocAucTiedScores) +{ + // tied scores are a single threshold (trapezoid, not a staircase) + VectorXf loss, tied(10); + tied << 0.2, 0.8, 0.8, 0.2, 0.8, 0.2, 0.8, 0.2, 0.2, 0.5; + EXPECT_NEAR(roc_auc_score(y, tied, loss), 0.78f, 1e-6f); + + // constant predictions carry no ranking information + VectorXf constant = VectorXf::Constant(10, 0.5f); + EXPECT_NEAR(roc_auc_score(y, constant, loss), 0.5f, 1e-6f); +} + +TEST_F(BinaryMetrics, RocAucEdgeCases) +{ + VectorXf loss, perfect(10), reversed(10); + perfect = y * 0.8f + VectorXf::Constant(10, 0.1f); + reversed = VectorXf::Constant(10, 1.0f) - perfect; + + EXPECT_NEAR(roc_auc_score(y, perfect, loss), 1.0f, 1e-6f); + EXPECT_NEAR(roc_auc_score(y, reversed, loss), 0.0f, 1e-6f); + + // undefined with a single class: we return 0.5 instead of throwing + VectorXf ones = VectorXf::Ones(10); + EXPECT_NEAR(roc_auc_score(ones, proba, loss), 0.5f, 1e-6f); +} + +TEST(Evaluation, ScorerBinaryNewMetrics) +{ + VectorXf y(4), yhat(4), loss_expected, loss; + y << 0.0, 1.0, 1.0, 0.0; + yhat << 0.1, 0.9, 0.2, 0.8; + + Scorer scorer("precision"); + ASSERT_NEAR(scorer.score(y, yhat, loss, {}), precision_score(y, yhat, loss_expected), 1e-6); + ASSERT_TRUE(loss.isApprox(loss_expected, 1e-6)); + + scorer.set_scorer("recall"); + ASSERT_NEAR(scorer.score(y, yhat, loss, {}), recall_score(y, yhat, loss_expected), 1e-6); + ASSERT_TRUE(loss.isApprox(loss_expected, 1e-6)); + + scorer.set_scorer("roc_auc"); + ASSERT_NEAR(scorer.score(y, yhat, loss, {}), roc_auc_score(y, yhat, loss_expected), 1e-6); + ASSERT_TRUE(loss.isApprox(loss_expected, 1e-6)); +} + +class MulticlassMetrics : public ::testing::Test { +protected: + VectorXf y; + ArrayXXf proba; + + void SetUp() override { + y.resize(6); proba.resize(6, 3); + + y << 0.0, 1.0, 2.0, 0.0, 1.0, 2.0; + proba << 0.7, 0.2, 0.1, + 0.3, 0.4, 0.3, + 0.2, 0.5, 0.3, // predicts 1, true 2 + 0.4, 0.35, 0.25, + 0.1, 0.3, 0.6, // predicts 2, true 1 + 0.1, 0.1, 0.8; + } +}; + +TEST_F(MulticlassMetrics, PrecisionRecall) +{ + VectorXf loss, misclassified(6); + misclassified << 0.0, 0.0, 1.0, 0.0, 1.0, 0.0; + + // per class precision = recall = {1, 0.5, 0.5} + EXPECT_NEAR(multi_precision_score(y, proba, loss), 2.0f/3.0f, 1e-6f); + ASSERT_TRUE(loss == misclassified); + + EXPECT_NEAR(multi_recall_score(y, proba, loss), 2.0f/3.0f, 1e-6f); + ASSERT_TRUE(loss == misclassified); +} + +TEST_F(MulticlassMetrics, RocAucAndAveragePrecision) +{ + VectorXf loss; + VectorXf expected_loss = multi_log_loss(y, proba); + + EXPECT_NEAR(multi_roc_auc_score(y, proba, loss), 0.8125f, 1e-6f); + ASSERT_TRUE(loss.isApprox(expected_loss, 1e-6f)); + + // weights change the one-vs-rest problems, since "rest" mixes classes + EXPECT_NEAR(multi_roc_auc_score(y, proba, loss, {1.0f, 2.0f, 3.0f}), + 0.7708333f, 1e-5f); + + EXPECT_NEAR(multi_average_precision_score(y, proba, loss), 0.75f, 1e-6f); + ASSERT_TRUE(loss.isApprox(expected_loss, 1e-6f)); +} + +TEST_F(MulticlassMetrics, AbsentClass) +{ + VectorXf loss, y_absent(6); + y_absent << 0.0, 1.0, 1.0, 0.0, 1.0, 0.0; // class 2 never occurs + + // ranking metrics skip classes that are absent from y + EXPECT_NEAR(multi_roc_auc_score(y_absent, proba, loss), 0.8055556f, 1e-5f); + EXPECT_NEAR(multi_average_precision_score(y_absent, proba, loss), 0.875f, 1e-5f); + + // precision/recall average over classes in y or in the predictions (class + // 2 is predicted once, so it counts with precision = recall = 0) + EXPECT_NEAR(multi_precision_score(y_absent, proba, loss), 2.0f/3.0f, 1e-6f); + EXPECT_NEAR(multi_recall_score(y_absent, proba, loss), 4.0f/9.0f, 1e-6f); +} + +TEST_F(MulticlassMetrics, PerfectPredictions) +{ + VectorXf loss; + ArrayXXf perfect = ArrayXXf::Constant(6, 3, 0.1f); + for (int i = 0; i < y.size(); ++i) + perfect(i, static_cast(y(i))) = 0.8f; + + EXPECT_NEAR(multi_precision_score(y, perfect, loss), 1.0f, 1e-6f); + EXPECT_NEAR(multi_recall_score(y, perfect, loss), 1.0f, 1e-6f); + EXPECT_NEAR(multi_roc_auc_score(y, perfect, loss), 1.0f, 1e-6f); + EXPECT_NEAR(multi_average_precision_score(y, perfect, loss), 1.0f, 1e-6f); +} + +TEST(Evaluation, ScorerMulticlassNewMetrics) +{ + VectorXf y(3), loss_expected, loss; + y << 0.0, 1.0, 2.0; + ArrayXXf proba(3, 3); + proba << 0.6, 0.3, 0.1, + 0.5, 0.3, 0.2, + 0.1, 0.2, 0.7; + + Scorer scorer("precision"); + ASSERT_NEAR(scorer.score(y, proba, loss, {}), multi_precision_score(y, proba, loss_expected), 1e-6); + ASSERT_TRUE(loss.isApprox(loss_expected, 1e-6)); + + scorer.set_scorer("recall"); + ASSERT_NEAR(scorer.score(y, proba, loss, {}), multi_recall_score(y, proba, loss_expected), 1e-6); + + scorer.set_scorer("roc_auc"); + ASSERT_NEAR(scorer.score(y, proba, loss, {}), multi_roc_auc_score(y, proba, loss_expected), 1e-6); + ASSERT_TRUE(loss.isApprox(loss_expected, 1e-6)); + + scorer.set_scorer("average_precision_score"); + ASSERT_NEAR(scorer.score(y, proba, loss, {}), multi_average_precision_score(y, proba, loss_expected), 1e-6); +} + TEST(Evaluation, MulticlassSoftmaxHasOneOutputPerClass) { ArrayXXf X(6, 2);