gext/fitting.py

@@ -61,14 +61,12 @@ class DiffFitting(AbstractFitting):
     def fit(self, vectors: List[np.ndarray], target: np.ndarray):
         """Given a set of vectors and a target return the fitting
         coefficients."""
+        if len(vectors) == 1:
+            raise ValueError("DiffFit does not work for one vector")
         target=target-vectors[-1]
         VECTORS=[]
-        print("lenvector", len(vectors))
-        if len(vectors)>1:
-            for i in range(2, len(vectors)+1):
-                print("lenvector", len(vectors))
-                VECTORS.append(vectors[i-2]-vectors[-1])
-                print(len(VECTORS))
+        for i in range(0, len(vectors)+1):
+            VECTORS.append(vectors[i]-vectors[-1])
         matrix = np.array(VECTORS).T
         a = matrix.T @ matrix
         b = matrix.T @ target
@@ -82,15 +80,15 @@ class DiffFitting(AbstractFitting):
             coefficients: np. ndarray) -> np.ndarray:
         """Given a set of vectors (or matrices) and the corresponding
         coefficients, build their linear combination."""
+        if len(vectors) == 1:
+            raise ValueError("DiffFit does not work for one vector")
         result = np.zeros(vectors[0].shape, dtype=np.float64)
         VECTORS_DiffFitting=[]
-        if len(vectors)>1:
-            for i in range(2,len(vectors)+1):
-                  VECTORS_DiffFitting.append(vectors[i-2]-vectors[-1])
-            for coeff, vector in zip(coefficients, vectors):
+        for i in range(0,len(vectors)+1):
+              VECTORS_DiffFitting.append(vectors[i]-vectors[-1])
+        for coeff, vector in zip(coefficients, VECTORS_DiffFitting):
             result += vector*coeff
         result=result+vectors[-1]
-                print(result.shape)
         return result
 
 class LeastSquare(AbstractFitting):

test_diff_fitting.py

+import os
+import sys
+import numpy as np
+
+import gext
+import gext.descriptors
+import gext.fitting
+import gext.grassmann
+from tests import utils
+
+SMALL = 1e-8
+THRESHOLD = 5e-2
+
+regularization = 0.0
+
+# load test data from json file
+data = utils.load_json(f"tests/urea.json")
+nelectrons = data["nelectrons"]
+natoms = data["trajectory"].shape[1]
+nbasis = data["overlaps"].shape[1]
+nframes = data["trajectory"].shape[0]
+
+# initialize an extrapolator
+extrapolator = gext.Extrapolator(nelectrons, nbasis, natoms,
+    nsteps=nframes, fitting_regularization=regularization,
+    fitting="diff")
+
+# load data in the extrapolator
+for (coords, coeff, overlap) in zip(data["trajectory"],
+        data["coefficients"], data["overlaps"]):
+    extrapolator.load_data(coords, coeff, overlap)
+
+descriptors = extrapolator.descriptors.get(10)
+target = descriptors[-1]
+
+fitting_calculator = extrapolator.fitting_calculator
+
+# check if things are reasonable
+for start in range(0, 8):
+    vectors = descriptors[start:-1]
+    fit_coefficients = fitting_calculator.fit(vectors, target)
+    fitted_target = fitting_calculator.linear_combination(vectors, fit_coefficients)
+    error = np.linalg.norm(target - fitted_target, ord=np.inf)
+    assert error < THRESHOLD
+
+# if we put the target in the vectors used for the fitting,
+# check that we get an error smaller than the regularization
+vectors = descriptors[:-1]
+vectors[0] = target
+fit_coefficients = fitting_calculator.fit(vectors, target)
+fitted_target = fitting_calculator.linear_combination(vectors, fit_coefficients)
+
+print(np.linalg.norm(target - fitted_target, ord=np.inf), max(SMALL, regularization))
+assert np.linalg.norm(target - fitted_target, ord=np.inf) < max(SMALL, regularization)