Pull branch 'devel' of github.com:SheffieldML/GPy into devel

GPy/core/model.py

@@ -13,6 +13,7 @@ import priors
 from ..util.linalg import jitchol
 from ..inference import optimization
 from .. import likelihoods
+import re
 
 class model(parameterised):
     def __init__(self):
@@ -245,7 +246,7 @@ class model(parameterised):
         for s in positive_strings:
             for i in self.grep_param_names(s):
                 if not (i in currently_constrained):
-                    to_make_positive.append(param_names[i])
+                    to_make_positive.append(re.escape(param_names[i]))
                     if warn:
                         print "Warning! constraining %s postive" % name
         if len(to_make_positive):

GPy/kern/__init__.py

@@ -2,5 +2,9 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, prod, prod_orthogonal, symmetric, coregionalise, rational_quadratic, fixed, rbfcos
+from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, periodic_exponential, periodic_Matern32, periodic_Matern52, prod, prod_orthogonal, symmetric, coregionalise, rational_quadratic, fixed, rbfcos, independent_outputs
+try:
+    from constructors import rbf_sympy, sympykern # these depend on sympy
+except:
+    pass
 from kern import kern

GPy/kern/constructors.py

@@ -25,6 +25,7 @@ from symmetric import symmetric as symmetric_part
 from coregionalise import coregionalise as coregionalise_part
 from rational_quadratic import rational_quadratic as rational_quadraticpart
 from rbfcos import rbfcos as rbfcospart
+from independent_outputs import independent_outputs as independent_output_part
 #TODO these s=constructors are not as clean as we'd like. Tidy the code up
 #using meta-classes to make the objects construct properly wthout them.
 
@@ -165,34 +166,40 @@ def Brownian(D,variance=1.):
     part = Brownianpart(D,variance)
     return kern(D, [part])
 
-import sympy as sp
+try:
-from sympykern import spkern
+    import sympy as sp
-from sympy.parsing.sympy_parser import parse_expr
+    from sympykern import spkern
+    from sympy.parsing.sympy_parser import parse_expr
+    sympy_available = True
+except ImportError:
+    sympy_available = False
 
-def rbf_sympy(D,ARD=False,variance=1., lengthscale=1.):
+if sympy_available:
-    """
+    def rbf_sympy(D,ARD=False,variance=1., lengthscale=1.):
-    Radial Basis Function covariance.
+        """
-    """
+        Radial Basis Function covariance.
-    X = [sp.var('x%i'%i) for i in range(D)]
+        """
-    Z = [sp.var('z%i'%i) for i in range(D)]
+        X = [sp.var('x%i'%i) for i in range(D)]
-    rbf_variance = sp.var('rbf_variance',positive=True)
+        Z = [sp.var('z%i'%i) for i in range(D)]
-    if ARD:
+        rbf_variance = sp.var('rbf_variance',positive=True)
-        rbf_lengthscales = [sp.var('rbf_lengthscale_%i'%i,positive=True) for i in range(D)]
+        if ARD:
-        dist_string = ' + '.join(['(x%i-z%i)**2/rbf_lengthscale_%i**2'%(i,i,i) for i in range(D)])
+            rbf_lengthscales = [sp.var('rbf_lengthscale_%i'%i,positive=True) for i in range(D)]
-        dist = parse_expr(dist_string)
+            dist_string = ' + '.join(['(x%i-z%i)**2/rbf_lengthscale_%i**2'%(i,i,i) for i in range(D)])
-        f =  rbf_variance*sp.exp(-dist/2.)
+            dist = parse_expr(dist_string)
-    else:
+            f =  rbf_variance*sp.exp(-dist/2.)
-        rbf_lengthscale = sp.var('rbf_lengthscale',positive=True)
+        else:
-        dist_string = ' + '.join(['(x%i-z%i)**2'%(i,i) for i in range(D)])
+            rbf_lengthscale = sp.var('rbf_lengthscale',positive=True)
-        dist = parse_expr(dist_string)
+            dist_string = ' + '.join(['(x%i-z%i)**2'%(i,i) for i in range(D)])
-        f =  rbf_variance*sp.exp(-dist/(2*rbf_lengthscale**2))
+            dist = parse_expr(dist_string)
-    return kern(D,[spkern(D,f)])
+            f =  rbf_variance*sp.exp(-dist/(2*rbf_lengthscale**2))
+        return kern(D,[spkern(D,f)])
 
