<commit_before>"""
Ridge regression
"""

# Author: <NAME> <<EMAIL>>
#         <NAME> <<EMAIL>>
#         <NAME> <<EMAIL>>
# License: Simplified BSD


from abc import ABCMeta, abstractmethod
import warnings

import numpy as np
from scipy import linalg
from scipy import sparse
from scipy.sparse import linalg as sp_linalg

from .base import LinearClassifierMixin, LinearModel
from ..base import RegressorMixin
from ..utils.extmath import safe_sparse_dot
from ..utils import safe_asarray
from ..preprocessing import LabelBinarizer
from ..grid_search import GridSearchCV


def ridge_regression(X, y, alpha, sample_weight=1.0, solver='auto',
                     max_iter=None, tol=1e-3):
    """Solve the ridge equation by the method of normal equations.

    Parameters
    ----------
    X : {array-like, sparse matrix, LinearOperator},
        shape = [n_samples, n_features]
        Training data

    y : array-like, shape = [n_samples] or [n_samples, n_targets]
        Target values

    max_iter : int, optional
        Maximum number of iterations for conjugate gradient solver.
        The default value is determined by scipy.sparse.linalg.

    sample_weight : float or numpy array of shape [n_samples]
        Individual weights for each sample

    solver : {'auto', 'dense_cholesky', 'lsqr', 'sparse_cg'}
        Solver to use in the computational routines:

        - 'auto' chooses the solver automatically based on the type of data.

        - 'dense_cholesky' uses the standard scipy.linalg.solve function to
          obtain a closed-form solution.

        - 'sparse_cg' uses the conjugate gradient solver as found in
          scipy.sparse.linalg.cg. As an iterative algorithm, this solver is
          more appropriate than 'dense_cholesky' for large-scale data
          (possibility to set `tol` and `max_iter`).

        - 'lsqr' uses the dedicated regularized least-squares routine
          scipy.sparse.linalg.lsqr. It is the fatest but may not be available
          in old scipy versions. It also uses an iterative procedure.

        All three solvers support both dense and sparse data.

    tol: float
        Precision of the solution.

    Returns
    -------
    coef: array, shape = [n_features] or [n_targets, n_features]
        Weight vector(s).

    Notes
    -----
    This function won't compute the intercept.
    """

    n_samples, n_features = X.shape
    has_sw = isinstance(sample_weight, np.ndarray) or sample_weight != 1.0

    if solver == 'auto':
        # cholesky if it's a dense array and cg in
        # any other case
        if hasattr(X, '__array__'):
            solver = 'dense_cholesky'
        else:
            solver = 'sparse_cg'

    elif solver == 'lsqr' and not hasattr(sp_linalg, 'lsqr'):
        warnings.warn("""lsqr not available on this machine, falling back
                      to sparse_cg.""")
        solver = 'sparse_cg'

    if has_sw:
        solver = 'dense_cholesky'

    if solver == 'sparse_cg':
        # gradient descent
        X1 = sp_linalg.aslinearoperator(X)
        if y.ndim == 1:
            y1 = np.reshape(y, (-1, 1))
        else:
            y1 = y
        coefs = np.empty((y1.shape[1], n_features))

        if n_features > n_samples:
            def mv(x):
                return X1.matvec(X1.rmatvec(x)) + alpha * x
        else:
            def mv(x):
                return X1.rmatvec(X1.matvec(x)) + alpha * x

        for i in range(y1.shape[1]):
            y_column = y1[:, i]
            if n_features > n_samples:
                # kernel ridge
                # w = X.T * inv(X X^t + alpha*Id) y
                C = sp_linalg.LinearOperator(
                    (n_samples, n_samples), matvec=mv, dtype=X.dtype)
                coef, info = sp_linalg.cg(C, y_column, tol=tol)
                coefs[i] = X1.rmatvec(coef)
            else:
                # ridge
                # w = inv(X^t X + alpha*Id) * X.T y
                y_column = X1.rmatvec(y_column)
                C = sp_linalg.LinearOperator(
                    (n_features, n_features), matvec=mv, dtype=X.dtype)
                coefs[i], info = sp_linalg.cg(C, y_column, maxiter=max_iter,
                                              tol=tol)
            if info != 0:
                raise ValueError("Failed with error code %d" % info)

        if y.ndim == 1:
            coefs = np.ravel(coefs)

        return coefs
    elif solver == "lsqr":
        if y.ndim == 1:
            y1 = np.reshape(y, (-1, 1))
        else:
            y1 = y
        coefs = np.empty((y1.shape[1], n_features))

        # According to the lsqr documentation, alpha = damp^2.
        sqrt_alpha = np.sqrt(alpha)

        for i in range(y1.shape[1]):
            y_column = y1[:, i]
            coefs[i] = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha,
                                      atol=tol, btol=tol, iter_lim=max_iter)[0]

        if y.ndim == 1:
            coefs = np.ravel(coefs)

        return coefs
    else:
        # normal equations (cholesky) method
        if n_features > n_samples or has_sw:
            # kernel ridge
            # w = X.T * inv(X X^t + alpha*Id) y
            K = safe_sparse_dot(X, X.T, dense_output=True)
            if has_sw:
                sw = np.sqrt(sample_weight)
                if y.ndim == 1:
                    y = y * sw
                else:
                    # Deal with multiple-output problems
                    y = y * sw[:, np.newaxis]
                K *= np.outer(sw, sw)
            K.flat[::n_samples + 1] += alpha
            dual_coef = linalg.solve(K, y,
                                     sym_pos=True, overwrite_a=True)
            coef = safe_sparse_dot(X.T, dual_coef, dense_output=True)
        else:
            # ridge
            # w = inv(X^t X + alpha*Id) * X.T y
            A = safe_sparse_dot(X.T, X, dense_output=True)
            A.flat[::n_features + 1] += alpha
            Xy = safe_sparse_dot(X.T, y, dense_output=True)
            coef = linalg.solve(A, Xy, sym_pos=True, overwrite_a=True)

        return coef.T
