Commit 8f2f7a17 authored by Reimar H Leike's avatar Reimar H Leike
Browse files

renamed and reinterpreted the arguments to be more intuitive

parent 853c9dfe
...@@ -32,7 +32,7 @@ def _gaussian_sf(x): ...@@ -32,7 +32,7 @@ def _gaussian_sf(x):
return 0.5*erfc(x/np.sqrt(2.)) return 0.5*erfc(x/np.sqrt(2.))
def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf): def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
ndim = start.shape[0] ndim = start.shape[0]
nlos = start.shape[1] nlos = start.shape[1]
inc = np.full(len(shp), 1, dtype=np.int64) inc = np.full(len(shp), 1, dtype=np.int64)
...@@ -85,19 +85,18 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf): ...@@ -85,19 +85,18 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
cdist *= corfac cdist *= corfac
wgt = np.diff(cdist) wgt = np.diff(cdist)
mdist = 0.5*(cdist[:-1]+cdist[1:]) mdist = 0.5*(cdist[:-1]+cdist[1:])
wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], erf) wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], sig[i], erf)
add = np.append(pos1, add) add = np.append(pos1, add)
add = np.cumsum(add) add = np.cumsum(add)
out[i] = (add, wgt) out[i] = (add, wgt)
return out return out
def apply_erf(wgt, dist, lo, mid, hi, erf): def apply_erf(wgt, dist, lo, mid, hi, sig, erf):
wgt = wgt.copy() wgt = wgt.copy()
mask = dist > hi mask = dist > hi
wgt[mask] = 0. wgt[mask] = 0.
mask = (dist > lo) & (dist <= hi) mask = (dist > lo) & (dist <= hi)
sig = (1/lo - 1/mid)/3
wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig) wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig)
return wgt return wgt
...@@ -118,22 +117,32 @@ class LOSResponse(LinearOperator): ...@@ -118,22 +117,32 @@ class LOSResponse(LinearOperator):
of sight. The first dimension must have as many entries as `domain` of sight. The first dimension must have as many entries as `domain`
has dimensions. The second dimensions must be identical for both arrays has dimensions. The second dimensions must be identical for both arrays
and indicated the total number of lines of sight. and indicated the total number of lines of sight.
sigmas_low, sigmas_up : numpy.ndarray(float) (optional) sigmas: numpy.ndarray(float) (optional)
sigmas_low is 1/parallax-1/(parallax+3*parallax_error), where the parallax If this is not None, the inverse of the lengths of the LOSs are assumed to be
error is assumed to be Gaussian distributed. Gaussian diastributed with these sigmas. The start point will remain the same,
sigmas_up is the distance at which the weight is truncated. but the endpoint is assumed to be unknown.
Should be at least 1/(parallax-3*parallax_error)-1/parallax, This is a typical statistical model for astrophysical parallaxes.
but could be higher. The LOS response then returns the expected integral
over the input given that the length of the LOS is unknown and therefore the
result is averaged over different endpoints.
default: None
truncation: float (optional)
Use only if the sigmas keyword argument is used!
This truncates the probability of the endpoint lying more sigmas away than
the truncation. Used to speed up computation and to avoid negative distances.
It should hold that 1./(1./length-sigma*truncation)>0 for all lengths of the
LOSs and all corresponding sigma of sigmas.
If unsure, leave blank. If unsure, leave blank.
default: 3.
Notes Notes
----- -----
`starts, `ends`, `sigmas_low`, and `sigmas_up` have to be identical on `starts, `ends`, `sigmas`, and `truncation` have to be identical on
every calling MPI task (i.e. the full LOS information has to be provided on every calling MPI task (i.e. the full LOS information has to be provided on
every task). every task).
""" """
def __init__(self, domain, starts, ends, sigmas_low=None, sigmas_up=None): def __init__(self, domain, starts, ends, sigmas=None, truncation=3.):
self._domain = DomainTuple.make(domain) self._domain = DomainTuple.make(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES self._capability = self.TIMES | self.ADJOINT_TIMES
...@@ -145,20 +154,15 @@ class LOSResponse(LinearOperator): ...@@ -145,20 +154,15 @@ class LOSResponse(LinearOperator):
starts = np.array(starts) starts = np.array(starts)
nlos = starts.shape[1] nlos = starts.shape[1]
ends = np.array(ends) ends = np.array(ends)
if sigmas_low is None: if sigmas is None:
sigmas_low = np.zeros(nlos, dtype=np.float32) sigmas = np.zeros(nlos, dtype=np.float32)
if sigmas_up is None: sigmas = np.array(sigmas)
sigmas_up = np.zeros(nlos, dtype=np.float32)
sigmas_low = np.array(sigmas_low)
sigmas_up = np.array(sigmas_up)
if starts.shape[0] != ndim: if starts.shape[0] != ndim:
raise TypeError("dimension mismatch") raise TypeError("dimension mismatch")
if nlos != sigmas_low.shape[0]: if nlos != sigmas.shape[0]:
raise TypeError("dimension mismatch") raise TypeError("dimension mismatch")
if starts.shape != ends.shape: if starts.shape != ends.shape:
raise TypeError("dimension mismatch") raise TypeError("dimension mismatch")
if sigmas_low.shape != sigmas_up.shape:
raise TypeError("dimension mismatch")
self._local_shape = dobj.local_shape(self.domain[0].shape) self._local_shape = dobj.local_shape(self.domain[0].shape)
local_zero_point = (np.array( local_zero_point = (np.array(
...@@ -168,7 +172,10 @@ class LOSResponse(LinearOperator): ...@@ -168,7 +172,10 @@ class LOSResponse(LinearOperator):
diffs = ends-starts diffs = ends-starts
difflen = np.linalg.norm(diffs, axis=0) difflen = np.linalg.norm(diffs, axis=0)
diffs /= difflen diffs /= difflen
real_ends = ends + sigmas_up*diffs real_distances = 1./(1./difflen - truncation*sigmas)
if np.any(real_distances<0):
raise ValueError("parallax error truncation to high: getting negative distances")
real_ends = starts + diffs*real_distances
lzp = local_zero_point.reshape((-1, 1)) lzp = local_zero_point.reshape((-1, 1))
dist = np.array(self.domain[0].distances).reshape((-1, 1)) dist = np.array(self.domain[0].distances).reshape((-1, 1))
localized_pixel_starts = (starts-lzp)/dist + 0.5 localized_pixel_starts = (starts-lzp)/dist + 0.5
...@@ -179,7 +186,10 @@ class LOSResponse(LinearOperator): ...@@ -179,7 +186,10 @@ class LOSResponse(LinearOperator):
localized_pixel_ends, localized_pixel_ends,
self._local_shape, self._local_shape,
np.array(self.domain[0].distances), np.array(self.domain[0].distances),
difflen-sigmas_low, difflen, difflen+sigmas_up, 1./(1./difflen+truncation*sigmas),
difflen,
1./(1./difflen-truncation*sigmas),
sigmas,
_gaussian_sf) _gaussian_sf)
boxsz = 16 boxsz = 16
......
...@@ -38,9 +38,10 @@ class Consistency_Tests(unittest.TestCase): ...@@ -38,9 +38,10 @@ class Consistency_Tests(unittest.TestCase):
def testLOSResponse(self, sp, dtype): def testLOSResponse(self, sp, dtype):
starts = np.random.randn(len(sp.shape), 10) starts = np.random.randn(len(sp.shape), 10)
ends = np.random.randn(len(sp.shape), 10) ends = np.random.randn(len(sp.shape), 10)
sigma_low = 1e-4*np.random.randn(10) sigmas = 1e-4*np.random.randn(10)**2
sigma_ups = 1e-5*np.random.randn(10) while np.any(1./(1./np.linalg.norm(ends-starts,axis=0)-2.*sigmas)<0):
op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups) sigmas /= 2.
op = ift.LOSResponse(sp, starts, ends, sigmas, 2.)
ift.extra.consistency_check(op, dtype, dtype) ift.extra.consistency_check(op, dtype, dtype)
@expand(product(_h_spaces + _p_spaces + _pow_spaces, @expand(product(_h_spaces + _p_spaces + _pow_spaces,
......
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