Added configuration file

this was done to solve the OpenMP problem on Windows/mac, but I think it
is useful in general. All unit tests pass except the sympy kern ones.

GPy/examples/dimensionality_reduction.py

@@ -26,7 +26,7 @@ def BGPLVM(seed=default_seed):
     lik = Gaussian(Y, normalize=True)
 
     k = GPy.kern.rbf_inv(Q, .5, np.ones(Q) * 2., ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
-    # k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+    # k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
     # k = GPy.kern.rbf(Q, ARD = False)  + GPy.kern.white(Q, 0.00001)
 
     m = GPy.models.BayesianGPLVM(lik, Q, kernel=k, num_inducing=num_inducing)

GPy/gpy_config.cfg

@@ -0,0 +1,7 @@
+# This is the configuration file for GPy
+
+[parallel]
+# Enable openmp support. This speeds up some computations, depending on the number
+# of cores available. Setting up a compiler with openmp support can be difficult on 
+# some platforms, hence this option.
+openmp=True

GPy/kern/parts/linear.py

@@ -7,6 +7,7 @@ import numpy as np
 from ...util.linalg import tdot
 from ...util.misc import fast_array_equal
 from scipy import weave
+from ...util.config import *
 
 class Linear(Kernpart):
     """
@@ -51,6 +52,26 @@ class Linear(Kernpart):
         self._Z, self._mu, self._S = np.empty(shape=(3, 1))
         self._X, self._X2, self._params = np.empty(shape=(3, 1))
 
+        # a set of optional args to pass to weave
+        weave_options_openmp = {'headers'           : ['<omp.h>'],
+                                'extra_compile_args': ['-fopenmp -O3'],
+                                'extra_link_args'   : ['-lgomp'],
+                                'libraries': ['gomp']}
+        weave_options_noopenmp = {'extra_compile_args': ['-O3']}
+
+
+        if config.getboolean('parallel', 'openmp'):
+            self.weave_options = weave_options_openmp
+            self.weave_support_code =  """
+            #include <omp.h>
+            #include <math.h>
+            """
+        else:
+            self.weave_options = weave_options_noopenmp
+            self.weave_support_code = """
+            #include <math.h>
+            """
+
     def _get_params(self):
         return self.variances
 
@@ -190,11 +211,17 @@ class Linear(Kernpart):
         #target_mu_dummy += (dL_dpsi2[:, :, :, None] * muAZZA).sum(1).sum(1)
         #target_S_dummy += (dL_dpsi2[:, :, :, None] * self.ZA[None, :, None, :] * self.ZA[None, None, :, :]).sum(1).sum(1)
 
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = "#pragma omp parallel for private(m,mm,q,qq,factor,tmp)"
+        else:
+            pragma_string = ''
+
         #Using weave, we can exploiut the symmetry of this problem:
         code = """
         int n, m, mm,q,qq;
         double factor,tmp;
-        #pragma omp parallel for private(m,mm,q,qq,factor,tmp)
+        %s
         for(n=0;n<N;n++){
           for(m=0;m<num_inducing;m++){
             for(mm=0;mm<=m;mm++){
@@ -218,19 +245,13 @@ class Linear(Kernpart):
             }
           }
         }
-        """
+        """ % pragma_string
-        support_code = """
-        #include <omp.h>
-        #include <math.h>
-        """
-        weave_options = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
-                         'extra_link_args'   : ['-lgomp']}
 
-        N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
+
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+        N,num_inducing,input_dim = int(mu.shape[0]),int(Z.shape[0]),int(mu.shape[1])
-                     arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
+        weave.inline(code, support_code=self.weave_support_code,
-                     type_converters=weave.converters.blitz,**weave_options)
+                    arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
+                    type_converters=weave.converters.blitz,**self.weave_options)
 
 
     def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target):
