add debug util module and try to debug sparsegp_mpi

GPy/core/sparse_gp_mpi.py

@@ -3,12 +3,7 @@
 
 import numpy as np
 from sparse_gp import SparseGP
-from parameterization.param import Param
-from ..inference.latent_function_inference import var_dtc
-from .. import likelihoods
-from parameterization.variational import VariationalPosterior
 from ..inference.latent_function_inference.var_dtc_parallel import update_gradients, VarDTC_minibatch
-from ..core.parameterization.parameter_core import OptimizationHandlable
 
 import logging
 logger = logging.getLogger("sparse gp mpi")
@@ -80,11 +75,15 @@ class SparseGP_MPI(SparseGP):
     #=====================================================
     
     @SparseGP.optimizer_array.setter
-    def optimizer_array(self, p):        
+    def optimizer_array(self, p):
         if self.mpi_comm != None:
             if self._IN_OPTIMIZATION_ and self.mpi_comm.rank==0:
                 self.mpi_comm.Bcast(np.int32(1),root=0)
             self.mpi_comm.Bcast(p, root=0)
+            
+        from ..util.debug import checkFinite
+        checkFinite(p, 'optimizer_array')
+        
         SparseGP.optimizer_array.fset(self,p)
         
     def optimize(self, optimizer=None, start=None, **kwargs):

GPy/inference/latent_function_inference/var_dtc_parallel.py

@@ -166,12 +166,16 @@ class VarDTC_minibatch(LatentFunctionInference):
         # Compute Common Components
         #======================================================================
         
+        from ...util.debug import checkFullRank
+        
         Kmm = kern.K(Z).copy()
         diag.add(Kmm, self.const_jitter)
+        checkFullRank(Kmm)
         Lm = jitchol(Kmm)
         
         LmInvPsi2LmInvT = backsub_both_sides(Lm,psi2_full,transpose='right')
         Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
+        checkFullRank(Lambda)
         LL = jitchol(Lambda)
         LL = np.dot(Lm,LL)
         b,_ = dtrtrs(LL, psi1Y_full.T)

GPy/util/debug.py

@@ -0,0 +1,36 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+"""
+The module for some general debug tools
+"""
+
+import numpy as np
+
+def checkFinite(arr, name=None):
+    if name is None:
+        name = 'Array with ID['+str(id(arr))+']'
+        
+    if np.any(np.logical_not(np.isfinite(arr))):
+        idx = np.where(np.logical_not(np.isfinite(arr)))[0]
+        print name+' at indices '+str(idx)+' have not finite values: '+str(arr[idx])+'!'
+        return False
+    return True
+
+def checkFullRank(m, tol=1e-10, name=None, force_check=False):
+    if name is None:
+        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:
+        print 'The size of '+name+'is too big to check (>=10000)!'
+        return True
+    
+    s = np.linalg.eigvals(m)
+    
+    if s.min()/s.max()<tol:
+        print name+' is close to singlar!'
+        print 'The eigen values of '+name+' is '+str(s)
+        return False
+    return True
+