direct_smoothing_operator.py 3.77 KB
Newer Older
1
2
# -*- coding: utf8 -*-

Martin Reinecke's avatar
Martin Reinecke committed
3
4
from __future__ import division
from builtins import range
5
6
import numpy as np

7
from ..endomorphic_operator import EndomorphicOperator
8
from ... import nifty_utilities as utilities
Martin Reinecke's avatar
Martin Reinecke committed
9
from ... import Field, DomainTuple
10
11


12
class DirectSmoothingOperator(EndomorphicOperator):
13
14
    def __init__(self, domain, sigma, log_distances=False,
                 default_spaces=None):
15
16
        super(DirectSmoothingOperator, self).__init__(default_spaces)

Martin Reinecke's avatar
Martin Reinecke committed
17
        self._domain = DomainTuple.make(domain)
18
19
20
21
        if len(self._domain) != 1:
            raise ValueError("DirectSmoothingOperator only accepts exactly one"
                             " space as input domain.")

Martin Reinecke's avatar
Martin Reinecke committed
22
        self._sigma = float(sigma)
23
24
25
26
        self._log_distances = log_distances
        self._effective_smoothing_width = 3.01

    def _times(self, x, spaces):
Martin Reinecke's avatar
Martin Reinecke committed
27
        if self._sigma == 0:
28
29
30
31
32
33
            return x.copy()

        # the domain of the smoothing operator contains exactly one space.
        # Hence, if spaces is None, but we passed LinearOperator's
        # _check_input_compatibility, we know that x is also solely defined
        # on that space
Martin Reinecke's avatar
Martin Reinecke committed
34
        return self._smooth(x, (0,) if spaces is None else spaces)
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

    # ---Mandatory properties and methods---
    @property
    def domain(self):
        return self._domain

    @property
    def self_adjoint(self):
        return True

    @property
    def unitary(self):
        return False

    # ---Added properties and methods---

Martin Reinecke's avatar
Martin Reinecke committed
51
    def _precompute(self, x):
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
        """ Does precomputations for Gaussian smoothing on a 1D irregular grid.

        Parameters
        ----------
        x: 1D floating point array or list containing the individual grid
            positions. Points must be given in ascending order.


        Returns
        -------
        ibegin: integer array of the same size as x
            ibegin[i] is the minimum grid index to consider when computing the
            smoothed value at grid index i
        nval: integer array of the same size as x
            nval[i] is the number of indices to consider when computing the
            smoothed value at grid index i.
        wgt: list with the same number of entries as x
            wgt[i] is an array with nval[i] entries containing the
            normalized smoothing weights.
        """

Martin Reinecke's avatar
Martin Reinecke committed
73
        dxmax = self._effective_smoothing_width*self._sigma
74
75
76
77
78
79
80

        x = np.asarray(x)

        ibegin = np.searchsorted(x, x-dxmax)
        nval = np.searchsorted(x, x+dxmax) - ibegin

        wgt = []
Martin Reinecke's avatar
Martin Reinecke committed
81
        expfac = 1. / (2. * self._sigma*self._sigma)
82
        for i in range(x.size):
Theo Steininger's avatar
Theo Steininger committed
83
            if nval[i] > 0:
84
85
86
87
88
89
                t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
                t = np.exp(-t*t*expfac)
                t *= 1./np.sum(t)
                wgt.append(t)
            else:
                wgt.append(np.array([]))
90
91
92

        return ibegin, nval, wgt

93
    def _smooth(self, x, spaces):
94
        # infer affected axes
Martin Reinecke's avatar
Martin Reinecke committed
95
        # we rely on the knowledge that `spaces` is a tuple with length 1.
Martin Reinecke's avatar
Martin Reinecke committed
96
        affected_axes = x.domain.axes[spaces[0]]
Martin Reinecke's avatar
Martin Reinecke committed
97
        if len(affected_axes) != 1:
98
99
            raise ValueError("By this implementation only one-dimensional "
                             "spaces can be smoothed directly.")
100
        axis = affected_axes[0]
101

102
        distances = x.domain[spaces[0]].get_k_length_array()
Martin Reinecke's avatar
Martin Reinecke committed
103
        if self._log_distances:
104
105
106
            distances = np.log(np.maximum(distances, 1e-15))

        ibegin, nval, wgt = self._precompute(distances)
Martin Reinecke's avatar
Martin Reinecke committed
107
        res = Field(x.domain, dtype=x.dtype)
108
109
110
111
112
        for sl in utilities.get_slice_list(x.val.shape, (axis,)):
            inp = x.val[sl]
            out = np.zeros(inp.shape[0], dtype=inp.dtype)
            for i in range(inp.shape[0]):
                out[ibegin[i]:ibegin[i]+nval[i]] += inp[i] * wgt[i][:]
Martin Reinecke's avatar
Martin Reinecke committed
113
114
            res.val[sl] = out
        return res