@@ -240,9 +261,15 @@ class Linear(Kernpart):
         #dummy_target += psi2_dZ.sum(0).sum(0)
 
         AZA = self.variances*self.ZAinner
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for private(n,mm,q)'
+        else:
+            pragma_string = ''
+
         code="""
         int n,m,mm,q;
-        #pragma omp parallel for private(n,mm,q)
+        %s
         for(m=0;m<num_inducing;m++){
           for(q=0;q<input_dim;q++){
             for(mm=0;mm<num_inducing;mm++){
@@ -252,22 +279,13 @@ class Linear(Kernpart):
             }
           }
         }
-        """
+        """ % pragma_string
-        support_code = """
-        #include <omp.h>
-        #include <math.h>
-        """
-        weave_options = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
-                         'extra_link_args'   : ['-lgomp']}
 
-        N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
+
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+        N,num_inducing,input_dim = int(mu.shape[0]),int(Z.shape[0]),int(mu.shape[1])
+        weave.inline(code, support_code=self.weave_support_code, 
                      arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
-                     type_converters=weave.converters.blitz,**weave_options)
+                     type_converters=weave.converters.blitz,**self.weave_options)
-
-
-
 
 
     #---------------------------------------#

GPy/kern/parts/rbf.py

@@ -7,6 +7,7 @@ import numpy as np
 from scipy import weave
 from ...util.linalg import tdot
 from ...util.misc import fast_array_equal
+from ...util.config import *
 
 class RBF(Kernpart):
     """
@@ -57,12 +58,27 @@ class RBF(Kernpart):
         self._X, self._X2, self._params = np.empty(shape=(3, 1))
 
         # a set of optional args to pass to weave
-        self.weave_options = {'headers'           : ['<omp.h>'],
+        weave_options_openmp = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
+                                'extra_compile_args': ['-fopenmp -O3'],
-                         'extra_link_args'   : ['-lgomp']}
+                                'extra_link_args'   : ['-lgomp'],
+                                'libraries': ['gomp']}
+        weave_options_noopenmp = {'extra_compile_args': ['-O3']}
 
 
 
+        if config.getboolean('parallel', 'openmp'):
+            self.weave_options = weave_options_openmp
+            self.weave_support_code =  """
+            #include <omp.h>
+            #include <math.h>
+            """
+        else:
+            self.weave_options = weave_options_noopenmp
+            self.weave_support_code = """
+            #include <math.h>
+            """
+
+
     def _get_params(self):
         return np.hstack((self.variance, self.lengthscale))
 
@@ -110,7 +126,7 @@ class RBF(Kernpart):
                   target(q+1) += var_len3(q)*tmp;
                 }
                 """
-                num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim
+                num_data, num_inducing, input_dim = int(X.shape[0]), int(X.shape[0]), int(self.input_dim)
                 weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
             else:
                 code = """
@@ -126,7 +142,7 @@ class RBF(Kernpart):
                   target(q+1) += var_len3(q)*tmp;
                 }
                 """
-                num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim
+                num_data, num_inducing, input_dim = int(X.shape[0]), int(X2.shape[0]), int(self.input_dim)
                 # [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)]
                 weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
         else:
@@ -287,10 +303,16 @@ class RBF(Kernpart):
             lengthscale2 = self.lengthscale2
         else:
             lengthscale2 = np.ones(input_dim) * self.lengthscale2
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for private(tmp)'
+        else:
+            pragma_string = ''
+
         code = """
         double tmp;
 
-        #pragma omp parallel for private(tmp)
+        %s
         for (int n=0; n<N; n++){
             for (int m=0; m<num_inducing; m++){
                for (int mm=0; mm<(m+1); mm++){
@@ -320,13 +342,20 @@ class RBF(Kernpart):
             }
         }
 
-        """
+        """ % pragma_string
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#include <omp.h>'
+        else:
+            pragma_string = ''
 
         support_code = """
-        #include <omp.h>
+        %s
         #include <math.h>
-        """
+        """ % pragma_string
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
+
+        N, num_inducing, input_dim = int(N), int(num_inducing), int(input_dim)
+        weave.inline(code, support_code=support_code,
                      arg_names=['N', 'num_inducing', 'input_dim', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'lengthscale2', '_psi2_denom', 'psi2_Zdist_sq', 'psi2_exponent', 'half_log_psi2_denom', 'psi2', 'variance_sq'],
                      type_converters=weave.converters.blitz, **self.weave_options)
 

GPy/kern/parts/rbf_inv.py

@@ -7,6 +7,8 @@ import numpy as np
 import hashlib
 from scipy import weave
 from ...util.linalg import tdot
+from ...util.config import *
+
 
 class RBFInv(RBF):
     """
@@ -58,11 +60,23 @@ class RBFInv(RBF):
         self._X, self._X2, self._params = np.empty(shape=(3, 1))
 
         # a set of optional args to pass to weave
-        self.weave_options = {'headers'           : ['<omp.h>'],
+        weave_options_openmp = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
+                                'extra_compile_args': ['-fopenmp -O3'],
-                         'extra_link_args'   : ['-lgomp']}
+                                'extra_link_args'   : ['-lgomp'],
-
+                                'libraries': ['gomp']}
+        weave_options_noopenmp = {'extra_compile_args': ['-O3']}
 
