def _pairwise_callable(X, Y, metric, **kwds):
    """Handle the callable case for pairwise_{distances,kernels}
    """
    try:
        X, Y = check_pairwise_arrays(X, Y)
    except ValueError:
        X, Y = check_pairwise_arrays(X, Y, dtype=object)  # try not to convert

    if X is Y:
        # Only calculate metric for upper triangle
        out = np.zeros((X.shape[0], Y.shape[0]), dtype='float')
        iterator = itertools.combinations(range(X.shape[0]), 2)
        for i, j in iterator:
            out[i, j] = metric(X[i], Y[j], **kwds)

        # Make symmetric
        # NB: out += out.T will produce incorrect results
        out = out + out.T

        # Calculate diagonal
        # NB: nonzero diagonals are allowed for both metrics and kernels
        for i in range(X.shape[0]):
            x = X[i]
            out[i, i] = metric(x, x, **kwds)

    else:
        # Calculate all cells
        out = np.empty((X.shape[0], Y.shape[0]), dtype='float')
        iterator = itertools.product(range(X.shape[0]), range(Y.shape[0]))
        for i, j in iterator:
            out[i, j] = metric(X[i], Y[j], **kwds)

    return out

def MaternKernel(X, Y=None, gamma = None, p = 0):
    """Compute the generalized normal kernel between X and Y.
    The generalized normal kernel is defined as::
        K(x, y) = exp(-gamma ||x-y||_1^beta)
    for each pair of rows x in X and y in Y.
    Parameters
    ----------
    X : array of shape (n_samples_X, n_features)
    Y : array of shape (n_samples_Y, n_features)
    gamma : float
    Returns
    -------
    kernel_matrix : array of shape (n_samples_X, n_samples_Y)
    """
    assert(p == int(p))

    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]

    r = manhattan_distances(X, Y)
    if p == 0:
        K = -gamma * r
        np.exp(K, K)    # exponentiate K in-place
    if p == 1:
        K = -gamma * r * math.sqrt(3)
        np.exp(K, K)    # exponentiate K in-place
        K *= (1+gamma * r * math.sqrt(3))
    if p == 1:
        K = -gamma * r * math.sqrt(5)
        np.exp(K, K)    # exponentiate K in-place
        K *= (1+gamma * r * math.sqrt(5) + 5./3. * (r*gamma)**2)
    return K

def poisson_kernel(X, Y=None, gamma=None, Sigma_inv = None):
    """
    Compute the poisson kernel between X and Y::
        K(x, y) = exp(-gamma ||x-mu||^2/mu)
        mu = centroid of X (=X if X.shape[0] == 1)
    for each pair of rows x in X and y in Y.
    Parameters
    ----------
    X : array of shape (n_samples_X, n_features)
    Y : array of shape (n_samples_Y, n_features)
    gamma : float
    Returns
    -------
    kernel_matrix : array of shape (n_samples_X, n_samples_Y)
    """
    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]
    if Sigma_inv is None:
        raise ValueError('Missing Sigma_inv')
    
    v = X - Y
    K = -0.5 * gamma * np.sqrt(v.dot(Sigma_inv).dot(v.T))
    np.exp(K, K)    # exponentiate K in-place
    return K

def trimmedrbf_kernel(X, Y=None, gamma=None, robust_gamma = None):
    """
    Compute the rbf (gaussian) kernel between X and Y::

        K(x, y) = exp(-gamma ||x-y||**2)

    for each pair of rows x in X and y in Y.

    Parameters
    ----------
    X : array of shape (n_samples_X, n_features)

    Y : array of shape (n_samples_Y, n_features)

    gamma : float

    Returns
    -------
    kernel_matrix : array of shape (n_samples_X, n_samples_Y)
    """
    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]

    K = euclidean_distances(X, Y, squared=True)
    print K
    print "SHape kernel" + str(np.where(np.sqrt(K) > robust_gamma)[0].shape)
    K[np.where(np.sqrt(K) > robust_gamma)] = robust_gamma**2
    
    K *= -gamma
    np.exp(K, K)    # exponentiate K in-place
    return K

def test_check_dense_matrices():
    # Ensure that pairwise array check works for dense matrices.
    # Check that if XB is None, XB is returned as reference to XA
    XA = np.resize(np.arange(40), (5, 8))
    XA_checked, XB_checked = check_pairwise_arrays(XA, None)
    assert XA_checked is XB_checked
    assert_array_equal(XA, XA_checked)

