From d8cb1ce25c29f7e56baaf52c4b7de47cf6cf3bd6 Mon Sep 17 00:00:00 2001
From: Daniel Varjas <dvarjas@gmail.com>
Date: Thu, 21 Sep 2023 16:33:06 +0200
Subject: [PATCH] fix error by eliminating small terms in model

---
 qsymm/model.py           | 18 ++++++++++++++++++
 qsymm/symmetry_finder.py |  9 +++++++--
 2 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/qsymm/model.py b/qsymm/model.py
index 2ec3b23..63b1458 100644
--- a/qsymm/model.py
+++ b/qsymm/model.py
@@ -705,6 +705,24 @@ class Model(UserDict):
             return all(allclose(self[key], other[key], rtol, atol, equal_nan)
                        for key in self.keys() | other.keys())
 
+    def eliminate_zeros(self, rtol=1e-05, atol=1e-08):
+        """Return a model with small terms removed. Tolerances are
+        in Frobenius matrix norm, relative tolerance compares to the
+        value with largest norm."""
+        if not issubclass(self.format, (np.ndarray, scipy.sparse.spmatrix)):
+            raise ValueError('Operation only supported for Models with '
+                             '`np.ndarray` or `scipy.sparse.spmatrix` data type.')
+        result = self.zeros_like()
+        # For empty model do nothing
+        if self.data == {}:
+            return result
+        # Write it explicitely so it works with sparse arrays
+        norm = lambda mat: np.sqrt(np.sum(np.abs(mat)**2))
+        max_norm = np.max([norm(val) for val in self.values()])
+        tol = max(atol, max_norm * rtol)
+        result.data = {key: copy(val) for key, val in self.items() if not norm(val) < tol}
+        return result
+
 
 class BlochModel(Model):
     def __init__(self, hamiltonian=None, locals=None, momenta=('k_x', 'k_y', 'k_z'),
diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py
index 6056201..8e3a0e2 100644
--- a/qsymm/symmetry_finder.py
+++ b/qsymm/symmetry_finder.py
@@ -488,9 +488,14 @@ def _find_unitary(model, Ps, g, sparse=False, checks=False):
                          '`np.ndarray` or `scipy.sparse.spmatrix` data type.')
     if g.U is not None:
         raise ValueError('g.U must be None.')
-    Rmodel = g.apply(model)
-    if set(model) != set(Rmodel):
+
+    # Remove potential small terms
+    model = model.eliminate_zeros()
+    Rmodel = g.apply(model).eliminate_zeros()
+    # After eliminating small entries, if g is a symmetry only the same keys should appear
+    if model.keys() != Rmodel.keys():
         return g
+
     HR, HL = [], []
     # Only test eigenvalues if all matrices are Hermitian
     ev_test = True
-- 
GitLab