r"""
        Get the gradients of this basis w.r.t.\ the mean and length scales.

        Parameters
        ----------
        x: ndarray
            (n, d) array of observations where n is the number of samples, and
            d is the dimensionality of x.
        mean: ndarray, optional
            array of shape (d,) frequency means (one for each dimension of X).
            If not input, this uses the value of the initial mean.
        lenscale: ndarray, optional
            array of shape (d,) length scales (one for each dimension of X). If
            not input, this uses the value of the initial length scale.

        Returns
        -------
        ndarray:
            shape (n, 4*nbases) where nbases is number of random rbf bases,
            again to the nearest larger two power. This is
            :math:`\partial \phi(\mathbf{x}) / \partial \boldsymbol\mu`
        ndarray:
            shape (n, 4*nbases) where nbases is number of random rbf bases,
            again to the nearest larger two power. This is
            :math:`\partial \phi(\mathbf{x}) / \partial \mathbf{l}`

def grad(self, X, mean=None, lenscale=None):
        r"""
        Get the gradients of this basis w.r.t.\ the mean and length scales.

        Parameters
        ----------
        x: ndarray
            (n, d) array of observations where n is the number of samples, and
            d is the dimensionality of x.
        mean: ndarray, optional
            array of shape (d,) frequency means (one for each dimension of X).
            If not input, this uses the value of the initial mean.
        lenscale: ndarray, optional
            array of shape (d,) length scales (one for each dimension of X). If
            not input, this uses the value of the initial length scale.

        Returns
        -------
        ndarray:
            shape (n, 4*nbases) where nbases is number of random rbf bases,
            again to the nearest larger two power. This is
            :math:`\partial \phi(\mathbf{x}) / \partial \boldsymbol\mu`
        ndarray:
            shape (n, 4*nbases) where nbases is number of random rbf bases,
            again to the nearest larger two power. This is
            :math:`\partial \phi(\mathbf{x}) / \partial \mathbf{l}`
        """
        d = X.shape[1]
        mean = self._check_dim(d, mean, paramind=0)
        lenscale = self._check_dim(d, lenscale, paramind=1)

        VX = self._makeVX(X / lenscale)
        mX = X.dot(mean)[:, np.newaxis]

        sinVXpmX = - np.sin(VX + mX)
        sinVXmmX = - np.sin(VX - mX)
        cosVXpmX = np.cos(VX + mX)
        cosVXmmX = np.cos(VX - mX)

        dPhi_len = []
        dPhi_mean = []
        for i, l in enumerate(lenscale):

            # Means
            dmX = X[:, [i]]
            dPhi_mean.append(np.hstack((dmX * sinVXpmX, dmX * cosVXpmX,
                                        -dmX * sinVXmmX, -dmX * cosVXmmX)) /
                             np.sqrt(2 * self.n))

            # Lenscales
            indlen = np.zeros(d)
            indlen[i] = 1. / l**2
            dVX = - self._makeVX(X * indlen)  # FIXME make this more efficient?
            dPhi_len.append(np.hstack((dVX * sinVXpmX, dVX * cosVXpmX,
                                       dVX * sinVXmmX, dVX * cosVXmmX)) /
                            np.sqrt(2 * self.n))

        dPhi_mean = np.dstack(dPhi_mean) if d != 1 else dPhi_mean[0]
        dPhi_len = np.dstack(dPhi_len) if d != 1 else dPhi_len[0]
        return dPhi_mean, dPhi_len