Repository: scikit-learn/scikit-learn
Pull request: [MRG+2] Allowing Gaussian process kernels on structured data - updated
Fixes #13936.

This PR is a from-scratch refresh of #13954 (with all previous modifications incorporated) based on the latest master branch.

Made the Gaussian process regressor and classifier compatible with generic kernels on arbitrary data. Specifically:
1. Added a `requires_vector_input` property to all the built-in GP kernels in `gaussian_process/kernels.py`.
2. Added a `GenericKernelMixin` base classes to overload the `requires_vector_input` property for kernels that can work on structured and/or generic data.
3. Let the GP regressor and classifier use different `check_array` and `check_X_y` parameters for vector/generic input. I.e., do not force the X and Y array to be at least 2D and numeric if the kernel can operate on generic input. Updated docstrings accordingly.
4. Provided a minimal example of GP regression and classification using a sequence kernel on variable-length strings in `plot_gpr_on_structured_data.py`.


Problem statement:
GaussianProcessRegressor unnecessarily forcing training samples to be array_like
The current interface precludes the application of the GPR module to cases where the training samples are sequences, trees, graphs etc. In fact, GPR should not care about whether its input are fixed-length feature vectors or not --- that should be handled by the kernel.

Expected usage scenario:
```
k = StructuredDataKernel()
gp = GaussianProcessRegressor( k )
X = [ object1, objects2, ... ]
Y = [ object3, objects4, ... ]
y = [ scalar_values... ]
gp.fit( X, y ) 
gp.predict( Y )
```
Error information (`sklearn/utils/validation.py", line 448, in check_array`):
```
ValueError: setting an array element with a sequence.
```
Expect outcome:
The learning workflow should proceed just fine as long as k returns valid inner product matrices between X and itself or Y.

Hints from issue discussion:
If you'd like to submit a pull request implementing that support (perhaps
with an example demonstrating a custom kernel) I think it would be well
received

OK. I'm working on it.

Commit message: documentation update

Patch to review:
diff --git a/doc/whats_new/v0.22.rst b/doc/whats_new/v0.22.rst
index 525ae77311064..4d5ee03853e16 100644
--- a/doc/whats_new/v0.22.rst
+++ b/doc/whats_new/v0.22.rst
@@ -405,6 +405,15 @@ Changelog
   estimator's constructor but not stored as attributes on the instance.
   :pr:`14464` by `Joel Nothman`_.
 
+- |Feature| Gaussian process models on structured data: :class:`gaussian_process.GaussianProcessRegressor`
+  and :class:`gaussian_process.GaussianProcessClassifier` can now accept a list
+  of generic objects (e.g. strings, trees, graphs, etc.) as the ``X`` argument
+  to their training/prediction methods.
+  A user-defined kernel should be provided for computing the kernel matrix among
+  the generic objects, and should inherit from :class:`gaussian_process.kernels.GenericKernelMixin`
+  to notify the GPR/GPC model that it handles non-vectorial samples.
+  :pr:`15557` by :user:`Yu-Hang Tang <yhtang>`.
+
 :mod:`sklearn.impute`
 .....................
 
diff --git a/examples/gaussian_process/plot_gpr_on_structured_data.py b/examples/gaussian_process/plot_gpr_on_structured_data.py
new file mode 100644
index 0000000000000..160ff0bb4bc45
--- /dev/null
+++ b/examples/gaussian_process/plot_gpr_on_structured_data.py
@@ -0,0 +1,172 @@
+"""
+==========================================================================
+Gaussian processes on discrete data structures
+==========================================================================
+
+This example illustrates the use of Gaussian processes for regression and
+classification tasks on data that are not in fixed-length feature vector form.
+This is achieved through the use of kernel functions that operates directly
+on discrete structures such as variable-length sequences, trees, and graphs.
+
+Specifically, here the input variables are some gene sequences stored as
+variable-length strings consisting of letters 'A', 'T', 'C', and 'G',
+while the output variables are floating point numbers and True/False labels
+in the regression and classification tasks, respectively.
+
+A kernel between the gene sequences is defined using R-convolution [1]_ by
+integrating a binary letter-wise kernel over all pairs of letters among a pair
+of strings.
+
+This example will generate three figures.
+
+In the first figure, we visualize the value of the kernel, i.e. the similarity
+of the sequences, using a colormap. Brighter color here indicates higher
+similarity.
+
+In the second figure, we show some regression result on a dataset of 6
+sequences. Here we use the 1st, 2nd, 4th, and 5th sequences as the training set
+to make predictions on the 3rd and 6th sequences.
+
+In the third figure, we demonstrate a classification model by training on 6
+sequences and make predictions on another 5 sequences. The ground truth here is
+simply  whether there is at least one 'A' in the sequence. Here the model makes
+four correct classifications and fails on one.
+
+.. [1] Haussler, D. (1999). Convolution kernels on discrete structures
+(Vol. 646). Technical report, Department of Computer Science, University of
+California at Santa Cruz.
+"""
+print(__doc__)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
+from sklearn.gaussian_process.kernels import GenericKernelMixin
+from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process import GaussianProcessClassifier
+from sklearn.base import clone
+
+
+class SequenceKernel(GenericKernelMixin, Kernel):
+    '''
+    A minimal (but valid) convolutional kernel for sequences of variable
+    lengths.'''
+    def __init__(self,
+                 baseline_similarity=0.5,
+                 baseline_similarity_bounds=(1e-5, 1)):
+        self.baseline_similarity = baseline_similarity
+        self.baseline_similarity_bounds = baseline_similarity_bounds
+
+    @property
+    def hyperparameter_baseline_similarity(self):
+        return Hyperparameter("baseline_similarity",
+                              "numeric",
+                              self.baseline_similarity_bounds)
+
+    def _f(self, s1, s2):
+        '''
+        kernel value between a pair of sequences
+        '''
+        return sum([1.0 if c1 == c2 else self.baseline_similarity
+                   for c1 in s1
+                   for c2 in s2])
+
+    def _g(self, s1, s2):
+        '''
+        kernel derivative between a pair of sequences
+        '''
+        return sum([0.0 if c1 == c2 else 1.0
+                    for c1 in s1
+                    for c2 in s2])
+
+    def __call__(self, X, Y=None, eval_gradient=False):
+        if Y is None:
+            Y = X
+
+        if eval_gradient:
+            return (np.array([[self._f(x, y) for y in Y] for x in X]),
+                    np.array([[[self._g(x, y)] for y in Y] for x in X]))
+        else:
+            return np.array([[self._f(x, y) for y in Y] for x in X])
+
+    def diag(self, X):
+        return np.array([self._f(x, x) for x in X])
+
+    def is_stationary(self):
+        return False
+
+    def clone_with_theta(self, theta):
+        cloned = clone(self)
+        cloned.theta = theta
+        return cloned
+
+
+kernel = SequenceKernel()
+
+X = np.array(['AGCT', 'AGC', 'AACT', 'TAA', 'AAA', 'GAACA'])
+y = np.array([1.0, 1.0, 2.0, 2.0, 3.0, 3.0])
+
+'''
+Sequence similarity matrix under the kernel
+===========================================
+'''
+
+K = kernel(X)
+D = kernel.diag(X)
+
+plt.figure(figsize=(8, 5))
+plt.imshow(np.diag(D**-0.5).dot(K).dot(np.diag(D**-0.5)))
+plt.xticks(np.arange(len(X)), X)
+plt.yticks(np.arange(len(X)), X)
+plt.title('Sequence similarity under the kernel')
+
+'''
+Regression
+==========
+'''
+
+training_idx = [0, 1, 3, 4]
+gp = GaussianProcessRegressor(kernel)
+gp.fit(X[training_idx], y[training_idx])
+
+plt.figure(figsize=(8, 5))
+plt.bar(np.arange(len(X)), gp.predict(X), color='b', label='prediction')
+plt.bar(training_idx, y[training_idx], width=0.2, color='r',
+        alpha=0.5, label='training')
+plt.xticks(np.arange(len(X)), X)
+plt.title('Regression on sequences')
+plt.legend()
+
+'''
+Classification
+==============
+'''
+
+X = np.array(['AGCT', 'CGA', 'TAAC', 'TCG', 'CTTT', 'TGCT'])
+# whether there are 'A's in the sequence
+t = np.array([True, True, True, False, False, False])
+
+gp = GaussianProcessClassifier(kernel)
+gp.fit(X, t)
+
+seqs_test = ['AAA', 'ATAG', 'CTC', 'CT', 'C']
+clss_test = [True, True, False, False, False]
+
+plt.figure(figsize=(8, 5))
+plt.scatter(np.arange(len(X)), [1.0 if c else -1.0 for c in t],
+            s=100, marker='o', edgecolor='none', facecolor=(1, 0.75, 0),
+            label='training')
+plt.scatter(len(X) + np.arange(len(seqs_test)),
+            [1.0 if c else -1.0 for c in clss_test],
+            s=100, marker='o', edgecolor='none', facecolor='r', label='truth')
+plt.scatter(len(X) + np.arange(len(seqs_test)),
+            [1.0 if c else -1.0 for c in gp.predict(seqs_test)],
+            s=100, marker='x', edgecolor=(0, 1.0, 0.3), linewidth=2,
+            label='prediction')
+plt.xticks(np.arange(len(X) + len(seqs_test)),
+           np.concatenate((X, seqs_test)))
+plt.yticks([-1, 1], [False, True])
+plt.title('Classification on sequences')
+plt.legend()
+
+plt.show()
diff --git a/sklearn/base.py b/sklearn/base.py
index 7ededd7a70548..9cb3efcef94af 100644
--- a/sklearn/base.py
+++ b/sklearn/base.py
@@ -383,8 +383,9 @@ def score(self, X, y, sample_weight=None):
         ----------
         X : array-like of shape (n_samples, n_features)
             Test samples. For some estimators this may be a
