<commit_before>"""
Different loss functions for gradient boosting are defined here.
"""
from __future__ import division, print_function, absolute_import
import numbers
import warnings
from collections import defaultdict

import numpy
import pandas
from scipy import sparse
from scipy.special import expit
from sklearn.utils.validation import check_random_state
from sklearn.base import BaseEstimator

from .commonutils import compute_knn_indices_of_signal, check_sample_weight, check_uniform_label
from .metrics_utils import bin_to_group_indices, compute_bin_indices, compute_group_weights, \
    group_indices_to_groups_matrix


__author__ = '<NAME>'


def compute_positions(y_pred, sample_weight):
    """
    For each event computes it position among other events by prediction.
    position = (weighted) part of elements with lower predictions => position belongs to [0, 1]

    This function is very close to `scipy.stats.rankdata`, but supports weights.
    """
    order = numpy.argsort(y_pred)
    ordered_weights = sample_weight[order]
    ordered_weights /= float(numpy.sum(ordered_weights))
    efficiencies = (numpy.cumsum(ordered_weights) - 0.5 * ordered_weights)
    return efficiencies[numpy.argsort(order)]


class AbstractLossFunction(BaseEstimator):
    """
    This is base class for loss functions used in `hep_ml`.
    Main differences compared to `scikit-learn` loss functions:
    1. losses are stateful, and may require fitting of training data before usage.
    2. thus, when computing gradients, hessians, one shall provide predictions of all events.
    3. losses are object that shall be passed as estimators to gradient boosting (see examples).
    4. only two-class case is supported, and different classes may have different role and meaning.
       (so no class-switching like possible as this was done in sklearn).
    """

    def fit(self, X, y, sample_weight):
        """ This method is optional, it is called before all the others."""
        return self

    def negative_gradient(self, y_pred):
        """The y_pred should contain all the events passed to `fit` method,
        moreover, the order should be the same"""
        raise NotImplementedError()

    def __call__(self, y_pred):
        """The y_pred should contain all the events passed to `fit` method,
        moreover, the order should be the same"""
        raise NotImplementedError()

    def prepare_tree_params(self, y_pred):
        """Prepares patameters for regression tree that minimizes MSE
        :param y_pred: contains predictions for all the events passed to `fit` method,
        moreover, the order should be the same
        :return: tuple (residual, sample_weight) with target and weight to be used in decision tree
        """
        return self.negative_gradient(y_pred), numpy.ones(len(y_pred))

    def prepare_new_leaves_values(self, terminal_regions, leaf_values,
                                  X, y, y_pred, sample_weight, update_mask, residual):
        """
        Method for pruning. Loss function can prepare better values for leaves

        :param terminal_regions: indices of terminal regions of each event.
        :param leaf_values: numpy.array, current mapping of leaf indices to prediction values.
        :param X: data (same as passed in fit, ignored usually)
        :param y: labels (same as passed in fit)
        :param y_pred: predictions before adding new tree.
        :param sample_weight: weights od samples (same as passed in fit)
        :param update_mask: which events to use during update?
        :param residual: computed value of negative gradient (before adding tree)
        :return: numpy.array with new prediction values for all leaves.
        """
        return leaf_values


class HessianLossFunction(AbstractLossFunction):
    def __init__(self, regularization=5.):
        """Loss function with diagonal hessian, provides uses Newton-Raphson step to update trees.
        :param regularization: float, penalty for leaves with few events,
            corresponds roughly to the number of added events of both classes to each leaf.
        """
        self.regularization = regularization

    def fit(self, X, y, sample_weight):
        self.regularization_ = self.regularization * numpy.mean(sample_weight)
        return self

    def hessian(self, y_pred):
        raise NotImplementedError('Override this method in loss function.')

    def prepare_tree_params(self, y_pred):
        grad = self.negative_gradient(y_pred)
        hess = self.hessian(y_pred) + 0.01
        return grad / hess, hess

    def prepare_new_leaves_values(self, terminal_regions, leaf_values,
                                  X, y, y_pred, sample_weight, update_mask, residual):
        """ This expression comes from optimization of second-order approximation of loss function."""
        minlength = len(leaf_values)
        nominators = numpy.bincount(terminal_regions, weights=residual, minlength=minlength)
        denominators = numpy.bincount(terminal_regions, weights=self.hessian(y_pred), minlength=minlength)
        return nominators / (denominators + self.regularization_)


