Added unit-tested string and dictionary methods for GPs

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

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))