From 424107740754ed5d28c40043cae774f6f5807bee Mon Sep 17 00:00:00 2001 From: sy3394 Date: Mon, 20 Feb 2023 15:50:09 +0200 Subject: [PATCH 01/10] first commit --- lyncs_quda/clover_field.py | 6 +- lyncs_quda/gauge_field.py | 12 +--- lyncs_quda/lattice_field.py | 129 +++++++++++++++++++----------------- lyncs_quda/lib.py | 4 +- test/test_field.py | 4 +- 5 files changed, 79 insertions(+), 76 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index ce89a84..cb747fa 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -70,7 +70,7 @@ def __init__( self._fmunu.global_lattice, dofs=(idof,), dtype=prec, - device=self._fmunu.device, + device=self._fmunu.device_id, empty=True, ) self._clover = new(idof) @@ -205,7 +205,7 @@ def clover_field(self): if not self._direct: lib.computeClover(self.quda_field, self._fmunu.quda_field, self.coeff) self._direct = True - return self._clover.field + return self._clover @property def inverse_field(self): @@ -213,7 +213,7 @@ def inverse_field(self): self.clover_field lib.cloverInvert(self.quda_field, self.computeTrLog) self._inverse = True - return self._cloverInv.field + return self._cloverInv @property def trLog(self): diff --git a/lyncs_quda/gauge_field.py b/lyncs_quda/gauge_field.py index d965791..ad2d4b5 100644 --- a/lyncs_quda/gauge_field.py +++ b/lyncs_quda/gauge_field.py @@ -65,9 +65,8 @@ def momentum(lattice, **kwargs): class GaugeField(LatticeField): "Mimics the quda::GaugeField object" - @LatticeField.field.setter - def field(self, field): - LatticeField.field.fset(self, field) + def _check_field(self, field): + super()._check_field(field)#probably ok if self.reconstruct == "INVALID": raise TypeError(f"Unrecognized field dofs {self.dofs}") @@ -130,11 +129,6 @@ def copy(self, other=None, out=None, **kwargs): return super().copy(other, out, **kwargs) - def __array_finalize__(self, obj): - "Support for __array_finalize__ standard" - # need to reset QUDA object when meta data of its Python wrapper is changed - self._quda = None - def _prepare(self, field, **kwargs): field = super()._prepare(field, **kwargs) is_momentum = kwargs.get("is_momentum", self.is_momentum) @@ -501,7 +495,7 @@ def project(self, tol=None): if tol is None: tol = numpy.finfo(self.dtype).eps - assert self.device == cupy.cuda.runtime.getDevice() + assert self.device_id == cupy.cuda.runtime.getDevice() fails = cupy.zeros((1,), dtype="int32") lib.projectSU3(self.quda_field, tol, to_pointer(fails.data.ptr, "int *")) # return fails.get()[0] # shouldn't we reduce? diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 539e031..a48c9eb 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -106,7 +106,7 @@ def backend(device=True): cupy.cuda.runtime.setDevice(lib.device_id) -class LatticeField(numpy.lib.mixins.NDArrayOperatorsMixin): +class LatticeField: #(numpy.lib.mixins.NDArrayOperatorsMixin): "Mimics the quda::LatticeField object" @classmethod @@ -164,26 +164,20 @@ def new(self, empty=True, **kwargs): self.global_lattice, dofs=kwargs.pop("dofs", self.dofs), dtype=kwargs.pop("dtype", self.dtype), - device=kwargs.pop("device", self.device), + device=kwargs.pop("device", self.device_id), comm=kwargs.pop("comm", self.comm), empty=empty, **kwargs, ) - out.__array_finalize__(self) return out - def __array_finalize__(self, obj): - "Support for __array_finalize__ standard" - if self.comm is None: - self.comm = obj.comm - def copy(self, other=None, out=None, **kwargs): "Returns out, a copy+=kwargs, of other if given, else of self" # src: other if given; otherwise self # dst: out (created anew, if not explicitly given) # ASSUME: if out != None, out <=> other/self+=kwargs # - # IF other and out are both given, this behaves like a classmethod + # If other and out are both given, this behaves like a classmethod # where out&other are casted into type(self) # check=False => here any output is accepted @@ -210,8 +204,8 @@ def equivalent(self, other, **kwargs): dtype = kwargs.get("dtype", self.dtype) if other.dtype != dtype: return False - device = kwargs.get("device", self.device) - if other.device != device: + device_id = kwargs.get("device", self.device_id) + if other.device_id != device_id: return False dofs = kwargs.get("dofs", self.dofs) if other.dofs != dofs: @@ -230,7 +224,6 @@ def _prepare(self, field, copy=False, check=False, **kwargs): if copy: return self.copy(other=field, **kwargs) raise ValueError("The given field is not appropriate") - field.__array_finalize__(self) return field def prepare(self, fields, **kwargs): @@ -258,8 +251,34 @@ def prepare_in(self, fields, **kwargs): kwargs.setdefault("copy", True) return self.prepare(fields, **kwargs) + _children = {} + + def __new__(cls, field, **kwargs): + #? does this get called when creating a new instance via view casting, as it is __new__ of a parent of the below child + # not sure as the doc says: 'view casting and new-from-template, the equivalent of ndarray.__new__(MySubClass,... is called, at the C level' + # if not called, this is ok; + if not isinstance(field, (numpy.ndarray, cupy.ndarray)): + raise TypeError( + f"Supporting only numpy or cupy for field, got {type(field)}" + ) + mod = numpy if isinstance(field, numpy.ndarray) else cupy + parent = mod.ndarray #type(field) + child = cls._children.setdefault(parent, type(cls.__name__+"ext", (cls, parent), {})) #inspect. + obj = mod.asarray(field).view(type=child) + + return obj + + #field check should be performed + def __array_finalize__(self, obj): + "Support for __array_finalize__ standard" + # called after __new__ in all cases of instance creation + # self: newly created instance + # obj: input instance + self._check_field(obj) + if obj is None: return + self.__init__(obj, comm=getattr(obj, "comm", None)) + def __init__(self, field, comm=None, **kwargs): - self.field = field if comm is False: self.comm = comm elif comm is None: @@ -269,21 +288,10 @@ def __init__(self, field, comm=None, **kwargs): self._quda = None self.activate() - def activate(self): - "Activates the current field. To be called before using the object in quda" - "to make sure the communicator is set for MPI" - # if self.comm is None, but #ranks>1, this will mess thigns up - lib.set_comm(self.comm) - - @property - def field(self): - "The underlaying lattice field" - return self._field - - @field.setter - def field(self, field): - if isinstance(field, LatticeField): - field = field.field + def _check_field(self, field=None): + if field is None: + field = self + #? the following can be removed? if not isinstance(field, (numpy.ndarray, cupy.ndarray)): raise TypeError( f"Supporting only numpy or cupy for field, got {type(field)}" @@ -292,29 +300,33 @@ def field(self, field): raise TypeError("Field is stored on a different device than the quda lib") if len(field.shape) < 4: raise ValueError("A lattice field should not have shape smaller than 4") - self._field = field - def get(self): + def activate(self): + "Activates the current field. To be called before using the object in quda" + "to make sure the communicator is set for MPI" + # if self.comm is None, but #ranks>1, this will mess thigns up + lib.set_comm(self.comm) + + def get(self, *args, **kwargs): "Returns the field as numpy array" - return self.__array__() + return self.__array__(*args, **kwargs) def __array__(self, *args, **kwargs): - out = self.field if self.device is not None: - out = out.get() - return out.__array__(*args, **kwargs) + return super().get(*args, **kwargs) + return super().__array__(*args, **kwargs) def complex_view(self): "Returns a complex view of the field" if self.iscomplex: - return self.field - return self.field.view(get_complex_dtype(self.dtype)) + return self + return self.view(get_complex_dtype(self.dtype)) def float_view(self): "Returns a complex view of the field" if not self.iscomplex: - return self.field - return self.field.view(get_float_dtype(self.dtype)) + return self + return self.view(get_float_dtype(self.dtype)) def default_view(self): "Returns the default view of the field including reshaping" @@ -323,26 +335,25 @@ def default_view(self): @property def backend(self): "The backend of the field: cupy or numpy" - if isinstance(self.field, cupy.ndarray): + if isinstance(self, cupy.ndarray): return cupy return numpy @property def device(self): "Device id of the field (None if not on GPUs)" - if isinstance(self.field, cupy.ndarray): - return self.field.device.id + if isinstance(self, cupy.ndarray): + return super().device return None @property - def shape(self): - "Shape of the field" - return self.field.shape - + def device_id(self): + return getattr(self.device, "id", None) + @property def location(self): "Memory location of the field (CPU or CUDA)" - return "CPU" if isinstance(self.field, numpy.ndarray) else "CUDA" + return "CPU" if isinstance(self, numpy.ndarray) else "CUDA" @property def quda_location(self): @@ -381,20 +392,15 @@ def dofs(self): "Shape of the per-site degrees of freedom" return self.shape[: -self.ndims] - @property - def dtype(self): - "Field data type" - return self.field.dtype - @property def iscomplex(self): "Whether the field dtype is complex" - return self.backend.iscomplexobj(self.field) + return self.backend.iscomplexobj(self) @property def isreal(self): "Whether the field dtype is real" - return self.backend.isrealobj(self.field) + return self.backend.isrealobj(self) @property def precision(self): @@ -424,9 +430,9 @@ def pad(self): @property def ptr(self): "Memory pointer" - if isinstance(self.field, numpy.ndarray): - return self.field.__array_interface__["data"][0] - return self.field.data.ptr + if isinstance(self, numpy.ndarray): + return self.__array_interface__["data"][0] + return self.data.ptr @property def quda_params(self): @@ -473,9 +479,8 @@ def gpu_field(self): def __array_ufunc__(self, ufunc, method, *args, **kwargs): out = kwargs.get("out", (self,))[0] - prepare = ( - lambda arg: out.prepare_in(arg).field + lambda arg: out.prepare_in(arg).view(self.backend.ndarray) if isinstance(arg, (LatticeField, cupy.ndarray, numpy.ndarray)) else arg ) @@ -488,8 +493,10 @@ def __array_ufunc__(self, ufunc, method, *args, **kwargs): kwargs[key] = prepare(val) fnc = getattr(ufunc, method) - - return self.prepare_out(fnc(*args, **kwargs), check=False) + result = fnc(*args, **kwargs) + if not isinstance(result, self.backend.ndarray): + return result + return self.prepare_out(result, check=False) def __bool__(self): - return bool(self.field.all()) + return bool(self.all()) diff --git a/lyncs_quda/lib.py b/lyncs_quda/lib.py index e3b6265..c8da08f 100644 --- a/lyncs_quda/lib.py +++ b/lyncs_quda/lib.py @@ -18,7 +18,7 @@ from lyncs_cppyy.ll import addressof, to_pointer from lyncs_utils import static_property, lazy_import from . import __path__ -from .config import QUDA_MPI, QUDA_GITVERSION, QUDA_PRECISION, QUDA_RECONSTRUCT +from .config import QUDA_MPI, GITVERSION, QUDA_PRECISION, QUDA_RECONSTRUCT cupy = lazy_import("cupy") @@ -38,7 +38,7 @@ def __init__(self, *args, **kwargs): self._device_id = 0 self._comm = None if not self.tune_dir: - self.tune_dir = user_data_dir("quda", "lyncs") + "/" + QUDA_GITVERSION + self.tune_dir = user_data_dir("quda", "lyncs") + "/" + GITVERSION super().__init__(*args, **kwargs) @property diff --git a/test/test_field.py b/test/test_field.py index 3e0fd0b..0176077 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -2,7 +2,7 @@ from pytest import raises import numpy as np import cupy as cp - +import inspect shape = (4, 3, 4, 4, 4, 4) @@ -22,6 +22,7 @@ def test_precision(): def test_numpy(): field = LatticeField(np.zeros(shape)) + print(inspect.getmodule(field)) assert field.location == "CPU" assert field.dims == shape[-4:] assert field.dofs == shape[:-4] @@ -55,6 +56,7 @@ def test_cupy(): assert field.precision == "double" assert not field + print(type(field)) assert field + 0 == field assert field + 1 != field assert field - 0 == field From 4b70ed1658a29e79433efaaf37f901ca4e5c5d67 Mon Sep 17 00:00:00 2001 From: sy3394 Date: Fri, 24 Feb 2023 15:35:47 +0200 Subject: [PATCH 02/10] 2nd commit --- lyncs_quda/clover_field.py | 19 +++++++++++-- lyncs_quda/gauge_field.py | 35 +++++++++++++++++++---- lyncs_quda/lattice_field.py | 57 +++++++++++++++++++------------------ lyncs_quda/lib.py | 8 ++++-- test/test_dirac.py | 36 +++++++++++------------ test/test_field.py | 9 +++++- test/test_gauge_field.py | 12 +++++--- test/test_spinor.py | 8 +++--- 8 files changed, 119 insertions(+), 65 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index cb747fa..ac376fe 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -36,6 +36,21 @@ class CloverField(LatticeField): * so that sigma_mu,nu = i[g_mu, g_nu], F_mu,nu = (Q_mu,nu - Q_nu,mu)/8 (1/2 is missing from sigma_mu,nu) """ + def __new__(cls, gauge, **kwargs): + #TODO: get dofs and local dims from kwargs, instead of getting them + # from self.shape assuming that it has the form (dofs, local_dims) + if not isinstance(field, (numpy.ndarray, cupy.ndarray)): + raise TypeError( + f"Supporting only numpy or cupy for field, got {type(field)}" + ) + parent = type(field) + child = cls._children.setdefault(parent, type(cls.__name__+"ext", (cls, parent), {})) + obj = inspect.getmodule(parent).asarray(field).view(type=child) + + #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) + #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) + return obj + def __init__( self, fmunu, @@ -53,7 +68,7 @@ def __init__( if not isinstance(fmunu, GaugeField): fmunu = GaugeField(fmunu) self._fmunu = fmunu.compute_fmunu() - super().__init__(self._fmunu.field, comm=self._fmunu.comm) + super().__init__(self._fmunu, comm=self._fmunu.comm) # QUDA clover field inherently works with real's not with complex's (c.f., include/clover_field_order.h) @@ -116,7 +131,7 @@ def default_view(self): shape = (2,) # even-odd shape += (self.dofs[0] // N, -1, N) - return self.field.view().reshape(shape) + return self.float_view().reshape(shape) @property def twisted(self): diff --git a/lyncs_quda/gauge_field.py b/lyncs_quda/gauge_field.py index ad2d4b5..b99d626 100644 --- a/lyncs_quda/gauge_field.py +++ b/lyncs_quda/gauge_field.py @@ -16,7 +16,7 @@ from time import time from math import sqrt from collections import defaultdict -import numpy +import numpy, inspect from lyncs_cppyy import make_shared, lib as tmp, to_pointer, array_to_pointers from lyncs_utils import prod, isiterable from .lib import lib, cupy @@ -65,9 +65,26 @@ def momentum(lattice, **kwargs): class GaugeField(LatticeField): "Mimics the quda::GaugeField object" - def _check_field(self, field): - super()._check_field(field)#probably ok - if self.reconstruct == "INVALID": + l_children = {} + + def new__(cls, field, **kwargs): + if not isinstance(field, (numpy.ndarray, cupy.ndarray)): + raise TypeError( + f"Supporting only numpy or cupy for field, got {type(field)}" + ) + parent = type(field) + child = cls._children.setdefault(parent, type(cls.__name__+"ext",(cls, parent), {})) + print("new 0",parent,child.__bases__,child.__mro__) + obj = inspect.getmodule(parent).asarray(field).view(type=child) + print("new") + + return obj + + def _check_field(self, field=None): + super()._check_field(field) + if field is None or not isinstance(field, LatticeField): + field = self + if field.reconstruct == "INVALID": raise TypeError(f"Unrecognized field dofs {self.dofs}") def new(self, reconstruct=None, geometry=None, **kwargs): @@ -126,12 +143,15 @@ def copy(self, other=None, out=None, **kwargs): dst = self if out is None else out if src.is_momentum != dst.is_momentum: kwargs.update({"is_momentum": True}) + print("changed") + print("gauge cp") return super().copy(other, out, **kwargs) def _prepare(self, field, **kwargs): field = super()._prepare(field, **kwargs) is_momentum = kwargs.get("is_momentum", self.is_momentum) + print("gageu _prep") field.is_momentum = is_momentum return field @@ -237,9 +257,12 @@ def is_momentum(self): @is_momentum.setter def is_momentum(self, value): + print("momo setter",value,self._quda) self._is_momentum = value if self._quda is not None: - self._quda.link_type = self.quda_link_type + self._quda = None + #self._quda.link_type = self.quda_link_type + #print(self.quda_link_type, self._quda.LinkType()) @property def t_boundary(self): @@ -288,8 +311,10 @@ def quda_params(self): def quda_field(self): "Returns an instance of quda::(cpu/cuda)GaugeField for QUDA_(CPU/CUDA)_FIELD_LOCATION" self.activate() + print("gauge_quda",self._quda) if self._quda is None: self._quda = make_shared(lib.GaugeField.Create(self.quda_params)) + print("after gagueg",self._quda.LinkType()) return self._quda def is_native(self): diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index a48c9eb..03f2fdd 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -8,7 +8,7 @@ from array import array from contextlib import contextmanager -import numpy +import numpy, inspect from lyncs_cppyy import nullptr from lyncs_utils import prod from .enums import QudaPrecision @@ -106,7 +106,7 @@ def backend(device=True): cupy.cuda.runtime.setDevice(lib.device_id) -class LatticeField: #(numpy.lib.mixins.NDArrayOperatorsMixin): +class LatticeField(numpy.lib.mixins.NDArrayOperatorsMixin): "Mimics the quda::LatticeField object" @classmethod @@ -182,14 +182,13 @@ def copy(self, other=None, out=None, **kwargs): # check=False => here any output is accepted out = self.prepare_out(out, check=False, **kwargs) - if other is None: other = self # we prepare other without copying because we do copy here! other = out.prepare_in(other, copy=False, check=False, **kwargs) try: out.quda_field.copy(other.quda_field) - except: # NotImplementedError: #raised if self is LatticeField# at least, serial version calls exit(1) from qudaError, which is not catched by this + except: # NotImplementedError: #raised if self is LatticeField# at least, serial version calls exit(1) from qudaError, which is not caught by this # As last resort trying to copy elementwise out.default_view()[:] = other.default_view() @@ -254,28 +253,33 @@ def prepare_in(self, fields, **kwargs): _children = {} def __new__(cls, field, **kwargs): - #? does this get called when creating a new instance via view casting, as it is __new__ of a parent of the below child - # not sure as the doc says: 'view casting and new-from-template, the equivalent of ndarray.__new__(MySubClass,... is called, at the C level' - # if not called, this is ok; + #TODO: get dofs and local dims from kwargs, instead of getting them + # from self.shape assuming that it has the form (dofs, local_dims) if not isinstance(field, (numpy.ndarray, cupy.ndarray)): raise TypeError( f"Supporting only numpy or cupy for field, got {type(field)}" ) - mod = numpy if isinstance(field, numpy.ndarray) else cupy - parent = mod.ndarray #type(field) - child = cls._children.setdefault(parent, type(cls.__name__+"ext", (cls, parent), {})) #inspect. - obj = mod.asarray(field).view(type=child) - + parent = type(field) + child = cls._children.setdefault(parent, type(cls.__name__+"ext",(cls, parent), {"__array_ufunc__":cls.__array_ufunc__})) + obj = inspect.getmodule(parent).asarray(field).view(type=child) + #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) + #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) + return obj #field check should be performed - def __array_finalize__(self, obj): + def __array_finalize__(self, obj, **kwargs): "Support for __array_finalize__ standard" - # called after __new__ in all cases of instance creation + # Note: this is called when creating a temporary, possibly causing + # some issues. Apparently, this temporary is flattened and loose + # shape information # self: newly created instance # obj: input instance self._check_field(obj) - if obj is None: return + if obj is None: return # can be removed; we are not bypassing ndarray.__new__ + # This will require some care when we use attr for dims and dofs in _check_field + #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) + #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) self.__init__(obj, comm=getattr(obj, "comm", None)) def __init__(self, field, comm=None, **kwargs): @@ -291,16 +295,11 @@ def __init__(self, field, comm=None, **kwargs): def _check_field(self, field=None): if field is None: field = self - #? the following can be removed? - if not isinstance(field, (numpy.ndarray, cupy.ndarray)): - raise TypeError( - f"Supporting only numpy or cupy for field, got {type(field)}" - ) if isinstance(field, cupy.ndarray) and field.device.id != lib.device_id: raise TypeError("Field is stored on a different device than the quda lib") if len(field.shape) < 4: raise ValueError("A lattice field should not have shape smaller than 4") - + def activate(self): "Activates the current field. To be called before using the object in quda" "to make sure the communicator is set for MPI" @@ -319,14 +318,15 @@ def __array__(self, *args, **kwargs): def complex_view(self): "Returns a complex view of the field" if self.iscomplex: - return self - return self.view(get_complex_dtype(self.dtype)) + return self.view(type=self.backend.ndarray) + return self.view(get_complex_dtype(self.dtype), self.backend.ndarray) def float_view(self): "Returns a complex view of the field" + #don't need to upcast if we keep dofs and dims as attributes if not self.iscomplex: - return self - return self.view(get_float_dtype(self.dtype)) + return self.view(type=self.backend.ndarray) + return self.view(get_float_dtype(self.dtype), self.backend.ndarray) def default_view(self): "Returns the default view of the field including reshaping" @@ -368,7 +368,7 @@ def ndims(self): @property def dims(self): "Shape of the local lattice dimensions" - return self.shape[-self.ndims :] + return self.shape[-self.ndims :] #self._dims #self.shape[-self.ndims :] @property def local_lattice(self): @@ -390,7 +390,7 @@ def quda_dims(self): @property def dofs(self): "Shape of the per-site degrees of freedom" - return self.shape[: -self.ndims] + return self.shape[: -self.ndims] #self._dofs #self.shape[: -self.ndims] @property def iscomplex(self): @@ -480,7 +480,7 @@ def gpu_field(self): def __array_ufunc__(self, ufunc, method, *args, **kwargs): out = kwargs.get("out", (self,))[0] prepare = ( - lambda arg: out.prepare_in(arg).view(self.backend.ndarray) + lambda arg: out.prepare_in(arg).view(type=self.backend.ndarray) if isinstance(arg, (LatticeField, cupy.ndarray, numpy.ndarray)) else arg ) @@ -496,6 +496,7 @@ def __array_ufunc__(self, ufunc, method, *args, **kwargs): result = fnc(*args, **kwargs) if not isinstance(result, self.backend.ndarray): return result + print("ufunc",self.shape,ufunc,type(result),result.shape,flush=True) return self.prepare_out(result, check=False) def __bool__(self): diff --git a/lyncs_quda/lib.py b/lyncs_quda/lib.py index c8da08f..1c9aa8a 100644 --- a/lyncs_quda/lib.py +++ b/lyncs_quda/lib.py @@ -18,7 +18,7 @@ from lyncs_cppyy.ll import addressof, to_pointer from lyncs_utils import static_property, lazy_import from . import __path__ -from .config import QUDA_MPI, GITVERSION, QUDA_PRECISION, QUDA_RECONSTRUCT +from .config import QUDA_MPI, QUDA_GITVERSION, QUDA_PRECISION, QUDA_RECONSTRUCT cupy = lazy_import("cupy") @@ -38,7 +38,7 @@ def __init__(self, *args, **kwargs): self._device_id = 0 self._comm = None if not self.tune_dir: - self.tune_dir = user_data_dir("quda", "lyncs") + "/" + GITVERSION + self.tune_dir = user_data_dir("quda", "lyncs") + "/" + QUDA_GITVERSION super().__init__(*args, **kwargs) @property @@ -292,7 +292,9 @@ def __del__(self): @fixture(scope="session") def fixlib(): "A fixture to guarantee that in pytest lib is finalized at the end" - if not lib.initialized: + if QUDA_MPI and MPI.COMM_WORLD.Get_size() > 1: + pass + elif not lib.initialized: lib.init_quda() yield lib if lib.initialized: diff --git a/test/test_dirac.py b/test/test_dirac.py index af931b3..2b11050 100644 --- a/test/test_dirac.py +++ b/test/test_dirac.py @@ -65,25 +65,25 @@ def test_zero(lib, lattice, device, gamma, mu, dtype=None): sf.uniform() kappa = random() dirac = gf.Dirac(kappa=kappa) - assert (dirac.M(sf).field == sf.field).all() - assert (dirac.Mdag(sf).field == sf.field).all() - assert (dirac.MdagM(sf).field == sf.field).all() - assert (dirac.MMdag(sf).field == sf.field).all() + assert (dirac.M(sf) == sf).all() + assert (dirac.Mdag(sf) == sf).all() + assert (dirac.MdagM(sf) == sf).all() + assert (dirac.MMdag(sf) == sf).all() dirac = gf.Dirac(kappa=kappa, mu=mu) - sfmu = (2 * kappa * mu) * 1j * sf.gamma5().field - assert np.allclose(dirac.M(sf).field, sf.field + sfmu) - assert np.allclose(dirac.Mdag(sf).field, sf.field - sfmu) - assert np.allclose(dirac.MdagM(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field) - assert np.allclose(dirac.MMdag(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field) + sfmu = (2 * kappa * mu) * 1j * sf.gamma5() + assert np.allclose(dirac.M(sf), sf + sfmu) + assert np.allclose(dirac.Mdag(sf), sf - sfmu) + assert np.allclose(dirac.MdagM(sf), (1 + (2 * kappa * mu) ** 2) * sf) + assert np.allclose(dirac.MMdag(sf), (1 + (2 * kappa * mu) ** 2) * sf) csw = random() dirac = gf.Dirac(kappa=kappa, mu=mu, csw=csw) - sfmu = (2 * kappa * mu) * 1j * sf.gamma5().field - assert np.allclose(dirac.M(sf).field, sf.field + sfmu) - assert np.allclose(dirac.Mdag(sf).field, sf.field - sfmu) - assert np.allclose(dirac.MdagM(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field) - assert np.allclose(dirac.MMdag(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field) + sfmu = (2 * kappa * mu) * 1j * sf.gamma5() + assert np.allclose(dirac.M(sf), sf + sfmu) + assert np.allclose(dirac.Mdag(sf), sf - sfmu) + assert np.allclose(dirac.MdagM(sf), (1 + (2 * kappa * mu) ** 2) * sf) + assert np.allclose(dirac.MMdag(sf), (1 + (2 * kappa * mu) ** 2) * sf) # @dtype_loop # enables dtype #Double precision multigrid has not been enabled @@ -99,10 +99,10 @@ def test_coarse_zero(lib, lattice, device, dtype=None): sf.uniform() dirac = gf.Dirac(clover=gf2) - assert (dirac.M(sf).field == sf.field).all() - assert (dirac.Mdag(sf).field == sf.field).all() - assert (dirac.MdagM(sf).field == sf.field).all() - assert (dirac.MMdag(sf).field == sf.field).all() + assert (dirac.M(sf) == sf).all() + assert (dirac.Mdag(sf) == sf).all() + assert (dirac.MdagM(sf) == sf).all() + assert (dirac.MMdag(sf) == sf).all() # @dtype_loop # enables dtype diff --git a/test/test_field.py b/test/test_field.py index 0176077..08bb1db 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -38,6 +38,9 @@ def test_numpy(): assert isinstance(field + 0, type(field)) field2 = field.copy() + print("GGGGGGG") + assert field+field2==0 + print("GGGGGGGGGGG@") field += 0 assert field2 == field field -= 0 @@ -47,7 +50,6 @@ def test_numpy(): field /= 1 assert field2 == field - def test_cupy(): field = LatticeField(cp.zeros(shape)) assert field.location == "CUDA" @@ -67,6 +69,11 @@ def test_cupy(): field2 = field.copy() field += 0 + print("G") + assert cp.add(field,field2)==0 + print("K") + assert cp.add(field,0)==0 + print("KK") assert field2 == field field -= 0 assert field2 == field diff --git a/test/test_gauge_field.py b/test/test_gauge_field.py index 980f427..289a9dd 100644 --- a/test/test_gauge_field.py +++ b/test/test_gauge_field.py @@ -60,7 +60,9 @@ def test_zero(lib, lattice, device, dtype): assert gf.project() == 4 * np.prod(lattice) gf.zero() - + print("NEWWWWW") + assert gf+gf==0 + print("NEWWWWEWEWE") gf2 = gf.new() gf2.gaussian() assert gf.dot(gf2) == 0 @@ -77,7 +79,9 @@ def test_zero(lib, lattice, device, dtype): gf3 = momentum(lattice, dtype=dtype, device=device) gf3.zero() + print("ss!!!!!!!!!!!!!!!!!",type(gf),gf.is_momentum) assert gf + gf3 == 0 + print("SSSS",gf.is_momentum) assert gf3 + gf == 0 @@ -134,7 +138,7 @@ def test_random(lib, lattice, device, dtype): gf2 = gf.copy() assert gf == gf2 - assert isclose(gf.norm2(), (gf.field**2).sum(), rel_tol=1e-6) + assert isclose(gf.norm2(), (gf**2).sum(), rel_tol=1e-6) @dtype_loop # enables dtype @@ -147,7 +151,7 @@ def test_exponential(lib, lattice, device, dtype): gf.unity() mom.copy(out=gf) - assert np.allclose(gf.field, 0) + assert np.allclose(gf, 0) # gf.is_momentum = False assert gf == 0 @@ -161,7 +165,7 @@ def test_exponential(lib, lattice, device, dtype): mom.gaussian(epsilon=0) gf2 = mom.exponentiate() - assert np.allclose(gf.field, gf2.field) + assert np.allclose(gf, gf2) assert gf2 == gf gf.gaussian() diff --git a/test/test_spinor.py b/test/test_spinor.py index 2c7d3b1..9ad4a21 100644 --- a/test/test_spinor.py +++ b/test/test_spinor.py @@ -51,13 +51,13 @@ def test_params(lib, lattice, device, dtype): def test_init(lib, lattice, device, dtype): sf = spinor(lattice, dtype=dtype, device=device) sf.zero() - assert (sf.field == 0).all() + assert (sf == 0).all() assert sf.norm1() == 0 assert sf.norm2() == 0 sf.uniform() - assert np.isclose(sf.field.mean(), 0.5 + 0.5j, atol=0.1) + assert np.isclose(sf.mean(), 0.5 + 0.5j, atol=0.1) sf.gaussian() - assert np.isclose(sf.field.mean(), 0.0, atol=0.1) + assert np.isclose(sf.mean(), 0.0, atol=0.1) # @dtype_loop # enables dtype @@ -68,4 +68,4 @@ def test_gamma5(lib, lattice, device, gamma, dtype=None): sf = spinor(lattice, dtype=dtype, device=device, gamma_basis=gamma) sf.uniform() sf2 = sf.gamma5().apply_gamma5() - assert (sf.field == sf2.field).all() + assert (sf == sf2).all() From adf5f13244aaeb0aa8ce85770eaf78cdf7d08306 Mon Sep 17 00:00:00 2001 From: sy3394 Date: Fri, 3 Mar 2023 18:42:29 +0200 Subject: [PATCH 03/10] lib, field, gauge, and clover are working --- lyncs_quda/clover_field.py | 146 ++++++++++++++++++------------------ lyncs_quda/gauge_field.py | 38 +++------- lyncs_quda/lattice_field.py | 13 ++-- lyncs_quda/spinor_field.py | 1 + test/test_clover_field.py | 6 +- test/test_field.py | 8 -- test/test_gauge_field.py | 36 +++++++-- 7 files changed, 129 insertions(+), 119 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index ac376fe..5d0a07e 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -6,7 +6,7 @@ "CloverField", ] -import numpy +import numpy, inspect from cppyy.gbl.std import vector from lyncs_cppyy import make_shared, to_pointer @@ -34,26 +34,48 @@ class CloverField(LatticeField): * Only rho is mutable. To change other params, a new instance should be created * QUDA convention for clover field := 1+i ( kappa csw )/4 sigma_mu,nu F_mu,nu (<-sigma_mu,nu: spinor tensor) * so that sigma_mu,nu = i[g_mu, g_nu], F_mu,nu = (Q_mu,nu - Q_nu,mu)/8 (1/2 is missing from sigma_mu,nu) + * Apparently, an input to QUDA clover object, coeff = kappa*csw + * wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp) """ - def __new__(cls, gauge, **kwargs): + _children = {} + + def __new__(cls, fmunu, **kwargs): #TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) - if not isinstance(field, (numpy.ndarray, cupy.ndarray)): + if not isinstance(fmunu, (numpy.ndarray, cupy.ndarray)): raise TypeError( - f"Supporting only numpy or cupy for field, got {type(field)}" + f"Supporting only numpy or cupy for field, got {type(fmunu)}" ) - parent = type(field) + is_clover = False + print("clv new0",type(fmunu),fmunu.shape,isinstance(fmunu,GaugeField),isinstance(fmunu,LatticeField),isinstance(fmunu, (numpy.ndarray, cupy.ndarray))) + if not isinstance(fmunu, GaugeField): + if kwargs.get("is_clover", False): + is_clover = True + parent = type(fmunu) + field = fmunu + print("it is clover") + else: + fmunu = GaugeField(fmunu) + print("try",fmunu.shape) + + if not is_clover: # not copying from a clover-field array + print("clv new",type(fmunu)) + idof = int((fmunu.ncol * fmunu.ndims) ** 2 / 2) + prec = fmunu.dtype + field = fmunu.backend.empty((idof,) + fmunu.dims, dtype=prec) + parent = fmunu.backend.ndarray child = cls._children.setdefault(parent, type(cls.__name__+"ext", (cls, parent), {})) obj = inspect.getmodule(parent).asarray(field).view(type=child) #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) + return obj def __init__( self, - fmunu, + obj, coeff=0.0, twisted=False, mu2=0, @@ -61,70 +83,50 @@ def __init__( eps2=0, rho=0, computeTrLog=False, + **kwargs ): - # ASSUME: coeff = kappa*csw wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp) - # ? better to store fmunu.quda_field to _fmunu -> import gauge_tensor to be used in some methods - # ? better to put clover into self.field -> need walk-around to make copy() work - if not isinstance(fmunu, GaugeField): - fmunu = GaugeField(fmunu) - self._fmunu = fmunu.compute_fmunu() - super().__init__(self._fmunu, comm=self._fmunu.comm) - - # QUDA clover field inherently works with real's not with complex's (c.f., include/clover_field_order.h) - - idof = int((self._fmunu.ncol * self._fmunu.ndims) ** 2 / 2) - prec = self._fmunu.precision - self._direct = ( - False # Here, it is a flag to indicate whether the field has been computed - ) - self._inverse = ( - False # Here, it is a flag to indicate whether the field has been computed - ) - - new = lambda idof: LatticeField.create( - self._fmunu.global_lattice, - dofs=(idof,), - dtype=prec, - device=self._fmunu.device_id, - empty=True, - ) - self._clover = new(idof) - self._cloverInv = new(idof) - self.coeff = coeff - self._twisted = twisted - self._twist_flavor = tf - self._mu2 = mu2 - self._eps2 = eps2 - self._rho = rho - self.computeTrLog = computeTrLog - - # shape, dofs, dtype, iscomlex, isreal are overwriten to report their values for the clover field, instead of _fmunu - - @property - def shape(self): - "Shape of the clover field" - return self._clover.shape - - @property - def dofs(self): - "Shape of the per-site degrees of freedom" - return self._clover.dofs - - @property - def dtype(self): - "Clover field data type" - return self._clover.dtype - - @property - def iscomplex(self): - "Whether the clover field dtype is complex" - return self._clover.iscomplex - - @property - def isreal(self): - "Whether the clover field dtype is real" - return self._clover.isreal - + # WARNING: ndarray object is not supposed to be view-casted to CloverField object + # except in __new__, for which __init__ will be called subsequently, + # as the result won't come with meta info such as 'coeff' or 'mu2' + + super().__init__(obj, getattr(obj, "comm", None)) + print("clv init",type(obj)) + if isinstance(obj, GaugeField): + # explicit construction + # QUDA clover field inherently works with real's not with complex's (c.f., include/clover_field_order.h) + self._cloverInv = LatticeField.create( + self.global_lattice, + dofs=self.dofs, + dtype=self.dtype, + device=self.device_id, + empty=True, + ) + self._fmunu = obj.compute_fmunu() + self._direct = ( + False # Here, it is a flag to indicate whether the field has been computed + ) + self._inverse = ( + False # Here, it is a flag to indicate whether the field has been computed + ) + self.coeff = coeff + self._twisted = twisted + self._twist_flavor = tf + self._mu2 = mu2 + self._eps2 = eps2 + self._rho = rho + self.computeTrLog = computeTrLog + elif isinstance(obj, CloverField): + # upcasting to ndarray or new from template + self.__dict__.update(obj.__dict__) + elif isinstance(obj, self.backend.ndarray): + pass + else: + raise ValueError("The input is expected to be ndarray or LatticeField object") + + def _prepare(self, field, copy=False, check=False, **kwargs): + # When CloverField object prepares its input, the input is assumed to be of CloverField + return super()._prepare(field, copy=copy, check=check, is_clover=True, **kwargs) + # naming suggestion: native_view? default_* saved for dofs+lattice? def default_view(self): N = 1 if self.order == "FLOAT2" else 4 @@ -193,7 +195,7 @@ def quda_params(self): params = lib.CloverFieldParam() lib.copy_struct(params, super().quda_params) params.inverse = True - params.clover = to_pointer(self._clover.ptr) + params.clover = to_pointer(self.ptr) params.cloverInv = to_pointer(self._cloverInv.ptr) params.coeff = self.coeff params.twisted = self.twisted @@ -220,7 +222,7 @@ def clover_field(self): if not self._direct: lib.computeClover(self.quda_field, self._fmunu.quda_field, self.coeff) self._direct = True - return self._clover + return self.view(self.backend.ndarray) @property def inverse_field(self): @@ -228,7 +230,7 @@ def inverse_field(self): self.clover_field lib.cloverInvert(self.quda_field, self.computeTrLog) self._inverse = True - return self._cloverInv + return self._cloverInv.view(self.backend.ndarray) @property def trLog(self): diff --git a/lyncs_quda/gauge_field.py b/lyncs_quda/gauge_field.py index b99d626..d9b5052 100644 --- a/lyncs_quda/gauge_field.py +++ b/lyncs_quda/gauge_field.py @@ -65,27 +65,21 @@ def momentum(lattice, **kwargs): class GaugeField(LatticeField): "Mimics the quda::GaugeField object" - l_children = {} - - def new__(cls, field, **kwargs): - if not isinstance(field, (numpy.ndarray, cupy.ndarray)): - raise TypeError( - f"Supporting only numpy or cupy for field, got {type(field)}" - ) - parent = type(field) - child = cls._children.setdefault(parent, type(cls.__name__+"ext",(cls, parent), {})) - print("new 0",parent,child.__bases__,child.__mro__) - obj = inspect.getmodule(parent).asarray(field).view(type=child) - print("new") - - return obj - + _children = {} + def _check_field(self, field=None): + # NOTE: + # - For now, we store dofs and local lattice info in the shape of the array + # - The canonical form for GaugeField is [dofs, local_lattice] + # where len(dofs) == 1,2 and len(local_lattice) = 4 super()._check_field(field) - if field is None or not isinstance(field, LatticeField): + if field is None: field = self - if field.reconstruct == "INVALID": - raise TypeError(f"Unrecognized field dofs {self.dofs}") + dofs = field.shape[:-self.ndims] + pdofs = prod(dofs) if dofs[0] not in self._geometry_values[1] or dofs[0]==1 else prod(dofs[1:]) + pdofs *= 2**(self.iscomplex) + if not (pdofs in (12, 8, 10) or sqrt(pdofs / 2).is_integer()): + raise TypeError(f"Unrecognized field dofs {dofs}") def new(self, reconstruct=None, geometry=None, **kwargs): "Returns a new empty field based on the current" @@ -143,15 +137,12 @@ def copy(self, other=None, out=None, **kwargs): dst = self if out is None else out if src.is_momentum != dst.is_momentum: kwargs.update({"is_momentum": True}) - print("changed") - print("gauge cp") return super().copy(other, out, **kwargs) def _prepare(self, field, **kwargs): field = super()._prepare(field, **kwargs) is_momentum = kwargs.get("is_momentum", self.is_momentum) - print("gageu _prep") field.is_momentum = is_momentum return field @@ -257,12 +248,9 @@ def is_momentum(self): @is_momentum.setter def is_momentum(self, value): - print("momo setter",value,self._quda) self._is_momentum = value if self._quda is not None: self._quda = None - #self._quda.link_type = self.quda_link_type - #print(self.quda_link_type, self._quda.LinkType()) @property def t_boundary(self): @@ -311,10 +299,8 @@ def quda_params(self): def quda_field(self): "Returns an instance of quda::(cpu/cuda)GaugeField for QUDA_(CPU/CUDA)_FIELD_LOCATION" self.activate() - print("gauge_quda",self._quda) if self._quda is None: self._quda = make_shared(lib.GaugeField.Create(self.quda_params)) - print("after gagueg",self._quda.LinkType()) return self._quda def is_native(self): diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 03f2fdd..0425fe3 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -179,7 +179,11 @@ def copy(self, other=None, out=None, **kwargs): # # If other and out are both given, this behaves like a classmethod # where out&other are casted into type(self) - + # TODO: + # * check if this dose not cause any bugs if it overwrites ndarray.copy + # - For now, we cast self to ndarray before performing ndarray methods like flatten + # - the second arg should be "order='C'" to match the signiture? + # check=False => here any output is accepted out = self.prepare_out(out, check=False, **kwargs) if other is None: @@ -189,7 +193,7 @@ def copy(self, other=None, out=None, **kwargs): try: out.quda_field.copy(other.quda_field) except: # NotImplementedError: #raised if self is LatticeField# at least, serial version calls exit(1) from qudaError, which is not caught by this - # As last resort trying to copy elementwise + # As a last resort, trying to copy elementwise out.default_view()[:] = other.default_view() return out @@ -218,7 +222,7 @@ def _prepare(self, field, copy=False, check=False, **kwargs): return self.new(**kwargs) cls = type(self) if not isinstance(field, cls): - field = cls(field) + field = cls(field, **kwargs) if check and not self.equivalent(field, **kwargs): if copy: return self.copy(other=field, **kwargs) @@ -268,7 +272,7 @@ def __new__(cls, field, **kwargs): return obj #field check should be performed - def __array_finalize__(self, obj, **kwargs): + def __array_finalize__(self, obj): "Support for __array_finalize__ standard" # Note: this is called when creating a temporary, possibly causing # some issues. Apparently, this temporary is flattened and loose @@ -496,7 +500,6 @@ def __array_ufunc__(self, ufunc, method, *args, **kwargs): result = fnc(*args, **kwargs) if not isinstance(result, self.backend.ndarray): return result - print("ufunc",self.shape,ufunc,type(result),result.shape,flush=True) return self.prepare_out(result, check=False) def __bool__(self): diff --git a/lyncs_quda/spinor_field.py b/lyncs_quda/spinor_field.py index 0187053..95842ec 100644 --- a/lyncs_quda/spinor_field.py +++ b/lyncs_quda/spinor_field.py @@ -37,6 +37,7 @@ def spinor_coarse(lattice, dofs=24, **kwargs): class SpinorField(LatticeField): "Mimics the quda::ColorSpinorField object" + _children = {} gammas = ["DEGRAND_ROSSI", "UKQCD", "CHIRAL"] @classmethod diff --git a/test/test_clover_field.py b/test/test_clover_field.py index 967f80c..084b63b 100644 --- a/test/test_clover_field.py +++ b/test/test_clover_field.py @@ -31,7 +31,6 @@ def test_default(lattice): @lattice_loop # enables lattice def test_params(lib, lattice, device, dtype): gf = gauge(lattice, dtype=dtype, device=device) - # gf.zero() clv = CloverField(gf, computeTrLog=True, coeff=0) params = clv.quda_params @@ -117,7 +116,10 @@ def test_unit(lib, lattice, device, dtype): gf, coeff=1.0, tf="SINGLET", twisted=True, mu2=mu2, computeTrLog=True ) - assert (clv.field == 0).all() + clv.fill(0) + assert clv == 0 + print("GGGGGGGGGGGG") + assert clv + 0 == 0 idof = int(((clv.ncol * clv.nspin) ** 2 / 2)) if dtype == "float64": tmp = np.zeros((idof,) + lattice, dtype=dtype).reshape((2, 2, 36, -1)) diff --git a/test/test_field.py b/test/test_field.py index 08bb1db..b664ba8 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -38,9 +38,6 @@ def test_numpy(): assert isinstance(field + 0, type(field)) field2 = field.copy() - print("GGGGGGG") - assert field+field2==0 - print("GGGGGGGGGGG@") field += 0 assert field2 == field field -= 0 @@ -69,11 +66,6 @@ def test_cupy(): field2 = field.copy() field += 0 - print("G") - assert cp.add(field,field2)==0 - print("K") - assert cp.add(field,0)==0 - print("KK") assert field2 == field field -= 0 assert field2 == field diff --git a/test/test_gauge_field.py b/test/test_gauge_field.py index 289a9dd..88ebd3f 100644 --- a/test/test_gauge_field.py +++ b/test/test_gauge_field.py @@ -10,7 +10,7 @@ epsilon_loop, ) from lyncs_cppyy.ll import addressof -from lyncs_utils import isclose, allclose +from lyncs_utils import isclose#, allclose @lattice_loop @@ -60,9 +60,6 @@ def test_zero(lib, lattice, device, dtype): assert gf.project() == 4 * np.prod(lattice) gf.zero() - print("NEWWWWW") - assert gf+gf==0 - print("NEWWWWEWEWE") gf2 = gf.new() gf2.gaussian() assert gf.dot(gf2) == 0 @@ -79,9 +76,7 @@ def test_zero(lib, lattice, device, dtype): gf3 = momentum(lattice, dtype=dtype, device=device) gf3.zero() - print("ss!!!!!!!!!!!!!!!!!",type(gf),gf.is_momentum) assert gf + gf3 == 0 - print("SSSS",gf.is_momentum) assert gf3 + gf == 0 @@ -239,6 +234,35 @@ def test_force(lib, lattice, device, epsilon): zeros = getattr(gf, path + "_field")(coeffs=0, force=True) assert zeros == 0 +from lyncs_utils import isiterable +from collections.abc import Mapping +def values(dct): + "Calls values, if available, or dict.values" + try: + return dct.values() + except AttributeError: + return dict.values(dct) +def allclose(left, right, **kwargs): + if isinstance(left, cp.ndarray) and isinstance(right, cp.ndarray): + return np.allclose(left,right) + if isinstance(left, cp.ndarray) and not isinstance(right, cp.ndarray): + left = [left] * len(right) + if not isinstance(left, cp.ndarray) and isinstance(right, cp.ndarray): + right = [right] * len(left) + if len(left) != len(right): + return False + if isinstance(left, Mapping) or isinstance(right, Mapping): + if not isinstance(left, Mapping): + pairs = zip(left, values(right)) + elif not isinstance(right, Mapping): + pairs = zip(values(left), right) + else: + if set(keys(left)) != set(keys(right)): + return False + pairs = dictzip(left, right, values_only=True) + else: + pairs = zip(left, right) + return all((allclose(*pair, **kwargs) for pair in pairs)) # @dtype_loop # enables dtype @device_loop # enables device From 29f835190af06c24fc67203ab2329b718a95534f Mon Sep 17 00:00:00 2001 From: sy3394 Date: Fri, 3 Mar 2023 19:01:03 +0200 Subject: [PATCH 04/10] some cleanup; working version execept for evenodd.py due to segfault --- lyncs_quda/clover_field.py | 5 ----- lyncs_quda/testing.py | 1 + test/test_clover_field.py | 1 - test/test_dirac.py | 3 +++ test/test_solver.py | 2 +- 5 files changed, 5 insertions(+), 7 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index 5d0a07e..b0e4dc6 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -48,19 +48,15 @@ def __new__(cls, fmunu, **kwargs): f"Supporting only numpy or cupy for field, got {type(fmunu)}" ) is_clover = False - print("clv new0",type(fmunu),fmunu.shape,isinstance(fmunu,GaugeField),isinstance(fmunu,LatticeField),isinstance(fmunu, (numpy.ndarray, cupy.ndarray))) if not isinstance(fmunu, GaugeField): if kwargs.get("is_clover", False): is_clover = True parent = type(fmunu) field = fmunu - print("it is clover") else: fmunu = GaugeField(fmunu) - print("try",fmunu.shape) if not is_clover: # not copying from a clover-field array - print("clv new",type(fmunu)) idof = int((fmunu.ncol * fmunu.ndims) ** 2 / 2) prec = fmunu.dtype field = fmunu.backend.empty((idof,) + fmunu.dims, dtype=prec) @@ -90,7 +86,6 @@ def __init__( # as the result won't come with meta info such as 'coeff' or 'mu2' super().__init__(obj, getattr(obj, "comm", None)) - print("clv init",type(obj)) if isinstance(obj, GaugeField): # explicit construction # QUDA clover field inherently works with real's not with complex's (c.f., include/clover_field_order.h) diff --git a/lyncs_quda/testing.py b/lyncs_quda/testing.py index 99a25b4..df89fec 100644 --- a/lyncs_quda/testing.py +++ b/lyncs_quda/testing.py @@ -9,6 +9,7 @@ "parallel_loop", "mtype_loop", "mark_mpi", + "epsilon_loop", ] from random import random diff --git a/test/test_clover_field.py b/test/test_clover_field.py index 084b63b..6019237 100644 --- a/test/test_clover_field.py +++ b/test/test_clover_field.py @@ -118,7 +118,6 @@ def test_unit(lib, lattice, device, dtype): clv.fill(0) assert clv == 0 - print("GGGGGGGGGGGG") assert clv + 0 == 0 idof = int(((clv.ncol * clv.nspin) ** 2 / 2)) if dtype == "float64": diff --git a/test/test_dirac.py b/test/test_dirac.py index 2b11050..6917b01 100644 --- a/test/test_dirac.py +++ b/test/test_dirac.py @@ -1,7 +1,9 @@ from random import random import numpy as np +from lyncs_utils import isclose from lyncs_quda import ( gauge, + momentum, spinor, gauge_coarse, gauge_scalar, @@ -15,6 +17,7 @@ dtype_loop, mu_loop, gamma_loop, + epsilon_loop, ) diff --git a/test/test_solver.py b/test/test_solver.py index 2fc4137..2aa01c3 100644 --- a/test/test_solver.py +++ b/test/test_solver.py @@ -20,6 +20,6 @@ def test_solve_random(lib, lattice, device, gamma, dtype=None): mat = dirac.M out = mat.solve(rhs, delta=1e-4) # this value allowed convergence for all cases res = mat(out) - res.field -= rhs.field + res -= rhs res = res.norm() / rhs.norm() assert res < 1e-9 From 1c54f6706f16a3c44f867f842f977b829ecf3d8d Mon Sep 17 00:00:00 2001 From: sy3394 Date: Tue, 7 Mar 2023 19:51:11 +0200 Subject: [PATCH 05/10] minor change for better performance --- lyncs_quda/clover_field.py | 14 +++++++++++--- lyncs_quda/dirac.py | 4 ++-- lyncs_quda/lattice_field.py | 15 +++++++++++---- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index b0e4dc6..3252704 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -43,6 +43,8 @@ class CloverField(LatticeField): def __new__(cls, fmunu, **kwargs): #TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) + if isinstance(fmunu, CloverField): + return fmunu if not isinstance(fmunu, (numpy.ndarray, cupy.ndarray)): raise TypeError( f"Supporting only numpy or cupy for field, got {type(fmunu)}" @@ -55,14 +57,20 @@ def __new__(cls, fmunu, **kwargs): field = fmunu else: fmunu = GaugeField(fmunu) - + + # Set parent if not is_clover: # not copying from a clover-field array idof = int((fmunu.ncol * fmunu.ndims) ** 2 / 2) prec = fmunu.dtype field = fmunu.backend.empty((idof,) + fmunu.dims, dtype=prec) parent = fmunu.backend.ndarray - child = cls._children.setdefault(parent, type(cls.__name__+"ext", (cls, parent), {})) - obj = inspect.getmodule(parent).asarray(field).view(type=child) + # Set child + if parent in cls._children.keys(): + child = cls._children.get(parent) + else: + child = type(cls.__name__+"ext",(cls, parent), {}) + cls._children.update({parent: child}) + obj = field.view(type=child) #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) diff --git a/lyncs_quda/dirac.py b/lyncs_quda/dirac.py index f35b02a..5de9d98 100644 --- a/lyncs_quda/dirac.py +++ b/lyncs_quda/dirac.py @@ -242,9 +242,9 @@ def action(self, phi, **params): out = 0 if not self.full and "CLOVER" in self.type: if self.even: - out -= 2 * self.clover.trLog[1] - else: out -= 2 * self.clover.trLog[0] + else: + out -= 2 * self.clover.trLog[1] parity = None if not self.full: diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 0425fe3..97d9b99 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -8,7 +8,7 @@ from array import array from contextlib import contextmanager -import numpy, inspect +import numpy from lyncs_cppyy import nullptr from lyncs_utils import prod from .enums import QudaPrecision @@ -259,13 +259,20 @@ def prepare_in(self, fields, **kwargs): def __new__(cls, field, **kwargs): #TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) + + if isinstance(field, LatticeField): + return field if not isinstance(field, (numpy.ndarray, cupy.ndarray)): raise TypeError( f"Supporting only numpy or cupy for field, got {type(field)}" ) - parent = type(field) - child = cls._children.setdefault(parent, type(cls.__name__+"ext",(cls, parent), {"__array_ufunc__":cls.__array_ufunc__})) - obj = inspect.getmodule(parent).asarray(field).view(type=child) + parent = type(field) + if parent in cls._children.keys(): + child = cls._children.get(parent) + else: + child = type(cls.__name__+"ext",(cls, parent), {}) + cls._children.update({parent: child}) + obj = field.view(type=child) #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) From 5f10bad93a2617b26a5d737384a449aa2ff494c6 Mon Sep 17 00:00:00 2001 From: sy3394 Date: Wed, 8 Mar 2023 16:48:19 +0200 Subject: [PATCH 06/10] applied suggestions --- lyncs_quda/clover_field.py | 2 +- lyncs_quda/gauge_field.py | 4 +--- lyncs_quda/lattice_field.py | 15 ++++++--------- lyncs_quda/spinor_field.py | 1 - test/test_field.py | 5 ++--- 5 files changed, 10 insertions(+), 17 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index 3252704..5d0f2f6 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -6,7 +6,7 @@ "CloverField", ] -import numpy, inspect +import numpy from cppyy.gbl.std import vector from lyncs_cppyy import make_shared, to_pointer diff --git a/lyncs_quda/gauge_field.py b/lyncs_quda/gauge_field.py index d9b5052..3ac8d99 100644 --- a/lyncs_quda/gauge_field.py +++ b/lyncs_quda/gauge_field.py @@ -16,7 +16,7 @@ from time import time from math import sqrt from collections import defaultdict -import numpy, inspect +import numpy from lyncs_cppyy import make_shared, lib as tmp, to_pointer, array_to_pointers from lyncs_utils import prod, isiterable from .lib import lib, cupy @@ -65,8 +65,6 @@ def momentum(lattice, **kwargs): class GaugeField(LatticeField): "Mimics the quda::GaugeField object" - _children = {} - def _check_field(self, field=None): # NOTE: # - For now, we store dofs and local lattice info in the shape of the array diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 97d9b99..372197f 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -260,19 +260,16 @@ def __new__(cls, field, **kwargs): #TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) - if isinstance(field, LatticeField): + if isinstance(field, cls): return field if not isinstance(field, (numpy.ndarray, cupy.ndarray)): raise TypeError( f"Supporting only numpy or cupy for field, got {type(field)}" ) parent = type(field) - if parent in cls._children.keys(): - child = cls._children.get(parent) - else: - child = type(cls.__name__+"ext",(cls, parent), {}) - cls._children.update({parent: child}) - obj = field.view(type=child) + if (cls, parent) not in cls._children: + cls._children[(cls,parent)] = type(cls.__name__+"ext",(cls, parent), {}) + obj = field.view(type=cls._children[(cls,parent)]) #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) @@ -379,7 +376,7 @@ def ndims(self): @property def dims(self): "Shape of the local lattice dimensions" - return self.shape[-self.ndims :] #self._dims #self.shape[-self.ndims :] + return self.shape[-self.ndims :] @property def local_lattice(self): @@ -401,7 +398,7 @@ def quda_dims(self): @property def dofs(self): "Shape of the per-site degrees of freedom" - return self.shape[: -self.ndims] #self._dofs #self.shape[: -self.ndims] + return self.shape[: -self.ndims] @property def iscomplex(self): diff --git a/lyncs_quda/spinor_field.py b/lyncs_quda/spinor_field.py index 95842ec..0187053 100644 --- a/lyncs_quda/spinor_field.py +++ b/lyncs_quda/spinor_field.py @@ -37,7 +37,6 @@ def spinor_coarse(lattice, dofs=24, **kwargs): class SpinorField(LatticeField): "Mimics the quda::ColorSpinorField object" - _children = {} gammas = ["DEGRAND_ROSSI", "UKQCD", "CHIRAL"] @classmethod diff --git a/test/test_field.py b/test/test_field.py index b664ba8..868eea9 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -2,7 +2,6 @@ from pytest import raises import numpy as np import cupy as cp -import inspect shape = (4, 3, 4, 4, 4, 4) @@ -22,7 +21,7 @@ def test_precision(): def test_numpy(): field = LatticeField(np.zeros(shape)) - print(inspect.getmodule(field)) + assert field.location == "CPU" assert field.dims == shape[-4:] assert field.dofs == shape[:-4] @@ -47,6 +46,7 @@ def test_numpy(): field /= 1 assert field2 == field + def test_cupy(): field = LatticeField(cp.zeros(shape)) assert field.location == "CUDA" @@ -55,7 +55,6 @@ def test_cupy(): assert field.precision == "double" assert not field - print(type(field)) assert field + 0 == field assert field + 1 != field assert field - 0 == field From 8220802cf3e833ce693e72aa8d7370242a7adb6a Mon Sep 17 00:00:00 2001 From: sy3394 Date: Wed, 8 Mar 2023 17:34:33 +0200 Subject: [PATCH 07/10] reduced redundancy in clover_field.py --- lyncs_quda/clover_field.py | 19 +------------------ lyncs_quda/lattice_field.py | 2 +- 2 files changed, 2 insertions(+), 19 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index 5d0f2f6..4dfe3fa 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -28,8 +28,6 @@ class CloverField(LatticeField): This is designed as an intermediary to QUDA CloverField class so that it should have 1-to-1 correspondence to an QUDA instance. Note: - * This class stores the corresponding gauge field in its "field" attribute - to make except-clause of copy() work * direct & inverse fields are both allocated upon initialization * Only rho is mutable. To change other params, a new instance should be created * QUDA convention for clover field := 1+i ( kappa csw )/4 sigma_mu,nu F_mu,nu (<-sigma_mu,nu: spinor tensor) @@ -38,8 +36,6 @@ class CloverField(LatticeField): * wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp) """ - _children = {} - def __new__(cls, fmunu, **kwargs): #TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) @@ -53,29 +49,16 @@ def __new__(cls, fmunu, **kwargs): if not isinstance(fmunu, GaugeField): if kwargs.get("is_clover", False): is_clover = True - parent = type(fmunu) field = fmunu else: fmunu = GaugeField(fmunu) - # Set parent if not is_clover: # not copying from a clover-field array idof = int((fmunu.ncol * fmunu.ndims) ** 2 / 2) prec = fmunu.dtype field = fmunu.backend.empty((idof,) + fmunu.dims, dtype=prec) - parent = fmunu.backend.ndarray - # Set child - if parent in cls._children.keys(): - child = cls._children.get(parent) - else: - child = type(cls.__name__+"ext",(cls, parent), {}) - cls._children.update({parent: child}) - obj = field.view(type=child) - - #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) - #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) - return obj + return super().__new__(cls, field, **kwargs) def __init__( self, diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 0208011..11e205d 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -267,7 +267,7 @@ def __new__(cls, field, **kwargs): raise TypeError( f"Supporting only numpy or cupy for field, got {type(field)}" ) - parent = type(field) + parent = numpy.ndarray if isinstance(field, numpy.ndarray) else cupy.ndarray if (cls, parent) not in cls._children: cls._children[(cls,parent)] = type(cls.__name__+"ext",(cls, parent), {}) obj = field.view(type=cls._children[(cls,parent)]) From 730730633a693f39283c32a5efaedb875c1cd282 Mon Sep 17 00:00:00 2001 From: Simone Bacchio Date: Wed, 8 Mar 2023 15:35:33 +0000 Subject: [PATCH 08/10] Applying black formatting (from Github Action) --- lyncs_quda/clover_field.py | 30 ++++++++++++++---------------- lyncs_quda/gauge_field.py | 10 +++++++--- lyncs_quda/lattice_field.py | 31 ++++++++++++++++--------------- test/test_field.py | 3 ++- test/test_gauge_field.py | 10 ++++++++-- 5 files changed, 47 insertions(+), 37 deletions(-) diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py index 4dfe3fa..ae4838b 100644 --- a/lyncs_quda/clover_field.py +++ b/lyncs_quda/clover_field.py @@ -32,12 +32,12 @@ class CloverField(LatticeField): * Only rho is mutable. To change other params, a new instance should be created * QUDA convention for clover field := 1+i ( kappa csw )/4 sigma_mu,nu F_mu,nu (<-sigma_mu,nu: spinor tensor) * so that sigma_mu,nu = i[g_mu, g_nu], F_mu,nu = (Q_mu,nu - Q_nu,mu)/8 (1/2 is missing from sigma_mu,nu) - * Apparently, an input to QUDA clover object, coeff = kappa*csw - * wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp) + * Apparently, an input to QUDA clover object, coeff = kappa*csw + * wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp) """ def __new__(cls, fmunu, **kwargs): - #TODO: get dofs and local dims from kwargs, instead of getting them + # TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) if isinstance(fmunu, CloverField): return fmunu @@ -52,14 +52,14 @@ def __new__(cls, fmunu, **kwargs): field = fmunu else: fmunu = GaugeField(fmunu) - - if not is_clover: # not copying from a clover-field array + + if not is_clover: # not copying from a clover-field array idof = int((fmunu.ncol * fmunu.ndims) ** 2 / 2) prec = fmunu.dtype field = fmunu.backend.empty((idof,) + fmunu.dims, dtype=prec) - + return super().__new__(cls, field, **kwargs) - + def __init__( self, obj, @@ -70,7 +70,7 @@ def __init__( eps2=0, rho=0, computeTrLog=False, - **kwargs + **kwargs, ): # WARNING: ndarray object is not supposed to be view-casted to CloverField object # except in __new__, for which __init__ will be called subsequently, @@ -88,12 +88,8 @@ def __init__( empty=True, ) self._fmunu = obj.compute_fmunu() - self._direct = ( - False # Here, it is a flag to indicate whether the field has been computed - ) - self._inverse = ( - False # Here, it is a flag to indicate whether the field has been computed - ) + self._direct = False # Here, it is a flag to indicate whether the field has been computed + self._inverse = False # Here, it is a flag to indicate whether the field has been computed self.coeff = coeff self._twisted = twisted self._twist_flavor = tf @@ -107,12 +103,14 @@ def __init__( elif isinstance(obj, self.backend.ndarray): pass else: - raise ValueError("The input is expected to be ndarray or LatticeField object") + raise ValueError( + "The input is expected to be ndarray or LatticeField object" + ) def _prepare(self, field, copy=False, check=False, **kwargs): # When CloverField object prepares its input, the input is assumed to be of CloverField return super()._prepare(field, copy=copy, check=check, is_clover=True, **kwargs) - + # naming suggestion: native_view? default_* saved for dofs+lattice? def default_view(self): N = 1 if self.order == "FLOAT2" else 4 diff --git a/lyncs_quda/gauge_field.py b/lyncs_quda/gauge_field.py index 20a4934..25f629a 100644 --- a/lyncs_quda/gauge_field.py +++ b/lyncs_quda/gauge_field.py @@ -82,9 +82,13 @@ def _check_field(self, field=None): super()._check_field(field) if field is None: field = self - dofs = field.shape[:-self.ndims] - pdofs = prod(dofs) if dofs[0] not in self._geometry_values[1] or dofs[0]==1 else prod(dofs[1:]) - pdofs *= 2**(self.iscomplex) + dofs = field.shape[: -self.ndims] + pdofs = ( + prod(dofs) + if dofs[0] not in self._geometry_values[1] or dofs[0] == 1 + else prod(dofs[1:]) + ) + pdofs *= 2 ** (self.iscomplex) if not (pdofs in (12, 8, 10) or sqrt(pdofs / 2).is_integer()): raise TypeError(f"Unrecognized field dofs {dofs}") diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 11e205d..8c441ae 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -184,7 +184,7 @@ def copy(self, other=None, out=None, **kwargs): # * check if this dose not cause any bugs if it overwrites ndarray.copy # - For now, we cast self to ndarray before performing ndarray methods like flatten # - the second arg should be "order='C'" to match the signiture? - + # check=False => here any output is accepted out = self.prepare_out(out, check=False, **kwargs) if other is None: @@ -256,9 +256,9 @@ def prepare_in(self, fields, **kwargs): return self.prepare(fields, **kwargs) _children = {} - + def __new__(cls, field, **kwargs): - #TODO: get dofs and local dims from kwargs, instead of getting them + # TODO: get dofs and local dims from kwargs, instead of getting them # from self.shape assuming that it has the form (dofs, local_dims) if isinstance(field, cls): @@ -269,14 +269,14 @@ def __new__(cls, field, **kwargs): ) parent = numpy.ndarray if isinstance(field, numpy.ndarray) else cupy.ndarray if (cls, parent) not in cls._children: - cls._children[(cls,parent)] = type(cls.__name__+"ext",(cls, parent), {}) - obj = field.view(type=cls._children[(cls,parent)]) - #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) - #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) - + cls._children[(cls, parent)] = type(cls.__name__ + "ext", (cls, parent), {}) + obj = field.view(type=cls._children[(cls, parent)]) + # self._dims = kwargs.get("dims", self.shape[-self.ndims :]) + # self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) + return obj - #field check should be performed + # field check should be performed def __array_finalize__(self, obj): "Support for __array_finalize__ standard" # Note: this is called when creating a temporary, possibly causing @@ -285,10 +285,11 @@ def __array_finalize__(self, obj): # self: newly created instance # obj: input instance self._check_field(obj) - if obj is None: return # can be removed; we are not bypassing ndarray.__new__ + if obj is None: + return # can be removed; we are not bypassing ndarray.__new__ # This will require some care when we use attr for dims and dofs in _check_field - #self._dims = kwargs.get("dims", self.shape[-self.ndims :]) - #self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) + # self._dims = kwargs.get("dims", self.shape[-self.ndims :]) + # self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) self.__init__(obj, comm=getattr(obj, "comm", None)) def __init__(self, field, comm=None, **kwargs): @@ -308,7 +309,7 @@ def _check_field(self, field=None): raise TypeError("Field is stored on a different device than the quda lib") if len(field.shape) < 4: raise ValueError("A lattice field should not have shape smaller than 4") - + def activate(self): "Activates the current field. To be called before using the object in quda" "to make sure the communicator is set for MPI" @@ -332,7 +333,7 @@ def complex_view(self): def float_view(self): "Returns a complex view of the field" - #don't need to upcast if we keep dofs and dims as attributes + # don't need to upcast if we keep dofs and dims as attributes if not self.iscomplex: return self.view(type=self.backend.ndarray) return self.view(get_float_dtype(self.dtype), self.backend.ndarray) @@ -358,7 +359,7 @@ def device(self): @property def device_id(self): return getattr(self.device, "id", None) - + @property def location(self): "Memory location of the field (CPU or CUDA)" diff --git a/test/test_field.py b/test/test_field.py index 868eea9..a600455 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -2,6 +2,7 @@ from pytest import raises import numpy as np import cupy as cp + shape = (4, 3, 4, 4, 4, 4) @@ -46,7 +47,7 @@ def test_numpy(): field /= 1 assert field2 == field - + def test_cupy(): field = LatticeField(cp.zeros(shape)) assert field.location == "CUDA" diff --git a/test/test_gauge_field.py b/test/test_gauge_field.py index 88ebd3f..210d364 100644 --- a/test/test_gauge_field.py +++ b/test/test_gauge_field.py @@ -10,7 +10,7 @@ epsilon_loop, ) from lyncs_cppyy.ll import addressof -from lyncs_utils import isclose#, allclose +from lyncs_utils import isclose # , allclose @lattice_loop @@ -234,17 +234,22 @@ def test_force(lib, lattice, device, epsilon): zeros = getattr(gf, path + "_field")(coeffs=0, force=True) assert zeros == 0 + from lyncs_utils import isiterable from collections.abc import Mapping + + def values(dct): "Calls values, if available, or dict.values" try: return dct.values() except AttributeError: return dict.values(dct) + + def allclose(left, right, **kwargs): if isinstance(left, cp.ndarray) and isinstance(right, cp.ndarray): - return np.allclose(left,right) + return np.allclose(left, right) if isinstance(left, cp.ndarray) and not isinstance(right, cp.ndarray): left = [left] * len(right) if not isinstance(left, cp.ndarray) and isinstance(right, cp.ndarray): @@ -264,6 +269,7 @@ def allclose(left, right, **kwargs): pairs = zip(left, right) return all((allclose(*pair, **kwargs) for pair in pairs)) + # @dtype_loop # enables dtype @device_loop # enables device @lattice_loop # enables lattice From b8c3a4f327b62a5e53968ddd941593142b310844 Mon Sep 17 00:00:00 2001 From: Simone Bacchio Date: Wed, 8 Mar 2023 20:47:04 +0200 Subject: [PATCH 09/10] Defining child --- lyncs_quda/lattice_field.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py index 8c441ae..24476c6 100644 --- a/lyncs_quda/lattice_field.py +++ b/lyncs_quda/lattice_field.py @@ -268,9 +268,10 @@ def __new__(cls, field, **kwargs): f"Supporting only numpy or cupy for field, got {type(field)}" ) parent = numpy.ndarray if isinstance(field, numpy.ndarray) else cupy.ndarray - if (cls, parent) not in cls._children: - cls._children[(cls, parent)] = type(cls.__name__ + "ext", (cls, parent), {}) - obj = field.view(type=cls._children[(cls, parent)]) + child = (cls, parent) + if child not in cls._children: + cls._children[child] = type(cls.__name__ + "ext", child, {}) + obj = field.view(type=cls._children[child]) # self._dims = kwargs.get("dims", self.shape[-self.ndims :]) # self._dofs = kwargs.get("dofs", field.shape[: -self.ndims]) From 2750c3789360876e58b71166a4802b9e70fa8e6c Mon Sep 17 00:00:00 2001 From: Simone Bacchio Date: Fri, 17 Mar 2023 09:37:12 +0200 Subject: [PATCH 10/10] Removing usage of allclose --- test/test_gauge_field.py | 45 +++++----------------------------------- 1 file changed, 5 insertions(+), 40 deletions(-) diff --git a/test/test_gauge_field.py b/test/test_gauge_field.py index 210d364..63c2fca 100644 --- a/test/test_gauge_field.py +++ b/test/test_gauge_field.py @@ -10,7 +10,7 @@ epsilon_loop, ) from lyncs_cppyy.ll import addressof -from lyncs_utils import isclose # , allclose +from lyncs_utils import isclose @lattice_loop @@ -235,41 +235,6 @@ def test_force(lib, lattice, device, epsilon): assert zeros == 0 -from lyncs_utils import isiterable -from collections.abc import Mapping - - -def values(dct): - "Calls values, if available, or dict.values" - try: - return dct.values() - except AttributeError: - return dict.values(dct) - - -def allclose(left, right, **kwargs): - if isinstance(left, cp.ndarray) and isinstance(right, cp.ndarray): - return np.allclose(left, right) - if isinstance(left, cp.ndarray) and not isinstance(right, cp.ndarray): - left = [left] * len(right) - if not isinstance(left, cp.ndarray) and isinstance(right, cp.ndarray): - right = [right] * len(left) - if len(left) != len(right): - return False - if isinstance(left, Mapping) or isinstance(right, Mapping): - if not isinstance(left, Mapping): - pairs = zip(left, values(right)) - elif not isinstance(right, Mapping): - pairs = zip(values(left), right) - else: - if set(keys(left)) != set(keys(right)): - return False - pairs = dictzip(left, right, values_only=True) - else: - pairs = zip(left, right) - return all((allclose(*pair, **kwargs) for pair in pairs)) - - # @dtype_loop # enables dtype @device_loop # enables device @lattice_loop # enables lattice @@ -287,9 +252,9 @@ def test_paths_wins(lib, lattice, device): rparts = gf.plaquette_field(sum_paths=False, insertion=mom, left=False) lparts = tuple(zip(*lparts)) rparts = tuple(zip(*rparts)) - assert allclose(lparts[0], rparts[0]) - assert allclose(lparts[1], out) - assert allclose(rparts[1], out) + assert all(np.allclose(left, right) for left, right in zip(lparts[0], rparts[0])) + assert all(np.allclose(left, out) for left in lparts[1]) + assert all(np.allclose(left, out) for left in rparts[1]) mom = momentum(lattice, dtype=dtype, device=device) mom.gaussian() @@ -298,7 +263,7 @@ def test_paths_wins(lib, lattice, device): rparts = gf.plaquette_field(sum_paths=False, insertion=mom, left=False) lparts = tuple(zip(*lparts)) rparts = tuple(zip(*rparts)) - assert allclose(lparts[0], rparts[0]) + assert all(np.allclose(left, right) for left, right in zip(lparts[0], rparts[0])) ltracs = tuple(loop.reduce() for loop in lparts[1]) rtracs = tuple(loop.reduce() for loop in rparts[1]) assert np.allclose(ltracs, rtracs)