merge mu's changes into devel

GPy/kern/__init__.py

@@ -3,7 +3,7 @@ from _src.rbf import RBF
 from _src.linear import Linear, LinearFull
 from _src.static import Bias, White
 from _src.brownian import Brownian
-from _src.stationary import Exponential, Matern32, Matern52, ExpQuad, RatQuad, Cosine
+from _src.stationary import Exponential, OU, Matern32, Matern52, ExpQuad, RatQuad, Cosine
 from _src.mlp import MLP
 from _src.periodic import PeriodicExponential, PeriodicMatern32, PeriodicMatern52
 from _src.independent_outputs import IndependentOutputs, Hierarchical
@@ -13,6 +13,8 @@ from _src.ODE_UY import ODE_UY
 from _src.ODE_UYC import ODE_UYC
 from _src.ODE_st import ODE_st
 from _src.ODE_t import ODE_t
+from _src.poly import Poly
+
 # TODO: put this in an init file somewhere
 #I'm commenting this out because the files were not added. JH. Remember to add the files before commiting
 try:

GPy/kern/_src/add.py

@@ -170,7 +170,4 @@ class Add(CombinationKernel):
         return self
 
     def input_sensitivity(self):
-        in_sen = np.zeros(self.input_dim)
-        for i, p in enumerate(self.parts):
-            in_sen[p.active_dims] += p.input_sensitivity()
-        return in_sen
+        return reduce(np.add, [k.input_sensitivity() for k in self.parts])

GPy/kern/_src/independent_outputs.py

@@ -32,7 +32,7 @@ def index_to_slices(index):
     [ret[ind_i].append(slice(*indexes_i)) for ind_i,indexes_i in zip(ind[switchpoints[:-1]],zip(switchpoints,switchpoints[1:]))]
     return ret
 
-class IndependentOutputs(Kern):
+class IndependentOutputs(CombinationKernel):
     """
     A kernel which can represent several independent functions.  this kernel
     'switches off' parts of the matrix where the output indexes are different.
@@ -180,6 +180,9 @@ class Hierarchical(CombinationKernel):
     def Kdiag(self,X):
         return np.diag(self.K(X))
 
+    def gradients_X(self, dL_dK, X, X2=None):
+        raise NotImplementedError
+
     def update_gradients_full(self,dL_dK,X,X2=None):
         slices = [index_to_slices(X[:,i]) for i in self.extra_dims]
         if X2 is None:

GPy/kern/_src/kern.py

@@ -34,36 +34,24 @@ class Kern(Parameterized):
 
             is the active_dimensions of inputs X we will work on.
             All kernels will get sliced Xes as inputs, if active_dims is not None
+            Only positive integers are allowed in active_dims!
             if active_dims is None, slicing is switched off and all X will be passed through as given.
 
         :param int input_dim: the number of input dimensions to the function
-        :param array-like|slice|None active_dims: list of indices on which dimensions this kernel works on, or none if no slicing
+        :param array-like|None active_dims: list of indices on which dimensions this kernel works on, or none if no slicing
 
         Do not instantiate.
         """
         super(Kern, self).__init__(name=name, *a, **kw)
-        try:
-            self.input_dim = int(input_dim)
-            self.active_dims = active_dims# if active_dims is not None else slice(0, input_dim, 1)
-        except TypeError:
-            # input_dim is something else then an integer
-            self.input_dim = input_dim
-            if active_dims is not None:
-                print "WARNING: given input_dim={} is not an integer and active_dims={} is given, switching off slicing"
-            self.active_dims = None
+        self.input_dim = int(input_dim)
+
+        if active_dims is None:
+            active_dims = np.arange(input_dim)
+
+        self.active_dims = np.array(active_dims, dtype=int)
+
+        assert self.active_dims.size == self.input_dim, "input_dim={} does not match len(active_dim)={}, active_dims={}".format(self.input_dim, self.active_dims.size, self.active_dims)
 
