direct_smoothing_operator.py 7.05 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
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
        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.")

        self._sigma = sigma
        self._log_distances = log_distances
        self._effective_smoothing_width = 3.01

    def _times(self, x, spaces):
        if self.sigma == 0:
            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
        if spaces is None:
            spaces = (0,)

        return self._smooth(x, spaces)

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

    @property
    def sigma(self):
        return self._sigma

    @property
    def log_distances(self):
        return self._log_distances
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

    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:
87
            dxmax = self._effective_smoothing_width*sigma
88
89
90
91
92
93
94
95
96

        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
97
            if nval[i] > 0:
98
99
100
101
102
103
                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([]))
104
105
106
107
108
109
110
111
112
113
114

        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
115
        for i in range(startindex, endindex):
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
            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
        if axis < 0:
            axis += nd
        if (axis >= nd):
            raise ValueError(
                "axis must be less than arr.ndim; axis=%d, rank=%d."
                % (axis, nd))
        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
139
        indlist = np.asarray(list(range(nd)))
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
        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

176
    def _smooth(self, x, spaces):
177
178
179
180
181
182
183
184
185
186
        # 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
187
        distance_array = x.domain[spaces[0]].get_distance_array()
188

Martin Reinecke's avatar
Martin Reinecke committed
189
        #MR FIXME: this causes calls of log(0.) which should probably be avoided
190
        if self.log_distances:
191
            np.log(np.maximum(distance_array,1e-15), out=distance_array)
192

Martin Reinecke's avatar
more    
Martin Reinecke committed
193
194
195
196
        augmented_data = x.val
        augmented_distance_array = distance_array
        true_start = 0
        true_end = x.shape[affected_axis]
197
198
199
200
201
202
203
204
205
206
207

        # perform the convolution along the affected axes
        # currently only one axis is supported
        data_axis = affected_axes[0]

        local_result = self._apply_along_axis(
                              data_axis,
                              augmented_data,
                              startindex=true_start,
                              endindex=true_end,
                              distances=augmented_distance_array,
208
                              smooth_length=self.sigma,
209
                              smoothing_width=self._effective_smoothing_width)
210

Martin Reinecke's avatar
more    
Martin Reinecke committed
211
        return local_result