Repository: uf3/uf3 (language: python)

Issue:
Fix to elusive bug in singlepoint training.
Single-point training (when the dataset only has 1 energy value) is often useful in active learning workflows involving UF3. This pull request has modified `least_squares.py` to allow single-point training. In addition, two new test cases have been added to `test_least_squares.py` related to this fix.

First, in `dataframe_to_tuples()`, I have removed the "not enough samples" exception. Having gone through the code, I think that this was put in place just to avoid the elusive bug that is solved by this pull request, which is described below.

If UF3 is trained using exactly one data point, the `fit()` and `fit_from_file()` methods raise a `RunTimeWarning` to warn about a divide-by-zero calculation, which arises because the standard deviation of a single point is zero. In the previous commit, there were already efforts to handle this divide-by-zero error in the `fit()` method but not in the `fit_from_file()` method:

    # from fit()
            try:
                energy_weight = 1 / len(y_e) / np.std(y_e)
                force_weight = 1 / len(y_f) / np.std(y_f)
            except (ZeroDivisionError, FloatingPointError):
                energy_weight = 1.0
                force_weight = 1 / len(y_f)

However, Numpy raises only a `RunTimeWarning`, not an error or an exception, so this issue is never handled. The result is that the fitted coefficients look something like:

    array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan,  0.,  0.,  0.])

This pull request fixes this by catching this `RunTimeWarning` in the try-except. By default, a warning is not treated as an exception unless a filter is added to elevate the warning to an exception status. That is the purpose of this line:

    warnings.filterwarnings("error", append=True) # to catch divide by zero warnings

After the try-except block, this filter is immediately removed so that it does not impact any other part of the code. This is the purpose of the line:

    warnings.filters.pop() # undo the filter



Merge singlepoint_training into master
Merge `singlepoint_training` branch to master.
Fix to elusive bug in singlepoint training.
Single-point training (when the dataset only has 1 energy value) is often useful in active learning workflows involving UF3. This pull request has modified `least_squares.py` to allow single-point training. In addition, two new test cases have been added to `test_least_squares.py` related to this fix.

First, in `dataframe_to_tuples()`, I have removed the "not enough samples" exception. Having gone through the code, I think that this was put in place just to avoid the elusive bug that is solved by this pull request, which is described below.

If UF3 is trained using exactly one data point, the `fit()` and `fit_from_file()` methods raise a `RunTimeWarning` to warn about a divide-by-zero calculation, which arises because the standard deviation of a single point is zero. In the previous commit, there were already efforts to handle this divide-by-zero error in the `fit()` method but not in the `fit_from_file()` method:

    # from fit()
            try:
                energy_weight = 1 / len(y_e) / np.std(y_e)
                force_weight = 1 / len(y_f) / np.std(y_f)
            except (ZeroDivisionError, FloatingPointError):
                energy_weight = 1.0
                force_weight = 1 / len(y_f)

However, Numpy raises only a `RunTimeWarning`, not an error or an exception, so this issue is never handled. The result is that the fitted coefficients look something like:

    array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan,  0.,  0.,  0.])

This pull request fixes this by catching this `RunTimeWarning` in the try-except. By default, a warning is not treated as an exception unless a filter is added to elevate the warning to an exception status. That is the purpose of this line:

    warnings.filterwarnings("error", append=True) # to catch divide by zero warnings

After the try-except block, this filter is immediately removed so that it does not impact any other part of the code. This is the purpose of the line:

    warnings.filters.pop() # undo the filter





Discussion / hints:

Fix patch:
diff --git a/uf3/regression/least_squares.py b/uf3/regression/least_squares.py
index 4c0ff9a8..916373cb 100644
--- a/uf3/regression/least_squares.py
+++ b/uf3/regression/least_squares.py
@@ -300,12 +300,10 @@ def fit(self,
                                   self.col_idx)
         gram_e, ord_e = batched_moore_penrose(x_e, y_e, batch_size=batch_size)
         if x_f is not None:
