Replace hard-coded np.allclose/isclose and math.isclose (for comp…

…lex expression) (#4164)

* replace hard-coded math.isclose

* add todo tag

* replace more isclose in code

* remove TODO tag

* replace more in code

* fix seemingly wrong quote position

* avoid import when it's used only once or very few

* replace last batch

* revert change to isotropic check

* replace some numpy isclose allclose

* remove debug tag

* remove some hard coded np allclose

* revert some change on very simple evals

* avoid unnecessary compare with zero

* revert simple comparison

* avoid minus zero

* revert some simple expressions

* simplify sci notation

* revert simple comparisons

* avoid 1.0e-x as it's already float

* revert simple compare

* use sci not

* revert simple

* use abs as we don't need always float

* fix round usage

* all close

* revert as i'm not sure about the shape broadcasting

* avoid import from numpy

* clean up math import, reduce namespace cluster

* simplify all close

* sci notation

* simplify import of math

* simplify assert all close

* simplify int(len(a) / b) to len(a) // b

dev_scripts/chemenv/explicit_permutations_plane_algorithm.py

@@ -140,7 +140,7 @@
         for icsm, csm in enumerate(csms):
             found = False
             for csm2 in csms_with_recorded_permutation:
-                if np.isclose(csm, csm2, rtol=0.0, atol=1.0e-6):
+                if np.isclose(csm, csm2, rtol=0.0, atol=1e-6):
                     found = True
                     break
             if not found:

src/pymatgen/analysis/bond_valence.py

@@ -3,10 +3,10 @@
 from __future__ import annotations
 
 import functools
+import math
 import operator
 import os
 from collections import defaultdict
-from math import exp, sqrt
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -55,8 +55,8 @@ def calculate_bv_sum(site, nn_list, scale_factor=1.0):
             r2 = BV_PARAMS[el2]["r"]
             c1 = BV_PARAMS[el1]["c"]
             c2 = BV_PARAMS[el2]["c"]
-            R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
-            vij = exp((R - nn.nn_distance * scale_factor) / 0.31)
+            R = r1 + r2 - r1 * r2 * (math.sqrt(c1) - math.sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
+            vij = math.exp((R - nn.nn_distance * scale_factor) / 0.31)
             bv_sum += vij * (1 if el1.X < el2.X else -1)
     return bv_sum
 
@@ -91,8 +91,8 @@ def calculate_bv_sum_unordered(site, nn_list, scale_factor=1):
                     r2 = BV_PARAMS[el2]["r"]
                     c1 = BV_PARAMS[el1]["c"]
                     c2 = BV_PARAMS[el2]["c"]
-                    R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
-                    vij = exp((R - nn.nn_distance * scale_factor) / 0.31)
+                    R = r1 + r2 - r1 * r2 * (math.sqrt(c1) - math.sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
+                    vij = math.exp((R - nn.nn_distance * scale_factor) / 0.31)
                     bv_sum += occu1 * occu2 * vij * (1 if el1.X < el2.X else -1)
     return bv_sum
 
@@ -173,7 +173,7 @@ def _calc_site_probabilities(self, site, nn):
                 sigma = data["std"]
                 # Calculate posterior probability. Note that constant
                 # factors are ignored. They have no effect on the results.
-                prob[sp.oxi_state] = exp(-((bv_sum - u) ** 2) / 2 / (sigma**2)) / sigma * PRIOR_PROB[sp]
+                prob[sp.oxi_state] = math.exp(-((bv_sum - u) ** 2) / 2 / (sigma**2)) / sigma * PRIOR_PROB[sp]
         # Normalize the probabilities
         try:
             prob = {k: v / sum(prob.values()) for k, v in prob.items()}
@@ -194,7 +194,7 @@ def _calc_site_probabilities_unordered(self, site, nn):
                     sigma = data["std"]
                     # Calculate posterior probability. Note that constant
                     # factors are ignored. They have no effect on the results.
-                    prob[el][sp.oxi_state] = exp(-((bv_sum - u) ** 2) / 2 / (sigma**2)) / sigma * PRIOR_PROB[sp]
+                    prob[el][sp.oxi_state] = math.exp(-((bv_sum - u) ** 2) / 2 / (sigma**2)) / sigma * PRIOR_PROB[sp]
             # Normalize the probabilities
             try:
                 prob[el] = {k: v / sum(prob[el].values()) for k, v in prob[el].items()}
@@ -263,7 +263,7 @@ def get_valences(self, structure: Structure):
                     # Retain probabilities that are at least 1/100 of highest prob.
                     filtered = list(
                         filter(
-                            lambda v: prob[elem.symbol][v] > 0.001 * prob[elem.symbol][val[0]],
+                            lambda v: prob[elem.symbol][v] > 1e-3 * prob[elem.symbol][val[0]],
                             val,
                         )
                     )

src/pymatgen/analysis/chemenv/coordination_environments/voronoi.py

@@ -173,7 +173,7 @@ def setup_voronoi_list(self, indices, voronoi_cutoff):
 
                     min_dist = min([min_dist, distances[ridge_point2]])
                     for iii, sss in enumerate(self.structure):
-                        if neighbors[ridge_point2].is_periodic_image(sss, tolerance=1.0e-6):
+                        if neighbors[ridge_point2].is_periodic_image(sss, tolerance=1e-6):
                             idx = iii
                             break
                     results2.append(

src/pymatgen/analysis/chemenv/utils/math_utils.py

@@ -2,9 +2,9 @@
 
 from __future__ import annotations
 
+import math
 import operator
 from functools import reduce
-from math import sqrt
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -56,7 +56,7 @@ def prime_factors(n: int) -> list[int]:
         list of all prime factors of the given natural n.
     """
     idx = 2
-    while idx <= sqrt(n):
+    while idx <= math.sqrt(n):
         if n % idx == 0:
             lst = prime_factors(n // idx)
             lst.append(idx)

src/pymatgen/analysis/diffraction/neutron.py

@@ -3,8 +3,8 @@
 from __future__ import annotations
 
 import json
+import math
 import os
-from math import asin, cos, degrees, pi, radians, sin
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -96,7 +96,7 @@ def get_pattern(self, structure: Structure, scaled=True, two_theta_range=(0, 90)
         min_r, max_r = (
             (0, 2 / wavelength)
             if two_theta_range is None
-            else [2 * sin(radians(t / 2)) / wavelength for t in two_theta_range]
+            else [2 * math.sin(math.radians(t / 2)) / wavelength for t in two_theta_range]
         )
 
         # Obtain crystallographic reciprocal lattice points within range
@@ -137,12 +137,12 @@ def get_pattern(self, structure: Structure, scaled=True, two_theta_range=(0, 90)
 
         for hkl, g_hkl, ind, _ in sorted(recip_pts, key=lambda i: (i[1], -i[0][0], -i[0][1], -i[0][2])):
             # Force miller indices to be integers
-            hkl = [int(round(i)) for i in hkl]
+            hkl = [round(i) for i in hkl]
             if g_hkl != 0:
                 d_hkl = 1 / g_hkl
 
                 # Bragg condition
-                theta = asin(wavelength * g_hkl / 2)
+                theta = math.asin(wavelength * g_hkl / 2)
 
                 # s = sin(theta) / wavelength = 1 / 2d = |ghkl| / 2 (d =
                 # 1/|ghkl|)
@@ -158,15 +158,15 @@ def get_pattern(self, structure: Structure, scaled=True, two_theta_range=(0, 90)
                 # Structure factor = sum of atomic scattering factors (with
                 # position factor exp(2j * pi * g.r and occupancies).
                 # Vectorized computation.
-                f_hkl = np.sum(coeffs * occus * np.exp(2j * pi * g_dot_r) * dw_correction)
+                f_hkl = np.sum(coeffs * occus * np.exp(2j * np.pi * g_dot_r) * dw_correction)
 
                 # Lorentz polarization correction for hkl
-                lorentz_factor = 1 / (sin(theta) ** 2 * cos(theta))
+                lorentz_factor = 1 / (math.sin(theta) ** 2 * math.cos(theta))
 
                 # Intensity for hkl is modulus square of structure factor
                 i_hkl = (f_hkl * f_hkl.conjugate()).real
 
-                two_theta = degrees(2 * theta)
+                two_theta = math.degrees(2 * theta)
 
                 if is_hex:
                     # Use Miller-Bravais indices for hexagonal lattices

src/pymatgen/analysis/diffraction/xrd.py

@@ -3,8 +3,8 @@
 from __future__ import annotations
 
 import json
+import math
 import os
-from math import asin, cos, degrees, pi, radians, sin
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -158,7 +158,7 @@ def get_pattern(self, structure: Structure, scaled=True, two_theta_range=(0, 90)
         min_r, max_r = (
             (0, 2 / wavelength)
             if two_theta_range is None
-            else [2 * sin(radians(t / 2)) / wavelength for t in two_theta_range]
+            else [2 * math.sin(math.radians(t / 2)) / wavelength for t in two_theta_range]
         )
 
         # Obtain crystallographic reciprocal lattice points within range
@@ -201,10 +201,10 @@ def get_pattern(self, structure: Structure, scaled=True, two_theta_range=(0, 90)
 
         for hkl, g_hkl, ind, _ in sorted(recip_pts, key=lambda i: (i[1], -i[0][0], -i[0][1], -i[0][2])):
             # Force miller indices to be integers
-            hkl = [int(round(i)) for i in hkl]
+            hkl = [round(i) for i in hkl]
             if g_hkl != 0:
                 # Bragg condition
-                theta = asin(wavelength * g_hkl / 2)
+                theta = math.asin(wavelength * g_hkl / 2)
 
                 # s = sin(theta) / wavelength = 1 / 2d = |ghkl| / 2 (d =
                 # 1/|ghkl|)
@@ -235,15 +235,15 @@ def get_pattern(self, structure: Structure, scaled=True, two_theta_range=(0, 90)
                 # Structure factor = sum of atomic scattering factors (with
                 # position factor exp(2j * pi * g.r and occupancies).
                 # Vectorized computation.
-                f_hkl = np.sum(fs * occus * np.exp(2j * pi * g_dot_r) * dw_correction)
+                f_hkl = np.sum(fs * occus * np.exp(2j * np.pi * g_dot_r) * dw_correction)
 
                 # Lorentz polarization correction for hkl
-                lorentz_factor = (1 + cos(2 * theta) ** 2) / (sin(theta) ** 2 * cos(theta))
+                lorentz_factor = (1 + math.cos(2 * theta) ** 2) / (math.sin(theta) ** 2 * math.cos(theta))
 
                 # Intensity for hkl is modulus square of structure factor
                 i_hkl = (f_hkl * f_hkl.conjugate()).real
 
-                two_theta = degrees(2 * theta)
+                two_theta = math.degrees(2 * theta)
 
                 if is_hex:
                     # Use Miller-Bravais indices for hexagonal lattices

src/pymatgen/analysis/elasticity/elastic.py

@@ -214,7 +214,7 @@ def directional_poisson_ratio(self, n: ArrayLike, m: ArrayLike, tol: float = 1e-
             tol (float): tolerance for testing of orthogonality
         """
         n, m = get_uvec(n), get_uvec(m)
-        if not np.abs(np.dot(n, m)) < tol:
+        if np.abs(np.dot(n, m)) >= tol:
             raise ValueError("n and m must be orthogonal")
         v = self.compliance_tensor.einsum_sequence([n] * 2 + [m] * 2)
         v *= -1 / self.compliance_tensor.einsum_sequence([n] * 4)
@@ -907,7 +907,7 @@ def find_eq_stress(strains, stresses, tol: float = 1e-10):
     eq_stress = stress_array[np.all(abs(strain_array) < tol, axis=(1, 2))]
 
     if eq_stress.size != 0:
-        all_same = (abs(eq_stress - eq_stress[0]) < 1e-8).all()
+        all_same = np.allclose(eq_stress, eq_stress[0], atol=1e-8, rtol=0)
         if len(eq_stress) > 1 and not all_same:
             raise ValueError(
                 "Multiple stresses found for equilibrium strain"

src/pymatgen/analysis/interface_reactions.py

@@ -150,6 +150,7 @@ def get_kinks(self) -> list[tuple[int, float, float, Reaction, float]]:
         critical_comp = self.pd.get_critical_compositions(self.comp1, self.comp2)
         x_kink, energy_kink, react_kink, energy_per_rxt_formula = [], [], [], []
 
+        # TODO: perhaps a bad idea to use full equality to compare coords
         if (c1_coord == c2_coord).all():
             x_kink = [0, 1]
             energy_kink = [self._get_energy(x) for x in x_kink]

src/pymatgen/analysis/local_env.py

@@ -127,8 +127,8 @@ def nearest_key(sorted_vals: list[int], skey: int) -> int:
                 continue
 
             el = site.specie.symbol
-            oxi_state = int(round(site.specie.oxi_state))
-            coord_no = int(round(vnn.get_cn(self._structure, idx)))
+            oxi_state = round(site.specie.oxi_state)
+            coord_no = round(vnn.get_cn(self._structure, idx))
             try:
                 tab_oxi_states = sorted(map(int, _ION_RADII[el]))
                 oxi_state = nearest_key(tab_oxi_states, oxi_state)
@@ -2888,9 +2888,9 @@ def get_order_parameters(
         if tol < 0.0:
             raise ValueError("Negative tolerance for weighted solid angle!")
 
-        left_of_unity = 1 - 1.0e-12
+        left_of_unity = 1 - 1e-12
         # The following threshold has to be adapted to non-Angstrom units.
-        very_small = 1.0e-12
+        very_small = 1e-12
         fac_bcc = 1 / math.exp(-0.5)
 
         # Find central site and its neighbors.
@@ -3330,7 +3330,7 @@ def get_order_parameters(
                     for j in range(n_neighbors):
                         ops[idx] += sum(qsp_theta[idx][j])
                         tmp_norm += float(sum(norms[idx][j]))
-                    ops[idx] = ops[idx] / tmp_norm if tmp_norm > 1.0e-12 else None  # type: ignore[operator]
+                    ops[idx] = ops[idx] / tmp_norm if tmp_norm > 1e-12 else None  # type: ignore[operator]
 
                 elif typ in {
                     "T",
@@ -3357,7 +3357,7 @@ def get_order_parameters(
                         for j in range(n_neighbors):
                             for k in range(len(qsp_theta[idx][j])):
                                 qsp_theta[idx][j][k] = (
-                                    qsp_theta[idx][j][k] / norms[idx][j][k] if norms[idx][j][k] > 1.0e-12 else 0.0
+                                    qsp_theta[idx][j][k] / norms[idx][j][k] if norms[idx][j][k] > 1e-12 else 0.0
                                 )
                             ops[idx] = max(qsp_theta[idx][j]) if j == 0 else max(ops[idx], *qsp_theta[idx][j])
 
@@ -3436,7 +3436,7 @@ class BrunnerNNReciprocal(NearNeighbors):
     largest reciprocal gap in interatomic distances.
     """
 
-    def __init__(self, tol: float = 1.0e-4, cutoff=8.0) -> None:
+    def __init__(self, tol: float = 1e-4, cutoff=8.0) -> None:
         """
         Args:
             tol (float): tolerance parameter for bond determination
@@ -3511,7 +3511,7 @@ class BrunnerNNRelative(NearNeighbors):
     of largest relative gap in interatomic distances.
     """
 
-    def __init__(self, tol: float = 1.0e-4, cutoff=8.0) -> None:
+    def __init__(self, tol: float = 1e-4, cutoff=8.0) -> None:
         """
         Args:
             tol (float): tolerance parameter for bond determination
@@ -3587,7 +3587,7 @@ class BrunnerNNReal(NearNeighbors):
     largest gap in interatomic distances.
     """
 
-    def __init__(self, tol: float = 1.0e-4, cutoff=8.0) -> None:
+    def __init__(self, tol: float = 1e-4, cutoff=8.0) -> None:
         """
         Args:
             tol (float): tolerance parameter for bond determination
@@ -3748,7 +3748,7 @@ def get_nn_info(self, structure: Structure, n: int):
         # calculate mean fictive ionic radius
         mefir = _get_mean_fictive_ionic_radius(firs)
 
-        # # iteratively solve MEFIR; follows equation 4 in Hoppe's EconN paper
+        # iteratively solve MEFIR; follows equation 4 in Hoppe's EconN paper
         prev_mefir = float("inf")
         while abs(prev_mefir - mefir) > 1e-4:
             # this is guaranteed to converge

src/pymatgen/analysis/magnetism/jahnteller.py

@@ -2,6 +2,7 @@
 
 from __future__ import annotations
 
+import math
 import os
 import warnings
 from typing import TYPE_CHECKING, Literal, cast
@@ -444,9 +445,9 @@ def _estimate_spin_state(
         # WARNING! this heuristic has not been robustly tested or benchmarked
         # using 'diff*0.25' as arbitrary measure, if known magmom is
         # too far away from expected value, we don't try to classify it
-        if known_magmom > mu_so_high or abs(mu_so_high - known_magmom) < diff * 0.25:
+        if known_magmom > mu_so_high or math.isclose(mu_so_high, known_magmom, abs_tol=diff * 0.25, rel_tol=0):
             return "high"
-        if known_magmom < mu_so_low or abs(mu_so_low - known_magmom) < diff * 0.25:
+        if known_magmom < mu_so_low or math.isclose(mu_so_low, known_magmom, abs_tol=diff * 0.25, rel_tol=0):
             return "low"
         return "unknown"
 

src/pymatgen/analysis/phase_diagram.py

@@ -1014,7 +1014,9 @@ def get_transition_chempots(self, element):
 
         clean_pots = []
         for c in sorted(critical_chempots):
-            if len(clean_pots) == 0 or abs(c - clean_pots[-1]) > PhaseDiagram.numerical_tol:
+            if len(clean_pots) == 0 or not math.isclose(
+                c, clean_pots[-1], abs_tol=PhaseDiagram.numerical_tol, rel_tol=0
+            ):
                 clean_pots.append(c)
         clean_pots.reverse()
         return tuple(clean_pots)
@@ -1996,7 +1998,11 @@ def fmt(fl):
 
                     x = coeffs[-1]
 
-                    if all(c >= -tol for c in coeffs) and (abs(sum(coeffs[:-1]) - 1) < tol) and (tol < x < 1 - tol):
+                    if (
+                        all(c >= -tol for c in coeffs)
+                        and (math.isclose(sum(coeffs[:-1]), 1, abs_tol=tol, rel_tol=0))
+                        and (tol < x < 1 - tol)
+                    ):
                         c1 = x / r1.num_atoms
                         c2 = (1 - x) / r2.num_atoms
                         factor = 1 / (c1 + c2)

src/pymatgen/analysis/piezo.py

@@ -38,7 +38,7 @@ def __new__(cls, input_array: ArrayLike, tol: float = 1e-3) -> Self:
                 representing the piezo tensor
         """
         obj = super().__new__(cls, input_array, check_rank=3)
-        if not (obj - np.transpose(obj, (0, 2, 1)) < tol).all():
+        if not np.allclose(obj, np.transpose(obj, (0, 2, 1)), atol=tol, rtol=0):
             warnings.warn("Input piezo tensor does not satisfy standard symmetries")
         return obj.view(cls)
 

src/pymatgen/analysis/piezo_sensitivity.py

@@ -181,7 +181,7 @@ def __init__(self, structure: Structure, ist, pointops, tol: float = 1e-3):
         self.IST_operations: list[list[list]] = []
 
         obj = self.ist
-        if not (obj - np.transpose(obj, (0, 1, 3, 2)) < tol).all():
+        if not np.allclose(obj, np.transpose(obj, (0, 1, 3, 2)), atol=tol, rtol=0):
             warnings.warn("Input internal strain tensor does not satisfy standard symmetries")
 
     def get_IST_operations(self, opstol=1e-3) -> list[list[list]]:

src/pymatgen/analysis/quasirrho.py

@@ -10,7 +10,7 @@
 
 from __future__ import annotations
 
-from math import isclose
+import math
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -221,7 +221,7 @@ def _get_quasirrho_thermo(
         linear = True
         for coord in coords[1:]:
             theta = abs(np.dot(coord - coords[0], v0) / np.linalg.norm(coord - coords[0]) / np.linalg.norm(v0))
-            if not isclose(theta, 1, abs_tol=1e-4):
+            if not math.isclose(theta, 1, abs_tol=1e-4):
                 linear = False
 
         # Rotational component of Entropy and Energy