Commit 2e0569ec authored by Malthe Kjær Bisbo's avatar Malthe Kjær Bisbo
Browse files

eta included

parent b18a63a1
......@@ -23,6 +23,8 @@ world = MPI.COMM_WORLD
import traceback
import sys
from time import time
class GOFEE():
"""GOFEE global structure search method.
......@@ -124,11 +126,6 @@ class GOFEE():
print('Using calculator from supplied structure(s)')
self.calc = calc
if gpr is not None:
self.gpr = gpr
else:
self.gpr = GPR()
if startgenerator is None:
assert structures is not None
self.startgenerator = None
......@@ -153,9 +150,6 @@ class GOFEE():
rattle_range=4)
self.candidate_generator = CandidateGenerator([1.0],[rattle])
# Initialize population
self.population = population(population_size=population_size, gpr=self.gpr, similarity2equal=0.9999)
# Define parallel communication
self.comm = world.Dup() # Important to avoid mpi-problems from call to ase.parallel in BFGS
self.master = self.comm.rank == 0
......@@ -180,6 +174,14 @@ class GOFEE():
if restart is None or not isfile(restart):
self.initialize()
if gpr is not None:
self.gpr = gpr
else:
self.gpr = GPR(template_structure=self.structures[0])
# Initialize population
self.population = population(population_size=population_size, gpr=self.gpr, similarity2equal=0.9999)
else:
self.read()
self.comm.barrier()
......@@ -226,6 +228,9 @@ class GOFEE():
self.structures.append(a)
def evaluate_initial_structures(self):
""" Evaluate energies and forces of all initial structures
(self.structures) that have not yet been evaluated.
"""
structures_init = []
for a in self.structures:
calc = a.get_calculator()
......@@ -323,13 +328,11 @@ class GOFEE():
def func1():
return [self.generate_candidate() for i in task_split[self.comm.rank]]
candidates = parallel_function_eval(self.comm, func1)
def func2():
return [self.surrogate_relaxation(candidates[i], Fmax=0.1, steps=200, kappa=self.kappa)
for i in task_split[self.comm.rank]]
relaxed_candidates = parallel_function_eval(self.comm, func2)
relaxed_candidates = self.certainty_filter(relaxed_candidates)
relaxed_candidates = self.population.pop_MLrelaxed + relaxed_candidates
if self.master:
......@@ -352,7 +355,6 @@ class GOFEE():
"""
calc = self.gpr.get_calculator(kappa)
a_relaxed = relax(a, calc, dmax_cov=self.dmax_cov, Fmax=Fmax, steps_max=steps)
# Evaluate uncertainty
E, Estd = self.gpr.predict_energy(a_relaxed, eval_std=True)
......
This diff is collapsed.
......@@ -112,7 +112,7 @@ class Fingerprint(object):
self.initialized = False
def initialize_from_atoms(self, atoms):
self.nsigma = 4
self.nsigma = 6
self.volume = atoms.get_volume()
self.pbc = atoms.get_pbc()
......
......@@ -71,16 +71,22 @@ class GPR():
prior:
"""
def __init__(self, descriptor=None, kernel='double', prior=None, n_restarts_optimizer=1):
def __init__(self, descriptor=None, kernel='double', prior=None, n_restarts_optimizer=1, template_structure=None):
if descriptor is None:
self.descriptor = Fingerprint()
else:
self.descriptor = descriptor
Nsplit_eta = None
if template_structure is not None:
self.descriptor.initialize_from_atoms(template_structure)
if hasattr(self.descriptor, 'use_angular'):
if self.descriptor.use_angular:
Nsplit_eta = self.descriptor.Nelements_2body
if kernel is 'single':
self.kernel = gauss_kernel()
self.kernel = gauss_kernel(Nsplit_eta=Nsplit_eta)
elif kernel is 'double':
self.kernel = double_gauss_kernel()
self.kernel = double_gauss_kernel(Nsplit_eta=Nsplit_eta)
else:
self.kernel = kernel
......@@ -102,7 +108,6 @@ class GPR():
if eval_std:
# Lines 5 and 6 in GPML
vk = np.dot(self.K_inv, k)
#K0 = self.kernel.kernel(x,x)
E_std = np.sqrt(self.K0 - np.dot(k, vk))
return E, E_std
else:
......@@ -123,11 +128,9 @@ class GPR():
if eval_std:
k = self.kernel.kernel_vector(x, self.X)
vk = np.dot(self.K_inv, k)
#K0 = self.kernel.kernel_value(x,x)
g = self.K0 - np.dot(k.T, vk)
assert g >= 0
F_std = 1/np.sqrt(g) * np.dot(k_ddr.T, vk)
#F_std = -1/np.sqrt(g) * np.dot(k_ddr.T, vk)
return F.reshape((-1,3)), F_std.reshape(-1,3)
else:
return F.reshape(-1,3)
......@@ -153,6 +156,7 @@ class GPR():
def optimize_hyperparameters(self, atoms_list=None, comm=None):
if self.n_restarts_optimizer == 0:
self.train(atoms_list)
return
if atoms_list is not None:
self.memory.save_data(atoms_list)
......@@ -207,7 +211,6 @@ class GPR():
L = cholesky(K, lower=True)
alpha = cho_solve((L, True), self.Y)
#log_likelihood = -0.5 * np.einsum("i,i", Y, alpha)
lml = -0.5 * np.dot(self.Y, alpha)
lml -= np.sum(np.log(np.diag(L)))
lml -= K.shape[0]/2 * np.log(2*np.pi)
......
......@@ -3,6 +3,8 @@ from abc import ABC, abstractmethod
from scipy.spatial.distance import pdist, cdist, squareform
from time import time
class kernel(ABC):
def __init__(self):
self._theta = None
......@@ -40,7 +42,6 @@ class kernel(ABC):
y = y.reshape((1,-1))
nx = len(x)
ny = y.shape[0]
f0 = self.kernel(x,y)
f_jac = np.zeros((ny,nx))
for i in range(nx):
x_up = np.copy(x)
......@@ -50,6 +51,7 @@ class kernel(ABC):
f_up = self.kernel(x_up,y)
f_down = self.kernel(x_down,y)
f_jac[:,i] = (f_up - f_down)/dx
return f_jac
......@@ -71,13 +73,18 @@ class kernel(ABC):
K_up = self(X, eval_gradient=False)
self.theta = theta_down
K_down = self(X, eval_gradient=False)
dK_dTheta[i,:,:] = (K_up - K_down)/dx
self.theta = theta
return dK_dTheta
class gauss_kernel(kernel):
def __init__(self, amplitude=10.0, length_scale=10.0, amplitude_bounds=(1e0, 1e3), length_scale_bounds=(1e-1, 1e1)):
def __init__(self, amplitude=100.0, amplitude_bounds=(1e0, 1e5),
length_scale=10.0, length_scale_bounds=(1e0, 1e3),
noise=1e-5, noise_bounds=(1e-5,1e-5),
eta=1, eta_bounds=(0.1,10),
Nsplit_eta=None):
self.amplitude = amplitude
self.length_scale = length_scale
......@@ -86,6 +93,22 @@ class gauss_kernel(kernel):
self.theta_bounds = np.array([amplitude_bounds, length_scale_bounds])
self.amplitude = amplitude
self.length_scale = length_scale
self.noise = noise
self.eta = eta
self.Nsplit_eta = Nsplit_eta
self.amplitude_bounds = amplitude_bounds
self.length_scale_bounds = length_scale_bounds
self.noise_bounds = noise_bounds
if self.Nsplit_eta is None:
self.eta_bounds = (eta,eta)
else:
self.eta_bounds = eta_bounds
self.theta_bounds = np.log(np.array([amplitude_bounds, length_scale_bounds, noise_bounds, self.eta_bounds]))
def __call__(self, X, eval_gradient=False):
if np.ndim(X) == 1:
X = X.reshape((1,-1))
......@@ -98,17 +121,32 @@ class gauss_kernel(kernel):
return K, K_gradient
else:
return K
def kernel(self, X,Y=None):
def __call__(self, X, eval_gradient=False):
K = self.kernel(X, with_noise=True)
if eval_gradient:
K_gradient = self.kernel_hyperparameter_gradient(X)
return K, K_gradient
else:
return K
def kernel(self, X, Y=None, with_noise=False):
if with_noise:
assert Y is None
if np.ndim(X) == 1:
X = X.reshape((1,-1))
if Y is None:
Y = X
elif np.ndim(Y) == 1:
Y = Y.reshape((1,-1))
X = self.apply_eta(X)
Y = self.apply_eta(Y)
d = cdist(X / self.length_scale,
Y / self.length_scale, metric='sqeuclidean')
K = self.amplitude * np.exp(-0.5 * d)
if with_noise:
K = self.amplitude * (np.exp(-0.5 * d) + self.noise*np.eye(X.shape[0]))
else:
K = self.amplitude * np.exp(-0.5 * d)
return K
def kernel_value(self, x,y):
......@@ -120,22 +158,29 @@ class gauss_kernel(kernel):
return K
def kernel_jacobian(self, X,Y):
""" Jacobian of the kernel with respect to X
"""
X = self.apply_eta(X)
Y = self.apply_eta(Y)
if np.ndim(X) == 1:
X = X.reshape((1,-1))
if np.ndim(Y) == 1:
Y = Y.reshape((1,-1))
K = self.kernel(X,Y).T
dK_dd = -1./(2*self.length_scale**2)*K
d = cdist(X / self.length_scale,
Y / self.length_scale, metric='sqeuclidean')
dK_dd = -self.amplitude * 1/(2*self.length_scale**2) * np.exp(-0.5 * d)
X = self.apply_eta(X)
Y = self.apply_eta(Y)
dd_df = 2*(X - Y)
dk_df = np.multiply(dK_dd, dd_df)
dk_df = np.multiply(dK_dd.T, dd_df)
return dk_df
@property
def theta(self):
"""Returns the log-transformed hyperparameters of the kernel.
"""
self._theta = np.array([self.amplitude, self.length_scale])
self._theta = np.array([self.amplitude, self.length_scale, self.noise, self.eta])
return np.log(self._theta)
@theta.setter
......@@ -147,28 +192,61 @@ class gauss_kernel(kernel):
self._theta = np.exp(theta)
self.amplitude = self._theta[0]
self.length_scale = self._theta[1]
self.noise = self._theta[2]
self.eta = self._theta[3]
def apply_eta(self, X):
Xeta = np.copy(X)
if self.Nsplit_eta is not None:
if np.ndim(X) == 1:
Xeta[self.Nsplit_eta:] *= self.eta
else:
Xeta[:,self.Nsplit_eta:] *= self.eta
return Xeta
def dK_da(self, X):
if np.ndim(X) == 1:
X = X.reshape((1,-1))
d = cdist(X / self.length_scale,
X / self.length_scale, metric='sqeuclidean')
dK_da = self.amplitude * np.exp(-0.5 * d)
dK_da = self.amplitude * (np.exp(-0.5 * d) + self.noise*np.eye(X.shape[0]))
return dK_da
def dK_dl(self, X):
if np.ndim(X) == 1:
X = X.reshape((1,-1))
d = cdist(X / self.length_scale,
X / self.length_scale, metric='sqeuclidean')
dK_dl = self.amplitude * d * np.exp(-0.5 * d)
return dK_dl
def dK_dn(self, X):
dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
return dK_dn
def dK_deta(self, X, dx=1e-5):
N_data = X.shape[0]
theta = np.copy(self.theta)
dK_deta = np.zeros((N_data, N_data))
theta_up = np.copy(theta)
theta_down = np.copy(theta)
theta_up[-1] += 0.5*dx
theta_down[-1] -= 0.5*dx
self.theta = theta_up
K_up = self(X, eval_gradient=False)
self.theta = theta_down
K_down = self(X, eval_gradient=False)
dK_dTheta = (K_up - K_down)/dx
self.theta = theta
return dK_dTheta
def kernel_hyperparameter_gradient(self, X):
"""Calculates the derivative of the kernel with respect to the
log transformed hyperparameters.
"""
return np.array([self.dK_da(X), self.dK_dl(X)])
dK_deta = self.dK_deta(X)
X = self.apply_eta(X)
return np.array([self.dK_da(X), self.dK_dl(X), self.dK_dn(X), dK_deta])
class double_gauss_kernel(kernel):
......@@ -210,13 +288,14 @@ class double_gauss_kernel(kernel):
def kernel(self, X, Y=None, with_noise=False):
if with_noise:
assert Y is None
X = self.apply_eta(X)
if np.ndim(X) == 1:
X = X.reshape((1,-1))
if Y is None:
Y = X
elif np.ndim(Y) == 1:
Y = Y.reshape((1,-1))
X = self.apply_eta(X)
Y = self.apply_eta(Y)
d1 = cdist(X / self.length_scale1,
Y / self.length_scale1, metric='sqeuclidean')
d2 = cdist(X / self.length_scale2,
......@@ -251,6 +330,8 @@ class double_gauss_kernel(kernel):
dK1_dd1 = -1/(2*self.length_scale1**2) * np.exp(-0.5 * d1)
dK2_dd2 = -1/(2*self.length_scale2**2) * np.exp(-0.5 * d2)
dK_dd = self.amplitude * (dK1_dd1 + self.weight*dK2_dd2)
X = self.apply_eta(X)
Y = self.apply_eta(Y)
dd_df = 2*(X - Y)
dk_df = np.multiply(dK_dd.T, dd_df)
......@@ -339,5 +420,6 @@ class double_gauss_kernel(kernel):
"""Calculates the derivative of the kernel with respect to the
log transformed hyperparameters.
"""
dK_deta = self.dK_deta(X)
X = self.apply_eta(X)
return np.array([self.dK_da(X), self.dK_dl1(X), self.dK_dl2(X), self.dK_dw(X), self.dK_dn(X), self.dK_deta(X)])
return np.array([self.dK_da(X), self.dK_dl1(X), self.dK_dl2(X), self.dK_dw(X), self.dK_dn(X), dK_deta])
......@@ -20,7 +20,8 @@ class test_gpr(unittest.TestCase):
#self.kernel = gauss_kernel()
a = read('structures.traj', index='0')
#self.gpr = gpr()
self.gpr = gpr(kernel='double')
#self.gpr = gpr(kernel='double', template_structure=a)
self.gpr = gpr(kernel='single', template_structure=a)
def tearDown(self):
print('tearDown\n')
......@@ -33,7 +34,7 @@ class test_gpr(unittest.TestCase):
_, lml_ddTheta = self.gpr.neg_log_marginal_likelihood(eval_gradient=True)
lml_ddTheta_numeric = self.gpr.numerical_neg_lml()
np.testing.assert_almost_equal(lml_ddTheta, lml_ddTheta_numeric)
np.testing.assert_almost_equal(lml_ddTheta, lml_ddTheta_numeric, decimal=6)
def test_forces(self):
traj = read('structures.traj', index=':50')
......@@ -45,7 +46,7 @@ class test_gpr(unittest.TestCase):
a = traj_predict[0]
F = self.gpr.predict_forces(a)
F_numeric = self.gpr.numerical_forces(a)
np.testing.assert_almost_equal(F, F_numeric)
np.testing.assert_almost_equal(F, F_numeric, decimal=6)
def test_forces_std(self):
traj = read('structures.traj', index=':50')
......@@ -57,7 +58,7 @@ class test_gpr(unittest.TestCase):
a = traj_predict[0]
_, Fstd = self.gpr.predict_forces(a, eval_std=True)
_, Fstd_numeric = self.gpr.numerical_forces(a, eval_std=True)
np.testing.assert_almost_equal(Fstd, Fstd_numeric)
np.testing.assert_almost_equal(Fstd, Fstd_numeric, decimal=6)
if __name__ == '__main__':
unittest.main()
......
......@@ -17,7 +17,8 @@ class test_kernel(unittest.TestCase):
#self.kernel = gauss_kernel()
#self.kernel_old = C(10, (1e1, 1e6)) * RBF(10, (1,1000))
self.kernel = double_gauss_kernel(Nsplit_eta=2)
#self.kernel = double_gauss_kernel(Nsplit_eta=2, eta=10)
self.kernel = gauss_kernel(Nsplit_eta=2, eta=10)
def tearDown(self):
print('tearDown\n')
......@@ -41,9 +42,9 @@ class test_kernel(unittest.TestCase):
X = np.c_[x1,x2].T
K_ddtheta = self.kernel.kernel_hyperparameter_gradient(X)
K_ddtheta_num = self.kernel.numerical_hyperparameter_gradient(X)
print(K_ddtheta, '\n\n')
print(K_ddtheta_num, '\n\n')
print(K_ddtheta-K_ddtheta_num, '\n')
#print(K_ddtheta, '\n\n')
#print(K_ddtheta_num, '\n\n')
#print(K_ddtheta-K_ddtheta_num, '\n')
np.testing.assert_almost_equal(K_ddtheta, K_ddtheta_num)
if __name__ == '__main__':
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment