diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 298607b6..bde249c8 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/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) diff --git a/GPy/gpy_config.cfg b/GPy/gpy_config.cfg new file mode 100644 index 00000000..8683f96c --- /dev/null +++ b/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 diff --git a/GPy/kern/parts/linear.py b/GPy/kern/parts/linear.py index ffcbcf5e..ab96bb31 100644 --- a/GPy/kern/parts/linear.py +++ b/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' : [''], + '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 + #include + """ + else: + self.weave_options = weave_options_noopenmp + self.weave_support_code = """ + #include + """ + 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 - #include - """ - weave_options = {'headers' : [''], - 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], - 'extra_link_args' : ['-lgomp']} + """ % pragma_string - N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1] - weave.inline(code, support_code=support_code, libraries=['gomp'], - arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'], - type_converters=weave.converters.blitz,**weave_options) + + 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','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 - #include - """ - weave_options = {'headers' : [''], - 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], - 'extra_link_args' : ['-lgomp']} + """ % pragma_string - 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) #---------------------------------------# diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index 855e2b71..585d687f 100644 --- a/GPy/kern/parts/rbf.py +++ b/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' : [''], - 'extra_compile_args': ['-fopenmp -O3'], # -march=native'], - 'extra_link_args' : ['-lgomp']} + weave_options_openmp = {'headers' : [''], + '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 + #include + """ + else: + self.weave_options = weave_options_noopenmp + self.weave_support_code = """ + #include + """ + + 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 + %s #include - """ - weave.inline(code, support_code=support_code, libraries=['gomp'], + """ % pragma_string + + 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) diff --git a/GPy/kern/parts/rbf_inv.py b/GPy/kern/parts/rbf_inv.py index 0433e96c..1cc05aaa 100644 --- a/GPy/kern/parts/rbf_inv.py +++ b/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' : [''], - 'extra_compile_args': ['-fopenmp -O3'], # -march=native'], - 'extra_link_args' : ['-lgomp']} - + weave_options_openmp = {'headers' : [''], + '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 + #include + """ + else: + self.weave_options = weave_options_noopenmp + self.weave_support_code = """ + #include + """ 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 - num_inducing = Zhat.shape[0] + N, input_dim = int(mu.shape[0]), int(mu.shape[1]) + 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 - #include - """ - weave.inline(code, support_code=support_code, libraries=['gomp'], + weave.inline(code, support_code=self.weave_support_code, 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) diff --git a/GPy/util/config.py b/GPy/util/config.py new file mode 100644 index 00000000..d2ed7543 --- /dev/null +++ b/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)) diff --git a/GPy/util/misc.py b/GPy/util/misc.py index 5866ecf9..d3f23b75 100644 --- a/GPy/util/misc.py +++ b/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 + %s #include - """ + """ % pragma_string - weave_options = {'headers' : [''], - 'extra_compile_args': ['-fopenmp -O3'], - 'extra_link_args' : ['-lgomp']} + weave_options_openmp = {'headers' : [''], + '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) - # libraries=['gomp'], **weave_options) + type_converters=weave.converters.blitz, **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) - #libraries=['gomp'], **weave_options) + type_converters=weave.converters.blitz, **weave_options) else: value = np.array_equal(A,B)