def mean_absolute_error(y_true, y_pred):
    """
    Mean absolute error and its standard deviation.
    
    If you need only mean absolute error, use 
    :func:`sklearn.metrics.mean_absolute_error`
    
    Parameters
    ----------
    y_true : array, shape(n_samples,)
        Ground truth scores
    y_pred : array, shape(n_samples,)
        Predicted scores

    Returns
    -------
    mean : float
        mean of squared errors
    stdev : float
        standard deviation of squared errors
    """

    # check inputs
    assert_all_finite(y_true)
    y_true = as_float_array(y_true)
    assert_all_finite(y_pred)
    y_pred = as_float_array(y_pred)
    check_consistent_length(y_true, y_pred)

    # calculate errors
    errs = np.abs(y_true - y_pred)
    mean = np.nanmean(errs)
    stdev = np.nanstd(errs)

    return mean, stdev

def test_np_matrix():
    """
    Confirm that input validation code does not return np.matrix
    """
    X = np.arange(12).reshape(3, 4)

    assert_false(isinstance(as_float_array(X), np.matrix))
    assert_false(isinstance(as_float_array(np.matrix(X)), np.matrix))
    assert_false(isinstance(as_float_array(sp.csc_matrix(X)), np.matrix))

    assert_false(isinstance(atleast2d_or_csr(X), np.matrix))
    assert_false(isinstance(atleast2d_or_csr(np.matrix(X)), np.matrix))
    assert_false(isinstance(atleast2d_or_csr(sp.csc_matrix(X)), np.matrix))

    assert_false(isinstance(atleast2d_or_csc(X), np.matrix))
    assert_false(isinstance(atleast2d_or_csc(np.matrix(X)), np.matrix))
    assert_false(isinstance(atleast2d_or_csc(sp.csr_matrix(X)), np.matrix))

    assert_false(isinstance(safe_asarray(X), np.matrix))
    assert_false(isinstance(safe_asarray(np.matrix(X)), np.matrix))
    assert_false(isinstance(safe_asarray(sp.lil_matrix(X)), np.matrix))

    assert_true(atleast2d_or_csr(X, copy=False) is X)
    assert_false(atleast2d_or_csr(X, copy=True) is X)
    assert_true(atleast2d_or_csc(X, copy=False) is X)
    assert_false(atleast2d_or_csc(X, copy=True) is X)

def test_np_matrix():
    # Confirm that input validation code does not return np.matrix
    X = np.arange(12).reshape(3, 4)

    assert_false(isinstance(as_float_array(X), np.matrix))
    assert_false(isinstance(as_float_array(np.matrix(X)), np.matrix))
    assert_false(isinstance(as_float_array(sp.csc_matrix(X)), np.matrix))

def transform_with_scaler( Y, scaler=None, wrt_X=[] ):
    Y = as_float_array( Y )
    if len(wrt_X) and not scaler:
        wrt_X = as_float_array( wrt_X )
        scaler = get_scaler( wrt_X )
    with_mean = scaler.mean_                                      
    with_stdv = scaler.std_                                       
    Z = scaler.transform(Y)                               
    return Z

def item_finder_statistics(y_true, y_pred):
    """
    Full Statistics of prediction performance

    * n_samples
    * mean_absolute_error: mean, stdev
    * mean_squared_error: mean, rmse, stdev
    * predicted: mean, stdev
    * true: mean, stdev

    Parameters
    ----------
    y_true : array, shape=(n_samples,)
        Ground truth scores
    y_pred : array, shape=(n_samples,)
        Predicted scores

    Returns
    -------
    stats : dict
        Full statistics of prediction performance
    """

    # check inputs
    assert_all_finite(y_true)
    if not is_binary_score(y_true):
        raise ValueError('True scores must be binary')
    y_true = as_float_array(y_true)
    assert_all_finite(y_pred)
    y_pred = as_float_array(y_pred)
    check_consistent_length(y_true, y_pred)

    # calc statistics
    stats = {}

    # dataset size
    stats['n_samples'] = y_true.size

    # descriptive statistics of ground truth scores
    stats['true'] = {'mean': np.mean(y_true), 'stdev': np.std(y_true)}

    # descriptive statistics of ground predicted scores
    stats['predicted'] = {'mean': np.mean(y_pred), 'stdev': np.std(y_pred)}

    # statistics at least 0 and 1 must be contained in a score array
    if is_binary_score(y_true, allow_uniform=False):

        # AUC (area under the curve)
        stats['area under the curve'] = skm.roc_auc_score(y_true, y_pred)

    return stats

