last modification of berny & geometric for executing fitting: diff &...

last modification of berny & geometric for executing fitting: diff & leastsquare, descriptor: distance & coulomb & H_core & overlap, tangent: fixed, one before the last

gext/descriptors.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 from scipy.spatial.distance import pdist
+from scipy.linalg import eigh
 
 class BaseDescriptor:
 
@@ -46,3 +47,17 @@ class FlattenMatrix(BaseDescriptor):
     def compute(self, matrix: np.ndarray) -> np.ndarray:
         """Compute the descriptor by flattening the matrix."""
         return matrix.flatten()
+
+class Coreeignvector(BaseDescriptor):
+
+
+    supported_options = {}
+
+    def __init__(self, **kwargs):
+        pass
+
+    def compute(self, main_matrix: np.ndarray, basis_matrix: np.ndarray) -> np.ndarray:
+        """Compute the descriptor by flattening the matrix."""
+        w,v=eigh(main_matrix, basis_matrix)
+
+        return w,v

gext/grassmann.py

@@ -17,6 +17,10 @@ def log(c: np.ndarray, c0: np.ndarray) -> np.ndarray:
     l = (np.identity(c.shape[0]) - c0 @ c0.T) @ cstar
     u, s, vt = np.linalg.svd(l, full_matrices=False)
     arcsin_s = np.diag(np.arcsin(s))
+    print('s', np.isnan(c))
+
+    print('s', np.isnan(c0))
+
     return u @ arcsin_s @ vt
 
 def exp(gamma: np.ndarray, c0: np.ndarray) -> np.ndarray:

gext/main.py

@@ -5,8 +5,9 @@ import numpy as np
 
 from . import grassmann
 from .fitting import LeastSquare, QuasiTimeReversible,DiffFitting
-from .descriptors import Distance, Coulomb, FlattenMatrix
+from .descriptors import Distance, Coulomb, FlattenMatrix, Coreeignvector
 from .buffer import CircularBuffer
+from scipy.linalg import eigh
 
 class Extrapolator:
 
@@ -22,6 +23,7 @@ class Extrapolator:
         "allow_partially_filled": True,
         "store_overlap": True,
         "tangent": "fixed",
+        "method": "lagrange_interpolation",
     }
 
     def __init__(self, nelectrons: int, nbasis: int, natoms: int, **kwargs):
@@ -50,7 +52,7 @@ class Extrapolator:
         self.options = {}
         descriptor_options = {}
         fitting_options = {}
-
+        method_options = {}
         # set specified options
         for key, value in kwargs.items():
             if key in self.supported_options:
@@ -59,6 +61,8 @@ class Extrapolator:
                 descriptor_options[key[11:]] = value
             elif key.startswith("fitting_"):
                 fitting_options[key[8:]] = value
+            #elif key.startswith("method_"):
+                #method_options[key[7:]] = value
             else:
                 raise ValueError(f"Unsupported option: {key}")
 
@@ -74,7 +78,7 @@ class Extrapolator:
 
         if isinstance(self.options["tangent"], int):
             if self.options["tangent"] < -self.options["nsteps"] or \
-                    self.options["tangent"] >= self.options["nsteps"]:
+                    self.options["tangent"] > self.options["nsteps"]:
                 raise ValueError("Unsupported tangent")
         else:
             if self.options["tangent"] != "fixed":
@@ -86,6 +90,14 @@ class Extrapolator:
             self.descriptor_calculator = Coulomb()
         elif self.options["descriptor"] == "H_core":
             self.descriptor_calculator = FlattenMatrix()
+        elif self.options["descriptor"] == "H_core_and_overlap":
+            self.descriptor_calculator = FlattenMatrix()
+        elif self.options["descriptor"]=="overlap":
+            self.descriptor_calculator = FlattenMatrix()
+        elif self.options["descriptor"]=="core_eign_vector" or "log_core_eign_vector":
+            self.descriptor_calculator=Coreeignvector()
+        elif self.options["descriptor"]=="index":
+            pass
         else:
             raise ValueError("Unsupported descriptor")
         
@@ -115,8 +127,16 @@ class Extrapolator:
             self._set_tangent(coeff)
 
         # push the new data to the corresponding vectors
+        if self.options["descriptor"]=="core_eign_vector":
+            self.descriptors.push(self._compute_descriptor(descriptor_input, overlap).flatten())
+        elif self.options["descriptor"]=="log_core_eign_vector":
+            self.descriptors.push(self._grassmann_log(self._normalize(self._compute_descriptor(descriptor_input, overlap), overlap)).flatten())
+        elif self.options["descriptor"]=="index":
+            self.descriptors.push(descriptor_input)
+        else: 
             self.descriptors.push(self._compute_descriptor(descriptor_input))
         
+
         if self.options["store_overlap"]:
             self.overlaps.push(overlap)
 
@@ -125,7 +145,7 @@ class Extrapolator:
         c_guess = self.guess_coefficients(coords, overlap)
         return c_guess @ c_guess.T
 
