diff --git a/flare/gp.py b/flare/gp.py index 6b3e3b560..dcefa8aa9 100644 --- a/flare/gp.py +++ b/flare/gp.py @@ -8,7 +8,9 @@ from flare.gp_algebra import get_ky_mat, get_ky_and_hyp, \ get_like_from_ky_mat, get_like_grad_from_mats, get_neg_likelihood, \ get_neg_like_grad, get_ky_and_hyp_par - +from flare.kernels import str_to_kernel +import json as json +from flare.util import NumpyEncoder class GaussianProcess: """ Gaussian Process Regression Model. @@ -185,12 +187,6 @@ def predict(self, x_t: AtomicEnvironment, d: int) -> [float, float]: pred_var = self_kern - \ np.matmul(np.matmul(k_v, self.ky_mat_inv), k_v) - # # get predictive variance (possibly slow) - # v_vec = solve_triangular(self.l_mat, k_v, lower=True) - # self_kern = self.kernel(x_t, x_t, self.bodies, d, d, self.hyps, - # self.cutoffs) - # pred_var = self_kern - np.matmul(v_vec, v_vec) - return pred_mean, pred_var def predict_local_energy(self, x_t: AtomicEnvironment) -> float: @@ -292,8 +288,8 @@ def set_L_alpha(self): self.ky_mat_inv = ky_mat_inv self.l_mat_inv = l_mat_inv - self.like = like - self.like_grad = like_grad + self.likelihood = like + self.likelihood_gradient = like_grad def update_L_alpha(self): """ @@ -331,3 +327,78 @@ def update_L_alpha(self): self.alpha = alpha self.ky_mat_inv = ky_mat_inv self.l_mat_inv = l_mat_inv + + def __str__(self): + + thestr = "GaussianProcess Object\n" + thestr += 'Kernel: {}\n'.format(self.kernel_name) + thestr += "Training points: {}\n".format(len(self.training_data)) + thestr += 'Cutoffs: {}\n'.format(self.cutoffs) + thestr += 'Model Likelihood: {}\n'.format(self.likelihood) + + thestr += 'Hyperparameters: \n' + if self.hyp_labels is None: + # Put unlabeled hyperparameters on one line + thestr = thestr[:-1] + thestr += str(self.hyps)+'\n' + else: + for hyp, label in zip(self.hyps, self.hyp_labels): + thestr += "{}: {}\n".format(label, hyp) + + return thestr + + def as_dict(self): + + out_dict = dict(vars(self)) + + out_dict['training_data'] = [env.as_dict() for env in + self.training_data] + # Remove the callables + del out_dict['kernel'] + del out_dict['kernel_grad'] + + return out_dict + + @staticmethod + def from_dict(dictionary): + + force_kernel, grad = str_to_kernel(dictionary['kernel_name'], + include_grad=True) + + if dictionary['energy_kernel'] is not None: + energy_kernel = str_to_kernel(dictionary['energy_kernel']) + else: + energy_kernel = None + + if dictionary['energy_force_kernel'] is not None: + energy_force_kernel = str_to_kernel(dictionary[ + 'energy_force_kernel']) + else: + energy_force_kernel = None + + new_gp = GaussianProcess(kernel=force_kernel, + kernel_grad=grad, + energy_kernel=energy_kernel, + energy_force_kernel=energy_force_kernel, + cutoffs=np.array(dictionary['cutoffs']), + hyps=np.array(dictionary['hyps']), + hyp_labels=dictionary['hyp_labels'], + par=dictionary['par'], + maxiter=dictionary['maxiter'], + opt_algorithm=dictionary['algo']) + + # Save time by attempting to load in computed attributes + new_gp.l_mat = np.array(dictionary.get('l_mat', None)) + new_gp.l_mat_inv = np.array(dictionary.get('l_mat_inv', None)) + new_gp.alpha = np.array(dictionary.get('alpha', None)) + new_gp.ky_mat = np.array(dictionary.get('ky_mat', None)) + new_gp.ky_mat_inv = np.array(dictionary.get('ky_mat_inv', None)) + + new_gp.training_data = [AtomicEnvironment.from_dict(env) for env in + dictionary['training_data']] + new_gp.training_labels = dictionary['training_labels'] + + new_gp.likelihood = dictionary['likelihood'] + new_gp.likelihood_gradient = dictionary['likelihood_gradient'] + + return new_gp \ No newline at end of file diff --git a/tests/test_gp.py b/tests/test_gp.py index edd40f5f5..e091e35d3 100644 --- a/tests/test_gp.py +++ b/tests/test_gp.py @@ -35,7 +35,7 @@ def get_random_structure(cell, unique_species, noa): # set the scope to module so it will only be setup once @pytest.fixture(scope='module') -def two_body_gp(): +def two_body_gp()->GaussianProcess: """Returns a GP instance with a two-body numba-based kernel""" print("\nSetting up...\n") @@ -55,6 +55,7 @@ def two_body_gp(): GaussianProcess(kernel=en.three_body, kernel_grad=en.three_body_grad, hyps=np.array([1, 1, 1]), + hyp_labels=['Length', 'Signal Var.', 'Noise Var.'], cutoffs=cutoffs) gaussian.update_db(test_structure, forces) @@ -79,7 +80,7 @@ def params(): @pytest.fixture(scope='module') -def test_point(): +def test_point()->AtomicEnvironment: """Create test point for kernel to compare against""" # params cell = np.eye(3) @@ -212,3 +213,50 @@ def test_update_L_alpha(): ky_mat_from_set = np.copy(gp_model.ky_mat) assert (np.all(np.absolute(ky_mat_from_update-ky_mat_from_set)) < 1e-6) + +def test_representation_method(two_body_gp): + + the_str = str(two_body_gp) + assert 'GaussianProcess Object' in the_str + assert 'Kernel: three_body' in the_str + assert 'Cutoffs: [0.8 0.8]' in the_str + assert 'Model Likelihood: ' in the_str + assert 'Length: 1' in the_str + assert 'Signal Var.: 1' in the_str + assert "Noise Var.: 1" in the_str + + +def test_serialization_method(two_body_gp,test_point): + """ + Serialize and then un-serialize a GP and ensure that no info was lost. + Compare one calculation to ensure predictions work correctly. + :param two_body_gp: + :return: + """ + old_gp_dict = two_body_gp.as_dict() + new_gp = GaussianProcess.from_dict(old_gp_dict) + new_gp_dict = new_gp.as_dict() + + assert len(new_gp_dict) == len(old_gp_dict) + + for k1, k2 in zip(sorted(new_gp_dict.keys()), sorted(old_gp_dict.keys())): + + x = new_gp_dict[k1] + y = new_gp_dict[k2] + + if isinstance(x, np.ndarray): + assert np.equal(x, y).all() + elif hasattr(x, '__len__'): + + if isinstance(x[0], np.ndarray): + assert np.equal(x, y).all() + + else: + for xx, yy in zip(x, y): + assert xx == yy + else: + assert x == y + + for d in [0, 1, 2]: + assert np.all(two_body_gp.predict(x_t=test_point, d=d) == + new_gp.predict(x_t=test_point, d=d)) \ No newline at end of file