-            precomputed kernel matrix instead, shape = (n_samples,
-            n_samples_fitted], where n_samples_fitted is the number of
+            precomputed kernel matrix or a list of generic objects instead,
+            shape = (n_samples, n_samples_fitted],
+            where n_samples_fitted is the number of
             samples used in the fitting for the estimator.
 
         y : array-like of shape (n_samples,) or (n_samples, n_outputs)
diff --git a/sklearn/gaussian_process/_gpc.py b/sklearn/gaussian_process/_gpc.py
index 072cf80dba250..f0812633160ef 100644
--- a/sklearn/gaussian_process/_gpc.py
+++ b/sklearn/gaussian_process/_gpc.py
@@ -115,8 +115,10 @@ def optimizer(obj_func, initial_theta, bounds):
 
     Attributes
     ----------
-    X_train_ : array-like of shape (n_samples, n_features)
-        Feature values in training data (also required for prediction)
+    X_train_ : sequence of length n_samples
+        Feature vectors or other representations of training data (also
+        required for prediction). Could either be array-like with shape =
+        (n_samples, n_features) or a list of objects
 
     y_train_ : array-like of shape (n_samples,)
         Target values in training data (also required for prediction)
@@ -160,8 +162,10 @@ def fit(self, X, y):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
-            Training data
+        X : sequence of length n_samples
+            Feature vectors or other representations of training data.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         y : array-like of shape (n_samples,)
             Target values, must be binary
@@ -248,7 +252,10 @@ def predict(self, X):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
+        X : sequence of length n_samples
+            Query points where the GP is evaluated for classification.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         Returns
         -------
@@ -270,7 +277,10 @@ def predict_proba(self, X):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
+        X : sequence of length n_samples
+            Query points where the GP is evaluated for classification.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         Returns
         -------
@@ -602,8 +612,10 @@ def fit(self, X, y):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
-            Training data
+        X : sequence of length n_samples
+            Feature vectors or other representations of training data.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         y : array-like of shape (n_samples,)
             Target values, must be binary
@@ -612,7 +624,12 @@ def fit(self, X, y):
         -------
         self : returns an instance of self.
         """
-        X, y = check_X_y(X, y, multi_output=False)
+        if self.kernel is None or self.kernel.requires_vector_input:
+            X, y = check_X_y(X, y, multi_output=False,
+                             ensure_2d=True, dtype="numeric")
+        else:
+            X, y = check_X_y(X, y, multi_output=False,
+                             ensure_2d=False, dtype=None)
 
         self.base_estimator_ = _BinaryGaussianProcessClassifierLaplace(
             self.kernel, self.optimizer, self.n_restarts_optimizer,
@@ -656,7 +673,10 @@ def predict(self, X):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
+        X : sequence of length n_samples
+            Query points where the GP is evaluated for classification.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         Returns
         -------
@@ -664,7 +684,12 @@ def predict(self, X):
             Predicted target values for X, values are from ``classes_``
         """
         check_is_fitted(self)
-        X = check_array(X)
+
+        if self.kernel is None or self.kernel.requires_vector_input:
+            X = check_array(X, ensure_2d=True, dtype="numeric")
+        else:
+            X = check_array(X, ensure_2d=False, dtype=None)
+
         return self.base_estimator_.predict(X)
 
     def predict_proba(self, X):
@@ -672,7 +697,10 @@ def predict_proba(self, X):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
+        X : sequence of length n_samples
+            Query points where the GP is evaluated for classification.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         Returns
         -------
@@ -686,7 +714,12 @@ def predict_proba(self, X):
             raise ValueError("one_vs_one multi-class mode does not support "
                              "predicting probability estimates. Use "
                              "one_vs_rest mode instead.")
-        X = check_array(X)
+
+        if self.kernel is None or self.kernel.requires_vector_input:
+            X = check_array(X, ensure_2d=True, dtype="numeric")
+        else:
+            X = check_array(X, ensure_2d=False, dtype=None)
+
         return self.base_estimator_.predict_proba(X)
 
     @property
diff --git a/sklearn/gaussian_process/_gpr.py b/sklearn/gaussian_process/_gpr.py
index a2be69abff794..1373444a2c9ae 100644
--- a/sklearn/gaussian_process/_gpr.py
+++ b/sklearn/gaussian_process/_gpr.py
@@ -114,8 +114,10 @@ def optimizer(obj_func, initial_theta, bounds):
 
     Attributes
     ----------
-    X_train_ : array-like of shape (n_samples, n_features)
-        Feature values in training data (also required for prediction)
+    X_train_ : sequence of length n_samples
+        Feature vectors or other representations of training data (also
+        required for prediction). Could either be array-like with shape =
+        (n_samples, n_features) or a list of objects
 
     y_train_ : array-like of shape (n_samples,) or (n_samples, n_targets)
         Target values in training data (also required for prediction)
@@ -164,8 +166,10 @@ def fit(self, X, y):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
-            Training data
+        X : sequence of length n_samples
+            Feature vectors or other representations of training data.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         y : array-like of shape (n_samples,) or (n_samples, n_targets)
             Target values
@@ -182,7 +186,12 @@ def fit(self, X, y):
 
         self._rng = check_random_state(self.random_state)
 
-        X, y = check_X_y(X, y, multi_output=True, y_numeric=True)
+        if self.kernel_.requires_vector_input:
+            X, y = check_X_y(X, y, multi_output=True, y_numeric=True,
+                             ensure_2d=True, dtype="numeric")
+        else:
+            X, y = check_X_y(X, y, multi_output=True, y_numeric=True,
+                             ensure_2d=False, dtype=None)
 
         # Normalize target value
         if self.normalize_y:
