Skip to content
Snippets Groups Projects
Commit 2d5d5431 authored by Michele Nottoli's avatar Michele Nottoli
Browse files

Regression.

parent f08e50f5
No related branches found
No related tags found
No related merge requests found
Pipeline #2107 failed
...@@ -31,3 +31,14 @@ class Coulomb(Distance): ...@@ -31,3 +31,14 @@ class Coulomb(Distance):
def compute(self, coords: np.ndarray) -> np.ndarray: def compute(self, coords: np.ndarray) -> np.ndarray:
"""Compute the Coulomb matrix as a descriptor.""" """Compute the Coulomb matrix as a descriptor."""
return 1.0/super().compute(coords) return 1.0/super().compute(coords)
class Flatten():
def __init__(self, **kwargs):
pass
def set_options(self, **kwargs):
pass
def compute(self, matrix: np.ndarray) -> np.ndarray:
return matrix.flatten()
...@@ -28,27 +28,35 @@ class AbstractFitting(abc.ABC): ...@@ -28,27 +28,35 @@ class AbstractFitting(abc.ABC):
self.options[option] = default_value self.options[option] = default_value
@abc.abstractmethod @abc.abstractmethod
def fit(self, vectors: List[np.ndarray], target:np.ndarray): def train(self, descriptors: List[np.ndarray], gammas: List[np.ndarray]):
"""Base method for computing new fitting coefficients.""" """Base method for computing new fitting coefficients."""
return np.zeros(0) return np.zeros(0)
def linear_combination(self, vectors: List[np.ndarray], @abc.abstractmethod
coefficients: np. ndarray) -> np.ndarray: def extrapolate(self, descriptor: np.ndarray):
"""Given a set of vectors (or matrices) and the corresponding """Extrapolate to a new descriptor:"""
coefficients, build their linear combination.""" return np.zeros(0)
result = np.zeros(vectors[0].shape, dtype=np.float64)
for coeff, vector in zip(coefficients, vectors):
result += vector*coeff
return result
class DiffFitting(AbstractFitting): class DiffFitting(AbstractFitting):
"""least square minimization fitting of the diffvectors""" """Least square minimization fitting with a reference."""
supported_options = { supported_options = {
"regularization": 0.0, "normalize": False,
"regularization": 1e-4,
"ref": -1
} }
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.gammas = []
self.ref = np.zeros(0)
self.A = np.zeros(0)
self.DT = np.zeros(0)
self.trained = False
def set_options(self, **kwargs): def set_options(self, **kwargs):
"""Set options for least square minimization""" """Set options for least square minimization"""
super().set_options(**kwargs) super().set_options(**kwargs)
...@@ -57,41 +65,66 @@ class DiffFitting(AbstractFitting): ...@@ -57,41 +65,66 @@ class DiffFitting(AbstractFitting):
or self.options["regularization"] > 100: or self.options["regularization"] > 100:
raise ValueError("Unsupported value for regularization") raise ValueError("Unsupported value for regularization")
def fit(self, vectors: List[np.ndarray], target: np.ndarray): def train(self, descriptors: List[np.ndarray], gammas: List[np.ndarray]):
"""Given a set of vectors and a target return the fitting """Given a set of descriptors and a target return the fitting
coefficients.""" coefficients."""
if len(vectors) == 1: if len(descriptors) == 1:
raise ValueError("DiffFitting does not work for one vector") raise ValueError("DiffFitting does not work for one vector")
ref = len(vectors) - 1 ref = len(descriptors) - 1
target = target - vectors[ref]
diff_vectors = [] self.ref = descriptors[ref]
for i in range(0, len(vectors)):
if i != ref:
diff_vectors.append(vectors[i] - vectors[ref])
matrix = np.array(diff_vectors).T diff_gammas = []
a = matrix.T @ matrix diff_descriptors = []
b = matrix.T @ target for i in range(0, len(descriptors)):
if i != ref:
diff_descriptors.append(descriptors[i] - self.ref)
diff_gammas.append(gammas[i]) # we assume the ref is tangent
self.gammas = diff_gammas
D = np.array(diff_descriptors).T
A = D.T @ D
norm = 1.0
if self.options["normalize"]:
norm = np.linalg.norm(D)
if self.options["regularization"] > 0.0: if self.options["regularization"] > 0.0:
a += np.identity(len(b))*self.options["regularization"]**2 A += np.identity(len(A))*self.options["regularization"]**2*norm**2
coefficients = np.linalg.solve(a, b)
self.A = A
# Convert diff coefficients to normal coefficients (like in self.DT = D.T
# least square fitting) self.trained = True
coefficients = np.insert(coefficients, ref, 1.0 - np.sum(coefficients))
def extrapolate(self, target: np.ndarray, vectors=None):
return coefficients if not self.trained:
raise ValueError("The model is not trained")
b = self.DT @ (target - self.ref)
coefficients = np.linalg.solve(self.A, b)
if vectors is None:
vectors = self.gammas
result = np.zeros(vectors[0].shape)
for c, vector in zip(coefficients, vectors):
result += vector * c
return result
class LeastSquare(AbstractFitting): class LeastSquare(AbstractFitting):
"""Simple least square minimization fitting.""" """Simple least square minimization fitting."""
supported_options = { supported_options = {
"regularization": 1e-3, "normalize": False,
"regularization": 1e-4,
} }
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.gammas = []
self.A = np.zeros(0)
self.DT = np.zeros(0)
self.trained = False
def set_options(self, **kwargs): def set_options(self, **kwargs):
"""Set options for least square minimization""" """Set options for least square minimization"""
super().set_options(**kwargs) super().set_options(**kwargs)
...@@ -100,25 +133,54 @@ class LeastSquare(AbstractFitting): ...@@ -100,25 +133,54 @@ class LeastSquare(AbstractFitting):
or self.options["regularization"] > 100: or self.options["regularization"] > 100:
raise ValueError("Unsupported value for regularization") raise ValueError("Unsupported value for regularization")
def fit(self, vectors: List[np.ndarray], target: np.ndarray): def train(self, descriptors: List[np.ndarray], gammas: List[np.ndarray]):
"""Given a set of vectors and a target return the fitting self.gammas = gammas
coefficients."""
matrix = np.array(vectors).T D = np.array(descriptors).T
a = matrix.T @ matrix A = D.T @ D
b = matrix.T @ target norm = 1.0
if self.options["normalize"]:
norm = np.linalg.norm(D)
if self.options["regularization"] > 0.0: if self.options["regularization"] > 0.0:
a += np.identity(len(b))*self.options["regularization"]**2 A += np.identity(len(A))*self.options["regularization"]**2*norm**2
coefficients = np.linalg.solve(a, b)
return np.array(coefficients, dtype=np.float64) self.A = A
self.DT = D.T
self.trained = True
def extrapolate(self, target: np.ndarray, vectors=None):
if not self.trained:
raise ValueError("The model is not trained")
b = self.DT @ target
coefficients = np.linalg.solve(self.A, b)
if vectors is None:
vectors = self.gammas
result = np.zeros(vectors[0].shape)
for c, vector in zip(coefficients, vectors):
result += vector * c
return result
class QuasiTimeReversible(AbstractFitting): class QuasiTimeReversible(AbstractFitting):
"""Quasi time reversible fitting scheme.""" """Quasi time reversible fitting scheme."""
supported_options = { supported_options = {
"regularization": 1e-3, "normalize": False,
"regularization": 1e-4,
"full_second_order": False
} }
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.gammas = []
self.past_target = np.zeros(0)
self.A = np.zeros(0)
self.DT = np.zeros(0)
self.q = 0
self.trained = False
def set_options(self, **kwargs): def set_options(self, **kwargs):
"""Set options for quasi time reversible fitting""" """Set options for quasi time reversible fitting"""
super().set_options(**kwargs) super().set_options(**kwargs)
...@@ -127,44 +189,65 @@ class QuasiTimeReversible(AbstractFitting): ...@@ -127,44 +189,65 @@ class QuasiTimeReversible(AbstractFitting):
or self.options["regularization"] > 100: or self.options["regularization"] > 100:
raise ValueError("Unsupported value for regularization") raise ValueError("Unsupported value for regularization")
def fit(self, vectors: List[np.ndarray], target: np.ndarray): def train(self, descriptors: List[np.ndarray], gammas: List[np.ndarray]):
"""Given a set of vectors and a target return the fitting """Given a set of descriptors and a target return the fitting
coefficients in a quasi time reversible scheme.""" coefficients in a quasi time reversible scheme."""
past_target = vectors[0] self.past_target = descriptors[0]
matrix = np.array(vectors[1:]).T self.gammas = gammas
q = matrix.shape[1] D = np.array(descriptors[1:]).T
if q == 1:
time_reversible_matrix = matrix self.q = D.shape[1]
elif q%2 == 0: if self.q == 1:
time_reversible_matrix = matrix[:, :q//2] + matrix[:, :q//2-1:-1] time_reversible_D = D
elif self.q%2 == 0:
time_reversible_D = D[:, :self.q//2] + D[:, :self.q//2-1:-1]
else: else:
time_reversible_matrix = matrix[:, :q//2+1] + matrix[:, :q//2-1:-1] time_reversible_D = D[:, :self.q//2+1] + D[:, :self.q//2-1:-1]
norm = 1.0
if self.options["normalize"]:
norm = np.linalg.norm(time_reversible_D)
a = time_reversible_matrix.T @ time_reversible_matrix A = time_reversible_D.T @ time_reversible_D
b = time_reversible_matrix.T @ (target + past_target)
if self.options["regularization"] > 0.0: if self.options["regularization"] > 0.0:
a += np.identity(len(b))*self.options["regularization"]**2 A += np.identity(len(A))*self.options["regularization"]**2*norm**2
coefficients = np.linalg.solve(a, b)
self.A = A
self.DT = time_reversible_D.T
self.trained = True
def extrapolate(self, target: np.ndarray, vectors=None):
if q == 1: b = self.DT @ (target + self.past_target)
coefficients = np.linalg.solve(self.A, b)
if self.q == 1:
full_coefficients = np.concatenate(([-1.0], coefficients)) full_coefficients = np.concatenate(([-1.0], coefficients))
elif q%2 == 0: elif self.q%2 == 0:
full_coefficients = np.concatenate(([-1.0], coefficients, full_coefficients = np.concatenate(([-1.0], coefficients,
coefficients[::-1])) coefficients[::-1]))
else: else:
full_coefficients = np.concatenate(([-1.0], coefficients[:-1], full_coefficients = np.concatenate(([-1.0], coefficients[:-1],
2.0*coefficients[-1:], coefficients[-2::-1])) 2.0*coefficients[-1:], coefficients[-2::-1]))
return np.array(full_coefficients, dtype=np.float64)
if vectors is None:
vectors = self.gammas
result = np.zeros(vectors[0].shape)
for vector, c in zip(vectors, full_coefficients):
result += vector * c
return result
class PolynomialRegression(AbstractFitting): class PolynomialRegression(AbstractFitting):
"""Polynomial regression.""" """Polynomial regression."""
supported_options = { supported_options = {
"regularization": 1e-3, "normalize": False,
"regularization": 1e-4,
"ref": -1, "ref": -1,
"order": 1} "order": 1}
...@@ -178,8 +261,8 @@ class PolynomialRegression(AbstractFitting): ...@@ -178,8 +261,8 @@ class PolynomialRegression(AbstractFitting):
self.matrix = np.zeros(0, dtype=np.float64) self.matrix = np.zeros(0, dtype=np.float64)
self.gamma_shape = None self.gamma_shape = None
def fit(self, descriptor_list: List[np.ndarray], gamma_list: List[np.ndarray]): def train(self, descriptor_list: List[np.ndarray], gamma_list: List[np.ndarray]):
"""Given a set of vectors and a set of gammas, construct the """Given a set of descriptors and a set of gammas, construct the
transformation matrix.""" transformation matrix."""
if self.gamma_shape is None: if self.gamma_shape is None:
...@@ -214,13 +297,24 @@ class PolynomialRegression(AbstractFitting): ...@@ -214,13 +297,24 @@ class PolynomialRegression(AbstractFitting):
gammas = np.reshape(gamma_list, gammas = np.reshape(gamma_list,
(len(gamma_list), self.gamma_shape[0]*self.gamma_shape[1])).T (len(gamma_list), self.gamma_shape[0]*self.gamma_shape[1])).T
norm = 1.0
if self.options["normalize"]:
norm = np.linalg.norm(d)
if self.options["order"] >= 2: if self.options["order"] >= 2:
if self.options["full_second_order"]:
outer = np.outer(d, d)
d = np.vstack((d, outer.flatten()))
else:
d = np.vstack((d, d**2)) d = np.vstack((d, d**2))
a = d @ np.linalg.inv(d.T @ d + np.identity(d.shape[1])*self.options["regularization"]**2) A = d.T @ d
self.matrix = a @ gammas.T if self.options["regularization"] > 0.0:
A += np.identity(len(A))*self.options["regularization"]**2*norm**2
self.matrix = d @ np.linalg.inv(A) @ gammas.T
def apply(self, descriptor): def extrapolate(self, descriptor):
"""Apply the matrix to the current descriptor.""" """Apply the matrix to the current descriptor."""
if self.options["ref"]: if self.options["ref"]:
descriptor -= self.descriptor_ref descriptor -= self.descriptor_ref
......
...@@ -5,7 +5,7 @@ import numpy as np ...@@ -5,7 +5,7 @@ import numpy as np
from . import grassmann from . import grassmann
from .fitting import LeastSquare, QuasiTimeReversible, PolynomialRegression, DiffFitting from .fitting import LeastSquare, QuasiTimeReversible, PolynomialRegression, DiffFitting
from .descriptors import Distance, Coulomb from .descriptors import Distance, Coulomb, Flatten
from .buffer import CircularBuffer from .buffer import CircularBuffer
class Extrapolator: class Extrapolator:
...@@ -78,6 +78,8 @@ class Extrapolator: ...@@ -78,6 +78,8 @@ class Extrapolator:
self.descriptor_calculator = Distance() self.descriptor_calculator = Distance()
elif self.options["descriptor"] == "coulomb": elif self.options["descriptor"] == "coulomb":
self.descriptor_calculator = Coulomb() self.descriptor_calculator = Coulomb()
elif self.options["descriptor"] == "flatten":
self.descriptor_calculator = Flatten()
else: else:
raise ValueError("Unsupported descriptor") raise ValueError("Unsupported descriptor")
self.descriptor_calculator.set_options(**descriptor_options) self.descriptor_calculator.set_options(**descriptor_options)
...@@ -108,6 +110,9 @@ class Extrapolator: ...@@ -108,6 +110,9 @@ class Extrapolator:
if self.options["store_overlap"]: if self.options["store_overlap"]:
self.overlaps.push(overlap) self.overlaps.push(overlap)
if self.tangent is None:
self._set_tangent(coeff)
def guess(self, coords: np.ndarray, overlap=None) -> np.ndarray: def guess(self, coords: np.ndarray, overlap=None) -> np.ndarray:
"""Get a new electronic density matrix to be used as a guess.""" """Get a new electronic density matrix to be used as a guess."""
c_guess = self.guess_coefficients(coords, overlap) c_guess = self.guess_coefficients(coords, overlap)
...@@ -139,20 +144,14 @@ class Extrapolator: ...@@ -139,20 +144,14 @@ class Extrapolator:
descriptor = self._compute_descriptor(coords) descriptor = self._compute_descriptor(coords)
# use the descriptors to find the fitting coefficients # use the descriptors to find the fitting coefficients
if self.options["fitting"] == "polynomialregression": self.fitting_calculator.train(prev_descriptors, gammas)
self.fitting_calculator.fit(prev_descriptors, gammas) gamma = self.fitting_calculator.extrapolate(descriptor)
gamma = self.fitting_calculator.apply(descriptor)
else:
fit_coefficients = self._fit(prev_descriptors, descriptor)
# use the fitting coefficients and the previous gammas to
# extrapolate a new gamma
gamma = self.fitting_calculator.linear_combination(gammas, fit_coefficients)
if self.options["verbose"]: #if self.options["fitting"] in ["leastsquare"]:
fit_descriptor = self.fitting_calculator.linear_combination( # fit_descriptor = self.fitting_calculator.extrapolate(descriptor,
prev_descriptors, fit_coefficients) # vectors=prev_descriptors)
print("error on descriptor:", \ # print("error on descriptor:", \
np.linalg.norm(fit_descriptor - descriptor, ord=np.inf)) # np.linalg.norm(fit_descriptor - descriptor, ord=np.inf))
# if the overlap is not given, use the coefficients to fit # if the overlap is not given, use the coefficients to fit
# a new overlap # a new overlap
...@@ -160,7 +159,7 @@ class Extrapolator: ...@@ -160,7 +159,7 @@ class Extrapolator:
if self.options["fitting"] == "polynomialregression": if self.options["fitting"] == "polynomialregression":
raise ValueError("The option polynomial regression needs the overlap") raise ValueError("The option polynomial regression needs the overlap")
overlaps = self.overlaps.get(n) overlaps = self.overlaps.get(n)
overlap = self.fitting_calculator.linear_combination(overlaps, fit_coefficients) overlap = self.fitting_calculator.extrapolate(descriptor, vectors=overlaps)
inverse_sqrt_overlap = self._inverse_sqrt_overlap(overlap) inverse_sqrt_overlap = self._inverse_sqrt_overlap(overlap)
else: else:
inverse_sqrt_overlap = self._inverse_sqrt_overlap(overlap) inverse_sqrt_overlap = self._inverse_sqrt_overlap(overlap)
...@@ -174,11 +173,11 @@ class Extrapolator: ...@@ -174,11 +173,11 @@ class Extrapolator:
count = self.descriptors.count count = self.descriptors.count
n = min(self.options["nsteps"], count) n = min(self.options["nsteps"], count)
if "ref" in self.fitting_calculator.options: if "ref" in self.fitting_calculator.options:
ref = self.fitting_calculator.options["ref"] index = self.fitting_calculator.options["ref"]
else:
ref = 0
coefficients = self.coefficients.get(n) coefficients = self.coefficients.get(n)
return coefficients[ref] return coefficients[index]
else:
return self.tangent
def _crop_coeff(self, coeff) -> np.ndarray: def _crop_coeff(self, coeff) -> np.ndarray:
"""Crop the coefficient matrix to remove the virtual orbitals.""" """Crop the coefficient matrix to remove the virtual orbitals."""
...@@ -218,8 +217,3 @@ class Extrapolator: ...@@ -218,8 +217,3 @@ class Extrapolator:
def _compute_descriptor(self, coords) -> np.ndarray: def _compute_descriptor(self, coords) -> np.ndarray:
"""Given a set of coordinates compute the corresponding descriptor.""" """Given a set of coordinates compute the corresponding descriptor."""
return self.descriptor_calculator.compute(coords) return self.descriptor_calculator.compute(coords)
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)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment