Source code for nifty7.library.dynamic_operator

# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

import numpy as np

from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..operators.diagonal_operator import DiagonalOperator
from ..operators.field_zero_padder import FieldZeroPadder
from ..operators.harmonic_operators import FFTOperator
from ..operators.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import FieldAdapter, Realizer
from ..sugar import makeOp
from .light_cone_operator import LightConeOperator, _field_from_function


def _float_or_listoffloat(inp):
    return [float(x) for x in inp] if isinstance(inp, list) else float(inp)


def _make_dynamic_operator(target, harmonic_padding, sm_s0, sm_x0, cone, keys, causal,
                           minimum_phase, sigc=None, quant=None, codomain=None):
    if not isinstance(target, RGSpace):
        raise TypeError("RGSpace required")
    if not target.harmonic:
        raise TypeError("Target space must be harmonic")
    if not (isinstance(harmonic_padding, int) or harmonic_padding is None
            or all(isinstance(ii, int) for ii in harmonic_padding)):
        raise TypeError
    sm_s0 = float(sm_s0)
    sm_x0 = _float_or_listoffloat(sm_x0)
    cone = bool(cone)
    if not all(isinstance(ss, str) for ss in keys):
        raise TypeError
    causal, minimum_phase = bool(causal), bool(minimum_phase)
    if sigc is not None:
        sigc = _float_or_listoffloat(sigc)
    if quant is not None:
        quant = float(quant)
    if cone and (sigc is None or quant is None):
        raise RuntimeError

    if codomain is None:
        codomain = target.get_default_codomain()
    dom = DomainTuple.make(codomain)
    ops = {}
    FFT = FFTOperator(dom)
    Real = Realizer(dom)
    ops['FFT'] = FFT
    ops['Real'] = Real
    if harmonic_padding is None:
        CentralPadd = ScalingOperator(FFT.target, 1.)
    else:
        if isinstance(harmonic_padding, int):
            harmonic_padding = list((harmonic_padding,)*len(FFT.target.shape))
        elif len(harmonic_padding) != len(FFT.target.shape):
            raise (ValueError, "Shape mismatch!")
        shp = [
            sh + harmonic_padding[ii] for ii, sh in enumerate(FFT.target.shape)
        ]
        CentralPadd = FieldZeroPadder(FFT.target, shp, central=True)
    ops['central_padding'] = CentralPadd
    sdom = CentralPadd.target[0].get_default_codomain()
    FFTB = FFTOperator(sdom)(Realizer(sdom))

    m = FieldAdapter(sdom, keys[0])

    dists = m.target[0].distances
    if isinstance(sm_x0, float):
        sm_x0 = list((sm_x0,)*len(dists))
    elif len(sm_x0) != len(dists):
        raise (ValueError, "Shape mismatch!")

    def smoothness_prior_func(x):
        res = 1.
        for i in range(len(dists)):
            res = res + (x[i]/sm_x0[i]/dists[i])**2
        return sm_s0/res

    Sm = makeOp(_field_from_function(m.target, smoothness_prior_func))
    m = CentralPadd.adjoint(FFTB(Sm(m)))
    ops['smoothed_dynamics'] = m

    m = -m.ptw("log")
    if not minimum_phase:
        m = m.ptw("exp")
    if causal or minimum_phase:
        m = Real.adjoint(FFT.inverse(Realizer(FFT.target).adjoint(m)))
        kernel = makeOp(
            _field_from_function(FFT.domain, (lambda x: 1. + np.sign(x[0]))))
        m = kernel(m)

    if cone and len(m.target.shape) > 1:
        if isinstance(sigc, float):
            sigc = list((sigc,)*(len(m.target.shape) - 1))
        elif len(sigc) != len(m.target.shape) - 1:
            raise (ValueError, "Shape mismatch!")
        c = FieldAdapter(UnstructuredDomain(len(sigc)), keys[1])
        c = makeOp(Field(c.target, np.array(sigc)))(c)

        lightspeed = ScalingOperator(c.target, -0.5)(c).ptw("exp")
        scaling = np.array(m.target[0].distances[1:])/m.target[0].distances[0]
        scaling = DiagonalOperator(Field(c.target, scaling))
        ops['lightspeed'] = scaling(lightspeed)

        c = LightConeOperator(c.target, m.target, quant) @ c.ptw("exp")
        ops['light_cone'] = c
        m = c*m

    if causal or minimum_phase:
        m = FFT(Real(m))
    if minimum_phase:
        m = m.ptw("exp")
    return m, ops


