diff --git a/pyproject.toml b/pyproject.toml
index 0b47f13418e..5193cf7ddf4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -56,8 +56,8 @@ classifiers = [
 dependencies = [
     "joblib>=1",
     "matplotlib>=3.8",
-    "monty>=2024.7.29",
-    "networkx>=3",  # PR #4116
+    "monty>=2024.10.21",
+    "networkx>=2.7",  # PR #4116
     "palettable>=3.3.3",
     "pandas>=2",
     "plotly>=4.5.0",
@@ -69,7 +69,7 @@ dependencies = [
     # https://github.com/scipy/scipy/issues/21052
     "scipy>=1.14.1; platform_system == 'Windows'",
     "spglib>=2.5.0",
-    "sympy>=1.2",
+    "sympy>=1.3",  # PR #4116
     "tabulate>=0.9",
     "tqdm>=4.60",
     "uncertainties>=3.1.4",
@@ -87,40 +87,34 @@ Issues = "https://github.com/materialsproject/pymatgen/issues"
 Pypi = "https://pypi.org/project/pymatgen"
 
 [project.optional-dependencies]
+# PR4128: netcdf4 1.7.[0/1] yanked, 1.7.1.post[1/2]/1.7.2 cause CI error
+abinit = ["netcdf4>=1.6.5,!=1.7.1.post1,!=1.7.1.post2,!=1.7.2"]
 ase = ["ase>=3.23.0"]
-# tblite only support Python 3.12+ through conda-forge
-# https://github.com/tblite/tblite/issues/175
-tblite = ["tblite[ase]>=0.3.0; python_version<'3.12'"]
-vis = ["vtk>=6.0.0"]
-abinit = ["netcdf4>=1.7.1"]
-mlp = ["chgnet>=0.3.8", "matgl>=1.1.3"]
-electronic_structure = ["fdint>=2.0.2"]
 ci = ["pytest-cov>=4", "pytest-split>=0.8", "pytest>=8"]
 docs = ["invoke", "sphinx", "sphinx_markdown_builder", "sphinx_rtd_theme"]
+electronic_structure = ["fdint>=2.0.2"]
+mlp = ["chgnet>=0.3.8", "matgl>=1.1.3"]
+numba = ["numba>=0.55"]
+numpy-v1 = ["numpy>=1.25.0,<2"] # Test NP1 on Windows (quite buggy ATM)
 optional = [
-    "ase>=3.23.0",
+    "pymatgen[abinit,ase,mlp,tblite]",
     "beautifulsoup4",
     # BoltzTraP2 build fails on Windows GitHub runners
     "BoltzTraP2>=24.9.4 ; platform_system != 'Windows'",
     "chemview>=0.6",
-    "chgnet>=0.3.8",
     "f90nml>=1.1.2",
     "galore>=0.6.1",
     "h5py>=3.11.0",
+    "hiphive>=1.3.1",
     "jarvis-tools>=2020.7.14",
-    "matgl>=1.1.3",
     "matplotlib>=3.8",
-    "netCDF4>=1.6.5",
     "phonopy>=2.23",
     "seekpath>=2.0.1",
-    # tblite only support Python 3.12+ through conda-forge
-    # https://github.com/tblite/tblite/issues/175
-    "hiphive>=1.3.1",
-    "openbabel-wheel>=3.1.1.20",
-    "tblite[ase]>=0.3.0; platform_system=='Linux' and python_version<'3.12'",
 ]
-numba = ["numba>=0.55"]
-numpy-v1 = ["numpy>=1.25.0,<2"] # Test NP1 on Windows (quite buggy ATM)
+# tblite only support Python 3.12+ through conda-forge
+# https://github.com/tblite/tblite/issues/175
+tblite = [ "tblite[ase]>=0.3.0; platform_system=='Linux' and python_version<'3.12'"]
+vis = ["vtk>=6.0.0"]
 
 [project.scripts]
 pmg = "pymatgen.cli.pmg:main"
diff --git a/src/pymatgen/analysis/adsorption.py b/src/pymatgen/analysis/adsorption.py
index 3e953ec00e7..fcb96a935e8 100644
--- a/src/pymatgen/analysis/adsorption.py
+++ b/src/pymatgen/analysis/adsorption.py
@@ -9,8 +9,6 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from matplotlib import patches
-from matplotlib.path import Path
 from monty.serialization import loadfn
 from scipy.spatial import Delaunay
 
@@ -664,6 +662,10 @@ def plot_slab(
         decay (float): how the alpha-value decays along the z-axis
         inverse (bool): invert z axis to plot opposite surface
     """
+    # Expensive import (PR4128)
+    from matplotlib import patches
+    from matplotlib.path import Path
+
     orig_slab = slab.copy()
     slab = reorient_z(slab)
     orig_cell = slab.lattice.matrix.copy()
diff --git a/src/pymatgen/analysis/interfaces/coherent_interfaces.py b/src/pymatgen/analysis/interfaces/coherent_interfaces.py
index cce7bb2f2e0..0a0633e273a 100644
--- a/src/pymatgen/analysis/interfaces/coherent_interfaces.py
+++ b/src/pymatgen/analysis/interfaces/coherent_interfaces.py
@@ -6,7 +6,6 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from numpy.testing import assert_allclose
 from scipy.linalg import polar
 
 from pymatgen.analysis.elasticity.strain import Deformation
@@ -101,21 +100,15 @@ def _find_matches(self) -> None:
         for match in self.zsl_matches:
             xform = get_2d_transform(film_vectors, match.film_vectors)
             strain, _rot = polar(xform)
-            assert_allclose(
-                strain,
-                np.round(strain),
-                atol=1e-12,
-                err_msg="Film lattice vectors changed during ZSL match, check your ZSL Generator parameters",
-            )
+            if not np.allclose(strain, np.round(strain), atol=1e-12):
+                raise ValueError("Film lattice vectors changed during ZSL match, check your ZSL Generator parameters")
 
             xform = get_2d_transform(substrate_vectors, match.substrate_vectors)
             strain, _rot = polar(xform)
-            assert_allclose(
-                strain,
-                strain.astype(int),
-                atol=1e-12,
-                err_msg="Substrate lattice vectors changed during ZSL match, check your ZSL Generator parameters",
-            )
+            if not np.allclose(strain, strain.astype(int), atol=1e-12):
+                raise ValueError(
+                    "Substrate lattice vectors changed during ZSL match, check your ZSL Generator parameters"
+                )
 
     def _find_terminations(self):
         """Find all terminations."""
@@ -226,18 +219,12 @@ def get_interfaces(
             ).astype(int)
             film_sl_slab = film_slab.copy()
             film_sl_slab.make_supercell(super_film_transform)
-            assert_allclose(
-                film_sl_slab.lattice.matrix[2],
-                film_slab.lattice.matrix[2],
-                atol=1e-08,
-                err_msg="2D transformation affected C-axis for Film transformation",
-            )
-            assert_allclose(
-                film_sl_slab.lattice.matrix[:2],
-                match.film_sl_vectors,
-                atol=1e-08,
-                err_msg="Transformation didn't make proper supercell for film",
-            )
+            if not np.allclose(film_sl_slab.lattice.matrix[2], film_slab.lattice.matrix[2], atol=1e-08):
+                raise ValueError(
+                    "2D transformation affected C-axis for Film transformation",
+                )
+            if not np.allclose(film_sl_slab.lattice.matrix[:2], match.film_sl_vectors, atol=1e-08):
+                raise ValueError("Transformation didn't make proper supercell for film")
 
             # Build substrate superlattice
             super_sub_transform = np.round(
@@ -245,18 +232,10 @@ def get_interfaces(
             ).astype(int)
             sub_sl_slab = sub_slab.copy()
             sub_sl_slab.make_supercell(super_sub_transform)
-            assert_allclose(
-                sub_sl_slab.lattice.matrix[2],
-                sub_slab.lattice.matrix[2],
-                atol=1e-08,
-                err_msg="2D transformation affected C-axis for Film transformation",
-            )
-            assert_allclose(
-                sub_sl_slab.lattice.matrix[:2],
-                match.substrate_sl_vectors,
-                atol=1e-08,
-                err_msg="Transformation didn't make proper supercell for substrate",
-            )
+            if not np.allclose(sub_sl_slab.lattice.matrix[2], sub_slab.lattice.matrix[2], atol=1e-08):
+                raise ValueError("2D transformation affected C-axis for Film transformation")
+            if not np.allclose(sub_sl_slab.lattice.matrix[:2], match.substrate_sl_vectors, atol=1e-08):
+                raise ValueError("Transformation didn't make proper supercell for substrate")
 
             # Add extra info
             match_dict = match.as_dict()
diff --git a/src/pymatgen/core/__init__.py b/src/pymatgen/core/__init__.py
index 840b16e9b69..68343cac523 100644
--- a/src/pymatgen/core/__init__.py
+++ b/src/pymatgen/core/__init__.py
@@ -5,7 +5,7 @@
 import os
 import warnings
 from importlib.metadata import PackageNotFoundError, version
-from typing import Any
+from typing import TYPE_CHECKING
 
 from ruamel.yaml import YAML
 
@@ -17,10 +17,14 @@
 from pymatgen.core.structure import IMolecule, IStructure, Molecule, PeriodicNeighbor, SiteCollection, Structure
 from pymatgen.core.units import ArrayWithUnit, FloatWithUnit, Unit
 
+if TYPE_CHECKING:
+    from typing import Any
+
 __author__ = "Pymatgen Development Team"
 __email__ = "pymatgen@googlegroups.com"
 __maintainer__ = "Shyue Ping Ong, Matthew Horton, Janosh Riebesell"
 __maintainer_email__ = "shyuep@gmail.com"
+
 try:
     __version__ = version("pymatgen")
 except PackageNotFoundError:  # pragma: no cover
diff --git a/src/pymatgen/core/interface.py b/src/pymatgen/core/interface.py
index 823d009849f..4bb1d76eec8 100644
--- a/src/pymatgen/core/interface.py
+++ b/src/pymatgen/core/interface.py
@@ -14,7 +14,6 @@
 
 import numpy as np
 from monty.fractions import lcm
-from numpy.testing import assert_allclose
 from scipy.cluster.hierarchy import fcluster, linkage
 from scipy.spatial.distance import squareform
 
@@ -2784,10 +2783,16 @@ def from_slabs(
             substrate_slab = substrate_slab.get_orthogonal_c_slab()
         if isinstance(film_slab, Slab):
             film_slab = film_slab.get_orthogonal_c_slab()
-        assert_allclose(film_slab.lattice.alpha, 90, 0.1)
-        assert_allclose(film_slab.lattice.beta, 90, 0.1)
-        assert_allclose(substrate_slab.lattice.alpha, 90, 0.1)
-        assert_allclose(substrate_slab.lattice.beta, 90, 0.1)
+
+        if not math.isclose(film_slab.lattice.alpha, 90, rel_tol=0.1) or not math.isclose(
+            film_slab.lattice.beta, 90, rel_tol=0.1
+        ):
+            raise ValueError("The lattice angles alpha or beta of the film slab are not approximately 90 degrees.")
+
+        if not math.isclose(substrate_slab.lattice.alpha, 90, rel_tol=0.1) or not math.isclose(
+            substrate_slab.lattice.beta, 90, rel_tol=0.1
+        ):
+            raise ValueError("The lattice angles alpha or beta of the substrate slab are not approximately 90 degrees.")
 
         # Ensure sub is right-handed
         # IE sub has surface facing "up"
diff --git a/src/pymatgen/core/operations.py b/src/pymatgen/core/operations.py
index d395c2c3d3b..5c6ea895a12 100644
--- a/src/pymatgen/core/operations.py
+++ b/src/pymatgen/core/operations.py
@@ -449,7 +449,7 @@ def as_xyz_str(self) -> str:
         Only works for integer rotation matrices.
         """
         # Check for invalid rotation matrix
-        if not np.all(np.isclose(self.rotation_matrix, np.round(self.rotation_matrix))):
+        if not np.allclose(self.rotation_matrix, np.round(self.rotation_matrix)):
             warnings.warn("Rotation matrix should be integer")
 
         return transformation_to_string(
diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py
index 06a58f77a3c..a7bc74ac2d6 100644
--- a/src/pymatgen/core/trajectory.py
+++ b/src/pymatgen/core/trajectory.py
@@ -15,7 +15,6 @@
 from monty.json import MSONable
 
 from pymatgen.core.structure import Composition, DummySpecies, Element, Lattice, Molecule, Species, Structure
-from pymatgen.io.ase import AseAtomsAdaptor
 
 if TYPE_CHECKING:
     from collections.abc import Iterator
@@ -580,9 +579,12 @@ def from_file(cls, filename: str | Path, constant_lattice: bool = True, **kwargs
             try:
                 from ase.io.trajectory import Trajectory as AseTrajectory
 
+                from pymatgen.io.ase import AseAtomsAdaptor
+
                 ase_traj = AseTrajectory(filename)
                 # Periodic boundary conditions should be the same for all frames so just check the first
                 pbc = ase_traj[0].pbc
+
                 if any(pbc):
                     structures = [AseAtomsAdaptor.get_structure(atoms) for atoms in ase_traj]
                 else:
diff --git a/src/pymatgen/io/vasp/outputs.py b/src/pymatgen/io/vasp/outputs.py
index 921aa71a6a4..d8603e03b50 100644
--- a/src/pymatgen/io/vasp/outputs.py
+++ b/src/pymatgen/io/vasp/outputs.py
@@ -24,7 +24,6 @@
 from monty.json import MSONable, jsanitize
 from monty.os.path import zpath
 from monty.re import regrep
-from numpy.testing import assert_allclose
 from tqdm import tqdm
 
 from pymatgen.core import Composition, Element, Lattice, Structure
@@ -4975,8 +4974,8 @@ def __init__(
 
                     if i_spin == 0:
                         self.kpoints.append(kpoint)
-                    else:
-                        assert_allclose(self.kpoints[i_nk], kpoint)
+                    elif not np.allclose(self.kpoints[i_nk], kpoint, rtol=1e-7, atol=0):
+                        raise ValueError(f"kpoints of {i_nk=} mismatch")
 
                     if verbose:
                         print(f"kpoint {i_nk: 4} with {nplane: 5} plane waves at {kpoint}")
diff --git a/src/pymatgen/symmetry/settings.py b/src/pymatgen/symmetry/settings.py
index c8d66a89c7c..ae407b7a412 100644
--- a/src/pymatgen/symmetry/settings.py
+++ b/src/pymatgen/symmetry/settings.py
@@ -7,8 +7,6 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from sympy import Matrix
-from sympy.parsing.sympy_parser import parse_expr
 
 from pymatgen.core.lattice import Lattice
 from pymatgen.core.operations import MagSymmOp, SymmOp
@@ -100,6 +98,10 @@ def parse_transformation_string(
         Returns:
             tuple[list[list[float]] | np.ndarray, list[float]]: transformation matrix & vector
         """
+        # Import sympy is expensive (PR4128)
+        from sympy import Matrix
+        from sympy.parsing.sympy_parser import parse_expr
+
         try:
             a, b, c = np.eye(3)
             b_change, o_shift = transformation_string.split(";")
diff --git a/tests/files/performance/import_time_linux.json b/tests/files/performance/import_time_linux.json
new file mode 100644
index 00000000000..3492c93a840
--- /dev/null
+++ b/tests/files/performance/import_time_linux.json
@@ -0,0 +1,19 @@
+{
+    "from pymatgen.core.bonds import CovalentBond": 289.5851116666108,
+    "from pymatgen.core.composition import Composition": 292.8479909999548,
+    "from pymatgen.core.interface import Interface": 969.5693099999593,
+    "from pymatgen.core.ion import Ion": 291.07530133334575,
+    "from pymatgen.core.lattice import Lattice": 288.8340153333881,
+    "from pymatgen.core.libxcfunc import LibxcFunc": 293.4184753333587,
+    "from pymatgen.core.molecular_orbitals import MolecularOrbitals": 294.19796566658835,
+    "from pymatgen.core.operations import SymmOp": 296.4627546666634,
+    "from pymatgen.core.periodic_table import Element": 295.95872066662804,
+    "from pymatgen.core.sites import Site": 292.66485499999817,
+    "from pymatgen.core.spectrum import Spectrum": 486.72776566669046,
+    "from pymatgen.core.structure import Structure": 291.01618733333606,
+    "from pymatgen.core.surface import Slab": 301.90875833329756,
+    "from pymatgen.core.tensors import Tensor": 304.27744800003137,
+    "from pymatgen.core.trajectory import Trajectory": 300.45536066666045,
+    "from pymatgen.core.units import Unit": 305.4779056666348,
+    "from pymatgen.core.xcfunc import XcFunc": 309.1085626666275
+}
\ No newline at end of file
diff --git a/tests/files/performance/import_time_macos.json b/tests/files/performance/import_time_macos.json
new file mode 100644
index 00000000000..9ac1e63fc9d
--- /dev/null
+++ b/tests/files/performance/import_time_macos.json
@@ -0,0 +1,19 @@
+{
+    "from pymatgen.core.bonds import CovalentBond": 321.8502360000457,
+    "from pymatgen.core.composition import Composition": 292.35445800009074,
+    "from pymatgen.core.interface import Interface": 855.005861000033,
+    "from pymatgen.core.ion import Ion": 240.8930970000256,
+    "from pymatgen.core.lattice import Lattice": 329.09868066659936,
+    "from pymatgen.core.libxcfunc import LibxcFunc": 306.6966386666839,
+    "from pymatgen.core.molecular_orbitals import MolecularOrbitals": 281.78087466661356,
+    "from pymatgen.core.operations import SymmOp": 299.9741943333447,
+    "from pymatgen.core.periodic_table import Element": 293.6565829999533,
+    "from pymatgen.core.sites import Site": 280.3443330000543,
+    "from pymatgen.core.spectrum import Spectrum": 459.20266666666976,
+    "from pymatgen.core.structure import Structure": 265.4675833332476,
+    "from pymatgen.core.surface import Slab": 306.0919996667053,
+    "from pymatgen.core.tensors import Tensor": 310.54281933325,
+    "from pymatgen.core.trajectory import Trajectory": 335.25658333329983,
+    "from pymatgen.core.units import Unit": 294.03472200003006,
+    "from pymatgen.core.xcfunc import XcFunc": 309.3993196666058
+}
\ No newline at end of file
diff --git a/tests/files/performance/import_time_windows.json b/tests/files/performance/import_time_windows.json
new file mode 100644
index 00000000000..270b73a9d65
--- /dev/null
+++ b/tests/files/performance/import_time_windows.json
@@ -0,0 +1,19 @@
+{
+    "from pymatgen.core.bonds import CovalentBond": 443.7560000000455,
+    "from pymatgen.core.composition import Composition": 441.06553333335796,
+    "from pymatgen.core.interface import Interface": 1828.751033333295,
+    "from pymatgen.core.ion import Ion": 443.58053333333675,
+    "from pymatgen.core.lattice import Lattice": 445.729999999988,
+    "from pymatgen.core.libxcfunc import LibxcFunc": 459.24773333338936,
+    "from pymatgen.core.molecular_orbitals import MolecularOrbitals": 440.4825999999957,
+    "from pymatgen.core.operations import SymmOp": 440.62226666665083,
+    "from pymatgen.core.periodic_table import Element": 441.64050000002436,
+    "from pymatgen.core.sites import Site": 442.1802333333744,
+    "from pymatgen.core.spectrum import Spectrum": 737.3025000000174,
+    "from pymatgen.core.structure import Structure": 445.0546333332568,
+    "from pymatgen.core.surface import Slab": 463.0683333333157,
+    "from pymatgen.core.tensors import Tensor": 463.5761666666743,
+    "from pymatgen.core.trajectory import Trajectory": 443.9995333333779,
+    "from pymatgen.core.units import Unit": 446.352766666602,
+    "from pymatgen.core.xcfunc import XcFunc": 469.42599999996065
+}
\ No newline at end of file
diff --git a/tests/io/abinit/test_netcdf.py b/tests/io/abinit/test_netcdf.py
index 77c313c53c3..ddc0e95ec50 100644
--- a/tests/io/abinit/test_netcdf.py
+++ b/tests/io/abinit/test_netcdf.py
@@ -92,6 +92,8 @@ def test_read_fe(self):
                 ref_magmom_collinear = [-0.5069359730980665]
                 path = os.path.join(tmp_dir, "Fe_magmoms_collinear_GSR.nc")
 
+                # TODO: PR4128, EtsfReader would fail in Ubuntu CI with netCDF4 > 1.6.5
+                # Need someone with knowledge in netCDF4 to fix it
                 with EtsfReader(path) as data:
                     structure = data.read_structure()
                     assert structure.site_properties["magmom"] == ref_magmom_collinear
diff --git a/tests/performance/test_import_time.py b/tests/performance/test_import_time.py
new file mode 100644
index 00000000000..9c9f9d5daf1
--- /dev/null
+++ b/tests/performance/test_import_time.py
@@ -0,0 +1,130 @@
+"""
+Test the import time of several important modules.
+
+NOTE:
+    - Toggle the "GEN_REF_TIME" to generate reference import time.
+
+Last update: 2024-10-26 (PR 4128)
+
+Runner specs:
+Linux:   4 CPU
+Windows: 4 CPU
+macOS:   3 CPU (M1)
+"""
+
+from __future__ import annotations
+
+import json
+import os
+import subprocess
+import time
+import warnings
+from typing import TYPE_CHECKING
+
+import pytest
+
+from pymatgen.util.testing import TEST_FILES_DIR
+
+if TYPE_CHECKING:
+    from typing import Literal
+
+if not os.getenv("CI"):
+    pytest.skip("ref time only comparable in CI runner", allow_module_level=True)
+
+# NOTE: Toggle this to generate reference import time
+GEN_REF_TIME: bool = False
+
+MODULES_TO_TEST: tuple[str, ...] = (
+    "from pymatgen.core.bonds import CovalentBond",
+    "from pymatgen.core.composition import Composition",
+    "from pymatgen.core.interface import Interface",
+    "from pymatgen.core.ion import Ion",
+    "from pymatgen.core.lattice import Lattice",
+    "from pymatgen.core.libxcfunc import LibxcFunc",
+    "from pymatgen.core.molecular_orbitals import MolecularOrbitals",
+    "from pymatgen.core.operations import SymmOp",
+    "from pymatgen.core.periodic_table import Element",
+    "from pymatgen.core.sites import Site",
+    "from pymatgen.core.spectrum import Spectrum",
+    "from pymatgen.core.structure import Structure",
+    "from pymatgen.core.surface import Slab",
+    "from pymatgen.core.tensors import Tensor",
+    "from pymatgen.core.trajectory import Trajectory",
+    "from pymatgen.core.units import Unit",
+    "from pymatgen.core.xcfunc import XcFunc",
+)
+
+# Get runner OS and reference file
+RUNNER_OS: Literal["linux", "windows", "macos"] = os.getenv("RUNNER_OS", "").lower()  # type: ignore[assignment]
+assert RUNNER_OS in {"linux", "windows", "macos"}
+
+# Skip test on macOS due to inconsistent import time results
+if RUNNER_OS == "macos":
+    pytest.skip("Import time measurements are unstable on macOS", allow_module_level=True)
+
+REF_FILE: str = f"{TEST_FILES_DIR}/performance/import_time_{RUNNER_OS}.json"
+
+
+@pytest.mark.skipif(not GEN_REF_TIME, reason="Set GEN_REF_TIME to generate reference import time.")
+def test_get_ref_import_time() -> None:
+    """A dummy test that would always fail, used to generate copyable reference time."""
+    import_times: dict[str, float] = {
+        module_import_cmd: _measure_import_time_in_ms(module_import_cmd) for module_import_cmd in MODULES_TO_TEST
+    }
+
+    # Print a copyable JSON format for easy reference updating
+    print("\nCopyable import time JSON:")
+    print(json.dumps(import_times, indent=4))
+
+    pytest.fail("Reference import times generated. Copy from output to update JSON file.")
+
+
+@pytest.mark.skipif(GEN_REF_TIME, reason="Generating reference import time.")
+def test_import_time(grace_percent: float = 0.5, hard_percent: float = 1.0) -> None:
+    """Test the import time of core modules to avoid performance regression.
+
+    Args:
+        grace_percent (float): Maximum allowed percentage increase in import time
+            before a warning is raised.
+        hard_percent (float): Maximum allowed percentage increase in import time
+            before the test fails.
+    """
+
+    with open(REF_FILE, encoding="utf-8") as file:
+        ref_import_times: dict[str, float] = json.load(file)
+
+    for module_import_cmd, ref_time in ref_import_times.items():
+        current_time: float = _measure_import_time_in_ms(module_import_cmd)
+
+        # Calculate thresholds for grace and hard limits
+        grace_threshold = ref_time * (1 + grace_percent)
+        hard_threshold = ref_time * (1 + hard_percent)
+
+        if current_time > grace_threshold:
+            if current_time > hard_threshold:
+                pytest.fail(f"{module_import_cmd} import too slow at {current_time:.2f} ms! {hard_threshold=:.2f} ms")
+            else:
+                warnings.warn(
+                    f"{module_import_cmd} import slightly slower than reference: {grace_threshold=:.2f} ms",
+                    stacklevel=2,
+                )
+
+
+def _measure_import_time_in_ms(module_import_cmd: str, count: int = 3) -> float:
+    """Measure import time of a module in milliseconds across several runs.
+
+    Args:
+        module_import_cmd (str): The module import command.
+        count (int): Number of runs to average.
+
+    Returns:
+        float: Import time in milliseconds.
+    """
+    total_time = 0.0
+
+    for _ in range(count):
+        start_time = time.perf_counter_ns()
+        subprocess.run(["python", "-c", f"{module_import_cmd}"], check=True)
+        total_time += time.perf_counter_ns() - start_time
+
+    return (total_time / count) / 1e6