@@ -273,8 +282,10 @@ def predict(self, X, return_std=False, return_cov=False):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples, n_features)
-            Query points where the GP is evaluated
+        X : sequence of length n_samples
+            Query points where the GP is evaluated.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         return_std : bool, default: False
             If True, the standard-deviation of the predictive distribution at
@@ -302,7 +313,10 @@ def predict(self, X, return_std=False, return_cov=False):
                 "Not returning standard deviation of predictions when "
                 "returning full covariance.")
 
-        X = check_array(X)
+        if self.kernel is None or self.kernel.requires_vector_input:
+            X = check_array(X, ensure_2d=True, dtype="numeric")
+        else:
+            X = check_array(X, ensure_2d=False, dtype=None)
 
         if not hasattr(self, "X_train_"):  # Unfitted;predict based on GP prior
             if self.kernel is None:
@@ -357,8 +371,10 @@ def sample_y(self, X, n_samples=1, random_state=0):
 
         Parameters
         ----------
-        X : array-like of shape (n_samples_X, n_features)
-            Query points where the GP samples are evaluated
+        X : sequence of length n_samples
+            Query points where the GP is evaluated.
+            Could either be array-like with shape = (n_samples, n_features)
+            or a list of objects.
 
         n_samples : int, default: 1
             The number of samples drawn from the Gaussian process
diff --git a/sklearn/gaussian_process/kernels.py b/sklearn/gaussian_process/kernels.py
index 1634113a009f3..7601eb97a68ef 100644
--- a/sklearn/gaussian_process/kernels.py
+++ b/sklearn/gaussian_process/kernels.py
@@ -31,6 +31,7 @@
 
 from ..metrics.pairwise import pairwise_kernels
 from ..base import clone
+from ..utils.validation import _num_samples
 
 
 def _check_length_scale(X, length_scale):
@@ -352,7 +353,7 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples
             Left argument of the returned kernel k(X, Y)
 
         Returns
@@ -365,6 +366,13 @@ def diag(self, X):
     def is_stationary(self):
         """Returns whether the kernel is stationary. """
 
+    @property
+    def requires_vector_input(self):
+        """Returns whether the kernel is defined on fixed-length feature
+        vectors or generic objects. Defaults to True for backward
+        compatibility."""
+        return True
+
 
 class NormalizedKernelMixin:
     """Mixin for kernels which are normalized: k(X, X)=1.
@@ -381,7 +389,7 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples
             Left argument of the returned kernel k(X, Y)
 
         Returns
@@ -403,6 +411,19 @@ def is_stationary(self):
         return True
 
 
+class GenericKernelMixin:
+    """Mixin for kernels which operate on generic objects such as variable-
+    length sequences, trees, and graphs.
+
+    .. versionadded:: 0.22
+    """
+
+    @property
+    def requires_vector_input(self):
+        """whether the kernel works only on fixed-length feature vectors."""
+        return False
+
+
 class CompoundKernel(Kernel):
     """Kernel which is composed of a set of other kernels.
 
@@ -481,12 +502,15 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples_X
             Left argument of the returned kernel k(X, Y)
+            Could either be array-like with shape = (n_samples_X, n_features)
+            or a list of objects.
 
-        Y : array, shape (n_samples_Y, n_features), (optional, default=None)
+        Y : sequence of length n_samples_Y
             Right argument of the returned kernel k(X, Y). If None, k(X, X)
-            if evaluated instead.
+            is evaluated instead. Y could either be array-like with
+            shape = (n_samples_Y, n_features) or a list of objects.
 
         eval_gradient : bool (optional, default=False)
             Determines whether the gradient with respect to the kernel
@@ -524,6 +548,12 @@ def is_stationary(self):
         """Returns whether the kernel is stationary. """
         return np.all([kernel.is_stationary() for kernel in self.kernels])
 
+    @property
+    def requires_vector_input(self):
+        """Returns whether the kernel is defined on discrete structures. """
+        return np.any([kernel.requires_vector_input
+                       for kernel in self.kernels])
+
     def diag(self, X):
         """Returns the diagonal of the kernel k(X, X).
 
@@ -533,8 +563,9 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
-            Left argument of the returned kernel k(X, Y)
+        X : sequence of length n_samples_X
+            Argument to the kernel. Could either be array-like with
+            shape = (n_samples_X, n_features) or a list of objects.
 
         Returns
         -------
@@ -646,6 +677,12 @@ def is_stationary(self):
         """Returns whether the kernel is stationary. """
         return self.k1.is_stationary() and self.k2.is_stationary()
 
+    @property
+    def requires_vector_input(self):
+        """Returns whether the kernel is stationary. """
+        return (self.k1.requires_vector_input or
+                self.k2.requires_vector_input)
+
 
 class Sum(KernelOperator):
     """Sum-kernel k1 + k2 of two kernels k1 and k2.
@@ -670,12 +707,15 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples_X
             Left argument of the returned kernel k(X, Y)
+            Could either be array-like with shape = (n_samples_X, n_features)
+            or a list of objects.
 
-        Y : array, shape (n_samples_Y, n_features), (optional, default=None)
+        Y : sequence of length n_samples_Y
             Right argument of the returned kernel k(X, Y). If None, k(X, X)
-            if evaluated instead.
+            is evaluated instead. Y could either be array-like with
+            shape = (n_samples_Y, n_features) or a list of objects.
 
         eval_gradient : bool (optional, default=False)
             Determines whether the gradient with respect to the kernel
@@ -707,8 +747,9 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
-            Left argument of the returned kernel k(X, Y)
+        X : sequence of length n_samples_X
+            Argument to the kernel. Could either be array-like with
+            shape = (n_samples_X, n_features) or a list of objects.
 
         Returns
         -------
@@ -744,12 +785,15 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples_X
             Left argument of the returned kernel k(X, Y)
+            Could either be array-like with shape = (n_samples_X, n_features)
+            or a list of objects.
 
-        Y : array, shape (n_samples_Y, n_features), (optional, default=None)
+        Y : sequence of length n_samples_Y
             Right argument of the returned kernel k(X, Y). If None, k(X, X)
-            if evaluated instead.
+            is evaluated instead. Y could either be array-like with
+            shape = (n_samples_Y, n_features) or a list of objects.
 
         eval_gradient : bool (optional, default=False)
             Determines whether the gradient with respect to the kernel
@@ -782,8 +826,9 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
-            Left argument of the returned kernel k(X, Y)
+        X : sequence of length n_samples_X
+            Argument to the kernel. Could either be array-like with
+            shape = (n_samples_X, n_features) or a list of objects.
 
         Returns
         -------
@@ -896,12 +941,15 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples_X
             Left argument of the returned kernel k(X, Y)
+            Could either be array-like with shape = (n_samples_X, n_features)
+            or a list of objects.
 
-        Y : array, shape (n_samples_Y, n_features), (optional, default=None)
+        Y : sequence of length n_samples_Y
             Right argument of the returned kernel k(X, Y). If None, k(X, X)
-            if evaluated instead.
+            is evaluated instead. Y could either be array-like with
+            shape = (n_samples_Y, n_features) or a list of objects.
 
         eval_gradient : bool (optional, default=False)
             Determines whether the gradient with respect to the kernel
@@ -935,8 +983,9 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
-            Left argument of the returned kernel k(X, Y)
+        X : sequence of length n_samples_X
+            Argument to the kernel. Could either be array-like with
+            shape = (n_samples_X, n_features) or a list of objects.
 
         Returns
         -------