-        if self.active_dims is not None and self.input_dim is not None:
-            assert isinstance(self.active_dims, (slice, list, tuple, np.ndarray)), 'active_dims needs to be an array-like or slice object over dimensions, {} given'.format(self.active_dims.__class__)
-            if isinstance(self.active_dims, slice):
-                self.active_dims = slice(self.active_dims.start or 0, self.active_dims.stop or self.input_dim, self.active_dims.step or 1)
-                active_dim_size = int(np.round((self.active_dims.stop-self.active_dims.start)/self.active_dims.step))
-            elif isinstance(self.active_dims, np.ndarray):
-                #assert np.all(self.active_dims >= 0), 'active dimensions need to be positive. negative indexing is not allowed'
-                assert self.active_dims.ndim == 1, 'only flat indices allowed, given active_dims.shape={}, provide only indexes to the dimensions (columns) of the input'.format(self.active_dims.shape)
-                active_dim_size = self.active_dims.size
-            else:
-                active_dim_size = len(self.active_dims)
-            assert active_dim_size == self.input_dim, "input_dim={} does not match len(active_dim)={}, active_dims={}".format(self.input_dim, active_dim_size, self.active_dims)
         self._sliced_X = 0
         self.useGPU = self._support_GPU and useGPU
 
@@ -176,8 +164,8 @@ class Kern(Parameterized):
         """
         Shortcut for tensor `prod`.
         """
-        assert self.active_dims == range(self.input_dim), "Can only use kernels, which have their input_dims defined from 0"
-        assert other.active_dims == range(other.input_dim), "Can only use kernels, which have their input_dims defined from 0"
+        assert np.all(self.active_dims == range(self.input_dim)), "Can only use kernels, which have their input_dims defined from 0"
+        assert np.all(other.active_dims == range(other.input_dim)), "Can only use kernels, which have their input_dims defined from 0"
         other.active_dims += self.input_dim
         return self.prod(other)
 
@@ -202,12 +190,12 @@ class Kern(Parameterized):
         return Prod([self, other], name)
 
     def _check_input_dim(self, X):
-        assert X.shape[1] == self.input_dim, "You did not specify active_dims and X has wrong shape: X_dim={}, whereas input_dim={}".format(X.shape[1], self.input_dim)
-    
-    def _check_active_dims(self, X):
-        assert X.shape[1] >= len(np.r_[self.active_dims]), "At least {} dimensional X needed, X.shape={!s}".format(len(np.r_[self.active_dims]), X.shape)
+        assert X.shape[1] == self.input_dim, "{} did not specify active_dims and X has wrong shape: X_dim={}, whereas input_dim={}".format(self.name, X.shape[1], self.input_dim)
+
+    def _check_active_dims(self, X):
+        assert X.shape[1] >= len(self.active_dims), "At least {} dimensional X needed, X.shape={!s}".format(len(self.active_dims), X.shape)
+
 
-            
 class CombinationKernel(Kern):
     """
     Abstract super class for combination kernels.
@@ -222,9 +210,10 @@ class CombinationKernel(Kern):
 
         :param list kernels: List of kernels to combine (can be only one element)
         :param str name: name of the combination kernel
-        :param array-like|slice extra_dims: if needed extra dimensions for the combination kernel to work on
+        :param array-like extra_dims: if needed extra dimensions for the combination kernel to work on
         """
         assert all([isinstance(k, Kern) for k in kernels])
+        extra_dims = np.array(extra_dims, dtype=int)
         input_dim, active_dims = self.get_input_dim_active_dims(kernels, extra_dims)
         # initialize the kernel with the full input_dim
         super(CombinationKernel, self).__init__(input_dim, active_dims, name)
@@ -238,16 +227,18 @@ class CombinationKernel(Kern):
     def get_input_dim_active_dims(self, kernels, extra_dims = None):
         #active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int))
         #active_dims = np.array(np.concatenate((active_dims, extra_dims if extra_dims is not None else [])), dtype=int)
-        input_dim = np.array([k.input_dim for k in kernels])
-        if np.all(input_dim[0]==input_dim):
-            input_dim = input_dim[0]
-        active_dims = None
+        input_dim = reduce(max, (k.active_dims.max() for k in kernels)) + 1
+
+        if extra_dims is not None:
+            input_dim += extra_dims.size
+
+        active_dims = np.arange(input_dim)
         return input_dim, active_dims
 
     def input_sensitivity(self):
         raise NotImplementedError("Choose the kernel you want to get the sensitivity for. You need to override the default behaviour for getting the input sensitivity to be able to get the input sensitivity. For sum kernel it is the sum of all sensitivities, TODO: product kernel? Other kernels?, also TODO: shall we return all the sensitivities here in the combination kernel? So we can combine them however we want? This could lead to just plot all the sensitivities here...")
 
-    def _check_input_dim(self, X):
+    def _check_active_dims(self, X):
         return
 
     def _check_input_dim(self, X):

GPy/kern/_src/linear.py