def histogram_intersection_kernel(X, Y = None, alpha = None, beta = None):
    """
    Source: https://github.com/kuantkid/scikit-learn/commit/16c82d8f2fe763df7bfee9bbcc40016fb84affcf
    Author: kuantkid
    Date: Nov 20, 2012

    Compute the histogram intersection kernel(min kernel) 
    between X and Y::
        K(x, y) = \\sum_i^n min(|x_i|^\x07lpha, |y_i|^\x08eta)
    Parameters
    ----------
    X : array of shape (n_samples_1, n_features)
    Y : array of shape (n_samples_2, n_features)
    gamma : float
    Returns
    -------
    Gram matrix : array of shape (n_samples_1, n_samples_2)
    """
    (X, Y,) = pairwise.check_pairwise_arrays(X, Y)
    if alpha is not None:
        X = np.abs(X) ** alpha
    if beta is not None:
        Y = np.abs(Y) ** beta
    (n_samples_1, n_features,) = X.shape
    (n_samples_2, _,) = Y.shape
    K = np.zeros(shape=(n_samples_1, n_samples_2), dtype=np.float)
    for i in range(n_samples_1):
        K[i] = np.sum(np.minimum(X[i], Y), axis=1)

    return K

def GeneralizedNormalKernel(X, Y=None, gamma = None, beta = 1):
    """Compute the generalized normal kernel between X and Y.
    The generalized normal kernel is defined as::
        K(x, y) = exp(-gamma ||x-y||_1^beta)
    for each pair of rows x in X and y in Y.
    Parameters
    ----------
    X : array of shape (n_samples_X, n_features)
    Y : array of shape (n_samples_Y, n_features)
    gamma : float
    Returns
    -------
    kernel_matrix : array of shape (n_samples_X, n_samples_Y)
    """

    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 1.0 / X.shape[1]

    if beta == 1:
        K = -gamma * manhattan_distances(X, Y)
    else:
        K = -gamma * manhattan_distances(X, Y) ** beta
    np.exp(K, K)    # exponentiate K in-place
    return K

def test_check_sparse_arrays():
    """ Ensures that checks return valid sparse matrices. """
    rng = np.random.RandomState(0)
    XA = rng.random_sample((5, 4))
    XA_sparse = csr_matrix(XA)
    XB = rng.random_sample((5, 4))
    XB_sparse = csr_matrix(XB)
    XA_checked, XB_checked = check_pairwise_arrays(XA_sparse, XB_sparse)

    # compare their difference because testing csr matrices for
    # equality with '==' does not work as expected.
    assert_true(abs(XA_sparse - XA_checked).nnz == 0)
    assert_true(abs(XB_sparse - XB_checked).nnz == 0)

    XA_checked, XB_checked = check_pairwise_arrays(XA_sparse, XA_sparse)
    assert_true(XA_sparse == XB_checked)
    assert_true(abs(XA_sparse - XA_checked).nnz == 0)

def test_check_XB_returned():
    """ Ensure that if XA and XB are given correctly, they return as equal."""
    # Check that if XB is not None, it is returned equal.
    # Note that the second dimension of XB is the same as XA.
    XA = np.resize(np.arange(40), (5, 8))
    XB = np.resize(np.arange(32), (4, 8))
    XA_checked, XB_checked = check_pairwise_arrays(XA, XB)
    assert_array_equal(XA, XA_checked)
    assert_array_equal(XB, XB_checked)

def test_check_sparse_arrays():
    """ Ensures that checks return valid sparse matrices. """
    rng = np.random.RandomState(0)
    XA = rng.random_sample((5, 4))
    XA_sparse = csr_matrix(XA)
    XB = rng.random_sample((5, 4))
    XB_sparse = csr_matrix(XB)
    XA_checked, XB_checked = check_pairwise_arrays(XA_sparse, XB_sparse)
    assert_equal(XA_sparse, XA_checked)
    assert_equal(XB_sparse, XB_checked)

