Commit 59edf42b authored by theos's avatar theos

Fixed a bug in LinearOperator in _check_input_compatibility.

Modified power_analyze in Field for complex power spectra.
parent 93d877c3
......@@ -177,8 +177,19 @@ class Field(object):
return random_arguments
def power_analyze(self, spaces=None, log=False, nbin=None, binbounds=None):
# check if the spaces input is valid
# ---Powerspectral methods---
def power_analyze(self, spaces=None, log=False, nbin=None, binbounds=None,
real_signal=True):
# assert that all spaces in `self.domain` are either harmonic or
# power_space instances
for sp in self.domain:
if not sp.harmonic and not isinstance(sp, PowerSpace):
raise AttributeError(
"ERROR: Field has a space in `domain` which is neither "
"harmonic nor a PowerSpace.")
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
if len(self.domain) == 1:
......@@ -201,20 +212,48 @@ class Field(object):
raise ValueError(about._errors.cstring(
"ERROR: Conversion of only one space at a time is allowed."))
# create the target PowerSpace instance
# Create the target PowerSpace instance:
# If the associated signal-space field was real, we extract the
# hermitian and anti-hermitian parts of `self` and put them
# into the real and imaginary parts of the power spectrum.
# If it was complex, all the power is put into a real power spectrum.
distribution_strategy = \
self.val.get_axes_local_distribution_strategy(
self.domain_axes[space_index])
if real_signal:
power_dtype = np.dtype('complex')
else:
power_dtype = np.dtype('float')
harmonic_domain = self.domain[space_index]
power_domain = PowerSpace(harmonic_domain=harmonic_domain,
datamodel=distribution_strategy,
log=log, nbin=nbin, binbounds=binbounds)
log=log, nbin=nbin, binbounds=binbounds,
dtype=power_dtype)
# extract pindex and rho from power_domain and calculate the spectrum
# extract pindex and rho from power_domain
pindex = power_domain.pindex
rho = power_domain.rho
power_spectrum = self._calculate_power_spectrum(
if real_signal:
hermitian_part, anti_hermitian_part = \
harmonic_domain.hermitian_decomposition(
self.val,
axes=self.domain_axes[space_index])
[hermitian_power, anti_hermitian_power] = \
[self._calculate_power_spectrum(
x=part,
pindex=pindex,
rho=rho,
axes=self.domain_axes[space_index])
for part in [hermitian_part, anti_hermitian_part]]
power_spectrum = hermitian_power + 1j * anti_hermitian_power
else:
power_spectrum = self._calculate_power_spectrum(
x=self.val,
pindex=pindex,
rho=rho,
......@@ -274,9 +313,17 @@ class Field(object):
return result_obj
def power_synthesize(self):
# check that all spaces in self.domain are real or instances of power_space
# check if field is real- or complex-valued
def power_synthesize(self, spaces=None, real_signal=True):
# assert that all spaces in `self.domain` are eiher of signal-type or
# power_space instances
for sp in self.domain:
if sp.harmonic and not isinstance(sp, PowerSpace):
raise AttributeError(
"ERROR: Field has a space in `domain` which is neither "
"harmonic nor a PowerSpace.")
# synthesize random fields in harmonic domain using
# np.random.multivariate_normal(mean=[0,0], cov=[[0.5,0],[0,0.5]], size=shape)
# ---Properties---
......@@ -365,12 +412,11 @@ class Field(object):
if dtype is None:
dtype = self.dtype
x = distributed_data_object(x,
global_shape=self.shape,
dtype=dtype,
distribution_strategy=self.datamodel)
return x
return_x = distributed_data_object(global_shape=self.shape,
dtype=dtype,
distribution_strategy=self.datamodel)
return_x.set_full_data(x, copy=False)
return return_x
def copy(self, domain=None, dtype=None, field_type=None,
datamodel=None):
......
......@@ -92,17 +92,18 @@ def hermitianize(x, axes=None):
# make the point inversions
flipped_x = _hermitianize_inverter(x, axes=axes)
flipped_x = flipped_x.conjugate()
# check if x was already hermitian
if (x == flipped_x).all():
return x
# average x and flipped_x.
x = (x + flipped_x) / 2.
try:
x.hermitian = True
except(AttributeError):
pass
# x = (x + flipped_x) / 2.
result_x = x + flipped_x
result_x /= 2.
return x
# try:
# x.hermitian = True
# except(AttributeError):
# pass
return result_x
def _hermitianize_inverter(x, axes):
......
......@@ -82,7 +82,8 @@ class LinearOperator(object):
return y
def inverse_times(self, x, spaces=None, types=None):
spaces, types = self._check_input_compatibility(x, spaces, types)
spaces, types = self._check_input_compatibility(x, spaces, types,
inverse=True)
y = self._inverse_times(x, spaces, types)
if not self.implemented:
......@@ -93,7 +94,8 @@ class LinearOperator(object):
if self.unitary:
return self.inverse_times(x, spaces, types)
spaces, types = self._check_input_compatibility(x, spaces, types)
spaces, types = self._check_input_compatibility(x, spaces, types,
inverse=True)
if not self.implemented:
x = x.weight(spaces=spaces)
......@@ -142,7 +144,7 @@ class LinearOperator(object):
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'inverse_adjoint_times'."))
def _check_input_compatibility(self, x, spaces, types):
def _check_input_compatibility(self, x, spaces, types, inverse=False):
if not isinstance(x, Field):
raise ValueError(about._errors.cstring(
"ERROR: supplied object is not a `nifty.Field`."))
......@@ -160,42 +162,49 @@ class LinearOperator(object):
# 2. Case:
# The domains of self and x match completely.
if not inverse:
self_domain = self.domain
self_field_type = self.field_type
else:
self_domain = self.target
self_field_type = self.field_type_target
if spaces is None:
if self.domain != () and self.domain != x.domain:
if self_domain != () and self_domain != x.domain:
raise ValueError(about._errors.cstring(
"ERROR: The operator's and and field's domains don't "
"match."))
else:
if len(self.domain) > 1:
if len(self_domain) > 1:
raise ValueError(about._errors.cstring(
"ERROR: Specifying `spaces` for operators with multiple "
"domain spaces is not valid."))
elif len(spaces) != len(self.domain):
elif len(spaces) != len(self_domain):
raise ValueError(about._errors.cstring(
"ERROR: Length of `spaces` does not match the number of "
"spaces in the operator's domain."))
elif len(spaces) == 1:
if x.domain[spaces[0]] != self.domain[0]:
if x.domain[spaces[0]] != self_domain[0]:
raise ValueError(about._errors.cstring(
"ERROR: The operator's and and field's domains don't "
"match."))
if types is None:
if self.field_type != () and self.field_type != x.field_type:
if self_field_type != () and self_field_type != x.field_type:
raise ValueError(about._errors.cstring(
"ERROR: The operator's and and field's field_types don't "
"match."))
else:
if len(self.field_type) > 1:
if len(self_field_type) > 1:
raise ValueError(about._errors.cstring(
"ERROR: Specifying `types` for operators with multiple "
"field-types is not valid."))
elif len(types) != len(self.field_type):
elif len(types) != len(self_field_type):
raise ValueError(about._errors.cstring(
"ERROR: Length of `types` does not match the number of "
"the operator's field-types."))
elif len(types) == 1:
if x.field_type[types[0]] != self.field_type[0]:
if x.field_type[types[0]] != self_field_type[0]:
raise ValueError(about._errors.cstring(
"ERROR: The operator's and and field's field_type "
"don't match."))
......
......@@ -44,6 +44,8 @@ class PowerSpace(Space):
self._pindex = power_index['pindex']
self._kindex = power_index['kindex']
self._rho = power_index['rho']
self._pundex = power_index['pundex']
self._k_array = power_index['k_array']
def compute_k_array(self, distribution_strategy):
raise NotImplementedError(about._errors.cstring(
......@@ -130,3 +132,10 @@ class PowerSpace(Space):
def rho(self):
return self._rho
@property
def pundex(self):
return self._pundex
@property
def k_array(self):
return self._k_array
......@@ -203,6 +203,46 @@ class RGSpace(Space):
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)
flipped_x = flipped_x.conjugate()
# average x and flipped_x.
hermitian_part = x + flipped_x
hermitian_part /= 2.
# use subtraction since it is faster than flipping another time
anti_hermitian_part = (x-hermitian_part)/1j
return (hermitian_part, anti_hermitian_part)
def _hermitianize_inverter(self, x, axes):
# calculate the number of dimensions the input array has
dimensions = len(x.shape)
# prepare the slicing object which will be used for mirroring
slice_primitive = [slice(None), ] * dimensions
# copy the input data
y = x.copy()
if axes is None:
axes = xrange(dimensions)
# flip in the desired directions
for i in axes:
slice_picker = slice_primitive[:]
slice_picker[i] = slice(1, None, None)
slice_picker = tuple(slice_picker)
slice_inverter = slice_primitive[:]
slice_inverter[i] = slice(None, 0, -1)
slice_inverter = tuple(slice_inverter)
try:
y.set_data(to_key=slice_picker, data=y,
from_key=slice_inverter)
except(AttributeError):
y[slice_picker] = y[slice_inverter]
return y
# ---Mandatory properties and methods---
@property
......
......@@ -273,6 +273,9 @@ class Space(object):
raise NotImplementedError(about._errors.cstring(
"ERROR: There is no generic k_array for Space base class."))
def hermitian_decomposition(self, x, axes=None):
raise NotImplementedError
def __repr__(self):
string = ""
string += str(type(self)) + "\n"
......
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