Cacher is the code snippet organizer for pro developers

import GPy
import numpy as np


class HeteroNoise(GPy.kern.Kern):
    """
    heteroscedastic white noise as function of X
    """

    def __init__(self, input_dim, variance_function=1., active_dims=None):
        """
        interpolating noise function:
            variance_function=lambda X:np.interp(X[:, 0], x, yerr)
        """
        super(HeteroNoise, self).__init__(input_dim, active_dims, name='noise')

        if callable(variance_function):
            self.variance_function = variance_function
        if np.isscalar(variance_function):
            self.variance_function = lambda X: np.full_like(X[:, 0], variance_function)
        else:
            raise TypeError("'variance_function' should be a function!")

    def K(self, X, X2=None):
        if X2 is None:
            return np.diag(self.Kdiag(X))
        else:
            return np.zeros((X.shape[0], X2.shape[0]))

    def Kdiag(self, X):
        return self.variance_function(X)

import GPy
import numpy as np


class Noise(GPy.kern.Kern):
    def __init__(self, input_dim, variance=1., active_dims=None):
        super(Noise, self).__init__(input_dim, active_dims, name='noise')
        assert input_dim == 1, "For this kernel we assume input_dim=1"

        from paramz.transformations import Logexp
        self.variance = GPy.Param('variance', variance, Logexp())
        self.link_parameters(self.variance)

    def K(self, X, X2=None):
        if X2 is None:
            return np.diag(self.Kdiag(X)) * self.variance
        else:
            return np.zeros((X.shape[0], X2.shape[0]))

    def Kdiag(self, X):
        return np.interp(X[:, 0], gmag, cvar) * self.variance

    def update_gradients_full(self, dL_dK, X, X2):
        if X2 is None:
            dKdvar = np.diag(self.Kdiag(X))
            self.variance.gradient = np.sum(dKdvar * dL_dK)
        else:
            self.variance.gradient = 0


def rbf_kern_multiquadric(r, epsilon):
    return np.sqrt((r / epsilon)**2 + 1)


def rbf_norm(x1, x2):
    """
    x1, x2: shape (nsample, nfeature)

    >>> rbf_norm(rand(10, 2), rand(5, 2)).shape
    # (10, 5)
    """
    return np.sqrt(((x1[..., None, :] - x2)**2).sum(axis=-1))


def rbf_eval(xi, x, w, epsilon=None, kern=rbf_kern_multiquadric):
    """
    xi: shape (nsample, nfeature)
    x: shape (nnode, nfeature)
        nodes
    w: shape (nnode,)
        weights
    """
    if epsilon is None:
        epsilon = (x.max() - x.min()) / len(x)
    r = rbf_norm(xi, x)
    A = kern(r, epsilon)
    return (A * w).sum(axis=-1)


class RbfNoise(GPy.kern.Kern):
    def __init__(self, input_dim, x, w, p, rbf_kern=rbf_kern_multiquadric, active_dims=None):
        """
        x: rbf nodes
        w: weights of nodes
        p: fiducial
        """
        super(RbfNoise, self).__init__(input_dim, active_dims, name='rbfnoise')

        self.epsilon = (x.max() - x.min()) / len(x)
        self.rbf_kern = rbf_kern
        self.x = x
        self.w = GPy.Param('w', w)
        self.p = GPy.Param('p', p)
        self.link_parameters(self.w, self.p)

        if not np.isscalar(p):
            raise NotImplementedError

    def K(self, X, X2=None):
        if X2 is None:
            return np.diag(self.Kdiag(X))
        else:
            return np.zeros((X.shape[0], X2.shape[0]))

    def lnKdiag(self, X):
        lnKdiag = rbf_eval(X, self.x, self.w, epsilon=self.epsilon, kern=self.rbf_kern) + self.p
        return lnKdiag

    def Kdiag(self, X):
        lnKdiag = self.lnKdiag(X)
        return np.exp(lnKdiag)

    def update_gradients_full(self, dL_dK, X, X2):
        if X2 is None:
            r = rbf_norm(X, self.x)
            A = self.rbf_kern(r, self.epsilon)
            lnKdiag = (A * self.w).sum(axis=-1)
            Kdiag = np.exp(lnKdiag)
            dL_dKdiag = np.diagonal(dL_dK)

            self.w.gradient = np.sum(dL_dKdiag * Kdiag * A.T, axis=-1)
            self.p.gradient = np.sum(dL_dKdiag * Kdiag)
        else:
            self.w.gradient = np.zeros_like(self.w)
            self.p.gradient = 0