+        if config.getboolean('parallel', 'openmp'):
+            self.weave_options = weave_options_openmp
+            self.weave_support_code =  """
+            #include <omp.h>
+            #include <math.h>
+            """
+        else:
+            self.weave_options = weave_options_noopenmp
+            self.weave_support_code = """
+            #include <math.h>
+            """
 
     def _get_params(self):
         return np.hstack((self.variance, self.inv_lengthscale))
@@ -109,7 +123,7 @@ class RBFInv(RBF):
                   target(q+1) += var_len3(q)*tmp*(-len2(q));
                 }
                 """
-                num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim
+                num_data, num_inducing, input_dim = int(X.shape[0]), int(X.shape[0]), int(self.input_dim)
                 weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3', 'len2'], type_converters=weave.converters.blitz, **self.weave_options)
             else:
                 code = """
@@ -125,7 +139,7 @@ class RBFInv(RBF):
                   target(q+1) += var_len3(q)*tmp*(-len2(q));
                 }
                 """
-                num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim
+                num_data, num_inducing, input_dim = int(X.shape[0]), int(X2.shape[0]), int(self.input_dim)
                 # [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)]
                 weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3', 'len2'], type_converters=weave.converters.blitz, **self.weave_options)
         else:
@@ -133,7 +147,7 @@ class RBFInv(RBF):
 
     def dK_dX(self, dL_dK, X, X2, target):
         self._K_computations(X, X2)
-        if X2 is None:            
+        if X2 is None:
             _K_dist = 2*(X[:, None, :] - X[None, :, :])
         else:
             _K_dist = X[:, None, :] - X2[None, :, :] # don't cache this in _K_computations because it is high memory. If this function is being called, chances are we're not in the high memory arena.
@@ -263,8 +277,8 @@ class RBFInv(RBF):
             self._Z, self._mu, self._S = Z, mu, S
 
     def weave_psi2(self, mu, Zhat):
-        N, input_dim = mu.shape
+        N, input_dim = int(mu.shape[0]), int(mu.shape[1])
-        num_inducing = Zhat.shape[0]
+        num_inducing = int(Zhat.shape[0])
 
         mudist = np.empty((N, num_inducing, num_inducing, input_dim))
         mudist_sq = np.empty((N, num_inducing, num_inducing, input_dim))
@@ -279,10 +293,16 @@ class RBFInv(RBF):
             inv_lengthscale2 = self.inv_lengthscale2
         else:
             inv_lengthscale2 = np.ones(input_dim) * self.inv_lengthscale2
+
+        if config.getboolean('parallel', 'openmp'):
+            pragma_string = '#pragma omp parallel for private(tmp)'
+        else:
+            pragma_string = ''
+
         code = """
         double tmp;
 