[docs]def dynamic_operator(*, target, harmonic_padding, sm_s0, sm_x0, key, causal=True, minimum_phase=False): """Constructs an operator encoding the Green's function of a linear homogeneous dynamic system. When evaluated, this operator returns the Green's function representation in harmonic space. This result can be used as a convolution kernel to construct solutions of the homogeneous stochastic differential equation encoded in this operator. Note that if causal is True, the Green's function is convolved with a step function in time, where the temporal axis is the first axis of the space. In this case the resulting function only extends up to half the length of the first axis of the space to avoid boundary effects during convolution. If minimum_phase is true then the spectrum of the Green's function is used to construct a corresponding minimum phase filter. Parameters ---------- target : RGSpace The harmonic space in which the Green's function shall be constructed. harmonic_padding : None, int, list of int Amount of central padding in harmonic space in pixels. If None the field is not padded at all. sm_s0 : float Cutoff for dynamic smoothness prior. sm_x0 : float, list of float Scaling of dynamic smoothness along each axis. key : String key for dynamics encoding parameter. causal : boolean Whether or not the Green's function shall be causal in time. Default is True. minimum_phase: boolean Whether or not the Green's function shall be a minimum phase filter. Default is False. Returns ------- Operator The Operator encoding the dynamic Green's function in target space. Dictionary of Operator A collection of sub-chains of Operator which can be used for plotting and evaluation. Notes ----- The first axis of the domain is interpreted the time axis. """ dct = { 'target': target, 'harmonic_padding': harmonic_padding, 'sm_s0': sm_s0, 'sm_x0': sm_x0, 'keys': [key], 'causal': causal, 'cone': False, 'minimum_phase': minimum_phase, } return _make_dynamic_operator(**dct)
[docs]def dynamic_lightcone_operator(*, target, harmonic_padding, sm_s0, sm_x0, key, lightcone_key, sigc, quant, causal=True, minimum_phase=False): '''Extends the functionality of :func:`dynamic_operator` to a Green's function which is constrained to be within a light cone. The resulting Green's function is constrained to be within a light cone. This is achieved via convolution of the function with a light cone in space-time. Thereby the first axis of the space is set to be the teporal axis. Parameters ---------- target : RGSpace The harmonic space in which the Green's function shall be constructed. It needs to have at least two dimensions. harmonic_padding : None, int, list of int Amount of central padding in harmonic space in pixels. If None the field is not padded at all. sm_s0 : float Cutoff for dynamic smoothness prior. sm_x0 : float, list of float Scaling of dynamic smoothness along each axis. key : String Key for dynamics encoding parameter. lightcone_key: String Key for lightspeed paramteter. sigc : float, list of float Variance of lightspeed parameter. quant : float Quantization of the light cone in pixels. causal : boolean Whether or not the Green's function shall be causal in time. Default is True. minimum_phase: boolean Whether or not the Green's function shall be a minimum phase filter. Default is False. Returns ------- Operator The Operator encoding the dynamic Green's function in harmonic space. Dictionary of Operator A collection of sub-chains of Operator which can be used for plotting and evaluation. Notes ----- The first axis of the domain is interpreted the time axis. ''' if len(target.shape) < 2: raise ValueError("Space must be at least 2 dimensional!") dct = { 'target': target, 'harmonic_padding': harmonic_padding, 'sm_s0': sm_s0, 'sm_x0': sm_x0, 'keys': [key, lightcone_key], 'causal': causal, 'cone': True, 'minimum_phase': minimum_phase, 'sigc': sigc, 'quant': quant } return _make_dynamic_operator(**dct)