def item_finder_report(y_true, y_pred, disp=True):
    """
    Report brief summary of prediction performance

    * AUC
    * number of data
    * mean and standard dev. of true scores
    * mean and standard dev. of predicted scores

    Parameters
    ----------
    y_true : array, shape(n_samples,)
        Ground truth scores
    y_pred : array, shape(n_samples,)
        Predicted scores
    disp : bool, optional, default=True
        if True, print report

    Returns
    -------
    stats : dict
        belief summary of prediction performance
    """

    # check inputs
    assert_all_finite(y_true)
    if not is_binary_score(y_true):
        raise ValueError('True scores must be binary')
    y_true = as_float_array(y_true)
    assert_all_finite(y_pred)
    y_pred = as_float_array(y_pred)
    check_consistent_length(y_true, y_pred)

    # calc statistics
    stats = {
        'n_samples': y_true.size,
        'true': {'mean': np.mean(y_true), 'stdev': np.std(y_true)},
        'predicted': {'mean': np.mean(y_pred), 'stdev': np.std(y_pred)}}

    # statistics at least 0 and 1 must be contained in a score array
    if is_binary_score(y_true, allow_uniform=False):
        stats['area under the curve'] = skm.roc_auc_score(y_true, y_pred)

    # display statistics
    if disp:
        print(
            json.dumps(
                stats, sort_keys=True, indent=4, separators=(',', ': '),
                ensure_ascii=False), file=sys.stderr)

    return stats

    def fit(self, X, y):

        n_samples, self.n_features = X.shape
        self.n_outputs = y.shape[1]
        self._init_fit(X)

        self.hidden_activations_ = self._get_hidden_activations(X)

        if self.regularized:
            self._solve_regularized(as_float_array(y, copy=True))
        else:
            self._solve(as_float_array(y, copy=True))

        return self

def score_predictor_report(y_true, y_pred, disp=True):
    """
    Report brief summary of prediction performance
    
    * mean absolute error
    * root mean squared error
    * number of data
    * mean and standard dev. of true scores
    * mean and standard dev. of predicted scores

    Parameters
    ----------
    y_true : array, shape(n_samples,)
        Ground truth scores
    y_pred : array, shape(n_samples,)
        Predicted scores
    disp : bool, optional, default=True
        if True, print report

    Returns
    -------
    stats : dict
        belief summary of prediction performance
    """

    # check inputs
    assert_all_finite(y_true)
    y_true = as_float_array(y_true)
    assert_all_finite(y_pred)
    y_pred = as_float_array(y_pred)
    check_consistent_length(y_true, y_pred)

    # calc statistics
    stats = {
        'mean absolute error': skm.mean_absolute_error(y_true, y_pred),
        'root mean squared error':
            np.sqrt(np.maximum(skm.mean_squared_error(y_true, y_pred), 0.)),
        'n_samples': y_true.size,
        'true': {'mean': np.mean(y_true), 'stdev': np.std(y_true)},
        'predicted': {'mean': np.mean(y_pred), 'stdev': np.std(y_pred)}}

    # display statistics
    if disp:
        print(json.dumps(
            stats, sort_keys=True, indent=4, separators=(',', ': '),
            ensure_ascii=False),
            file=sys.stderr)

    return stats

