Coverage for /builds/ase/ase/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 prior: Prior class, as in ase.optimize.gpmin.prior 

19 Defaults to ZeroPrior 

20 

21 kernel: Kernel function for the regression, as in 

22 ase.optimize.gpmin.kernel 

23 Defaults to the Squared Exponential kernel with derivatives 

24 """ 

25 

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

27 if kernel is None: 

28 self.kernel = SquaredExponential() 

29 else: 

30 self.kernel = kernel 

31 

32 if prior is None: 

33 self.prior = ZeroPrior() 

34 else: 

35 self.prior = prior 

36 

37 def set_hyperparams(self, params): 

38 """Set hyperparameters of the regression. 

39 This is a list containing the parameters of the 

40 kernel and the regularization (noise) 

41 of the method as the last entry. 

42 """ 

43 self.hyperparams = params 

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

45 self.noise = params[-1] 

46 

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

48 """Produces a PES model from data. 

49 

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

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

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

53 

54 Parameters: 

55 

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

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

58 shape (nsamples, D+1) 

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

60 """ 

61 if noise is not None: 

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

63 

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

65 n = self.X.shape[0] 

66 D = self.X.shape[1] 

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

68 D * [self.noise])) 

69 

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

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

72 

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

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

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

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

77 check_finite=True) 

78 

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

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

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

82 f : prediction: [y, grady] 

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

84 of f. 

85 

86 Parameters: 

87 

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

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

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

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

92 """ 

93 n = self.X.shape[0] 

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

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

96 if get_variance: 

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

98 check_finite=False) 

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

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

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

102 V = variance - covariance 

103 return f, V 

104 return f 

105 

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

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

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

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

110 

111 Parameters: 

112 

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

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

115 in the training set- 

116 """ 

117 X, Y = args 

118 # Come back to this 

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

120 self.train(X, Y) 

121 y = Y.flatten() 

122 

123 # Compute log likelihood 

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

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

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

127 

128 # Gradient of the loglikelihood 

129 grad = self.kernel.gradient(X) 

130 

131 # vectorizing the derivative of the log likelihood 

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

133 for g in grad]) 

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

135 for g in grad]) 

136 

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

138 return -logP, -DlogP 

139 

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

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

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

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

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

145 optimal value at the end of execution 

146 

147 Parameters: 

148 

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

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

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

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

153 log-likelihood. 

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

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

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

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

158 

159 Returns: 

160 

161 result (dict) : 

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

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

164 False otherwise 

165 } 

166 """ 

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

168 arguments = (X, Y) 

169 if eps is not None: 

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

171 else: 

172 bounds = None 

173 

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

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

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

177 

178 if not result.success: 

179 converged = False 

180 else: 

181 converged = True 

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

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

184 self.set_hyperparams(self.hyperparams) 

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