# Source code for GPy.kern.src.coregionalize

# Copyright (c) 2012, James Hensman and Ricardo Andrade

from .kern import Kern
import numpy as np
from ...core.parameterization import Param
from paramz.transformations import Logexp
from ...util.config import config # for assesing whether to use cython

try:
from . import coregionalize_cython
use_coregionalize_cython = config.getboolean('cython', 'working')
except ImportError:
print('warning in coregionalize: failed to import cython module: falling back to numpy')
use_coregionalize_cython = False

[docs]class Coregionalize(Kern):
"""
Covariance function for intrinsic/linear coregionalization models

This covariance has the form:

.. math::
\mathbf{B} = \mathbf{W}\mathbf{W}^\intercal + \mathrm{diag}(kappa)

An intrinsic/linear coregionalization covariance function of the form:

.. math::
k_2(x, y)=\mathbf{B} k(x, y)

it is obtained as the tensor product between a covariance function
k(x, y) and B.

:param output_dim: number of outputs to coregionalize
:type output_dim: int
:param rank: number of columns of the W matrix (this parameter is ignored if parameter W is not None)
:type rank: 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  (output_dim, )

.. note: see coregionalization examples in GPy.examples.regression for some usage.
"""
def __init__(self, input_dim, output_dim, rank=1, W=None, kappa=None, active_dims=None, name='coregion'):
super(Coregionalize, self).__init__(input_dim, active_dims, name=name)
self.output_dim = output_dim
self.rank = rank
if self.rank>output_dim:
print("Warning: Unusual choice of rank, it should normally be less than the output_dim.")
if W is None:
W = 0.5*np.random.randn(self.output_dim, self.rank)/np.sqrt(self.rank)
else:
assert W.shape==(self.output_dim, self.rank)
self.W = Param('W', W)
if kappa is None:
kappa = 0.5*np.ones(self.output_dim)
else:
assert kappa.shape==(self.output_dim, )
self.kappa = Param('kappa', kappa, Logexp())

[docs]    def parameters_changed(self):
self.B = np.dot(self.W, self.W.T) + np.diag(self.kappa)

[docs]    def K(self, X, X2=None):
if use_coregionalize_cython:
return self._K_cython(X, X2)
else:
return self._K_numpy(X, X2)

def _K_numpy(self, X, X2=None):
index = np.asarray(X, dtype=np.int)
if X2 is None:
return self.B[index,index.T]
else:
index2 = np.asarray(X2, dtype=np.int)
return self.B[index,index2.T]

def _K_cython(self, X, X2=None):
if X2 is None:
return coregionalize_cython.K_symmetric(self.B, np.asarray(X, dtype=np.int64)[:,0])
return coregionalize_cython.K_asymmetric(self.B, np.asarray(X, dtype=np.int64)[:,0], np.asarray(X2, dtype=np.int64)[:,0])

[docs]    def Kdiag(self, X):
return np.diag(self.B)[np.asarray(X, dtype=np.int).flatten()]

[docs]    def update_gradients_full(self, dL_dK, X, X2=None):
index = np.asarray(X, dtype=np.int)
if X2 is None:
index2 = index
else:
index2 = np.asarray(X2, dtype=np.int)

#attempt to use cython for a nasty double indexing loop: fall back to numpy
if use_coregionalize_cython:
else:

dkappa = np.diag(dL_dK_small).copy()
dL_dK_small += dL_dK_small.T
dW = (self.W[:, None, :]*dL_dK_small[:, :, None]).sum(0)

index, index2 = index[:,0], index2[:,0]
dL_dK_small = np.zeros_like(self.B)
for i in range(self.output_dim):
tmp1 = dL_dK[index==i]
for j in range(self.output_dim):
dL_dK_small[j,i] = tmp1[:,index2==j].sum()
return dL_dK_small

index, index2 = np.int64(index[:,0]), np.int64(index2[:,0])