From eda70a19f7f7472d4ef8eb470193233a6083bb6e Mon Sep 17 00:00:00 2001 From: Martin Reinecke <martin@mpa-garching.mpg.de> Date: Wed, 4 Dec 2019 20:42:32 +0100 Subject: [PATCH] stage 1 --- .gitlab-ci.yml | 9 +- nifty6/__init__.py | 7 - nifty6/data_objects/distributed_do.py | 590 --------------------- nifty6/data_objects/numpy_do.py | 136 +---- nifty6/dobj.py | 17 +- nifty6/domain_tuple.py | 10 - nifty6/domains/domain.py | 13 - nifty6/domains/lm_space.py | 2 +- nifty6/domains/power_space.py | 19 +- nifty6/domains/rg_space.py | 13 +- nifty6/field.py | 77 +-- nifty6/internal_config.py | 25 - nifty6/library/los_response.py | 11 +- nifty6/logger.py | 10 +- nifty6/minimization/scipy_minimizer.py | 5 - nifty6/operators/diagonal_operator.py | 21 +- nifty6/operators/distributors.py | 28 +- nifty6/operators/field_zero_padder.py | 29 +- nifty6/operators/harmonic_operators.py | 86 +-- nifty6/operators/outer_product_operator.py | 9 +- nifty6/operators/regridding_operator.py | 18 +- nifty6/plot.py | 4 - nifty6/sugar.py | 15 +- test/test_field.py | 98 ++-- 24 files changed, 141 insertions(+), 1111 deletions(-) delete mode 100644 nifty6/data_objects/distributed_do.py delete mode 100644 nifty6/internal_config.py diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b5a74b32f..f05eb5312 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -39,17 +39,10 @@ test_serial: script: - pytest-3 -q --cov=nifty6 test - > - python3 -m coverage report --omit "*plot*,*distributed_do*" | tee coverage.txt + python3 -m coverage report --omit "*plot*" | tee coverage.txt - > grep TOTAL coverage.txt | awk '{ print "TOTAL: "$4; }' -test_mpi: - stage: test - variables: - OMPI_MCA_btl_vader_single_copy_mechanism: none - script: - - mpiexec -n 2 --bind-to none pytest-3 -q test - pages: stage: release script: diff --git a/nifty6/__init__.py b/nifty6/__init__.py index 1ddfd7966..9906e63bf 100644 --- a/nifty6/__init__.py +++ b/nifty6/__init__.py @@ -1,7 +1,5 @@ from .version import __version__ -from . import dobj - from .domains.domain import Domain from .domains.structured_domain import StructuredDomain from .domains.unstructured_domain import UnstructuredDomain @@ -93,10 +91,5 @@ from .linearization import Linearization from .operator_spectrum import operator_spectrum -from . import internal_config -_scheme = internal_config.parallelization_scheme() -if _scheme == "Samples": - from .minimization.metric_gaussian_kl_mpi import MetricGaussianKL_MPI - # We deliberately don't set __all__ here, because we don't want people to do a # "from nifty6 import *"; that would swamp the global namespace. diff --git a/nifty6/data_objects/distributed_do.py b/nifty6/data_objects/distributed_do.py deleted file mode 100644 index d289d6e27..000000000 --- a/nifty6/data_objects/distributed_do.py +++ /dev/null @@ -1,590 +0,0 @@ -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. -# -# Copyright(C) 2013-2019 Max-Planck-Society -# -# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. - -import sys -from functools import reduce - -import numpy as np -from mpi4py import MPI - -from .random import Random - -__all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", - "empty", "zeros", "ones", "empty_like", "vdot", "exp", - "log", "tanh", "sqrt", "from_object", "from_random", - "local_data", "ibegin", "ibegin_from_shape", "np_allreduce_sum", - "np_allreduce_min", "np_allreduce_max", - "distaxis", "from_local_data", "from_global_data", "to_global_data", - "redistribute", "default_distaxis", "is_numpy", "absmax", "norm", - "lock", "locked", "uniform_full", "transpose", "to_global_data_rw", - "ensure_not_distributed", "ensure_default_distributed", - "tanh", "conjugate", "sin", "cos", "tan", "log10", "log1p", "expm1", - "sinh", "cosh", "sinc", "absolute", "sign", "clip"] - -_comm = MPI.COMM_WORLD -ntask = _comm.Get_size() -rank = _comm.Get_rank() -master = (rank == 0) - - -def is_numpy(): - return False - - -def _shareSize(nwork, nshares, myshare): - return (nwork//nshares) + int(myshare < nwork % nshares) - - -def _shareRange(nwork, nshares, myshare): - nbase = nwork//nshares - additional = nwork % nshares - lo = myshare*nbase + min(myshare, additional) - hi = lo + nbase + int(myshare < additional) - return lo, hi - - -def local_shape(shape, distaxis=0): - if len(shape) == 0 or distaxis == -1: - return shape - shape2 = list(shape) - shape2[distaxis] = _shareSize(shape[distaxis], ntask, rank) - return tuple(shape2) - - -class data_object(object): - def __init__(self, shape, data, distaxis): - self._shape = tuple(shape) - if len(self._shape) == 0: - distaxis = -1 - if not isinstance(data, np.ndarray): - data = np.full((), data) - self._distaxis = distaxis - self._data = data - if local_shape(self._shape, self._distaxis) != self._data.shape: - raise ValueError("shape mismatch") - - def copy(self): - return data_object(self._shape, self._data.copy(), self._distaxis) - -# def _sanity_checks(self): -# # check whether the distaxis is consistent -# if self._distaxis < -1 or self._distaxis >= len(self._shape): -# raise ValueError -# itmp = np.array(self._distaxis) -# otmp = np.empty(ntask, dtype=np.int) -# _comm.Allgather(itmp, otmp) -# if np.any(otmp != self._distaxis): -# raise ValueError -# # check whether the global shape is consistent -# itmp = np.array(self._shape) -# otmp = np.empty((ntask, len(self._shape)), dtype=np.int) -# _comm.Allgather(itmp, otmp) -# for i in range(ntask): -# if np.any(otmp[i, :] != self._shape): -# raise ValueError -# # check shape of local data -# if self._distaxis < 0: -# if self._data.shape != self._shape: -# raise ValueError -# else: -# itmp = np.array(self._shape) -# itmp[self._distaxis] = _shareSize(self._shape[self._distaxis], -# ntask, rank) -# if np.any(self._data.shape != itmp): -# raise ValueError - - @property - def dtype(self): - return self._data.dtype - - @property - def shape(self): - return self._shape - - @property - def size(self): - return np.prod(self._shape) - - @property - def real(self): - return data_object(self._shape, self._data.real, self._distaxis) - - @property - def imag(self): - return data_object(self._shape, self._data.imag, self._distaxis) - - def conj(self): - return data_object(self._shape, self._data.conj(), self._distaxis) - - def conjugate(self): - return data_object(self._shape, self._data.conjugate(), self._distaxis) - - def _contraction_helper(self, op, mpiop, axis): - if axis is not None: - if len(axis) == len(self._data.shape): - axis = None - if axis is None: - res = np.array(getattr(self._data, op)()) - if (self._distaxis == -1): - return res[()] - res2 = np.empty((), dtype=res.dtype) - _comm.Allreduce(res, res2, mpiop) - return res2[()] - - if self._distaxis in axis: - res = getattr(self._data, op)(axis=axis) - res2 = np.empty_like(res) - _comm.Allreduce(res, res2, mpiop) - return from_global_data(res2, distaxis=0) - else: - # perform the contraction on the local data - res = getattr(self._data, op)(axis=axis) - if self._distaxis == -1: - return from_global_data(res, distaxis=0) - shp = list(res.shape) - shift = 0 - for ax in axis: - if ax < self._distaxis: - shift += 1 - shp[self._distaxis-shift] = self.shape[self._distaxis] - return from_local_data(shp, res, self._distaxis-shift) - - def sum(self, axis=None): - return self._contraction_helper("sum", MPI.SUM, axis) - - def prod(self, axis=None): - return self._contraction_helper("prod", MPI.PROD, axis) - -# def min(self, axis=None): -# return self._contraction_helper("min", MPI.MIN, axis) - -# def max(self, axis=None): -# return self._contraction_helper("max", MPI.MAX, axis) - - def mean(self, axis=None): - if axis is None: - sz = self.size - else: - sz = reduce(lambda x, y: x*y, [self.shape[i] for i in axis]) - return self.sum(axis)/sz - - def std(self, axis=None): - return np.sqrt(self.var(axis)) - - # FIXME: to be improved! - def var(self, axis=None): - if axis is not None and len(axis) != len(self.shape): - raise ValueError("functionality not yet supported") - return (abs(self-self.mean())**2).mean() - - def _binary_helper(self, other, op): - a = self - if isinstance(other, data_object): - b = other - if a._shape != b._shape: - raise ValueError("shapes are incompatible.") - if a._distaxis != b._distaxis: - raise ValueError("distributions are incompatible.") - a = a._data - b = b._data - elif np.isscalar(other): - a = a._data - b = other - else: - return NotImplemented - - tval = getattr(a, op)(b) - if tval is a: - return self - else: - return data_object(self._shape, tval, self._distaxis) - - def clip(self, min=None, max=None): - return data_object(self._shape, np.clip(self._data, min, max)) - - def __neg__(self): - return data_object(self._shape, -self._data, self._distaxis) - - def __abs__(self): - return data_object(self._shape, abs(self._data), self._distaxis) - - def all(self): - return self.sum() == self.size - - def any(self): - return self.sum() != 0 - - def fill(self, value): - self._data.fill(value) - - -for op in ["__add__", "__radd__", "__iadd__", - "__sub__", "__rsub__", "__isub__", - "__mul__", "__rmul__", "__imul__", - "__truediv__", "__rtruediv__", "__itruediv__", - "__floordiv__", "__rfloordiv__", "__ifloordiv__", - "__pow__", "__rpow__", "__ipow__", - "__lt__", "__le__", "__gt__", "__ge__", "__eq__", "__ne__"]: - def func(op): - def func2(self, other): - return self._binary_helper(other, op=op) - return func2 - setattr(data_object, op, func(op)) - - -def full(shape, fill_value, dtype=None, distaxis=0): - return data_object(shape, np.full(local_shape(shape, distaxis), - fill_value, dtype), distaxis) - - -def uniform_full(shape, fill_value, dtype=None, distaxis=0): - return data_object( - shape, np.broadcast_to(fill_value, local_shape(shape, distaxis)), - distaxis) - - -def empty(shape, dtype=None, distaxis=0): - return data_object(shape, np.empty(local_shape(shape, distaxis), - dtype), distaxis) - - -def zeros(shape, dtype=None, distaxis=0): - return data_object(shape, np.zeros(local_shape(shape, distaxis), dtype), - distaxis) - - -def ones(shape, dtype=None, distaxis=0): - return data_object(shape, np.ones(local_shape(shape, distaxis), dtype), - distaxis) - - -def empty_like(a, dtype=None): - return data_object(np.empty_like(a._data, dtype)) - - -def vdot(a, b): - tmp = np.array(np.vdot(a._data, b._data)) - if a._distaxis == -1: - return tmp[()] - res = np.empty((), dtype=tmp.dtype) - _comm.Allreduce(tmp, res, MPI.SUM) - return res[()] - - -def _math_helper(x, function, out): - function = getattr(np, function) - if out is not None: - function(x._data, out=out._data) - return out - else: - return data_object(x.shape, function(x._data), x._distaxis) - - -_current_module = sys.modules[__name__] - -for f in ["sqrt", "exp", "log", "tanh", "conjugate", "sin", "cos", "tan", - "sinh", "cosh", "sinc", "absolute", "sign", "log10", "log1p", - "expm1"]: - def func(f): - def func2(x, out=None): - return _math_helper(x, f, out) - return func2 - setattr(_current_module, f, func(f)) - - -def clip(x, a_min=None, a_max=None): - return data_object(x.shape, np.clip(x._data, a_min, a_max), x._distaxis) - - -def from_object(object, dtype, copy, set_locked): - if dtype is None: - dtype = object.dtype - dtypes_equal = dtype == object.dtype - if set_locked and dtypes_equal and locked(object): - return object - if not dtypes_equal and not copy: - raise ValueError("cannot change data type without copying") - if set_locked and not copy: - raise ValueError("cannot lock object without copying") - data = np.array(object._data, dtype=dtype, copy=copy) - if set_locked: - data.flags.writeable = False - return data_object(object._shape, data, distaxis=object._distaxis) - - -# This function draws all random numbers on all tasks, to produce the same -# array independent on the number of tasks -# MR FIXME: depending on what is really wanted/needed (i.e. same result -# independent of number of tasks, performance etc.) we need to adjust the -# algorithm. -def from_random(random_type, shape, dtype=np.float64, **kwargs): - generator_function = getattr(Random, random_type) - if len(shape) == 0: - ldat = generator_function(dtype=dtype, shape=shape, **kwargs) - ldat = _comm.bcast(ldat) - return from_local_data(shape, ldat, distaxis=-1) - for i in range(ntask): - lshape = list(shape) - lshape[0] = _shareSize(shape[0], ntask, i) - ldat = generator_function(dtype=dtype, shape=lshape, **kwargs) - if i == rank: - outdat = ldat - return from_local_data(shape, outdat, distaxis=0) - - -def local_data(arr): - return arr._data - - -def ibegin_from_shape(glob_shape, distaxis=0): - res = [0] * len(glob_shape) - if distaxis < 0: - return res - res[distaxis] = _shareRange(glob_shape[distaxis], ntask, rank)[0] - return tuple(res) - - -def ibegin(arr): - res = [0] * arr._data.ndim - res[arr._distaxis] = _shareRange(arr._shape[arr._distaxis], ntask, rank)[0] - return tuple(res) - - -def np_allreduce_sum(arr): - res = np.empty_like(arr) - _comm.Allreduce(arr, res, MPI.SUM) - return res - - -def np_allreduce_min(arr): - res = np.empty_like(arr) - _comm.Allreduce(arr, res, MPI.MIN) - return res - - -def np_allreduce_max(arr): - res = np.empty_like(arr) - _comm.Allreduce(arr, res, MPI.MAX) - return res - - -def distaxis(arr): - return arr._distaxis - - -def from_local_data(shape, arr, distaxis=0): - if arr.dtype.kind not in "fciub": - raise TypeError - return data_object(shape, arr, distaxis) - - -def from_global_data(arr, sum_up=False, distaxis=0): - if arr.dtype.kind not in "fciub": - raise TypeError - if sum_up: - arr = np_allreduce_sum(arr) - if arr.ndim == 0: - distaxis = -1 - if distaxis == -1: - return data_object(arr.shape, arr, distaxis) - lo, hi = _shareRange(arr.shape[distaxis], ntask, rank) - sl = [slice(None)]*len(arr.shape) - sl[distaxis] = slice(lo, hi) - return data_object(arr.shape, arr[tuple(sl)], distaxis) - - -def to_global_data(arr): - if arr._distaxis == -1: - return arr._data - tmp = redistribute(arr, dist=-1) - return tmp._data - - -def to_global_data_rw(arr): - if arr._distaxis == -1: - return arr._data.copy() - tmp = redistribute(arr, dist=-1) - return tmp._data - - -def redistribute(arr, dist=None, nodist=None): - if dist is not None: - if nodist is not None: - raise ValueError - if dist == arr._distaxis: - return arr - else: - if nodist is None: - raise ValueError - if arr._distaxis not in nodist: - return arr - dist = -1 - for i in range(len(arr.shape)): - if i not in nodist: - dist = i - break - - if arr._distaxis == -1: # all data available, just pick the proper subset - return from_global_data(arr._data, distaxis=dist) - if dist == -1: # gather all data on all tasks - tmp = np.moveaxis(arr._data, arr._distaxis, 0) - slabsize = np.prod(tmp.shape[1:])*tmp.itemsize - sz = np.empty(ntask, dtype=np.int) - for i in range(ntask): - sz[i] = slabsize*_shareSize(arr.shape[arr._distaxis], ntask, i) - disp = np.empty(ntask, dtype=np.int) - disp[0] = 0 - disp[1:] = np.cumsum(sz[:-1]) - tmp = np.require(tmp, requirements="C") - out = np.empty(arr.size, dtype=arr.dtype) - _comm.Allgatherv(tmp, [out, sz, disp, MPI.BYTE]) - shp = np.array(arr._shape) - shp[1:arr._distaxis+1] = shp[0:arr._distaxis] - shp[0] = arr.shape[arr._distaxis] - out = out.reshape(shp) - out = np.moveaxis(out, 0, arr._distaxis) - return from_global_data(out, distaxis=-1) - - # real redistribution via Alltoallv - ssz0 = arr._data.size//arr.shape[dist] - ssz = np.empty(ntask, dtype=np.int) - rszall = arr.size//arr.shape[dist]*_shareSize(arr.shape[dist], ntask, rank) - rbuf = np.empty(rszall, dtype=arr.dtype) - rsz0 = rszall//arr.shape[arr._distaxis] - rsz = np.empty(ntask, dtype=np.int) - if dist == 0: # shortcut possible - sbuf = np.ascontiguousarray(arr._data) - for i in range(ntask): - lo, hi = _shareRange(arr.shape[dist], ntask, i) - ssz[i] = ssz0*(hi-lo) - rsz[i] = rsz0*_shareSize(arr.shape[arr._distaxis], ntask, i) - else: - sbuf = np.empty(arr._data.size, dtype=arr.dtype) - sslice = [slice(None)]*arr._data.ndim - ofs = 0 - for i in range(ntask): - lo, hi = _shareRange(arr.shape[dist], ntask, i) - sslice[dist] = slice(lo, hi) - ssz[i] = ssz0*(hi-lo) - sbuf[ofs:ofs+ssz[i]] = arr._data[tuple(sslice)].flat - ofs += ssz[i] - rsz[i] = rsz0*_shareSize(arr.shape[arr._distaxis], ntask, i) - ssz *= arr._data.itemsize - rsz *= arr._data.itemsize - sdisp = np.append(0, np.cumsum(ssz[:-1])) - rdisp = np.append(0, np.cumsum(rsz[:-1])) - s_msg = [sbuf, (ssz, sdisp), MPI.BYTE] - r_msg = [rbuf, (rsz, rdisp), MPI.BYTE] - _comm.Alltoallv(s_msg, r_msg) - del sbuf # free memory - if arr._distaxis == 0: - rbuf = rbuf.reshape(local_shape(arr.shape, dist)) - arrnew = from_local_data(arr.shape, rbuf, distaxis=dist) - else: - arrnew = np.empty(local_shape(arr.shape, dist), dtype=arr.dtype) - rslice = [slice(None)]*arr._data.ndim - ofs = 0 - for i in range(ntask): - lo, hi = _shareRange(arr.shape[arr._distaxis], ntask, i) - rslice[arr._distaxis] = slice(lo, hi) - sz = rsz[i]//arr._data.itemsize - arrnew[tuple(rslice)].flat = rbuf[ofs:ofs+sz] - ofs += sz - arrnew = from_local_data(arr.shape, arrnew, distaxis=dist) - return arrnew - - -def transpose(arr): - if len(arr.shape) != 2 or arr._distaxis != 0: - raise ValueError("bad input") - ssz0 = arr._data.size//arr.shape[1] - ssz = np.empty(ntask, dtype=np.int) - rszall = arr.size//arr.shape[1]*_shareSize(arr.shape[1], ntask, rank) - rbuf = np.empty(rszall, dtype=arr.dtype) - rsz0 = rszall//arr.shape[0] - rsz = np.empty(ntask, dtype=np.int) - sbuf = np.empty(arr._data.size, dtype=arr.dtype) - ofs = 0 - for i in range(ntask): - lo, hi = _shareRange(arr.shape[1], ntask, i) - ssz[i] = ssz0*(hi-lo) - sbuf[ofs:ofs+ssz[i]] = arr._data[:, lo:hi].flat - ofs += ssz[i] - rsz[i] = rsz0*_shareSize(arr.shape[0], ntask, i) - ssz *= arr._data.itemsize - rsz *= arr._data.itemsize - sdisp = np.append(0, np.cumsum(ssz[:-1])) - rdisp = np.append(0, np.cumsum(rsz[:-1])) - s_msg = [sbuf, (ssz, sdisp), MPI.BYTE] - r_msg = [rbuf, (rsz, rdisp), MPI.BYTE] - _comm.Alltoallv(s_msg, r_msg) - del sbuf # free memory - sz2 = _shareSize(arr.shape[1], ntask, rank) - arrnew = np.empty((sz2, arr.shape[0]), dtype=arr.dtype) - ofs = 0 - for i in range(ntask): - lo, hi = _shareRange(arr.shape[0], ntask, i) - sz = rsz[i]//arr._data.itemsize - arrnew[:, lo:hi] = rbuf[ofs:ofs+sz].reshape(hi-lo, sz2).T - ofs += sz - return from_local_data((arr.shape[1], arr.shape[0]), arrnew, 0) - - -def default_distaxis(): - return 0 - - -def lock(arr): - arr._data.flags.writeable = False - - -def locked(arr): - return not arr._data.flags.writeable - - -def ensure_not_distributed(arr, axes): - if arr._distaxis in axes: - arr = redistribute(arr, nodist=axes) - return arr, arr._data - - -def ensure_default_distributed(arr): - if arr._distaxis != 0: - arr = redistribute(arr, dist=0) - return arr - - -def absmax(arr): - if arr._data.size == 0: - tmp = np.array(0, dtype=arr._data.dtype) - else: - tmp = np.asarray(np.linalg.norm(arr._data.reshape(-1), ord=np.inf)) - res = np.empty_like(tmp) - _comm.Allreduce(tmp, res, MPI.MAX) - return res[()] - - -def norm(arr, ord=2): - if ord == np.inf: - return absmax(arr) - tmp = np.asarray(np.linalg.norm(arr._data.reshape(-1), ord=ord) ** ord) - res = np.empty_like(tmp) - if len(arr._data.shape) == 0: - res = tmp - else: - _comm.Allreduce(tmp, res, MPI.SUM) - return res[()] ** (1./ord) diff --git a/nifty6/data_objects/numpy_do.py b/nifty6/data_objects/numpy_do.py index a058c6ed4..bb34d2121 100644 --- a/nifty6/data_objects/numpy_do.py +++ b/nifty6/data_objects/numpy_do.py @@ -18,141 +18,15 @@ # Data object module that uses simple numpy ndarrays. import numpy as np -from numpy import ndarray as data_object -from numpy import empty, empty_like, ones, zeros, full -from numpy import absolute, sign, clip, vdot -from numpy import sin, cos, sinh, cosh, tan, tanh -from numpy import exp, log, log10, sqrt, sinc, log1p, expm1 +#from numpy import ndarray as data_object +#from numpy import empty, empty_like, ones, zeros, full +#from numpy import absolute, sign, clip, vdot +#from numpy import sin, cos, sinh, cosh, tan, tanh +#from numpy import exp, log, log10, sqrt, sinc, log1p, expm1 from .random import Random -__all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", - "empty", "zeros", "ones", "empty_like", "vdot", "exp", - "log", "tanh", "sqrt", "from_object", "from_random", - "local_data", "ibegin", "ibegin_from_shape", "np_allreduce_sum", - "np_allreduce_min", "np_allreduce_max", - "distaxis", "from_local_data", "from_global_data", "to_global_data", - "redistribute", "default_distaxis", "is_numpy", "absmax", "norm", - "lock", "locked", "uniform_full", "to_global_data_rw", - "ensure_not_distributed", "ensure_default_distributed", - "clip", "sin", "cos", "tan", "sinh", "cosh", - "absolute", "sign", "sinc", "log10", "log1p", "expm1"] - -ntask = 1 -rank = 0 -master = True - - -def is_numpy(): - return True - - -def from_object(object, dtype, copy, set_locked): - if dtype is None: - dtype = object.dtype - dtypes_equal = dtype == object.dtype - if set_locked and dtypes_equal and locked(object): - return object - if not dtypes_equal and not copy: - raise ValueError("cannot change data type without copying") - if set_locked and not copy: - raise ValueError("cannot lock object without copying") - res = np.array(object, dtype=dtype, copy=copy) - if set_locked: - lock(res) - return res - def from_random(random_type, shape, dtype=np.float64, **kwargs): generator_function = getattr(Random, random_type) return generator_function(dtype=dtype, shape=shape, **kwargs) - - -def local_data(arr): - return arr - - -def ibegin_from_shape(glob_shape, distaxis=-1): - return (0,)*len(glob_shape) - - -def ibegin(arr): - return (0,)*arr.ndim - - -def np_allreduce_sum(arr): - return arr - - -def np_allreduce_min(arr): - return arr - - -def np_allreduce_max(arr): - return arr - - -def distaxis(arr): - return -1 - - -def from_local_data(shape, arr, distaxis=-1): - if tuple(shape) != arr.shape: - raise ValueError - if arr.dtype.kind not in "fciub": - raise TypeError - return arr - - -def from_global_data(arr, sum_up=False, distaxis=-1): - if arr.dtype.kind not in "fciub": - raise TypeError - return arr - - -def to_global_data(arr): - return arr - - -def to_global_data_rw(arr): - return arr.copy() - - -def redistribute(arr, dist=None, nodist=None): - return arr - - -def default_distaxis(): - return -1 - - -def local_shape(glob_shape, distaxis=-1): - return glob_shape - - -def lock(arr): - arr.flags.writeable = False - - -def locked(arr): - return not arr.flags.writeable - - -def uniform_full(shape, fill_value, dtype=None, distaxis=-1): - return np.broadcast_to(fill_value, shape) - - -def ensure_not_distributed(arr, axes): - return arr, arr - - -def ensure_default_distributed(arr): - return arr - - -def absmax(arr): - return np.linalg.norm(arr.reshape(-1), ord=np.inf) - - -def norm(arr, ord=2): - return np.linalg.norm(arr.reshape(-1), ord=ord) diff --git a/nifty6/dobj.py b/nifty6/dobj.py index 9f3d3dc7e..fed5881fd 100644 --- a/nifty6/dobj.py +++ b/nifty6/dobj.py @@ -15,19 +15,4 @@ # # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. -from . import internal_config - - -_scheme = internal_config.parallelization_scheme() - -if _scheme in ("Samples", "None"): - from .data_objects.numpy_do import * -else: - try: - from mpi4py import MPI - if MPI.COMM_WORLD.Get_size() == 1: - from .data_objects.numpy_do import * - else: - from .data_objects.distributed_do import * - except ImportError: - from .data_objects.numpy_do import * +from .data_objects.numpy_do import * diff --git a/nifty6/domain_tuple.py b/nifty6/domain_tuple.py index 9200b85f3..84df832ca 100644 --- a/nifty6/domain_tuple.py +++ b/nifty6/domain_tuple.py @@ -109,16 +109,6 @@ class DomainTuple(object): """ return self._shape - @property - def local_shape(self): - """tuple of int: number of pixels along each axis on the local task - - The shape of the array-like object required to store information - defined on part of the domain which is stored on the local MPI task. - """ - from .dobj import local_shape - return local_shape(self._shape) - @property def size(self): """int : total number of pixels. diff --git a/nifty6/domains/domain.py b/nifty6/domains/domain.py index c5c332eda..a5bed8306 100644 --- a/nifty6/domains/domain.py +++ b/nifty6/domains/domain.py @@ -84,19 +84,6 @@ class Domain(metaclass=NiftyMeta): """ raise NotImplementedError - @property - def local_shape(self): - """tuple of int: number of pixels along each axis on the local task, - mainly relevant for MPI. - - See :meth:`.shape()` for general explanation of property. - - The shape of the array-like object required to store information - defined on part of the domain which is stored on the local MPI task. - """ - from ..dobj import local_shape - return local_shape(self.shape) - @property def size(self): """int: total number of pixels. diff --git a/nifty6/domains/lm_space.py b/nifty6/domains/lm_space.py index 081106ee4..520d1c878 100644 --- a/nifty6/domains/lm_space.py +++ b/nifty6/domains/lm_space.py @@ -87,7 +87,7 @@ class LMSpace(StructuredDomain): for m in range(1, mmax+1): ldist[idx:idx+2*(lmax+1-m)] = tmp[2*m:] idx += 2*(lmax+1-m) - return Field.from_global_data(self, ldist) + return Field.from_arr(self, ldist) def get_unique_k_lengths(self): return np.arange(self.lmax+1, dtype=np.float64) diff --git a/nifty6/domains/power_space.py b/nifty6/domains/power_space.py index d200ce1b9..f5774e3fd 100644 --- a/nifty6/domains/power_space.py +++ b/nifty6/domains/power_space.py @@ -17,7 +17,6 @@ import numpy as np -from .. import dobj from .structured_domain import StructuredDomain @@ -176,26 +175,20 @@ class PowerSpace(StructuredDomain): tbb = 0.5*(tmp[:-1]+tmp[1:]) else: tbb = binbounds - locdat = np.searchsorted(tbb, k_length_array.local_data) - temp_pindex = dobj.from_local_data( - k_length_array.val.shape, locdat, - dobj.distaxis(k_length_array.val)) + temp_pindex = np.searchsorted(tbb, k_length_array.val) nbin = len(tbb)+1 - temp_rho = np.bincount(dobj.local_data(temp_pindex).ravel(), - minlength=nbin) - temp_rho = dobj.np_allreduce_sum(temp_rho) + temp_rho = np.bincount(temp_pindex.ravel(), minlength=nbin) if (temp_rho == 0).any(): raise ValueError("empty bins detected") # The explicit conversion to float64 is necessary because bincount # sometimes returns its result as an integer array, even when # floating-point weights are present ... - temp_k_lengths = np.bincount( - dobj.local_data(temp_pindex).ravel(), - weights=k_length_array.local_data.ravel(), + temp_k_lengths = np.bincount(temp_pindex.ravel(), + weights=k_length_array.val.ravel(), minlength=nbin).astype(np.float64, copy=False) - temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho + temp_k_lengths = temp_k_lengths / temp_rho temp_k_lengths.flags.writeable = False - dobj.lock(temp_pindex) + temp_pindex.flags.writeable = False temp_dvol = temp_rho*pdvol temp_dvol.flags.writeable = False self._powerIndexCache[key] = (binbounds, temp_pindex, diff --git a/nifty6/domains/rg_space.py b/nifty6/domains/rg_space.py index 344af55fb..3d0abbdcc 100644 --- a/nifty6/domains/rg_space.py +++ b/nifty6/domains/rg_space.py @@ -18,7 +18,6 @@ from functools import reduce import numpy as np -from .. import dobj from ..field import Field from .structured_domain import StructuredDomain @@ -95,18 +94,17 @@ class RGSpace(StructuredDomain): return self._dvol def _get_dist_array(self): - ibegin = dobj.ibegin_from_shape(self._shape) - res = np.arange(self.local_shape[0], dtype=np.float64) + ibegin[0] + res = np.arange(self.shape[0], dtype=np.float64) res = np.minimum(res, self.shape[0]-res)*self.distances[0] if len(self.shape) == 1: - return Field.from_local_data(self, res) + return Field.from_arr(self, res) res *= res for i in range(1, len(self.shape)): - tmp = np.arange(self.local_shape[i], dtype=np.float64) + ibegin[i] + tmp = np.arange(self.shape[i], dtype=np.float64) tmp = np.minimum(tmp, self.shape[i]-tmp)*self.distances[i] tmp *= tmp res = np.add.outer(res, tmp) - return Field.from_local_data(self, np.sqrt(res)) + return Field.from_arr(self, np.sqrt(res)) def get_k_length_array(self): if (not self.harmonic): @@ -133,8 +131,7 @@ class RGSpace(StructuredDomain): tmp[t2] = True return np.sqrt(np.nonzero(tmp)[0])*self.distances[0] else: # do it the hard way - # FIXME: this needs to improve for MPI. Maybe unique()/gather()? - tmp = self.get_k_length_array().to_global_data() + tmp = self.get_k_length_array().val tmp = np.unique(tmp) tol = 1e-12*tmp[-1] # remove all points that are closer than tol to their right diff --git a/nifty6/field.py b/nifty6/field.py index e7572c03f..cf706eba3 100644 --- a/nifty6/field.py +++ b/nifty6/field.py @@ -18,7 +18,7 @@ from functools import reduce import numpy as np -from . import dobj, utilities +from . import utilities from .domain_tuple import DomainTuple @@ -47,16 +47,16 @@ class Field(object): def __init__(self, domain, val): if not isinstance(domain, DomainTuple): raise TypeError("domain must be of type DomainTuple") - if type(val) is not dobj.data_object: + if not isinstance(val, np.ndarray): if np.isscalar(val): - val = dobj.full(domain.shape, val) + val = np.full(domain.shape, val) else: - raise TypeError("val must be of type dobj.data_object") + raise TypeError("val must be of type numpy.ndarray") if domain.shape != val.shape: raise ValueError("shape mismatch between val and domain") self._domain = domain self._val = val - dobj.lock(self._val) + self._val.flags.writeable = False @staticmethod def scalar(val): @@ -93,7 +93,7 @@ class Field(object): return Field(domain, val) @staticmethod - def from_global_data(domain, arr, sum_up=False): + def from_arr(domain, arr): """Returns a Field constructed from `domain` and `arr`. Parameters @@ -103,49 +103,8 @@ class Field(object): arr : numpy.ndarray The data content to be used for the new Field. Its shape must match the shape of `domain`. - sum_up : bool, optional - If True, the contents of `arr` are summed up over all MPI tasks - (if any), and the sum is used as data content. - If False, the contens of `arr` are used directly, and must be - identical on all MPI tasks. """ - return Field(DomainTuple.make(domain), - dobj.from_global_data(arr, sum_up)) - - @staticmethod - def from_local_data(domain, arr): - return Field(DomainTuple.make(domain), - dobj.from_local_data(domain.shape, arr)) - - def to_global_data(self): - """Returns an array containing the full data of the field. - - Returns - ------- - numpy.ndarray : array containing all field entries. - Its shape is identical to `self.shape`. - """ - return dobj.to_global_data(self._val) - - def to_global_data_rw(self): - """Returns a modifiable array containing the full data of the field. - - Returns - ------- - numpy.ndarray - Array containing all field entries, which can be modified. Its - shape is identical to `self.shape`. - """ - return dobj.to_global_data_rw(self._val) - - @property - def local_data(self): - """numpy.ndarray : locally residing field data - - Returns a handle to the part of the array data residing on the local - task (or to the entore array if MPI is not active). - """ - return dobj.local_data(self._val) + return Field(DomainTuple.make(domain), arr) def cast_domain(self, new_domain): """Returns a field with the same data, but a different domain @@ -181,6 +140,7 @@ class Field(object): Field The newly created Field. """ + from . import dobj domain = DomainTuple.make(domain) return Field(domain=domain, val=dobj.from_random(random_type, dtype=dtype, @@ -188,7 +148,7 @@ class Field(object): @property def val(self): - """dobj.data_object : the data object storing the field's entries. + """numpy.ndarray : the data object storing the field's entries. Notes ----- @@ -281,7 +241,7 @@ class Field(object): Field The weighted field. """ - aout = self.local_data.copy() + aout = self.val.copy() spaces = utilities.parse_spaces(spaces, len(self._domain)) @@ -295,15 +255,12 @@ class Field(object): new_shape[self._domain.axes[ind][0]: self._domain.axes[ind][-1]+1] = wgt.shape wgt = wgt.reshape(new_shape) - if dobj.distaxis(self._val) >= 0 and ind == 0: - # we need to distribute the weights along axis 0 - wgt = dobj.local_data(dobj.from_global_data(wgt)) aout *= wgt**power fct = fct**power if fct != 1.: aout *= fct - return Field.from_local_data(self._domain, aout) + return Field(self._domain, aout) def outer(self, x): """Computes the outer product of 'self' with x. @@ -351,7 +308,7 @@ class Field(object): spaces = utilities.parse_spaces(spaces, ndom) if len(spaces) == ndom: - return dobj.vdot(self._val, x._val) + return np.vdot(self._val, x._val) # If we arrive here, we have to do a partial dot product. # For the moment, do this the explicit, non-optimized way return (self.conjugate()*x).sum(spaces=spaces) @@ -369,7 +326,7 @@ class Field(object): float The L2-norm of the field values. """ - return dobj.norm(self._val, ord) + return np.linalg.norm(self._val.reshape(-1), ord=ord) def conjugate(self): """Returns the complex conjugate of the field. @@ -622,9 +579,9 @@ class Field(object): return 0.5*(1.+self.tanh()) def clip(self, min=None, max=None): - min = min.local_data if isinstance(min, Field) else min - max = max.local_data if isinstance(max, Field) else max - return Field(self._domain, dobj.clip(self._val, min, max)) + min = min.val if isinstance(min, Field) else min + max = max.val if isinstance(max, Field) else max + return Field(self._domain, np.clip(self._val, min, max)) def one_over(self): return 1/self @@ -667,6 +624,6 @@ for f in ["sqrt", "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "tanh", "absolute", "sinc", "sign", "log10", "log1p", "expm1"]: def func(f): def func2(self): - return Field(self._domain, getattr(dobj, f)(self.val)) + return Field(self._domain, getattr(np, f)(self.val)) return func2 setattr(Field, f, func(f)) diff --git a/nifty6/internal_config.py b/nifty6/internal_config.py deleted file mode 100644 index c6c157aa0..000000000 --- a/nifty6/internal_config.py +++ /dev/null @@ -1,25 +0,0 @@ -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. -# -# Copyright(C) 2019 Max-Planck-Society -# -# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. - - -# Internal configuration switches, typically for experimental features. -# Leave unchanged unless you know what you are doing! - -def parallelization_scheme(): - return "Standard" - # return "Samples" - # return "None" diff --git a/nifty6/library/los_response.py b/nifty6/library/los_response.py index f094c097a..11042f094 100644 --- a/nifty6/library/los_response.py +++ b/nifty6/library/los_response.py @@ -20,7 +20,6 @@ from scipy.sparse import coo_matrix from scipy.sparse.linalg import aslinearoperator from scipy.special import erfc -from .. import dobj from ..domain_tuple import DomainTuple from ..domains.rg_space import RGSpace from ..domains.unstructured_domain import UnstructuredDomain @@ -164,11 +163,6 @@ class LOSResponse(LinearOperator): if starts.shape != ends.shape: raise TypeError("dimension mismatch") - self._local_shape = dobj.local_shape(self.domain[0].shape) - local_zero_point = (np.array( - dobj.ibegin_from_shape(self.domain[0].shape)) * - np.array(self.domain[0].distances)) - diffs = ends-starts difflen = np.linalg.norm(diffs, axis=0) diffs /= difflen @@ -177,10 +171,9 @@ class LOSResponse(LinearOperator): raise ValueError("parallax error truncation to high: " "getting negative distances") real_ends = starts + diffs*real_distances - lzp = local_zero_point.reshape((-1, 1)) dist = np.array(self.domain[0].distances).reshape((-1, 1)) - localized_pixel_starts = (starts-lzp)/dist + 0.5 - localized_pixel_ends = (real_ends-lzp)/dist + 0.5 + localized_pixel_starts = starts/dist + 0.5 + localized_pixel_ends = real_ends/dist + 0.5 # get the shape of the local data slice w_i = _comp_traverse(localized_pixel_starts, diff --git a/nifty6/logger.py b/nifty6/logger.py index 29ba83215..1754d08c2 100644 --- a/nifty6/logger.py +++ b/nifty6/logger.py @@ -18,16 +18,12 @@ def _logger_init(): import logging - from . import dobj res = logging.getLogger('NIFTy6') res.setLevel(logging.DEBUG) res.propagate = False - if dobj.rank == 0: - ch = logging.StreamHandler() - ch.setLevel(logging.DEBUG) - res.addHandler(ch) - else: - res.addHandler(logging.NullHandler()) + ch = logging.StreamHandler() + ch.setLevel(logging.DEBUG) + res.addHandler(ch) return res diff --git a/nifty6/minimization/scipy_minimizer.py b/nifty6/minimization/scipy_minimizer.py index 64da960ba..21f7bbc31 100644 --- a/nifty6/minimization/scipy_minimizer.py +++ b/nifty6/minimization/scipy_minimizer.py @@ -17,7 +17,6 @@ import numpy as np -from .. import dobj from ..field import Field from ..logger import logger from ..multi_field import MultiField @@ -105,8 +104,6 @@ class _ScipyMinimizer(Minimizer): """ def __init__(self, method, options, need_hessp, bounds): - if not dobj.is_numpy(): - raise NotImplementedError self._method = method self._options = options self._need_hessp = need_hessp @@ -155,8 +152,6 @@ class _ScipyCG(Minimizer): gradient implementation and should not be used otherwise. """ def __init__(self, tol, maxiter): - if not dobj.is_numpy(): - raise NotImplementedError self._tol = tol self._maxiter = maxiter diff --git a/nifty6/operators/diagonal_operator.py b/nifty6/operators/diagonal_operator.py index 1cce2068b..8470428d7 100644 --- a/nifty6/operators/diagonal_operator.py +++ b/nifty6/operators/diagonal_operator.py @@ -17,7 +17,7 @@ import numpy as np -from .. import dobj, utilities +from .. import utilities from ..domain_tuple import DomainTuple from ..field import Field from .endomorphic_operator import EndomorphicOperator @@ -78,16 +78,12 @@ class DiagonalOperator(EndomorphicOperator): for space_index in self._spaces: active_axes += self._domain.axes[space_index] - if self._spaces[0] == 0: - self._ldiag = diagonal.local_data - else: - self._ldiag = diagonal.to_global_data() - locshape = dobj.local_shape(self._domain.shape, 0) + self._ldiag = diagonal.val self._reshaper = [shp if i in active_axes else 1 - for i, shp in enumerate(locshape)] + for i, shp in enumerate(self._domain.shape)] self._ldiag = self._ldiag.reshape(self._reshaper) else: - self._ldiag = diagonal.local_data + self._ldiag = diagonal.val self._fill_rest() def _fill_rest(self): @@ -95,8 +91,7 @@ class DiagonalOperator(EndomorphicOperator): self._complex = utilities.iscomplextype(self._ldiag.dtype) self._capability = self._all_ops if not self._complex: - lmin = self._ldiag.min() if self._ldiag.size > 0 else 1. - self._diagmin = dobj.np_allreduce_min(np.array(lmin))[()] + self._diagmin = self._ldiag.min() def _from_ldiag(self, spc, ldiag): res = DiagonalOperator.__new__(DiagonalOperator) @@ -160,10 +155,10 @@ class DiagonalOperator(EndomorphicOperator): (self._diagmin == 0. and from_inverse)): raise ValueError("operator not positive definite") if from_inverse: - res = samp.local_data/np.sqrt(self._ldiag) + res = samp.val/np.sqrt(self._ldiag) else: - res = samp.local_data*np.sqrt(self._ldiag) - return Field.from_local_data(self._domain, res) + res = samp.val*np.sqrt(self._ldiag) + return Field(self._domain, res) def draw_sample(self, from_inverse=False, dtype=np.float64): res = Field.from_random(random_type="normal", domain=self._domain, diff --git a/nifty6/operators/distributors.py b/nifty6/operators/distributors.py index c6e9ef243..9c06113da 100644 --- a/nifty6/operators/distributors.py +++ b/nifty6/operators/distributors.py @@ -17,7 +17,6 @@ import numpy as np -from .. import dobj from ..domain_tuple import DomainTuple from ..domains.dof_space import DOFSpace from ..domains.power_space import PowerSpace @@ -70,7 +69,7 @@ class DOFDistributor(LinearOperator): nbin = 0 else: nbin = ldat.max() - nbin = dobj.np_allreduce_max(np.array(nbin))[()] + 1 + nbin = nbin + 1 if partner.scalar_dvol is not None: wgt = np.bincount(dofdex.local_data.ravel(), minlength=nbin) wgt = wgt*partner.scalar_dvol @@ -82,7 +81,6 @@ class DOFDistributor(LinearOperator): # sometimes returns its result as an integer array, even when # floating-point weights are present ... wgt = wgt.astype(np.float64, copy=False) - wgt = dobj.np_allreduce_sum(wgt) if (wgt == 0).any(): raise ValueError("empty bins detected") @@ -95,42 +93,30 @@ class DOFDistributor(LinearOperator): self._domain = DomainTuple.make(dom) self._capability = self.TIMES | self.ADJOINT_TIMES - if dobj.default_distaxis() in self._domain.axes[self._space]: - dofdex = dobj.local_data(dofdex) - else: # dofdex must be available fully on every task - dofdex = dobj.to_global_data(dofdex) self._dofdex = dofdex.ravel() firstaxis = self._target.axes[self._space][0] lastaxis = self._target.axes[self._space][-1] - arrshape = dobj.local_shape(self._target.shape, 0) + arrshape = self._target.shape presize = np.prod(arrshape[0:firstaxis], dtype=np.int) postsize = np.prod(arrshape[lastaxis+1:], dtype=np.int) self._hshape = (presize, self._domain[self._space].shape[0], postsize) self._pshape = (presize, self._dofdex.size, postsize) def _adjoint_times(self, x): - arr = x.local_data + arr = x.val arr = arr.reshape(self._pshape) oarr = np.zeros(self._hshape, dtype=x.dtype) oarr = special_add_at(oarr, 1, self._dofdex, arr) - if dobj.distaxis(x.val) in x.domain.axes[self._space]: - oarr = oarr.reshape(self._domain.shape) - res = Field.from_global_data(self._domain, oarr, sum_up=True) - else: - oarr = oarr.reshape(self._domain.local_shape) - res = Field.from_local_data(self._domain, oarr) + oarr = oarr.reshape(self._domain.shape) + res = Field.from_arr(self._domain, oarr) return res def _times(self, x): - if dobj.distaxis(x.val) in x.domain.axes[self._space]: - arr = x.to_global_data() - else: - arr = x.local_data + arr = x.val arr = arr.reshape(self._hshape) oarr = np.empty(self._pshape, dtype=x.dtype) oarr[()] = arr[(slice(None), self._dofdex, slice(None))] - return Field.from_local_data( - self._target, oarr.reshape(self._target.local_shape)) + return Field(self._target, oarr.reshape(self._target.shape)) def apply(self, x, mode): self._check_input(x, mode) diff --git a/nifty6/operators/field_zero_padder.py b/nifty6/operators/field_zero_padder.py index eabc7d2f0..afd3e47a7 100644 --- a/nifty6/operators/field_zero_padder.py +++ b/nifty6/operators/field_zero_padder.py @@ -17,7 +17,7 @@ import numpy as np -from .. import dobj, utilities +from .. import utilities from ..domain_tuple import DomainTuple from ..domains.rg_space import RGSpace from ..field import Field @@ -76,41 +76,38 @@ class FieldZeroPadder(LinearOperator): idx = (slice(None),) * d - v, x = dobj.ensure_not_distributed(v, (d,)) - if mode == self.TIMES: - shp = list(x.shape) + shp = list(v.shape) shp[d] = tgtshp[d] - xnew = np.zeros(shp, dtype=x.dtype) + xnew = np.zeros(shp, dtype=v.dtype) if self._central: - Nyquist = x.shape[d]//2 + Nyquist = v.shape[d]//2 i1 = idx + (slice(0, Nyquist+1),) - xnew[i1] = x[i1] + xnew[i1] = v[i1] i1 = idx + (slice(None, -(Nyquist+1), -1),) - xnew[i1] = x[i1] -# if (x.shape[d] & 1) == 0: # even number of pixels + xnew[i1] = v[i1] +# if (v.shape[d] & 1) == 0: # even number of pixels # i1 = idx+(Nyquist,) # xnew[i1] *= 0.5 # i1 = idx+(-Nyquist,) # xnew[i1] *= 0.5 else: - xnew[idx + (slice(0, x.shape[d]),)] = x + xnew[idx + (slice(0, v.shape[d]),)] = x else: # ADJOINT_TIMES if self._central: shp = list(x.shape) shp[d] = tgtshp[d] - xnew = np.zeros(shp, dtype=x.dtype) + xnew = np.zeros(shp, dtype=v.dtype) Nyquist = xnew.shape[d]//2 i1 = idx + (slice(0, Nyquist+1),) - xnew[i1] = x[i1] + xnew[i1] = v[i1] i1 = idx + (slice(None, -(Nyquist+1), -1),) - xnew[i1] += x[i1] + xnew[i1] += v[i1] # if (xnew.shape[d] & 1) == 0: # even number of pixels # i1 = idx+(Nyquist,) # xnew[i1] *= 0.5 else: - xnew = x[idx + (slice(0, tgtshp[d]),)] + xnew = v[idx + (slice(0, tgtshp[d]),)] curshp[d] = xnew.shape[d] - v = dobj.from_local_data(curshp, xnew, distaxis=dobj.distaxis(v)) - return Field(self._tgt(mode), dobj.ensure_default_distributed(v)) + return Field(self._tgt(mode), xnew) diff --git a/nifty6/operators/harmonic_operators.py b/nifty6/operators/harmonic_operators.py index 95d34ac8e..bd4a5e337 100644 --- a/nifty6/operators/harmonic_operators.py +++ b/nifty6/operators/harmonic_operators.py @@ -17,7 +17,7 @@ import numpy as np -from .. import dobj, utilities, fft +from .. import utilities, fft from ..domain_tuple import DomainTuple from ..domains.gl_space import GLSpace from ..domains.lm_space import LMSpace @@ -81,37 +81,7 @@ class FFTOperator(LinearOperator): fct = ncells axes = x.domain.axes[self._space] tdom = self._tgt(mode) - oldax = dobj.distaxis(x.val) - if oldax not in axes: # straightforward, no redistribution needed - ldat = x.local_data - ldat = func(ldat, axes=axes) - tmp = dobj.from_local_data(x.val.shape, ldat, distaxis=oldax) - elif len(axes) < len(x.shape) or len(axes) == 1: - # we can use one FFT pass in between the redistributions - tmp = dobj.redistribute(x.val, nodist=axes) - newax = dobj.distaxis(tmp) - ldat = dobj.local_data(tmp) - ldat = func(ldat, axes=axes) - tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=newax) - tmp = dobj.redistribute(tmp, dist=oldax) - else: # two separate FFTs needed - rem_axes = tuple(i for i in axes if i != oldax) - tmp = x.val - ldat = dobj.local_data(tmp) - ldat = func(ldat, axes=rem_axes) - if oldax != 0: - raise ValueError("bad distribution") - ldat2 = ldat.reshape((ldat.shape[0], - np.prod(ldat.shape[1:]))) - shp2d = (x.val.shape[0], np.prod(x.val.shape[1:])) - tmp = dobj.from_local_data(shp2d, ldat2, distaxis=0) - tmp = dobj.transpose(tmp) - ldat2 = dobj.local_data(tmp) - ldat2 = func(ldat2, axes=(1,)) - tmp = dobj.from_local_data(tmp.shape, ldat2, distaxis=0) - tmp = dobj.transpose(tmp) - ldat2 = dobj.local_data(tmp).reshape(ldat.shape) - tmp = dobj.from_local_data(x.val.shape, ldat2, distaxis=0) + tmp = func(x.val, axes=axes) Tval = Field(tdom, tmp) if mode & (LinearOperator.TIMES | LinearOperator.ADJOINT_TIMES): fct *= self._domain[self._space].scalar_dvol @@ -178,45 +148,7 @@ class HartleyOperator(LinearOperator): def _apply_cartesian(self, x, mode): axes = x.domain.axes[self._space] tdom = self._tgt(mode) - oldax = dobj.distaxis(x.val) - if oldax not in axes: # straightforward, no redistribution needed - ldat = x.local_data - ldat = fft.hartley(ldat, axes=axes) - tmp = dobj.from_local_data(x.val.shape, ldat, distaxis=oldax) - elif len(axes) < len(x.shape) or len(axes) == 1: - # we can use one Hartley pass in between the redistributions - tmp = dobj.redistribute(x.val, nodist=axes) - newax = dobj.distaxis(tmp) - ldat = dobj.local_data(tmp) - ldat = fft.hartley(ldat, axes=axes) - tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=newax) - tmp = dobj.redistribute(tmp, dist=oldax) - else: # two separate, full FFTs needed - # ideal strategy for the moment would be: - # - do real-to-complex FFT on all local axes - # - fill up array - # - redistribute array - # - do complex-to-complex FFT on remaining axis - # - add re+im - # - redistribute back - rem_axes = tuple(i for i in axes if i != oldax) - tmp = x.val - ldat = dobj.local_data(tmp) - ldat = fft.fftn(ldat, axes=rem_axes) - if oldax != 0: - raise ValueError("bad distribution") - ldat2 = ldat.reshape((ldat.shape[0], - np.prod(ldat.shape[1:]))) - shp2d = (x.val.shape[0], np.prod(x.val.shape[1:])) - tmp = dobj.from_local_data(shp2d, ldat2, distaxis=0) - tmp = dobj.transpose(tmp) - ldat2 = dobj.local_data(tmp) - ldat2 = fft.fftn(ldat2, axes=(1,)) - ldat2 = ldat2.real+ldat2.imag - tmp = dobj.from_local_data(tmp.shape, ldat2, distaxis=0) - tmp = dobj.transpose(tmp) - ldat2 = dobj.local_data(tmp).reshape(ldat.shape) - tmp = dobj.from_local_data(x.val.shape, ldat2, distaxis=0) + tmp = fft.hartley(x.val, axes=axes) Tval = Field(tdom, tmp) if mode & (LinearOperator.TIMES | LinearOperator.ADJOINT_TIMES): fct = self._domain[self._space].scalar_dvol @@ -318,18 +250,14 @@ class SHTOperator(LinearOperator): axes = x.domain.axes[self._space] axis = axes[0] v = x.val - v, idat = dobj.ensure_not_distributed(v, (axis,)) - distaxis = dobj.distaxis(v) p2h = not x.domain[self._space].harmonic tdom = self._tgt(mode) func = self._slice_p2h if p2h else self._slice_h2p - odat = np.empty(dobj.local_shape(tdom.shape, distaxis=distaxis), - dtype=x.dtype) - for slice in utilities.get_slice_list(idat.shape, axes): - odat[slice] = func(idat[slice]) - odat = dobj.from_local_data(tdom.shape, odat, distaxis) - return Field(tdom, dobj.ensure_default_distributed(odat)) + odat = np.empty(tdom.shape, dtype=x.dtype) + for slice in utilities.get_slice_list(v.shape, axes): + odat[slice] = func(v[slice]) + return Field(tdom, odat) def _unpickleSHTOperator(*args): diff --git a/nifty6/operators/outer_product_operator.py b/nifty6/operators/outer_product_operator.py index 21412c4c3..89c311155 100644 --- a/nifty6/operators/outer_product_operator.py +++ b/nifty6/operators/outer_product_operator.py @@ -41,10 +41,9 @@ class OuterProduct(LinearOperator): def apply(self, x, mode): self._check_input(x, mode) if mode == self.TIMES: - return Field.from_global_data( + return Field( self._target, np.multiply.outer( - self._field.to_global_data(), x.to_global_data())) + self._field.val, x.val)) axes = len(self._field.shape) - return Field.from_global_data( - self._domain, np.tensordot( - self._field.to_global_data(), x.to_global_data(), axes)) + return Field( + self._domain, np.tensordot(self._field.val, x.val, axes)) diff --git a/nifty6/operators/regridding_operator.py b/nifty6/operators/regridding_operator.py index 84ef00e59..404882cbe 100644 --- a/nifty6/operators/regridding_operator.py +++ b/nifty6/operators/regridding_operator.py @@ -17,7 +17,6 @@ import numpy as np -from .. import dobj from ..domain_tuple import DomainTuple from ..domains.rg_space import RGSpace from ..field import Field @@ -81,18 +80,15 @@ class RegriddingOperator(LinearOperator): idx = (slice(None),) * d wgt = self._frac[d-d0].reshape((1,)*d + (-1,) + (1,)*(ndim-d-1)) - v, x = dobj.ensure_not_distributed(v, (d,)) - if mode == self.ADJOINT_TIMES: - shp = list(x.shape) + shp = list(v.shape) shp[d] = tgtshp[d] - xnew = np.zeros(shp, dtype=x.dtype) - xnew = special_add_at(xnew, d, self._bindex[d-d0], x*(1.-wgt)) - xnew = special_add_at(xnew, d, self._bindex[d-d0]+1, x*wgt) + xnew = np.zeros(shp, dtype=v.dtype) + xnew = special_add_at(xnew, d, self._bindex[d-d0], v*(1.-wgt)) + xnew = special_add_at(xnew, d, self._bindex[d-d0]+1, v*wgt) else: # TIMES - xnew = x[idx + (self._bindex[d-d0],)] * (1.-wgt) - xnew += x[idx + (self._bindex[d-d0]+1,)] * wgt + xnew = v[idx + (self._bindex[d-d0],)] * (1.-wgt) + xnew += v[idx + (self._bindex[d-d0]+1,)] * wgt curshp[d] = xnew.shape[d] - v = dobj.from_local_data(curshp, xnew, distaxis=dobj.distaxis(v)) - return Field(self._tgt(mode), dobj.ensure_default_distributed(v)) + return Field(self._tgt(mode), xnew) diff --git a/nifty6/plot.py b/nifty6/plot.py index c4f547704..11921756e 100644 --- a/nifty6/plot.py +++ b/nifty6/plot.py @@ -19,7 +19,6 @@ import os import numpy as np -from . import dobj from .domains.gl_space import GLSpace from .domains.hp_space import HPSpace from .domains.power_space import PowerSpace @@ -172,9 +171,6 @@ def _find_closest(A, target): def _makeplot(name, block=True, dpi=None): import matplotlib.pyplot as plt - if dobj.rank != 0: - plt.close() - return if name is None: plt.show(block=block) if block: diff --git a/nifty6/sugar.py b/nifty6/sugar.py index c60970bca..87b602647 100644 --- a/nifty6/sugar.py +++ b/nifty6/sugar.py @@ -20,7 +20,7 @@ from time import time import numpy as np -from . import dobj, utilities +from . import utilities from .domain_tuple import DomainTuple from .domains.power_space import PowerSpace from .field import Field @@ -61,7 +61,7 @@ def PS_field(pspace, func): """ if not isinstance(pspace, PowerSpace): raise TypeError - data = dobj.from_global_data(func(pspace.k_lengths)) + data = func(pspace.k_lengths) return Field(DomainTuple.make(pspace), data) @@ -280,7 +280,7 @@ def from_random(random_type, domain, dtype=np.float64, **kwargs): return Field.from_random(random_type, domain, dtype, **kwargs) -def from_global_data(domain, arr, sum_up=False): +def from_global_data(domain, arr): """Convenience function creating Fields/MultiFields from Numpy arrays or dicts of Numpy arrays. @@ -290,11 +290,6 @@ def from_global_data(domain, arr, sum_up=False): the intended domain of the output field arr : Numpy array if `domain` corresponds to a `DomainTuple`, dictionary of Numpy arrays if `domain` corresponds to a `MultiDomain` - sum_up : bool - Only meaningful if MPI is enabled - If `True`, the contents of the arrays on all tasks are added together, - otherwise it is assumed that the array on each task holds the correct - field values. Returns ------- @@ -303,7 +298,7 @@ def from_global_data(domain, arr, sum_up=False): """ if isinstance(domain, (dict, MultiDomain)): return MultiField.from_global_data(domain, arr, sum_up) - return Field.from_global_data(domain, arr, sum_up) + return Field.from_arr(domain, arr) def from_local_data(domain, arr): @@ -324,7 +319,7 @@ def from_local_data(domain, arr): """ if isinstance(domain, (dict, MultiDomain)): return MultiField.from_local_data(domain, arr) - return Field.from_local_data(domain, arr) + return Field.from_arr(domain, arr) def makeDomain(domain): diff --git a/test/test_field.py b/test/test_field.py index 9f95a311f..54b830c62 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -28,7 +28,7 @@ SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES] @pmp('domain', SPACE_COMBINATIONS) @pmp('attribute_desired_type', - [['domain', ift.DomainTuple], ['val', ift.dobj.data_object], + [['domain', ift.DomainTuple], ['val', np.ndarray], ['shape', tuple], ['size', (np.int, np.int64)]]) def test_return_types(domain, attribute_desired_type): attribute = attribute_desired_type[0] @@ -72,8 +72,8 @@ def test_power_synthesize_analyze(space1, space2): sc1.add(sp.sum(spaces=1)/fp2.sum()) sc2.add(sp.sum(spaces=0)/fp1.sum()) - assert_allclose(sc1.mean.local_data, fp1.local_data, rtol=0.2) - assert_allclose(sc2.mean.local_data, fp2.local_data, rtol=0.2) + assert_allclose(sc1.mean.val, fp1.val, rtol=0.2) + assert_allclose(sc2.mean.val, fp2.val, rtol=0.2) @pmp('space1', [ @@ -101,8 +101,8 @@ def test_DiagonalOperator_power_analyze2(space1, space2): sc1.add(sp.sum(spaces=1)/fp2.sum()) sc2.add(sp.sum(spaces=0)/fp1.sum()) - assert_allclose(sc1.mean.local_data, fp1.local_data, rtol=0.2) - assert_allclose(sc2.mean.local_data, fp2.local_data, rtol=0.2) + assert_allclose(sc1.mean.val, fp1.val, rtol=0.2) + assert_allclose(sc2.mean.val, fp2.val, rtol=0.2) @pmp('space', [ @@ -112,7 +112,7 @@ def test_DiagonalOperator_power_analyze2(space1, space2): ]) def test_norm(space): f = ift.Field.from_random("normal", domain=space, dtype=np.complex128) - gd = f.to_global_data().reshape(-1) + gd = f.val.reshape(-1) assert_allclose(f.norm(), np.linalg.norm(gd)) assert_allclose(f.norm(1), np.linalg.norm(gd, ord=1)) assert_allclose(f.norm(2), np.linalg.norm(gd, ord=2)) @@ -133,7 +133,7 @@ def test_vdot2(): x2 = ift.RGSpace((150,)) m = ift.Field.full((x1, x2), .5) res = m.vdot(m, spaces=1) - assert_allclose(res.local_data, 37.5) + assert_allclose(res.val, 37.5) def test_outer(): @@ -142,7 +142,7 @@ def test_outer(): m1 = ift.Field.full(x1, .5) m2 = ift.Field.full(x2, 3.) res = m1.outer(m2) - assert_allclose(res.to_global_data(), np.full((9, 3), 1.5)) + assert_allclose(res.val, np.full((9, 3), 1.5)) def test_sum(): @@ -152,62 +152,62 @@ def test_sum(): 2, 12, ), distances=(0.3,)) - m1 = ift.Field.from_global_data(ift.makeDomain(x1), np.arange(9)) + m1 = ift.Field(ift.makeDomain(x1), np.arange(9)) m2 = ift.Field.full(ift.makeDomain((x1, x2)), 0.45) res1 = m1.sum() res2 = m2.sum(spaces=1) assert_allclose(res1, 36) - assert_allclose(res2.to_global_data(), np.full(9, 2*12*0.45)) + assert_allclose(res2.val, np.full(9, 2*12*0.45)) def test_integrate(): x1 = ift.RGSpace((9,), distances=2.) x2 = ift.RGSpace((2, 12), distances=(0.3,)) - m1 = ift.Field.from_global_data(ift.makeDomain(x1), np.arange(9)) + m1 = ift.Field(ift.makeDomain(x1), np.arange(9)) m2 = ift.Field.full(ift.makeDomain((x1, x2)), 0.45) res1 = m1.integrate() res2 = m2.integrate(spaces=1) assert_allclose(res1, 36*2) - assert_allclose(res2.to_global_data(), np.full(9, 2*12*0.45*0.3**2)) + assert_allclose(res2.val, np.full(9, 2*12*0.45*0.3**2)) def test_dataconv(): s1 = ift.RGSpace((10,)) - ld = np.arange(ift.dobj.local_shape(s1.shape)[0]) + ld = np.arange(s1.shape[0]) gd = np.arange(s1.shape[0]) - assert_equal(ld, ift.from_local_data(s1, ld).local_data) - assert_equal(gd, ift.from_global_data(s1, gd).to_global_data()) + assert_equal(ld, ift.from_local_data(s1, ld).val) + assert_equal(gd, ift.from_global_data(s1, gd).val) def test_cast_domain(): s1 = ift.RGSpace((10,)) s2 = ift.RGSpace((10,), distances=20.) d = np.arange(s1.shape[0]) - d2 = ift.from_global_data(s1, d).cast_domain(s2).to_global_data() + d2 = ift.from_global_data(s1, d).cast_domain(s2).val assert_equal(d, d2) def test_empty_domain(): f = ift.Field.full((), 5) - assert_equal(f.local_data, 5) + assert_equal(f.val, 5) f = ift.Field.full(None, 5) - assert_equal(f.local_data, 5) + assert_equal(f.val, 5) def test_trivialities(): s1 = ift.RGSpace((10,)) f1 = ift.Field.full(s1, 27) - assert_equal(f1.clip(min=29).local_data, 29.) - assert_equal(f1.clip(max=25).local_data, 25.) - assert_equal(f1.local_data, f1.real.local_data) - assert_equal(f1.local_data, (+f1).local_data) + assert_equal(f1.clip(min=29).val, 29.) + assert_equal(f1.clip(max=25).val, 25.) + assert_equal(f1.val, f1.real.val) + assert_equal(f1.val, (+f1).val) f1 = ift.Field.full(s1, 27. + 3j) - assert_equal(f1.one_over().local_data, (1./f1).local_data) - assert_equal(f1.real.local_data, 27.) - assert_equal(f1.imag.local_data, 3.) + assert_equal(f1.one_over().val, (1./f1).val) + assert_equal(f1.real.val, 27.) + assert_equal(f1.imag.val, 3.) assert_equal(f1.sum(), f1.sum(0)) - assert_equal(f1.conjugate().local_data, - ift.Field.full(s1, 27. - 3j).local_data) + assert_equal(f1.conjugate().val, + ift.Field.full(s1, 27. - 3j).val) f1 = ift.from_global_data(s1, np.arange(10)) # assert_equal(f1.min(), 0) # assert_equal(f1.max(), 9) @@ -218,7 +218,7 @@ def test_weight(): s1 = ift.RGSpace((10,)) f = ift.Field.full(s1, 10.) f2 = f.weight(1) - assert_equal(f.weight(1).local_data, f2.local_data) + assert_equal(f.weight(1).val, f2.val) assert_equal(f.domain.total_volume(), 1) assert_equal(f.domain.total_volume(0), 1) assert_equal(f.domain.total_volume((0,)), 1) @@ -291,43 +291,43 @@ def test_err(): def test_stdfunc(): s = ift.RGSpace((200,)) f = ift.Field.full(s, 27) - assert_equal(f.local_data, 27) + assert_equal(f.val, 27) assert_equal(f.shape, (200,)) assert_equal(f.dtype, np.int) fx = ift.full(f.domain, 0) assert_equal(f.dtype, fx.dtype) assert_equal(f.shape, fx.shape) - assert_equal(fx.local_data, 0) + assert_equal(fx.val, 0) fx = ift.full(f.domain, 1) assert_equal(f.dtype, fx.dtype) assert_equal(f.shape, fx.shape) - assert_equal(fx.local_data, 1) + assert_equal(fx.val, 1) fx = ift.full(f.domain, 67.) assert_equal(f.shape, fx.shape) - assert_equal(fx.local_data, 67.) + assert_equal(fx.val, 67.) f = ift.Field.from_random("normal", s) f2 = ift.Field.from_random("normal", s) - assert_equal((f > f2).local_data, f.local_data > f2.local_data) - assert_equal((f >= f2).local_data, f.local_data >= f2.local_data) - assert_equal((f < f2).local_data, f.local_data < f2.local_data) - assert_equal((f <= f2).local_data, f.local_data <= f2.local_data) - assert_equal((f != f2).local_data, f.local_data != f2.local_data) - assert_equal((f == f2).local_data, f.local_data == f2.local_data) - assert_equal((f + f2).local_data, f.local_data + f2.local_data) - assert_equal((f - f2).local_data, f.local_data - f2.local_data) - assert_equal((f*f2).local_data, f.local_data*f2.local_data) - assert_equal((f/f2).local_data, f.local_data/f2.local_data) - assert_equal((-f).local_data, -(f.local_data)) - assert_equal(abs(f).local_data, abs(f.local_data)) + assert_equal((f > f2).val, f.val > f2.val) + assert_equal((f >= f2).val, f.val >= f2.val) + assert_equal((f < f2).val, f.val < f2.val) + assert_equal((f <= f2).val, f.val <= f2.val) + assert_equal((f != f2).val, f.val != f2.val) + assert_equal((f == f2).val, f.val == f2.val) + assert_equal((f + f2).val, f.val + f2.val) + assert_equal((f - f2).val, f.val - f2.val) + assert_equal((f*f2).val, f.val*f2.val) + assert_equal((f/f2).val, f.val/f2.val) + assert_equal((-f).val, -(f.val)) + assert_equal(abs(f).val, abs(f.val)) def test_emptydomain(): f = ift.Field.full((), 3.) assert_equal(f.sum(), 3.) assert_equal(f.prod(), 3.) - assert_equal(f.local_data, 3.) - assert_equal(f.local_data.shape, ()) - assert_equal(f.local_data.size, 1) + assert_equal(f.val, 3.) + assert_equal(f.val.shape, ()) + assert_equal(f.val.size, 1) assert_equal(f.vdot(f), 9.) @@ -342,7 +342,7 @@ def test_funcs(num, dom, func): f = ift.Field.full(dom, num) res = getattr(f, func)() res2 = getattr(np, func)(num) - assert_allclose(res.local_data, res2) + assert_allclose(res.val, res2) @pmp('rtype', ['normal', 'pm1', 'uniform']) @@ -356,4 +356,4 @@ def test_field_of_objects(): arr = np.array(['x', 'y', 'z']) sp = ift.RGSpace(3) with assert_raises(TypeError): - ift.Field.from_global_data(sp, arr) + ift.Field(sp, arr) -- GitLab