Improve GPT related functions.

pyquda_utils/core.py

@@ -3,6 +3,8 @@
 import numpy
 
 from pyquda import (
+    initGPU,
+    initQUDA,
     init,
     getCoordFromRank,
     getRankFromCoord,

pyquda_utils/gpt.py

@@ -48,7 +48,7 @@ def LatticeGaugeGPT(lattice: List[g.lattice], gen_simd_width: int, gauge: Lattic
         assert latt_info.size == gauge.latt_info.size
         for index in range(latt_info.Nd):
             gpt_shape = [i for sl in zip(gpt_simd, gpt_latt) for i in sl]
-            lattice[index].mview()[0] = (
+            lattice[index].mview()[0][:] = (
                 gauge.lexico()
                 .astype(f"<c{2 * gpt_prec}")
                 .reshape(*gpt_shape, Nc, Nc)
@@ -80,7 +80,7 @@ def LatticePropagatorGPT(lattice: g.lattice, gen_simd_width: int, propagator: La
     else:
         assert latt_info.size == propagator.latt_info.size
         gpt_shape = [i for sl in zip(gpt_simd, gpt_latt) for i in sl]
-        lattice.mview()[0] = (
+        lattice.mview()[0][:] = (
             propagator.lexico()
             .astype(f"<c{2 * gpt_prec}")
             .reshape(*gpt_shape, Ns, Ns, Nc, Nc)

pyquda_utils/milc_rhmc_param.py

@@ -0,0 +1,121 @@
+from math import log
+from typing import Dict, Tuple
+
+from pyquda.field import X, Y, Z, T
+from pyquda.action.abstract import LoopParam, RationalParam
+
+
+def symanzik_1loop_hisq(u_0: float, n_flavor: int):
+    beta_0 = 1
+    beta_1 = -beta_0 / (20 * u_0 * u_0) * (1.00 - (0.6264 - 1.1746 * n_flavor) * log(u_0))
+    beta_2 = beta_0 / (u_0 * u_0) * (0.0433 - 0.0156 * n_flavor) * log(u_0)
+    return LoopParam(
+        path=[
+            [X, Y, -X, -Y],
+            [X, Z, -X, -Z],
+            [X, T, -X, -T],
+            [Y, Z, -Y, -Z],
+            [Y, T, -Y, -T],
+            [Z, T, -Z, -T],
+            [X, X, Y, -X, -X, -Y],
+            [X, X, Z, -X, -X, -Z],
+            [X, X, T, -X, -X, -T],
+            [Y, Y, X, -Y, -Y, -X],
+            [Y, Y, Z, -Y, -Y, -Z],
+            [Y, Y, T, -Y, -Y, -T],
+            [Z, Z, X, -Z, -Z, -X],
+            [Z, Z, Y, -Z, -Z, -Y],
+            [Z, Z, T, -Z, -Z, -T],
+            [T, T, X, -T, -T, -X],
+            [T, T, Y, -T, -T, -Y],
+            [T, T, Z, -T, -T, -Z],
+            [X, Y, Z, -X, -Y, -Z],
+            [X, Y, -Z, -X, -Y, Z],
+            [X, -Y, Z, -X, Y, -Z],
+            [X, -Y, -Z, -X, Y, Z],
+            [X, Y, T, -X, -Y, -T],
+            [X, Y, -T, -X, -Y, T],
+            [X, -Y, T, -X, Y, -T],
+            [X, -Y, -T, -X, Y, T],
+            [X, Z, T, -X, -Z, -T],
+            [X, Z, -T, -X, -Z, T],
+            [X, -Z, T, -X, Z, -T],
+            [X, -Z, -T, -X, Z, T],
+            [Y, Z, T, -Y, -Z, -T],
+            [Y, Z, -T, -Y, -Z, T],
+            [Y, -Z, T, -Y, Z, -T],
+            [Y, -Z, -T, -Y, Z, T],
+        ],
+        coeff=[beta_0 for _ in range(6)] + [beta_1 for _ in range(12)] + [beta_2 for _ in range(16)],
+    )
+
+
+def poly4_rhmc_param(filename: str):
+    rhmc_params: Dict[Tuple[Tuple[float, ...], Tuple[int, ...]], RationalParam] = {}
+    with open(filename, "r") as f:
+        lines = [line.strip().split(" ") for line in f.readlines() if not line.startswith("#") and not line == "\n"]
+    n_pseudo = int(lines[0][1])
+    index = 1
+    for _ in range(n_pseudo):
+        while not lines[index][0] == "y_MD":
+            index += 1
+        index += 3
+        order = int(lines[index][1])
+        index += 1
+        norm_molecular_dynamics = float(lines[index][1])
+        index += 1
+        residue_molecular_dynamics = [float(lines[index + i][1]) for i in range(order)]
+        index += order + 1
+        offset_molecular_dynamics = [float(lines[index + i][1]) for i in range(order)]
+        index += order
+
+        n_flavor = tuple(int(lines[index][1 + i]) for i in range(4) if lines[index][1 + i] != "0")
+        mass = tuple(float(lines[index + 2][1 + i]) for i in range(4) if lines[index][1 + i] != "0")
+        print(mass, n_flavor)
+        index += 3
+        order = int(lines[index][1])
+        index += 1
+        norm_pseudo_fermion = float(lines[index][1])
+        index += 1
+        residue_pseudo_fermion = [float(lines[index + i][1]) for i in range(order)]
+        index += order + 1
+        offset_pseudo_fermion = [float(lines[index + i][1]) for i in range(order)]
+        index += order
+
+        index += 3
+        order = int(lines[index][1])
+        index += 1
+        norm_fermion_action = float(lines[index][1])
+        index += 1
+        residue_fermion_action = [float(lines[index + i][1]) for i in range(order)]
+        index += order + 1
+        offset_fermion_action = [float(lines[index + i][1]) for i in range(order)]
+        index += order
+
+        key = (mass, n_flavor)
+        if key not in rhmc_params:
+            rhmc_params[key] = RationalParam(
+                norm_molecular_dynamics,
+                residue_molecular_dynamics,
+                offset_molecular_dynamics,
+                norm_pseudo_fermion,
+                residue_pseudo_fermion,
+                offset_pseudo_fermion,
+                norm_fermion_action,
+                residue_fermion_action,
+                offset_fermion_action,
+            )
+        else:
+            assert rhmc_params[key] == RationalParam(
+                norm_molecular_dynamics,
+                residue_molecular_dynamics,
+                offset_molecular_dynamics,
+                norm_pseudo_fermion,
+                residue_pseudo_fermion,
+                offset_pseudo_fermion,
+                norm_fermion_action,
+                residue_fermion_action,
+                offset_fermion_action,
+            )
+
+    return rhmc_params

tests/test.shift.py

@@ -1,18 +1,19 @@
 from check_pyquda import weak_field
 
 from pyquda_utils import core, io
+from pyquda_utils.core import X, Y, Z, T
 
 core.init([1, 1, 1, 1], [4, 4, 4, 8], 1, 1.0, resource_path=".cache")
 
 gauge = io.readQIOGauge(weak_field)
 gauge.toDevice()
-gauge_shift = gauge.shift([4, 5, 6, 7])
+gauge_shift = gauge.shift([-X, -Y, -Z, -T])
 
-print(gauge.data[0, 0, 0, 0, 0, 0])
-print(gauge_shift.data[0, 1, 0, 0, 0, 0])
-print(gauge.data[1, 0, 0, 0, 0, 0])
-print(gauge_shift.data[1, 1, 0, 0, 1, 0])
-print(gauge.data[2, 0, 0, 0, 0, 0])
-print(gauge_shift.data[2, 1, 0, 1, 0, 0])
-print(gauge.data[3, 0, 0, 0, 0, 0])
-print(gauge_shift.data[3, 1, 1, 0, 0, 0])
+print(gauge[X].data[0, 0, 0, 0, 0])
+print(gauge_shift[X].data[1, 0, 0, 0, 0])
+print(gauge[Y].data[0, 0, 0, 0, 0])
+print(gauge_shift[Y].data[1, 0, 0, 1, 0])
+print(gauge[Z].data[0, 0, 0, 0, 0])
+print(gauge_shift[Z].data[1, 0, 1, 0, 0])
+print(gauge[T].data[0, 0, 0, 0, 0])
+print(gauge_shift[T].data[1, 1, 0, 0, 0])