Coverage for ase / optimize / gpmin / gp.py: 53.85%

1# fmt: off 

2 

3import numpy as np 

4from scipy.linalg import cho_factor, cho_solve, solve_triangular 

5from scipy.optimize import minimize 

6 

7from ase.optimize.gpmin.kernel import SquaredExponential 

8from ase.optimize.gpmin.prior import ZeroPrior 

9 

10 

11class GaussianProcess(): 

12 """Gaussian Process Regression 

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

14 ase.optimize.gpmin 

15 

16 Parameters 

17 ---------- 

18 

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

20 Defaults to ZeroPrior 

21 

22 kernel: Kernel function for the regression, as in 

23 ase.optimize.gpmin.kernel 

24 Defaults to the Squared Exponential kernel with derivatives 

25 """ 

26 

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

28 if kernel is None: 

29 self.kernel = SquaredExponential() 

30 else: 

31 self.kernel = kernel 

32 

33 if prior is None: 

34 self.prior = ZeroPrior() 

35 else: 

36 self.prior = prior 

37 

38 def set_hyperparams(self, params): 

39 """Set hyperparameters of the regression. 

40 This is a list containing the parameters of the 

41 kernel and the regularization (noise) 

42 of the method as the last entry. 

43 """ 

44 self.hyperparams = params 

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

46 self.noise = params[-1] 

47 

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

49 """Produces a PES model from data. 

50 

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

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

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

54 

55 Parameters 

56 ---------- 

57 

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

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

60 shape (nsamples, D+1) 

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

62 """ 

63 if noise is not None: 

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

65 

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

67 n = self.X.shape[0] 

68 D = self.X.shape[1] 

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

70 D * [self.noise])) 

71 

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

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

74 

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

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

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

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

79 check_finite=True) 

80 

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

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

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

84 f : prediction: [y, grady] 

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

86 of f. 

87 

88 Parameters 

89 ---------- 

90 

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

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

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

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

95 """ 

96 n = self.X.shape[0] 

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

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

99 if get_variance: 

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

101 check_finite=False) 

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

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

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

105 V = variance - covariance 

106 return f, V 

107 return f 

108 

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

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

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

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

113 

114 Parameters 

115 ---------- 

116 

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

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

119 in the training set- 

120 """ 

121 X, Y = args 

122 # Come back to this 

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

124 self.train(X, Y) 

125 y = Y.flatten() 

126 

127 # Compute log likelihood 

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

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

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

131 

132 # Gradient of the loglikelihood 

133 grad = self.kernel.gradient(X) 

134 

135 # vectorizing the derivative of the log likelihood 

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

137 for g in grad]) 

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

139 for g in grad]) 

140 

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

142 return -logP, -DlogP 

143 

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

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

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

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

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

149 optimal value at the end of execution 

150 

151 Parameters 

152 ---------- 

153 

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

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

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

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

158 log-likelihood. 

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

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

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

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

163 

164 Returns 

165 ------- 

166 

167 result (dict) : 

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

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

170 False otherwise 

171 } 

172 """ 

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

174 arguments = (X, Y) 

175 if eps is not None: 

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

177 else: 

178 bounds = None 

179 

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

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

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

183 

184 if not result.success: 

185 converged = False 

186 else: 

187 converged = True 

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

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

190 self.set_hyperparams(self.hyperparams) 

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