-            try:
-                energy_weight = 1 / np.sqrt(len(y_e)) / np.std(y_e)
-                force_weight = 1 / np.sqrt(len(y_f)) / np.std(y_f)
-            except (ZeroDivisionError, FloatingPointError):
-                energy_weight = 1.0
-                force_weight = 1 / np.sqrt(len(y_f))
+            energy_weight, force_weight = calc_E_F_weights(len(y_e),
+                                                           len(y_f),
+                                                           np.std(y_e),
+                                                           np.std(y_f))
             x_f, y_f = freeze_columns(x_f,
                                       y_f,
                                       self.mask,
@@ -403,8 +401,10 @@ def fit_from_file(self,
             gram_f += g_f
             ord_e += o_e
             ord_f += o_f
-        energy_weight = 1 / np.sqrt(e_variance.n) / e_variance.std
-        force_weight = 1 / np.sqrt(f_variance.n) / f_variance.std
+        energy_weight, force_weight = calc_E_F_weights(e_variance.n,
+                                                       f_variance.n,
+                                                       e_variance.std,
+                                                       f_variance.std)
         gram, ordinate = self.combine_weighted_gram(gram_e,
                                                     gram_f,
                                                     ord_e,
@@ -669,9 +669,6 @@ def dataframe_to_tuples(df_features,
         y (np.ndarray): target vector.
         w (np.ndarray): weight vector for machine learning.
     """
-    if len(df_features) <= 1:
-        raise ValueError(
-            "Not enough samples ({} provided)".format(len(df_features)))
     names = df_features.index.get_level_values(0)
     y_index = df_features.index.get_level_values(-1)
     energy_mask = (y_index == energy_key)
@@ -1122,3 +1119,28 @@ def find_pair_potential_well(coefficients, rounding_factor):
             # no actual well
             well_idx = peak_idx + 1
     return well_idx
+
+
+def calc_E_F_weights(n_e, n_f, std_e, std_f):
+    """
+    Calculates weights applied to energy and force components of the
+    least-squares problem (excluding kappa, which is applied in
+    self.combine_weighted_gram()).
+
+    Args:
+        n_e (int): number of energy samples.
+        n_f (int): number of force samples.
+        e_stddev (float): standard deviation of energy samples.
+        f_stddev (float): standard deviation of force samples.
+
+    Returns:
+        energy_weight (float): weight applied to energy components.
+        force_weight (float): weight applied to force components.
+    """
+    if std_e == 0:  # single point or really bad dataset
+        energy_weight = 1.0
+        force_weight = 1 / np.sqrt(n_f)
+    else:
+        energy_weight = 1 / np.sqrt(n_e) / std_e
+        force_weight = 1 / np.sqrt(n_f) / std_f
+    return energy_weight, force_weight
\ No newline at end of file


Test patch:
diff --git a/tests/data/singlepoint_fit/df_features_test_singlepoint_fit_from_file.h5 b/tests/data/singlepoint_fit/df_features_test_singlepoint_fit_from_file.h5
new file mode 100644
index 00000000..7a8653dc
Binary files /dev/null and b/tests/data/singlepoint_fit/df_features_test_singlepoint_fit_from_file.h5 differ
diff --git a/tests/test_least_squares.py b/tests/test_least_squares.py
index 3eab9052..30589b30 100644
--- a/tests/test_least_squares.py
+++ b/tests/test_least_squares.py
@@ -1,6 +1,10 @@
 import numpy as np
+import uf3
 from uf3.regression import least_squares
-from uf3.regression import regularize
+from uf3.data import composition
+from uf3.representation import bspline
+
+import os
 
 
 def simple_problem(n_features, n_samples, seed=0):
@@ -87,4 +91,81 @@ def test_frozen_coefficients():
 #     assert np.allclose([r[k] for k in args], [1, 1e-5, 0])
 #     args = ["c_2b", "c_3b", "r_1b", "r_2b", "r_3b"]  # d = 3
 #     r = least_squares.interpret_regularizer([3, -1], 3)
-#     assert np.allclose([r[k] for k in args], [1000, 0.1, 1e-5, 0, 0])
\ No newline at end of file
+#     assert np.allclose([r[k] for k in args], [1000, 0.1, 1e-5, 0, 0])
+
+def test_singlepoint_fit():
+    chemical_system = composition.ChemicalSystem(["Al"])
+    bspline_config = bspline.BSplineBasis(chemical_system)
+    n_features = sum(bspline_config.partition_sizes)
+    x_e, y_e, _ = simple_problem(n_features, 1, seed=0)  # single point
+    x_f, y_f, _ = simple_problem(n_features, 3, seed=1)
+    regularizer = np.eye(n_features) * 1e-6
+    model = least_squares.WeightedLinearModel(bspline_config,
+                                              regularizer=regularizer)
+    model.fit(x_e, y_e, x_f, y_f)
+    assert sum(~np.isfinite(model.coefficients)) == 0  # no nan or inf
+
+def test_singlepoint_fit_from_file():
+    chemical_system = composition.ChemicalSystem(["Al"])
+    bspline_config = bspline.BSplineBasis(chemical_system)
+    pkg_directory = os.path.dirname(os.path.dirname(uf3.__file__))
+    data_directory = os.path.join(pkg_directory, "tests/data")
+    features = os.path.join(data_directory, "singlepoint_fit",
+                            "df_features_test_singlepoint_fit_from_file.h5")
+    n_features = 19
+    regularizer = np.eye(n_features) * 1e-6
+    model = least_squares.WeightedLinearModel(bspline_config,
+                                              regularizer=regularizer)
+    model.fit_from_file(features,
+                        subset=['0_0'])
+    assert sum(~np.isfinite(model.coefficients)) == 0  # no nan or inf
+
+def test_loss_function():
+    chemical_system = composition.ChemicalSystem(["Al"], degree=2)
+    # Don't trim columns for this test
+    bspline_config = bspline.BSplineBasis(chemical_system,
+                                          leading_trim=0,
+                                          trailing_trim=0)
+    n_features = sum(bspline_config.partition_sizes)
+    n_samples_e = 30
+    n_samples_f = 500
+    x, y_true, _ = simple_problem(n_features, n_samples_e + n_samples_f, seed=0)
+    noise = np.random.normal(0, 0.1, n_samples_e + n_samples_f)
+    y = y_true + noise
+    x_e = x[:n_samples_e]
+    y_e = y[:n_samples_e]
+    x_f = x[n_samples_e:]
+    y_f = y[n_samples_e:]
+    regularizer = np.zeros((n_features, n_features))  # no regularizer for this test
+    model = least_squares.WeightedLinearModel(bspline_config,
+                                              regularizer=regularizer)
+    kappa = 0.25
+    model.fit(x_e, y_e, x_f, y_f, weight=kappa)
+
+    # Test if the fitted coefficients are the minimum to the loss function
+    e_weight = kappa / len(y_e) / np.var(y_e)
+    f_weight = (1 - kappa) / len(y_f) / np.var(y_f)
+    def loss_function(c):
+        e_loss = np.sum((y_e - np.dot(x_e, c)) ** 2) * e_weight
+        f_loss = np.sum((y_f - np.dot(x_f, c)) ** 2) * f_weight
+        return e_loss + f_loss
+
+    # loss evaluation for given coefficients
+    c_ref = model.coefficients
+    loss_ref = loss_function(c_ref)
+
+    # loss evaluation for perturbed coefficients
+    for i in range(len(c_ref)):
+        c = c_ref.copy()
+        c[i] += 1e-6
+        loss = loss_function(c)
+        assert loss > loss_ref
+        c[i] -= 2e-6
+        loss = loss_function(c)
+        assert loss > loss_ref
+
+    # loss evaluation for random perturbations
+    for i in range(10):
+        c = c_ref + np.random.normal(0, 1e-6, len(c_ref))
+        loss = loss_function(c)
+        assert loss > loss_ref