-def sympykern(D,k):
+    def sympykern(D,k):
-    """
+        """
-    A kernel from a symbolic sympy representation
+        A kernel from a symbolic sympy representation
-    """
+        """
-    return kern(D,[spkern(D,k)])
+        return kern(D,[spkern(D,k)])
+del sympy_available
 
 def periodic_exponential(D=1,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
     """
@@ -318,3 +325,14 @@ def rbfcos(D,variance=1.,frequencies=None,bandwidths=None,ARD=False):
     """
     part = rbfcospart(D,variance,frequencies,bandwidths,ARD)
     return kern(D,[part])
+
+def independent_outputs(k):
+    """
+    Construct a kernel with independent outputs from an existing kernel
+    """
+    for sl in k.input_slices:
+        assert (sl.start is None) and (sl.stop is None), "cannot adjust input slices! (TODO)"
+    parts = [independent_output_part(p) for p in k.parts]
+    return kern(k.D+1,parts)
+
+

GPy/kern/independent_outputs.py

@@ -0,0 +1,97 @@
+# Copyright (c) 2012, James Hesnsman
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+from kernpart import kernpart
+import numpy as np
+
+def index_to_slices(index):
+    """
+    take a numpy array of integers (index) and return a  nested list of slices such that the slices describe the start, stop points for each integer in the index. 
+
+    e.g.
+    >>> index = np.asarray([0,0,0,1,1,1,2,2,2])
+    returns
+    >>> [[slice(0,3,None)],[slice(3,6,None)],[slice(6,9,None)]]
+
+    or, a more complicated example
+    >>> index = np.asarray([0,0,1,1,0,2,2,2,1,1])
+    returns
+    >>> [[slice(0,2,None),slice(4,5,None)],[slice(2,4,None),slice(8,10,None)],[slice(5,8,None)]]
+    """
+
+    #contruct the return structure
+    ind = np.asarray(index,dtype=np.int64)
+    ret = [[] for i in range(ind.max()+1)]
+
+    #find the switchpoints
+    ind_ = np.hstack((ind,ind[0]+ind[-1]+1))
+    switchpoints = np.nonzero(ind_ - np.roll(ind_,+1))[0]
+
+    [ret[ind_i].append(slice(*indexes_i)) for ind_i,indexes_i in zip(ind[switchpoints[:-1]],zip(switchpoints,switchpoints[1:]))]
+    return ret
+
+class independent_outputs(kernpart):
+    """
+    A kernel part shich can reopresent several independent functions.
+    this kernel 'switches off' parts of the matrix where the output indexes are different.
+
+    The index of the functions is given by the last column in the input X
+    the rest of the columns of X are passed to the kernel for computation (in blocks).
+
+    """
+    def __init__(self,k):
+        self.D = k.D + 1
+        self.Nparam = k.Nparam
+        self.name = 'iops('+ k.name + ')'
+        self.k = k
+
+    def _get_params(self):
+        return self.k._get_params()
+
+    def _set_params(self,x):
+        self.k._set_params(x)
+        self.params = x
+
+    def _get_param_names(self):
+        return self.k._get_param_names()
+
+    def K(self,X,X2,target):
+        #Sort out the slices from the input data
+        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        if X2 is None:
+            X2,slices2 = X,slices
+        else:
+            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
+
+        [[[self.k.K(X[s],X2[s2],target[s,s2]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+
+    def Kdiag(self,X,target):
+        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        [[self.k.Kdiag(X[s],target[s]) for s in slices_i] for slices_i in slices]
+
+    def dK_dtheta(self,dL_dK,X,X2,target):
+        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        if X2 is None:
+            X2,slices2 = X,slices
+        else:
+            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
+        [[[self.k.dK_dtheta(dL_dK[s,s2],X[s],X2[s2],target) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+
+
+    def dK_dX(self,dL_dK,X,X2,target):
+        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        if X2 is None:
+            X2,slices2 = X,slices
+        else:
+            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
+        [[[self.k.dK_dX(dL_dK[s,s2],X[s],X2[s2],target[s,:-1]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+
+    def dKdiag_dX(self,dL_dKdiag,X,target):
+        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        [[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target[s,:-1]) for s in slices_i] for slices_i in slices]
+
+
+    def dKdiag_dtheta(self,dL_dKdiag,X,target):
+        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        [[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target) for s in slices_i] for slices_i in slices]