# Source code for openmdao.surrogate_models.response_surface


"""
Surrogate Model based on second order response surface equations.
"""

from numpy import zeros, einsum
from numpy.dual import lstsq
from openmdao.surrogate_models.surrogate_model import SurrogateModel
from six.moves import range

[docs]class ResponseSurface(SurrogateModel):
"""
Surrogate Model based on second order response surface equations.

Attributes
----------
betas : ndarray
Vector of response surface equation coefficients.
m : int
Number of training points.
n : int
Number of independent variables.
"""

[docs]    def __init__(self):
"""
Initialize all attributes.
"""
super(ResponseSurface, self).__init__()

self.m = 0  # number of training points
self.n = 0  # number of independents
# vector of response surface equation coefficients
self.betas = zeros(0)

[docs]    def train(self, x, y):
"""
Calculate response surface equation coefficients using least squares regression.

Parameters
----------
x : array-like
Training input locations
y : array-like
Model responses at given inputs.
"""
super(ResponseSurface, self).train(x, y)

m = self.m = x.shape
n = self.n = x.shape

X = zeros((m, ((n + 1) * (n + 2)) // 2))

# Modify X to include constant, squared terms and cross terms

# Constant Terms
X[:, 0] = 1.0

# Linear Terms
X[:, 1:n + 1] = x

X_offset = X[:, n + 1:]
for i in range(n):
# Z = einsum('i,ij->ij', X, Y) is equivalent to, but much faster and
# memory efficient than, diag(X).dot(Y) for vector X and 2D array Y.
# I.e. Z[i,j] = X[i]*Y[i,j]
X_offset[:, :n - i] = einsum('i,ij->ij', x[:, i], x[:, i:])
X_offset = X_offset[:, n - i:]

# Determine response surface equation coefficients (betas) using least
# squares
self.betas, rs, r, s = lstsq(X, y)

[docs]    def predict(self, x):
"""
Calculate predicted value of response based on the current response surface model.

Parameters
----------
x : array-like
Point at which the surrogate is evaluated.

Returns
-------
float
Predicted response.
"""
super(ResponseSurface, self).predict(x)

n = x.size

X = zeros(((self.n + 1) * (self.n + 2)) // 2)

# Modify X to include constant, squared terms and cross terms

# Constant Terms
X = 1.0

# Linear Terms
X[1:n + 1] = x

X_offset = X[n + 1:]
for i in range(n):
X_offset[:n - i] = x[i] * x[i:]
X_offset = X_offset[n - i:]

# Predict new_y using X and betas
return X.dot(self.betas)

[docs]    def linearize(self, x):
"""
Calculate the jacobian of the Kriging surface at the requested point.

Parameters
----------
x : array-like
Point at which the surrogate Jacobian is evaluated.

Returns
-------
ndarray
Jacobian of surrogate output wrt inputs.
"""
n = self.n
betas = self.betas

x = x.flat

jac = betas[1:n + 1, :].copy()
beta_offset = betas[n + 1:, :]
for i in range(n):
jac[i, :] += x[i:].dot(beta_offset[:n - i, :])
jac[i:, :] += x[i] * beta_offset[:n - i, :]
beta_offset = beta_offset[n - i:, :]

return jac.T