def k_modes(X, n_clusters, n_init=1, max_iter=5,
            verbose=False, tol=1e-4, random_state=None, copy_x=True, n_jobs=1):
    """K-modes clustering algorithm."""
    if n_init <= 0:
        raise ValueError("Invalid number of initializations."
                         " n_init=%d must be bigger than zero." % n_init)

    X = as_float_array(X, copy=copy_x)
    matrix_all_irm = _compute_all_irm(X, n_clusters)
    best_labels, best_modes, best_mirm = None, None, -np.inf

    if n_jobs == 1:

        for j in range(2,n_clusters+1):
            # For a single thread, less memory is needed if we just store one set
            # of the best results (as opposed to one set per run per thread).
            for it in range(n_init):
                # run a k-modes once
                labels, modes, mirm_sum = _kmodes_single(
                    X, j, matrix_all_irm, max_iter=max_iter,
                    verbose=verbose, tol=tol, random_state=random_state)
                # determine if these results are the best so far
                if mirm_sum >= best_mirm:
                    best_labels = labels.copy()
                    best_modes = modes.copy()
                    best_mirm = mirm_sum
    else:
        # TODO:
        pass

    return best_modes, best_labels, best_mirm

    def fit(self, X, y):
        """
        Fit the model using X, y as training data.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape [n_samples, n_features]
            Training vectors, where n_samples is the number of samples
            and n_features is the number of features.

        y : array-like of shape [n_samples, n_outputs]
            Target values (class labels in classification, real numbers in
            regression)

        Returns
        -------
        self : object

            Returns an instance of self.
        """
        # fit random hidden layer and compute the hidden layer activations
        self.hidden_activations_ = self.hidden_layer.fit_transform(X)

        # solve the regression from hidden activations to outputs
        self._fit_regression(as_float_array(y, copy=True))

        return self

    def fit(self, X, y=None, **params):
        """Fit the model with X.

        Parameters
        ----------
        X: array-like, shape (n_samples, n_features)
            Training data, where n_samples in the number of samples
            and n_features is the number of features.

        Returns
        -------
        self : object
            Returns the instance itself.
            
        Notes
        -----
        Calling multiple times will update the components
        """

        X = array2d(X)
        n_samples, n_features = X.shape
        X = as_float_array(X, copy=self.copy)

        if self.iteration != 0 and n_features != self.components_.shape[1]:
            raise ValueError("The dimensionality of the new data and the existing components_ does not match")

        # incrementally fit the model
        for i in range(0, X.shape[0]):
            self.partial_fit(X[i, :])

        return self

 def transform(self, X):
     X = as_float_array(X, copy=self.copy) 
     if self.mean_ is not None and self.std_ is not None:
         X -= self.mean_
         X /= self.std_
     X_whitend = np.dot(X, self.components_)
     return X_whitend

    def _fit(self, X):
        """Fit the model to the data X.

        Parameters
        ----------
        X: array-like, shape (n_samples, n_features)
            Training vector, where n_samples in the number of samples and
            n_features is the number of features.

        Returns
        -------
        X : ndarray, shape (n_samples, n_features)
            The input data, copied, centered and whitened when requested.
        """
        random_state = self._random_state
        X = np.atleast_2d(as_float_array(X))
        self._initialize(X, random_state)
        for it in range(self.n_iter):
            if it % 10 == 0:
                self._print_status(it)
            else:
                logger.info("<{}>".format(it))
            self._sample_topics(random_state)
        self._print_status(self.n_iter)
        self.components_ = self.nzw + self.eta
        self.components_ /= np.sum(self.components_, axis=1, keepdims=True)
        return self.ndz

 def fit(self, X):
     n_samples, self.n_features = X.shape
     self.n_outputs = X.shape[1]
     self._init_fit(X)
     self.hidden_activations_ = self._get_hidden_activations(X)
     self._regularized(as_float_array(X, copy=True))
     #self.coef_output_ = safe_sparse_dot(pinv2(self.hidden_activations_), X)
     return self

 def fit(self, X):
     X = as_float_array(X, copy=self.copy)
     self._mean = np.mean(X, axis=0)
     X -= self._mean
     sigma = np.dot(X.T,X) / X.shape[1]
     U, S, V = linalg.svd(sigma)
     tmp = np.dot(U, np.diag(1 / np.sqrt(S + self.regularization)))
     self._components = np.dot(tmp, U.T)
     return self

 def fit(self, X, y=None):
     X = array2d(X)
     X = as_float_array(X, copy = self.copy)
     print X.shape
     sigma = np.dot(X.T,X) / X.shape[1]
     U, S, V = linalg.svd(sigma)
     tmp = np.dot(U, np.diag(1/np.sqrt(S+self.regularization)))
     self.components_ = np.dot(tmp, U.T)
     return self

def import_lda(filename):
    lines = filename.read().splitlines()
    nr_lines = len(lines)
    labels = range(nr_lines)
    features = range(nr_lines)
    for i in range(nr_lines):
        labels[i], data = lines[i].split(',')
        features[i] = map(float, data.split())
    X = as_float_array(features)
    return labels, X

def test_as_float_array():
    """
    Test function for as_float_array
    """
    X = np.ones((3, 10), dtype=np.int32)
    X = X + np.arange(10, dtype=np.int32)
    # Checks that the return type is ok
    X2 = as_float_array(X, copy=False)
    np.testing.assert_equal(X2.dtype, np.float32)
    # Another test
    X = X.astype(np.int64)
    X2 = as_float_array(X, copy=True)
    # Checking that the array wasn't overwritten
    assert_true(as_float_array(X, False) is not X)
    # Checking that the new type is ok
    np.testing.assert_equal(X2.dtype, np.float64)
    # Here, X is of the right type, it shouldn't be modified
    X = np.ones((3, 2), dtype=np.float32)
    assert_true(as_float_array(X, copy=False) is X)

    def transform(self, X, y=None, copy=None):
        """Perform ZCA whitening

        Parameters
        ----------
        X : array-like with shape [n_samples, n_features]
            The data to whiten along the features axis.
        """
        check_is_fitted(self, 'mean_')
        X = as_float_array(X, copy=self.copy)
        return np.dot(X - self.mean_, self.whiten_.T)

def get_distance_metric(X, dtype='euclidean', normalize=True, nrows=10 ):
    X = as_float_array( X )
    dist = nn.DistanceMetric.get_metric('euclidean')
    d = dist.pairwise(X)
    print '-' * 80
    print introspection( dist )
    nrows=min(nrows,len(d))
    for row in d[:nrows]:
        print [ "%.2f" % x for x in row[:nrows] ]
    print '-' * 80
    return dist