"""
Robust location and covariance estimators.

Here are implemented estimators that are resistant to outliers.

"""
# Author: Virgile Fritsch <virgile.fritsch@inria.fr>
#
# License: BSD 3 clause
import warnings
import numbers
import numpy as np
from scipy import linalg
from scipy.stats import chi2

from . import empirical_covariance, EmpiricalCovariance
from ..utils.extmath import fast_logdet, pinvh
from ..utils import check_random_state, check_array


# Minimum Covariance Determinant
#   Implementing of an algorithm by Rousseeuw & Van Driessen described in
#   (A Fast Algorithm for the Minimum Covariance Determinant Estimator,
#   1999, American Statistical Association and the American Society
#   for Quality, TECHNOMETRICS)
# XXX Is this really a public function? It's not listed in the docs or
# exported by sklearn.covariance. Deprecate?
def c_step(X, n_support, remaining_iterations=30, initial_estimates=None,
           verbose=False, cov_computation_method=empirical_covariance,
           random_state=None):
    """C_step procedure described in [Rouseeuw1984]_ aiming at computing MCD.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        Data set in which we look for the n_support observations whose
        scatter matrix has minimum determinant.

    n_support : int, > n_samples / 2
        Number of observations to compute the robust estimates of location
        and covariance from.

    remaining_iterations : int, optional
        Number of iterations to perform.
        According to [Rouseeuw1999]_, two iterations are sufficient to get
        close to the minimum, and we never need more than 30 to reach
        convergence.

    initial_estimates : 2-tuple, optional
        Initial estimates of location and shape from which to run the c_step
        procedure:
        - initial_estimates[0]: an initial location estimate
        - initial_estimates[1]: an initial covariance estimate

    verbose : boolean, optional
        Verbose mode.

    random_state : integer or numpy.RandomState, optional
        The random generator used. If an integer is given, it fixes the
        seed. Defaults to the global numpy random number generator.

    cov_computation_method : callable, default empirical_covariance
        The function which will be used to compute the covariance.
        Must return shape (n_features, n_features)

    Returns
    -------
    location : array-like, shape (n_features,)
        Robust location estimates.

    covariance : array-like, shape (n_features, n_features)
        Robust covariance estimates.

    support : array-like, shape (n_samples,)
        A mask for the `n_support` observations whose scatter matrix has
        minimum determinant.

    References
    ----------
    .. [Rouseeuw1999] A Fast Algorithm for the Minimum Covariance Determinant
        Estimator, 1999, American Statistical Association and the American
        Society for Quality, TECHNOMETRICS

    """
    X = np.asarray(X)
    random_state = check_random_state(random_state)
    return _c_step(X, n_support, remaining_iterations=remaining_iterations,
                   initial_estimates=initial_estimates, verbose=verbose,
                   cov_computation_method=cov_computation_method,
                   random_state=random_state)


def _c_step(X, n_support, random_state, remaining_iterations=30,
            initial_estimates=None, verbose=False,
            cov_computation_method=empirical_covariance):
    n_samples, n_features = X.shape

    # Initialisation
    support = np.zeros(n_samples, dtype=bool)
    if initial_estimates is None:
        # compute initial robust estimates from a random subset
        support[random_state.permutation(n_samples)[:n_support]] = True
    else:
        # get initial robust estimates from the function parameters
        location = initial_estimates[0]
        covariance = initial_estimates[1]
        # run a special iteration for that case (to get an initial support)
        precision = pinvh(covariance)
        X_centered = X - location
        dist = (np.dot(X_centered, precision) * X_centered).sum(1)
        # compute new estimates
        support[np.argsort(dist)[:n_support]] = True

    X_support = X[support]
    location = X_support.mean(0)
    covariance = cov_computation_method(X_support)

    # Iterative procedure for Minimum Covariance Determinant computation
    det = fast_logdet(covariance)
    previous_det = np.inf
    while (det < previous_det) and (remaining_iterations > 0):
        # save old estimates values
        previous_location = location
        previous_covariance = covariance
        previous_det = det
        previous_support = support
        # compute a new support from the full data set mahalanobis distances
        precision = pinvh(covariance)
        X_centered = X - location
        dist = (np.dot(X_centered, precision) * X_centered).sum(axis=1)
        # compute new estimates
        support = np.zeros(n_samples, dtype=bool)
        support[np.argsort(dist)[:n_support]] = True
        X_support = X[support]
        location = X_support.mean(axis=0)
        covariance = cov_computation_method(X_support)
        det = fast_logdet(covariance)
        # update remaining iterations for early stopping
        remaining_iterations -= 1

    previous_dist = dist
    dist = (np.dot(X - location, precision) * (X - location)).sum(axis=1)
    # Catch computation errors
    if np.isinf(det):
        raise ValueError(
            "Singular covariance matrix. "
            "Please check that the covariance matrix corresponding "
            "to the dataset is full rank and that MinCovDet is used with "
            "Gaussian-distributed data (or at least data drawn from a "
            "unimodal, symmetric distribution.")
    # Check convergence
    if np.allclose(det, previous_det):
        # c_step procedure converged
        if verbose:
            print("Optimal couple (location, covariance) found before"
                  " ending iterations (%d left)" % (remaining_iterations))
        results = location, covariance, det, support, dist
    elif det > previous_det:
        # determinant has increased (should not happen)
        warnings.warn("Warning! det > previous_det (%.15f > %.15f)"
                      % (det, previous_det), RuntimeWarning)
        results = previous_location, previous_covariance, \
            previous_det, previous_support, previous_dist

    # Check early stopping
    if remaining_iterations == 0:
        if verbose:
            print('Maximum number of iterations reached')
        results = location, covariance, det, support, dist

    return results


