Merge pull request #155 from spglib/numpy20

Migrate to numpy 2.0

.github/dependabot.yml

@@ -3,9 +3,9 @@ updates:
   - package-ecosystem: "pip"
     directory: "/"
     schedule:
-      interval: "weekly"
+      interval: "monthly"
 
   - package-ecosystem: "github-actions"
     directory: "/"
     schedule:
-      interval: "weekly"
+      interval: "monthly"

.github/workflows/testing.yml

@@ -11,8 +11,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        # Python 3.11 is blocked by the following issue: https://github.com/h5py/h5py/issues/2146
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v4

.pre-commit-config.yaml

@@ -17,7 +17,7 @@ repos:
     - id: black
   # linter
   - repo: https://github.com/PyCQA/flake8
-    rev: 7.0.0
+    rev: 7.1.0
     hooks:
     - id: flake8
       exclude: ^docs/

docs/changelog.md

@@ -1,5 +1,9 @@
 # Change Log
 
+## v0.3.5 (18 Jun. 2024)
+
+- Make compatible with NumPy 2.0.0
+
 ## v0.3.3 (11 Jan. 2023)
 - Improve comparison between linear subspaces [[#53]](https://github.com/spglib/spgrep/pull/53)
     - Use Grassmann distance to measure linear subspaces: {func}`spgrep.utils.grassmann_distance`

docs/development/development.md

@@ -48,8 +48,7 @@ vim docs/changelog.md
 
 # Push with tag
 git tag <next-version>
-git push origin main
-git push origin --tags
+git push origin <next-version>
 ```
 
 ## Subpages

examples/lattice_vibration.py

@@ -37,7 +37,7 @@ def get_displacements_representation(
         perm_i = permutations[i]
         shifts[i] = positions @ Ri.T + vi[None, :] - positions[perm_i]
 
-    perm_rep = np.zeros((little_order, num_atoms, num_atoms), dtype=np.complex_)
+    perm_rep = np.zeros((little_order, num_atoms, num_atoms), dtype=np.complex128)
     for i, Ri in enumerate(little_rotations):
         for kappa in range(num_atoms):
             kappa2 = permutations[i, kappa]
@@ -48,7 +48,7 @@ def get_displacements_representation(
     # Rotation matrix in cartesian (order, 3, 3)
     A = np.transpose(lattice)  # column-wise lattice vectors
     Ainv = np.linalg.inv(A)
-    rotation_rep = np.array([A @ r @ Ainv for r in little_rotations], dtype=np.complex_)
+    rotation_rep = np.array([A @ r @ Ainv for r in little_rotations], dtype=np.complex128)
 
     rep = np.einsum("ipq,iab->ipaqb", perm_rep, rotation_rep, optimize="greedy")
     return rep.reshape(-1, num_atoms * 3, num_atoms * 3)

setup.py

@@ -23,16 +23,16 @@
 # What packages are optional?
 EXTRAS = {
     "dev": [
-        "pytest==8.2.2",
-        "pytest-cov==5.0.0",
+        "pytest",
+        "pytest-cov",
         "pre-commit",
         "black",
         "mypy",
         "flake8",
         "pyupgrade",
         "pydocstyle",
         "nbqa",
-        "phonopy==2.24.2",
+        "phonopy",
         # Jupyter notebook
         "notebook",
         "matplotlib",
@@ -93,6 +93,7 @@
         "Programming Language :: Python :: 3.8",
         "Programming Language :: Python :: 3.9",
         "Programming Language :: Python :: 3.10",
+        "Programming Language :: Python :: 3.11",
         "Topic :: Scientific/Engineering :: Mathematics",
         "Topic :: Scientific/Engineering :: Physics",
     ],

src/spgrep/corep.py

@@ -186,14 +186,14 @@ def enumerate_spinor_small_corepresentations(
         if indicator == 1:
             # Unitary matrix s.t. irrep @ U = conj_irrep @ U
             U = get_intertwiner(irrep, conj_irrep, atol, max_num_random_generations)
-            corep = np.zeros((order, dim, dim), dtype=np.complex_)
+            corep = np.zeros((order, dim, dim), dtype=np.complex128)
             corep[xsg_indices] = irrep
             corep[a0u] = (
                 np.conj(factor_system[a0_idx, xsg_indices])[:, None, None] * U[None, :, :] @ irrep
             )
         elif indicator == -1:
             U = get_intertwiner(irrep, conj_irrep, atol, max_num_random_generations)
-            corep = np.zeros((order, 2 * dim, 2 * dim), dtype=np.complex_)
+            corep = np.zeros((order, 2 * dim, 2 * dim), dtype=np.complex128)
 
             # [ [irrep, 0],
             #   [0, conj_irrep]]
@@ -202,7 +202,7 @@ def enumerate_spinor_small_corepresentations(
 
             # [ [0, -U],
             #   [U, 0] ]
-            corep_a0 = np.zeros((2 * dim, 2 * dim), dtype=np.complex_)
+            corep_a0 = np.zeros((2 * dim, 2 * dim), dtype=np.complex128)
             corep_a0[:dim, dim:] = -U
             corep_a0[dim:, :dim] = U
 
@@ -212,7 +212,7 @@ def enumerate_spinor_small_corepresentations(
                 @ corep[xsg_indices]
             )
         elif indicator == 0:
-            corep = np.zeros((order, 2 * dim, 2 * dim), dtype=np.complex_)
+            corep = np.zeros((order, 2 * dim, 2 * dim), dtype=np.complex128)
 
             # [ [irrep, 0],
             #   [0, conj_irrep]]
@@ -221,11 +221,11 @@ def enumerate_spinor_small_corepresentations(
 
             # [ [0, omega(a0, a0) irrep[a0 * a0]],
             #   [1, 0] ]
-            corep_a0 = np.zeros((2 * dim, 2 * dim), dtype=np.complex_)
+            corep_a0 = np.zeros((2 * dim, 2 * dim), dtype=np.complex128)
             corep_a0[:dim, dim:] = (
                 factor_system[a0_idx, a0_idx] * irrep[xsg_indices_mapping[table[a0_idx, a0_idx]]]
             )
-            corep_a0[dim:, :dim] = np.eye(dim, dtype=np.complex_)
+            corep_a0[dim:, :dim] = np.eye(dim, dtype=np.complex128)
 
             corep[a0u] = (
                 np.conj(factor_system[a0_idx, xsg_indices])[:, None, None]
@@ -277,12 +277,12 @@ def get_corep_spinor_factor_system(
             [0, -1],
             [1, 0],
         ],
-        dtype=np.complex_,
+        dtype=np.complex128,
     )
 
     # Assign a unitary or anti-unitary operator for each magnetic operation
     order = len(rotations)
-    unitary_rotations = np.zeros((order, 2, 2), dtype=np.complex_)
+    unitary_rotations = np.zeros((order, 2, 2), dtype=np.complex128)
     anti_linear = np.zeros((order,), dtype=np.bool_)
     for i, (rotation, time_reversal) in enumerate(zip(rotations, time_reversals)):
         unitary_rotation = get_spinor_unitary_rotation(lattice, rotation)
@@ -297,7 +297,7 @@ def get_corep_spinor_factor_system(
 
     # Factor system for spinor co-rep
     table = get_cayley_table(rotations, time_reversals)
-    corep_spinor_factor_system = np.zeros((order, order), dtype=np.complex_)
+    corep_spinor_factor_system = np.zeros((order, order), dtype=np.complex128)
     for i, (ui, ai) in enumerate(zip(unitary_rotations, anti_linear)):
         for j, uj in enumerate(unitary_rotations):
             if ai:

src/spgrep/group.py

@@ -131,7 +131,7 @@ def get_factor_system_from_little_group(
         # Never take modulus!
         residuals[i] = rotation.T @ kpoint - kpoint
 
-    factor_system = np.zeros((n, n), dtype=np.complex_)
+    factor_system = np.zeros((n, n), dtype=np.complex128)
     for i, residual in enumerate(residuals):
         for j, translation in enumerate(little_translations):
             factor_system[i, j] = np.exp(-2j * np.pi * np.dot(residual, translation))

src/spgrep/irreps.py

@@ -177,7 +177,7 @@ def enumerate_unitary_irreps(
     """
     order = rotations.shape[0]
     if factor_system is None:
-        factor_system = np.ones((order, order), dtype=np.complex_)
+        factor_system = np.ones((order, order), dtype=np.complex128)
 
     if method == "Neto":
         table = get_cayley_table(rotations)
@@ -264,10 +264,10 @@ def enumerate_unitary_irreps_from_regular_representation(
 
     # For (m, i), reg[m, i, :] has only one nonzero entry.
     # To reduce computational time, suppress reg to only nonzero elements
-    reg_nonzero = np.zeros((n, n), dtype=np.complex_)
+    reg_nonzero = np.zeros((n, n), dtype=np.complex128)
     lookup = np.zeros((n, n), dtype=int)
     for m, i in product(range(n), repeat=2):
-        idx = np.nonzero(reg[m, i, :])[0]
+        idx = np.nonzero(reg[m, i, :])[0].item()
         reg_nonzero[m, i] = reg[m, i, idx]
         lookup[m, i] = idx
 
@@ -277,7 +277,7 @@ def enumerate_unitary_irreps_from_regular_representation(
         hermite_random = rng.random((n, n)) + rng.random((n, n)) * 1j
         hermite_random += np.conj(hermite_random.T)
 
-        hermite_random_reordered = np.zeros((n, n, n), dtype=np.complex_)
+        hermite_random_reordered = np.zeros((n, n, n), dtype=np.complex128)
         meshi, meshj = np.meshgrid(range(n), range(n))
         # hermite_random_reordered[m, i, j] = hermite_random[lookup[m, i], lookup[m, j]]
         for m in range(n):
@@ -434,7 +434,7 @@ def enumerate_unitary_irreps_from_solvable_group_chain(
     """
     identity = get_identity_index(table)
     group = [identity]  # int -> GroupIdx
-    irreps = [np.ones((1, 1, 1), dtype=np.complex_)]
+    irreps = [np.ones((1, 1, 1), dtype=np.complex128)]
 
     # Extend subgroups from identity to whole
     for r in solvable_chain_generators[::-1]:
@@ -504,13 +504,13 @@ def enumerate_unitary_irreps_from_solvable_group_chain(
                     omegaq = omega * np.exp(2j * np.pi * q / p)
                     delta_r = intertwiner / omegaq  # Rep. matrix for r
                     delta_rm = [
-                        np.eye(intertwiner.shape[0], dtype=np.complex_)
+                        np.eye(intertwiner.shape[0], dtype=np.complex128)
                     ]  # delta_rm[m] is rep. matrix for r^m
                     for m in range(1, p):
                         # D(r^m) = D(r) @ D(r^{m-1}) / mu(r, r^{m-1})
                         delta_rm.append(delta_r @ delta_rm[m - 1] / factor_system[r, rm[m - 1]])
 
-                    next_irrep = np.zeros((len(group), dim, dim), dtype=np.complex_)
+                    next_irrep = np.zeros((len(group), dim, dim), dtype=np.complex128)
                     for m in range(p):
                         for s in subgroup:
                             idx = table[rm[m], s]
@@ -522,7 +522,7 @@ def enumerate_unitary_irreps_from_solvable_group_chain(
                     next_sub_irreps.append(next_irrep)
             else:
                 # Mutually inequivalent
-                next_irrep = np.zeros((len(group), dim * p, dim * p), dtype=np.complex_)
+                next_irrep = np.zeros((len(group), dim * p, dim * p), dtype=np.complex128)
                 for m in range(p):
                     for s in subgroup:
                         idx = table[rm[m], s]
@@ -613,7 +613,7 @@ def get_physically_irrep(
         real_irrep = np.real(np.einsum("il,klm,mj->kij", T, irrep, np.conj(T), optimize="greedy"))
 
     elif indicator in [-1, 0]:
-        real_irrep = np.empty((order, 2 * dim, 2 * dim), dtype=np.float_)
+        real_irrep = np.empty((order, 2 * dim, 2 * dim), dtype=np.float64)
         # [ [Re D(g),  Im D(g)]
         #   [-Im D(g), Re D(g)] ]
         real_irrep[:, :dim, :dim] = np.real(irrep)

src/spgrep/representation.py

@@ -59,7 +59,7 @@ def get_projective_regular_representation(
     n = len(rotations)
     table = get_cayley_table(rotations)
 
-    reg = np.zeros((n, n, n), dtype=np.complex_)
+    reg = np.zeros((n, n, n), dtype=np.complex128)
     for k, j in product(range(n), repeat=2):
         reg[k, table[k, j], j] = factor_system[k, j]
 
@@ -120,7 +120,7 @@ def get_character(representation: NDArrayComplex) -> NDArrayComplex:
     -------
     character: array, (order, )
     """
-    character = np.einsum("ijj->i", representation, optimize="greedy").astype(np.complex_)
+    character = np.einsum("ijj->i", representation, optimize="greedy").astype(np.complex128)
     return character
 
 
@@ -251,7 +251,7 @@ def is_representation(
     """Return true if given matrix function is a (projective) representation with given factor system."""
     order = rep.shape[0]
     if factor_system is None:
-        factor_system = np.ones((order, order), dtype=np.complex_)
+        factor_system = np.ones((order, order), dtype=np.complex128)
 
     for i, ri in enumerate(rep):
         for j, rj in enumerate(rep):
@@ -302,7 +302,7 @@ def check_spacegroup_representation(
     """Check definition of representation. This function works for primitive and conventional cell."""
     order = len(little_rotations)
     if spinor_factor_system is None:
-        spinor_factor_system = np.ones((order, order), dtype=np.complex_)
+        spinor_factor_system = np.ones((order, order), dtype=np.complex128)
 
     little_rotations_int = [ndarray2d_to_integer_tuple(rotation) for rotation in little_rotations]
 

src/spgrep/spinor.py

@@ -125,13 +125,13 @@ def get_spinor_factor_system(
     order = len(rotations)
 
     # Assign a SU(2) rotation for each O(3) operation
-    unitary_rotations = np.zeros((order, 2, 2), dtype=np.complex_)
+    unitary_rotations = np.zeros((order, 2, 2), dtype=np.complex128)
     for i, rotation in enumerate(rotations):
         unitary_rotations[i] = get_spinor_unitary_rotation(lattice, rotation)
 
     # Factor system from spin
     table = get_cayley_table(rotations)
-    spinor_factor_system = np.zeros((order, order), dtype=np.complex_)
+    spinor_factor_system = np.zeros((order, order), dtype=np.complex128)
     for (i, ui), (j, uj) in product(enumerate(unitary_rotations), repeat=2):
         # si @ sj = sk in O(3)
         k = table[i, j]

src/spgrep/utils.py

@@ -11,8 +11,8 @@
 
 NDArrayBool: TypeAlias = NDArray[np.bool_]
 NDArrayInt: TypeAlias = NDArray[np.int_]
-NDArrayFloat: TypeAlias = NDArray[np.float_]
-NDArrayComplex: TypeAlias = NDArray[np.complex_]
+NDArrayFloat: TypeAlias = NDArray[np.float64]
+NDArrayComplex: TypeAlias = NDArray[np.complex128]
 
 
 def is_integer_array(array: NDArrayFloat, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
@@ -54,7 +54,7 @@ def is_prime(n: int) -> bool:
     return True
 
 
-def nroot(z: np.complex_, n: int) -> np.complex_:
+def nroot(z: np.complex128, n: int) -> np.complex128:
     """Return `n`-th power root of `z` with the minimum angle."""
     root = z ** (1 / n)
     r = np.absolute(root)

tests/test_corep.py

@@ -57,7 +57,7 @@ def test_corep_spinor_factor_system(request, symmetry_and_lattice):
     for unitary_rotation in unitary_rotations:
         assert np.allclose(
             unitary_rotation @ np.conj(unitary_rotation).T,
-            np.eye(2, dtype=np.complex_),
+            np.eye(2, dtype=np.complex128),
         )
 
     # Cocycle condition

tests/test_pir.py

@@ -57,7 +57,7 @@ def check_spacegroup_pir(
                 m12 = phase * rep[idx]
             else:
                 half_dim = rep.shape[1] // 2
-                rep_additional = np.zeros((2 * half_dim, 2 * half_dim), dtype=np.complex_)
+                rep_additional = np.zeros((2 * half_dim, 2 * half_dim), dtype=np.complex128)
                 # [ [ cos(kt), -sin(kt) ],
                 #   [ sin(kt), cos(kt) ] ]
                 rep_additional[:half_dim, :half_dim] = np.cos(

tests/test_spinor.py

@@ -48,7 +48,7 @@ def test_spinor_factor_system_symmorphic(C3v, hexagonal_lattice):
     for unitary_rotation in unitary_rotations:
         assert np.allclose(
             unitary_rotation @ np.conj(unitary_rotation).T,
-            np.eye(2, dtype=np.complex_),
+            np.eye(2, dtype=np.complex128),
         )
 
     # Check factor system
@@ -129,7 +129,7 @@ def test_get_spacegroup_spinor_irreps_from_primitive_symmetry(kpoint, shape_expe
     for unitary_rotation in little_unitary_rotations:
         assert np.allclose(
             unitary_rotation @ np.conj(unitary_rotation).T,
-            np.eye(2, dtype=np.complex_),
+            np.eye(2, dtype=np.complex128),
         )
 
     # Check as representation
@@ -159,7 +159,7 @@ def test_get_crystallographic_pointgroup_spinor_irreps_from_symmetry(Oh):
     for unitary_rotation in unitary_rotations:
         assert np.allclose(
             unitary_rotation @ np.conj(unitary_rotation).T,
-            np.eye(2, dtype=np.complex_),
+            np.eye(2, dtype=np.complex128),
         )
 
     # Check as representation