Skip to content

GitLab

Commit 144b6fef authored by theos's avatar theos
Browse files

Huge step towards field on multiple spaces. 

Added functionality for multiple field_types per field.
parent d6800a9b
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 600 additions and 481 deletions
demos/demo_wf1.py
View file @ 144b6fef
......@@ -93,7 +93,7 @@ if __name__ == "__main__":
#temp_result = (D.inverse_times(m)-xi)
n_power = x_space.enforce_power(s.var()/x_space.get_dim())
n_power = x_space.enforce_power(s.var()/x_space.dim)
s_power = S.get_power()
s.plot(title="signal", save = 'plot_s.png')
......
field_types/base_field_type.py
View file @ 144b6fef
......@@ -21,7 +21,15 @@ class Field_type(object):
def dtype(self):
return self._dtype
def get_dof(self, split=False):
@property
def dof(self):
return self._dof
@property
def dof_split(self):
return self._dof_split
def _get_dof(self, split=False):
if issubclass(self.dtype.type, np.complexfloating):
multiplicator = 2
else:
......@@ -49,3 +57,6 @@ class Field_type(object):
def complement_cast(self, x, axis=None):
return x
def dot_contraction(self, x, axes):
raise NotImplementedError
field_types/field_array.py
View file @ 144b6fef
......@@ -5,4 +5,5 @@ from base_field_type import Field_type
class Field_array(Field_type):
pass
def dot_contraction(self, x, axes):
return x.sum(axis=axes)
lm/nifty_lm.py
View file @ 144b6fef
......@@ -173,7 +173,7 @@ class lm_space(point_space):
self.power_indices = lm_power_indices(
lmax=self.paradict['lmax'],
dim=self.get_dim(),
dim=self.dim,
allowed_distribution_strategies=LM_DISTRIBUTION_STRATEGIES)
@property
......@@ -210,12 +210,14 @@ class lm_space(point_space):
mmax=self.paradict['mmax'],
dtype=self.dtype)
def get_shape(self):
@property
def shape(self):
lmax = self.paradict['lmax']
mmax = self.paradict['mmax']
return (np.int((mmax + 1) * (lmax + 1) - ((mmax + 1) * mmax) // 2),)
def get_dof(self, split=False):
@property
def dof(self):
"""
Computes the number of degrees of freedom of the space, taking into
account symmetry constraints and complex-valuedness.
......@@ -234,12 +236,14 @@ class lm_space(point_space):
lmax = self.paradict['lmax']
mmax = self.paradict['mmax']
dof = np.int((lmax + 1) * (2 * mmax + 1) - (mmax + 1) * mmax)
if split:
return (dof, )
else:
return dof
return dof
@property
def dof_split(self):
return (self.dof,)
def get_meta_volume(self, split=False):
@property
def meta_volume(self):
"""
Calculates the meta volumes.
......@@ -265,12 +269,13 @@ class lm_space(point_space):
volume 1, the ones with :math:`m>0` have meta volume 2, sinnce they
each determine another component with negative :math:`m`.
"""
if not split:
return np.float(self.get_dof())
else:
mol = self.cast(1, dtype=np.float)
mol[self.paradict['lmax'] + 1:] = 2 # redundant: (l,m) and (l,-m)
return mol
return np.float(self.dof())
@property
def meta_volume_split(self):
mol = self.cast(1, dtype=np.float)
mol[self.paradict['lmax'] + 1:] = 2 # redundant: (l,m) and (l,-m)
return mol
def _complement_cast(self, x, axis=None, **kwargs):
if axis is None:
......@@ -441,11 +446,11 @@ class lm_space(point_space):
# x = 0
#
# elif arg['random'] == "pm1":
# x = random.pm1(dtype=self.dtype, shape=self.get_shape())
# x = random.pm1(dtype=self.dtype, shape=self.shape)
#
# elif arg['random'] == "gau":
# x = random.gau(dtype=self.dtype,
# shape=self.get_shape(),
# shape=self.shape,
# mean=arg['mean'],
# std=arg['std'])
......@@ -471,7 +476,7 @@ class lm_space(point_space):
# elif arg['random'] == "uni":
# x = random.uni(dtype=self.dtype,
# shape=self.get_shape(),
# shape=self.shape,
# vmin=arg['vmin'],
# vmax=arg['vmax'])
#
......@@ -481,37 +486,51 @@ class lm_space(point_space):
sample = self.cast(sample)
return sample
def calc_dot(self, x, y):
"""
Computes the discrete inner product of two given arrays of field
values.
Parameters
----------
x : numpy.ndarray
First array
y : numpy.ndarray
Second array
Returns
-------
dot : scalar
Inner product of the two arrays.
"""
x = self.cast(x)
y = self.cast(y)
# def calc_dot(self, x, y):
# """
# Computes the discrete inner product of two given arrays of field
# values.
#
# Parameters
# ----------
# x : numpy.ndarray
# First array
# y : numpy.ndarray
# Second array
#
# Returns
# -------
# dot : scalar
# Inner product of the two arrays.
# """
# x = self.cast(x)
# y = self.cast(y)
#
# lmax = self.paradict['lmax']
#
# x_low = x[:lmax + 1]
# x_high = x[lmax + 1:]
# y_low = y[:lmax + 1]
# y_high = y[lmax + 1:]
#
# dot = (x_low.real * y_low.real).sum()
# dot += 2 * (x_high.real * y_high.real).sum()
# dot += 2 * (x_high.imag * y_high.imag).sum()
# return dot
def dot_contraction(self, x, axes):
assert len(axes) == 1
axis = axes[0]
lmax = self.paradict['lmax']
x_low = x[:lmax + 1]
x_high = x[lmax + 1:]
y_low = y[:lmax + 1]
y_high = y[lmax + 1:]
# extract the low and high parts of x
extractor = ()
extractor += (slice(None),)*axis
low_extractor = extractor + (slice(None, lmax+1), )
high_extractor = extractor + (slice(lmax+1), )
dot = (x_low.real * y_low.real).sum()
dot += 2 * (x_high.real * y_high.real).sum()
dot += 2 * (x_high.imag * y_high.imag).sum()
return dot
result = x[low_extractor].sum(axes) + 2 * x[high_extractor].sum(axes)
return result
def calc_transform(self, x, codomain=None, **kwargs):
"""
......@@ -884,7 +903,7 @@ class lm_space(point_space):
fig.canvas.draw()
def getlm(self): # > compute all (l,m)
index = np.arange(self.get_dim())
index = np.arange(self.dim)
n = 2 * self.paradict['lmax'] + 1
m = np.ceil(
(n - np.sqrt(n**2 - 8 * (index - self.paradict['lmax']))) / 2
......@@ -1008,29 +1027,12 @@ class gl_space(point_space):
nlon=self.paradict['nlon'],
dtype=self.dtype)
def get_shape(self):
@property
def shape(self):
return (np.int((self.paradict['nlat'] * self.paradict['nlon'])),)
def get_dof(self, split=False):
"""
Computes the number of degrees of freedom of the space.
Returns
-------
dof : int
Number of degrees of freedom of the space.
Notes
-----
Since the :py:class:`gl_space` class only supports real-valued
fields, the number of degrees of freedom is the number of pixels.
"""
if split:
return self.get_shape()
else:
return self.get_dim()
def get_meta_volume(self, split=False):
@property
def meta_volume(self):
"""
Calculates the meta volumes.
......@@ -1055,11 +1057,12 @@ class gl_space(point_space):
For Gauss-Legendre pixelizations, the meta volumes are the pixel
sizes.
"""
if not split:
return np.float(4 * np.pi)
else:
mol = self.cast(1, dtype=np.float)
return self.calc_weight(mol, power=1)
return np.float(4 * np.pi)
@property
def meta_volume_split(self):
mol = self.cast(1, dtype=np.float)
return self.calc_weight(mol, power=1)
# TODO: Extend to binning/log
def enforce_power(self, spec, size=None, kindex=None):
......@@ -1171,14 +1174,14 @@ class gl_space(point_space):
arg = random.parse_arguments(self, **kwargs)
# if(arg is None):
# x = np.zeros(self.get_shape(), dtype=self.dtype)
# x = np.zeros(self.shape, dtype=self.dtype)
#
# elif(arg['random'] == "pm1"):
# x = random.pm1(dtype=self.dtype, shape=self.get_shape())
# x = random.pm1(dtype=self.dtype, shape=self.shape)
#
# elif(arg['random'] == "gau"):
# x = random.gau(dtype=self.dtype,
# shape=self.get_shape(),
# shape=self.shape,
# mean=arg['mean'],
# std=arg['std'])
#
......@@ -1204,7 +1207,7 @@ class gl_space(point_space):
# elif(arg['random'] == "uni"):
# x = random.uni(dtype=self.dtype,
# shape=self.get_shape(),
# shape=self.shape,
# vmin=arg['vmin'],
# vmax=arg['vmax'])
#
......@@ -1667,29 +1670,12 @@ class hp_space(point_space):
def copy(self):
return hp_space(nside=self.paradict['nside'])
def get_shape(self):
@property
def shape(self):
return (np.int(12 * self.paradict['nside']**2),)
def get_dof(self, split=False):
"""
Computes the number of degrees of freedom of the space.
Returns
-------
dof : int
Number of degrees of freedom of the space.
Notes
-----
Since the :py:class:`hp_space` class only supports real-valued
fields, the number of degrees of freedom is the number of pixels.
"""
if split:
return self.get_shape()
else:
return self.get_dim()
def get_meta_volume(self, split=False):
@property
def meta_volume(self):
"""
Calculates the meta volumes.
......@@ -1713,11 +1699,12 @@ class hp_space(point_space):
-----
For HEALpix discretizations, the meta volumes are the pixel sizes.
"""
if not split:
return np.float(4 * np.pi)
else:
mol = self.cast(1, dtype=np.float)
return self.calc_weight(mol, power=1)
return np.float(4 * np.pi)
@property
def meta_volume_split(self):
mol = self.cast(1, dtype=np.float)
return self.calc_weight(mol, power=1)
# TODO: Extend to binning/log
def enforce_power(self, spec, size=None, kindex=None):
......@@ -1822,13 +1809,13 @@ class hp_space(point_space):
arg = random.parse_arguments(self, **kwargs)
# if arg is None:
# x = np.zeros(self.get_shape(), dtype=self.dtype)
# x = np.zeros(self.shape, dtype=self.dtype)
#
# elif arg['random'] == "pm1":
# x = random.pm1(dtype=self.dtype, shape=self.get_shape())
# x = random.pm1(dtype=self.dtype, shape=self.shape)
#
# elif arg['random'] == "gau":
# x = random.gau(dtype=self.dtype, shape=self.get_shape(),
# x = random.gau(dtype=self.dtype, shape=self.shape,
# mean=arg['mean'],
# std=arg['std'])
......@@ -1849,7 +1836,7 @@ class hp_space(point_space):
# elif arg['random'] == "uni":
# x = random.uni(dtype=self.dtype, shape=self.get_shape(),
# x = random.uni(dtype=self.dtype, shape=self.shape,
# vmin=arg['vmin'],
# vmax=arg['vmax'])
#
......
nifty_core.py
View file @ 144b6fef
......@@ -259,7 +259,7 @@ class space(object):
return not self.__eq__(x)
def __len__(self):
return int(self.get_dim())
return int(self.dim)
def copy(self):
return space(para=self.para,
......@@ -277,11 +277,13 @@ class space(object):
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'apply_scalar_function'."))
def get_shape(self):
@property
def shape(self):
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'shape'."))
def get_dim(self):
@property
def dim(self):
"""
Computes the dimension of the space, i.e.\ the number of pixels.
......@@ -299,7 +301,8 @@ class space(object):
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'dim'."))
def get_dof(self):
@property
def dof(self):
"""
Computes the number of degrees of freedom of the space.
......@@ -311,6 +314,19 @@ class space(object):
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'dof'."))
@property
def dof_split(self):
"""
Computes the number of degrees of freedom of the space.
Returns
-------
dof : int
Number of degrees of freedom of the space.
"""
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'dof_split'."))
def complement_cast(self, x, axis=None):
return x
......@@ -501,7 +517,7 @@ class space(object):
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'norm'."))
def calc_dot(self, x, y):
def dot_contraction(self, x, axes):
"""
Computes the discrete inner product of two given arrays of field
values.
......@@ -828,10 +844,12 @@ class point_space(space):
def apply_scalar_function(self, x, function, inplace=False):
return x.apply_scalar_function(function, inplace=inplace)
def get_shape(self):
@property
def shape(self):
return (self.paradict['num'],)
def get_dim(self):
@property
def dim(self):
"""
Computes the dimension of the space, i.e.\ the number of points.
......@@ -846,9 +864,10 @@ class point_space(space):
dim : {int, numpy.ndarray}
Dimension(s) of the space.
"""
return np.prod(self.get_shape())
return np.prod(self.shape)
def get_dof(self, split=False):
@property
def dof(self):
"""
Computes the number of degrees of freedom of the space, i.e./ the
number of points for real-valued fields and twice that number for
......@@ -859,23 +878,38 @@ class point_space(space):
dof : int
Number of degrees of freedom of the space.
"""
if split:
dof = self.get_shape()
if issubclass(self.dtype.type, np.complexfloating):
dof = tuple(np.array(dof)*2)
else:
dof = self.get_dim()
if issubclass(self.dtype.type, np.complexfloating):
dof = dof * 2
dof = self.dim
if issubclass(self.dtype.type, np.complexfloating):
dof = dof * 2
return dof
def get_vol(self, split=False):
if split:
return self.distances
else:
return np.prod(self.distances)
@property
def dof_split(self):
"""
Computes the number of degrees of freedom of the space, i.e./ the
number of points for real-valued fields and twice that number for
complex-valued fields.
def get_meta_volume(self, split=False):
Returns
-------
dof : int
Number of degrees of freedom of the space.
"""
dof = self.shape
if issubclass(self.dtype.type, np.complexfloating):
dof = tuple(np.array(dof)*2)
return dof
@property
def vol(self, split=False):
return np.prod(self.distances)
@property
def vol_split(self):
return self.distances
@property
def meta_volume(self):
"""
Calculates the meta volumes.
......@@ -896,11 +930,12 @@ class point_space(space):
mol : {numpy.ndarray, float}
Meta volume of the pixels or the complete space.
"""
if not split:
return self.get_dim() * self.get_vol()
else:
mol = self.cast(1, dtype=np.dtype('float'))
return self.calc_weight(mol, power=1)
return self.dim() * self.vol()
@property
def meta_volume_split(self):
mol = self.cast(1, dtype=np.dtype('float'))
return self.calc_weight(mol, power=1)
def enforce_power(self, spec, **kwargs):
"""
......@@ -1028,110 +1063,110 @@ class point_space(space):
"""
return self.copy()
def get_random_values(self, **kwargs):
"""
Generates random field values according to the specifications given
by the parameters.
Returns
-------
x : numpy.ndarray
Valid field values.
Other parameters
----------------
random : string, *optional*
Specifies the probability distribution from which the random
numbers are to be drawn.
Supported distributions are:
- "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
- "gau" (normal distribution with zero-mean and a given
standard
deviation or variance)
- "syn" (synthesizes from a given power spectrum)
- "uni" (uniform distribution over [vmin,vmax[)
(default: None).
dev : float, *optional*
Standard deviation (default: 1).
var : float, *optional*
Variance, overriding `dev` if both are specified
(default: 1).
spec : {scalar, list, numpy.ndarray, nifty.field, function},
*optional*
Power spectrum (default: 1).
pindex : numpy.ndarray, *optional*
Indexing array giving the power spectrum index of each band
(default: None).
kindex : numpy.ndarray, *optional*
Scale of each band (default: None).
codomain : nifty.space, *optional*
A compatible codomain with power indices (default: None).
log : bool, *optional*
Flag specifying if the spectral binning is performed on
logarithmic
scale or not; if set, the number of used bins is set
automatically (if not given otherwise); by default no binning
is done (default: None).
nbin : integer, *optional*
Number of used spectral bins; if given `log` is set to
``False``;
integers below the minimum of 3 induce an automatic setting;
by default no binning is done (default: None).
binbounds : {list, array}, *optional*
User specific inner boundaries of the bins, which are preferred
over the above parameters; by default no binning is done
(default: None).
vmin : {scalar, list, ndarray, field}, *optional*
Lower limit of the uniform distribution if ``random == "uni"``
(default: 0).
vmin : float, *optional*
Lower limit for a uniform distribution (default: 0).
vmax : float, *optional*
Upper limit for a uniform distribution (default: 1).
"""
arg = random.parse_arguments(self, **kwargs)
if arg is None:
return self.cast(0)
# Prepare the empty distributed_data_object
sample = distributed_data_object(
global_shape=self.get_shape(),
dtype=self.dtype)
# Case 1: uniform distribution over {-1,+1}/{1,i,-1,-i}
if arg['random'] == 'pm1':
sample.apply_generator(lambda s: random.pm1(dtype=self.dtype,
shape=s))
# Case 2: normal distribution with zero-mean and a given standard
# deviation or variance
elif arg['random'] == 'gau':
std = arg['std']
if np.isscalar(std) or std is None:
processed_std = std
else:
try:
processed_std = sample.distributor. \
extract_local_data(std)
except(AttributeError):
processed_std = std
sample.apply_generator(lambda s: random.gau(dtype=self.dtype,
shape=s,
mean=arg['mean'],
std=processed_std))
# Case 3: uniform distribution
elif arg['random'] == 'uni':
sample.apply_generator(lambda s: random.uni(dtype=self.dtype,
shape=s,