direct_smoothing_operator.py 6.29 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
9


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

        self._domain = self._parse_domain(domain)
        if len(self._domain) != 1:
            raise ValueError("DirectSmoothingOperator only accepts exactly one"
                             " space as input domain.")

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

    def _times(self, x, spaces):
Martin Reinecke's avatar
Martin Reinecke committed
25
        if self._sigma == 0:
26
27
28
29
30
31
            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
32
        return self._smooth(x, (0,) if spaces is None else spaces)
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

    # ---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---

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    def _precompute(self, x, sigma, dxmax=None):
        """ 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.

        sigma: The sigma of the Gaussian with which the function living on x
            should be smoothed, in the same units as x.
        dxmax: (optional) The maximum distance up to which smoothing is
            performed, in the same units as x. Default is 3.01*sigma.

        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.
        """

        if dxmax is None:
76
            dxmax = self._effective_smoothing_width*sigma
77
78
79
80
81
82
83
84
85

        x = np.asarray(x)

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

        wgt = []
        expfac = 1. / (2. * sigma*sigma)
        for i in range(x.size):
Theo Steininger's avatar
Theo Steininger committed
86
            if nval[i] > 0:
87
88
89
90
91
92
                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([]))
93
94
95
96
97
98
99
100
101
102
103

        return ibegin, nval, wgt

    def _apply_kernel_along_array(self, power, startindex, endindex,
                                  distances, smooth_length, smoothing_width,
                                  ibegin, nval, wgt):

        if smooth_length == 0.0:
            return power[startindex:endindex]

        p_smooth = np.zeros(endindex-startindex, dtype=power.dtype)
Martin Reinecke's avatar
Martin Reinecke committed
104
        for i in range(startindex, endindex):
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
            imin = max(startindex, ibegin[i])
            imax = min(endindex, ibegin[i]+nval[i])
            p_smooth[imin:imax] += (power[i] *
                                    wgt[i][imin-ibegin[i]:imax-imin+ibegin[i]])

        return p_smooth

    def _apply_along_axis(self, axis, arr, startindex, endindex, distances,
                          smooth_length, smoothing_width):

        nd = arr.ndim
        ibegin, nval, wgt = self._precompute(
                distances, smooth_length, smooth_length*smoothing_width)

        ind = np.zeros(nd-1, dtype=np.int)
        i = np.zeros(nd, dtype=object)
        shape = arr.shape
Martin Reinecke's avatar
Martin Reinecke committed
122
        indlist = np.asarray(list(range(nd)))
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
        indlist = np.delete(indlist, axis)
        i[axis] = slice(None, None)
        outshape = np.asarray(shape).take(indlist)

        i.put(indlist, ind)

        Ntot = np.product(outshape)
        holdshape = outshape
        slicedArr = arr[tuple(i.tolist())]
        res = self._apply_kernel_along_array(
                    slicedArr, startindex, endindex, distances,
                    smooth_length, smoothing_width, ibegin, nval, wgt)

        outshape = np.asarray(arr.shape)
        outshape[axis] = endindex - startindex
        outarr = np.zeros(outshape, dtype=arr.dtype)
        outarr[tuple(i.tolist())] = res
        k = 1
        while k < Ntot:
            # increment the index
            ind[nd-1] += 1
            n = -1
            while (ind[n] >= holdshape[n]) and (n > (1-nd)):
                ind[n-1] += 1
                ind[n] = 0
                n -= 1
            i.put(indlist, ind)
            slicedArr = arr[tuple(i.tolist())]
            res = self._apply_kernel_along_array(
                    slicedArr, startindex, endindex, distances,
                    smooth_length, smoothing_width, ibegin, nval, wgt)
            outarr[tuple(i.tolist())] = res
            k += 1

        return outarr

159
    def _smooth(self, x, spaces):
160
161
162
163
164
165
166
167
        # infer affected axes
        # we rely on the knowledge, that `spaces` is a tuple with length 1.
        affected_axes = x.domain_axes[spaces[0]]
        if len(affected_axes) > 1:
            raise ValueError("By this implementation only one-dimensional "
                             "spaces can be smoothed directly.")
        affected_axis = affected_axes[0]

Martin Reinecke's avatar
more    
Martin Reinecke committed
168
        distance_array = x.domain[spaces[0]].get_distance_array()
Martin Reinecke's avatar
Martin Reinecke committed
169
        if self._log_distances:
Martin Reinecke's avatar
Martin Reinecke committed
170
            np.log(np.maximum(distance_array, 1e-15), out=distance_array)
171

Martin Reinecke's avatar
Martin Reinecke committed
172
        return self._apply_along_axis(
Martin Reinecke's avatar
Martin Reinecke committed
173
174
175
176
177
                              affected_axis,
                              x.val,
                              startindex=0,
                              endindex=x.shape[affected_axis],
                              distances=distance_array,
Martin Reinecke's avatar
Martin Reinecke committed
178
                              smooth_length=self._sigma,
179
                              smoothing_width=self._effective_smoothing_width)