direct_smoothing_operator.py 5.46 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---

Martin Reinecke's avatar
Martin Reinecke committed
49
    def _precompute(self, x):
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
        """ 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
71
        dxmax = self._effective_smoothing_width*self._sigma
72
73
74
75
76
77
78

        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
79
        expfac = 1. / (2. * self._sigma*self._sigma)
80
        for i in range(x.size):
Theo Steininger's avatar
Theo Steininger committed
81
            if nval[i] > 0:
82
83
84
85
86
87
                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([]))
88
89
90
91
92
93
94

        return ibegin, nval, wgt

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

        p_smooth = np.zeros(endindex-startindex, dtype=power.dtype)
Martin Reinecke's avatar
Martin Reinecke committed
95
        for i in range(startindex, endindex):
96
97
98
99
100
101
102
            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

Martin Reinecke's avatar
Martin Reinecke committed
103
    def _apply_along_axis(self, axis, arr, startindex, endindex, distances):
104
105

        nd = arr.ndim
Martin Reinecke's avatar
Martin Reinecke committed
106
        ibegin, nval, wgt = self._precompute(distances)
107
108
109
110

        ind = np.zeros(nd-1, dtype=np.int)
        i = np.zeros(nd, dtype=object)
        shape = arr.shape
Martin Reinecke's avatar
Martin Reinecke committed
111
        indlist = np.asarray(list(range(nd)))
112
113
114
115
116
117
118
119
120
121
        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(
Martin Reinecke's avatar
Martin Reinecke committed
122
                    slicedArr, startindex, endindex, ibegin, nval, wgt)
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139

        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(
Martin Reinecke's avatar
Martin Reinecke committed
140
                    slicedArr, startindex, endindex, ibegin, nval, wgt)
141
142
143
144
145
            outarr[tuple(i.tolist())] = res
            k += 1

        return outarr

146
    def _smooth(self, x, spaces):
147
        # infer affected axes
Martin Reinecke's avatar
Martin Reinecke committed
148
        # we rely on the knowledge that `spaces` is a tuple with length 1.
149
        affected_axes = x.domain_axes[spaces[0]]
Martin Reinecke's avatar
Martin Reinecke committed
150
        if len(affected_axes) != 1:
151
152
153
154
            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
155
        distance_array = x.domain[spaces[0]].get_distance_array()
Martin Reinecke's avatar
Martin Reinecke committed
156
        if self._log_distances:
Martin Reinecke's avatar
Martin Reinecke committed
157
            np.log(np.maximum(distance_array, 1e-15), out=distance_array)
158

Martin Reinecke's avatar
Martin Reinecke committed
159
        return self._apply_along_axis(
Martin Reinecke's avatar
Martin Reinecke committed
160
161
162
163
                              affected_axis,
                              x.val,
                              startindex=0,
                              endindex=x.shape[affected_axis],
Martin Reinecke's avatar
Martin Reinecke committed
164
                              distances=distance_array)