Commit 2517d766 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge master; currently broken

parent 5cc52382
# 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
# 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 <>.
# 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 import RGSpace
from 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:
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.
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.
`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)) *
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,
boxsz = 16
nlos = len(w_i)
npix =
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)),
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,
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)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment