Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make kriging surrogate picklable #154

Merged
merged 5 commits into from
Jun 13, 2019
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
make kriging poly, corr and mfk rho_regr options an enumerate
  • Loading branch information
relf committed Jun 13, 2019
commit 412487270b24a546be5427b9334e98e37ee83c59
40 changes: 20 additions & 20 deletions smt/extensions/mfk.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,9 @@ def _initialize(self):
super(MFK, self)._initialize()
declare = self.options.declare

declare('rho_regr', 'constant',types=FunctionType,\
values=('constant', 'linear', 'quadratic'), desc='regr. term')
declare('theta0', None, types=(list, np.ndarray), \
declare('rho_regr', 'constant', types=str,\
values=('constant', 'linear', 'quadratic'), desc='Regression function type for rho')
declare('theta0', types=(list, np.ndarray), \
desc='Initial hyperparameters')
declare('optim_var', False, types = bool, \
values = (True, False), \
Expand Down Expand Up @@ -145,12 +145,12 @@ def _new_train(self):


# Regression matrix and parameters
self.F_all[lvl] = self.options['poly'](self.X_norma)
self.F_all[lvl] = self._regression_types[self.options['poly']](self.X_norma)
self.p_all[lvl] = self.F_all[lvl].shape[1]

# Concatenate the autoregressive part for levels > 0
if lvl > 0:
F_rho = self.options['rho_regr'](self.X_norma)
F_rho = self._regression_types[self.options['rho_regr']](self.X_norma)
self.q_all[lvl] = F_rho.shape[1]
self.F_all[lvl] = np.hstack((F_rho*np.dot(self._predict_intermediate_values(self.X_norma, lvl, descale=False),
np.ones((1,self.q_all[lvl]))), self.F_all[lvl]))
Expand Down Expand Up @@ -194,7 +194,7 @@ def _new_train(self):
self._new_train()

def _componentwise_distance(self,dx,opt=0):
d = componentwise_distance(dx,self.options['corr'].__name__,
d = componentwise_distance(dx, self.options['corr'],
self.nx)
return d

Expand Down Expand Up @@ -225,8 +225,8 @@ def _predict_intermediate_values(self, X, lvl, descale = True):
if descale :
X = (X - self.X_mean) / self.X_std
## X = (X - self.X_mean[0]) / self.X_std[0]
f = self.options['poly'](X)
f0 = self.options['poly'](X)
f = self._regression_types[self.options['poly']](X)
f0 = self._regression_types[self.options['poly']](X)
dx = manhattan_distances(X, Y=self.X_norma_all[0], sum_over_features=False)
d = self._componentwise_distance(dx)
# Get regression function and correlation
Expand All @@ -236,7 +236,7 @@ def _predict_intermediate_values(self, X, lvl, descale = True):
beta = self.optimal_par[0]['beta']
Ft = solve_triangular(C, F, lower=True)
yt = solve_triangular(C, self.y_norma_all[0], lower=True)
r_ = self.options['corr'](self.optimal_theta[0], d).reshape(n_eval, self.nt_all[0])
r_ = self._correlation_types[self.options['corr']](self.optimal_theta[0], d).reshape(n_eval, self.nt_all[0])
gamma = self.optimal_par[0]['gamma']


Expand All @@ -247,10 +247,10 @@ def _predict_intermediate_values(self, X, lvl, descale = True):
for i in range(1,lvl):
F = self.F_all[i]
C = self.optimal_par[i]['C']
g = self.options['rho_regr'](X)
g = self._regression_types[self.options['rho_regr']](X)
dx = manhattan_distances(X, Y=self.X_norma_all[i], sum_over_features=False)
d = self._componentwise_distance(dx)
r_ = self.options['corr'](self.optimal_theta[i], d).reshape(n_eval, self.nt_all[i])
r_ = self._correlation_types[self.options['corr']](self.optimal_theta[i], d).reshape(n_eval, self.nt_all[i])
f = np.vstack((g.T*mu[:,i-1], f0.T))
Ft = solve_triangular(C, F, lower=True)
yt = solve_triangular(C, self.y_norma_all[i], lower=True)
Expand Down Expand Up @@ -339,7 +339,7 @@ def predict_variances_all_levels(self, X):
beta = self.optimal_par[0]['beta']
Ft = solve_triangular(C, F, lower=True)
yt = solve_triangular(C, self.y_norma_all[0], lower=True)
r_ = self.options['corr'](self.optimal_theta[0], d).reshape(n_eval, self.nt_all[0])
r_ = self._correlation_types[self.options['corr']](self.optimal_theta[0], d).reshape(n_eval, self.nt_all[0])
gamma = self.optimal_par[0]['gamma']

# Scaled predictor
Expand All @@ -360,7 +360,7 @@ def predict_variances_all_levels(self, X):
for i in range(1,nlevel):
F = self.F_all[i]
C = self.optimal_par[i]['C']
g = self.options['rho_regr'](X)
g = self._regression_types[self.options['rho_regr']](X)
dx = manhattan_distances(X, Y=self.X_norma_all[i], sum_over_features=False)
d = self._componentwise_distance(dx)
r_ = self.options['corr'](self.optimal_theta[i], d).reshape(n_eval, self.nt_all[i])
Expand Down Expand Up @@ -424,28 +424,28 @@ def _predict_derivatives(self, x, kx):

dy_dx = np.zeros((n_eval, lvl))

if self.options['corr'].__name__ != 'squar_exp':
if self.options['corr'] != 'squar_exp':
raise ValueError(
'The derivative is only available for square exponential kernel')
if self.options['poly'].__name__ == 'constant':
if self.options['poly'] == 'constant':
df = np.zeros([n_eval,1])
elif self.options['poly'].__name__ == 'linear':
elif self.options['poly'] == 'linear':
df = np.zeros((n_eval, self.nx + 1))
df[:,1:] = 1
else:
raise ValueError(
'The derivative is only available for ordinary kriging or '+
'universal kriging using a linear trend')
df0 = copy.deepcopy(df)
if self.options['rho_regr'].__name__ != 'constant' :
if self.options['rho_regr'] != 'constant' :
raise ValueError(
'The derivative is only available for regression rho constant')
# Get pairwise componentwise L1-distances to the input training set
dx = manhattan_distances(x, Y=self.X_norma_all[0], sum_over_features=
False)
d = self._componentwise_distance(dx)
# Compute the correlation function
r_ = self.options['corr'](self.optimal_theta[0], d).reshape(n_eval,self.nt_all[0])
r_ = self._correlation_types[self.options['corr']](self.optimal_theta[0], d).reshape(n_eval,self.nt_all[0])

# Beta and gamma = R^-1(y-FBeta)
beta = self.optimal_par[0]['beta']
Expand All @@ -463,10 +463,10 @@ def _predict_derivatives(self, x, kx):
for i in range(1,lvl):
F = self.F_all[i]
C = self.optimal_par[i]['C']
g = self.options['rho_regr'](x)
g = self._regression_types[self.options['rho_regr']](x)
dx = manhattan_distances(x, Y=self.X_norma_all[i], sum_over_features=False)
d = self._componentwise_distance(dx)
r_ = self.options['corr'](self.optimal_theta[i], d).reshape(n_eval, self.nt_all[i])
r_ = self._correlation_types[self.options['corr']](self.optimal_theta[i], d).reshape(n_eval, self.nt_all[i])
df = np.vstack((g.T*dy_dx[:,i-1], df0.T))

Ft = solve_triangular(C, F, lower=True)
Expand Down
2 changes: 1 addition & 1 deletion smt/surrogate_models/gekpls.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,6 @@ def _compute_pls(self,X,y):

def _componentwise_distance(self,dx,opt=0):

d = componentwise_distance_PLS(dx,self.options['corr'].__name__,
d = componentwise_distance_PLS(dx,self.options['corr'],
self.options['n_comp'],self.coeff_pls)
return d
2 changes: 1 addition & 1 deletion smt/surrogate_models/kpls.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,6 @@ def _compute_pls(self,X,y):
return X,y

def _componentwise_distance(self,dx,opt=0):
d = componentwise_distance_PLS(dx,self.options['corr'].__name__,
d = componentwise_distance_PLS(dx, self.options['corr'],
self.options['n_comp'],self.coeff_pls)
return d
4 changes: 2 additions & 2 deletions smt/surrogate_models/kplsk.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,9 +35,9 @@ def _compute_pls(self,X,y):
def _componentwise_distance(self,dx,opt=0):
if opt == 0:
# Kriging step
d = componentwise_distance(dx,self.options['corr'].__name__,self.nx)
d = componentwise_distance(dx,self.options['corr'],self.nx)
else:
# KPLS step
d = componentwise_distance_PLS(dx,self.options['corr'].__name__,
d = componentwise_distance_PLS(dx,self.options['corr'],
self.options['n_comp'],self.coeff_pls)
return d
13 changes: 6 additions & 7 deletions smt/surrogate_models/krg.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
'''
"""
Author: Dr. Mohamed A. Bouhlel <[email protected]>

This package is distributed under New BSD license.
'''
"""

from __future__ import division
import warnings
Expand All @@ -15,13 +15,12 @@
The kriging class.
"""

class KRG(KrgBased):

class KRG(KrgBased):
def _initialize(self):
super(KRG, self)._initialize()
self.name = 'Kriging'
self.name = "Kriging"

def _componentwise_distance(self,dx,opt=0):
d = componentwise_distance(dx,self.options['corr'].__name__,
self.nx)
def _componentwise_distance(self, dx, opt=0):
d = componentwise_distance(dx, self.options["corr"], self.nx)
return d
62 changes: 15 additions & 47 deletions smt/surrogate_models/krg_based.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,6 @@
TODO:
- fail_iteration and nb_iter_max to remove from options
- define outputs['sol'] = self.sol

- debug _train: self_pkl = pickle.dumps(obj)
cPickle.PicklingError: Can't pickle <type 'function'>: attribute lookup __builtin__.function failed

"""

from __future__ import division
Expand Down Expand Up @@ -50,11 +46,11 @@ def _initialize(self):
super(KrgBased, self)._initialize()
declare = self.options.declare
supports = self.supports
declare('poly', 'constant',types=FunctionType,values=('constant', 'linear', 'quadratic'),
desc='regr. term')
declare('corr', 'squar_exp', types=FunctionType,values=('abs_exp', 'squar_exp'),
desc='type of corr. func.')
declare('data_dir', values=None, types=str,
declare('poly', 'constant', values=('constant', 'linear', 'quadratic'),
desc='Regression function type')
declare('corr', 'squar_exp', values=('abs_exp', 'squar_exp'),
desc='Correlation function type')
declare('data_dir', types=str,
desc='Directory for loading / saving cached data; None means do not save or load')
declare('theta0', [1e-2], types=(list, np.ndarray), desc='Initial hyperparameters')
self.name = 'KrigingBased'
Expand Down Expand Up @@ -92,7 +88,7 @@ def _new_train(self):
raise Exception("Multiple input features cannot have the same value.")

# Regression matrix and parameters
self.F = self.options['poly'](self.X_norma)
self.F = self._regression_types[self.options['poly']](self.X_norma)
n_samples_F = self.F.shape[0]
if self.F.ndim > 1:
p = self.F.shape[1]
Expand All @@ -112,16 +108,13 @@ def _train(self):
"""
Train the model
"""
"""
inputs = {'self': self}
with cached_operation(inputs, self.options['data_dir']) as outputs:
if outputs:
self.sol = outputs['sol']
else:
self._new_train()
#outputs['sol'] = self.sol
"""
self._new_train()

def _reduced_likelihood_function(self, theta):

Expand Down Expand Up @@ -183,7 +176,7 @@ def _reduced_likelihood_function(self, theta):
theta = tmp_var[:-1]
noise = tmp_var[-1]

r = self.options['corr'](theta, self.D).reshape(-1,1)
r = self._correlation_types[self.options['corr']](theta, self.D).reshape(-1,1)

R = np.eye(self.nt) * (1. + nugget+ noise)
R[self.ij[:, 0], self.ij[:, 1]] = r[:,0]
Expand Down Expand Up @@ -276,10 +269,10 @@ def _predict_values(self, x):
False)
d = self._componentwise_distance(dx)
# Compute the correlation function
r = self.options['corr'](self.optimal_theta, d).reshape(n_eval,self.nt)
r = self._correlation_types[self.options['corr']](self.optimal_theta, d).reshape(n_eval,self.nt)
y = np.zeros(n_eval)
# Compute the regression function
f = self.options['poly'](x)
f = self._regression_types[self.options['poly']](x)
# Scaled predictor
y_ = np.dot(f, self.optimal_par['beta']) + np.dot(r,
self.optimal_par['gamma'])
Expand Down Expand Up @@ -314,14 +307,14 @@ def _predict_derivatives(self, x, kx):
False)
d = self._componentwise_distance(dx)
# Compute the correlation function
r = self.options['corr'](self.optimal_theta, d).reshape(n_eval,self.nt)
r = self._correlation_types[self.options['corr']](self.optimal_theta, d).reshape(n_eval,self.nt)

if self.options['corr'].__name__ != 'squar_exp':
if self.options['corr'] != 'squar_exp':
raise ValueError(
'The derivative is only available for square exponential kernel')
if self.options['poly'].__name__ == 'constant':
if self.options['poly'] == 'constant':
df = np.array([0])
elif self.options['poly'].__name__ == 'linear':
elif self.options['poly'] == 'linear':
df = np.zeros((self.nx + 1, self.nx))
df[1:,:] = 1
else:
Expand Down Expand Up @@ -493,7 +486,7 @@ def minus_reduced_likelihood_function(log10t):
if exit_function:
return best_optimal_rlf_value, best_optimal_par, best_optimal_theta

if self.options['corr'].__name__ == 'squar_exp':
if self.options['corr'] == 'squar_exp':
self.options['theta0'] = (best_optimal_theta*self.coeff_pls**2).sum(1)
else:
self.options['theta0'] = (best_optimal_theta*np.abs(self.coeff_pls)).sum(1)
Expand All @@ -508,24 +501,7 @@ def _check_param(self):
"""
This function check some parameters of the model.
"""
# Check regression model
if not callable(self.options['poly']):
if self.options['poly'] in self._regression_types:
self.options['poly'] = self._regression_types[
self.options['poly']]
else:
raise ValueError("regr should be one of %s or callable, "
"%s was given." % (self._regression_types.keys(),
self.options['poly']))
if self.name == 'MFK' and not callable(self.options['rho_regr']):
if self.options['rho_regr'] in self._regression_types:
self.options['rho_regr'] = self._regression_types[
self.options['rho_regr']]
else:
raise ValueError("rho_regr should be one of %s or callable, "
"%s was given." % (self._regression_types.keys(),
self.options['rho_regr']))

# FIXME: _check_param should be overriden in corresponding subclasses
if self.name in ['KPLS', 'KPLSK', 'GEKPLS']:
d = self.options['n_comp']
else:
Expand All @@ -538,14 +514,6 @@ def _check_param(self):
raise ValueError('the number of dim %s should be equal to the length of theta0 %s.' %
(d, len(self.options['theta0'])))

if not callable(self.options['corr']):
if self.options['corr'] in self._correlation_types:
self.options['corr'] = self._correlation_types[self.options['corr']]
else:
raise ValueError("corr should be one of %s or callable, "
"%s was given."% (self._correlation_types.keys(),
self.options['corr']))

if self.supports['training_derivatives']:
if not(1 in self.training_points[None]):
raise Exception('Derivative values are needed for using the GEKPLS model.')
Expand Down