<commit_msg>BUG: fix sample weights in ridge<commit_after>"""
Ridge regression
"""

# Author: <NAME> <<EMAIL>>
#         <NAME> <<EMAIL>>
#         <NAME> <<EMAIL>>
# License: Simplified BSD


from abc import ABCMeta, abstractmethod
import warnings

import numpy as np
from scipy import linalg
from scipy import sparse
from scipy.sparse import linalg as sp_linalg

from .base import LinearClassifierMixin, LinearModel
from ..base import RegressorMixin
from ..utils.extmath import safe_sparse_dot
from ..utils import safe_asarray
from ..preprocessing import LabelBinarizer
from ..grid_search import GridSearchCV


def ridge_regression(X, y, alpha, sample_weight=1.0, solver='auto',
                     max_iter=None, tol=1e-3):
    """Solve the ridge equation by the method of normal equations.

    Parameters
    ----------
    X : {array-like, sparse matrix, LinearOperator},
        shape = [n_samples, n_features]
        Training data

    y : array-like, shape = [n_samples] or [n_samples, n_targets]
        Target values

    max_iter : int, optional
        Maximum number of iterations for conjugate gradient solver.
        The default value is determined by scipy.sparse.linalg.

    sample_weight : float or numpy array of shape [n_samples]
        Individual weights for each sample

    solver : {'auto', 'dense_cholesky', 'lsqr', 'sparse_cg'}
        Solver to use in the computational routines:

        - 'auto' chooses the solver automatically based on the type of data.

        - 'dense_cholesky' uses the standard scipy.linalg.solve function to
          obtain a closed-form solution.

        - 'sparse_cg' uses the conjugate gradient solver as found in
          scipy.sparse.linalg.cg. As an iterative algorithm, this solver is
          more appropriate than 'dense_cholesky' for large-scale data
          (possibility to set `tol` and `max_iter`).

        - 'lsqr' uses the dedicated regularized least-squares routine
          scipy.sparse.linalg.lsqr. It is the fatest but may not be available
          in old scipy versions. It also uses an iterative procedure.

        All three solvers support both dense and sparse data.

    tol: float
        Precision of the solution.

    Returns
    -------
    coef: array, shape = [n_features] or [n_targets, n_features]
        Weight vector(s).

    Notes
    -----
    This function won't compute the intercept.
    """

    n_samples, n_features = X.shape
    has_sw = isinstance(sample_weight, np.ndarray) or sample_weight != 1.0

    if solver == 'auto':
        # cholesky if it's a dense array and cg in
        # any other case
        if hasattr(X, '__array__'):
            solver = 'dense_cholesky'
        else:
            solver = 'sparse_cg'

    elif solver == 'lsqr' and not hasattr(sp_linalg, 'lsqr'):
        warnings.warn("""lsqr not available on this machine, falling back
                      to sparse_cg.""")
        solver = 'sparse_cg'

    if has_sw:
        solver = 'dense_cholesky'

    if solver == 'sparse_cg':
        # gradient descent
        X1 = sp_linalg.aslinearoperator(X)
        if y.ndim == 1:
            y1 = np.reshape(y, (-1, 1))
        else:
            y1 = y
        coefs = np.empty((y1.shape[1], n_features))

        if n_features > n_samples:
            def mv(x):
                return X1.matvec(X1.rmatvec(x)) + alpha * x
        else:
            def mv(x):
                return X1.rmatvec(X1.matvec(x)) + alpha * x

        for i in range(y1.shape[1]):
            y_column = y1[:, i]
            if n_features > n_samples:
                # kernel ridge
                # w = X.T * inv(X X^t + alpha*Id) y
                C = sp_linalg.LinearOperator(
                    (n_samples, n_samples), matvec=mv, dtype=X.dtype)
                coef, info = sp_linalg.cg(C, y_column, tol=tol)
                coefs[i] = X1.rmatvec(coef)
            else:
                # ridge
                # w = inv(X^t X + alpha*Id) * X.T y
                y_column = X1.rmatvec(y_column)
                C = sp_linalg.LinearOperator(
                    (n_features, n_features), matvec=mv, dtype=X.dtype)
                coefs[i], info = sp_linalg.cg(C, y_column, maxiter=max_iter,
                                              tol=tol)
            if info != 0:
                raise ValueError("Failed with error code %d" % info)

        if y.ndim == 1:
            coefs = np.ravel(coefs)

        return coefs
    elif solver == "lsqr":
        if y.ndim == 1:
            y1 = np.reshape(y, (-1, 1))
        else:
            y1 = y
        coefs = np.empty((y1.shape[1], n_features))

        # According to the lsqr documentation, alpha = damp^2.
        sqrt_alpha = np.sqrt(alpha)

        for i in range(y1.shape[1]):
            y_column = y1[:, i]
            coefs[i] = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha,
                                      atol=tol, btol=tol, iter_lim=max_iter)[0]

        if y.ndim == 1:
            coefs = np.ravel(coefs)

        return coefs
    else:
        # normal equations (cholesky) method
        if n_features > n_samples or has_sw:
            # kernel ridge
            # w = X.T * inv(X X^t + alpha*Id) y
            K = safe_sparse_dot(X, X.T, dense_output=True)
            if has_sw:
                # We are doing a little danse with the sample weights to
                # avoid copying the original X, which could be big
                sw = np.sqrt(sample_weight)
                if y.ndim == 1:
                    y = y * sw
                else:
                    # Deal with multiple-output problems
                    y = y * sw[:, np.newaxis]
                K *= np.outer(sw, sw)
            K.flat[::n_samples + 1] += alpha
            dual_coef = linalg.solve(K, y,
                                     sym_pos=True, overwrite_a=True)
            if has_sw:
                if dual_coef.ndim == 1:
                    dual_coef *= sw
                else:
                    # Deal with multiple-output problems
                    dual_coef *= sw[:, np.newaxis]
            coef = safe_sparse_dot(X.T, dual_coef, dense_output=True)
        else:
            # ridge
            # w = inv(X^t X + alpha*Id) * X.T y
            A = safe_sparse_dot(X.T, X, dense_output=True)
            A.flat[::n_features + 1] += alpha
            Xy = safe_sparse_dot(X.T, y, dense_output=True)
            coef = linalg.solve(A, Xy, sym_pos=True, overwrite_a=True)

        return coef.T

{'fraction_non_alphanumeric': 0.09879000159210317, 'fraction_numerical': 0.006607228148384015, 'mean_word_length': 1.7263454861111112, 'pattern_counts': {'":': 2, '<': 21, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '750752', 'n_tokens_mistral': 4134, 'n_tokens_neox': 3860, 'n_words': 1128}