def test_check_tuple_input():
    # Ensures that checks return valid tuples.
    rng = np.random.RandomState(0)
    XA = rng.random_sample((5, 4))
    XA_tuples = tuplify(XA)
    XB = rng.random_sample((5, 4))
    XB_tuples = tuplify(XB)
    XA_checked, XB_checked = check_pairwise_arrays(XA_tuples, XB_tuples)
    assert_array_equal(XA_tuples, XA_checked)
    assert_array_equal(XB_tuples, XB_checked)

def test_check_preserve_type():
    # Ensures that type float32 is preserved.
    XA = np.resize(np.arange(40), (5, 8)).astype(np.float32)
    XB = np.resize(np.arange(40), (5, 8)).astype(np.float32)

    XA_checked, XB_checked = check_pairwise_arrays(XA, None)
    assert_equal(XA_checked.dtype, np.float32)

    # both float32
    XA_checked, XB_checked = check_pairwise_arrays(XA, XB)
    assert_equal(XA_checked.dtype, np.float32)
    assert_equal(XB_checked.dtype, np.float32)

    # mismatched A
    XA_checked, XB_checked = check_pairwise_arrays(XA.astype(np.float), XB)
    assert_equal(XA_checked.dtype, np.float)
    assert_equal(XB_checked.dtype, np.float)

    # mismatched B
    XA_checked, XB_checked = check_pairwise_arrays(XA, XB.astype(np.float))
    assert_equal(XA_checked.dtype, np.float)
    assert_equal(XB_checked.dtype, np.float)

    def _cosine_distances_prenorm(self, X, Y):
        """
        Return cosine distances based on a prenormalized vectors.

        It allows for much faster computation of cosine distances.
        """
        if not self.prenorm:
            raise Exception(
                'Vectors must be prenormalized!')
        if Y is None:
            Y = X
        X, Y = smp.check_pairwise_arrays(X, Y)
        sims = X.dot(Y.T)

        if scipy.sparse.issparse(sims):
            sims = sims.todense()

        return 1 - sims

def first_periodic_kernel(X, Y=None, gamma=None, period=None):
    # TODO: Add mathematical form of the kernel in the docstring
    """Compute the first periodic kernel between *X* and *Y*.

    Parameters
    ----------
    X : array of shape (n_samples_X, n_features)

    Y : array of shape (n_samples_Y, n_features)

    gamma : float, default None
        If None, default to 1.0 / n_samples_X

    period : float, default None
        If None, default to 2 * pi.

        This parameter should not be default as
        wrong estimation lead to poor learning score.

    Returns
    -------
    kernel_matrix : array of shape (n_samples_X, n_samples_Y)
    """
    X, Y = check_pairwise_arrays(X, Y)
    if gamma is None:
        gamma = 0.8

    if period is None:
        period = 2. * pi

    a = -log(gamma) / period
    b = 2 * pi / period
    c = sqrt(pi / a) * (exp(- b ** 2 / (4 * a)) + 1)
    K = euclidean_distances(X, Y, squared=True)

    # TODO: Optimize to avoid temporary?
    return exp(-a * K) * (1 + cos(b * sqrt(K))) / c

def cosine_similarity(X, Y=None):
    """Compute cosine similarity between samples in X and Y.

    Cosine similarity, or the cosine kernel, computes similarity as the
    normalized dot product of X and Y:

        K(X, Y) = <X, Y> / (||X||*||Y||)

    On L2-normalized data, this function is equivalent to linear_kernel.

    Parameters
    ----------
    X : array_like, sparse matrix
        with shape (n_samples_X, n_features).

    Y : array_like, sparse matrix (optional)
        with shape (n_samples_Y, n_features).

    Returns
    -------
    kernel matrix : array_like
        An array with shape (n_samples_X, n_samples_Y).
    """
    # to avoid recursive import

    X, Y = check_pairwise_arrays(X, Y)

    X_normalized = normalize(X, copy=True)
    if X is Y:
        Y_normalized = X_normalized
    else:
        Y_normalized = normalize(Y, copy=True)

    K = linear_kernel(X_normalized, Y_normalized)

    return K

def _prep_X_Y_for_cython(X, Y):
    X, Y = check_pairwise_arrays(X, Y)
    X, Y = X.astype(np.double, order='C'), Y.astype(np.double, order='C').T  # transposing Y here!
    res = np.zeros((X.shape[0], Y.shape[1]), dtype=X.dtype)
    return X, Y, res

