diff --git a/GPy/util/Tango.py b/GPy/util/Tango.py
new file mode 100644
index 00000000..fec0a3d1
--- /dev/null
+++ b/GPy/util/Tango.py
@@ -0,0 +1,166 @@
+import matplotlib as mpl
+
+import pylab as pb
+import sys
+#sys.path.append('/home/james/mlprojects/sitran_cluster/')
+#from switch_pylab_backend import *
+
+
+#this stuff isn;t really Tango related: maybe it could be moved out? TODO
+def removeRightTicks(ax=None):
+    ax = ax or pb.gca()
+    for i, line in enumerate(ax.get_yticklines()):
+        if i%2 == 1:   # odd indices
+            line.set_visible(False) 
+def removeUpperTicks(ax=None):
+    ax = ax or pb.gca()
+    for i, line in enumerate(ax.get_xticklines()):
+        if i%2 == 1:   # odd indices
+            line.set_visible(False) 
+def fewerXticks(ax=None,divideby=2):
+    ax = ax or pb.gca()
+    ax.set_xticks(ax.get_xticks()[::divideby])
+
+
+coloursHex = {\
+"Aluminium6":"#2e3436",\
+"Aluminium5":"#555753",\
+"Aluminium4":"#888a85",\
+"Aluminium3":"#babdb6",\
+"Aluminium2":"#d3d7cf",\
+"Aluminium1":"#eeeeec",\
+"lightPurple":"#ad7fa8",\
+"mediumPurple":"#75507b",\
+"darkPurple":"#5c3566",\
+"lightBlue":"#729fcf",\
+"mediumBlue":"#3465a4",\
+"darkBlue": "#204a87",\
+"lightGreen":"#8ae234",\
+"mediumGreen":"#73d216",\
+"darkGreen":"#4e9a06",\
+"lightChocolate":"#e9b96e",\
+"mediumChocolate":"#c17d11",\
+"darkChocolate":"#8f5902",\
+"lightRed":"#ef2929",\
+"mediumRed":"#cc0000",\
+"darkRed":"#a40000",\
+"lightOrange":"#fcaf3e",\
+"mediumOrange":"#f57900",\
+"darkOrange":"#ce5c00",\
+"lightButter":"#fce94f",\
+"mediumButter":"#edd400",\
+"darkButter":"#c4a000"}
+
+darkList = [coloursHex['darkBlue'],coloursHex['darkRed'],coloursHex['darkGreen'], coloursHex['darkOrange'], coloursHex['darkButter'], coloursHex['darkPurple'], coloursHex['darkChocolate'], coloursHex['Aluminium6']]
+mediumList = [coloursHex['mediumBlue'], coloursHex['mediumRed'],coloursHex['mediumGreen'], coloursHex['mediumOrange'], coloursHex['mediumButter'], coloursHex['mediumPurple'], coloursHex['mediumChocolate'], coloursHex['Aluminium5']]
+lightList = [coloursHex['lightBlue'], coloursHex['lightRed'],coloursHex['lightGreen'], coloursHex['lightOrange'], coloursHex['lightButter'], coloursHex['lightPurple'], coloursHex['lightChocolate'], coloursHex['Aluminium4']]
+
+def currentDark():
+    return darkList[-1]
+def currentMedium():
+    return mediumList[-1]
+def currentLight():
+    return lightList[-1]
+
+def nextDark():
+    darkList.append(darkList.pop(0))
+    return darkList[-1]
+def nextMedium():
+    mediumList.append(mediumList.pop(0))
+    return mediumList[-1]
+def nextLight():
+    lightList.append(lightList.pop(0))
+    return lightList[-1]
+
+def reset():
+    while not darkList[0]==coloursHex['darkBlue']:
+        darkList.append(darkList.pop(0))
+    while not mediumList[0]==coloursHex['mediumBlue']:
+        mediumList.append(mediumList.pop(0))
+    while not lightList[0]==coloursHex['lightBlue']:
+        lightList.append(lightList.pop(0))
+
+def setLightFigures():
+    mpl.rcParams['axes.edgecolor']=coloursHex['Aluminium6']
+    mpl.rcParams['axes.facecolor']=coloursHex['Aluminium2']
+    mpl.rcParams['axes.labelcolor']=coloursHex['Aluminium6']
+    mpl.rcParams['figure.edgecolor']=coloursHex['Aluminium6']
+    mpl.rcParams['figure.facecolor']=coloursHex['Aluminium2']
+    mpl.rcParams['grid.color']=coloursHex['Aluminium6']
+    mpl.rcParams['savefig.edgecolor']=coloursHex['Aluminium2']
+    mpl.rcParams['savefig.facecolor']=coloursHex['Aluminium2']
+    mpl.rcParams['text.color']=coloursHex['Aluminium6']
+    mpl.rcParams['xtick.color']=coloursHex['Aluminium6']
+    mpl.rcParams['ytick.color']=coloursHex['Aluminium6']
+
+def setDarkFigures():
+    mpl.rcParams['axes.edgecolor']=coloursHex['Aluminium2']
+    mpl.rcParams['axes.facecolor']=coloursHex['Aluminium6']
+    mpl.rcParams['axes.labelcolor']=coloursHex['Aluminium2']
+    mpl.rcParams['figure.edgecolor']=coloursHex['Aluminium2']
+    mpl.rcParams['figure.facecolor']=coloursHex['Aluminium6']
+    mpl.rcParams['grid.color']=coloursHex['Aluminium2']
+    mpl.rcParams['savefig.edgecolor']=coloursHex['Aluminium6']
+    mpl.rcParams['savefig.facecolor']=coloursHex['Aluminium6']
+    mpl.rcParams['text.color']=coloursHex['Aluminium2']
+    mpl.rcParams['xtick.color']=coloursHex['Aluminium2']
+    mpl.rcParams['ytick.color']=coloursHex['Aluminium2']
+
+def hex2rgb(hexcolor):
+    hexcolor = [hexcolor[1+2*i:1+2*(i+1)] for i in range(3)]
+    r,g,b = [int(n,16) for n in hexcolor]
+    return (r,g,b)
+
+coloursRGB = dict([(k,hex2rgb(i)) for k,i in coloursHex.items()])
+
+cdict_RB = {'red' :((0.,coloursRGB['mediumRed'][0]/256.,coloursRGB['mediumRed'][0]/256.),
+                     (.5,coloursRGB['mediumPurple'][0]/256.,coloursRGB['mediumPurple'][0]/256.),
+                     (1.,coloursRGB['mediumBlue'][0]/256.,coloursRGB['mediumBlue'][0]/256.)),
+            'green':((0.,coloursRGB['mediumRed'][1]/256.,coloursRGB['mediumRed'][1]/256.),
+                     (.5,coloursRGB['mediumPurple'][1]/256.,coloursRGB['mediumPurple'][1]/256.),
+                     (1.,coloursRGB['mediumBlue'][1]/256.,coloursRGB['mediumBlue'][1]/256.)),
+            'blue':((0.,coloursRGB['mediumRed'][2]/256.,coloursRGB['mediumRed'][2]/256.),
+                      (.5,coloursRGB['mediumPurple'][2]/256.,coloursRGB['mediumPurple'][2]/256.),
+                      (1.,coloursRGB['mediumBlue'][2]/256.,coloursRGB['mediumBlue'][2]/256.))}
+cmap_RB = mpl.colors.LinearSegmentedColormap('TangoRedBlue',cdict_RB,256)
+
+
+cdict_BGR = {'red' :((0.,coloursRGB['mediumBlue'][0]/256.,coloursRGB['mediumBlue'][0]/256.),
+                     (.5,coloursRGB['mediumGreen'][0]/256.,coloursRGB['mediumGreen'][0]/256.),
+                     (1.,coloursRGB['mediumRed'][0]/256.,coloursRGB['mediumRed'][0]/256.)),
+            'green':((0.,coloursRGB['mediumBlue'][1]/256.,coloursRGB['mediumBlue'][1]/256.),
+                     (.5,coloursRGB['mediumGreen'][1]/256.,coloursRGB['mediumGreen'][1]/256.),
+                     (1.,coloursRGB['mediumRed'][1]/256.,coloursRGB['mediumRed'][1]/256.)),
+            'blue':((0.,coloursRGB['mediumBlue'][2]/256.,coloursRGB['mediumBlue'][2]/256.),
+                      (.5,coloursRGB['mediumGreen'][2]/256.,coloursRGB['mediumGreen'][2]/256.),
+                      (1.,coloursRGB['mediumRed'][2]/256.,coloursRGB['mediumRed'][2]/256.))}
+cmap_BGR = mpl.colors.LinearSegmentedColormap('TangoRedBlue',cdict_BGR,256)
+
+cdict_Alu = {'red' :((0./5,coloursRGB['Aluminium1'][0]/256.,coloursRGB['Aluminium1'][0]/256.),
+                     (1./5,coloursRGB['Aluminium2'][0]/256.,coloursRGB['Aluminium2'][0]/256.),
+                     (2./5,coloursRGB['Aluminium3'][0]/256.,coloursRGB['Aluminium3'][0]/256.),
+                     (3./5,coloursRGB['Aluminium4'][0]/256.,coloursRGB['Aluminium4'][0]/256.),
+                     (4./5,coloursRGB['Aluminium5'][0]/256.,coloursRGB['Aluminium5'][0]/256.),
+                     (5./5,coloursRGB['Aluminium6'][0]/256.,coloursRGB['Aluminium6'][0]/256.)),
+           'green' :((0./5,coloursRGB['Aluminium1'][1]/256.,coloursRGB['Aluminium1'][1]/256.),
+                     (1./5,coloursRGB['Aluminium2'][1]/256.,coloursRGB['Aluminium2'][1]/256.),
+                     (2./5,coloursRGB['Aluminium3'][1]/256.,coloursRGB['Aluminium3'][1]/256.),
+                     (3./5,coloursRGB['Aluminium4'][1]/256.,coloursRGB['Aluminium4'][1]/256.),
+                     (4./5,coloursRGB['Aluminium5'][1]/256.,coloursRGB['Aluminium5'][1]/256.),
+                     (5./5,coloursRGB['Aluminium6'][1]/256.,coloursRGB['Aluminium6'][1]/256.)),
+            'blue' :((0./5,coloursRGB['Aluminium1'][2]/256.,coloursRGB['Aluminium1'][2]/256.),
+                     (1./5,coloursRGB['Aluminium2'][2]/256.,coloursRGB['Aluminium2'][2]/256.),
+                     (2./5,coloursRGB['Aluminium3'][2]/256.,coloursRGB['Aluminium3'][2]/256.),
+                     (3./5,coloursRGB['Aluminium4'][2]/256.,coloursRGB['Aluminium4'][2]/256.),
+                     (4./5,coloursRGB['Aluminium5'][2]/256.,coloursRGB['Aluminium5'][2]/256.),
+                     (5./5,coloursRGB['Aluminium6'][2]/256.,coloursRGB['Aluminium6'][2]/256.))}
+cmap_Alu = mpl.colors.LinearSegmentedColormap('TangoAluminium',cdict_Alu,256)
+
+if __name__=='__main__':
+    import pylab as pb
+    pb.figure()
+    pb.pcolor(pb.rand(10,10),cmap=cmap_RB)
+    pb.colorbar()
+    pb.show()
+
+    
diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py
new file mode 100644
index 00000000..9e495be2
--- /dev/null
+++ b/GPy/util/__init__.py
@@ -0,0 +1,7 @@
+import linalg
+import misc
+import plot
+import squashers
+import Tango
+import misc
+import warping_functions
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
new file mode 100644
index 00000000..3e51b7df
--- /dev/null
+++ b/GPy/util/linalg.py
@@ -0,0 +1,150 @@
+import numpy as np
+from scipy import linalg, optimize
+import pylab as pb
+import Tango
+import sys
+import re
+import pdb
+import cPickle
+import types
+import scipy.lib.lapack.flapack
+import scipy as sp
+
+def mdot(*args):
+   """Multiply all the arguments using matrix product rules.
+   The output is equivalent to multiplying the arguments one by one
+   from left to right using dot().
+   Precedence can be controlled by creating tuples of arguments,
+   for instance mdot(a,((b,c),d)) multiplies a (a*((b*c)*d)).
+   Note that this means the output of dot(a,b) and mdot(a,b) will differ if
+   a or b is a pure tuple of numbers.
+   """
+   if len(args)==1:
+       return args[0]
+   elif len(args)==2:
+       return _mdot_r(args[0],args[1])
+   else:
+       return _mdot_r(args[:-1],args[-1])
+
+def _mdot_r(a,b):
+   """Recursive helper for mdot"""
+   if type(a)==types.TupleType:
+       if len(a)>1:
+           a = mdot(*a)
+       else:
+           a = a[0]
+   if type(b)==types.TupleType:
+       if len(b)>1:
+           b = mdot(*b)
+       else:
+           b = b[0]
+   return np.dot(a,b)
+
+def jitchol(A,maxtries=5):
+    """
+    Arguments
+    ---------
+    A : An almost pd square matrix
+
+    Returns
+    -------
+    cholesky(K)
+
+    Notes
+    -----
+    Adds jitter to K, to enforce positive-definiteness
+    if stuff breaks, please check:
+    np.allclose(sp.linalg.cholesky(XXT, lower = True), np.triu(sp.linalg.cho_factor(XXT)[0]).T)
+    """
+    try:
+        return linalg.cholesky(A, lower = True)
+    except linalg.LinAlgError:
+        diagA = np.diag(A)
+        if np.any(diagA<0.):
+            raise linalg.LinAlgError, "not pd: negative diagonal elements"
+        jitter= diagA.mean()*1e-6
+        for i in range(1,maxtries+1):
+            try:
+                print 'Warning: adding jitter of '+str(jitter)
+                return linalg.cholesky(A+np.eye(A.shape[0])*jitter, lower = True)
+            except:
+                jitter *= 10
+
+        raise linalg.LinAlgError,"not positive definite, even with jitter."
+
+def pdinv(A):
+    """
+    Arguments
+    ---------
+    :param A: A DxD pd numpy array
+
+    Returns
+    -------
+    inv : the inverse of A
+    hld: 0.5* the log of the determinant of A
+    """
+    L = jitchol(A)
+    hld = np.sum(np.log(np.diag(L)))
+
+    inv = sp.lib.lapack.flapack.dpotri(L)[0]
+    # inv = linalg.flapack.dpotri(L,lower = 1)[0]
+    inv = np.tril(inv)+np.tril(inv,-1).T
+
+    return inv, hld
+
+
+def chol_inv(L):
+    """
+    Inverts a Cholesky lower triangular matrix
+
+    :param L: lower triangular matrix
+    :rtype: inverse of L
+
+    """
+
+    return linalg.flapack.dtrtri(L, lower = True)[0]
+
+
+def multiple_pdinv(A):
+    """
+    Arguments
+    ---------
+    :param A: A DxDxN numpy array (each A[:,:,i] is pd)
+
+    Returns
+    -------
+    invs : the inverses of A
+    hld: 0.5* the log of the determinants of A
+    """
+    N = A.shape[-1]
+    chols = [jitchol(A[:,:,i]) for i in range(N)]
+    halflogdets = [np.sum(np.log(np.diag(L[0]))) for L in chols]
+    invs = [linalg.flapack.dpotri(L[0],True)[0] for L in chols]
+    invs = [np.triu(I)+np.triu(I,1).T for I in invs]
+    return np.dstack(invs),np.array(halflogdets)
+
+
+def PCA(Y, Q):
+    """
+    Principal component analysis: maximum likelihood solution by SVD
+
+    Arguments
+    ---------
+    :param Y: NxD np.array of data
+    :param Q: int, dimension of projection
+
+    Returns
+    -------
+    X - NxQ np.array of dimensionality reduced data
+    W - QxD mapping from X to Y
+    """
+    if not np.allclose(Y.mean(axis=0), 0.0):
+        print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"
+        Y -= Y.mean(axis=0)
+
+    Z = linalg.svd(Y, full_matrices = False)
+    [X, W] = [Z[0][:,0:Q], np.dot(np.diag(Z[1]), Z[2]).T[:,0:Q]]
+    v = X.std(axis=0)
+    X /= v;
+    W *= v;
+    return X, W.T
diff --git a/GPy/util/misc.py b/GPy/util/misc.py
new file mode 100644
index 00000000..0725a973
--- /dev/null
+++ b/GPy/util/misc.py
@@ -0,0 +1,56 @@
+import numpy as np
+
+def linear_grid(D, n = 100, min_max = (-100, 100)):
+    """
+    Creates a D-dimensional grid of n linearly spaced points
+
+    Parameters:
+
+    D:        dimension of the grid
+    n:        number of points
+    min_max:  (min, max) list
+
+
+    """
+
+    g = np.linspace(min_max[0], min_max[1], n)
+    G = np.ones((n, D))
+
+    return G*g[:,None]
+
+def kmm_init(X, m = 10):
+    """
+    This is the same initialization algorithm that is used
+    in Kmeans++. It's quite simple and very useful to initialize
+    the locations of the inducing points in sparse GPs.
+    
+    :param X: data
+    :param m: number of inducing points
+    """
+
+    # compute the distances
+    XXT = np.dot(X, X.T)
+    D = (-2.*XXT + np.diag(XXT)[:,np.newaxis] + np.diag(XXT)[np.newaxis,:])
+
+    # select the first point
+    s = np.random.permutation(X.shape[0])[0]
+    inducing = [s]
+    prob = D[s]/D[s].sum()
+
+    for z in range(m-1):
+        s = np.random.multinomial(1, prob.flatten()).argmax()
+        inducing.append(s)
+        prob = D[s]/D[s].sum()
+
+    inducing = np.array(inducing)
+    return X[inducing]
+
+if __name__ == '__main__':
+    import pylab as plt
+    X = np.linspace(1,10, 100)[:, None]
+    X = X[np.random.permutation(X.shape[0])[:20]]
+    inducing = kmm_init(X, m = 5)
+    plt.figure()
+    plt.plot(X.flatten(), np.ones((X.shape[0],)), 'x')
+    plt.plot(inducing, 0.5* np.ones((len(inducing),)), 'o')
+    plt.ylim((0.0, 10.0))
diff --git a/GPy/util/plot.py b/GPy/util/plot.py
new file mode 100644
index 00000000..3bd7d8a0
--- /dev/null
+++ b/GPy/util/plot.py
@@ -0,0 +1,73 @@
+import Tango
+import pylab as pb
+import numpy as np
+
+def gpplot(x,mu,var,edgecol=Tango.coloursHex['darkBlue'],fillcol=Tango.coloursHex['lightBlue'],axes=None,**kwargs):
+    if axes is None:
+        axes = pb.gca()
+    mu = mu.flatten()
+    x = x.flatten()
+
+    #here's the mean
+    axes.plot(x,mu,color=edgecol,linewidth=2)
+
+    #ensure variance is a vector
+    if len(var.shape)>1:
+        err = 2*np.sqrt(np.diag(var))
+    else:
+        err = 2*np.sqrt(var)
+
+    #here's the 2*std box
+    kwargs['linewidth']=0.5
+    if not 'alpha' in kwargs.keys():
+        kwargs['alpha'] = 0.3
+    axes.fill(np.hstack((x,x[::-1])),np.hstack((mu+err,mu[::-1]-err[::-1])),color=fillcol,**kwargs)
+
+    #this is the edge:
+    axes.plot(x,mu+err,color=edgecol,linewidth=0.2)
+    axes.plot(x,mu-err,color=edgecol,linewidth=0.2)
+
+def removeRightTicks(ax=None):
+    ax = ax or pb.gca()
+    for i, line in enumerate(ax.get_yticklines()):
+        if i%2 == 1:   # odd indices
+            line.set_visible(False)
+def removeUpperTicks(ax=None):
+    ax = ax or pb.gca()
+    for i, line in enumerate(ax.get_xticklines()):
+        if i%2 == 1:   # odd indices
+            line.set_visible(False)
+def fewerXticks(ax=None,divideby=2):
+    ax = ax or pb.gca()
+    ax.set_xticks(ax.get_xticks()[::divideby])
+
+def align_subplots(N,M,xlim=None, ylim=None):
+    """make all of the subplots have the same limits, turn off unnecessary ticks"""
+    #find sensible xlim,ylim
+    if xlim is None:
+        xlim = [np.inf,-np.inf]
+        for i in range(N*M):
+            pb.subplot(N,M,i+1)
+            xlim[0] = min(xlim[0],pb.xlim()[0])
+            xlim[1] = max(xlim[1],pb.xlim()[1])
+    if ylim is None:
+        ylim = [np.inf,-np.inf]
+        for i in range(N*M):
+            pb.subplot(N,M,i+1)
+            ylim[0] = min(ylim[0],pb.ylim()[0])
+            ylim[1] = max(ylim[1],pb.ylim()[1])
+
+    for i in range(N*M):
+        pb.subplot(N,M,i+1)
+        pb.xlim(xlim)
+        pb.ylim(ylim)
+        if (i)%M:
+            pb.yticks([])
+        else:
+            removeRightTicks()
+        if i<(M*(N-1)):
+            pb.xticks([])
+        else:
+            removeUpperTicks()
+
+
diff --git a/GPy/util/squashers.py b/GPy/util/squashers.py
new file mode 100644
index 00000000..b338c343
--- /dev/null
+++ b/GPy/util/squashers.py
@@ -0,0 +1,14 @@
+import numpy as np
+
+def sigmoid(x):
+    return 1./(1.+np.exp(-x))
+
+def softmax(x):
+    ex = np.exp(x-x.max(1)[:,None])
+    return ex/ex.sum(1)[:,np.newaxis]
+
+def single_softmax(x):
+    ex = np.exp(x)
+    return ex/ex.sum()
+
+
diff --git a/GPy/util/warping_functions.py b/GPy/util/warping_functions.py
new file mode 100644
index 00000000..35450b59
--- /dev/null
+++ b/GPy/util/warping_functions.py
@@ -0,0 +1,153 @@
+import numpy as np
+import scipy as sp
+import pylab as plt
+
+class WarpingFunction(object):
+    """
+    abstract function for warping
+    z = f(y)
+    """
+
+    def __init__(self):
+        raise NotImplementedError
+
+    def f(self,y,psi):
+        """function transformation
+        y is a list of values (GP training data) of shpape [N,1]
+        """
+        raise NotImplementedError
+
+    def fgrad_y(self,y,psi):
+        """gradient of f w.r.t to y"""
+        raise NotImplementedError
+
+    def fgrad_y_psi(self,y,psi):
+        """gradient of f w.r.t to y"""
+        raise NotImplementedError
+
+    def f_inv(self,z,psi):
+        """inverse function transformation"""
+        raise NotImplementedError
+
+    def get_param_names(self):
+        raise NotImplementedError
+
+    def plot(self, psi, xmin, xmax):
+        y = np.arange(xmin, xmax, 0.01)
+        f_y = self.f(y, psi)
+        plt.figure()
+        plt.plot(y, f_y)
+        plt.xlabel('y')
+        plt.ylabel('f(y)')
+        plt.title('warping function')
+
+class TanhWarpingFunction(WarpingFunction):
+
+    def __init__(self,n_terms=3):
+        """n_terms specifies the number of tanh terms to be used"""
+        self.n_terms = n_terms
+        self.num_parameters = 3 * self.n_terms
+
+    def f(self,y,psi):
+        """transform y with f using parameter vector psi
+        psi = [[a,b,c]]
+        f = \sum_{terms} a * tanh(b*(y+c))
+        """
+
+        #1. check that number of params is consistent
+        assert psi.shape[0] == self.n_terms, 'inconsistent parameter dimensions'
+        assert psi.shape[1] == 3, 'inconsistent parameter dimensions'
+
+        #2. exponentiate the a and b (positive!)
+        mpsi = psi.copy()
+
+        #3. transform data
+        z = y.copy()
+        for i in range(len(mpsi)):
+            a,b,c = mpsi[i]
+            z += a*np.tanh(b*(y+c))
+        return z
+
+
+    def f_inv(self, y, psi, iterations = 10):
+        """
+        calculate the numerical inverse of f
+
+        == input ==
+        iterations: number of N.R. iterations
+
+        """
+
+        y = y.copy()
+        z = np.ones_like(y)
+
+        for i in range(iterations):
+            z -= (self.f(z, psi) - y)/self.fgrad_y(z,psi)
+
+        return z
+
+
+    def fgrad_y(self, y, psi, return_precalc = False):
+        """
+        gradient of f w.r.t to y ([N x 1])
+        returns: Nx1 vector of derivatives, unless return_precalc is true,
+        then it also returns the precomputed stuff
+        """
+
+        mpsi = psi.copy()
+
+        # vectorized version
+
+        # S = (mpsi[:,1]*(y + mpsi[:,2])).T
+        S = (mpsi[:,1]*(y[:,:,None] + mpsi[:,2])).T
+        R = np.tanh(S)
+        D = 1-R**2
+
+        # GRAD = (1+(mpsi[:,0:1]*mpsi[:,1:2]*D).sum(axis=0))[:,np.newaxis]
+        GRAD = (1+(mpsi[:,0:1][:,:,None]*mpsi[:,1:2][:,:,None]*D).sum(axis=0)).T
+
+        if return_precalc:
+            # return GRAD,S.sum(axis=1),R.sum(axis=1),D.sum(axis=1)
+            return GRAD, S, R, D
+
+
+        return GRAD
+
+
+    def fgrad_y_psi(self, y, psi, return_covar_chain = False):
+        """
+        gradient of f w.r.t to y and psi
+
+        returns: NxIx3 tensor of partial derivatives
+
+        """
+
+        # 1. exponentiate the a and b (positive!)
+        mpsi = psi.copy()
+        w, s, r, d = self.fgrad_y(y, psi, return_precalc = True)
+
+        gradients = np.zeros((y.shape[0], y.shape[1], len(mpsi), 3))
+        for i in range(len(mpsi)):
+            a,b,c = mpsi[i]
+            gradients[:,:,i,0] = (b*(1.0/np.cosh(s[i]))**2).T
+            gradients[:,:,i,1] = a*(d[i] - 2.0*s[i]*r[i]*(1.0/np.cosh(s[i]))**2).T
+            gradients[:,:,i,2] = (-2.0*a*(b**2)*r[i]*((1.0/np.cosh(s[i]))**2)).T
+
+
+        if return_covar_chain:
+            covar_grad_chain = np.zeros((y.shape[0], y.shape[1], len(mpsi), 3))
+
+            for i in range(len(mpsi)):
+                a,b,c = mpsi[i]
+                covar_grad_chain[:, :, i, 0] = (r[i]).T
+                covar_grad_chain[:, :, i, 1] = (a*(y + c) * ((1.0/np.cosh(s[i]))**2).T)
+                covar_grad_chain[:, :, i, 2] = a*b*((1.0/np.cosh(s[i]))**2).T
+
+            return gradients, covar_grad_chain
+
+        return gradients
+
+    def get_param_names(self):
+        variables = ['a', 'b', 'c']
+        names = sum([['warp_tanh_%s_t%i' % (variables[n],q) for n in range(3)] for q in range(self.n_terms)],[])
+        return names