diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 8a334278..f9d7a99a 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -140,6 +140,33 @@ def white(input_dim,variance=1.): part = parts.white.White(input_dim,variance) return kern(input_dim, [part]) +def eq_ode1(output_dim, W=None, rank=1, kappa=None, length_scale=1., decay=None, delay=None): + """Covariance function for first order differential equation driven by an exponentiated quadratic covariance. + + This outputs of this kernel have the form + .. math:: + \frac{\text{d}y_j}{\text{d}t} = \sum_{i=1}^R w_{j,i} f_i(t-\delta_j) +\sqrt{\kappa_j}g_j(t) - d_jy_j(t) + + where :math:`R` is the rank of the system, :math:`w_{j,i}` is the sensitivity of the :math:`j`th output to the :math:`i`th latent function, :math:`d_j` is the decay rate of the :math:`j`th output and :math:`f_i(t)` and :math:`g_i(t)` are independent latent Gaussian processes goverened by an exponentiated quadratic covariance. + + :param output_dim: number of outputs driven by latent function. + :type output_dim: int + :param W: sensitivities of each output to the latent driving function. + :type W: ndarray (output_dim x rank). + :param rank: If rank is greater than 1 then there are assumed to be a total of rank latent forces independently driving the system, each with identical covariance. + :type rank: int + :param decay: decay rates for the first order system. + :type decay: array of length output_dim. + :param delay: delay between latent force and output response. + :type delay: array of length output_dim. + :param kappa: diagonal term that allows each latent output to have an independent component to the response. + :type kappa: array of length output_dim. + + .. Note: see first order differential equation examples in GPy.examples.regression for some usage. + """ + part = parts.eq_ode1.Eq_ode1(output_dim, W, rank, kappa, length_scale, decay, delay) + return kern(2, [part]) + def exponential(input_dim,variance=1., lengthscale=None, ARD=False): """ @@ -346,7 +373,7 @@ def symmetric(k): k_.parts = [symmetric.Symmetric(p) for p in k.parts] return k_ -def coregionalize(num_outputs,W_columns=1, W=None, kappa=None): +def coregionalize(output_dim,W_columns=1, W=None, kappa=None): """ Coregionlization matrix B, of the form: .. math:: @@ -358,18 +385,18 @@ def coregionalize(num_outputs,W_columns=1, W=None, kappa=None): it is obtainded as the tensor product between a kernel k(x,y) and B. - :param num_outputs: the number of outputs to corregionalize - :type num_outputs: int + :param output_dim: the number of outputs to corregionalize + :type output_dim: int :param W_columns: number of columns of the W matrix (this parameter is ignored if parameter W is not None) :type W_colunns: int :param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalisation matrix B :type W: numpy array of dimensionality (num_outpus, W_columns) :param kappa: a vector which allows the outputs to behave independently - :type kappa: numpy array of dimensionality (num_outputs,) + :type kappa: numpy array of dimensionality (output_dim,) :rtype: kernel object """ - p = parts.coregionalize.Coregionalize(num_outputs,W_columns,W,kappa) + p = parts.coregionalize.Coregionalize(output_dim,W_columns,W,kappa) return kern(1,[p]) @@ -429,12 +456,12 @@ def hierarchical(k): _parts = [parts.hierarchical.Hierarchical(k.parts)] return kern(k.input_dim+len(k.parts),_parts) -def build_lcm(input_dim, num_outputs, kernel_list = [], W_columns=1,W=None,kappa=None): +def build_lcm(input_dim, output_dim, kernel_list = [], W_columns=1,W=None,kappa=None): """ Builds a kernel of a linear coregionalization model :input_dim: Input dimensionality - :num_outputs: Number of outputs + :output_dim: Number of outputs :kernel_list: List of coregionalized kernels, each element in the list will be multiplied by a different corregionalization matrix :type kernel_list: list of GPy kernels :param W_columns: number tuples of the corregionalization parameters 'coregion_W' @@ -448,11 +475,11 @@ def build_lcm(input_dim, num_outputs, kernel_list = [], W_columns=1,W=None,kappa k.input_dim = input_dim warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.") - k_coreg = coregionalize(num_outputs,W_columns,W,kappa) + k_coreg = coregionalize(output_dim,W_columns,W,kappa) kernel = kernel_list[0]**k_coreg.copy() for k in kernel_list[1:]: - k_coreg = coregionalize(num_outputs,W_columns,W,kappa) + k_coreg = coregionalize(output_dim,W_columns,W,kappa) kernel += k**k_coreg.copy() return kernel diff --git a/GPy/kern/parts/__init__.py b/GPy/kern/parts/__init__.py index 643483f5..fff24d04 100644 --- a/GPy/kern/parts/__init__.py +++ b/GPy/kern/parts/__init__.py @@ -2,10 +2,11 @@ import bias import Brownian import coregionalize import exponential +import eq_ode1 import finite_dimensional import fixed import gibbs -#import hetero #hetero.py is not commited: omitting for now. JH. +import hetero #hetero.py is not commited: omitting for now. JH. import hierarchical import independent_outputs import linear diff --git a/GPy/kern/parts/coregionalize.py b/GPy/kern/parts/coregionalize.py index 4dc8d909..8e960038 100644 --- a/GPy/kern/parts/coregionalize.py +++ b/GPy/kern/parts/coregionalize.py @@ -13,7 +13,7 @@ class Coregionalize(Kernpart): This covariance has the form .. math:: - \mathbf{B} = \mathbf{W}\mathbf{W}^\top + kappa \mathbf{I} + \mathbf{B} = \mathbf{W}\mathbf{W}^\top + \text{diag}(kappa) An intrinsic/linear coregionalization covariance function of the form .. math:: @@ -22,33 +22,33 @@ class Coregionalize(Kernpart): it is obtained as the tensor product between a covariance function k(x,y) and B. - :param num_outputs: number of outputs to coregionalize - :type num_outputs: int + :param output_dim: number of outputs to coregionalize + :type output_dim: int :param W_columns: number of columns of the W matrix (this parameter is ignored if parameter W is not None) :type W_colunns: int :param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalization matrix B :type W: numpy array of dimensionality (num_outpus, W_columns) :param kappa: a vector which allows the outputs to behave independently - :type kappa: numpy array of dimensionality (num_outputs,) + :type kappa: numpy array of dimensionality (output_dim,) .. Note: see coregionalization examples in GPy.examples.regression for some usage. """ - def __init__(self,num_outputs,W_columns=1, W=None, kappa=None): + def __init__(self, output_dim, rank=1, W=None, kappa=None): self.input_dim = 1 self.name = 'coregion' - self.num_outputs = num_outputs - self.W_columns = W_columns + self.output_dim = output_dim + self.rank = rank if W is None: - self.W = 0.5*np.random.randn(self.num_outputs,self.W_columns)/np.sqrt(self.W_columns) + self.W = 0.5*np.random.randn(self.output_dim,self.rank)/np.sqrt(self.rank) else: - assert W.shape==(self.num_outputs,self.W_columns) + assert W.shape==(self.output_dim,self.rank) self.W = W if kappa is None: - kappa = 0.5*np.ones(self.num_outputs) + kappa = 0.5*np.ones(self.output_dim) else: - assert kappa.shape==(self.num_outputs,) + assert kappa.shape==(self.output_dim,) self.kappa = kappa - self.num_params = self.num_outputs*(self.W_columns + 1) + self.num_params = self.output_dim*(self.rank + 1) self._set_params(np.hstack([self.W.flatten(),self.kappa])) def _get_params(self): @@ -56,12 +56,12 @@ class Coregionalize(Kernpart): def _set_params(self,x): assert x.size == self.num_params - self.kappa = x[-self.num_outputs:] - self.W = x[:-self.num_outputs].reshape(self.num_outputs,self.W_columns) + self.kappa = x[-self.output_dim:] + self.W = x[:-self.output_dim].reshape(self.output_dim,self.rank) self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa) def _get_param_names(self): - return sum([['W%i_%i'%(i,j) for j in range(self.W_columns)] for i in range(self.num_outputs)],[]) + ['kappa_%i'%i for i in range(self.num_outputs)] + return sum([['W%i_%i'%(i,j) for j in range(self.rank)] for i in range(self.output_dim)],[]) + ['kappa_%i'%i for i in range(self.output_dim)] def K(self,index,index2,target): index = np.asarray(index,dtype=np.int) @@ -79,26 +79,26 @@ class Coregionalize(Kernpart): if index2 is None: code=""" for(int i=0;i