Commit 96bdce1d authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'NIFTy_5' into better_pytest

parents 14d98173 34af9f25
......@@ -34,18 +34,22 @@ build_docker_from_cache:
- docker build -t $CONTAINER_TEST_IMAGE .
- docker push $CONTAINER_TEST_IMAGE
test:
test_serial:
stage: test
variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none
script:
- mpiexec -n 2 --bind-to none pytest-3 -q test
- pytest-3 -q --cov=nifty5 test
- >
python3 -m coverage report --omit "*plotting*,*distributed_do*"
- >
python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
test_mpi:
stage: test
variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none
script:
- mpiexec -n 2 --bind-to none pytest-3 -q test
pages:
stage: release
before_script:
......
......@@ -84,7 +84,8 @@ MPI support is added via:
### Running the tests
In oder to run the tests one needs two additional packages:
To run the tests, additional packages are required:
sudo apt-get install python3-coverage python3-parameterized python3-pytest python3-pytest-cov
Afterwards the tests (including a coverage report) can be run using the
......
......@@ -53,8 +53,8 @@ if __name__ == '__main__':
A = ift.create_power_operator(harmonic_space, sqrtpspec)
# Set up a sky model and instrumental response
sky = ift.positive_tanh(HT(A))
# Set up a sky operator and instrumental response
sky = ift.sigmoid(HT(A))
GR = ift.GeometryRemover(position_space)
R = GR
......
......@@ -76,7 +76,7 @@ if __name__ == '__main__':
a = ift.PS_field(p_space, sqrtpspec)
A = pd(a)
# Define sky model
# Define sky operator
sky = ift.exp(HT(ift.makeOp(A)))
M = ift.DiagonalOperator(exposure)
......
......@@ -47,7 +47,7 @@ if __name__ == '__main__':
position_space = ift.RGSpace([128, 128])
# Set up an amplitude model for the field
# Set up an amplitude operator for the field
# The parameters mean:
# 64 spectral bins
#
......@@ -60,9 +60,9 @@ if __name__ == '__main__':
# 0.5 = low variance of power-law slope
# 0.4 = y-intercept mean
# 0.3 = relatively high y-intercept variance
A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
A = ift.AmplitudeOperator(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
# Build the model for a correlated signal
# Build the operator for a correlated signal
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = A.target[0]
......@@ -76,7 +76,7 @@ if __name__ == '__main__':
# correlated_field = ift.CorrelatedField(position_space, A)
# Apply a nonlinearity
signal = ift.positive_tanh(correlated_field)
signal = ift.sigmoid(correlated_field)
# Build the line-of-sight response and define signal response
LOS_starts, LOS_ends = random_los(100) if mode == 1 else radial_los(100)
......@@ -98,7 +98,7 @@ if __name__ == '__main__':
name='Newton', tol=1e-7, iteration_limit=35)
minimizer = ift.NewtonCG(ic_newton)
# Set up model likelihood and information Hamiltonian
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
H = ift.Hamiltonian(likelihood, ic_sampling)
......
......@@ -15,7 +15,7 @@ From such a perspective,
- IFT problems largely consist of the combination of several high dimensional
*minimization* problems.
- Within NIFTy, *models* are used to define the characteristic equations and
- Within NIFTy, *operators* are used to define the characteristic equations and
properties of the problems.
- The equations are built mostly from the application of *linear operators*,
but there may also be nonlinear functions involved.
......@@ -99,13 +99,13 @@ Combinations of domains
The fundamental classes described above are often sufficient to specify the
domain of a field. In some cases, however, it will be necessary to define the
field on a product of elementary domains instead of a single one.
More sophisticated models also require a set of several such fields.
More sophisticated operators also require a set of several such fields.
Some examples are:
- sky emission depending on location and energy. This could be represented by
a product of an :class:`HPSpace` (for location) with an :class:`RGSpace`
(for energy).
- a polarised field, which could be modeled as a product of any structured
- a polarized field, which could be modeled as a product of any structured
domain (representing location) with a four-element
:class:`UnstructuredDomain` holding Stokes I, Q, U and V components.
- a model for the sky emission, which holds both the current realization
......@@ -136,7 +136,7 @@ A :class:`Field` object consists of the following components:
- a data type (e.g. numpy.float64)
- an array containing the actual values
Usually, the array is stored in the for of a ``numpy.ndarray``, but for very
Usually, the array is stored in the form of a ``numpy.ndarray``, but for very
resource-intensive tasks NIFTy also provides an alternative storage method to
be used with distributed memory processing.
......@@ -269,59 +269,59 @@ The properties :attr:`~LinearOperator.adjoint` and
were the original operator's adjoint or inverse, respectively.
Models
======
Operators
=========
Model classes (represented by NIFTy5's abstract :class:`Model` class) are used to construct
Operator classes (represented by NIFTy5's abstract :class:`Operator` class) are used to construct
the equations of a specific inference problem.
Most models are defined via a position, which is a :class:`MultiField` object,
Most operators are defined via a position, which is a :class:`MultiField` object,
their value at this position, which is again a :class:`MultiField` object and a Jacobian derivative,
which is a :class:`LinearOperator` and is needed for the minimization procedure.
Using the existing basic model classes one can construct more complicated models, as
Using the existing basic operator classes one can construct more complicated operators, as
NIFTy allows for easy and self-consinstent combination via point-wise multiplication,
addition and subtraction. The model resulting from these operations then automatically
addition and subtraction. The operator resulting from these operations then automatically
contains the correct Jacobians, positions and values.
Notably, :class:`Constant` and :class:`Variable` allow for an easy way to turn
inference of specific quantities on and off.
The basic model classes also allow for more complex operations on models such as
The basic operator classes also allow for more complex operations on operators such as
the application of :class:`LinearOperators` or local non-linearities.
As an example one may consider the following combination of ``x``, which is a model of type
:class:`Variable` and ``y``, which is a model of type :class:`Constant`::
As an example one may consider the following combination of ``x``, which is an operator of type
:class:`Variable` and ``y``, which is an operator of type :class:`Constant`::
z = x*x + y
``z`` will then be a model with the following properties::
``z`` will then be an operator with the following properties::
z.value = x.value*x.value + y.value
z.position = Union(x.position, y.position)
z.jacobian = 2*makeOp(x.value)
Basic models
Basic operators
------------
# FIXME All this is outdated!
Basic model classes provided by NIFTy are
Basic operator classes provided by NIFTy are
- :class:`Constant` contains a constant value and has a zero valued Jacobian.
Like other models, it has a position, but its value does not depend on it.
Like other operators, it has a position, but its value does not depend on it.
- :class:`Variable` returns the position as its value, its derivative is one.
- :class:`LinearModel` applies a :class:`LinearOperator` on the model.
- :class:`LocalModel` applies a non-linearity locally on the model.
- :class:`MultiModel` combines various models into one. In this case the position,
value and Jacobian are combined into corresponding :class:`MultiFields` and operators.
Advanced models
---------------
Advanced operators
------------------
NIFTy also provides a library of more sophisticated models which are used for more
NIFTy also provides a library of more sophisticated operators which are used for more
specific inference problems. Currently these are:
- :class:`AmplitudeModel`, which returns a smooth power spectrum.
- :class:`PointModel`, which models point sources which follow a inverse gamma distribution.
- :class:`SmoothSkyModel`, which models a diffuse lognormal field. It takes an amplitude model
- :class:`AmplitudeOperator`, which returns a smooth power spectrum.
- :class:`InverseGammaOperator`, which models point sources which follow a inverse gamma distribution.
- :class:`CorrelatedField`, which models a diffuse log-normal field. It takes an amplitude operator
to specify the correlation structure of the field.
......
......@@ -174,7 +174,7 @@ NIFTy takes advantage of this formulation in several ways:
1) All prior degrees of freedom have unit covariance which improves the condition number of operators which need to be inverted.
2) The amplitude operator can be regarded as part of the response, :math:`{R'=R\,A}`. In general, more sophisticated responses can be constructed out of the composition of simpler operators.
3) The response can be non-linear, e.g. :math:`{R'(s)=R \exp(A\,\xi)}`, see demos/getting_started_2.py.
4) The amplitude operator can be made dependent on unknowns as well, e.g. :math:`A=A(\tau)= F\, \widehat{e^\tau}` represents an amplitude model with a positive definite, unknown spectrum defined in Fourier domain. The amplitude field :math:`{\tau}` would get its own amplitude model, with a cepstrum (spectrum of a log spectrum) defined in quefrency space (harmonic space of a logarithmically binned harmonic space) to regularize its degrees of freedom by imposing some (user-defined level of) spectral smoothness.
4) The amplitude operator can be made dependent on unknowns as well, e.g. :math:`A=A(\tau)= F\, \widehat{e^\tau}` represents an amplitude operator with a positive definite, unknown spectrum defined in Fourier domain. The amplitude field :math:`{\tau}` would get its own amplitude operator, with a cepstrum (spectrum of a log spectrum) defined in quefrency space (harmonic space of a logarithmically binned harmonic space) to regularize its degrees of freedom by imposing some (user-defined level of) spectral smoothness.
5) NIFTy can calculate the gradient of the information Hamiltonian and the Fischer information metric with respect to all unknown parameters, here :math:`{\xi}` and :math:`{\tau}`, by automatic differentiation. The gradients are used for MAP and HMCF estimates, and the Fischer matrix is required in addition to the gradient by Metric Gaussian Variational Inference (MGVI), which is available in NIFTy as well. MGVI is an implicit operator extension of Automatic Differentiation Variational Inference (ADVI).
The reconstruction of a non-Gaussian signal with unknown covarinance from a non-trivial (tomographic) response is demonstrated in demos/getting_started_3.py. Here, the uncertainty of the field and the power spectrum of its generating process are probed via posterior samples provided by the MGVI algorithm.
......
......@@ -73,8 +73,8 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
from .library.amplitude_model import AmplitudeModel
from .library.inverse_gamma_model import InverseGammaModel
from .library.amplitude_operator import AmplitudeOperator
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.wiener_filter_curvature import WienerFilterCurvature
......
......@@ -31,7 +31,8 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
"redistribute", "default_distaxis", "is_numpy", "absmax", "norm",
"lock", "locked", "uniform_full", "transpose", "to_global_data_rw",
"ensure_not_distributed", "ensure_default_distributed",
"clipped_exp"]
"tanh", "conjugate", "sin", "cos", "tan",
"sinh", "cosh", "sinc", "absolute", "sign", "clip"]
_comm = MPI.COMM_WORLD
ntask = _comm.Get_size()
......@@ -211,6 +212,9 @@ class data_object(object):
else:
return data_object(self._shape, tval, self._distaxis)
def clip(self, min=None, max=None):
return data_object(self._shape, np.clip(self._data, min, max))
def __neg__(self):
return data_object(self._shape, -self._data, self._distaxis)
......@@ -292,7 +296,8 @@ def _math_helper(x, function, out):
_current_module = sys.modules[__name__]
for f in ["sqrt", "exp", "log", "tanh", "conjugate"]:
for f in ["sqrt", "exp", "log", "tanh", "conjugate", "sin", "cos", "tan",
"sinh", "cosh", "sinc", "absolute", "sign"]:
def func(f):
def func2(x, out=None):
return _math_helper(x, f, out)
......@@ -300,8 +305,8 @@ for f in ["sqrt", "exp", "log", "tanh", "conjugate"]:
setattr(_current_module, f, func(f))
def clipped_exp(a):
return data_object(x.shape, np.exp(np.clip(x.data, -300, 300), x.distaxis))
def clip(x, a_min=None, a_max=None):
return data_object(x.shape, np.clip(x.data, a_min, a_max), x.distaxis)
def from_object(object, dtype, copy, set_locked):
......
......@@ -21,7 +21,8 @@ import numpy as np
from numpy import empty, empty_like, exp, full, log
from numpy import ndarray as data_object
from numpy import ones, sqrt, tanh, vdot, zeros
from numpy import sin, cos, tan, sinh, cosh, sinc
from numpy import absolute, sign, clip
from .random import Random
__all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
......@@ -33,7 +34,8 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
"redistribute", "default_distaxis", "is_numpy", "absmax", "norm",
"lock", "locked", "uniform_full", "to_global_data_rw",
"ensure_not_distributed", "ensure_default_distributed",
"clipped_exp"]
"clip", "sin", "cos", "tan", "sinh",
"cosh", "absolute", "sign", "sinc"]
ntask = 1
rank = 0
......@@ -149,7 +151,3 @@ def absmax(arr):
def norm(arr, ord=2):
return np.linalg.norm(arr.reshape(-1), ord=ord)
def clipped_exp(arr):
return np.exp(np.clip(arr, -300, 300))
......@@ -628,11 +628,14 @@ class Field(object):
def flexible_addsub(self, other, neg):
return self-other if neg else self+other
def positive_tanh(self):
def sigmoid(self):
return 0.5*(1.+self.tanh())
def clipped_exp(self):
return Field(self._domain, dobj.clipped_exp(self._val))
def clip(self, min=None, max=None):
return Field(self._domain, dobj.clip(self._val, min, max))
def one_over(self):
return 1/self
def _binary_op(self, other, op):
# if other is a field, make sure that the domains match
......@@ -669,7 +672,9 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
return func2
setattr(Field, op, func(op))
for f in ["sqrt", "exp", "log", "tanh"]:
for f in ["sqrt", "exp", "log", "tanh",
"sin", "cos", "tan", "cosh", "sinh",
"absolute", "sinc", "sign"]:
def func(f):
def func2(self):
return Field(self._domain, getattr(dobj, f)(self.val))
......
......@@ -75,14 +75,14 @@ def make_adjust_variances(a,
def do_adjust_variances(position,
amplitude_model,
amplitude_operator,
minimizer,
xi_key='xi',
samples=[]):
h_space = position[xi_key].domain[0]
pd = PowerDistributor(h_space, amplitude_model.target[0])
a = pd(amplitude_model)
pd = PowerDistributor(h_space, amplitude_operator.target[0])
a = pd(amplitude_operator)
xi = ducktape(None, position.domain, xi_key)
ham = make_adjust_variances(a, xi, position, samples=samples)
......
......@@ -62,7 +62,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
return Field.from_global_data(domain, cepstrum_field)
def CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
def _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
'''
Parameters
----------
......@@ -88,7 +88,7 @@ def CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
return res
def SlopeModel(logk_space, sm, sv, im, iv):
def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
'''
Parameters
----------
......@@ -109,8 +109,8 @@ def SlopeModel(logk_space, sm, sv, im, iv):
return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
def AmplitudeModel(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True):
def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True):
'''
Computes a smooth power spectrum.
Output is defined on a PowerSpace.
......@@ -137,9 +137,9 @@ def AmplitudeModel(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
et = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = et.domain[0]
smooth = CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
smooth = _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
smooth = smooth.ducktape(keys[0])
linear = SlopeModel(logk_space, sm, sv, im, iv)
linear = _SlopePowerSpectrum(logk_space, sm, sv, im, iv)
linear = linear.ducktape(keys[1])
fac = ScalingOperator(0.5, smooth.target)
......
......@@ -16,7 +16,6 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..domain_tuple import DomainTuple
from ..multi_domain import MultiDomain
from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
......@@ -24,7 +23,7 @@ from ..operators.simple_linear_operators import ducktape
from ..operators.scaling_operator import ScalingOperator
def CorrelatedField(s_space, amplitude_model, name='xi'):
def CorrelatedField(s_space, amplitude_operator, name='xi'):
'''
Function for construction of correlated fields
......@@ -32,23 +31,26 @@ def CorrelatedField(s_space, amplitude_model, name='xi'):
----------
s_space : Domain
Field domain
amplitude_model: Operator
model for correlation structure
amplitude_operator: Operator
operator for correlation structure
name : string
MultiField component name
'''
h_space = s_space.get_default_codomain()
ht = HarmonicTransformOperator(h_space, s_space)
p_space = amplitude_model.target[0]
p_space = amplitude_operator.target[0]
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_model)
A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol
vol = ScalingOperator(vol**(-0.5), h_space)
return ht(vol(A)*ducktape(h_space, None, name))
def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
amplitude_model_energy, name="xi"):
def MfCorrelatedField(s_space_spatial,
s_space_energy,
amplitude_operator_spatial,
amplitude_operator_energy,
name="xi"):
'''
Method for construction of correlated multi-frequency fields
'''
......@@ -59,8 +61,8 @@ def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
ht2 = HarmonicTransformOperator(ht1.target, space=1)
ht = ht2(ht1)
p_space_spatial = amplitude_model_spatial.target[0]
p_space_energy = amplitude_model_energy.target[0]
p_space_spatial = amplitude_operator_spatial.target[0]
p_space_energy = amplitude_operator_energy.target[0]
pd_spatial = PowerDistributor(h_space, p_space_spatial, 0)
pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1)
......@@ -69,8 +71,8 @@ def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
dom_distr_spatial = ContractionOperator(pd.domain, 1).adjoint
dom_distr_energy = ContractionOperator(pd.domain, 0).adjoint
a_spatial = dom_distr_spatial(amplitude_model_spatial)
a_energy = dom_distr_energy(amplitude_model_energy)
a_spatial = dom_distr_spatial(amplitude_operator_spatial)
a_energy = dom_distr_energy(amplitude_operator_energy)
a = a_spatial*a_energy
A = pd(a)
return ht(A*ducktape(h_space, None, name))
......@@ -25,9 +25,9 @@ from ..operators.operator import Operator
from ..sugar import makeOp
class InverseGammaModel(Operator):
class InverseGammaOperator(Operator):
def __init__(self, domain, alpha, q, delta=0.001):
"""Model which transforms a Gaussian into an inverse gamma distribution.
"""Operator which transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
......
......@@ -28,14 +28,14 @@ from ..field import Field
from ..operators.linear_operator import LinearOperator
def _gaussian_error_function(x):
return 0.5*erfc(x*np.sqrt(2.))
def _gaussian_sf(x):
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]
nlos = start.shape[1]
inc = np.full(len(shp), 1)
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]
......@@ -59,7 +59,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
dmin += 1e-7
dmax -= 1e-7
if dmin >= dmax: # no intersection
out[i] = (np.full(0, 0), np.full(0, 0.))
out[i] = (np.full(0, 0, dtype=np.int64), np.full(0, 0.))
continue
# determine coordinates of first cell crossing
c_first = np.ceil(start[:, i]+dir*dmin)
......@@ -75,7 +75,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
tmp = np.arange(start=c_first[j], stop=dmax,
step=abs(1./dir[j]))
cdist = np.append(cdist, tmp)
add = np.append(add, np.full(len(tmp), step))
add = np.append(add, np.full(len(tmp), step, dtype=np.int64))
idx = np.argsort(cdist)
cdist = cdist[idx]
add = add[idx]
......@@ -85,21 +85,19 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
cdist *= corfac
wgt = np.diff(cdist)
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.cumsum(add)
out[i] = (add, wgt)
return out
def apply_erf(wgt, dist, lo, mid, hi, erf):
def apply_erf(wgt, dist, lo, mid, hi, sig, erf):
wgt = wgt.copy()
mask = dist > hi
wgt[mask] = 0.
mask = (dist > mid) & (dist <= hi)
wgt[mask] *= erf((dist[mask]-mid)/(hi-mid))
mask = (dist <= mid) & (dist > lo)
wgt[mask] *= erf((dist[mask]-mid)/(mid-lo))
mask = (dist > lo) & (dist <= hi)
wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig)
return wgt
......@@ -119,17 +117,32 @@ class LOSResponse(LinearOperator):
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.
sigmas_low, sigmas_up : numpy.ndarray(float) (optional)
For expert use. If unsure, leave blank.
sigmas: numpy.ndarray(float) (optional)
If this is not None, the inverse of the lengths of the LOSs are assumed
to be Gaussian distributed with these sigmas. The start point will
remain the same, but the endpoint is assumed to be unknown.
This is a typical statistical model for astrophysical parallaxes.
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.
default: 3.
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 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._capability = self.TIMES | self.ADJOINT_TIMES
......@@ -141,20 +154,15 @@ class LOSResponse(LinearOperator):
starts = np.array(starts)
nlos = starts.shape[1]
ends = np.array(ends)
if sigmas_low is None:
sigmas_low = np.zeros(nlos, dtype=np.float32)
if sigmas_up is None:
sigmas_up = np.zeros(nlos, dtype=np.float32)
sigmas_low = np.array(sigmas_low)
sigmas_up = np.array(sigmas_up)
if sigmas is None:
sigmas = np.zeros(nlos, dtype=np.float32)
sigmas = np.array(sigmas)
if starts.shape[0] != ndim:
raise TypeError("dimension mismatch")
if nlos != sigmas_low.shape[0]:
if nlos != sigmas.shape[0]:
raise TypeError("dimension mismatch")
if starts.shape != ends.shape:
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)
local_zero_point = (np.array(
......@@ -164,7 +172,11 @@ class LOSResponse(LinearOperator):
diffs = ends-starts
difflen = np.linalg.norm(diffs, axis=0)
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))
dist = np.array(self.domain[0].distances).reshape((-1, 1))
localized_pixel_starts = (starts-lzp)/dist + 0.5
......@@ -175,8 +187,11 @@ class LOSResponse(LinearOperator):
localized_pixel_ends,
self._local_shape,
np.array(self.domain[0].distances),