class AdaLossFunction(HessianLossFunction):
    """ AdaLossFunction is the same as Exponential Loss Function (exploss)"""
    # TODO write better update
    def fit(self, X, y, sample_weight):
        self.y = y
        self.sample_weight = sample_weight
        self.y_signed = 2 * y - 1
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)
        return self

    def __call__(self, y_pred):
        return numpy.sum(self.sample_weight * numpy.exp(- self.y_signed * y_pred))

    def negative_gradient(self, y_pred):
        return self.y_signed * self.sample_weight * numpy.exp(- self.y_signed * y_pred)

    def hessian(self, y_pred):
        return self.sample_weight * numpy.exp(- self.y_signed * y_pred)


class BinomialDevianceLossFunction(HessianLossFunction):
    """
    Binomial deviance, aka Logistic loss function (logloss).
    """

    def fit(self, X, y, sample_weight):
        self.y = y
        self.sample_weight = sample_weight
        self.y_signed = 2 * y - 1
        self.adjusted_regularization = numpy.mean(sample_weight) * self.regularization
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)
        return self

    def __call__(self, y_pred):
        return numpy.sum(self.sample_weight * numpy.logaddexp(0, - self.y_signed * y_pred))

    def negative_gradient(self, y_pred):
        return self.y_signed * self.sample_weight * expit(- self.y_signed * y_pred)

    def hessian(self, y_pred):
        expits = expit(self.y_signed * y_pred)
        return self.sample_weight * expits * (1 - expits)


class CompositeLossFunction(HessianLossFunction):
    """
    This is exploss for bck and logloss for signal with proper constants.
    Such kind of loss functions is very useful to optimize AMS.
    """

    def fit(self, X, y, sample_weight):
        self.y = y
        self.sample_weight = sample_weight
        self.y_signed = 2 * y - 1
        self.sig_w = (y == 1) * self.sample_weight
        self.bck_w = (y == 0) * self.sample_weight
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)
        return self

    def __call__(self, y_pred):
        result = numpy.sum(self.sig_w * numpy.logaddexp(0, -y_pred))
        result += numpy.sum(self.bck_w * numpy.exp(0.5 * y_pred))
        return result

    def negative_gradient(self, y_pred):
        result = self.sig_w * expit(- y_pred)
        result -= 0.5 * self.bck_w * numpy.exp(0.5 * y_pred)
        return result

    def hessian(self, y_pred):
        expits = expit(- y_pred)
        return self.sig_w * expits * (1 - expits) + self.bck_w * 0.25 * numpy.exp(0.5 * y_pred)


