Commit 30a97cc9 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

make power_synthesize easier to understand

parent fbd91b74
Pipeline #18295 failed with stage
in 5 minutes and 35 seconds
......@@ -2,7 +2,6 @@
import numpy as np
import nifty2go as ift
from nifty2go import plotting
if __name__ == "__main__":
signal_to_noise = 1.5 # The signal to noise ratio
......@@ -56,8 +55,7 @@ if __name__ == "__main__":
mock_power = ift.Field(domain=(power_space_1, power_space_2),
val=np.outer(mock_power_1.val, mock_power_2.val))
diagonal = mock_power.power_synthesize(spaces=(0, 1), mean=1, std=0,
real_signal=False)**2
diagonal = mock_power.power_synthesize_special(spaces=(0, 1))**2
diagonal = diagonal.real
S = ift.DiagonalOperator(domain=(harmonic_space_1, harmonic_space_2),
......@@ -109,7 +107,7 @@ if __name__ == "__main__":
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
ift.plotting.plot(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), name='uncertainty.pdf',zmin=0.,zmax=3.,title="Uncertainty map",xlabel="x")
ift.plotting.plot(ift.Field(plot_space, val=mock_signal.val.real), name='mock_signal.pdf')
ift.plotting.plot(ift.Field(plot_space, val=data.val.real), name='data.pdf')
ift.plotting.plot(ift.Field(plot_space, val=m.val.real), name='map.pdf')
ift.plotting.plot(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), name='uncertainty.pdf',zmin=0.,zmax=3.,title="Uncertainty map",colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=mock_signal.val.real), name='mock_signal.pdf',colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=data.val.real), name='data.pdf',colormap="Planck-like")
ift.plotting.plot(ift.Field(plot_space, val=m.val.real), name='map.pdf',colormap="Planck-like")
......@@ -57,7 +57,7 @@ if __name__ == "__main__":
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) #|\label{code:wf_variance_weighting}|
exit()
# Plotting #|\label{code:wf_plotting}|
ift.plotting.plot(variance,name="uncertainty.pdf",xlabel='Pixel index', ylabel='Pixel index')
ift.plotting.plot(mock_signal,name="mock_signal.pdf",xlabel='Pixel index', ylabel='Pixel index')
......
......@@ -18,13 +18,8 @@
from __future__ import division
from builtins import range
import numpy as np
from .domain_object import DomainObject
from .spaces.power_space import PowerSpace
from . import nifty_utilities as utilities
from .random import Random
from functools import reduce
......@@ -161,8 +156,7 @@ class Field(object):
"""
generator_function = getattr(Random, random_type)
return Field(domain=domain,
val=generator_function(dtype=dtype,
return Field(domain=domain, val=generator_function(dtype=dtype,
shape=utilities.domains2shape(domain), **kwargs))
# ---Powerspectral methods---
......@@ -232,23 +226,17 @@ class Field(object):
raise ValueError("No space for analysis specified.")
if keep_phase_information:
parts = [self.real.copy(), self.imag.copy()]
parts = [self.real*self.real, self.imag*self.imag]
else:
parts = [self]
parts = [abs(part)**2 for part in parts]
parts = [self.real*self.real + self.imag*self.imag]
for space_index in spaces:
parts = [self._single_power_analyze(
work_field=part,
space_index=space_index,
binbounds=binbounds)
parts = [self._single_power_analyze(work_field=part,
space_index=space_index,
binbounds=binbounds)
for part in parts]
if keep_phase_information:
return parts[0] + 1j*parts[1]
else:
return parts[0]
return parts[0] + 1j*parts[1] if keep_phase_information else parts[0]
@staticmethod
def _single_power_analyze(work_field, space_index, binbounds):
......@@ -278,7 +266,6 @@ class Field(object):
@staticmethod
def _calculate_power_spectrum(field_val, pdomain, axes=None):
pindex = pdomain.pindex
if axes is not None:
pindex = Field._shape_up_pindex(pindex, field_val.shape, axes)
......@@ -287,7 +274,7 @@ class Field(object):
axis=axes)
rho = pdomain.rho
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
new_rho_shape = [1] * len(power_spectrum.shape)
new_rho_shape[axes[0]] = len(rho)
rho = rho.reshape(new_rho_shape)
power_spectrum /= rho
......@@ -303,8 +290,26 @@ class Field(object):
result_obj[()] = pindex.reshape(semiscaled_local_shape)
return result_obj
def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
mean=0., std=1.):
def _prep_powersynth(self, spaces):
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
spaces = range(len(self.domain))
# create the result domain
result_domain = list(self.domain)
for i in spaces:
if not isinstance(self.domain[i], PowerSpace):
raise ValueError("A PowerSpace is needed for field "
"synthetization.")
result_domain[i] = self.domain[i].harmonic_partner
spec = np.sqrt(self.val)
for i in spaces:
spec = self._spec_to_rescaler(spec, i)
return (result_domain, spec)
def power_synthesize(self, spaces=None, real_power=True, real_signal=True):
""" Yields a sampled field with `self`**2 as its power spectrum.
This method draws a Gaussian random field in the harmonic partner
......@@ -322,14 +327,6 @@ class Field(object):
real_signal : boolean *optional*
True will result in a purely real signal-space field
(default : True).
mean : float *optional*
The mean of the Gaussian noise field which is used for the Field
synthetization (default : None).
if mean==None : mean will be set to 0
std : float *optional*
The standard deviation of the Gaussian noise field which is used
for the Field synthetization (default : None).
if std==None : std will be set to 1
Returns
-------
......@@ -353,51 +350,36 @@ class Field(object):
"""
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
spaces = range(len(self.domain))
for i in spaces:
if not isinstance(self.domain[i], PowerSpace):
raise ValueError("A PowerSpace is needed for field "
"synthetization.")
# create the result domain
result_domain = list(self.domain)
for i in spaces:
result_domain[i] = self.domain[i].harmonic_partner
result_domain, spec = self._prep_powersynth(spaces)
# create random samples: one or two, depending on whether the
# power spectrum is real or complex
result_list = [self.__class__.from_random(
'normal',
mean=mean,
std=std,
domain=result_domain,
dtype=np.complex)
for x in range(1 if real_power else 2)]
# from now on extract the values from the random fields for further
# processing without killing the fields.
# if the signal-space field should be real, hermitianize the field
# components
spec = np.sqrt(self.val)
for power_space_index in spaces:
spec = self._spec_to_rescaler(spec, power_space_index)
result = [self.from_random('normal', mean=0., std=1.,
domain=result_domain,
dtype=np.float if real_signal
else np.complex)
for x in range(1 if real_power else 2)]
# MR: dummy call - will be removed soon
if real_signal:
self.from_random('normal', mean=0., std=1.,
domain=result_domain, dtype=np.float)
# apply the rescaler to the random fields
result_list[0] *= spec.real
result[0] *= spec.real
if not real_power:
result_list[1] *= spec.imag
result[1] *= spec.imag
if real_signal:
result_list = [i.real*np.sqrt(2.) for i in result_list]
return result[0] if real_power else result[0] + 1j*result[1]
if real_power:
return result_list[0]
else:
return result_list[0] + 1j*result_list[1]
def power_synthesize_special(self, spaces=None):
result_domain, spec = self._prep_powersynth(spaces)
# MR: dummy call - will be removed soon
self.from_random('normal', mean=0., std=1.,
domain=result_domain, dtype=np.complex)
return spec.real
def _spec_to_rescaler(self, spec, power_space_index):
power_space = self.domain[power_space_index]
......@@ -535,7 +517,8 @@ class Field(object):
fct *= wgt
else:
new_shape = np.ones(len(self.shape), dtype=np.int)
new_shape[self.domain_axes[ind][0]:self.domain_axes[ind][-1]+1]=wgt.shape
new_shape[self.domain_axes[ind][0]:
self.domain_axes[ind][-1]+1] = wgt.shape
wgt = wgt.reshape(new_shape)
new_field *= wgt**power
fct = fct**power
......
......@@ -104,7 +104,7 @@ class DiagonalOperator(EndomorphicOperator):
def _adjoint_inverse_times(self, x, spaces):
return self._times_helper(x, spaces,
operation=lambda z:
z.conjugate().__rtruediv__)
z.conjugate().__rtruediv__)
def diagonal(self, bare=False, copy=True):
""" Returns the diagonal of the Operator.
......@@ -219,4 +219,4 @@ class DiagonalOperator(EndomorphicOperator):
reshaped_local_diagonal = np.reshape(self._diagonal.val, reshaper)
# here the actual multiplication takes place
return Field(x.domain,val=operation(reshaped_local_diagonal)(x.val))
return Field(x.domain, val=operation(reshaped_local_diagonal)(x.val))
......@@ -40,8 +40,10 @@ class Random(object):
else:
x = np.random.normal(loc=0., scale=1., size=shape)
x = x.astype(dtype, copy=False)
x *= std
x += mean
if std != 1.:
x *= std
if mean != 0.:
x += mean
return x
@staticmethod
......
......@@ -58,10 +58,9 @@ def create_power_operator(domain, power_spectrum, dtype=None):
raise TypeError("power_spectrum must be callable")
power_domain = PowerSpace(domain)
fp = Field(power_domain,
val=power_spectrum(power_domain.kindex), dtype=dtype)
# MR FIXME: why generate a non-random random field?
f = fp.power_synthesize(mean=1, std=0, real_signal=False)
fp = Field(power_domain, val=power_spectrum(power_domain.kindex),
dtype=dtype)
f = fp.power_synthesize_special()
if not issubclass(fp.dtype.type, np.complexfloating):
f = f.real
......
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