1import numpy as np 

2from scipy.optimize import minimize 

3from scipy.linalg import solve_triangular, cho_factor, cho_solve 

4from ase.optimize.gpmin.kernel import SquaredExponential 

5from ase.optimize.gpmin.prior import ZeroPrior 

6 

7 

8class GaussianProcess(): 

9 """Gaussian Process Regression 

10 It is recommended to be used with other Priors and Kernels from 

11 ase.optimize.gpmin 

12 

13 Parameters: 

14 

15 prior: Prior class, as in ase.optimize.gpmin.prior 

16 Defaults to ZeroPrior 

17 

18 kernel: Kernel function for the regression, as in 

19 ase.optimize.gpmin.kernel 

20 Defaults to the Squared Exponential kernel with derivatives 

21 """ 

22 def __init__(self, prior=None, kernel=None): 

23 if kernel is None: 

24 self.kernel = SquaredExponential() 

25 else: 

26 self.kernel = kernel 

27 

28 if prior is None: 

29 self.prior = ZeroPrior() 

30 else: 

31 self.prior = prior 

32 

33 def set_hyperparams(self, params): 

34 """Set hyperparameters of the regression. 

35 This is a list containing the parameters of the 

36 kernel and the regularization (noise) 

37 of the method as the last entry. 

38 """ 

39 self.hyperparams = params 

40 self.kernel.set_params(params[:-1]) 

41 self.noise = params[-1] 

42 

43 def train(self, X, Y, noise=None): 

44 """Produces a PES model from data. 

45 

46 Given a set of observations, X, Y, compute the K matrix 

47 of the Kernel given the data (and its cholesky factorization) 

48 This method should be executed whenever more data is added. 

49 

50 Parameters: 

51 

52 X: observations (i.e. positions). numpy array with shape: nsamples x D 

53 Y: targets (i.e. energy and forces). numpy array with 

54 shape (nsamples, D+1) 

55 noise: Noise parameter in the case it needs to be restated. 

56 """ 

57 if noise is not None: 

58 self.noise = noise # Set noise attribute to a different value 

59 

60 self.X = X.copy() # Store the data in an attribute 

61 n = self.X.shape[0] 

62 D = self.X.shape[1] 

63 regularization = np.array(n * ([self.noise * self.kernel.l] + 

64 D * [self.noise])) 

65 

66 K = self.kernel.kernel_matrix(X) # Compute the kernel matrix 

67 K[range(K.shape[0]), range(K.shape[0])] += regularization**2 

68 

69 self.m = self.prior.prior(X) 

70 self.a = Y.flatten() - self.m 

71 self.L, self.lower = cho_factor(K, lower=True, check_finite=True) 

72 cho_solve((self.L, self.lower), self.a, overwrite_b=True, 

73 check_finite=True) 

74 

75 def predict(self, x, get_variance=False): 

76 """Given a trained Gaussian Process, it predicts the value and the 

77 uncertainty at point x. It returns f and V: 

78 f : prediction: [y, grady] 

79 V : Covariance matrix. Its diagonal is the variance of each component 

80 of f. 

81 

82 Parameters: 

83 

84 x (1D np.array): The position at which the prediction is computed 

85 get_variance (bool): if False, only the prediction f is returned 

86 if True, the prediction f and the variance V are 

87 returned: Note V is O(D*nsample2) 

88 """ 

89 n = self.X.shape[0] 

90 k = self.kernel.kernel_vector(x, self.X, n) 

91 f = self.prior.prior(x) + np.dot(k, self.a) 

92 if get_variance: 

93 v = solve_triangular(self.L, k.T.copy(), lower=True, 

94 check_finite=False) 

95 variance = self.kernel.kernel(x, x) 

96 # covariance = np.matmul(v.T, v) 

97 covariance = np.tensordot(v, v, axes=(0, 0)) 

98 V = variance - covariance 

99 return f, V 

100 return f 

101 

102 def neg_log_likelihood(self, params, *args): 

103 """Negative logarithm of the marginal likelihood and its derivative. 

104 It has been built in the form that suits the best its optimization, 

105 with the scipy minimize module, to find the optimal hyperparameters. 

106 

107 Parameters: 

108 

109 l: The scale for which we compute the marginal likelihood 

110 *args: Should be a tuple containing the inputs and targets 

111 in the training set- 

112 """ 

113 X, Y = args 

114 # Come back to this 

115 self.kernel.set_params(np.array([params[0], params[1], self.noise])) 

116 self.train(X, Y) 

117 y = Y.flatten() 

118 

119 # Compute log likelihood 

120 logP = (-0.5 * np.dot(y - self.m, self.a) - 

121 np.sum(np.log(np.diag(self.L))) - 

122 X.shape[0] * 0.5 * np.log(2 * np.pi)) 

123 

124 # Gradient of the loglikelihood 

125 grad = self.kernel.gradient(X) 

126 

127 # vectorizing the derivative of the log likelihood 

128 D_P_input = np.array([np.dot(np.outer(self.a, self.a), g) 

129 for g in grad]) 

130 D_complexity = np.array([cho_solve((self.L, self.lower), g) 

131 for g in grad]) 

132 

133 DlogP = 0.5 * np.trace(D_P_input - D_complexity, axis1=1, axis2=2) 

134 return -logP, -DlogP 

135 

136 def fit_hyperparameters(self, X, Y, tol=1e-2, eps=None): 

137 """Given a set of observations, X, Y; optimize the scale 

138 of the Gaussian Process maximizing the marginal log-likelihood. 

139 This method calls TRAIN there is no need to call the TRAIN method 

140 again. The method also sets the parameters of the Kernel to their 

141 optimal value at the end of execution 

142 

143 Parameters: 

144 

145 X: observations(i.e. positions). numpy array with shape: nsamples x D 

146 Y: targets (i.e. energy and forces). 

147 numpy array with shape (nsamples, D+1) 

148 tol: tolerance on the maximum component of the gradient of the 

149 log-likelihood. 

150 (See scipy's L-BFGS-B documentation: 

151 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) 

152 eps: include bounds to the hyperparameters as a +- a percentage 

153 if eps is None there are no bounds in the optimization 

154 

155 Returns: 

156 

157 result (dict) : 

158 result = {'hyperparameters': (numpy.array) New hyperparameters, 

159 'converged': (bool) True if it converged, 

160 False otherwise 

161 } 

162 """ 

163 params = np.copy(self.hyperparams)[:2] 

164 arguments = (X, Y) 

165 if eps is not None: 

166 bounds = [((1 - eps) * p, (1 + eps) * p) for p in params] 

167 else: 

168 bounds = None 

169 

170 result = minimize(self.neg_log_likelihood, params, args=arguments, 

171 method='L-BFGS-B', jac=True, bounds=bounds, 

172 options={'gtol': tol, 'ftol': 0.01 * tol}) 

173 

174 if not result.success: 

175 converged = False 

176 else: 

177 converged = True 

178 self.hyperparams = np.array([result.x.copy()[0], 

179 result.x.copy()[1], self.noise]) 

180 self.set_hyperparams(self.hyperparams) 

181 return {'hyperparameters': self.hyperparams, 'converged': converged}