@@ -952,8 +1001,14 @@ def is_stationary(self):
         """Returns whether the kernel is stationary. """
         return self.kernel.is_stationary()
 
+    @property
+    def requires_vector_input(self):
+        """Returns whether the kernel is defined on discrete structures. """
+        return self.kernel.requires_vector_input
+
 
-class ConstantKernel(StationaryKernelMixin, Kernel):
+class ConstantKernel(StationaryKernelMixin, GenericKernelMixin,
+                     Kernel):
     """Constant kernel.
 
     Can be used as part of a product-kernel where it scales the magnitude of
@@ -988,12 +1043,15 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples_X
             Left argument of the returned kernel k(X, Y)
+            Could either be array-like with shape = (n_samples_X, n_features)
+            or a list of objects.
 
-        Y : array, shape (n_samples_Y, n_features), (optional, default=None)
+        Y : sequence of length n_samples_Y
             Right argument of the returned kernel k(X, Y). If None, k(X, X)
-            if evaluated instead.
+            is evaluated instead. Y could either be array-like with
+            shape = (n_samples_Y, n_features) or a list of objects.
 
         eval_gradient : bool (optional, default=False)
             Determines whether the gradient with respect to the kernel
@@ -1009,21 +1067,20 @@ def __call__(self, X, Y=None, eval_gradient=False):
             hyperparameter of the kernel. Only returned when eval_gradient
             is True.
         """
-        X = np.atleast_2d(X)
         if Y is None:
             Y = X
         elif eval_gradient:
             raise ValueError("Gradient can only be evaluated when Y is None.")
 