def pairwise_distances(X, Y=None, metric="euclidean", n_jobs=1, **kwds):
    """ Compute the distance matrix from a vector array X and optional Y.

    This method takes either a vector array or a distance matrix, and returns
    a distance matrix. If the input is a vector array, the distances are
    computed. If the input is a distances matrix, it is returned instead.

    This method provides a safe way to take a distance matrix as input, while
    preserving compatibility with many other algorithms that take a vector
    array.

    If Y is given (default is None), then the returned matrix is the pairwise
    distance between the arrays from both X and Y.

    Valid values for metric are:

    - From scikit-learn: ['cityblock', 'cosine', 'euclidean', 'l1', 'l2',
      'manhattan']. These metrics support sparse matrix inputs.

    - From scipy.spatial.distance: ['braycurtis', 'canberra', 'chebyshev',
      'correlation', 'dice', 'hamming', 'jaccard', 'kulsinski', 'mahalanobis',
      'matching', 'minkowski', 'rogerstanimoto', 'russellrao', 'seuclidean',
      'sokalmichener', 'sokalsneath', 'sqeuclidean', 'yule']
      See the documentation for scipy.spatial.distance for details on these
      metrics. These metrics do not support sparse matrix inputs.

    Note that in the case of 'cityblock', 'cosine' and 'euclidean' (which are
    valid scipy.spatial.distance metrics), the scikit-learn implementation
    will be used, which is faster and has support for sparse matrices (except
    for 'cityblock'). For a verbose description of the metrics from
    scikit-learn, see the __doc__ of the sklearn.pairwise.distance_metrics
    function.

    Read more in the :ref:`User Guide <metrics>`.

    Parameters
    ----------
    X : array [n_samples_a, n_samples_a] if metric == "precomputed", or, \
             [n_samples_a, n_features] otherwise
        Array of pairwise distances between samples, or a feature array.

    Y : array [n_samples_b, n_features], optional
        An optional second feature array. Only allowed if metric != "precomputed".

    metric : string, or callable
        The metric to use when calculating distance between instances in a
        feature array. If metric is a string, it must be one of the options
        allowed by scipy.spatial.distance.pdist for its metric parameter, or
        a metric listed in pairwise.PAIRWISE_DISTANCE_FUNCTIONS.
        If metric is "precomputed", X is assumed to be a distance matrix.
        Alternatively, if metric is a callable function, it is called on each
        pair of instances (rows) and the resulting value recorded. The callable
        should take two arrays from X as input and return a value indicating
        the distance between them.

    n_jobs : int
        The number of jobs to use for the computation. This works by breaking
        down the pairwise matrix into n_jobs even slices and computing them in
        parallel.

        If -1 all CPUs are used. If 1 is given, no parallel computing code is
        used at all, which is useful for debugging. For n_jobs below -1,
        (n_cpus + 1 + n_jobs) are used. Thus for n_jobs = -2, all CPUs but one
        are used.

    `**kwds` : optional keyword parameters
        Any further parameters are passed directly to the distance function.
        If using a scipy.spatial.distance metric, the parameters are still
        metric dependent. See the scipy docs for usage examples.

    Returns
    -------
    D : array [n_samples_a, n_samples_a] or [n_samples_a, n_samples_b]
        A distance matrix D such that D_{i, j} is the distance between the
        ith and jth vectors of the given matrix X, if Y is None.
        If Y is not None, then D_{i, j} is the distance between the ith array
        from X and the jth array from Y.

    """
    if (metric not in _VALID_METRICS and
            not callable(metric) and metric != "precomputed"):
        raise ValueError("Unknown metric %s. "
                         "Valid metrics are %s, or 'precomputed', or a "
                         "callable" % (metric, _VALID_METRICS))

    if metric == "precomputed":
        X, _ = check_pairwise_arrays(X, Y, precomputed=True)
        return X
    elif metric in PAIRWISE_DISTANCE_FUNCTIONS:
        func = PAIRWISE_DISTANCE_FUNCTIONS[metric]
    elif callable(metric):
        func = partial(_pairwise_callable, metric=metric, **kwds)
    else:
        if issparse(X) or issparse(Y):
            raise TypeError("scipy distance metrics do not"
                            " support sparse matrices.")

        dtype = bool if metric in PAIRWISE_BOOLEAN_FUNCTIONS else None

        X, Y = check_pairwise_arrays(X, Y, dtype=dtype)
#.........这里部分代码省略.........