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

1import numpy as np 

2from scipy.linalg import cho_factor, cho_solve, solve_triangular 

3from scipy.optimize import minimize 

4 

5from ase.optimize.gpmin.kernel import SquaredExponential 

6from ase.optimize.gpmin.prior import ZeroPrior 

7 

8 

9class GaussianProcess(): 

10 """Gaussian Process Regression 

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

12 ase.optimize.gpmin 

13 

14 Parameters: 

15 

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

17 Defaults to ZeroPrior 

18 

19 kernel: Kernel function for the regression, as in 

20 ase.optimize.gpmin.kernel 

21 Defaults to the Squared Exponential kernel with derivatives 

22 """ 

23 

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

25 if kernel is None: 

26 self.kernel = SquaredExponential() 

27 else: 

28 self.kernel = kernel 

29 

30 if prior is None: 

31 self.prior = ZeroPrior() 

32 else: 

33 self.prior = prior 

34 

35 def set_hyperparams(self, params): 

36 """Set hyperparameters of the regression. 

37 This is a list containing the parameters of the 

38 kernel and the regularization (noise) 

39 of the method as the last entry. 

40 """ 

41 self.hyperparams = params 

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

43 self.noise = params[-1] 

44 

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

46 """Produces a PES model from data. 

47 

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

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

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

51 

52 Parameters: 

53 

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

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

56 shape (nsamples, D+1) 

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

58 """ 

59 if noise is not None: 

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

61 

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

63 n = self.X.shape[0] 

64 D = self.X.shape[1] 

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

66 D * [self.noise])) 

67 

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

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

70 

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

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

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

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

75 check_finite=True) 

76 

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

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

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

80 f : prediction: [y, grady] 

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

82 of f. 

83 

84 Parameters: 

85 

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

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

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

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

90 """ 

91 n = self.X.shape[0] 

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

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

94 if get_variance: 

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

96 check_finite=False) 

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

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

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

100 V = variance - covariance 

101 return f, V 

102 return f 

103 

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

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

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

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

108 

109 Parameters: 

110 

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

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

113 in the training set- 

114 """ 

115 X, Y = args 

116 # Come back to this 

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

118 self.train(X, Y) 

119 y = Y.flatten() 

120 

121 # Compute log likelihood 

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

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

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

125 

126 # Gradient of the loglikelihood 

127 grad = self.kernel.gradient(X) 

128 

129 # vectorizing the derivative of the log likelihood 

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

131 for g in grad]) 

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

133 for g in grad]) 

134 

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

136 return -logP, -DlogP 

137 

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

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

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

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

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

143 optimal value at the end of execution 

144 

145 Parameters: 

146 

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

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

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

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

151 log-likelihood. 

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

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

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

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

156 

157 Returns: 

158 

159 result (dict) : 

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

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

162 False otherwise 

163 } 

164 """ 

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

166 arguments = (X, Y) 

167 if eps is not None: 

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

169 else: 

170 bounds = None 

171 

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

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

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

175 

176 if not result.success: 

177 converged = False 

178 else: 

179 converged = True 

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

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

182 self.set_hyperparams(self.hyperparams) 

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