-        K = np.full((X.shape[0], Y.shape[0]), self.constant_value,
+        K = np.full((_num_samples(X), _num_samples(Y)), self.constant_value,
                     dtype=np.array(self.constant_value).dtype)
         if eval_gradient:
             if not self.hyperparameter_constant_value.fixed:
-                return (K, np.full((X.shape[0], X.shape[0], 1),
+                return (K, np.full((_num_samples(X), _num_samples(X), 1),
                                    self.constant_value,
                                    dtype=np.array(self.constant_value).dtype))
             else:
-                return K, np.empty((X.shape[0], X.shape[0], 0))
+                return K, np.empty((_num_samples(X), _num_samples(X), 0))
         else:
             return K
 
@@ -1036,22 +1093,24 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
-            Left argument of the returned kernel k(X, Y)
+        X : sequence of length n_samples_X
+            Argument to the kernel. Could either be array-like with
+            shape = (n_samples_X, n_features) or a list of objects.
 
         Returns
         -------
         K_diag : array, shape (n_samples_X,)
             Diagonal of kernel k(X, X)
         """
-        return np.full(X.shape[0], self.constant_value,
+        return np.full(_num_samples(X), self.constant_value,
                        dtype=np.array(self.constant_value).dtype)
 
     def __repr__(self):
         return "{0:.3g}**2".format(np.sqrt(self.constant_value))
 
 
-class WhiteKernel(StationaryKernelMixin, Kernel):
+class WhiteKernel(StationaryKernelMixin, GenericKernelMixin,
+                  Kernel):
     """White kernel.
 
     The main use-case of this kernel is as part of a sum-kernel where it
@@ -1085,12 +1144,15 @@ def __call__(self, X, Y=None, eval_gradient=False):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
+        X : sequence of length n_samples_X
             Left argument of the returned kernel k(X, Y)
+            Could either be array-like with shape = (n_samples_X, n_features)
+            or a list of objects.
 
-        Y : array, shape (n_samples_Y, n_features), (optional, default=None)
+        Y : sequence of length n_samples_Y
             Right argument of the returned kernel k(X, Y). If None, k(X, X)
-            if evaluated instead.
+            is evaluated instead. Y could either be array-like with
+            shape = (n_samples_Y, n_features) or a list of objects.
 
         eval_gradient : bool (optional, default=False)
             Determines whether the gradient with respect to the kernel
@@ -1106,22 +1168,21 @@ def __call__(self, X, Y=None, eval_gradient=False):
             hyperparameter of the kernel. Only returned when eval_gradient
             is True.
         """
-        X = np.atleast_2d(X)
         if Y is not None and eval_gradient:
             raise ValueError("Gradient can only be evaluated when Y is None.")
 
         if Y is None:
-            K = self.noise_level * np.eye(X.shape[0])
+            K = self.noise_level * np.eye(_num_samples(X))
             if eval_gradient:
                 if not self.hyperparameter_noise_level.fixed:
                     return (K, self.noise_level
-                            * np.eye(X.shape[0])[:, :, np.newaxis])
+                            * np.eye(_num_samples(X))[:, :, np.newaxis])
                 else:
-                    return K, np.empty((X.shape[0], X.shape[0], 0))
+                    return K, np.empty((_num_samples(X), _num_samples(X), 0))
             else:
                 return K
         else:
-            return np.zeros((X.shape[0], Y.shape[0]))
+            return np.zeros((_num_samples(X), _num_samples(Y)))
 
     def diag(self, X):
         """Returns the diagonal of the kernel k(X, X).
@@ -1132,15 +1193,16 @@ def diag(self, X):
 
         Parameters
         ----------
-        X : array, shape (n_samples_X, n_features)
-            Left argument of the returned kernel k(X, Y)
+        X : sequence of length n_samples_X
+            Argument to the kernel. Could either be array-like with
+            shape = (n_samples_X, n_features) or a list of objects.
 
         Returns
         -------
         K_diag : array, shape (n_samples_X,)
             Diagonal of kernel k(X, X)
         """
-        return np.full(X.shape[0], self.noise_level,
+        return np.full(_num_samples(X), self.noise_level,
                        dtype=np.array(self.noise_level).dtype)
 
     def __repr__(self):
diff --git a/sklearn/gaussian_process/tests/_mini_sequence_kernel.py b/sklearn/gaussian_process/tests/_mini_sequence_kernel.py
new file mode 100644
index 0000000000000..5c833d38750b6
--- /dev/null
+++ b/sklearn/gaussian_process/tests/_mini_sequence_kernel.py
@@ -0,0 +1,50 @@
+from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
+from sklearn.gaussian_process.kernels import GenericKernelMixin
+from sklearn.gaussian_process.kernels import StationaryKernelMixin
+import numpy as np
+from sklearn.base import clone
+
+
+class MiniSeqKernel(GenericKernelMixin,
+                    StationaryKernelMixin,
+                    Kernel):
+    '''
+    a mimimal (but valid) convolutional kernel for sequences of variable length
+    '''
+    def __init__(self,
+                 baseline_similarity=0.5,
+                 baseline_similarity_bounds=(1e-5, 1)):
+        self.baseline_similarity = baseline_similarity
+        self.baseline_similarity_bounds = baseline_similarity_bounds
+
+    @property
+    def hyperparameter_baseline_similarity(self):
+        return Hyperparameter("baseline_similarity",
+                              "numeric",
+                              self.baseline_similarity_bounds)
+
+    def _f(self, s1, s2):
+        return sum([1.0 if c1 == c2 else self.baseline_similarity
+                   for c1 in s1
+                   for c2 in s2])
+
+    def _g(self, s1, s2):
+        return sum([0.0 if c1 == c2 else 1.0 for c1 in s1 for c2 in s2])
+
+    def __call__(self, X, Y=None, eval_gradient=False):
+        if Y is None:
+            Y = X
+
+        if eval_gradient:
+            return (np.array([[self._f(x, y) for y in Y] for x in X]),
+                    np.array([[[self._g(x, y)] for y in Y] for x in X]))
+        else:
+            return np.array([[self._f(x, y) for y in Y] for x in X])
+
+    def diag(self, X):
+        return np.array([self._f(x, x) for x in X])
+
+    def clone_with_theta(self, theta):
+        cloned = clone(self)
+        cloned.theta = theta
+        return cloned
diff --git a/sklearn/gaussian_process/tests/test_gpc.py b/sklearn/gaussian_process/tests/test_gpc.py
index aec5cc147223f..72d550231f4ea 100644
--- a/sklearn/gaussian_process/tests/test_gpc.py
+++ b/sklearn/gaussian_process/tests/test_gpc.py
@@ -11,12 +11,15 @@
 
 from sklearn.gaussian_process import GaussianProcessClassifier
 from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
+from sklearn.gaussian_process.tests._mini_sequence_kernel import MiniSeqKernel
 
 from sklearn.utils._testing import assert_almost_equal, assert_array_equal
 
 
 def f(x):
     return np.sin(x)
+
+
 X = np.atleast_2d(np.linspace(0, 10, 30)).T
 X2 = np.atleast_2d([2., 4., 5.5, 6.5, 7.5]).T
 y = np.array(f(X).ravel() > 0, dtype=int)
@@ -44,12 +47,22 @@ def test_predict_consistent(kernel):
                        gpc.predict_proba(X)[:, 1] >= 0.5)
 
 
+def test_predict_consistent_structured():
+    # Check binary predict decision has also predicted probability above 0.5.
+    X = ['A', 'AB', 'B']
+    y = np.array([True, False, True])
+    kernel = MiniSeqKernel(baseline_similarity_bounds='fixed')
+    gpc = GaussianProcessClassifier(kernel=kernel).fit(X, y)
+    assert_array_equal(gpc.predict(X),
+                       gpc.predict_proba(X)[:, 1] >= 0.5)
+
+
 @pytest.mark.parametrize('kernel', non_fixed_kernels)
 def test_lml_improving(kernel):
     # Test that hyperparameter-tuning improves log-marginal likelihood.
     gpc = GaussianProcessClassifier(kernel=kernel).fit(X, y)
     assert (gpc.log_marginal_likelihood(gpc.kernel_.theta) >
-                   gpc.log_marginal_likelihood(kernel.theta))
+            gpc.log_marginal_likelihood(kernel.theta))
 
 
 @pytest.mark.parametrize('kernel', kernels)
@@ -139,7 +152,7 @@ def optimizer(obj_func, initial_theta, bounds):
     gpc.fit(X, y_mc)
     # Checks that optimizer improved marginal likelihood
     assert (gpc.log_marginal_likelihood(gpc.kernel_.theta) >
-                   gpc.log_marginal_likelihood(kernel.theta))
+            gpc.log_marginal_likelihood(kernel.theta))
 
 
 @pytest.mark.parametrize('kernel', kernels)
diff --git a/sklearn/gaussian_process/tests/test_gpr.py b/sklearn/gaussian_process/tests/test_gpr.py
index 64b177ce17c48..eb4bc6dec1761 100644
--- a/sklearn/gaussian_process/tests/test_gpr.py
+++ b/sklearn/gaussian_process/tests/test_gpr.py
@@ -13,6 +13,7 @@
 from sklearn.gaussian_process.kernels \
     import RBF, ConstantKernel as C, WhiteKernel
 from sklearn.gaussian_process.kernels import DotProduct
+from sklearn.gaussian_process.tests._mini_sequence_kernel import MiniSeqKernel
 
 from sklearn.utils._testing \
     import (assert_array_less,
@@ -53,12 +54,26 @@ def test_gpr_interpolation(kernel):
     assert_almost_equal(np.diag(y_cov), 0.)
 
 
+def test_gpr_interpolation_structured():
+    # Test the interpolating property for different kernels.
+    kernel = MiniSeqKernel(baseline_similarity_bounds='fixed')
+    X = ['A', 'B', 'C']
+    y = np.array([1, 2, 3])
+    gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)
+    y_pred, y_cov = gpr.predict(X, return_cov=True)
+
+    assert_almost_equal(kernel(X, eval_gradient=True)[1].ravel(),
+                        (1 - np.eye(len(X))).ravel())
+    assert_almost_equal(y_pred, y)
+    assert_almost_equal(np.diag(y_cov), 0.)
+
+
 @pytest.mark.parametrize('kernel', non_fixed_kernels)
 def test_lml_improving(kernel):
     # Test that hyperparameter-tuning improves log-marginal likelihood.
     gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)
     assert (gpr.log_marginal_likelihood(gpr.kernel_.theta) >
-                   gpr.log_marginal_likelihood(kernel.theta))
+            gpr.log_marginal_likelihood(kernel.theta))
 
 
 @pytest.mark.parametrize('kernel', kernels)
@@ -66,7 +81,7 @@ def test_lml_precomputed(kernel):
     # Test that lml of optimized kernel is stored correctly.
     gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)
     assert (gpr.log_marginal_likelihood(gpr.kernel_.theta) ==
-                 gpr.log_marginal_likelihood())
+            gpr.log_marginal_likelihood())
 
 
 @pytest.mark.parametrize('kernel', kernels)
@@ -179,7 +194,7 @@ def test_anisotropic_kernel():
     kernel = RBF([1.0, 1.0])
     gpr = GaussianProcessRegressor(kernel=kernel).fit(X, y)
     assert (np.exp(gpr.kernel_.theta[1]) >
-                   np.exp(gpr.kernel_.theta[0]) * 5)
+            np.exp(gpr.kernel_.theta[0]) * 5)
 
 
 def test_random_starts():
@@ -297,7 +312,7 @@ def optimizer(obj_func, initial_theta, bounds):
     gpr.fit(X, y)
     # Checks that optimizer improved marginal likelihood
     assert (gpr.log_marginal_likelihood(gpr.kernel_.theta) >
-                   gpr.log_marginal_likelihood(gpr.kernel.theta))
+            gpr.log_marginal_likelihood(gpr.kernel.theta))
 
 
 def test_gpr_correct_error_message():
diff --git a/sklearn/gaussian_process/tests/test_kernels.py b/sklearn/gaussian_process/tests/test_kernels.py
index ecb636d13103b..aa176f33aa50c 100644
--- a/sklearn/gaussian_process/tests/test_kernels.py
+++ b/sklearn/gaussian_process/tests/test_kernels.py
@@ -14,22 +14,22 @@
 from sklearn.gaussian_process.kernels \
     import (RBF, Matern, RationalQuadratic, ExpSineSquared, DotProduct,
             ConstantKernel, WhiteKernel, PairwiseKernel, KernelOperator,
-            Exponentiation, Kernel)
+            Exponentiation, Kernel, CompoundKernel)
 from sklearn.base import clone
 
 from sklearn.utils._testing import (assert_almost_equal, assert_array_equal,
-                                   assert_array_almost_equal,
-                                   assert_raise_message)
+                                    assert_array_almost_equal,
+                                    assert_raise_message)
 
 
 X = np.random.RandomState(0).normal(0, 1, (5, 2))
 Y = np.random.RandomState(0).normal(0, 1, (6, 2))
 
-kernel_white = RBF(length_scale=2.0) + WhiteKernel(noise_level=3.0)
+kernel_rbf_plus_white = RBF(length_scale=2.0) + WhiteKernel(noise_level=3.0)
 kernels = [RBF(length_scale=2.0), RBF(length_scale_bounds=(0.5, 2.0)),
            ConstantKernel(constant_value=10.0),
            2.0 * RBF(length_scale=0.33, length_scale_bounds="fixed"),
-           2.0 * RBF(length_scale=0.5), kernel_white,
+           2.0 * RBF(length_scale=0.5), kernel_rbf_plus_white,
            2.0 * RBF(length_scale=[0.5, 2.0]),
            2.0 * Matern(length_scale=0.33, length_scale_bounds="fixed"),
            2.0 * Matern(length_scale=0.5, nu=0.5),
@@ -92,8 +92,7 @@ def test_kernel_theta(kernel):
     # Check that values returned in theta are consistent with
     # hyperparameter values (being their logarithms)
     for i, hyperparameter in enumerate(kernel.hyperparameters):
-        assert (theta[i] ==
-                     np.log(getattr(kernel, hyperparameter.name)))
+        assert (theta[i] == np.log(getattr(kernel, hyperparameter.name)))
 
     # Fixed kernel parameters must be excluded from theta and gradient.
     for i, hyperparameter in enumerate(kernel.hyperparameters):
@@ -129,7 +128,7 @@ def test_kernel_theta(kernel):
 @pytest.mark.parametrize('kernel',
                          [kernel for kernel in kernels
                           # Identity is not satisfied on diagonal
-                          if kernel != kernel_white])
+                          if kernel != kernel_rbf_plus_white])
 def test_auto_vs_cross(kernel):
     # Auto-correlation and cross-correlation should be consistent.
     K_auto = kernel(X)
@@ -186,6 +185,27 @@ def test_kernel_stationary(kernel):
     assert_almost_equal(K[0, 0], np.diag(K))
 
 
+@pytest.mark.parametrize('kernel',  kernels)
+def test_kernel_input_type(kernel):
+    # Test whether kernels is for vectors or structured data
+    if isinstance(kernel, Exponentiation):
+        assert(kernel.requires_vector_input ==
+               kernel.kernel.requires_vector_input)
+    if isinstance(kernel, KernelOperator):
+        assert(kernel.requires_vector_input ==
+               (kernel.k1.requires_vector_input or
+                kernel.k2.requires_vector_input))
+
+
+def test_compound_kernel_input_type():
+    kernel = CompoundKernel([WhiteKernel(noise_level=3.0)])
+    assert(kernel.requires_vector_input is not True)
+
+    kernel = CompoundKernel([WhiteKernel(noise_level=3.0),
+                             RBF(length_scale=2.0)])
+    assert(kernel.requires_vector_input)
+
+
 def check_hyperparameters_equal(kernel1, kernel2):
     # Check that hyperparameters of two kernels are equal
     for attr in set(dir(kernel1) + dir(kernel2)):
@@ -236,8 +256,7 @@ def test_kernel_clone_after_set_params(kernel):
             params['length_scale_bounds'] = bounds * 2
         kernel_cloned.set_params(**params)
         kernel_cloned_clone = clone(kernel_cloned)
-        assert (kernel_cloned_clone.get_params() ==
-                     kernel_cloned.get_params())
+        assert (kernel_cloned_clone.get_params() == kernel_cloned.get_params())
         assert id(kernel_cloned_clone) != id(kernel_cloned)
         check_hyperparameters_equal(kernel_cloned, kernel_cloned_clone)
 
@@ -266,7 +285,7 @@ def test_kernel_versus_pairwise(kernel):
     # Check that GP kernels can also be used as pairwise kernels.
 
     # Test auto-kernel
-    if kernel != kernel_white:
+    if kernel != kernel_rbf_plus_white:
         # For WhiteKernel: k(X) != k(X,X). This is assumed by
         # pairwise_kernels
         K1 = kernel(X)


Code review comments:
diff_hunk: @@ -0,0 +1,172 @@
+"""
+==========================================================================
+Gaussian processes on discrete data structures
+==========================================================================
+
+This example illustrates the use of Gaussian processes for regression and
+classification tasks on data that are not in fixed-length feature vector form.
+This is achieved through the use of kernel functions that operates directly
+on discrete structures such as variable-length sequences, trees, and graphs.
+
+Specifically, here the input variables are some gene sequences stored as
+variable-length strings consisting of letters 'A', 'T', 'C', and 'G',
+while the output variables are floating point numbers and True/False labels
+in the regression and classification tasks, respectively.
+
+A kernel between the gene sequences is defined using R-convolution [1]_ by
+integrating a binary letter-wise kernel over all pairs of letters among a pair
+of strings.
+
+This example will generate three figures.
+
+In the first figure, we visualize the value of the kernel, i.e. the similarity
+of the sequences, using a colormap. Brighter color here indicates higher
+similarity.
+
+In the second figure, we show some regression result on a dataset of 6
+sequences. Here we use the 1st, 2nd, 4th, and 5th sequences as the training set
+to make predictions on the 3rd and 6th sequences.
+
+In the third figure, we demonstrate a classification model by training on 6
+sequences and make predictions on another 5 sequences. The ground truth here is
+simply  whether there is at least one 'A' in the sequence. Here the model makes
+four correct classifications and fails on one.
+
+.. [1] Haussler, D. (1999). Convolution kernels on discrete structures
+(Vol. 646). Technical report, Department of Computer Science, University of
+California at Santa Cruz.
+"""
+print(__doc__)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
+from sklearn.gaussian_process.kernels import GenericKernelMixin
+from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process import GaussianProcessClassifier
+from sklearn.base import clone
+
+
+class SequenceKernel(GenericKernelMixin, Kernel):
+    '''
+    A minimal (but valid) convolutional kernel for sequences of variable
+    lengths.'''
+    def __init__(self,
+                 baseline_similarity=0.5,
+                 baseline_similarity_bounds=(1e-5, 1)):
+        self.baseline_similarity = baseline_similarity
+        self.baseline_similarity_bounds = baseline_similarity_bounds
+
+    @property
+    def hyperparameter_baseline_similarity(self):
+        return Hyperparameter("baseline_similarity",
+                              "numeric",
+                              self.baseline_similarity_bounds)
+
+    def _f(self, s1, s2):
+        '''
+        kernel value between a pair of sequences
+        '''
+        return sum([1.0 if c1 == c2 else self.baseline_similarity
+                   for c1 in s1
+                   for c2 in s2])
+
+    def _g(self, s1, s2):
+        '''
+        kernel derivative between a pair of sequences
+        '''
+        return sum([0.0 if c1 == c2 else 1.0
+                    for c1 in s1
+                    for c2 in s2])
+
+    def __call__(self, X, Y=None, eval_gradient=False):
+        if Y is None:
+            Y = X
+
+        if eval_gradient:
+            return (np.array([[self._f(x, y) for y in Y] for x in X]),
+                    np.array([[[self._g(x, y)] for y in Y] for x in X]))
+        else:
+            return np.array([[self._f(x, y) for y in Y] for x in X])
+
+    def diag(self, X):
+        return np.array([self._f(x, x) for x in X])
+
+    def is_stationary(self):
+        return False
+
+    def clone_with_theta(self, theta):
+        cloned = clone(self)
+        cloned.theta = theta
+        return cloned
+
+
+kernel = SequenceKernel()
+
+X = np.array(['AGCT', 'AGC', 'AACT', 'TAA', 'AAA', 'GAACA'])
+y = np.array([1.0, 1.0, 2.0, 2.0, 3.0, 3.0])
+
+'''
+Sequence similarity matrix under the kernel
+===========================================
+'''
+
+K = kernel(X)
+D = kernel.diag(X)
+
+plt.figure(figsize=(8, 5))
+plt.imshow(np.diag(D**-0.5).dot(K).dot(np.diag(D**-0.5)))
+plt.xticks(np.arange(len(X)), X)
+plt.yticks(np.arange(len(X)), X)
+plt.title('Sequence similarity under the kernel')
+
+'''
+Regression
+==========
+'''
+
+training_idx = [0, 1, 3, 4]
+gp = GaussianProcessRegressor(kernel)
+gp.fit(X[training_idx], y[training_idx])
+
+plt.figure(figsize=(8, 5))
+plt.bar(np.arange(len(X)), gp.predict(X), color='b', label='prediction')
+plt.bar(training_idx, y[training_idx], width=0.2, color='r',
+        alpha=0.5, label='training')
+plt.xticks(np.arange(len(X)), X)
+plt.title('Regression on sequences')
+plt.legend()
+
+'''
+Classification
+==============
+'''
+
+X = np.array(['AGCT', 'CGA', 'TAAC', 'TCG', 'CTTT', 'TGCT'])
+# whether there are 'A's in the sequence
+t = np.array([True, True, True, False, False, False])
+
+gp = GaussianProcessClassifier(kernel)
+gp.fit(X, t)
+
+seqs_test = ['AAA', 'ATAG', 'CTC', 'CT', 'C']
path: examples/gaussian_process/plot_gpr_on_structured_data.py
text: @user1:
The variable names are no clear. I propose to rename `X` into `X_train`, `seqs_test` into `X_test`, `t` into `y_train`, and `clss_test` into `y_test`.

@author:
Fixed in 57da4d657ff5ee9ef6e580574b5ad2ed896c5916

diff_hunk: @@ -114,8 +114,10 @@ def optimizer(obj_func, initial_theta, bounds):
 
     Attributes
     ----------
-    X_train_ : array-like of shape (n_samples, n_features)
-        Feature values in training data (also required for prediction)
+    X_train_ : sequence of length n_samples
+        Feature vectors or other representations of training data (also
+        required for prediction). Could either be array-like with shape =
+        (n_samples, n_features) or a list of objects
path: sklearn/gaussian_process/_gpr.py
text: @user1:
```suggestion
        (n_samples, n_features) or a list of objects.
```

diff_hunk: @@ -403,6 +411,19 @@ def is_stationary(self):
         return True
 
 
+class GenericKernelMixin:
+    """Mixin for kernels which operate on generic objects such as variable-
+    length sequences, trees, and graphs.
+
+    .. versionadded:: 0.22
+    """
+
+    @property
+    def requires_vector_input(self):
+        """whether the kernel works only on fixed-length feature vectors."""
path: sklearn/gaussian_process/kernels.py
text: @user1:
```suggestion
        """Whether the kernel works only on fixed-length feature vectors."""
```

diff_hunk: @@ -0,0 +1,50 @@
+from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
+from sklearn.gaussian_process.kernels import GenericKernelMixin
+from sklearn.gaussian_process.kernels import StationaryKernelMixin
+import numpy as np
+from sklearn.base import clone
+
+
+class MiniSeqKernel(GenericKernelMixin,
+                    StationaryKernelMixin,
+                    Kernel):
+    '''
+    a mimimal (but valid) convolutional kernel for sequences of variable length
path: sklearn/gaussian_process/tests/_mini_sequence_kernel.py
text: @user1:
```suggestion
    A mimimal (but valid) convolutional kernel for sequences of variable length.
```

diff_hunk: @@ -186,6 +185,27 @@ def test_kernel_stationary(kernel):
     assert_almost_equal(K[0, 0], np.diag(K))
 
 
+@pytest.mark.parametrize('kernel',  kernels)
+def test_kernel_input_type(kernel):
+    # Test whether kernels is for vectors or structured data
+    if isinstance(kernel, Exponentiation):
+        assert(kernel.requires_vector_input ==
+               kernel.kernel.requires_vector_input)
+    if isinstance(kernel, KernelOperator):
+        assert(kernel.requires_vector_input ==
+               (kernel.k1.requires_vector_input or
+                kernel.k2.requires_vector_input))
+
+
+def test_compound_kernel_input_type():
+    kernel = CompoundKernel([WhiteKernel(noise_level=3.0)])
+    assert(kernel.requires_vector_input is not True)
+
+    kernel = CompoundKernel([WhiteKernel(noise_level=3.0),
+                             RBF(length_scale=2.0)])
+    assert(kernel.requires_vector_input)
path: sklearn/gaussian_process/tests/test_kernels.py
text: @user1:
```suggestion
    assert kernel.requires_vector_input
```

diff_hunk: @@ -186,6 +185,27 @@ def test_kernel_stationary(kernel):
     assert_almost_equal(K[0, 0], np.diag(K))
 
 
+@pytest.mark.parametrize('kernel',  kernels)
+def test_kernel_input_type(kernel):
+    # Test whether kernels is for vectors or structured data
+    if isinstance(kernel, Exponentiation):
+        assert(kernel.requires_vector_input ==
+               kernel.kernel.requires_vector_input)
+    if isinstance(kernel, KernelOperator):
+        assert(kernel.requires_vector_input ==
+               (kernel.k1.requires_vector_input or
+                kernel.k2.requires_vector_input))
+
+
+def test_compound_kernel_input_type():
+    kernel = CompoundKernel([WhiteKernel(noise_level=3.0)])
+    assert(kernel.requires_vector_input is not True)
path: sklearn/gaussian_process/tests/test_kernels.py
text: @user1:
```suggestion
    assert not kernel.requires_vector_input
```

diff_hunk: @@ -0,0 +1,172 @@
+"""
+==========================================================================
+Gaussian processes on discrete data structures
+==========================================================================
+
+This example illustrates the use of Gaussian processes for regression and
+classification tasks on data that are not in fixed-length feature vector form.
+This is achieved through the use of kernel functions that operates directly
+on discrete structures such as variable-length sequences, trees, and graphs.
+
+Specifically, here the input variables are some gene sequences stored as
+variable-length strings consisting of letters 'A', 'T', 'C', and 'G',
+while the output variables are floating point numbers and True/False labels
+in the regression and classification tasks, respectively.
+
+A kernel between the gene sequences is defined using R-convolution [1]_ by
+integrating a binary letter-wise kernel over all pairs of letters among a pair
+of strings.
+
+This example will generate three figures.
+
+In the first figure, we visualize the value of the kernel, i.e. the similarity
+of the sequences, using a colormap. Brighter color here indicates higher
+similarity.
+
+In the second figure, we show some regression result on a dataset of 6
+sequences. Here we use the 1st, 2nd, 4th, and 5th sequences as the training set
+to make predictions on the 3rd and 6th sequences.
+
+In the third figure, we demonstrate a classification model by training on 6
+sequences and make predictions on another 5 sequences. The ground truth here is
+simply  whether there is at least one 'A' in the sequence. Here the model makes
+four correct classifications and fails on one.
+
+.. [1] Haussler, D. (1999). Convolution kernels on discrete structures
+(Vol. 646). Technical report, Department of Computer Science, University of
+California at Santa Cruz.
+"""
+print(__doc__)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
+from sklearn.gaussian_process.kernels import GenericKernelMixin
+from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process import GaussianProcessClassifier
+from sklearn.base import clone
+
+
+class SequenceKernel(GenericKernelMixin, Kernel):
+    '''
+    A minimal (but valid) convolutional kernel for sequences of variable
+    lengths.'''
+    def __init__(self,
+                 baseline_similarity=0.5,
+                 baseline_similarity_bounds=(1e-5, 1)):
+        self.baseline_similarity = baseline_similarity
+        self.baseline_similarity_bounds = baseline_similarity_bounds
+
+    @property
+    def hyperparameter_baseline_similarity(self):
+        return Hyperparameter("baseline_similarity",
+                              "numeric",
+                              self.baseline_similarity_bounds)
+
+    def _f(self, s1, s2):
+        '''
+        kernel value between a pair of sequences
+        '''
+        return sum([1.0 if c1 == c2 else self.baseline_similarity
+                   for c1 in s1
+                   for c2 in s2])
+
+    def _g(self, s1, s2):
+        '''
+        kernel derivative between a pair of sequences
+        '''
+        return sum([0.0 if c1 == c2 else 1.0
+                    for c1 in s1
+                    for c2 in s2])
+
+    def __call__(self, X, Y=None, eval_gradient=False):
+        if Y is None:
+            Y = X
+
+        if eval_gradient:
+            return (np.array([[self._f(x, y) for y in Y] for x in X]),
+                    np.array([[[self._g(x, y)] for y in Y] for x in X]))
+        else:
+            return np.array([[self._f(x, y) for y in Y] for x in X])
+
+    def diag(self, X):
+        return np.array([self._f(x, x) for x in X])
+
+    def is_stationary(self):
+        return False
+
+    def clone_with_theta(self, theta):
+        cloned = clone(self)
+        cloned.theta = theta
+        return cloned
+
+
+kernel = SequenceKernel()
+
+X = np.array(['AGCT', 'AGC', 'AACT', 'TAA', 'AAA', 'GAACA'])
+y = np.array([1.0, 1.0, 2.0, 2.0, 3.0, 3.0])
+
+'''
+Sequence similarity matrix under the kernel
+===========================================
+'''
+
+K = kernel(X)
+D = kernel.diag(X)
+
+plt.figure(figsize=(8, 5))
+plt.imshow(np.diag(D**-0.5).dot(K).dot(np.diag(D**-0.5)))
+plt.xticks(np.arange(len(X)), X)
+plt.yticks(np.arange(len(X)), X)
+plt.title('Sequence similarity under the kernel')
+
+'''
+Regression
+==========
+'''
+
+training_idx = [0, 1, 3, 4]
+gp = GaussianProcessRegressor(kernel)
+gp.fit(X[training_idx], y[training_idx])
+
+plt.figure(figsize=(8, 5))
+plt.bar(np.arange(len(X)), gp.predict(X), color='b', label='prediction')
+plt.bar(training_idx, y[training_idx], width=0.2, color='r',
+        alpha=0.5, label='training')
path: examples/gaussian_process/plot_gpr_on_structured_data.py
text: @user1:
nitpick: with `alpha=0.5`, the legend is not the same color as the plot. I would use `alpha=1`.
```suggestion
        alpha=1, label='training')
```

diff_hunk: @@ -0,0 +1,172 @@
+"""
+==========================================================================
+Gaussian processes on discrete data structures
+==========================================================================
+
+This example illustrates the use of Gaussian processes for regression and
+classification tasks on data that are not in fixed-length feature vector form.
+This is achieved through the use of kernel functions that operates directly
+on discrete structures such as variable-length sequences, trees, and graphs.
+
+Specifically, here the input variables are some gene sequences stored as
+variable-length strings consisting of letters 'A', 'T', 'C', and 'G',
+while the output variables are floating point numbers and True/False labels
+in the regression and classification tasks, respectively.
+
+A kernel between the gene sequences is defined using R-convolution [1]_ by
+integrating a binary letter-wise kernel over all pairs of letters among a pair
+of strings.
+
+This example will generate three figures.
+
+In the first figure, we visualize the value of the kernel, i.e. the similarity
+of the sequences, using a colormap. Brighter color here indicates higher
+similarity.
+
+In the second figure, we show some regression result on a dataset of 6
+sequences. Here we use the 1st, 2nd, 4th, and 5th sequences as the training set
+to make predictions on the 3rd and 6th sequences.
+
+In the third figure, we demonstrate a classification model by training on 6
+sequences and make predictions on another 5 sequences. The ground truth here is
+simply  whether there is at least one 'A' in the sequence. Here the model makes
+four correct classifications and fails on one.
+
+.. [1] Haussler, D. (1999). Convolution kernels on discrete structures
+(Vol. 646). Technical report, Department of Computer Science, University of
+California at Santa Cruz.
+"""
+print(__doc__)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
+from sklearn.gaussian_process.kernels import GenericKernelMixin
+from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process import GaussianProcessClassifier
+from sklearn.base import clone
+
+
+class SequenceKernel(GenericKernelMixin, Kernel):
+    '''
+    A minimal (but valid) convolutional kernel for sequences of variable
+    lengths.'''
+    def __init__(self,
+                 baseline_similarity=0.5,
+                 baseline_similarity_bounds=(1e-5, 1)):
+        self.baseline_similarity = baseline_similarity
+        self.baseline_similarity_bounds = baseline_similarity_bounds
+
+    @property
+    def hyperparameter_baseline_similarity(self):
+        return Hyperparameter("baseline_similarity",
+                              "numeric",
+                              self.baseline_similarity_bounds)
+
+    def _f(self, s1, s2):
+        '''
+        kernel value between a pair of sequences
+        '''
+        return sum([1.0 if c1 == c2 else self.baseline_similarity
+                   for c1 in s1
+                   for c2 in s2])
+
+    def _g(self, s1, s2):
+        '''
+        kernel derivative between a pair of sequences
+        '''
+        return sum([0.0 if c1 == c2 else 1.0
+                    for c1 in s1
+                    for c2 in s2])
+
+    def __call__(self, X, Y=None, eval_gradient=False):
+        if Y is None:
+            Y = X
+
+        if eval_gradient:
+            return (np.array([[self._f(x, y) for y in Y] for x in X]),
+                    np.array([[[self._g(x, y)] for y in Y] for x in X]))
+        else:
+            return np.array([[self._f(x, y) for y in Y] for x in X])
+
+    def diag(self, X):
+        return np.array([self._f(x, x) for x in X])
+
+    def is_stationary(self):
+        return False
+
+    def clone_with_theta(self, theta):
+        cloned = clone(self)
+        cloned.theta = theta
+        return cloned
+
+
+kernel = SequenceKernel()
+
+X = np.array(['AGCT', 'AGC', 'AACT', 'TAA', 'AAA', 'GAACA'])
+y = np.array([1.0, 1.0, 2.0, 2.0, 3.0, 3.0])
+
+'''
+Sequence similarity matrix under the kernel
+===========================================
+'''
+
+K = kernel(X)
+D = kernel.diag(X)
+
+plt.figure(figsize=(8, 5))
+plt.imshow(np.diag(D**-0.5).dot(K).dot(np.diag(D**-0.5)))
+plt.xticks(np.arange(len(X)), X)
+plt.yticks(np.arange(len(X)), X)
+plt.title('Sequence similarity under the kernel')
+
+'''
+Regression
+==========
+'''
+
+training_idx = [0, 1, 3, 4]
+gp = GaussianProcessRegressor(kernel)
path: examples/gaussian_process/plot_gpr_on_structured_data.py
text: @user1:
```suggestion
gp = GaussianProcessRegressor(kernel=kernel)
```

diff_hunk: @@ -383,8 +383,9 @@ def score(self, X, y, sample_weight=None):
         ----------
         X : array-like of shape (n_samples, n_features)
             Test samples. For some estimators this may be a
-            precomputed kernel matrix instead, shape = (n_samples,
-            n_samples_fitted], where n_samples_fitted is the number of
+            precomputed kernel matrix or a list of generic objects instead,
+            shape = (n_samples, n_samples_fitted],
path: sklearn/base.py
text: @user1:
```suggestion
            shape = (n_samples, n_samples_fitted),
```

diff_hunk: @@ -115,8 +115,10 @@ def optimizer(obj_func, initial_theta, bounds):
 
     Attributes
     ----------
-    X_train_ : array-like of shape (n_samples, n_features)
-        Feature values in training data (also required for prediction)
+    X_train_ : sequence of length n_samples
+        Feature vectors or other representations of training data (also
+        required for prediction). Could either be array-like with shape =
+        (n_samples, n_features) or a list of objects
path: sklearn/gaussian_process/_gpc.py
text: @user1:
```suggestion
        (n_samples, n_features) or a list of objects.
```