fix the SSGPLVM with MPI

GPy/core/parameterization/parameter_core.py

@@ -425,9 +425,6 @@ class Indexable(Nameable, Observable):
     def _connect_fixes(self):
         from ties_and_remappings import Tie
         self._ensure_fixes()
-#         for c, ind in self.constraints.iteritems():
-#             if c == __fixed__ or isinstance(c,Tie):
-#                 self._fixes_[ind] = FIXED
         [np.put(self._fixes_, ind, FIXED) for c, ind in self.constraints.iteritems()
             if c == __fixed__ or isinstance(c,Tie)]
         if np.all(self._fixes_): self._fixes_ = None

GPy/core/parameterization/ties_and_remappings.py

@@ -59,7 +59,7 @@ class Tie(Remapping):
         uvals = np.unique(vals)
         if len(uvals)==1:
             #all of the tied things are at the same value
-            if (self.value==uvals[0]).all():
+            if np.all(self.value==uvals[0]):
                 return # DO NOT DO ANY CHANGES IF THE TIED PART IS NOT CHANGED!
             self.value[...] = uvals[0]
         elif len(uvals)==2:
@@ -72,7 +72,7 @@ class Tie(Remapping):
     def parameters_changed(self):
         #ensure all out parameters have the correct value, as specified by our mapping
         index = self._highest_parent_.constraints[self]
-        if (self._highest_parent_.param_array[index]==self.value).all():
+        if np.all(self._highest_parent_.param_array[index]==self.value):
             return # STOP TRIGGER THE UPDATE LOOP MULTIPLE TIMES!!!
         self._highest_parent_.param_array[index] = self.mapping()
         [p.notify_observers(which=self) for p in self.tied_parameters]

GPy/core/parameterization/variational.py

@@ -150,6 +150,9 @@ class SpikeAndSlabPosterior(VariationalPosterior):
             n.parameters[dc['mean']._parent_index_] = dc['mean']
             n.parameters[dc['variance']._parent_index_] = dc['variance']
             n.parameters[dc['binary_prob']._parent_index_] = dc['binary_prob']
+            n._gradient_array_ = None
+            oversize = self.size - self.mean.size - self.variance.size
+            n.size = n.mean.size + n.variance.size + oversize
             n.ndim = n.mean.ndim
             n.shape = n.mean.shape
             n.num_data = n.mean.shape[0]

GPy/inference/latent_function_inference/var_dtc_parallel.py

@@ -328,7 +328,7 @@ def update_gradients(model, mpi_comm=None):
         X = model.X
     else:
         Y = model.Y_local
-        X = model.X_local
+        X = model.X[model.N_range[0]:model.N_range[1]]
 
     model._log_marginal_likelihood, dL_dKmm, model.posterior = model.inference_method.inference_likelihood(model.kern, X, model.Z, model.likelihood, Y)
     
@@ -350,7 +350,7 @@ def update_gradients(model, mpi_comm=None):
             if mpi_comm ==None:
                 X_slice = model.X[n_range[0]:n_range[1]]
             else:
-                X_slice = model.X[model.Y_range[0]+n_range[0]:model.Y_range[0]+n_range[1]]
+                X_slice = model.X[model.N_range[0]+n_range[0]:model.N_range[0]+n_range[1]]
             
             #gradients w.r.t. kernel
             model.kern.update_gradients_expectations(variational_posterior=X_slice, Z=model.Z, dL_dpsi0=grad_dict['dL_dpsi0'], dL_dpsi1=grad_dict['dL_dpsi1'], dL_dpsi2=grad_dict['dL_dpsi2'])
@@ -396,7 +396,7 @@ def update_gradients(model, mpi_comm=None):
         KL_div_all = np.array(KL_div)
         mpi_comm.Allreduce([np.float64(KL_div), MPI.DOUBLE], [KL_div_all, MPI.DOUBLE])
         KL_div = KL_div_all
