diff --git a/gext/descriptors.py b/gext/descriptors.py
index 838446f5afbb73dd3eecf25691882dd8c05f02ae..728e353e1b187ba0336a58aca3b62e2ca4653c06 100644
--- a/gext/descriptors.py
+++ b/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
diff --git a/gext/grassmann.py b/gext/grassmann.py
index f99d06cc5a56fb1bc57e3d1de14d0292f0a9a6cd..0fc84196b2b8aa2ecd698620fd33c7ab8bdaeed9 100644
--- a/gext/grassmann.py
+++ b/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:
diff --git a/gext/main.py b/gext/main.py
index 4a16d517ca12ad82c4fbebfc534e3defa7a0436e..8cf5333717f5313c7364fd6aa90f1063023ccc9b 100644
--- a/gext/main.py
+++ b/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,9 +90,17 @@ 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")
-
+
self.descriptor_calculator.set_options(**descriptor_options)
if self.options["fitting"] == "leastsquare":
@@ -115,7 +127,15 @@ class Extrapolator:
self._set_tangent(coeff)
# push the new data to the corresponding vectors
- self.descriptors.push(self._compute_descriptor(descriptor_input))
+ 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,27 +165,38 @@ class Extrapolator:
# use the descriptors to find the fitting coefficients
prev_descriptors= self.descriptors.get(n)
- descriptor = self._compute_descriptor(descriptor_input)
- fit_coefficients = self._fit(prev_descriptors, descriptor)
- print(fit_coefficients)
+ 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)
- # 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]))
-
- gamma = self.fitting_calculator.linear_combination(gammas, fit_coefficients)
- if self.options["verbose"]:
- fit_descriptor = self.fitting_calculator.linear_combination(
- prev_descriptors, fit_coefficients)
- print("error on descriptor:", \
- np.linalg.norm(fit_descriptor - descriptor, ord=np.inf))
-
- # if the overlap is not given, use the coefficients to fit
- # a new overlap
+ if self.options["method"]=="linear_interpolation":
+ fit_coefficients = self._fit(prev_descriptors, descriptor)
+
+ # use the fitting coefficients and the previous gammas to
+ # extrapolate a new gamma
+ gammas=[]
+ for i in range(len(coeffs)):
+ gammas.append(self._grassmann_log(coeffs[i]))
+
+ gamma = self.fitting_calculator.linear_combination(gammas, fit_coefficients)
+
+ if self.options["verbose"]:
+ fit_descriptor = self.fitting_calculator.linear_combination(
+ prev_descriptors, fit_coefficients)
+ 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:
overlaps = self.overlaps.get(n)
overlap = self.fitting_calculator.linear_combination(overlaps, fit_coefficients)
@@ -173,10 +204,11 @@ 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)
- return inverse_sqrt_overlap @ c_guess
-
+ 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."""
if self.options["tangent"] == "fixed":
@@ -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."""
- return self.descriptor_calculator.compute(descriptor_input)
+ 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