class RankBoostLossFunction(HessianLossFunction):
    def __init__(self, request_column, messup_penalty='square', regularization=0.1, update_terations=1):
        self.update_terations = update_terations
        self.messup_penalty = messup_penalty
        self.request_column = request_column
        HessianLossFunction.__init__(self, regularization=regularization)

    def fit(self, X, y, sample_weight):
        self.queries = X[self.request_column]
        self.y = y
        self.possible_queries, normed_queries = numpy.unique(self.queries, return_inverse=True)
        self.possible_ranks, normed_ranks = numpy.unique(self.y, return_inverse=True)

        self.lookups = [normed_ranks, normed_queries * len(self.possible_ranks) + normed_ranks]
        self.minlengths = [len(self.possible_ranks), len(self.possible_ranks) * len(self.possible_queries)]
        self.rank_penalties = numpy.zeros([len(self.possible_ranks), len(self.possible_ranks)], dtype=float)
        for r1 in self.possible_ranks:
            for r2 in self.possible_ranks:
                if r1 < r2:
                    if self.messup_penalty == 'square':
                        self.rank_penalties[r1, r2] = (r2 - r1) ** 2
                    elif self.messup_penalty == 'linear':
                        self.rank_penalties[r1, r2] = r2 - r1
                    else:
                        raise NotImplementedError()

        self.penalty_matrices = []
        self.penalty_matrices.append(self.rank_penalties / 10.)
        self.penalty_matrices.append(sparse.block_diag([self.rank_penalties] * len(self.possible_queries)))
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)

    def __call__(self, y_pred):
        """
        loss is defined as  w_ij exp(pred_i - pred_j),
        w_ij is zero if label_i <= label_j
        All the other labels are:

        w_ij = (alpha + beta * [query_i = query_j]) rank_penalty_{type_i, type_j}
        rank_penalty_{ij} is zero if i <= j

        :param y_pred: predictions of shape [n_samples]
        :return: value of loss, float
        """
        y_pred -= y_pred.mean()
        pos_exponent = numpy.exp(y_pred)
        neg_exponent = numpy.exp(-y_pred)
        result = 0.
        for lookup, length, penalty_matrix in zip(self.lookups, self.minlengths, self.penalty_matrices):
            pos_stats = numpy.bincount(lookup, weights=pos_exponent)
            neg_stats = numpy.bincount(lookup, weights=neg_exponent)
            result += pos_stats.T.dot(penalty_matrix.dot(neg_stats))
        # assert numpy.shape(result) == tuple()
        return result

    def negative_gradient(self, y_pred):
        y_pred -= y_pred.mean()
        pos_exponent = numpy.exp(y_pred)
        neg_exponent = numpy.exp(-y_pred)
        gradient = numpy.zeros(len(y_pred), dtype=float)
        for lookup, length, penalty_matrix in zip(self.lookups, self.minlengths, self.penalty_matrices):
            pos_stats = numpy.bincount(lookup, weights=pos_exponent)
            neg_stats = numpy.bincount(lookup, weights=neg_exponent)
            gradient += pos_exponent * penalty_matrix.dot(neg_stats)[lookup]
            gradient -= neg_exponent * penalty_matrix.T.dot(pos_stats)[lookup]
        return - gradient

    def hessian(self, y_pred):
        y_pred -= y_pred.mean()
        pos_exponent = numpy.exp(y_pred)
        neg_exponent = numpy.exp(-y_pred)
        result = numpy.zeros(len(y_pred), dtype=float)<commit_msg>added per-query normalization to RankBoost<commit_after>"""
Different loss functions for gradient boosting are defined here.
"""
from __future__ import division, print_function, absolute_import
import numbers
import warnings
from collections import defaultdict

import numpy
import pandas
from scipy import sparse
from scipy.special import expit
from sklearn.utils.validation import check_random_state
from sklearn.base import BaseEstimator

from .commonutils import compute_knn_indices_of_signal, check_sample_weight, check_uniform_label
from .metrics_utils import bin_to_group_indices, compute_bin_indices, compute_group_weights, \
    group_indices_to_groups_matrix


__author__ = '<NAME>'


def compute_positions(y_pred, sample_weight):
    """
    For each event computes it position among other events by prediction.
    position = (weighted) part of elements with lower predictions => position belongs to [0, 1]

    This function is very close to `scipy.stats.rankdata`, but supports weights.
    """
    order = numpy.argsort(y_pred)
    ordered_weights = sample_weight[order]
    ordered_weights /= float(numpy.sum(ordered_weights))
    efficiencies = (numpy.cumsum(ordered_weights) - 0.5 * ordered_weights)
    return efficiencies[numpy.argsort(order)]


class AbstractLossFunction(BaseEstimator):
    """
    This is base class for loss functions used in `hep_ml`.
    Main differences compared to `scikit-learn` loss functions:
    1. losses are stateful, and may require fitting of training data before usage.
    2. thus, when computing gradients, hessians, one shall provide predictions of all events.
    3. losses are object that shall be passed as estimators to gradient boosting (see examples).
    4. only two-class case is supported, and different classes may have different role and meaning.
       (so no class-switching like possible as this was done in sklearn).
    """

    def fit(self, X, y, sample_weight):
        """ This method is optional, it is called before all the others."""
        return self

    def negative_gradient(self, y_pred):
        """The y_pred should contain all the events passed to `fit` method,
        moreover, the order should be the same"""
        raise NotImplementedError()

    def __call__(self, y_pred):
        """The y_pred should contain all the events passed to `fit` method,
        moreover, the order should be the same"""
        raise NotImplementedError()

    def prepare_tree_params(self, y_pred):
        """Prepares patameters for regression tree that minimizes MSE
        :param y_pred: contains predictions for all the events passed to `fit` method,
        moreover, the order should be the same
        :return: tuple (residual, sample_weight) with target and weight to be used in decision tree
        """
        return self.negative_gradient(y_pred), numpy.ones(len(y_pred))

    def prepare_new_leaves_values(self, terminal_regions, leaf_values,
                                  X, y, y_pred, sample_weight, update_mask, residual):
        """
        Method for pruning. Loss function can prepare better values for leaves

        :param terminal_regions: indices of terminal regions of each event.
        :param leaf_values: numpy.array, current mapping of leaf indices to prediction values.
        :param X: data (same as passed in fit, ignored usually)
        :param y: labels (same as passed in fit)
        :param y_pred: predictions before adding new tree.
        :param sample_weight: weights od samples (same as passed in fit)
        :param update_mask: which events to use during update?
        :param residual: computed value of negative gradient (before adding tree)
        :return: numpy.array with new prediction values for all leaves.
        """
        return leaf_values