-        [mpi_comm.Allgatherv([pp.copy(), MPI.DOUBLE], [pa, (model.Y_list*pa.shape[-1], None), MPI.DOUBLE]) for pp,pa in zip(model.get_X_gradients(X),model.get_X_gradients(model.X))]
+        [mpi_comm.Allgatherv([pp.copy(), MPI.DOUBLE], [pa, (model.N_list*pa.shape[-1], None), MPI.DOUBLE]) for pp,pa in zip(model.get_X_gradients(X),model.get_X_gradients(model.X))]
         from ...models import SSGPLVM
         if isinstance(model, SSGPLVM):
             grad_pi = np.array(model.variational_prior.pi.gradient)

GPy/kern/_src/psi_comp/ssrbf_psi_gpucomp.py

@@ -339,7 +339,7 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
                              'grad_l_gpu'               :gpuarray.empty((Q,),np.float64, order='F'),
                              'grad_mu_gpu'              :gpuarray.empty((N,Q,),np.float64, order='F'),
                              'grad_S_gpu'               :gpuarray.empty((N,Q,),np.float64, order='F'),
-                             'grad_gamma_gpu'               :gpuarray.empty((N,Q,),np.float64, order='F'),
+                             'grad_gamma_gpu'           :gpuarray.empty((N,Q,),np.float64, order='F'),
                              }
         else:
             assert N==self.gpuCache['mu_gpu'].shape[0]
@@ -356,8 +356,6 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
         self.gpuCache['S_gpu'].set(np.asfortranarray(S))
         self.gpuCache['gamma_gpu'].set(np.asfortranarray(gamma))
         N,Q = self.gpuCache['S_gpu'].shape
-        # t=self.g_compDenom(self.gpuCache['log_denom1_gpu'],self.gpuCache['log_denom2_gpu'],self.gpuCache['l_gpu'],self.gpuCache['S_gpu'], np.int32(N), np.int32(Q), block=(self.threadnum,1,1), grid=(self.blocknum,1),time_kernel=True)
-        # print 'g_compDenom '+str(t)
         self.g_compDenom.prepared_call((self.blocknum,1),(self.threadnum,1,1), self.gpuCache['log_denom1_gpu'].gpudata,self.gpuCache['log_denom2_gpu'].gpudata,self.gpuCache['log_gamma_gpu'].gpudata,self.gpuCache['log_gamma1_gpu'].gpudata,self.gpuCache['gamma_gpu'].gpudata,self.gpuCache['l_gpu'].gpudata,self.gpuCache['S_gpu'].gpudata, np.int32(N), np.int32(Q))
         
     def reset_derivative(self):
@@ -393,7 +391,6 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
         Z_gpu = self.gpuCache['Z_gpu']
         mu_gpu = self.gpuCache['mu_gpu']
         S_gpu = self.gpuCache['S_gpu']
-        gamma_gpu = self.gpuCache['gamma_gpu']
         log_denom1_gpu = self.gpuCache['log_denom1_gpu']
         log_denom2_gpu = self.gpuCache['log_denom2_gpu']
         log_gamma_gpu = self.gpuCache['log_gamma_gpu']
@@ -403,10 +400,6 @@ class PSICOMP_SSRBF_GPU(PSICOMP_RBF):
         psi0[:] = variance
         self.g_psi1computations.prepared_call((self.blocknum,1),(self.threadnum,1,1),psi1_gpu.gpudata, log_denom1_gpu.gpudata, log_gamma_gpu.gpudata, log_gamma1_gpu.gpudata, np.float64(variance),l_gpu.gpudata,Z_gpu.gpudata,mu_gpu.gpudata,S_gpu.gpudata, np.int32(N), np.int32(M), np.int32(Q))
         self.g_psi2computations.prepared_call((self.blocknum,1),(self.threadnum,1,1),psi2_gpu.gpudata, psi2n_gpu.gpudata, log_denom2_gpu.gpudata, log_gamma_gpu.gpudata, log_gamma1_gpu.gpudata, np.float64(variance),l_gpu.gpudata,Z_gpu.gpudata,mu_gpu.gpudata,S_gpu.gpudata, np.int32(N), np.int32(M), np.int32(Q))
