Commit 7ae27a71 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into 'mpitests'

Master

See merge request !171
parents e7d60d48 e1e13bd9
Pipeline #15052 passed with stage
in 9 minutes and 15 seconds
...@@ -97,10 +97,8 @@ if __name__ == "__main__": ...@@ -97,10 +97,8 @@ if __name__ == "__main__":
# The information source # The information source
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"], realized_power = log(sh.power_analyze(binbounds=p_space.binbounds))
nbin=p_space.config["nbin"])) data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"]))
d_data = d.val.get_full_data().real d_data = d.val.get_full_data().real
if rank == 0: if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html') pl.plot([go.Heatmap(z=d_data)], filename='data.html')
......
import numpy as np
from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
SmoothingOperator, DiagonalOperator, create_power_operator
from nifty.library import WienerFilterCurvature
#import plotly.offline as pl
#import plotly.graph_objs as go
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
if __name__ == "__main__":
distribution_strategy = 'fftw'
# Setting up physical constants
# total length of Interval or Volume the field lives on, e.g. in meters
L = 2.
# typical distance over which the field is correlated (in same unit as L)
correlation_length = 0.1
# variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
field_variance = 2.
# smoothing length of response (in same unit as L)
response_sigma = 0.1
# defining resolution (pixels per dimension)
N_pixels = 512
# Setting up derived constants
k_0 = 1./correlation_length
# note that field_variance**2 = a*k_0/4. for this analytic form of power
# spectrum
a = field_variance**2/k_0*4.
pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
pixel_length = L/N_pixels
# Setting up the geometry
s_space = RGSpace([N_pixels, N_pixels], distances=pixel_length)
fft = FFTOperator(s_space, domain_dtype=np.float, target_dtype=np.complex)
h_space = fft.target[0]
inverse_fft = FFTOperator(h_space, target=s_space,
domain_dtype=np.complex, target_dtype=np.float)
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
sp = Field(p_space, val=pow_spec,
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
ss = fft.inverse_times(sh)
R = SmoothingOperator(s_space, sigma=response_sigma)
R_harmonic = ComposedOperator([inverse_fft, R], default_spaces=[0, 0])
signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
n = Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
mean=0)
d = R(ss) + n
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(d))
wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic)
m = wiener_curvature.inverse_times(j)
m_s = inverse_fft(m)
import numpy as np
from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
DiagonalOperator, ResponseOperator, plotting,\
create_power_operator
from nifty.library import WienerFilterCurvature
if __name__ == "__main__":
distribution_strategy = 'not'
# Setting up variable parameters
# Typical distance over which the field is correlated
correlation_length = 0.01
# Variance of field in position space sqrt(<|s_x|^2>)
field_variance = 2.
# smoothing length of response (in same unit as L)
response_sigma = 0.1
# The signal to noise ratio
signal_to_noise = 0.7
# note that field_variance**2 = a*k_0/4. for this analytic form of power
# spectrum
def power_spectrum(k):
a = 4 * correlation_length * field_variance**2
return a / (1 + k * correlation_length) ** 4
# Setting up the geometry
# Total side-length of the domain
L = 2.
# Grid resolution (pixels per axis)
N_pixels = 512
signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
fft = FFTOperator(harmonic_space, target=signal_space,
domain_dtype=np.complex, target_dtype=np.float)
power_space = PowerSpace(harmonic_space,
distribution_strategy=distribution_strategy)
# Creating the mock data
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum,
distribution_strategy=distribution_strategy)
mock_power = Field(power_space, val=power_spectrum,
distribution_strategy=distribution_strategy)
np.random.seed(43)
mock_harmonic = mock_power.power_synthesize(real_signal=True)
mock_signal = fft(mock_harmonic)
R = ResponseOperator(signal_space, sigma=(response_sigma,))
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
N = DiagonalOperator(data_domain,
diagonal=mock_signal.var()/signal_to_noise,
bare=True)
noise = Field.from_random(domain=data_domain,
random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0)
data = R(mock_signal) + noise
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic)
m = wiener_curvature.inverse_times(j)
m_s = fft(m)
plotter = plotting.RG2DPlotter()
plotter.title = 'mock_signal.html';
plotter(mock_signal)
plotter.title = 'data.html'
plotter(Field(signal_space,
val=data.val.get_full_data().reshape(signal_space.shape)))
plotter.title = 'map.html'; plotter(m_s)
\ No newline at end of file
...@@ -18,6 +18,7 @@ ...@@ -18,6 +18,7 @@
from __future__ import division from __future__ import division
import ast
import itertools import itertools
import numpy as np import numpy as np
...@@ -270,7 +271,7 @@ class Field(Loggable, Versionable, object): ...@@ -270,7 +271,7 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods--- # ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=False, nbin=None, def power_analyze(self, spaces=None, logarithmic=None, nbin=None,
binbounds=None, keep_phase_information=False): binbounds=None, keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`. """ Computes the square root power spectrum for a subspace of `self`.
...@@ -287,14 +288,15 @@ class Field(Loggable, Versionable, object): ...@@ -287,14 +288,15 @@ class Field(Loggable, Versionable, object):
(default : None). (default : None).
logarithmic : boolean *optional* logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning. True if the output PowerSpace should use logarithmic binning.
{default : False} {default : None}
nbin : int *optional* nbin : int *optional*
The number of bins the resulting PowerSpace shall have The number of bins the resulting PowerSpace shall have
(default : None). (default : None).
if nbin==None : maximum number of bins is used if nbin==None : maximum number of bins is used
binbounds : array-like *optional* binbounds : array-like *optional*
Inner bounds of the bins (default : None). Inner bounds of the bins (default : None).
if binbounds==None : bins are inferred. Overwrites nbins and log Overrides nbin and logarithmic.
if binbounds==None : bins are inferred.
keep_phase_information : boolean, *optional* keep_phase_information : boolean, *optional*
If False, return a real-valued result containing the power spectrum If False, return a real-valued result containing the power spectrum
of the input Field. of the input Field.
...@@ -397,14 +399,9 @@ class Field(Loggable, Versionable, object): ...@@ -397,14 +399,9 @@ class Field(Loggable, Versionable, object):
logarithmic=logarithmic, nbin=nbin, logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds) binbounds=binbounds)
# extract pindex and rho from power_domain
pindex = power_domain.pindex
rho = power_domain.rho
power_spectrum = cls._calculate_power_spectrum( power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val, field_val=work_field.val,
pindex=pindex, pdomain=power_domain,
rho=rho,
axes=work_field.domain_axes[space_index]) axes=work_field.domain_axes[space_index])
# create the result field and put power_spectrum into it # create the result field and put power_spectrum into it
...@@ -421,8 +418,11 @@ class Field(Loggable, Versionable, object): ...@@ -421,8 +418,11 @@ class Field(Loggable, Versionable, object):
return result_field return result_field
@classmethod @classmethod
def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None): def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
pindex = pdomain.pindex
# MR FIXME: how about iterating over slices, instead of replicating
# pindex? Would save memory and probably isn't slower.
if axes is not None: if axes is not None:
pindex = cls._shape_up_pindex( pindex = cls._shape_up_pindex(
pindex=pindex, pindex=pindex,
...@@ -431,6 +431,7 @@ class Field(Loggable, Versionable, object): ...@@ -431,6 +431,7 @@ class Field(Loggable, Versionable, object):
axes=axes) axes=axes)
power_spectrum = pindex.bincount(weights=field_val, power_spectrum = pindex.bincount(weights=field_val,
axis=axes) axis=axes)
rho = pdomain.rho
if axes is not None: 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) new_rho_shape[axes[0]] = len(rho)
...@@ -755,7 +756,7 @@ class Field(Loggable, Versionable, object): ...@@ -755,7 +756,7 @@ class Field(Loggable, Versionable, object):
Returns Returns
------- -------
out : tuple out : tuple
The output object. The tuple contains the dimansions of the spaces The output object. The tuple contains the dimensions of the spaces
in domain. in domain.
See Also See Also
...@@ -1519,7 +1520,8 @@ class Field(Loggable, Versionable, object): ...@@ -1519,7 +1520,8 @@ class Field(Loggable, Versionable, object):
temp_domain.append(repository.get('s_' + str(i), hdf5_group)) temp_domain.append(repository.get('s_' + str(i), hdf5_group))
new_field.domain = tuple(temp_domain) new_field.domain = tuple(temp_domain)
exec('new_field.domain_axes = ' + hdf5_group.attrs['domain_axes']) new_field.domain_axes = ast.literal_eval(
hdf5_group.attrs['domain_axes'])
try: try:
new_field._val = repository.get('val', hdf5_group) new_field._val = repository.get('val', hdf5_group)
......
...@@ -103,14 +103,12 @@ class CriticalPowerEnergy(Energy): ...@@ -103,14 +103,12 @@ class CriticalPowerEnergy(Energy):
posterior_sample = generate_posterior_sample( posterior_sample = generate_posterior_sample(
self.m, self.D) self.m, self.D)
projected_sample = posterior_sample.power_analyze( projected_sample = posterior_sample.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"], binbounds=self.position.domain[0].binbounds)
nbin=self.position.domain[0].config["nbin"])
w += (projected_sample) * self.rho w += (projected_sample) * self.rho
w /= float(self.samples) w /= float(self.samples)
else: else:
w = self.m.power_analyze( w = self.m.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"], binbounds=self.position.domain[0].binbounds)
nbin=self.position.domain[0].config["nbin"])
w *= self.rho w *= self.rho
self._w = w self._w = w
return self._w return self._w
......
...@@ -44,9 +44,6 @@ class LineSearch(Loggable, object): ...@@ -44,9 +44,6 @@ class LineSearch(Loggable, object):
__metaclass__ = abc.ABCMeta __metaclass__ = abc.ABCMeta
def __init__(self): def __init__(self):
self.line_energy = None self.line_energy = None
self.f_k_minus_1 = None self.f_k_minus_1 = None
self.preferred_initial_step_size = None self.preferred_initial_step_size = None
......
...@@ -373,6 +373,10 @@ class MPIFFT(Transform): ...@@ -373,6 +373,10 @@ class MPIFFT(Transform):
original_shape = inp.shape original_shape = inp.shape
inp = inp.reshape(inp.shape[0], 1) inp = inp.reshape(inp.shape[0], 1)
axes = (0, ) axes = (0, )
if original_shape[0]%2!=0:
raise AttributeError("MPI-FFTs of onedimensional arrays "
"with odd length are currently not supported due to a "
"bug in FFTW. Please use a grid with even length.")
if current_info is None: if current_info is None:
transform_shape = list(inp.shape) transform_shape = list(inp.shape)
......
...@@ -38,7 +38,7 @@ class InvertibleOperatorMixin(object): ...@@ -38,7 +38,7 @@ class InvertibleOperatorMixin(object):
(default: ConjugateGradient) (default: ConjugateGradient)
preconditioner : LinearOperator preconditioner : LinearOperator
Preconditioner that is used by ConjugateGraduent if no minimizer was Preconditioner that is used by ConjugateGradient if no minimizer was
given. given.
Attributes Attributes
......
...@@ -33,8 +33,7 @@ class ResponseOperator(LinearOperator): ...@@ -33,8 +33,7 @@ class ResponseOperator(LinearOperator):
domain : tuple of DomainObjects, i.e. Spaces and FieldTypes domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
The domain on which the Operator's input Field lives. The domain on which the Operator's input Field lives.
target : tuple of DomainObjects, i.e. Spaces and FieldTypes target : tuple of DomainObjects, i.e. Spaces and FieldTypes
The domain in which the outcome of the operator lives. As the Operator The domain in which the outcome of the operator lives.
is endomorphic this is the same as its domain.
unitary : boolean unitary : boolean
Indicates whether the Operator is unitary or not. Indicates whether the Operator is unitary or not.
......
...@@ -21,7 +21,6 @@ class FFTSmoothingOperator(SmoothingOperator): ...@@ -21,7 +21,6 @@ class FFTSmoothingOperator(SmoothingOperator):
# transform to the (global-)default codomain and perform all remaining # transform to the (global-)default codomain and perform all remaining
# steps therein # steps therein
transformator = self._get_transformator(x.dtype) transformator = self._get_transformator(x.dtype)
transformed_x = transformator(x, spaces=spaces) transformed_x = transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]] codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]] coaxes = transformed_x.domain_axes[spaces[0]]
......
...@@ -21,3 +21,6 @@ class HealpixPlotter(PlotterBase): ...@@ -21,3 +21,6 @@ class HealpixPlotter(PlotterBase):
def _initialize_figure(self): def _initialize_figure(self):
return Figure2D(plots=None) return Figure2D(plots=None)
def _parse_data(self, data, field, spaces):
return data
...@@ -37,7 +37,8 @@ class Prober(object): ...@@ -37,7 +37,8 @@ class Prober(object):
""" """
def __init__(self, domain=None, distribution_strategy=None, probe_count=8, def __init__(self, domain=None, distribution_strategy=None, probe_count=8,
random_type='pm1', compute_variance=False): random_type='pm1', probe_dtype=np.float,
compute_variance=False):
self._domain = utilities.parse_domain(domain) self._domain = utilities.parse_domain(domain)
self._distribution_strategy = \ self._distribution_strategy = \
...@@ -45,6 +46,7 @@ class Prober(object): ...@@ -45,6 +46,7 @@ class Prober(object):
self._probe_count = self._parse_probe_count(probe_count) self._probe_count = self._parse_probe_count(probe_count)
self._random_type = self._parse_random_type(random_type) self._random_type = self._parse_random_type(random_type)
self.compute_variance = bool(compute_variance) self.compute_variance = bool(compute_variance)
self.probe_dtype = np.dtype(probe_dtype)
# ---Properties--- # ---Properties---
...@@ -104,6 +106,7 @@ class Prober(object): ...@@ -104,6 +106,7 @@ class Prober(object):
""" a random-probe generator """ """ a random-probe generator """
f = Field.from_random(random_type=self.random_type, f = Field.from_random(random_type=self.random_type,
domain=self.domain, domain=self.domain,
dtype=self.probe_dtype,
distribution_strategy=self.distribution_strategy) distribution_strategy=self.distribution_strategy)
uid = np.random.randint(1e18) uid = np.random.randint(1e18)
return (uid, f) return (uid, f)
......
...@@ -22,7 +22,7 @@ import numpy as np ...@@ -22,7 +22,7 @@ import numpy as np
from nifty.spaces.space import Space from nifty.spaces.space import Space
from d2o import arange from d2o import arange, distributed_data_object
class LMSpace(Space): class LMSpace(Space):
...@@ -127,6 +127,24 @@ class LMSpace(Space): ...@@ -127,6 +127,24 @@ class LMSpace(Space):
return x.copy() return x.copy()
def get_distance_array(self, distribution_strategy): def get_distance_array(self, distribution_strategy):
if distribution_strategy == 'not': # short cut
lmax = self.lmax
ldist = np.empty((self.dim,), dtype=np.float64)
ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
tmp = np.empty((2*lmax+2), dtype=np.float64)
tmp[0::2] = np.arange(lmax+1)
tmp[1::2] = np.arange(lmax+1)
idx = lmax+1
for l in range(1, lmax+1):
ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
idx += 2*(lmax+1-l)
dists = distributed_data_object(
global_shape=self.shape,
dtype=np.float,
distribution_strategy=distribution_strategy)
dists.set_local_data(ldist)
return dists
dists = arange(start=0, stop=self.shape[0], dists = arange(start=0, stop=self.shape[0],
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
...@@ -136,6 +154,12 @@ class LMSpace(Space): ...@@ -136,6 +154,12 @@ class LMSpace(Space):
return dists return dists
def get_unique_distances(self):
return np.arange(self.lmax+1, dtype=np.float64)
def get_natural_binbounds(self):
return np.arange(self.lmax, dtype=np.float64) + 0.5
@staticmethod @staticmethod
def _distance_array_helper(index_array, lmax): def _distance_array_helper(index_array, lmax):
u = 2*lmax + 1 u = 2*lmax + 1
......
...@@ -17,4 +17,3 @@ ...@@ -17,4 +17,3 @@
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
from power_space import PowerSpace from power_space import PowerSpace
from power_index_factory import PowerIndexFactory
\ No newline at end of file
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from power_indices import PowerIndices
class _PowerIndexFactory(object):
def __init__(self):
self.power_indices_storage = {}
def get_power_index(self, domain, distribution_strategy,
logarithmic=False, nbin=None, binbounds=None):
key = (domain, distribution_strategy)
if key not in self.power_indices_storage:
self.power_indices_storage[key] = \
PowerIndices(domain, distribution_strategy,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
power_indices = self.power_indices_storage[key]
power_index = power_indices.get_index_dict(logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
return power_index
PowerIndexFactory = _PowerIndexFactory()
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from d2o import distributed_data_object,\
STRATEGIES as DISTRIBUTION_STRATEGIES