...
 
Commits (19)
...@@ -90,7 +90,7 @@ if __name__ == '__main__': ...@@ -90,7 +90,7 @@ if __name__ == '__main__':
# Build the line-of-sight response and define signal response # Build the line-of-sight response and define signal response
LOS_starts, LOS_ends = random_los(100) if mode == 0 else radial_los(100) LOS_starts, LOS_ends = random_los(100) if mode == 0 else radial_los(100)
R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends) R = ift.LOSResponse2(position_space, starts=LOS_starts, ends=LOS_ends)
signal_response = R(signal) signal_response = R(signal)
# Specify noise # Specify noise
......
...@@ -76,6 +76,7 @@ from .plot import Plot ...@@ -76,6 +76,7 @@ from .plot import Plot
from .library.smooth_linear_amplitude import SLAmplitude, CepstrumOperator from .library.smooth_linear_amplitude import SLAmplitude, CepstrumOperator
from .library.inverse_gamma_operator import InverseGammaOperator from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse from .library.los_response import LOSResponse
from .library.los_response2 import LOSResponse2
from .library.dynamic_operator import (dynamic_operator, from .library.dynamic_operator import (dynamic_operator,
dynamic_lightcone_operator) dynamic_lightcone_operator)
from .library.light_cone_operator import LightConeOperator from .library.light_cone_operator import LightConeOperator
......
# 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/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import aslinearoperator
from scipy.special import erfc
from .. import dobj
from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..operators.linear_operator import LinearOperator
from .. import dobj
from array import array
class Pixellister(object):
def __init__(self, start, end, inc):
self._Rmax = 0.5*np.sqrt(start.shape[0]) # cell center to corner
self._start = start
self._dir = end - start
self._dirlen = np.linalg.norm(self._dir, axis=0)
self._dir *= 1./self._dirlen.reshape((1, -1))
self._inc = inc
def _losdist(self, idx, pos):
pos2 = pos.reshape((-1, 1))
didx = self._dir[:, idx]
t1 = pos2 - self._start[:, idx]
dotpr = np.sum(didx*t1, axis=0)
dotpr = np.maximum(0., np.minimum(self._dirlen[idx], dotpr))
return np.linalg.norm(t1 - didx*dotpr, axis=0)
def build_pixellist(self, idx, lo, hi, pixels):
diff = hi-lo
mid = (lo+hi-1.)*0.5
dist = self._losdist(idx, mid)
if np.all(diff == 1):
for i in range(len(dist)):
if dist[i] < self._Rmax:
pixels[0][idx[i]].append(np.sum(self._inc*lo))
pixels[1][idx[i]].append(dist[i])
else:
t1 = mid-lo
idx2 = idx[dist < np.linalg.norm(t1) + self._Rmax]
del dist
if len(idx2) > 0:
dim = np.argmax(diff)
imid = lo[dim]+diff[dim]//2
lmid = lo.copy()
lmid[dim] = imid
hmid = hi.copy()
hmid[dim] = imid
self.build_pixellist(idx2, lo, hmid, pixels)
self.build_pixellist(idx2, lmid, hi, pixels)
def _comp_traverse(start, end, shp):
nlos = start.shape[1]
ndim = start.shape[0]
inc = np.full(len(shp), 1, dtype=np.int64)
for i in range(-2, -len(shp)-1, -1):
inc[i] = inc[i+1]*shp[i+1]
pl = Pixellister(start, end, inc)
pixels = [[array("l") for _ in range(nlos)],
[array("d") for _ in range(nlos)]]
pl.build_pixellist(np.arange(nlos), np.array(shp)*0, np.array(shp), pixels)
return [(np.array(p), np.exp(-(np.array(w)**2)))
for p, w in zip(pixels[0], pixels[1])]
class LOSResponse2(LinearOperator):
"""Line-of-sight response operator
This operator transforms from a single RGSpace to an UnstructuredDomain
with as many entries as there were lines of sight passed to the
constructor. Adjoint application is also provided.
Parameters
----------
domain : RGSpace or DomainTuple
The operator's input domain. This must be a single RGSpace, with
all distances the same
starts, ends : numpy.ndarray(float) with two dimensions
Arrays containing the start and end points of the individual lines
of sight. The first dimension must have as many entries as `domain`
has dimensions. The second dimensions must be identical for both arrays
and indicated the total number of lines of sight.
Notes
-----
`starts and `ends` have to be identical on
every calling MPI task (i.e. the full LOS information has to be provided on
every task).
"""
def __init__(self, domain, starts, ends):
super(LOSResponse2, self).__init__()
self._domain = DomainTuple.make(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
if ((not isinstance(self.domain[0], RGSpace)) or
(len(self._domain) != 1)):
raise TypeError("The domain must be exactly one RGSpace instance.")
# if not np.all(self.domain[0].distances == self.domain[0].distances[0]):
# raise ValueError("domain distances must all be equal")
ndim = len(self.domain[0].shape)
starts = np.array(starts)
nlos = starts.shape[1]
ends = np.array(ends)
if starts.shape[0] != ndim:
raise TypeError("dimension mismatch")
if starts.shape != ends.shape:
raise TypeError("dimension mismatch")
self._local_shape = dobj.local_shape(self.domain[0].shape)
local_zero_point = (np.array(
dobj.ibegin_from_shape(self.domain[0].shape)) *
np.array(self.domain[0].distances))
diffs = ends-starts
difflen = np.linalg.norm(diffs, axis=0)
diffs /= difflen
lzp = local_zero_point.reshape((-1, 1))
dist = np.array(self.domain[0].distances).reshape((-1, 1))
localized_pixel_starts = (starts-lzp)/dist + 0.5
localized_pixel_ends = (ends-lzp)/dist + 0.5
# get the shape of the local data slice
w_i = _comp_traverse(localized_pixel_starts,
localized_pixel_ends,
self._local_shape)
boxsz = 16
nlos = len(w_i)
npix = np.prod(self._local_shape)
ntot = 0
for i in w_i:
ntot += len(i[1])
pri = np.empty(ntot, dtype=np.float64)
ilos = np.empty(ntot, dtype=np.int32)
iarr = np.empty(ntot, dtype=np.int32)
xwgt = np.empty(ntot, dtype=np.float32)
ofs = 0
cnt = 0
for i in w_i:
nval = len(i[1])
ilos[ofs:ofs+nval] = cnt
iarr[ofs:ofs+nval] = i[0]
xwgt[ofs:ofs+nval] = i[1]
fullidx = np.unravel_index(i[0], self._local_shape)
tmp = np.zeros(nval, dtype=np.float64)
fct = 1.
for j in range(ndim):
tmp += (fullidx[j]//boxsz)*fct
fct *= self._local_shape[j]
tmp += cnt/float(nlos)
tmp += iarr[ofs:ofs+nval]/(float(nlos)*float(npix))
pri[ofs:ofs+nval] = tmp
ofs += nval
cnt += 1
xtmp = np.argsort(pri)
ilos = ilos[xtmp]
iarr = iarr[xtmp]
xwgt = xwgt[xtmp]
self._smat = aslinearoperator(
coo_matrix((xwgt, (ilos, iarr)),
shape=(nlos, np.prod(self._local_shape))))
self._target = DomainTuple.make(UnstructuredDomain(nlos))
def apply(self, x, mode):
self._check_input(x, mode)
if mode == self.TIMES:
result_arr = self._smat.matvec(x.local_data.reshape(-1))
return Field.from_global_data(self._target, result_arr,
sum_up=True)
local_input_data = x.to_global_data().reshape(-1)
res = self._smat.rmatvec(local_input_data).reshape(self._local_shape)
return Field.from_local_data(self._domain, res)
...@@ -48,6 +48,14 @@ def testLOSResponse(sp, dtype): ...@@ -48,6 +48,14 @@ def testLOSResponse(sp, dtype):
ift.extra.consistency_check(op, dtype, dtype) ift.extra.consistency_check(op, dtype, dtype)
@pmp('sp', _p_RG_spaces)
def testLOSResponse2(sp, dtype):
starts = np.random.randn(len(sp.shape), 10)
ends = np.random.randn(len(sp.shape), 10)
op = ift.LOSResponse2(sp, starts, ends)
ift.extra.consistency_check(op, dtype, dtype)
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces) @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
def testOperatorCombinations(sp, dtype): def testOperatorCombinations(sp, dtype):
a = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype)) a = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
......