-        # t = self.g_psi1computations(psi1_gpu, log_denom1_gpu, np.float64(variance),l_gpu,Z_gpu,mu_gpu,S_gpu, np.int32(N), np.int32(M), np.int32(Q), block=(self.threadnum,1,1), grid=(self.blocknum,1),time_kernel=True)
-        # print 'g_psi1computations '+str(t)
-        # t = self.g_psi2computations(psi2_gpu, psi2n_gpu, log_denom2_gpu, np.float64(variance),l_gpu,Z_gpu,mu_gpu,S_gpu, np.int32(N), np.int32(M), np.int32(Q), block=(self.threadnum,1,1), grid=(self.blocknum,1),time_kernel=True)
-        # print 'g_psi2computations '+str(t)
         
         if self.GPU_direct:
             return psi0, psi1_gpu, psi2_gpu

GPy/models/ss_gplvm.py

@@ -45,7 +45,6 @@ class SSGPLVM(SparseGP):
         gamma[:] = 0.5 + 0.1 * np.random.randn(X.shape[0], input_dim)
         gamma[gamma>1.-1e-9] = 1.-1e-9
         gamma[gamma<1e-9] = 1e-9
-        gamma[:] = 0.5
                 
         if Z is None:
             Z = np.random.permutation(X.copy())[:num_inducing]
@@ -73,19 +72,18 @@ class SSGPLVM(SparseGP):
         self.add_parameter(self.X, index=0)
         self.add_parameter(self.variational_prior)
                 
-        if self.group_spike:
-            [self.X.gamma[:,i].tie('tieGamma'+str(i)) for i in xrange(self.X.gamma.shape[1])] # Tie columns together
-        
         if mpi_comm != None:
             from ..util.mpi import divide_data
-            Y_start, Y_end, Y_list = divide_data(Y.shape[0], mpi_comm)
+            N_start, N_end, N_list = divide_data(Y.shape[0], mpi_comm)
-            self.Y_local = self.Y[Y_start:Y_end]
+            self.N_range = (N_start, N_end)
-            self.X_local = self.X[Y_start:Y_end]
+            self.N_list = np.array(N_list)
-            self.Y_range = (Y_start, Y_end)
+            self.Y_local = self.Y[N_start:N_end]
-            self.Y_list = np.array(Y_list)
+            print 'MPI RANK: '+str(self.mpi_comm.rank)+' with datasize: '+str(self.N_range)
-            print self.mpi_comm.rank, self.Y_range
             mpi_comm.Bcast(self.param_array, root=0)
         
+        if self.group_spike:
+            [self.X.gamma[:,i].tie('tieGamma'+str(i)) for i in xrange(self.X.gamma.shape[1])] # Tie columns together
+        
     def set_X_gradients(self, X, X_grad):
         """Set the gradients of the posterior distribution of X in its specific form."""
         X.mean.gradient, X.variance.gradient, X.binary_prob.gradient = X_grad
@@ -124,9 +122,10 @@ class SSGPLVM(SparseGP):
         dc = super(SSGPLVM, self).__getstate__()
         dc['mpi_comm'] = None
         if self.mpi_comm != None:
+            del dc['N_range']
+            del dc['N_list']
             del dc['Y_local']
             del dc['X_local']
-            del dc['Y_range']
         return dc
  
     def __setstate__(self, state):

GPy/util/gpu_init.py

@@ -14,3 +14,9 @@ try:
     initSuccess = True
 except:
     initSuccess = False
+    
+def initGPU(gpuid=None):
+    if gpuid==None:
+        return pycuda.tools.make_default_context()
+    else:
+        return pycuda.driver.Device(gpuid).make_context()