-        #pragma omp parallel for private(tmp)
+        %s
         for (int n=0; n<N; n++){
             for (int m=0; m<num_inducing; m++){
                for (int mm=0; mm<(m+1); mm++){
@@ -312,13 +332,9 @@ class RBFInv(RBF):
             }
         }
 
-        """
+        """ % pragma_string
 
-        support_code = """
+        weave.inline(code, support_code=self.weave_support_code,
-        #include <omp.h>
-        #include <math.h>
-        """
-        weave.inline(code, support_code=support_code, libraries=['gomp'],
                      arg_names=['N', 'num_inducing', 'input_dim', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'inv_lengthscale2', '_psi2_denom', 'psi2_Zdist_sq', 'psi2_exponent', 'half_log_psi2_denom', 'psi2', 'variance_sq'],
                      type_converters=weave.converters.blitz, **self.weave_options)
 

GPy/util/config.py

@@ -0,0 +1,17 @@
+#
+# This loads the configuration
+#
+import ConfigParser
+import os
+config = ConfigParser.ConfigParser()
+
+user_file = os.path.join(os.getenv('HOME'),'.gpy_config.cfg')
+default_file = os.path.join('..','gpy_config.cfg')
+
+# 1. check if the user has a ~/.gpy_config.cfg
+if os.path.isfile(user_file):
+    config.read(user_file)
+else:
+    # 2. if not, use the default one
+    path = os.path.dirname(__file__)
+    config.read(os.path.join(path,default_file))

GPy/util/misc.py

@@ -3,6 +3,7 @@
 
 import numpy as np
 from scipy import weave
+from config import *
 
 def opt_wrapper(m, **kwargs):
     """
@@ -57,11 +58,18 @@ def kmm_init(X, m = 10):
     return X[inducing]
 
 def fast_array_equal(A, B):
+
+
+    if config.getboolean('parallel', 'openmp'):
+        pragma_string = '#pragma omp parallel for private(i, j)'
+    else:
+        pragma_string = ''
+
     code2="""
     int i, j;
     return_val = 1;
 
-    // #pragma omp parallel for private(i, j)
+    %s
     for(i=0;i<N;i++){
        for(j=0;j<D;j++){
           if(A(i, j) != B(i, j)){
@@ -70,13 +78,18 @@ def fast_array_equal(A, B):
           }
        }
     }
-    """
+    """ % pragma_string
+
+    if config.getboolean('parallel', 'openmp'):
+        pragma_string = '#pragma omp parallel for private(i, j, z)'
+    else:
+        pragma_string = ''
 
     code3="""
     int i, j, z;
     return_val = 1;
 
-    // #pragma omp parallel for private(i, j, z)
+    %s
     for(i=0;i<N;i++){
        for(j=0;j<D;j++){
          for(z=0;z<Q;z++){
@@ -87,20 +100,33 @@ def fast_array_equal(A, B):
           }
        }
     }
-    """
+    """ % pragma_string
+
+    if config.getboolean('parallel', 'openmp'):
+        pragma_string = '#include <omp.h>'
+    else:
+        pragma_string = ''
 
     support_code = """
-    // #include <omp.h>
+    %s
     #include <math.h>
-    """
+    """ % pragma_string
 
-    weave_options = {'headers'           : ['<omp.h>'],
-                     'extra_compile_args': ['-fopenmp -O3'],
-                     'extra_link_args'   : ['-lgomp']}
 
+    weave_options_openmp = {'headers'           : ['<omp.h>'],
+                            'extra_compile_args': ['-fopenmp -O3'],
+                            'extra_link_args'   : ['-lgomp'],
+                            'libraries': ['gomp']}
+    weave_options_noopenmp = {'extra_compile_args': ['-O3']}
+
+    if config.getboolean('parallel', 'openmp'):
+        weave_options = weave_options_openmp
+    else:
+        weave_options = weave_options_noopenmp
 
     value = False
 
+
     if (A == None) and (B == None):
         return True
     elif ((A == None) and (B != None)) or ((A != None) and (B == None)):
@@ -110,14 +136,12 @@ def fast_array_equal(A, B):
             N, D = [int(i) for i in A.shape]
             value = weave.inline(code2, support_code=support_code,
                                  arg_names=['A', 'B', 'N', 'D'],
-                                 type_converters=weave.converters.blitz)
+                                 type_converters=weave.converters.blitz, **weave_options)
-            # libraries=['gomp'], **weave_options)
         elif A.ndim == 3:
             N, D, Q = [int(i) for i in A.shape]
             value = weave.inline(code3, support_code=support_code,
                                  arg_names=['A', 'B', 'N', 'D', 'Q'],
-                                 type_converters=weave.converters.blitz)
+                                 type_converters=weave.converters.blitz, **weave_options)
-            #libraries=['gomp'], **weave_options)
         else:
             value = np.array_equal(A,B)