class HessianLossFunction(AbstractLossFunction):
    def __init__(self, regularization=5.):
        """Loss function with diagonal hessian, provides uses Newton-Raphson step to update trees.
        :param regularization: float, penalty for leaves with few events,
            corresponds roughly to the number of added events of both classes to each leaf.
        """
        self.regularization = regularization

    def fit(self, X, y, sample_weight):
        self.regularization_ = self.regularization * numpy.mean(sample_weight)
        return self

    def hessian(self, y_pred):
        raise NotImplementedError('Override this method in loss function.')

    def prepare_tree_params(self, y_pred):
        grad = self.negative_gradient(y_pred)
        hess = self.hessian(y_pred) + 0.01
        return grad / hess, hess

    def prepare_new_leaves_values(self, terminal_regions, leaf_values,
                                  X, y, y_pred, sample_weight, update_mask, residual):
        """ This expression comes from optimization of second-order approximation of loss function."""
        minlength = len(leaf_values)
        nominators = numpy.bincount(terminal_regions, weights=residual, minlength=minlength)
        denominators = numpy.bincount(terminal_regions, weights=self.hessian(y_pred), minlength=minlength)
        return nominators / (denominators + self.regularization_)


class AdaLossFunction(HessianLossFunction):
    """ AdaLossFunction is the same as Exponential Loss Function (exploss)"""
    # TODO write better update
    def fit(self, X, y, sample_weight):
        self.y = y
        self.sample_weight = sample_weight
        self.y_signed = 2 * y - 1
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)
        return self

    def __call__(self, y_pred):
        return numpy.sum(self.sample_weight * numpy.exp(- self.y_signed * y_pred))

    def negative_gradient(self, y_pred):
        return self.y_signed * self.sample_weight * numpy.exp(- self.y_signed * y_pred)

    def hessian(self, y_pred):
        return self.sample_weight * numpy.exp(- self.y_signed * y_pred)


class BinomialDevianceLossFunction(HessianLossFunction):
    """
    Binomial deviance, aka Logistic loss function (logloss).
    """

    def fit(self, X, y, sample_weight):
        self.y = y
        self.sample_weight = sample_weight
        self.y_signed = 2 * y - 1
        self.adjusted_regularization = numpy.mean(sample_weight) * self.regularization
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)
        return self

    def __call__(self, y_pred):
        return numpy.sum(self.sample_weight * numpy.logaddexp(0, - self.y_signed * y_pred))

    def negative_gradient(self, y_pred):
        return self.y_signed * self.sample_weight * expit(- self.y_signed * y_pred)

    def hessian(self, y_pred):
        expits = expit(self.y_signed * y_pred)
        return self.sample_weight * expits * (1 - expits)


class CompositeLossFunction(HessianLossFunction):
    """
    This is exploss for bck and logloss for signal with proper constants.
    Such kind of loss functions is very useful to optimize AMS.
    """

    def fit(self, X, y, sample_weight):
        self.y = y
        self.sample_weight = sample_weight
        self.y_signed = 2 * y - 1
        self.sig_w = (y == 1) * self.sample_weight
        self.bck_w = (y == 0) * self.sample_weight
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)
        return self

    def __call__(self, y_pred):
        result = numpy.sum(self.sig_w * numpy.logaddexp(0, -y_pred))
        result += numpy.sum(self.bck_w * numpy.exp(0.5 * y_pred))
        return result

    def negative_gradient(self, y_pred):
        result = self.sig_w * expit(- y_pred)
        result -= 0.5 * self.bck_w * numpy.exp(0.5 * y_pred)
        return result

    def hessian(self, y_pred):
        expits = expit(- y_pred)
        return self.sig_w * expits * (1 - expits) + self.bck_w * 0.25 * numpy.exp(0.5 * y_pred)