def select_candidates(X, n_support, n_trials, select=1, n_iter=30,
                      verbose=False,
                      cov_computation_method=empirical_covariance,
                      random_state=None):
    """Finds the best pure subset of observations to compute MCD from it.

    The purpose of this function is to find the best sets of n_support
    observations with respect to a minimization of their covariance
    matrix determinant. Equivalently, it removes n_samples-n_support
    observations to construct what we call a pure data set (i.e. not
    containing outliers). The list of the observations of the pure
    data set is referred to as the `support`.

    Starting from a random support, the pure data set is found by the
    c_step procedure introduced by Rousseeuw and Van Driessen in
    [Rouseeuw1999]_.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        Data (sub)set in which we look for the n_support purest observations.

    n_support : int, [(n + p + 1)/2] < n_support < n
        The number of samples the pure data set must contain.

    select : int, int > 0
        Number of best candidates results to return.

    n_trials : int, nb_trials > 0 or 2-tuple
        Number of different initial sets of observations from which to
        run the algorithm.
        Instead of giving a number of trials to perform, one can provide a
        list of initial estimates that will be used to iteratively run
        c_step procedures. In this case:
        - n_trials[0]: array-like, shape (n_trials, n_features)
          is the list of `n_trials` initial location estimates
        - n_trials[1]: array-like, shape (n_trials, n_features, n_features)
          is the list of `n_trials` initial covariances estimates

    n_iter : int, nb_iter > 0
        Maximum number of iterations for the c_step procedure.
        (2 is enough to be close to the final solution. "Never" exceeds 20).

    random_state : integer or numpy.RandomState, default None
        The random generator used. If an integer is given, it fixes the
        seed. Defaults to the global numpy random number generator.

    cov_computation_method : callable, default empirical_covariance
        The function which will be used to compute the covariance.
        Must return shape (n_features, n_features)

    verbose : boolean, default False
        Control the output verbosity.

    See Also
    ---------
    c_step

    Returns
    -------
    best_locations : array-like, shape (select, n_features)
        The `select` location estimates computed from the `select` best
        supports found in the data set (`X`).

    best_covariances : array-like, shape (select, n_features, n_features)
        The `select` covariance estimates computed from the `select`
        best supports found in the data set (`X`).

    best_supports : array-like, shape (select, n_samples)
        The `select` best supports found in the data set (`X`).

    References
    ----------
    .. [Rouseeuw1999] A Fast Algorithm for the Minimum Covariance Determinant
        Estimator, 1999, American Statistical Association and the American
        Society for Quality, TECHNOMETRICS

    """
    random_state = check_random_state(random_state)
    n_samples, n_features = X.shape

    if isinstance(n_trials, numbers.Integral):
        run_from_estimates = False
    elif isinstance(n_trials, tuple):
        run_from_estimates = True
        estimates_list = n_trials
        n_trials = estimates_list[0].shape[0]
    else:
        raise TypeError("Invalid 'n_trials' parameter, expected tuple or "
                        " integer, got %s (%s)" % (n_trials, type(n_trials)))

    # compute `n_trials` location and shape estimates candidates in the subset
    all_estimates = []
    if not run_from_estimates:
        # perform `n_trials` computations from random initial supports
        for j in range(n_trials):
            all_estimates.append(
                _c_step(
                    X, n_support, remaining_iterations=n_iter, verbose=verbose,
                    cov_computation_method=cov_computation_method,
                    random_state=random_state))
    else:
        # perform computations from every given initial estimates
        for j in range(n_trials):
            initial_estimates = (estimates_list[0][j], estimates_list[1][j])
            all_estimates.append(_c_step(
                X, n_support, remaining_iterations=n_iter,
                initial_estimates=initial_estimates, verbose=verbose,
                cov_computation_method=cov_computation_method,
                random_state=random_state))
    all_locs_sub, all_covs_sub, all_dets_sub, all_supports_sub, all_ds_sub = \
        zip(*all_estimates)
    # find the `n_best` best results among the `n_trials` ones
    index_best = np.argsort(all_dets_sub)[:select]