diff --git a/.gitignore b/.gitignore index 73431343..60866848 100644 --- a/.gitignore +++ b/.gitignore @@ -36,3 +36,6 @@ nosetests.xml #vim *.swp + +#bfgs optimiser leaves this lying around +iterate.dat diff --git a/GPy/core/model.py b/GPy/core/model.py index 6e0dea60..46164f14 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -227,8 +227,8 @@ class model(parameterised): start = self.extract_param() optimizer = optimization.get_optimizer(optimizer) - opt = optimizer(start, f_fp=f_fp,f=f,fp=fp, **kwargs) - opt.run() + opt = optimizer(start, model = self, **kwargs) + opt.run(f_fp=f_fp, f=f, fp=fp) self.optimization_runs.append(opt) self.expand_param(opt.x_opt) diff --git a/GPy/inference/optimization.py b/GPy/inference/optimization.py index c96d2a3b..4cf56b69 100644 --- a/GPy/inference/optimization.py +++ b/GPy/inference/optimization.py @@ -19,15 +19,12 @@ class Optimizer(): :param messages: print messages from the optimizer? :type messages: (True | False) :param max_f_eval: maximum number of function evaluations - + :rtype: optimizer object. - - """ - def __init__(self, x_init, f_fp, f, fp , messages=False, max_f_eval=1e4, ftol=None, gtol=None, xtol=None): + + """ + def __init__(self, x_init, messages=False, model = None, max_f_eval=1e4, ftol=None, gtol=None, xtol=None): self.opt_name = None - self.f_fp = f_fp - self.f = f - self.fp = fp self.x_init = x_init self.messages = messages self.f_opt = None @@ -40,14 +37,15 @@ class Optimizer(): self.xtol = xtol self.gtol = gtol self.ftol = ftol - - def run(self): + self.model = model + + def run(self, **kwargs): start = dt.datetime.now() - self.opt() + self.opt(**kwargs) end = dt.datetime.now() self.time = str(end-start) - def opt(self): + def opt(self, f_fp = None, f = None, fp = None): raise NotImplementedError, "this needs to be implemented to use the optimizer class" def plot(self): @@ -71,7 +69,7 @@ class opt_tnc(Optimizer): Optimizer.__init__(self, *args, **kwargs) self.opt_name = "TNC (Scipy implementation)" - def opt(self): + def opt(self, f_fp = None, f = None, fp = None): """ Run the TNC optimizer @@ -79,7 +77,7 @@ class opt_tnc(Optimizer): tnc_rcstrings = ['Local minimum', 'Converged', 'XConverged', 'Maximum number of f evaluations reached', 'Line search failed', 'Function is constant'] - assert self.f_fp != None, "TNC requires f_fp" + assert f_fp != None, "TNC requires f_fp" opt_dict = {} if self.xtol is not None: @@ -89,10 +87,10 @@ class opt_tnc(Optimizer): if self.gtol is not None: opt_dict['pgtol'] = self.gtol - opt_result = optimize.fmin_tnc(self.f_fp, self.x_init, messages = self.messages, + opt_result = optimize.fmin_tnc(f_fp, self.x_init, messages = self.messages, maxfun = self.max_f_eval, **opt_dict) self.x_opt = opt_result[0] - self.f_opt = self.f_fp(self.x_opt)[0] + self.f_opt = f_fp(self.x_opt)[0] self.funct_eval = opt_result[1] self.status = tnc_rcstrings[opt_result[2]] @@ -101,14 +99,14 @@ class opt_lbfgsb(Optimizer): Optimizer.__init__(self, *args, **kwargs) self.opt_name = "L-BFGS-B (Scipy implementation)" - def opt(self): + def opt(self, f_fp = None, f = None, fp = None): """ Run the optimizer """ rcstrings = ['Converged', 'Maximum number of f evaluations reached', 'Error'] - assert self.f_fp != None, "BFGS requires f_fp" + assert f_fp != None, "BFGS requires f_fp" if self.messages: iprint = 1 @@ -123,10 +121,10 @@ class opt_lbfgsb(Optimizer): if self.gtol is not None: opt_dict['pgtol'] = self.gtol - opt_result = optimize.fmin_l_bfgs_b(self.f_fp, self.x_init, iprint = iprint, + opt_result = optimize.fmin_l_bfgs_b(f_fp, self.x_init, iprint = iprint, maxfun = self.max_f_eval, **opt_dict) self.x_opt = opt_result[0] - self.f_opt = self.f_fp(self.x_opt)[0] + self.f_opt = f_fp(self.x_opt)[0] self.funct_eval = opt_result[2]['funcalls'] self.status = rcstrings[opt_result[2]['warnflag']] @@ -135,7 +133,7 @@ class opt_simplex(Optimizer): Optimizer.__init__(self, *args, **kwargs) self.opt_name = "Nelder-Mead simplex routine (via Scipy)" - def opt(self): + def opt(self, f_fp = None, f = None, fp = None): """ The simplex optimizer does not require gradients. """ @@ -150,7 +148,7 @@ class opt_simplex(Optimizer): if self.gtol is not None: print "WARNING: simplex doesn't have an gtol arg, so I'm going to ignore it" - opt_result = optimize.fmin(self.f, self.x_init, (), disp = self.messages, + opt_result = optimize.fmin(f, self.x_init, (), disp = self.messages, maxfun = self.max_f_eval, full_output=True, **opt_dict) self.x_opt = opt_result[0] diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index be3c902f..ead4d0fc 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -2,5 +2,5 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) -from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, rbf_ARD, spline, Brownian, linear_ARD +from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, rbf_ARD, spline, Brownian, linear_ARD, rbf_sympy from kern import kern diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 99b0aab3..7bbac967 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -21,6 +21,7 @@ from Brownian import Brownian as Brownianpart #TODO these s=constructors are not as clean as we'd like. Tidy the code up #using meta-classes to make the objects construct properly wthout them. + def rbf(D,variance=1., lengthscale=1.): """ Construct an RBF kernel @@ -170,3 +171,20 @@ def Brownian(D,variance=1.): """ part = Brownianpart(D,variance) return kern(D, [part]) + +import sympy as sp +from sympykern import spkern +from sympy.parsing.sympy_parser import parse_expr + +def rbf_sympy(D,variance=1., lengthscale=1.): + """ + Radial Basis Function covariance. + """ + X = [sp.var('x%i'%i) for i in range(D)] + Z = [sp.var('z%i'%i) for i in range(D)] + rbf_variance = sp.var('rbf_variance',positive=True) + rbf_lengthscale = sp.var('rbf_lengthscale',positive=True) + dist_string = ' + '.join(['(x%i-z%i)**2'%(i,i) for i in range(D)]) + dist = parse_expr(dist_string) + f = rbf_variance*sp.exp(-dist/(2*rbf_lengthscale**2)) + return kern(D,[spkern(D,f,np.array([variance,lengthscale]))]) diff --git a/GPy/kern/sympy_helpers.cpp b/GPy/kern/sympy_helpers.cpp new file mode 100644 index 00000000..2af4737a --- /dev/null +++ b/GPy/kern/sympy_helpers.cpp @@ -0,0 +1,10 @@ +#include +double DiracDelta(double x){ + if((x<0.000001) & (x>-0.000001))//go on, laught at my c++ skills + return 1.0; + else + return 0.0; +}; +double DiracDelta(double x,int foo){ + return 0.0; +}; diff --git a/GPy/kern/sympy_helpers.h b/GPy/kern/sympy_helpers.h new file mode 100644 index 00000000..29244eca --- /dev/null +++ b/GPy/kern/sympy_helpers.h @@ -0,0 +1,3 @@ +#include +double DiracDelta(double x); +double DiracDelta(double x, int foo); diff --git a/GPy/kern/sympykern.py b/GPy/kern/sympykern.py new file mode 100644 index 00000000..4d912dc8 --- /dev/null +++ b/GPy/kern/sympykern.py @@ -0,0 +1,286 @@ +import numpy as np +import sympy as sp +from sympy.utilities.codegen import codegen +from sympy.core.cache import clear_cache +from scipy import weave +import re +import os +import sys +current_dir = os.path.dirname(os.path.abspath(os.path.dirname(__file__))) +import tempfile +import pdb +from kernpart import kernpart + +class spkern(kernpart): + """ + A kernel object, where all the hard work in done by sympy. + + :param k: the covariance function + :type k: a positive definite sympy function of x1, z1, x2, z2... + + To construct a new sympy kernel, you'll need to define: + - a kernel function using a sympy object. Ensure that the kernel is of the form k(x,z). + - that's it! we'll extract the variables from the function k. + + Note: + - to handle multiple inputs, call them x1, z1, etc + - to handle multpile correlated outputs, you'll need to define each covariance function and 'cross' variance function. TODO + """ + def __init__(self,D,k,param=None): + self.name='sympykern' + self._sp_k = k + sp_vars = [e for e in k.atoms() if e.is_Symbol] + self._sp_x= sorted([e for e in sp_vars if e.name[0]=='x'],key=lambda x:int(x.name[1:])) + self._sp_z= sorted([e for e in sp_vars if e.name[0]=='z'],key=lambda z:int(z.name[1:])) + assert all([x.name=='x%i'%i for i,x in enumerate(self._sp_x)]) + assert all([z.name=='z%i'%i for i,z in enumerate(self._sp_z)]) + assert len(self._sp_x)==len(self._sp_z) + self.D = len(self._sp_x) + assert self.D == D + self._sp_theta = sorted([e for e in sp_vars if not (e.name[0]=='x' or e.name[0]=='z')],key=lambda e:e.name) + self.Nparam = len(self._sp_theta) + + #deal with param + if param is None: + param = np.ones(self.Nparam) + assert param.size==self.Nparam + self.set_param(param) + + #Differentiate! + self._sp_dk_dtheta = [sp.diff(k,theta) for theta in self._sp_theta] + self._sp_dk_dx = [sp.diff(k,xi) for xi in self._sp_x] + #self._sp_dk_dz = [sp.diff(k,zi) for zi in self._sp_z] + + #self.compute_psi_stats() + self._gen_code() + + self.weave_kwargs = {\ + 'support_code':self._function_code,\ + 'include_dirs':[tempfile.gettempdir(), os.path.join(current_dir,'kern/')],\ + 'headers':['"sympy_helpers.h"'],\ + 'sources':[os.path.join(current_dir,"kern/sympy_helpers.cpp")],\ + #'extra_compile_args':['-ftree-vectorize', '-mssse3', '-ftree-vectorizer-verbose=5'],\ + 'extra_compile_args':[],\ + 'extra_link_args':['-lgomp'],\ + 'verbose':True} + + def __add__(self,other): + return spkern(self._sp_k+other._sp_k) + + def compute_psi_stats(self): + #define some normal distributions + mus = [sp.var('mu%i'%i,real=True) for i in range(self.D)] + Ss = [sp.var('S%i'%i,positive=True) for i in range(self.D)] + normals = [(2*sp.pi*Si)**(-0.5)*sp.exp(-0.5*(xi-mui)**2/Si) for xi, mui, Si in zip(self._sp_x, mus, Ss)] + + #do some integration! + #self._sp_psi0 = ?? + self._sp_psi1 = self._sp_k + for i in range(self.D): + print 'perfoming integrals %i of %i'%(i+1,2*self.D) + sys.stdout.flush() + self._sp_psi1 *= normals[i] + self._sp_psi1 = sp.integrate(self._sp_psi1,(self._sp_x[i],-sp.oo,sp.oo)) + clear_cache() + self._sp_psi1 = self._sp_psi1.simplify() + + #and here's psi2 (eek!) + zprime = [sp.Symbol('zp%i'%i) for i in range(self.D)] + self._sp_psi2 = self._sp_k.copy()*self._sp_k.copy().subs(zip(self._sp_z,zprime)) + for i in range(self.D): + print 'perfoming integrals %i of %i'%(self.D+i+1,2*self.D) + sys.stdout.flush() + self._sp_psi2 *= normals[i] + self._sp_psi2 = sp.integrate(self._sp_psi2,(self._sp_x[i],-sp.oo,sp.oo)) + clear_cache() + self._sp_psi2 = self._sp_psi2.simplify() + + + def _gen_code(self): + #generate c functions from sympy objects + (foo_c,self._function_code),(foo_h,self._function_header) = \ + codegen([('k',self._sp_k)] \ + + [('dk_d%s'%x.name,dx) for x,dx in zip(self._sp_x,self._sp_dk_dx)]\ + #+ [('dk_d%s'%z.name,dz) for z,dz in zip(self._sp_z,self._sp_dk_dz)]\ + + [('dk_d%s'%theta.name,dtheta) for theta,dtheta in zip(self._sp_theta,self._sp_dk_dtheta)]\ + ,"C",'foobar',argument_sequence=self._sp_x+self._sp_z+self._sp_theta) + #put the header file where we can find it + f = file(os.path.join(tempfile.gettempdir(),'foobar.h'),'w') + f.write(self._function_header) + f.close() + + #get rid of derivatives of DiracDelta + self._function_code = re.sub('DiracDelta\(.+?,.+?\)','0.0',self._function_code) + + #Here's some code to do the looping for K + arglist = ", ".join(["X[i*D+%s]"%x.name[1:] for x in self._sp_x]\ + + ["Z[j*D+%s]"%z.name[1:] for z in self._sp_z]\ + + ["param[%i]"%i for i in range(self.Nparam)]) + + self._K_code =\ + """ + int i; + int j; + int N = target_array->dimensions[0]; + int M = target_array->dimensions[1]; + int D = X_array->dimensions[1]; + #pragma omp parallel for private(j) + for (i=0;idimensions[0]; + int D = X_array->dimensions[1]; + #pragma omp parallel for + for (i=0;idimensions[0]; + int M = partial_array->dimensions[1]; + int D = X_array->dimensions[1]; + #pragma omp parallel for private(j) + for (i=0;idimensions[0]; + int D = X_array->dimensions[1]; + for (i=0;idimensions[0]; + int M = partial_array->dimensions[1]; + int D = X_array->dimensions[1]; + #pragma omp parallel for private(j) + for (i=0;idimensions[0]; + int M = partial_array->dimensions[1]; + int D = X_array->dimensions[1]; + for (i=0;idimensions[0]; + int M = 0; + int D = X_array->dimensions[1]; + for (i=0;i