class RankBoostLossFunction(HessianLossFunction):
    def __init__(self, request_column, messup_penalty='square', regularization=0.1, update_terations=1):
        self.update_terations = update_terations
        self.messup_penalty = messup_penalty
        self.request_column = request_column
        HessianLossFunction.__init__(self, regularization=regularization)

    def fit(self, X, y, sample_weight):
        self.queries = X[self.request_column]
        self.y = y
        self.possible_queries, normed_queries = numpy.unique(self.queries, return_inverse=True)
        self.possible_ranks, normed_ranks = numpy.unique(self.y, return_inverse=True)

        self.lookups = [normed_ranks, normed_queries * len(self.possible_ranks) + normed_ranks]
        self.minlengths = [len(self.possible_ranks), len(self.possible_ranks) * len(self.possible_queries)]
        self.rank_penalties = numpy.zeros([len(self.possible_ranks), len(self.possible_ranks)], dtype=float)
        for r1 in self.possible_ranks:
            for r2 in self.possible_ranks:
                if r1 < r2:
                    if self.messup_penalty == 'square':
                        self.rank_penalties[r1, r2] = (r2 - r1) ** 2
                    elif self.messup_penalty == 'linear':
                        self.rank_penalties[r1, r2] = r2 - r1
                    else:
                        raise NotImplementedError()

        self.penalty_matrices = []
        self.penalty_matrices.append(self.rank_penalties / numpy.sqrt(1 + len(y)))
        n_queries = numpy.bincount(normed_queries)
        assert len(n_queries) == len(self.possible_queries)
        self.penalty_matrices.append(sparse.block_diag([self.rank_penalties * 1. / numpy.sqrt(1 + nq) for nq in n_queries]))
        HessianLossFunction.fit(self, X, y, sample_weight=sample_weight)

    def __call__(self, y_pred):
        """
        loss is defined as  w_ij exp(pred_i - pred_j),
        w_ij is zero if label_i <= label_j
        All the other labels are:

        w_ij = (alpha + beta * [query_i = query_j]) rank_penalty_{type_i, type_j}
        rank_penalty_{ij} is zero if i <= j

        :param y_pred: predictions of shape [n_samples]
        :return: value of loss, float
        """
        y_pred -= y_pred.mean()
        pos_exponent = numpy.exp(y_pred)
        neg_exponent = numpy.exp(-y_pred)
        result = 0.
        for lookup, length, penalty_matrix in zip(self.lookups, self.minlengths, self.penalty_matrices):
            pos_stats = numpy.bincount(lookup, weights=pos_exponent)
            neg_stats = numpy.bincount(lookup, weights=neg_exponent)
            result += pos_stats.T.dot(penalty_matrix.dot(neg_stats))
        return result

    def negative_gradient(self, y_pred):
        y_pred -= y_pred.mean()
        pos_exponent = numpy.exp(y_pred)
        neg_exponent = numpy.exp(-y_pred)
        gradient = numpy.zeros(len(y_pred), dtype=float)
        for lookup, length, penalty_matrix in zip(self.lookups, self.minlengths, self.penalty_matrices):
            pos_stats = numpy.bincount(lookup, weights=pos_exponent)
            neg_stats = numpy.bincount(lookup, weights=neg_exponent)
            gradient += pos_exponent * penalty_matrix.dot(neg_stats)[lookup]
            gradient -= neg_exponent * penalty_matrix.T.dot(pos_stats)[lookup]
        return - gradient

    def hessian(self, y_pred):
        y_pred -= y_pred.mean()
        pos_exponent = numpy.exp(y_pred)
        neg_exponent = numpy.exp(-y_pred)
        result = numpy.zeros(len(y_pred), dtype=float)
{'fraction_non_alphanumeric': 0.08800636830001769, 'fraction_numerical': 0.004820449318945693, 'mean_word_length': 2.9553961868112646, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3660398', 'n_tokens_mistral': 6904, 'n_tokens_neox': 6594, 'n_words': 1882}