Commit 3de1629c authored by Malthe Kjær Bisbo's avatar Malthe Kjær Bisbo
Browse files

eta added in kernel

parent 8e764b0f
from __future__ import print_function
# ******NOTICE***************
# optimize.py module by Travis E. Oliphant
#
# You may copy and use this module as you see fit with no
# guarantee implied provided you keep this notice in all copies.
# *****END NOTICE************
import time
import numpy as np import numpy as np
from numpy import eye, absolute, sqrt, isinf
from ase.utils.linesearch import LineSearch
from ase.optimize.bfgslinesearch import BFGSLineSearch from ase.optimize.bfgslinesearch import BFGSLineSearch
from ase.constraints import FixAtoms from ase.constraints import FixAtoms
from ase.utils import basestring
from ase.data import covalent_radii from ase.data import covalent_radii
from ase.ga.utilities import get_mic_distance from ase.ga.utilities import get_mic_distance
from utils import check_valid_bondlengths from utils import check_valid_bondlengths
# These have been copied from Numeric's MLab.py
# I don't think they made the transition to scipy_core
# Modified from scipy_optimize
abs = absolute
pymin = min
pymax = max
__version__ = '0.1'
class BFGSLineSearch_constrained(BFGSLineSearch): class BFGSLineSearch_constrained(BFGSLineSearch):
def __init__(self, atoms, pos_init=None, restart=None, logfile='-', maxstep=.2, def __init__(self, atoms, pos_init=None, restart=None, logfile='-', maxstep=.2,
...@@ -118,3 +96,22 @@ class BFGSLineSearch_constrained(BFGSLineSearch): ...@@ -118,3 +96,22 @@ class BFGSLineSearch_constrained(BFGSLineSearch):
if not valid_bondlengths: if not valid_bondlengths:
valid_bondlengths = False valid_bondlengths = False
return valid_bondlengths return valid_bondlengths
def relax(structure, calc, Fmax=0.05, steps_max=200, dmax_cov=None):
a = structure.copy()
# Set calculator
a.set_calculator(calc)
pos_init = a.get_positions()
# Catch if linesearch fails
try:
dyn = BFGSLineSearch_constrained(a,
logfile=None,
pos_init=pos_init,
dmax_cov=dmax_cov)
dyn.run(fmax = Fmax, steps = steps_max)
except Exception as err:
print('Error in surrogate-relaxation:', err, flush=True)
traceback.print_exc()
traceback.print_exc(file=sys.stderr)
return a
\ No newline at end of file
...@@ -10,41 +10,19 @@ from ase.io import read, write, Trajectory ...@@ -10,41 +10,19 @@ from ase.io import read, write, Trajectory
from ase.calculators.singlepoint import SinglePointCalculator from ase.calculators.singlepoint import SinglePointCalculator
from ase.calculators.dftb import Dftb from ase.calculators.dftb import Dftb
from bfgslinesearch_constrained import BFGSLineSearch_constrained
from ase.ga.relax_attaches import VariansBreak
from parallel_utils import split, parallel_function_eval from parallel_utils import split, parallel_function_eval
from bfgslinesearch_constrained import relax
from candidate_operations.candidate_generation import CandidateGenerator from candidate_operations.candidate_generation import CandidateGenerator
from candidate_operations.basic_mutations import RattleMutation from candidate_operations.basic_mutations import RattleMutation
from mpi4py import MPI from mpi4py import MPI
world = MPI.COMM_WORLD world = MPI.COMM_WORLD
import traceback import traceback
import sys import sys
def relax(structure, calc, Fmax=0.05, steps_max=200, dmax_cov=None):
a = structure.copy()
# Set calculator
a.set_calculator(calc)
pos_init = a.get_positions()
# Catch if linesearch fails
try:
dyn = BFGSLineSearch_constrained(a,
logfile=None,
pos_init=pos_init,
dmax_cov=dmax_cov)
dyn.run(fmax = Fmax, steps = steps_max)
except Exception as err:
print('Error in surrogate-relaxation:', err, flush=True)
traceback.print_exc()
traceback.print_exc(file=sys.stderr)
return a
class GOFEE(): class GOFEE():
"""GOFEE global structure search method. """GOFEE global structure search method.
...@@ -254,9 +232,6 @@ class GOFEE(): ...@@ -254,9 +232,6 @@ class GOFEE():
a = rattle.get_new_candidate([a]) a = rattle.get_new_candidate([a])
a = self.evaluate(a) a = self.evaluate(a)
structures_init.append(a) structures_init.append(a)
# Potentially do a few relaxation steps.
### missing code ###
for a in structures_init: for a in structures_init:
self.write(a) self.write(a)
...@@ -350,9 +325,7 @@ class GOFEE(): ...@@ -350,9 +325,7 @@ class GOFEE():
relaxed_candidates = parallel_function_eval(self.comm, func2) relaxed_candidates = parallel_function_eval(self.comm, func2)
relaxed_candidates = self.certainty_filter(relaxed_candidates) relaxed_candidates = self.certainty_filter(relaxed_candidates)
#if (self.NsearchIter % self.Nuse_pop_as_candidates) == 0: relaxed_candidates = self.population.pop_MLrelaxed + relaxed_candidates
if self.steps % 3 == 0:
relaxed_candidates = self.population.pop_MLrelaxed + relaxed_candidates
if self.master: if self.master:
Epred = np.array([a.info['key_value_pairs']['Epred'] for a in relaxed_candidates]) Epred = np.array([a.info['key_value_pairs']['Epred'] for a in relaxed_candidates])
......
...@@ -80,7 +80,7 @@ class GPR(): ...@@ -80,7 +80,7 @@ class GPR():
if kernel is 'single': if kernel is 'single':
self.kernel = gauss_kernel() self.kernel = gauss_kernel()
elif kernel is 'double': elif kernel is 'double':
self.kernel = double_gauss_kernel() self.kernel = double_gauss_kernel(Nsplit_eta=self.descriptor.Nelements_2body)
else: else:
self.kernel = kernel self.kernel = kernel
......
...@@ -8,16 +8,16 @@ class kernel(ABC): ...@@ -8,16 +8,16 @@ class kernel(ABC):
self._theta = None self._theta = None
@abstractmethod @abstractmethod
def kernel(self): def __call__(self):
pass pass
# @abstractmethod @abstractmethod
# def kernel_vector(self): def kernel(self):
# pass pass
# @abstractmethod @abstractmethod
# def kernel_matrix(self): def kernel_vector(self):
# pass pass
@abstractmethod @abstractmethod
def kernel_jacobian(self): def kernel_jacobian(self):
...@@ -72,6 +72,7 @@ class kernel(ABC): ...@@ -72,6 +72,7 @@ class kernel(ABC):
self.theta = theta_down self.theta = theta_down
K_down = self(X, eval_gradient=False) K_down = self(X, eval_gradient=False)
dK_dTheta[i,:,:] = (K_up - K_down)/dx dK_dTheta[i,:,:] = (K_up - K_down)/dx
self.theta = theta
return dK_dTheta return dK_dTheta
...@@ -98,10 +99,12 @@ class gauss_kernel(kernel): ...@@ -98,10 +99,12 @@ class gauss_kernel(kernel):
else: else:
return K return K
def kernel(self, X,Y): def kernel(self, X,Y=None):
if np.ndim(X) == 1: if np.ndim(X) == 1:
X = X.reshape((1,-1)) X = X.reshape((1,-1))
if np.ndim(Y) == 1: if Y is None:
Y = X
elif np.ndim(Y) == 1:
Y = Y.reshape((1,-1)) Y = Y.reshape((1,-1))
d = cdist(X / self.length_scale, d = cdist(X / self.length_scale,
Y / self.length_scale, metric='sqeuclidean') Y / self.length_scale, metric='sqeuclidean')
...@@ -116,14 +119,6 @@ class gauss_kernel(kernel): ...@@ -116,14 +119,6 @@ class gauss_kernel(kernel):
K = self.kernel(x,Y).reshape(-1) K = self.kernel(x,Y).reshape(-1)
return K return K
def kernel_matrix(self, X,Y=None):
if Y is None:
d = cdist(X / self.length_scale, X / self.length_scale, metric='sqeuclidean')
else:
d = cdist(X / self.length_scale, Y / self.length_scale, metric='sqeuclidean')
K = self.amplitude * np.exp(-0.5 * d)
return K
def kernel_jacobian(self, X,Y): def kernel_jacobian(self, X,Y):
if np.ndim(X) == 1: if np.ndim(X) == 1:
X = X.reshape((1,-1)) X = X.reshape((1,-1))
...@@ -181,37 +176,41 @@ class double_gauss_kernel(kernel): ...@@ -181,37 +176,41 @@ class double_gauss_kernel(kernel):
length_scale1=10.0, length_scale1_bounds=(1e0, 1e3), length_scale1=10.0, length_scale1_bounds=(1e0, 1e3),
length_scale2=10.0, length_scale2_bounds=(1e0, 1e3), length_scale2=10.0, length_scale2_bounds=(1e0, 1e3),
weight=0.01, weight_bounds=(0.01,0.01), weight=0.01, weight_bounds=(0.01,0.01),
noise=1e-5, noise_bounds=(1e-5,1e-5)): noise=1e-5, noise_bounds=(1e-5,1e-5),
eta=1, eta_bounds=(0.1,10),
Nsplit_eta=None):
self.amplitude = amplitude self.amplitude = amplitude
self.length_scale1 = length_scale1 self.length_scale1 = length_scale1
self.length_scale2 = length_scale2 self.length_scale2 = length_scale2
self.weight = weight self.weight = weight
self.noise = noise self.noise = noise
self.eta = eta
self.Nsplit_eta = Nsplit_eta
self.amplitude_bounds = amplitude_bounds self.amplitude_bounds = amplitude_bounds
self.length_scale1_bounds = length_scale1_bounds self.length_scale1_bounds = length_scale1_bounds
self.length_scale2_bounds = length_scale2_bounds self.length_scale2_bounds = length_scale2_bounds
self.weight_bounds = weight_bounds self.weight_bounds = weight_bounds
self.noise_bounds = noise_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_scale1_bounds, length_scale2_bounds, weight_bounds, noise_bounds])) self.theta_bounds = np.log(np.array([amplitude_bounds, length_scale1_bounds, length_scale2_bounds, weight_bounds, noise_bounds, self.eta_bounds]))
def __call__(self, X, eval_gradient=False): def __call__(self, X, eval_gradient=False):
if np.ndim(X) == 1: K = self.kernel(X, with_noise=True)
X = X.reshape((1,-1))
d1 = cdist(X / self.length_scale1,
X / self.length_scale1, metric='sqeuclidean')
d2 = cdist(X / self.length_scale2,
X / self.length_scale2, metric='sqeuclidean')
K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2) + self.noise*np.eye(X.shape[0]))
if eval_gradient: if eval_gradient:
K_gradient = self.kernel_hyperparameter_gradient(X) K_gradient = self.kernel_hyperparameter_gradient(X)
return K, K_gradient return K, K_gradient
else: else:
return K return K
def kernel(self, X,Y=None): 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: if np.ndim(X) == 1:
X = X.reshape((1,-1)) X = X.reshape((1,-1))
if Y is None: if Y is None:
...@@ -222,7 +221,10 @@ class double_gauss_kernel(kernel): ...@@ -222,7 +221,10 @@ class double_gauss_kernel(kernel):
Y / self.length_scale1, metric='sqeuclidean') Y / self.length_scale1, metric='sqeuclidean')
d2 = cdist(X / self.length_scale2, d2 = cdist(X / self.length_scale2,
Y / self.length_scale2, metric='sqeuclidean') Y / self.length_scale2, metric='sqeuclidean')
K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2)) if with_noise:
K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2) + self.noise*np.eye(X.shape[0]))
else:
K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2))
return K return K
def kernel_value(self, x,y): def kernel_value(self, x,y):
...@@ -233,17 +235,11 @@ class double_gauss_kernel(kernel): ...@@ -233,17 +235,11 @@ class double_gauss_kernel(kernel):
K = self.kernel(x,Y).reshape(-1) K = self.kernel(x,Y).reshape(-1)
return K return K
def kernel_matrix(self, X,Y=None):
if Y is None:
d = cdist(X / self.length_scale, X / self.length_scale, metric='sqeuclidean')
else:
d = cdist(X / self.length_scale, Y / self.length_scale, metric='sqeuclidean')
K = np.exp(-0.5 * d)
return K
def kernel_jacobian(self, X,Y): def kernel_jacobian(self, X,Y):
""" Jacobian of the kernel with respect to X """ Jacobian of the kernel with respect to X
""" """
X = self.apply_eta(X)
Y = self.apply_eta(Y)
if np.ndim(X) == 1: if np.ndim(X) == 1:
X = X.reshape((1,-1)) X = X.reshape((1,-1))
if np.ndim(Y) == 1: if np.ndim(Y) == 1:
...@@ -264,7 +260,7 @@ class double_gauss_kernel(kernel): ...@@ -264,7 +260,7 @@ class double_gauss_kernel(kernel):
def theta(self): def theta(self):
"""Returns the log-transformed hyperparameters of the kernel. """Returns the log-transformed hyperparameters of the kernel.
""" """
self._theta = np.array([self.amplitude, self.length_scale1, self.length_scale2, self.weight, self.noise]) self._theta = np.array([self.amplitude, self.length_scale1, self.length_scale2, self.weight, self.noise, self.eta])
return np.log(self._theta) return np.log(self._theta)
@theta.setter @theta.setter
...@@ -279,6 +275,16 @@ class double_gauss_kernel(kernel): ...@@ -279,6 +275,16 @@ class double_gauss_kernel(kernel):
self.length_scale2 = self._theta[2] self.length_scale2 = self._theta[2]
self.weight = self._theta[3] self.weight = self._theta[3]
self.noise = self._theta[4] self.noise = self._theta[4]
self.eta = self._theta[5]
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): def dK_da(self, X):
d1 = cdist(X / self.length_scale1, d1 = cdist(X / self.length_scale1,
...@@ -310,8 +316,28 @@ class double_gauss_kernel(kernel): ...@@ -310,8 +316,28 @@ class double_gauss_kernel(kernel):
dK_dn = self.amplitude * self.noise * np.eye(X.shape[0]) dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
return dK_dn 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): def kernel_hyperparameter_gradient(self, X):
"""Calculates the derivative of the kernel with respect to the """Calculates the derivative of the kernel with respect to the
log transformed hyperparameters. log transformed hyperparameters.
""" """
return np.array([self.dK_da(X), self.dK_dl1(X), self.dK_dl2(X), self.dK_dw(X), self.dK_dn(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)])
No preview for this file type
import numpy as np import numpy as np
import unittest import unittest
from gpr import GPR as gpr from gpr import GPR as gpr
from bfgslinesearch_zlim import BFGSLineSearch_zlim from bfgslinesearch_constrained import BFGSLineSearch_constrained
from ase.io import read from ase.io import read
# old gpr
from kernels import RBF, ConstantKernel as C, WhiteKernel
from descriptor.fingerprint import Fingerprint
from delta_functions_multi.delta import delta as deltaFunc
from GPR import GPR
from custom_calculators import krr_calculator
def relax(structure, calc, Fmax=0.05, steps=200, dmax_cov=None, label='test'):
a = structure.copy()
# Set calculator
a.set_calculator(calc)
pos_init = a.get_positions()
# Catch if linesearch fails
try:
dyn = BFGSLineSearch_zlim(a,
logfile=label+'.log',
trajectory=label+'.traj',
pos_init=pos_init,
dmax_cov=dmax_cov)
dyn.run(fmax = Fmax, steps = steps)
except Exception as err:
print('Error in surrogate-relaxation:', err)
traceback.print_exc()
traceback.print_exc(file=sys.stderr)
return a
def initialize_old_gpr(atoms):
### Set up feature ###
# Radial part
Rc1 = 6
binwidth1 = 0.2
sigma1 = 0.2
# Angular part
Rc2 = 4
Nbins2 = 30
sigma2 = 0.2
gamma = 2
# Radial/angular weighting
eta = 20
use_angular = True
# Initialize feature
featureCalculator = Fingerprint(Rc1=Rc1, Rc2=Rc2, binwidth1=binwidth1, Nbins2=Nbins2, sigma1=sigma1, sigma2=sigma2, gamma=gamma, eta=eta, use_angular=use_angular)
#kernel = C(10, (1e1, 1e6)) * RBF(10, (1,1000))
kernel = C(100, (1e1, 1e5)) * (RBF(10, (1,1000)) + C(0.01, (0.01, 0.01)) * RBF(10, (1,1000)) + WhiteKernel(1e-5, (1e-5,1e-5)))
delta = deltaFunc(atoms=atoms, rcut=6)
gpr = GPR(kernel=kernel,
featureCalculator=featureCalculator,
delta_function=delta,
bias_func=None,
optimize=False,
n_restarts_optimizer=1)
return gpr
def get_E_with_std(traj, gpr):
E = []
F = []
for a in traj:
e = gpr.predict_energy(a, )
class test_gpr(unittest.TestCase): class test_gpr(unittest.TestCase):
@classmethod @classmethod
...@@ -85,67 +19,16 @@ class test_gpr(unittest.TestCase): ...@@ -85,67 +19,16 @@ class test_gpr(unittest.TestCase):
print('setUp') print('setUp')
#self.kernel = gauss_kernel() #self.kernel = gauss_kernel()
a = read('structures.traj', index='0') a = read('structures.traj', index='0')
self.gpr_old = initialize_old_gpr(a)
#self.gpr = gpr() #self.gpr = gpr()
self.gpr = gpr(kernel='double') self.gpr = gpr(kernel='double')
def tearDown(self): def tearDown(self):
print('tearDown\n') print('tearDown\n')
def test_compare_training_with_old(self):
traj = read('structures.traj', index=':50')
traj_train = traj[:40]
traj_predict = traj[40:]
self.gpr_old.train(traj_train, optimize=False)
self.gpr.train(traj_train)
np.testing.assert_almost_equal(self.gpr.alpha, self.gpr_old.alpha)
E_old = np.array([self.gpr_old.predict_energy(a, return_error=True)[:2] for a in traj_predict])
E = np.array([self.gpr.predict_energy(a, eval_std=True) for a in traj_predict])
np.testing.assert_almost_equal(E, E_old)
F_old = np.array([self.gpr_old.predict_force(a).reshape((-1,3)) for a in traj_predict])
F = np.array([self.gpr.predict_forces(a) for a in traj_predict])
np.testing.assert_almost_equal(F, F_old)
# The sign of the std are different in new and old GPR.
Fstd_old = -np.array([self.gpr_old.predict_force(a, return_error=True)[1].reshape((-1,3)) for a in traj_predict])
Fstd = np.array([self.gpr.predict_forces(a, eval_std=True)[1] for a in traj_predict])
np.testing.assert_almost_equal(Fstd, Fstd_old)
def test_compare_lml_with_old(self):
traj = read('structures.traj', index=':50')
traj_train = traj[:40]
traj_predict = traj[40:]
self.gpr_old.train(traj_train)
self.gpr.train(traj_train)
lml_old = self.gpr_old.log_marginal_likelihood_value_
lml_new = -self.gpr.neg_log_marginal_likelihood(eval_gradient=False)
np.testing.assert_almost_equal(lml_new, lml_old)
def test_compare_lml_gradient_with_old(self):
traj = read('structures.traj', index=':50')
traj_train = traj[:40]
traj_predict = traj[40:]
self.gpr_old.train(traj_train)
self.gpr.train(traj_train)
_, lml_ddTheta_old = self.gpr_old.log_marginal_likelihood(self.gpr_old.kernel.theta, eval_gradient=True)
_, lml_ddTheta = self.gpr.neg_log_marginal_likelihood(eval_gradient=True)
np.testing.assert_almost_equal(-lml_ddTheta, lml_ddTheta_old)
def test_lml_gradient(self): def test_lml_gradient(self):
traj = read('structures.traj', index=':50') traj = read('structures.traj', index=':50')
traj_train = traj[:40] traj_train = traj[:40]
traj_predict = traj[40:] traj_predict = traj[40:]