@@ -12,6 +12,7 @@ from ...core.parameterization.transformations import Logexp
 from ...util.caching import Cache_this
 from ...core.parameterization import variational
 from psi_comp import linear_psi_comp
+from ...util.config import *
 
 class Linear(Kern):
     """
@@ -224,12 +225,23 @@ class Linear(Kern):
         AZZA = ZA.T[:, None, :, None] * ZA[None, :, None, :]
         AZZA = AZZA + AZZA.swapaxes(1, 2)
         AZZA_2 = AZZA/2.
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for private(m,mm,q,qq,factor,tmp)'
+            header_string = '#include <omp.h>'
+            weave_options = {'headers'           : ['<omp.h>'],
+                                'extra_compile_args': ['-fopenmp -O3'],
+                                'extra_link_args'   : ['-lgomp'],
+                                'libraries': ['gomp']}
+        else:
+            pragma_string = ''
+            header_string = ''
+            weave_options = {'extra_compile_args': ['-O3']}
 
         #Using weave, we can exploit the symmetry of this problem:
         code = """
         int n, m, mm,q,qq;
         double factor,tmp;
-        #pragma omp parallel for private(m,mm,q,qq,factor,tmp)
+        %s
         for(n=0;n<N;n++){
           for(m=0;m<num_inducing;m++){
             for(mm=0;mm<=m;mm++){
@@ -253,26 +265,36 @@ class Linear(Kern):
             }
           }
         }
-        """
+        """ % pragma_string
         support_code = """
-        #include <omp.h>
+        %s
         #include <math.h>
-        """
-        weave_options = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
-                         'extra_link_args'   : ['-lgomp']}
+        """ % header_string
         mu = vp.mean
         N,num_inducing,input_dim,mu = mu.shape[0],Z.shape[0],mu.shape[1],param_to_array(mu)
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+        weave.inline(code, support_code=support_code,
                      arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
 
 
     def _weave_dpsi2_dZ(self, dL_dpsi2, Z, vp, target):
         AZA = self.variances*self._ZAinner(vp, Z)
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for private(n,mm,q)'
+            header_string = '#include <omp.h>'
+            weave_options =  {'headers'           : ['<omp.h>'],
+                              'extra_compile_args': ['-fopenmp -O3'],
+                              'extra_link_args'   : ['-lgomp'],
+                              'libraries': ['gomp']}
+        else:
+            pragma_string = ''
+            header_string = ''
+            weave_options = {'extra_compile_args': ['-O3']}
+
         code="""
         int n,m,mm,q;
-        #pragma omp parallel for private(n,mm,q)
+        %s
         for(m=0;m<num_inducing;m++){
           for(q=0;q<input_dim;q++){
             for(mm=0;mm<num_inducing;mm++){
@@ -282,18 +304,15 @@ class Linear(Kern):
             }
           }
         }
-        """
+        """ % pragma_string
         support_code = """
-        #include <omp.h>
+        %s
         #include <math.h>
-        """
-        weave_options = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
-                         'extra_link_args'   : ['-lgomp']}
+        """ % header_string
 
         N,num_inducing,input_dim = vp.mean.shape[0],Z.shape[0],vp.mean.shape[1]
         mu = param_to_array(vp.mean)
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+        weave.inline(code, support_code=support_code,
                      arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
 

GPy/kern/_src/poly.py

@@ -0,0 +1,42 @@
+# Copyright (c) 2014, James Hensman
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import numpy as np
+from kern import Kern
+from ...util.misc import param_to_array
+from ...core.parameterization import Param
+from ...core.parameterization.transformations import Logexp
+class Poly(Kern):
+    """
+    Polynomial kernel
+    """
+
+    def __init__(self, input_dim, variance=1., order=3., active_dims=None, name='poly'):
+        super(Poly, self).__init__(input_dim, active_dims, name)
+        self.variance = Param('variance', variance, Logexp())
+        self.add_parameter(self.variance)
+        self.order=order
+
+    def K(self, X, X2=None):
+        return (self._dot_product(X, X2) + 1.)**self.order * self.variance
+
+    def _dot_product(self, X, X2=None):
+        if X2 is None:
+            return np.dot(X, X.T)
+        else:
+            return np.dot(X, X2.T)
+
+    def Kdiag(self, X):
+        return self.variance*(np.square(X).sum(1) + 1.)**self.order
+
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        self.variance.gradient = np.sum(dL_dK * (self._dot_product(X, X2) + 1.)**self.order)
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        raise NotImplementedError
+
+    def gradients_X(self, dL_dK, X, X2=None):
+        raise NotImplementedError
+
+    def gradients_X_diag(self, dL_dKdiag, X):
+        raise NotImplementedError

GPy/kern/_src/rbf.py

@@ -10,6 +10,7 @@ from GPy.util.caching import Cache_this
 from ...core.parameterization import variational
 from psi_comp import ssrbf_psi_comp
 from psi_comp.ssrbf_psi_gpucomp import PSICOMP_SSRBF
+from ...util.config import *
 
 class RBF(Stationary):
     """
@@ -231,6 +232,16 @@ class RBF(Stationary):
 
     @Cache_this(limit=1)
     def _psi2computations(self, Z, vp):
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for private(tmp, exponent_tmp)'
+            header_string = '#include <omp.h>'
+            libraries = ['gomp']
+        else:
+            pragma_string = ''
+            header_string = ''
+            libraries = []
+
         mu, S = vp.mean, vp.variance
 
         N, Q = mu.shape
@@ -253,8 +264,7 @@ class RBF(Stationary):
         variance_sq = float(np.square(self.variance))
         code = """
         double tmp, exponent_tmp;
-
-        #pragma omp parallel for private(tmp, exponent_tmp)
+        %s 
         for (int n=0; n<N; n++)
         {
             for (int m=0; m<M; m++)
@@ -278,20 +288,20 @@ class RBF(Stationary):
                         tmp = -Zdist_sq(m,mm,q) - tmp - half_log_denom(n,q);
                         exponent_tmp += tmp;
                     }
-                    //compute psi2 by exponontiating
+                    //compute psi2 by exponentiating
                     psi2(n,m,mm) = variance_sq * exp(exponent_tmp);
                     psi2(n,mm,m) = psi2(n,m,mm);
                 }
             }
         }
-        """
+        """ % pragma_string
 
         support_code = """
-        #include <omp.h>
+        %s
         #include <math.h>
-        """
+        """ % header_string
         mu = param_to_array(mu)
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+        weave.inline(code, support_code=support_code, libraries=libraries,
                      arg_names=['N', 'M', 'Q', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'denom_l2', 'Zdist_sq', 'half_log_denom', 'psi2', 'variance_sq'],
                      type_converters=weave.converters.blitz, **self.weave_options)
 
@@ -303,12 +313,20 @@ class RBF(Stationary):
         #return 2.*np.einsum( 'ijk,ijk,ijkl,il->l', dL_dpsi2, psi2, Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2, 1./(2.*S + l2))*self.lengthscale
 
         result = np.zeros(self.input_dim)
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for reduction(+:tmp)'
+            header_string = '#include <omp.h>'
+            libraries = ['gomp']
+        else:
+            pragma_string = ''
+            header_string = ''
+            libraries = []
         code = """
         double tmp;
         for(int q=0; q<Q; q++)
         {
             tmp = 0.0;
-            #pragma omp parallel for reduction(+:tmp)
+            %s
             for(int n=0; n<N; n++)
             {
                 for(int m=0; m<M; m++)
@@ -326,16 +344,16 @@ class RBF(Stationary):
             result(q) = tmp;
         }
 
-        """
+        """ % pragma_string
         support_code = """
-        #include <omp.h>
+        %s
         #include <math.h>
-        """
+        """ % header_string
         N,Q = S.shape
         M = psi2.shape[-1]
 
         S = param_to_array(S)
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+        weave.inline(code, support_code=support_code, libraries=libraries,
                      arg_names=['psi2', 'dL_dpsi2', 'N', 'M', 'Q', 'mudist_sq', 'l2', 'Zdist_sq', 'S', 'result'],
                      type_converters=weave.converters.blitz, **self.weave_options)
 

GPy/kern/_src/stationary.py

@@ -192,6 +192,27 @@ class Exponential(Stationary):
     def dK_dr(self, r):
         return -0.5*self.K_of_r(r)
 
+
+class OU(Stationary):
+    """
+    OU kernel:
+
+    .. math::
+
+       k(r) = \\sigma^2 \exp(- r) \\ \\ \\ \\  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
+
+    """
+
+    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='OU'):
+        super(OU, self).__init__(input_dim, variance, lengthscale, ARD, active_dims, name)
+
+    def K_of_r(self, r):
+        return self.variance * np.exp(-r)
+
+    def dK_dr(self,r):
+        return -1.*self.variance*np.exp(-r)
+
+
 class Matern32(Stationary):
     """
     Matern 3/2 kernel: