Merge pull request #326 from SheffieldML/kern

[kernel] fix #218 and #325

.coveragerc

@@ -2,7 +2,7 @@
 [run]
 branch = True
 source = GPy
-omit = ./GPy/testing/*.py, travis_tests.py, setup.py, ./GPy/__version__.py
+omit = ./GPy/testing/*.py, travis_tests.py, setup.py, ./GPy/__version__.py, ./GPy/plotting/*
 
 [report]
 # Regexes for lines to exclude from consideration

GPy/core/parameterization/__init__.py

@@ -3,7 +3,7 @@
 
 from .param import Param
 from .parameterized import Parameterized
-from paramz import transformations
+from . import transformations
 
 from paramz.core import lists_and_dicts, index_operations, observable_array, observable
 from paramz import ties_and_remappings, ObsAr

GPy/kern/src/add.py

@@ -19,8 +19,8 @@ class Add(CombinationKernel):
             if isinstance(kern, Add):
                 del subkerns[i]
                 for part in kern.parts[::-1]:
-                    kern.unlink_parameter(part)
+                    #kern.unlink_parameter(part)
-                    subkerns.insert(i, part)
+                    subkerns.insert(i, part.copy())
         super(Add, self).__init__(subkerns, name)
         self._exact_psicomp = self._check_exact_psicomp()
 
@@ -241,16 +241,20 @@ class Add(CombinationKernel):
             [np.add(target_grads[i],grads[i],target_grads[i]) for i in range(len(grads))]
         return target_grads
 
-    def add(self, other):
+    #def add(self, other):
-        if isinstance(other, Add):
+    #    parts = self.parts
-            other_params = other.parameters[:]
+    #    if 0:#isinstance(other, Add):
-            for p in other_params:
+    #        #other_params = other.parameters[:]
-                other.unlink_parameter(p)
+    #        for p in other.parts[:]:
-            self.link_parameters(*other_params)
+    #            other.unlink_parameter(p)
-        else:
+    #        parts.extend(other.parts)
-            self.link_parameter(other)
+    #        #self.link_parameters(*other_params)
-        self.input_dim, self._all_dims_active = self.get_input_dim_active_dims(self.parts)
+    #        
-        return self
+    #    else:
+    #        #self.link_parameter(other)
+    #        parts.append(other)
+    #    #self.input_dim, self._all_dims_active = self.get_input_dim_active_dims(parts)
+    #    return Add([p for p in parts], self.name)
 
     def input_sensitivity(self, summarize=True):
         if summarize:

GPy/kern/src/kern.py

@@ -48,11 +48,12 @@ class Kern(Parameterized):
 
         if active_dims is None:
             active_dims = np.arange(input_dim)
-
+        
-        self.active_dims = active_dims
+        self.active_dims = np.asarray(active_dims, np.int_)
-        self._all_dims_active = np.atleast_1d(active_dims).astype(int)
+        
-
+        self._all_dims_active = np.atleast_1d(self.active_dims).astype(int)
-        assert self._all_dims_active.size == self.input_dim, "input_dim={} does not match len(active_dim)={}, _all_dims_active={}".format(self.input_dim, self._all_dims_active.size, self._all_dims_active)
+        
+        assert self.active_dims.size == self.input_dim, "input_dim={} does not match len(active_dim)={}".format(self.input_dim, self._all_dims_active.size)
 
         self._sliced_X = 0
         self.useGPU = self._support_GPU and useGPU
@@ -322,10 +323,20 @@ class CombinationKernel(Kern):
         :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)
+        extra_dims = np.asarray(extra_dims, dtype=int)
-        input_dim, active_dims = self.get_input_dim_active_dims(kernels, extra_dims)
+        
+        active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int))
+        
+        input_dim = active_dims.size
+        if extra_dims is not None:
+            input_dim += extra_dims.size
+
         # initialize the kernel with the full input_dim
         super(CombinationKernel, self).__init__(input_dim, active_dims, name)
+
+        effective_input_dim = reduce(max, (k._all_dims_active.max() for k in kernels)) + 1
+        self._all_dims_active = np.array(np.concatenate((np.arange(effective_input_dim), extra_dims if extra_dims is not None else [])), dtype=int)
+        
         self.extra_dims = extra_dims
         self.link_parameters(*kernels)
 
@@ -333,16 +344,8 @@ class CombinationKernel(Kern):
     def parts(self):
         return self.parameters
 
-    def get_input_dim_active_dims(self, kernels, extra_dims = None):
+    def _set_all_dims_ative(self):
-        self.active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int))
+        self._all_dims_active = np.atleast_1d(self.active_dims).astype(int)        
-        #_all_dims_active = np.array(np.concatenate((_all_dims_active, extra_dims if extra_dims is not None else [])), dtype=int)
-        input_dim = reduce(max, (k._all_dims_active.max() for k in kernels)) + 1
-
-        if extra_dims is not None:
-            input_dim += extra_dims.size
-
-        _all_dims_active = np.arange(input_dim)
-        return input_dim, _all_dims_active
 
     def input_sensitivity(self, summarize=True):
         """

GPy/plotting/gpy_plot/latent_plots.py

@@ -147,6 +147,7 @@ def _plot_magnification(self, canvas, which_indices, Xgrid,
     def plot_function(x):
         Xtest_full = np.zeros((x.shape[0], Xgrid.shape[1]))
         Xtest_full[:, which_indices] = x
+
         mf = self.predict_magnification(Xtest_full, kern=kern, mean=mean, covariance=covariance)
         return mf.reshape(resolution, resolution).T
     imshow_kwargs = update_not_existing_kwargs(imshow_kwargs, pl().defaults.magnification)
@@ -215,7 +216,12 @@ def _plot_latent(self, canvas, which_indices, Xgrid,
     def plot_function(x):
         Xtest_full = np.zeros((x.shape[0], Xgrid.shape[1]))
         Xtest_full[:, which_indices] = x
-        mf = np.log(self.predict(Xtest_full, kern=kern)[1])
+        mf = self.predict(Xtest_full, kern=kern)[1]
+        if mf.shape[1]==self.output_dim:
+            mf = mf.sum(-1)
+        else:
+            mf *= self.output_dim
+        mf = np.log(mf)
         return mf.reshape(resolution, resolution).T
 
     imshow_kwargs = update_not_existing_kwargs(imshow_kwargs, pl().defaults.latent)

GPy/testing/kernel_tests.py

@@ -344,11 +344,14 @@ class KernelTestsMiscellaneous(unittest.TestCase):
         N, D = 100, 10
         self.X = np.linspace(-np.pi, +np.pi, N)[:,None] * np.random.uniform(-10,10,D)
         self.rbf = GPy.kern.RBF(2, active_dims=np.arange(0,4,2))
+        self.rbf.randomize()
         self.linear = GPy.kern.Linear(2, active_dims=(3,9))
+        self.linear.randomize()
         self.matern = GPy.kern.Matern32(3, active_dims=np.array([1,7,9]))
+        self.matern.randomize()
         self.sumkern = self.rbf + self.linear
         self.sumkern += self.matern
-        self.sumkern.randomize()
+        #self.sumkern.randomize()
 
     def test_which_parts(self):
         self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=[self.linear, self.matern]), self.linear.K(self.X)+self.matern.K(self.X)))
@@ -358,6 +361,21 @@ class KernelTestsMiscellaneous(unittest.TestCase):
     def test_active_dims(self):
         np.testing.assert_array_equal(self.sumkern.active_dims, [0,1,2,3,7,9])
         np.testing.assert_array_equal(self.sumkern._all_dims_active, range(10))
+        tmp = self.linear+self.rbf
+        np.testing.assert_array_equal(tmp.active_dims, [0,2,3,9])
+        np.testing.assert_array_equal(tmp._all_dims_active, range(10))
+        tmp = self.matern+self.rbf
+        np.testing.assert_array_equal(tmp.active_dims, [0,1,2,7,9])
+        np.testing.assert_array_equal(tmp._all_dims_active, range(10))
+        tmp = self.matern+self.rbf*self.linear
+        np.testing.assert_array_equal(tmp.active_dims, [0,1,2,3,7,9])
+        np.testing.assert_array_equal(tmp._all_dims_active, range(10))
+        tmp = self.matern+self.rbf+self.linear
+        np.testing.assert_array_equal(tmp.active_dims, [0,1,2,3,7,9])
+        np.testing.assert_array_equal(tmp._all_dims_active, range(10))
+        tmp = self.matern*self.rbf*self.linear
+        np.testing.assert_array_equal(tmp.active_dims, [0,1,2,3,7,9])
+        np.testing.assert_array_equal(tmp._all_dims_active, range(10))
 
 class KernelTestsNonContinuous(unittest.TestCase):
     def setUp(self):

GPy/testing/plotting_tests.py

@@ -100,7 +100,7 @@ def _image_comparison(baseline_images, extensions=['pdf','svg','png'], tol=11):
             fig.axes[0].set_axis_off()
             fig.set_frameon(False)
             fig.canvas.draw()
-            fig.savefig(os.path.join(result_dir, "{}.{}".format(base, ext)), transparent=True, edgecolor='none', facecolor='none')
+            fig.savefig(os.path.join(result_dir, "{}.{}".format(base, ext)), transparent=True, edgecolor='none', facecolor='none', bbox='tight')
     for num, base in zip(plt.get_fignums(), baseline_images):
         for ext in extensions:
             #plt.close(num)
@@ -116,7 +116,7 @@ def _image_comparison(baseline_images, extensions=['pdf','svg','png'], tol=11):
 def test_figure():
     np.random.seed(1239847)
     from GPy.plotting import plotting_library as pl
-    import matplotlib
+    #import matplotlib
     matplotlib.rcParams.update(matplotlib.rcParamsDefault)
     matplotlib.rcParams[u'figure.figsize'] = (4,3)
     matplotlib.rcParams[u'text.usetex'] = False
@@ -160,7 +160,7 @@ def test_figure():
 
 def test_kernel():
     np.random.seed(1239847)
-    import matplotlib
+    #import matplotlib
     matplotlib.rcParams.update(matplotlib.rcParamsDefault)
     matplotlib.rcParams[u'figure.figsize'] = (4,3)
     matplotlib.rcParams[u'text.usetex'] = False

GPy/testing/util_tests.py

@@ -0,0 +1,49 @@
+#===============================================================================
+# Copyright (c) 2016, Max Zwiessele
+# All rights reserved.
+# 
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+# 
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+# 
+# * Neither the name of GPy.testing.util_tests nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+# 
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#===============================================================================
+
+import unittest, numpy as np
+
+class TestDebug(unittest.TestCase):
+    def test_checkFinite(self):
+        from GPy.util.debug import checkFinite
+        array = np.random.normal(0, 1, 100).reshape(25,4)
+        self.assertTrue(checkFinite(array, name='test'))
+        
+        array[np.random.binomial(1, .3, array.shape).astype(bool)] = np.nan
+        self.assertFalse(checkFinite(array))
+
+    def test_checkFullRank(self):
+        from GPy.util.debug import checkFullRank
+        from GPy.util.linalg import tdot
+        array = np.random.normal(0, 1, 100).reshape(25,4)
+        self.assertFalse(checkFullRank(tdot(array), name='test'))
+        
+        array = np.random.normal(0, 1, (25,25))
+        self.assertTrue(checkFullRank(tdot(array)))

GPy/util/__init__.py

@@ -15,4 +15,4 @@ from . import diag
 from . import initialization
 from . import multioutput
 from . import parallel
-
+from . import functions

GPy/util/debug.py

@@ -22,7 +22,7 @@ def checkFullRank(m, tol=1e-10, name=None, force_check=False):
         name = 'Matrix with ID['+str(id(m))+']'
     assert len(m.shape)==2 and m.shape[0]==m.shape[1], 'The input of checkFullRank has to be a square matrix!'
 
-    if not force_check and m.shape[0]>=10000:
+    if not force_check and m.shape[0]>=10000: # pragma: no cover
         print('The size of '+name+'is too big to check (>=10000)!')
         return True
 

GPy/util/functions.py

@@ -1,27 +1,31 @@
 # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 import numpy as np
-from scipy.special import erf, erfc, erfcx
+from scipy import special
+from scipy.special import erfcx
 import sys
 epsilon = sys.float_info.epsilon
 lim_val = -np.log(epsilon) 
 
-def logisticln(x):
+def logisticln(x): # pragma: no cover
     return np.where(x<lim_val, np.where(x>-lim_val, -np.log(1+np.exp(-x)), -x), -np.log(1+epsilon))
 
-def logistic(x):
+def logistic(x): # pragma: no cover
-    return np.where(x<lim_val, np.where(x>-lim_val, 1/(1+np.exp(-x)), epsilon/(epsilon+1)), 1/(1+epsilon))
+    return special.expit(x)
+    #return np.where(x<lim_val, np.where(x>-lim_val, 1/(1+np.exp(-x)), epsilon/(epsilon+1)), 1/(1+epsilon))
 
-def normcdf(x):
+def normcdf(x): # pragma: no cover
-    g=0.5*erfc(-x/np.sqrt(2))
+    return special.ndtr(x)
-    return np.where(g==0, epsilon, np.where(g==1, 1-epsilon, g)) 
+    #g=0.5*erfc(-x/np.sqrt(2))
+    #return np.where(g==0, epsilon, np.where(g==1, 1-epsilon, g)) 
 
-def normcdfln(x):
+def normcdfln(x): # pragma: no cover
-    return np.where(x < 0, -.5*x*x + np.log(.5) + np.log(erfcx(-x/np.sqrt(2))), np.log(normcdf(x)))
+    return special.log_ndtr(x)
+    #return np.where(x < 0, -.5*x*x + np.log(.5) + np.log(erfcx(-x/np.sqrt(2))), np.log(normcdf(x)))
 
-def clip_exp(x):
+def clip_exp(x): # pragma: no cover
     return np.where(x<lim_val, np.where(x>-lim_val, np.exp(x), epsilon), 1/epsilon)
 
-def differfln(x0, x1):
+def differfln(x0, x1): # pragma: no cover
     # this is a, hopefully!, a numerically more stable variant of log(erf(x0)-erf(x1)) = log(erfc(x1)-erfc(x0)).
     return np.where(x0>x1, -x1*x1 + np.log(erfcx(x1)-np.exp(-x0**2+x1**2)*erfcx(x0)), -x0*x0 + np.log(np.exp(-x1**2+x0**2)*erfcx(x1) - erfcx(x0)))