-    def guess_coefficients(self, descriptor_input: np.ndarray, overlap=None) -> np.ndarray:
+    def guess_coefficients(self, descriptor_input: np.ndarray, overlap) -> np.ndarray:
         """Get a new coefficient matrix to be used as a guess."""
 
         # check if we have enough data points to perform an extrapolation
@@ -145,13 +165,22 @@ class Extrapolator:
 
         # use the descriptors to find the fitting coefficients
         prev_descriptors= self.descriptors.get(n)
+        if self.options["descriptor"]=="core_eign_vector":
+            descriptor =self._compute_descriptor(descriptor_input, overlap).flatten()
+        elif self.options["descriptor"]=="log_core_eign_vector":
+            descriptor =self._grassmann_log(self._normalize(self._compute_descriptor(descriptor_input, overlap), overlap)).flatten()
+        elif self.options["descriptor"]=="index":
+            descriptor =descriptor_input
+        else:
             descriptor =self._compute_descriptor(descriptor_input)
+
+        coeffs=self.coeffs.get(n)
+
+        if self.options["method"]=="linear_interpolation":
             fit_coefficients = self._fit(prev_descriptors, descriptor)
-        print(fit_coefficients)
 
             # use the fitting coefficients and the previous gammas to
             # extrapolate a new gamma
-        coeffs=self.coeffs.get(n)
             gammas=[]
             for i in range(len(coeffs)):
                 gammas.append(self._grassmann_log(coeffs[i]))
@@ -164,6 +193,8 @@ class Extrapolator:
                 print("error on descriptor:", \
                     np.linalg.norm(fit_descriptor - descriptor, ord=np.inf))
 
+            # use the overlap and gamma to find a new set of coefficients
+            c_guess = self._grassmann_exp(gamma)
             # if the overlap is not given, use the coefficients to fit
             # a new overlap
         if overlap is None:
@@ -173,9 +204,10 @@ class Extrapolator:
         else:
             inverse_sqrt_overlap = self._inverse_sqrt_overlap(overlap)
 
-        # use the overlap and gamma to find a new set of coefficients
-        c_guess = self._grassmann_exp(gamma)
+        if self.options["method"]=="linear_interpolation":
             return inverse_sqrt_overlap @ c_guess
+        else:
+            return inverse_sqrt_overlap @ self._grassmann_exp(self._lagrange_polynomial(prev_descriptors, coeffs, descriptor))
     
     def _get_tangent(self) -> np.ndarray:
         """Get the tangent point."""
@@ -219,11 +251,39 @@ class Extrapolator:
         q, s, vt = np.linalg.svd(overlap, full_matrices=False)
         return q @ np.diag(1.0/np.sqrt(s)) @ vt
 
-    def _compute_descriptor(self, descriptor_input) -> np.ndarray:
+    def _compute_descriptor(self, descriptor_input, overlap=None) -> np.ndarray:
         """Given an input compute the corresponding descriptor."""
+        if overlap is not None:
+            #w,v= self.descriptor_calculator.compute(descriptor_input, overlap)
+            w,v = eigh(descriptor_input, overlap)
+            print('check', np.max(descriptor_input@ v - overlap@ v@ np.diag(w)))
+            print('check2', np.max(self._normalize(v,overlap).T@ overlap@ self._normalize(v, overlap)-np.eye(len(self._normalize(v,overlap).T@ overlap@ self._normalize(v, overlap)))))
+            return self._normalize(v, overlap)
+        else:
             return self.descriptor_calculator.compute(descriptor_input)
 
     def _fit(self, prev_descriptors, descriptor) -> np.ndarray:
         """Fit the current descriptor using previous descriptors and
         the specified fitting scheme."""
         return self.fitting_calculator.fit(prev_descriptors, descriptor)
+
+    def _lagrange_polynomial(self, prev_descriptors, coeffs, descriptor):
+        r = range(len(prev_descriptors))
+        print('r', r)
+        print('descriptor', descriptor)
+        print('prev_descriptors', prev_descriptors)
+        t=np.zeros((len(prev_descriptors)))
+        t[0]=0
+        for k in range(1,len(prev_descriptors)):
+            t[k] = np.linalg.norm(prev_descriptors[k]-prev_descriptors[k-1])+t[k-1]
+        a = [self._grassmann_log(coeffs[i]) / self._product(t[i] - t[j] for j in r if j != i) for i in r]    
+        return sum(a[i] * self._product([ np.linalg.norm(descriptor-prev_descriptors[k])+t[k] - t[j] for j in r if j != i]) for i in r)
+        if self.options['descriptor']=="index":
+            a = [self._grassmann_log(coeffs[i]) / self._product(prev_descriptors[i] - prev_descriptors[j] for j in r if j != i) for i in r]
+            return sum(a[i] * self._product([(descriptor - prev_descriptors[j]) for j in r if j != i]) for i in r)
+       
+    def _product(self, a):
+        p = 1.
+        for i in a:
+            p *= i
+        return p