regridding_operator.py 3.77 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
# 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17
18
19
20

import numpy as np

from ..domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..domains.rg_space import RGSpace
22
from ..field import Field
Martin Reinecke's avatar
Martin Reinecke committed
23
from ..utilities import infer_space, special_add_at
24
25
26
27
from .linear_operator import LinearOperator


class RegriddingOperator(LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
28
    """Linearly interpolates an RGSpace to an RGSpace with coarser resolution.
Philipp Arras's avatar
Docs    
Philipp Arras committed
29
30
31
32
33
34
35
36
37
38
39

    Parameters
    ----------
    domain : Domain, DomainTuple or tuple of Domain
        domain[space] needs to be an :class:`RGSpace`.
    new_shape : tuple of int
        Shape of the space which domain[space] is replaced by. Each entry must
        be smaller or equal to the respective entry in `domain[space].shape`.
    space : int
        Index of space in `domain` on which the operator shall act.
        Default is 0.
40
    """
Martin Reinecke's avatar
Martin Reinecke committed
41
    def __init__(self, domain, new_shape, space=0):
42
        self._domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
43
44
        self._space = infer_space(self._domain, space)
        dom = self._domain[self._space]
45

Martin Reinecke's avatar
Martin Reinecke committed
46
47
48
49
50
51
        if not isinstance(dom, RGSpace):
            raise TypeError("RGSpace required")
        if len(new_shape) != len(dom.shape):
            raise ValueError("Shape mismatch")
        if any([a > b for a, b in zip(new_shape, dom.shape)]):
            raise ValueError("New shape must not be larger than old shape")
Philipp Arras's avatar
Docs    
Philipp Arras committed
52
53
        if any([ii <= 0 for ii in new_shape]):
            raise ValueError('New shape must not be zero or negative.')
54

Martin Reinecke's avatar
Martin Reinecke committed
55
56
        newdist = tuple(dom.distances[i]*dom.shape[i]/new_shape[i]
                        for i in range(len(dom.shape)))
57

Martin Reinecke's avatar
Martin Reinecke committed
58
59
60
61
62
        tgt = RGSpace(new_shape, newdist)
        self._target = list(self._domain)
        self._target[self._space] = tgt
        self._target = DomainTuple.make(self._target)
        self._capability = self.TIMES | self.ADJOINT_TIMES
63

Martin Reinecke's avatar
Martin Reinecke committed
64
65
66
67
68
        ndim = len(new_shape)
        self._bindex = [None] * ndim
        self._frac = [None] * ndim
        for d in range(ndim):
            tmp = np.arange(new_shape[d])*(newdist[d]/dom.distances[d])
Martin Reinecke's avatar
Martin Reinecke committed
69
            self._bindex[d] = np.minimum(dom.shape[d]-2, tmp.astype(np.int))
Martin Reinecke's avatar
bug fix    
Martin Reinecke committed
70
            self._frac[d] = tmp-self._bindex[d]
71
72
73

    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
74
        v = x.val
Martin Reinecke's avatar
Martin Reinecke committed
75
76
        ndim = len(self.target.shape)
        curshp = list(self._dom(mode).shape)
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
77
        tgtshp = self._tgt(mode).shape
Martin Reinecke's avatar
Martin Reinecke committed
78
79
80
81
        d0 = self._target.axes[self._space][0]
        for d in self._target.axes[self._space]:
            idx = (slice(None),) * d
            wgt = self._frac[d-d0].reshape((1,)*d + (-1,) + (1,)*(ndim-d-1))
82

Martin Reinecke's avatar
Martin Reinecke committed
83
            if mode == self.ADJOINT_TIMES:
Martin Reinecke's avatar
stage 1    
Martin Reinecke committed
84
                shp = list(v.shape)
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
85
                shp[d] = tgtshp[d]
Martin Reinecke's avatar
stage 1    
Martin Reinecke committed
86
87
88
                xnew = np.zeros(shp, dtype=v.dtype)
                xnew = special_add_at(xnew, d, self._bindex[d-d0], v*(1.-wgt))
                xnew = special_add_at(xnew, d, self._bindex[d-d0]+1, v*wgt)
Martin Reinecke's avatar
Martin Reinecke committed
89
            else:  # TIMES
Martin Reinecke's avatar
stage 1    
Martin Reinecke committed
90
91
                xnew = v[idx + (self._bindex[d-d0],)] * (1.-wgt)
                xnew += v[idx + (self._bindex[d-d0]+1,)] * wgt
92

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
93
            curshp[d] = xnew.shape[d]
Martin Reinecke's avatar
Martin Reinecke committed
94
            v = xnew
Martin Reinecke's avatar
stage 1    
Martin Reinecke committed
95
        return Field(self._tgt(mode), xnew)