Commit f23c4d72 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'master' into byebye_fixed_point_voodoo

parents 1ce8d473 1d10be46
Pipeline #15141 passed with stage
in 8 minutes and 7 seconds
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 @@
from __future__ import division
import ast
import itertools
import numpy as np
......@@ -465,7 +466,7 @@ class Field(Loggable, Versionable, object):
return result_obj
def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
mean=None, std=None):
mean=None, std=None, distribution_strategy=None):
""" Yields a sampled field with `self`**2 as its power spectrum.
This method draws a Gaussian random field in the harmonic partner
......@@ -540,13 +541,16 @@ class Field(Loggable, Versionable, object):
else:
result_list = [None, None]
if distribution_strategy is None:
distribution_strategy = gc['default_distribution_strategy']
result_list = [self.__class__.from_random(
'normal',
mean=mean,
std=std,
domain=result_domain,
dtype=np.complex,
distribution_strategy=self.distribution_strategy)
distribution_strategy=distribution_strategy)
for x in result_list]
# from now on extract the values from the random fields for further
......@@ -1491,7 +1495,8 @@ class Field(Loggable, Versionable, object):
temp_domain.append(repository.get('s_' + str(i), hdf5_group))
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:
new_field._val = repository.get('val', hdf5_group)
......
......@@ -44,9 +44,6 @@ class LineSearch(Loggable, object):
__metaclass__ = abc.ABCMeta
def __init__(self):
self.line_energy = None
self.f_k_minus_1 = None
self.preferred_initial_step_size = None
......
......@@ -373,6 +373,10 @@ class MPIFFT(Transform):
original_shape = inp.shape
inp = inp.reshape(inp.shape[0], 1)
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:
transform_shape = list(inp.shape)
......
......@@ -71,7 +71,8 @@ class InvertibleOperatorMixin(object):
def _times(self, x, spaces, x0=None):
if x0 is None:
x0 = Field(self.target, val=0., dtype=x.dtype)
x0 = Field(self.target, val=0., dtype=x.dtype,
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.inverse_times,
b=x,
......@@ -80,7 +81,8 @@ class InvertibleOperatorMixin(object):
def _adjoint_times(self, x, spaces, x0=None):
if x0 is None:
x0 = Field(self.domain, val=0., dtype=x.dtype)
x0 = Field(self.domain, val=0., dtype=x.dtype,
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.adjoint_inverse_times,
b=x,
......@@ -89,7 +91,8 @@ class InvertibleOperatorMixin(object):
def _inverse_times(self, x, spaces, x0=None):
if x0 is None:
x0 = Field(self.domain, val=0., dtype=x.dtype)
x0 = Field(self.domain, val=0., dtype=x.dtype,
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.times,
b=x,
......@@ -98,7 +101,8 @@ class InvertibleOperatorMixin(object):
def _adjoint_inverse_times(self, x, spaces, x0=None):
if x0 is None:
x0 = Field(self.target, val=0., dtype=x.dtype)
x0 = Field(self.target, val=0., dtype=x.dtype,
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.adjoint_times,
b=x,
......
......@@ -21,7 +21,6 @@ class FFTSmoothingOperator(SmoothingOperator):
# transform to the (global-)default codomain and perform all remaining
# steps therein
transformator = self._get_transformator(x.dtype)
transformed_x = transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]]
......
......@@ -21,3 +21,6 @@ class HealpixPlotter(PlotterBase):
def _initialize_figure(self):
return Figure2D(plots=None)
def _parse_data(self, data, field, spaces):
return data
......@@ -37,7 +37,8 @@ class Prober(object):
"""
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._distribution_strategy = \
......@@ -45,6 +46,7 @@ class Prober(object):
self._probe_count = self._parse_probe_count(probe_count)
self._random_type = self._parse_random_type(random_type)
self.compute_variance = bool(compute_variance)
self.probe_dtype = np.dtype(probe_dtype)
# ---Properties---
......@@ -104,6 +106,7 @@ class Prober(object):
""" a random-probe generator """
f = Field.from_random(random_type=self.random_type,
domain=self.domain,
dtype=self.probe_dtype,
distribution_strategy=self.distribution_strategy)
uid = np.random.randint(1e18)
return (uid, f)
......
......@@ -22,7 +22,7 @@ import numpy as np
from nifty.spaces.space import Space
from d2o import arange
from d2o import arange, distributed_data_object
class LMSpace(Space):
......@@ -138,8 +138,10 @@ class LMSpace(Space):
for l in range(1, lmax+1):
ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
idx += 2*(lmax+1-l)
dists = arange(start=0, stop=self.shape[0],
distribution_strategy=distribution_strategy)
dists = distributed_data_object(
global_shape=self.shape,
dtype=np.float,
distribution_strategy=distribution_strategy)
dists.set_local_data(ldist)
return dists
......
......@@ -16,6 +16,7 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import ast
import numpy as np
from d2o import distributed_data_object
......@@ -311,6 +312,6 @@ class PowerSpace(Space):
@classmethod
def _from_hdf5(cls, hdf5_group, repository):
hp = repository.get('harmonic_partner', hdf5_group)
exec("bb = " + hdf5_group.attrs['binbounds'])
bb = ast.literal_eval(hdf5_group.attrs['binbounds'])
ds = hdf5_group.attrs['distribution_strategy']
return PowerSpace(hp, ds, binbounds=bb)
......@@ -63,7 +63,7 @@ class FFTOperatorTests(unittest.TestCase):
assert_equal(res[zc1 * (dim1 // 2), zc2 * (dim2 // 2)], 0.)
@expand(product(["numpy", "fftw", "fftw_mpi"],
[10, 11], [False, True], [False, True],
[12, ], [False, True], [False, True],
[0.1, 1, 3.7],
[np.float64, np.complex128, np.float32, np.complex64]))
def test_fft1D(self, module, dim1, zc1, zc2, d, itp):
......@@ -86,7 +86,7 @@ class FFTOperatorTests(unittest.TestCase):
rtol=tol, atol=tol)
@expand(product(["numpy", "fftw", "fftw_mpi"],
[10, 11], [9, 12], [False, True],
[12, 15], [9, 12], [False, True],
[False, True], [False, True], [False, True], [0.1, 1, 3.7],
[0.4, 1, 2.7],
[np.float64, np.complex128, np.float32, np.complex64]))
......
......@@ -21,6 +21,7 @@ import unittest
from numpy.testing import assert_equal
from keepers import Repository
from test.common import expand, generate_spaces
from nifty import Field
from nose.plugins.skip import SkipTest
import os
......@@ -33,6 +34,12 @@ class SpaceSerializationTests(unittest.TestCase):
raise SkipTest
repo = Repository('test.h5')
repo.add(space, 'space')
field = Field(space,val=42.)
repo.add(field, 'field')
repo.commit()
assert_equal(space, repo.get('space'))
os.remove('test.h5')
assert_equal(field, repo.get('field'))
try:
os.remove('test.h5')
except OSError:
pass
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