Skip to content

GitLab

Commit 42abacbd authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanups

parent 6e043763
Pipeline #21899 passed with stage
in 4 minutes and 19 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 85 additions and 126 deletions
nifty/data_objects/distributed_do.py
View file @ 42abacbd
......@@ -297,6 +297,9 @@ def from_object(object, dtype=None, copy=True):
# 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)
for i in range(ntask):
......@@ -387,8 +390,6 @@ def redistribute(arr, dist=None, nodist=None):
return from_global_data(out, distaxis=-1)
# real redistribution via Alltoallv
# temporary slow, but simple solution for comparison purposes:
# return redistribute(redistribute(arr,dist=-1),dist=dist)
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)
......
nifty/data_objects/random.py
View file @ 42abacbd
......@@ -24,38 +24,31 @@ class Random(object):
@staticmethod
def pm1(dtype, shape):
if issubclass(dtype, (complex, np.complexfloating)):
x = np.array([1 + 0j, 0 + 1j, -1 + 0j, 0 - 1j], dtype=dtype)
x = np.array([1+0j, 0+1j, -1+0j, 0-1j], dtype=dtype)
x = x[np.random.randint(4, size=shape)]
else:
x = 2 * np.random.randint(2, size=shape) - 1
return x.astype(dtype)
x = 2*np.random.randint(2, size=shape) - 1
return x.astype(dtype, copy=False)
@staticmethod
def normal(dtype, shape, mean=0., std=1.):
if issubclass(dtype, (complex, np.complexfloating)):
x = np.empty(shape, dtype=dtype)
x.real = np.random.normal(loc=mean.real, scale=std*np.sqrt(0.5),
size=shape)
x.imag = np.random.normal(loc=mean.imag, scale=std*np.sqrt(0.5),
size=shape)
x.real = np.random.normal(mean.real, std*np.sqrt(0.5), shape)
x.imag = np.random.normal(mean.imag, std*np.sqrt(0.5), shape)
else:
x = np.random.normal(loc=mean, scale=std, size=shape)
x = x.astype(dtype, copy=False)
x = np.random.normal(mean, std, shape).astype(dtype, copy=False)
return x
@staticmethod
def uniform(dtype, shape, low=0., high=1.):
if issubclass(dtype, (complex, np.complexfloating)):
x = np.empty(shape, dtype=dtype)
x.real = (high - low) * np.random.random(shape) + low
x.imag = (high - low) * np.random.random(shape) + low
x.real = np.random.uniform(low, high, shape)
x.imag = np.random.uniform(low, high, shape)
elif dtype in [np.dtype('int8'), np.dtype('int16'), np.dtype('int32'),
np.dtype('int64')]:
x = np.random.random_integers(min(low, high),
high=max(low, high),
size=shape)
x = np.random.random.randint(low, high+1, shape)
else:
x = (high - low) * np.random.random(shape) + low
x = np.random.uniform(low, high, shape)
return x.astype(dtype, copy=False)
nifty/field.py
View file @ 42abacbd
......@@ -36,15 +36,12 @@ class Field(object):
Parameters
----------
domain : DomainObject
One of the space types NIFTY supports. RGSpace, GLSpace, HPSpace,
LMSpace or PowerSpace. It might also be a FieldArray, which is
an unstructured domain.
domain : None, DomainTuple, tuple of DomainObjects, or single DomainObject
val : scalar, numpy.ndarray, Field
val : None, Field, data_object, or scalar
The values the array should contain after init. A scalar input will
fill the whole array with this scalar. If an array is provided the
array's dimensions must match the domain's.
fill the whole array with this scalar. If a data_object is provided,
its dimensions must match the domain's.
dtype : type
A numpy.type. Most common are float and complex.
......@@ -53,24 +50,14 @@ class Field(object):
Attributes
----------
val : numpy.ndarray
val : data_object
domain : DomainTuple
See Parameters.
dtype : type
Contains the datatype stored in the Field.
Raise
-----
TypeError
Raised if
*the given domain contains something that is not a DomainObject
instance
*val is an array that has a different dimension than the domain
"""
# ---Initialization methods---
def __init__(self, domain=None, val=None, dtype=None, copy=False):
self.domain = self._parse_domain(domain=domain, val=val)
......@@ -188,8 +175,6 @@ class Field(object):
def fill(self, fill_value):
self._val.fill(fill_value)
# ---Properties---
@property
def val(self):
""" Returns the data object associated with this Field.
......@@ -234,8 +219,6 @@ class Field(object):
""" The imaginary part of the field (data is not copied)."""
return Field(self.domain, self.val.imag)
# ---Special unary/binary operations---
def copy(self):
""" Returns a full copy of the Field.
......@@ -456,8 +439,6 @@ class Field(object):
raise ValueError("domains are incompatible.")
dobj.local_data(self.val)[()] = dobj.local_data(other.val)[()]
# ---General binary methods---
def _binary_helper(self, other, op):
# if other is a field, make sure that the domains match
if isinstance(other, Field):
......
nifty/operators/fft_operator.py
View file @ 42abacbd
......@@ -72,6 +72,7 @@ class FFTOperator(LinearOperator):
if "domain" or "target" are not of the proper type.
"""
# MR FIXME: target should only be a single DomainObject, not the full tuple
def __init__(self, domain, target=None, space=None):
super(FFTOperator, self).__init__()
......
nifty/sugar.py
View file @ 42abacbd
......@@ -39,7 +39,6 @@ def PS_field(pspace, func, dtype=None):
def _single_power_analyze(field, idx, binbounds):
from .operators.power_projection_operator import PowerProjectionOperator
power_domain = PowerSpace(field.domain[idx], binbounds)
ppo = PowerProjectionOperator(field.domain, power_domain, idx)
return ppo(field.weight(-1))
......@@ -76,11 +75,6 @@ def power_analyze(field, spaces=None, binbounds=None,
of power_analyze with keep_phase_information=False.
(default : False).
Raise
-----
TypeError
Raised if any of the input field's domains is not harmonic
Returns
-------
out : Field
......@@ -88,14 +82,11 @@ def power_analyze(field, spaces=None, binbounds=None,
the power spectrum of 'field'.
"""
# check if all spaces in `field.domain` are either harmonic or
# power_space instances
for sp in field.domain:
if not sp.harmonic and not isinstance(sp, PowerSpace):
dobj.mprint("WARNING: Field has a space in `domain` which is "
"neither harmonic nor a PowerSpace.")
# check if the `spaces` input is valid
if spaces is None:
spaces = range(len(field.domain))
else:
......@@ -120,15 +111,10 @@ def power_analyze(field, spaces=None, binbounds=None,
def _compute_spec(field, spaces):
from .operators.power_projection_operator import PowerProjectionOperator
if spaces is None:
spaces = range(len(field.domain))
else:
spaces = utilities.cast_iseq_to_tuple(spaces)
spaces = range(len(field.domain)) if spaces is None \
else utilities.cast_iseq_to_tuple(spaces)
# create the result domain
result_domain = list(field.domain)
spec = sqrt(field)
for i in spaces:
result_domain[i] = field.domain[i].harmonic_partner
......@@ -174,7 +160,6 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
Raises
------
ValueError : If domain specified by `spaces` is not a PowerSpace.
"""
spec = _compute_spec(field, spaces)
......@@ -187,12 +172,11 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
else np.complex)
for x in range(1 if real_power else 2)]
# MR: dummy call - will be removed soon
# MR FIXME: dummy call - will be removed soon
if real_signal:
field.from_random('normal', mean=0., std=1.,
domain=spec.domain, dtype=np.float)
# apply the rescaler to the random fields
result[0] *= spec.real
if not real_power:
result[1] *= spec.imag
......@@ -203,7 +187,7 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
def power_synthesize_special(field, spaces=None):
spec = _compute_spec(field, spaces)
# MR: dummy call - will be removed soon
# MR FIXME: dummy call - will be removed soon
field.from_random('normal', mean=0., std=1.,
domain=spec.domain, dtype=np.complex)
......@@ -223,8 +207,7 @@ def create_power_field(domain, power_spectrum, dtype=None):
else:
power_domain = PowerSpace(domain)
fp = PS_field(power_domain, power_spectrum, dtype)
P = PowerProjectionOperator(domain, power_domain)
f = P.adjoint_times(fp)
f = PowerProjectionOperator(domain, power_domain).adjoint_times(fp)
if not issubclass(fp.dtype.type, np.complexfloating):
f = f.real
......@@ -263,8 +246,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
space = 0
space = int(space)
return DiagonalOperator(
create_power_field(domain[space],
power_spectrum, dtype).weight(1),
create_power_field(domain[space], power_spectrum, dtype).weight(1),
domain=domain,
spaces=space)
......@@ -330,5 +312,6 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
interdomain[i] = codomain[i]
fft_op_list += [FFTOperator(domain=domain, target=interdomain,
space=i)]
# MR FIXME: this looks slightly fishy...
domain = interdomain
return ComposedOperator(fft_op_list)
nifty/utilities.py
View file @ 42abacbd
......@@ -124,34 +124,6 @@ class NiftyMeta(_DocStringInheritor, abc.ABCMeta):
pass
def _fill_upper_half(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
np.subtract(tmp[slice2].real, tmp[slice2].imag, out=res[slice1])
for i, ax in enumerate(axes[:-1]):
dim1 = [slice(None)]*ax + [slice(0, 1)]
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half(tmp[dim1], res[dim1], axes2)
def _fill_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
np.add(tmp.real, tmp.imag, out=res[slice1])
_fill_upper_half(tmp, res, axes)
return res
def hartley(a, axes=None):
# Check if the axes provided are valid given the shape
if axes is not None and \
......@@ -162,36 +134,35 @@ def hartley(a, axes=None):
from pyfftw.interfaces.numpy_fft import rfftn
tmp = rfftn(a, axes=axes)
return _fill_array(tmp, np.empty_like(a), axes)
def _fill_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
np.add(tmp.real, tmp.imag, out=res[slice1])
def _fill_upper_half(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
np.subtract(tmp[slice2].real, tmp[slice2].imag, out=res[slice1])
for i, ax in enumerate(axes[:-1]):
dim1 = [slice(None)]*ax + [slice(0, 1)]
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half(tmp[dim1], res[dim1], axes2)
_fill_upper_half(tmp, res, axes)
return res
def _fill_upper_half_complex(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
#np.conjugate(tmp[slice2], out=res[slice1])
res[slice1] = np.conjugate(tmp[slice2])
for i, ax in enumerate(axes[:-1]):
dim1 = [slice(None)]*ax + [slice(0, 1)]
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half_complex(tmp[dim1], res[dim1], axes2)
def _fill_complex_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
res[slice1] = tmp
_fill_upper_half_complex(tmp, res, axes)
return res
return _fill_array(tmp, np.empty_like(a), axes)
# Do a real-to-complex forward FFT and return the _full_ output array
......@@ -205,7 +176,36 @@ def my_fftn_r2c(a, axes=None):
from pyfftw.interfaces.numpy_fft import rfftn
tmp = rfftn(a, axes=axes)
return _fill_complex_array(tmp, np.empty_like(a,dtype=tmp.dtype), axes)
def _fill_complex_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
res[slice1] = tmp
def _fill_upper_half_complex(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
# np.conjugate(tmp[slice2], out=res[slice1])
res[slice1] = np.conjugate(tmp[slice2])
for i, ax in enumerate(axes[:-1]):
dim1 = [slice(None)]*ax + [slice(0, 1)]
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half_complex(tmp[dim1], res[dim1], axes2)
_fill_upper_half_complex(tmp, res, axes)
return res
return _fill_complex_array(tmp, np.empty_like(a, dtype=tmp.dtype), axes)
def general_axpy(a, x, y, out):
......
nifty/version.py
View file @ 42abacbd
......@@ -20,4 +20,4 @@
# 1) we don't load dependencies by storing it in __init__.py
# 2) we can import it in setup.py for the same reason
# 3) we can import it into your module module
__version__ = '3.1.0'
__version__ = '3.9.0'
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment