Skip to content

GitLab

N
NIFTy
Commit fd98cba2 authored by Theo Steininger's avatar Theo Steininger
Browse files
Options

Integrated real harmonic RGSpace representation deeply into NIFTy.

parent 39b6628e
Pipeline #16544 failed with stage
    in 6 minutes and 17 seconds
    Hide whitespace changes
    Inline Side-by-side
    Showing with 158 additions and 583 deletions
    demos/wiener_filter_via_curvature.py
    View file @ fd98cba2
    ......@@ -2,13 +2,14 @@ import numpy as np
    from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
    DiagonalOperator, ResponseOperator, plotting,\
    create_power_operator
    create_power_operator, nifty_configuration
    from nifty.library import WienerFilterCurvature
    if __name__ == "__main__":
    distribution_strategy = 'not'
    nifty_configuration['default_distribution_strategy'] = 'fftw'
    nifty_configuration['harmonic_rg_base'] = 'real'
    # Setting up variable parameters
    ......@@ -36,19 +37,17 @@ if __name__ == "__main__":
    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)
    fft = FFTOperator(harmonic_space, target=signal_space)
    power_space = PowerSpace(harmonic_space)
    # Creating the mock data
    S = create_power_operator(harmonic_space, power_spectrum=power_spectrum,
    distribution_strategy=distribution_strategy)
    S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
    mock_power = Field(power_space, val=power_spectrum,
    distribution_strategy=distribution_strategy)
    mock_power = Field(power_space, val=power_spectrum)
    np.random.seed(43)
    mock_harmonic = mock_power.power_synthesize(real_signal=True)
    if nifty_configuration['harmonic_rg_base'] == 'real':
    mock_harmonic = mock_harmonic.real
    mock_signal = fft(mock_harmonic)
    R = ResponseOperator(signal_space, sigma=(response_sigma,))
    ......@@ -74,9 +73,10 @@ if __name__ == "__main__":
    plotter = plotting.RG2DPlotter()
    plotter.path = 'mock_signal.html'
    plotter(mock_signal)
    plotter(mock_signal.real)
    plotter.path = 'data.html'
    plotter(Field(signal_space,
    val=data.val.get_full_data().reshape(signal_space.shape)))
    plotter(Field(
    signal_space,
    val=data.val.get_full_data().real.reshape(signal_space.shape)))
    plotter.path = 'map.html'
    plotter(m_s)
    plotter(m_s.real)
    demos/wiener_filter_via_curvature_real.py deleted 100644 → 0
    View file @ 39b6628e
    import numpy as np
    from nifty import RGSpace, PowerSpace, Field, RealFFTOperator,\
    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.05
    # 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.01
    # 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 = RealFFTOperator.get_default_codomain(signal_space)
    fft = RealFFTOperator(harmonic_space, target=signal_space)
    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_harmonic = mock_harmonic.real + mock_harmonic.imag
    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.path = 'mock_signal.html'
    plotter(mock_signal)
    plotter.path = 'data.html'
    plotter(Field(signal_space,
    val=data.val.get_full_data().reshape(signal_space.shape)))
    plotter.path = 'map.html'
    plotter(m_s)
    nifty/config/nifty_config.py
    View file @ fd98cba2
    ......@@ -70,11 +70,18 @@ variable_default_distribution_strategy = keepers.Variable(
    if z == 'fftw' else True),
    genus='str')
    variable_harmonic_rg_base = keepers.Variable(
    'harmonic_rg_base',
    ['real', 'complex'],
    lambda z: z in ['real', 'complex'],
    genus='str')
    nifty_configuration = keepers.get_Configuration(
    name='NIFTy',
    variables=[variable_fft_module,
    variable_default_field_dtype,
    variable_default_distribution_strategy],
    variable_default_distribution_strategy,
    variable_harmonic_rg_base],
    file_name='NIFTy.conf',
    search_paths=[os.path.expanduser('~') + "/.config/nifty/",
    os.path.expanduser('~') + "/.config/",
    ......
    nifty/field.py
    View file @ fd98cba2
    ......@@ -19,7 +19,6 @@
    from __future__ import division
    import ast
    import itertools
    import numpy as np
    from keepers import Versionable,\
    ......
    nifty/operators/fft_operator/__init__.py
    View file @ fd98cba2
    ......@@ -18,4 +18,3 @@
    from transformations import *
    from .fft_operator import FFTOperator
    from .real_fft_operator import RealFFTOperator
    nifty/operators/fft_operator/fft_operator.py
    View file @ fd98cba2
    ......@@ -142,17 +142,11 @@ class FFTOperator(LinearOperator):
    backward_class, self.target[0], self.domain[0], module=module)
    # Store the dtype information
    if domain_dtype is None:
    self.logger.info("Setting domain_dtype to np.complex.")
    self.domain_dtype = np.complex
    else:
    self.domain_dtype = np.dtype(domain_dtype)
    self.domain_dtype = \
    None if domain_dtype is None else np.dtype(domain_dtype)
    if target_dtype is None:
    self.logger.info("Setting target_dtype to np.complex.")
    self.target_dtype = np.complex
    else:
    self.target_dtype = np.dtype(target_dtype)
    self.target_dtype = \
    None if target_dtype is None else np.dtype(target_dtype)
    def _times(self, x, spaces):
    spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
    ......@@ -172,8 +166,10 @@ class FFTOperator(LinearOperator):
    result_domain = list(x.domain)
    result_domain[spaces[0]] = self.target[0]
    result_dtype = \
    new_val.dtype if self.target_dtype is None else self.target_dtype
    result_field = x.copy_empty(domain=result_domain,
    dtype=self.target_dtype)
    dtype=result_dtype)
    result_field.set_val(new_val=new_val, copy=True)
    return result_field
    ......@@ -196,8 +192,11 @@ class FFTOperator(LinearOperator):
    result_domain = list(x.domain)
    result_domain[spaces[0]] = self.domain[0]
    result_dtype = \
    new_val.dtype if self.domain_dtype is None else self.domain_dtype
    result_field = x.copy_empty(domain=result_domain,
    dtype=self.domain_dtype)
    dtype=result_dtype)
    result_field.set_val(new_val=new_val, copy=True)
    return result_field
    ......
    nifty/operators/fft_operator/real_fft_operator.py deleted 100644 → 0
    View file @ 39b6628e
    # 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
    import nifty.nifty_utilities as utilities
    from nifty.spaces import RGSpace,\
    GLSpace,\
    HPSpace,\
    LMSpace
    from nifty.operators.linear_operator import LinearOperator
    from transformations import RGRGTransformation,\
    LMGLTransformation,\
    LMHPTransformation,\
    GLLMTransformation,\
    HPLMTransformation,\
    TransformationCache
    class RealFFTOperator(LinearOperator):
    """Transforms between a pair of position and harmonic domains.
    Built-in domain pairs are
    - a harmonic and a non-harmonic RGSpace (with matching distances)
    - a HPSpace and a LMSpace
    - a GLSpace and a LMSpace
    Within a domain pair, both orderings are possible.
    The operator provides a "times" and an "adjoint_times" operation.
    For a pair of RGSpaces, the "adjoint_times" operation is equivalent to
    "inverse_times"; for the sphere-related domains this is not the case, since
    the operator matrix is not square.
    In contrast to the FFTOperator, RealFFTOperator accepts and returns
    real-valued fields only. For the harmonic-space counterpart of a
    real-valued field living on an RGSpace, the sum of real
    and imaginary components is stored. Since the full complex field has
    Hermitian symmetry, this is sufficient to reconstruct the full field
    whenever needed (e.g. during the transform back to position space).
    Parameters
    ----------
    domain: Space or single-element tuple of Spaces
    The domain of the data that is input by "times" and output by
    "adjoint_times".
    target: Space or single-element tuple of Spaces (optional)
    The domain of the data that is output by "times" and input by
    "adjoint_times".
    If omitted, a co-domain will be chosen automatically.
    Whenever "domain" is an RGSpace, the codomain (and its parameters) are
    uniquely determined.
    For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
    co-domain is chosen that should work satisfactorily in most situations,
    but for full control, the user should explicitly specify a codomain.
    module: String (optional)
    Software module employed for carrying out the transform operations.
    For RGSpace pairs this can be "scalar" or "mpi", where "scalar" is
    always available (using pyfftw if available, else numpy.fft), and "mpi"
    requires pyfftw and offers MPI parallelization.
    For sphere-related domains, only "pyHealpix" is
    available. If omitted, "fftw" is selected for RGSpaces if available,
    else "numpy"; on the sphere the default is "pyHealpix".
    Attributes
    ----------
    domain: Tuple of Spaces (with one entry)
    The domain of the data that is input by "times" and output by
    "adjoint_times".
    target: Tuple of Spaces (with one entry)
    The domain of the data that is output by "times" and input by
    "adjoint_times".
    unitary: bool
    Returns True if the operator is unitary (currently only the case if
    the domain and codomain are RGSpaces), else False.
    Raises
    ------
    ValueError:
    if "domain" or "target" are not of the proper type.
    """
    default_codomain_dictionary = {RGSpace: RGSpace,
    HPSpace: LMSpace,
    GLSpace: LMSpace,
    LMSpace: GLSpace,
    }
    transformation_dictionary = {(RGSpace, RGSpace): RGRGTransformation,
    (HPSpace, LMSpace): HPLMTransformation,
    (GLSpace, LMSpace): GLLMTransformation,
    (LMSpace, HPSpace): LMHPTransformation,
    (LMSpace, GLSpace): LMGLTransformation
    }
    def __init__(self, domain, target=None, module=None,
    default_spaces=None):
    super(RealFFTOperator, self).__init__(default_spaces)
    # Initialize domain and target
    self._domain = self._parse_domain(domain)
    if len(self.domain) != 1:
    raise ValueError("TransformationOperator accepts only exactly one "
    "space as input domain.")
    if target is None:
    target = (self.get_default_codomain(self.domain[0]), )
    self._target = self._parse_domain(target)
    if len(self.target) != 1:
    raise ValueError("TransformationOperator accepts only exactly one "
    "space as output target.")
    # Create transformation instances
    forward_class = self.transformation_dictionary[
    (self.domain[0].__class__, self.target[0].__class__)]
    backward_class = self.transformation_dictionary[
    (self.target[0].__class__, self.domain[0].__class__)]
    self._forward_transformation = TransformationCache.create(
    forward_class, self.domain[0], self.target[0], module=module)
    self._backward_transformation = TransformationCache.create(
    backward_class, self.target[0], self.domain[0], module=module)
    def _prep(self, x, spaces, dom):
    assert issubclass(x.dtype.type,np.floating), \
    "Argument must be real-valued"
    spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
    if spaces is None:
    # this case means that x lives on only one space, which is
    # identical to the space in the domain of `self`. Otherwise the
    # input check of LinearOperator would have failed.
    axes = x.domain_axes[0]
    else:
    axes = x.domain_axes[spaces[0]]
    if spaces is None:
    result_domain = dom
    else:
    result_domain = list(x.domain)
    result_domain[spaces[0]] = dom[0]
    result_field = x.copy_empty(domain=result_domain, dtype=x.dtype)
    return spaces, axes, result_field
    def _times(self, x, spaces):
    spaces, axes, result_field = self._prep(x, spaces, self.target)
    if type(self._domain[0]) != RGSpace:
    new_val = self._forward_transformation.transform(x.val, axes=axes)
    result_field.set_val(new_val=new_val, copy=True)
    return result_field
    if self._target[0].harmonic: # going to harmonic space
    new_val = self._forward_transformation.transform(x.val, axes=axes)
    result_field.set_val(new_val=new_val.real+new_val.imag)
    else:
    tval = self._domain[0].hermitianize_inverter(x.val, axes)
    tval = 0.5*((x.val+tval)+1j*(x.val-tval))
    new_val = self._forward_transformation.transform(tval, axes=axes)
    result_field.set_val(new_val=new_val.real)
    return result_field
    def _adjoint_times(self, x, spaces):
    spaces, axes, result_field = self._prep(x, spaces, self.domain)
    if type(self._domain[0]) != RGSpace:
    new_val = self._backward_transformation.transform(x.val, axes=axes)
    result_field.set_val(new_val=new_val, copy=True)
    return result_field
    if self._domain[0].harmonic: # going to harmonic space
    new_val = self._backward_transformation.transform(x.val, axes=axes)
    result_field.set_val(new_val=new_val.real+new_val.imag)
    else:
    tval = self._target[0].hermitianize_inverter(x.val, axes)
    tval = 0.5*((x.val+tval)+1j*(x.val-tval))
    new_val = self._backward_transformation.transform(tval, axes=axes)
    result_field.set_val(new_val=new_val.real)
    return result_field
    # ---Mandatory properties and methods---
    @property
    def domain(self):
    return self._domain
    @property
    def target(self):
    return self._target
    @property
    def unitary(self):
    return (self._forward_transformation.unitary and
    self._backward_transformation.unitary)
    # ---Added properties and methods---
    @classmethod
    def get_default_codomain(cls, domain):
    """Returns a codomain to the given domain.
    Parameters
    ----------
    domain: Space
    An instance of RGSpace, HPSpace, GLSpace or LMSpace.
    Returns
    -------
    target: Space
    A (more or less perfect) counterpart to "domain" with respect
    to a FFT operation.
    Whenever "domain" is an RGSpace, the codomain (and its parameters)
    are uniquely determined.
    For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
    co-domain is chosen that should work satisfactorily in most
    situations. For full control however, the user should not rely on
    this method.
    Raises
    ------
    ValueError:
    if no default codomain is defined for "domain".
    """
    domain_class = domain.__class__
    try:
    codomain_class = cls.default_codomain_dictionary[domain_class]
    except KeyError:
    raise ValueError("Unknown domain")
    try:
    transform_class = cls.transformation_dictionary[(domain_class,
    codomain_class)]
    except KeyError:
    raise ValueError(
    "No transformation for domain-codomain pair found.")
    return transform_class.get_codomain(domain)
    nifty/operators/fft_operator/transformations/rg_transforms.py
    View file @ fd98cba2
    ......@@ -260,7 +260,7 @@ class MPIFFT(Transform):
    p()
    if p.has_output:
    result = p.output_array
    result = p.output_array.copy()
    if result.shape != val.shape:
    raise ValueError("Output shape is different than input shape. "
    "Maybe fftw tries to optimize the "
    ......
    nifty/operators/fft_operator/transformations/rgrgtransformation.py
    View file @ fd98cba2
    ......@@ -43,6 +43,8 @@ class RGRGTransformation(Transformation):
    else:
    raise ValueError('Unsupported FFT module:' + module)
    self.harmonic_base = nifty_configuration['harmonic_rg_base']
    # ---Mandatory properties and methods---
    @property
    ......@@ -144,7 +146,31 @@ class RGRGTransformation(Transformation):
    val = self._transform.domain.weight(val, power=1, axes=axes)
    # Perform the transformation
    Tval = self._transform.transform(val, axes, **kwargs)
    if self.harmonic_base == 'complex':
    Tval = self._transform.transform(val, axes, **kwargs)
    else:
    if issubclass(val.dtype.type, np.complexfloating):
    Tval_real = self._transform.transform(val.real, axes,
    **kwargs)
    Tval_imag = self._transform.transform(val.imag, axes,
    **kwargs)
    if self.codomain.harmonic:
    Tval_real.data.real += Tval_real.data.imag
    Tval_real.data.imag = \
    Tval_imag.data.real + Tval_imag.data.imag
    else:
    Tval_real.data.real -= Tval_real.data.imag
    Tval_real.data.imag = \
    Tval_imag.data.real - Tval_imag.data.imag
    Tval = Tval_real
    else:
    Tval = self._transform.transform(val, axes, **kwargs)
    if self.codomain.harmonic:
    Tval.data.real += Tval.data.imag
    else:
    Tval.data.real -= Tval.data.imag
    Tval = Tval.real
    if not self._transform.codomain.harmonic:
    # correct for inverse fft.
    ......
    nifty/spaces/rg_space/rg_space.py
    View file @ fd98cba2
    ......@@ -36,7 +36,7 @@ from d2o import distributed_data_object,\
    STRATEGIES as DISTRIBUTION_STRATEGIES
    from nifty.spaces.space import Space
    from nifty.config import nifty_configuration
    class RGSpace(Space):
    """
    ......@@ -122,33 +122,36 @@ class RGSpace(Space):
    # return fixed_points
    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()
    # flip in the desired directions
    for k in range(len(axes)):
    i = axes[k]
    slice_picker = slice_primitive[:]
    slice_inverter = slice_primitive[:]
    if (not self.zerocenter[k]) or self.shape[k] % 2 == 0:
    slice_picker[i] = slice(1, None, None)
    slice_inverter[i] = slice(None, 0, -1)
    else:
    slice_picker[i] = slice(None)
    slice_inverter[i] = slice(None, None, -1)
    slice_picker = tuple(slice_picker)
    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
    if nifty_configuration['harmonic_rg_base'] == 'real':
    return x
    else:
    # 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()
    # flip in the desired directions
    for k in range(len(axes)):
    i = axes[k]
    slice_picker = slice_primitive[:]
    slice_inverter = slice_primitive[:]
    if (not self.zerocenter[k]) or self.shape[k] % 2 == 0:
    slice_picker[i] = slice(1, None, None)
    slice_inverter[i] = slice(None, 0, -1)
    else:
    slice_picker[i] = slice(None)
    slice_inverter[i] = slice(None, None, -1)
    slice_picker = tuple(slice_picker)
    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---
    ......
    nifty/sugar.py
    View file @ fd98cba2
    ......@@ -16,6 +16,8 @@
    # 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 nifty import Space,\
    PowerSpace,\
    Field,\
    ......@@ -71,9 +73,10 @@ def create_power_operator(domain, power_spectrum, dtype=None,
    distribution_strategy='not')
    f = fp.power_synthesize(mean=1, std=0, real_signal=False,
    distribution_strategy=distribution_strategy)
    # MR FIXME: we need the real part here. Could this also be achieved
    # by setting real_signal=True in the call above?
    f = f.real
    if not issubclass(fp.dtype.type, np.complexfloating):