There is a maintenance of MPCDF Gitlab on Thursday, April 22st 2020, 9:00 am CEST - Expect some service interruptions during this time

Commit 30d072b9 authored by theos's avatar theos

Revised SmoothOperator implementation.

parent 3454501e
......@@ -207,7 +207,7 @@ class Field(object):
if not self.domain[space_index].harmonic:
raise ValueError(about._errors.cstring(
"ERROR: Conversion of only one space at a time is allowed."))
"ERROR: The analyzed space must be harmonic."))
# Create the target PowerSpace instance:
# If the associated signal-space field was real, we extract the
......
......@@ -53,45 +53,54 @@ class SmoothOperator(EndomorphicOperator):
return self._sigma
def _smooth_helper(self, x, spaces, types, inverse=False):
# copy for doing the actual smoothing
smooth_out = x.copy()
if spaces is not None and self.sigma != 0:
if self.sigma == 0:
return x.copy()
# the domain of the smoothing operator contains exactly one space.
# Hence, if spaces is None, but we passed LinearOperator's
# _check_input_compatibility, we know that x is also solely defined
# on that space
if spaces is None:
spaces = (0,)
else:
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
axes = x.domain_axes[spaces[0]]
Transformator = FFTOperator(x.domain[spaces[0]])
transform = FFTOperator(x.domain[spaces[0]])
# transform to the (global-)default codomain and perform all remaining
# steps therein
transformed_x = Transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]]
# transform
smooth_out = transform(smooth_out, spaces=spaces[0])
# create the kernel using the knowledge of codomain about itself
axes_local_distribution_strategy = \
transformed_x.val.get_axes_local_distribution_strategy(axes=coaxes)
# create the kernel
space_obj = smooth_out.domain[spaces[0]]
kernel = space_obj.distance_array(
x.val.get_axes_local_distribution_strategy(axes=axes))
kernel = kernel.apply_scalar_function(
x.domain[spaces[0]].codomain_smoothing_function(self.sigma))
kernel = codomain.distance_array(
distribution_strategy=axes_local_distribution_strategy)
kernel.apply_scalar_function(
codomain.get_smoothing_kernel_function(self.sigma),
inplace=True)
# local data
local_val = smooth_out.val.get_local_data(copy=False)
# now, apply the kernel to transformed_x
# this is done node-locally utilizing numpys reshaping in order to
# apply the kernel to the correct axes
local_transformed_x = transformed_x.val.get_local_data(copy=False)
local_kernel = kernel.get_local_data(copy=False)
# extract local kernel and reshape
local_kernel = kernel.get_local_data(copy=False)
new_shape = np.ones(len(local_val.shape), dtype=np.int)
for space_axis, val_axis in zip(range(len(space_obj.shape)), axes):
new_shape[val_axis] = local_kernel.shape[space_axis]
local_kernel = local_kernel.reshape(new_shape)
reshaper = [transformed_x.shape[i] if i in coaxes else 1
for i in xrange(len(transformed_x.shape))]
local_kernel = np.reshape(local_kernel, reshaper)
# multiply kernel
if inverse:
local_val /= kernel
else:
local_val *= kernel
# apply the kernel
if inverse:
local_transformed_x /= local_kernel
else:
local_transformed_x *= local_kernel
smooth_out.val.set_local_data(local_val, copy=False)
transformed_x.val.set_local_data(local_transformed_x, copy=False)
# inverse transform
smooth_out = transform.inverse_times(smooth_out, spaces=spaces[0])
result = Transformator.inverse_times(transformed_x, spaces=spaces)
return smooth_out
return result
......@@ -153,16 +153,9 @@ class GLSpace(Space):
return result_x
def _distance_array_helper(self, qr_tuple):
numerator = np.sqrt(np.sin(qr_tuple[1])**2 +
(np.sin(qr_tuple[0]) * np.cos(qr_tuple[1]))**2)
denominator = np.cos(qr_tuple[0]) * np.cos(qr_tuple[1])
return np.arctan(numerator / denominator)
def distance_array(self, distribution_strategy):
dists = arange(
start = 0, stop = self.shape[0], dtype=np.float128,
start=0, stop=self.shape[0], dtype=np.float128,
distribution_strategy=distribution_strategy
)
......@@ -172,12 +165,15 @@ class GLSpace(Space):
return dists
def _distance_array_helper(self, qr_tuple):
numerator = np.sqrt(np.sin(qr_tuple[1])**2 +
(np.sin(qr_tuple[0]) * np.cos(qr_tuple[1]))**2)
denominator = np.cos(qr_tuple[0]) * np.cos(qr_tuple[1])
def codomain_smoothing_function(self, sigma, target):
if sigma is None:
sigma = np.sqrt(2) * np.pi / (target.lmax + 1)
return np.arctan(numerator / denominator)
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma**2)
def get_smoothing_kernel_function(self, sigma):
raise NotImplementedError
# ---Added properties and methods---
......
......@@ -126,39 +126,6 @@ class HPSpace(Space):
self._nside = self._parse_nside(nside)
def distance_array(self, distribution_strategy):
"""
Calculates distance from center to all the points on the sphere
Parameters
----------
distribution_strategy: Result d2o's distribution strategy
Returns
-------
dists: distributed_data_object
"""
dists = arange(
start=0, stop = self.shape[0], dtype=np.float128,
distribution_strategy=distribution_strategy
)
# setting the center to fixed value
center_vec = (1, 0, 0)
dists = dists.apply_scalar_function(
lambda x: np.arccos(np.dot(hp.pix2vec(self.nside, int(x)),
center_vec))
)
return dists
def codomain_smoothing_function(self, sigma, target):
if sigma is None:
sigma = np.sqrt(2) * np.pi / (target.lmax + 1)
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma**2)
# ---Mandatory properties and methods---
@property
......@@ -192,6 +159,36 @@ class HPSpace(Space):
return result_x
def distance_array(self, distribution_strategy):
"""
Calculates distance from center to all the points on the sphere
Parameters
----------
distribution_strategy: Result d2o's distribution strategy
Returns
-------
dists: distributed_data_object
"""
dists = arange(
start=0, stop=self.shape[0], dtype=np.float128,
distribution_strategy=distribution_strategy
)
# setting the center to fixed value
center_vec = (1, 0, 0)
dists = dists.apply_scalar_function(
lambda x: np.arccos(np.dot(hp.pix2vec(self.nside, int(x)),
center_vec))
)
return dists
def get_smoothing_kernel_function(self, sigma, target):
raise NotImplementedError
# ---Added properties and methods---
@property
......
......@@ -110,10 +110,14 @@ class LMSpace(Space):
self._lmax = self._parse_lmax(lmax)
self._mmax = self._parse_mmax(mmax)
def compute_k_array(self, distribution_strategy):
# return a d2o with the l-value at the individual pixels
# e.g. for 0<=l<=2: [0,1,2, 1,1,2,2, 2,2]
pass
def distance_array(self, distribution_strategy):
raise NotImplementedError
def get_smoothing_kernel_function(self, sigma):
if sigma is None:
sigma = np.sqrt(2) * np.pi / (self.lmax + 1)
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma**2)
# ---Mandatory properties and methods---
......
......@@ -48,10 +48,6 @@ class PowerSpace(Space):
self._pundex = power_index['pundex']
self._k_array = power_index['k_array']
def compute_k_array(self, distribution_strategy):
raise NotImplementedError(about._errors.cstring(
"ERROR: There is no k_array implementation for PowerSpace."))
def pre_cast(self, x, axes=None):
if callable(x):
return x(self.kindex)
......@@ -109,6 +105,15 @@ class PowerSpace(Space):
return result_x
def distance_array(self, distribution_strategy):
raise NotImplementedError(about._errors.cstring(
"ERROR: There is no distance_array implementation for PowerSpace.")
)
def get_smoothing_kernel_function(self, sigma):
raise NotImplementedError(about._errors.cstring(
"ERROR: There is no smoothing function for PowerSpace."))
# ---Added properties and methods---
@property
......
......@@ -152,64 +152,6 @@ class RGSpace(Space):
self._distances = self._parse_distances(distances)
self._zerocenter = self._parse_zerocenter(zerocenter)
def distance_array(self, distribution_strategy):
"""
Calculates an n-dimensional array with its entries being the
lengths of the k-vectors from the zero point of the grid.
Parameters
----------
None : All information is taken from the parent object.
Returns
-------
nkdict : distributed_data_object
"""
shape = self.shape
# prepare the distributed_data_object
nkdict = distributed_data_object(
global_shape=shape,
dtype=np.float128,
distribution_strategy=distribution_strategy)
if distribution_strategy in DISTRIBUTION_STRATEGIES['slicing']:
# get the node's individual slice of the first dimension
slice_of_first_dimension = slice(
*nkdict.distributor.local_slice[0:2])
elif distribution_strategy in DISTRIBUTION_STRATEGIES['not']:
slice_of_first_dimension = slice(0, shape[0])
else:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported distribution strategy"))
dists = self._distance_array_helper(slice_of_first_dimension)
nkdict.set_local_data(dists)
return nkdict
def _distance_array_helper(self, slice_of_first_dimension):
dk = self.distances
shape = self.shape
inds = []
for a in shape:
inds += [slice(0, a)]
cords = np.ogrid[inds]
dists = ((np.float128(0) + cords[0] - shape[0] // 2) * dk[0])**2
# apply zerocenterQ shift
if self.zerocenter[0] == False:
dists = np.fft.fftshift(dists)
# only save the individual slice
dists = dists[slice_of_first_dimension]
for ii in range(1, len(shape)):
temp = ((cords[ii] - shape[ii] // 2) * dk[ii])**2
if self.zerocenter[ii] == False:
temp = np.fft.fftshift(temp)
dists = dists + temp
dists = np.sqrt(dists)
return dists
def hermitian_decomposition(self, x, axes=None):
# compute the hermitian part
flipped_x = self._hermitianize_inverter(x, axes=axes)
......@@ -284,6 +226,70 @@ class RGSpace(Space):
result_x = x*weight
return result_x
def distance_array(self, distribution_strategy):
"""
Calculates an n-dimensional array with its entries being the
lengths of the k-vectors from the zero point of the grid.
Parameters
----------
None : All information is taken from the parent object.
Returns
-------
nkdict : distributed_data_object
"""
shape = self.shape
# prepare the distributed_data_object
nkdict = distributed_data_object(
global_shape=shape,
dtype=np.float128,
distribution_strategy=distribution_strategy)
if distribution_strategy in DISTRIBUTION_STRATEGIES['slicing']:
# get the node's individual slice of the first dimension
slice_of_first_dimension = slice(
*nkdict.distributor.local_slice[0:2])
elif distribution_strategy in DISTRIBUTION_STRATEGIES['not']:
slice_of_first_dimension = slice(0, shape[0])
else:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported distribution strategy"))
dists = self._distance_array_helper(slice_of_first_dimension)
nkdict.set_local_data(dists)
return nkdict
def _distance_array_helper(self, slice_of_first_dimension):
dk = self.distances
shape = self.shape
inds = []
for a in shape:
inds += [slice(0, a)]
cords = np.ogrid[inds]
dists = ((np.float128(0) + cords[0] - shape[0] // 2) * dk[0])**2
# apply zerocenterQ shift
if self.zerocenter[0] == False:
dists = np.fft.fftshift(dists)
# only save the individual slice
dists = dists[slice_of_first_dimension]
for ii in range(1, len(shape)):
temp = ((cords[ii] - shape[ii] // 2) * dk[ii])**2
if self.zerocenter[ii] == False:
temp = np.fft.fftshift(temp)
dists = dists + temp
dists = np.sqrt(dists)
return dists
def get_smoothing_kernel_function(self, sigma):
if sigma is None:
sigma = np.sqrt(2) * np.max(self.distances)
return lambda x: np.exp(-2. * np.pi**2 * x**2 * sigma**2)
# ---Added properties and methods---
@property
......@@ -317,8 +323,3 @@ class RGSpace(Space):
temp[:] = zerocenter
return tuple(temp)
def codomain_smoothing_function(self, sigma):
if sigma is None:
sigma = np.sqrt(2) * np.max(self.distances)
return lambda x: np.exp(-2. * np.pi**2 * x**2 * sigma**2)
......@@ -223,6 +223,12 @@ class Space(object):
else:
return False
def pre_cast(self, x, axes=None):
return x
def post_cast(self, x, axes=None):
return x
@abc.abstractproperty
def harmonic(self):
raise NotImplementedError
......@@ -266,20 +272,13 @@ class Space(object):
"""
raise NotImplementedError
def pre_cast(self, x, axes=None):
return x
def post_cast(self, x, axes=None):
return x
def distance_array(self, distribution_strategy):
raise NotImplementedError(about._errors.cstring(
"ERROR: There is no generic distance_array for Space base class."))
def codomain_smoothing_function(self, sigma, target):
def get_smoothing_kernel_function(self, sigma):
raise NotImplementedError(about._errors.cstring(
"ERROR: There is no generic smoothing fuction for Space base class."
))
"ERROR: There is no generic smoothing function for Space class."))
def hermitian_decomposition(self, x, axes=None):
raise NotImplementedError
......
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