Commit d90436e4 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_5' into docstrings_torsten

parents 0df19473 928858f4
...@@ -10,6 +10,7 @@ setup.cfg ...@@ -10,6 +10,7 @@ setup.cfg
.document .document
.svn/ .svn/
*.csv *.csv
.pytest_cache/
# from https://github.com/github/gitignore/blob/master/Python.gitignore # from https://github.com/github/gitignore/blob/master/Python.gitignore
......
...@@ -34,18 +34,22 @@ build_docker_from_cache: ...@@ -34,18 +34,22 @@ build_docker_from_cache:
- docker build -t $CONTAINER_TEST_IMAGE . - docker build -t $CONTAINER_TEST_IMAGE .
- docker push $CONTAINER_TEST_IMAGE - docker push $CONTAINER_TEST_IMAGE
test: test_serial:
stage: test stage: test
variables:
OMPI_MCA_btl_vader_single_copy_mechanism: none
script: script:
- mpiexec -n 2 --bind-to none pytest-3 -q test
- pytest-3 -q --cov=nifty5 test - pytest-3 -q --cov=nifty5 test
- > - >
python3 -m coverage report --omit "*plotting*,*distributed_do*" python3 -m coverage report --omit "*plotting*,*distributed_do*"
- > - >
python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }' 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: pages:
stage: release stage: release
before_script: before_script:
......
...@@ -7,9 +7,9 @@ RUN apt-get update && apt-get install -y \ ...@@ -7,9 +7,9 @@ RUN apt-get update && apt-get install -y \
libfftw3-dev \ libfftw3-dev \
python3 python3-pip python3-dev python3-future python3-scipy cython3 \ python3 python3-pip python3-dev python3-future python3-scipy cython3 \
# Documentation build dependencies # Documentation build dependencies
python3-sphinx python3-sphinx-rtd-theme python3-numpydoc \ python3-sphinx python3-sphinx-rtd-theme \
# Testing dependencies # Testing dependencies
python3-coverage python3-parameterized python3-pytest python3-pytest-cov \ python3-coverage python3-pytest python3-pytest-cov \
# Optional NIFTy dependencies # Optional NIFTy dependencies
openmpi-bin libopenmpi-dev python3-mpi4py \ openmpi-bin libopenmpi-dev python3-mpi4py \
# Packages needed for NIFTy # Packages needed for NIFTy
......
rm -rf docs/build docs/source/mod
sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty5 sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty5
sphinx-build -b html docs/source/ docs/build/ sphinx-build -b html docs/source/ docs/build/
.. currentmodule:: nifty5
============= =============
Code Overview Code Overview
...@@ -37,9 +36,12 @@ Domains ...@@ -37,9 +36,12 @@ Domains
Abstract base class Abstract base class
------------------- -------------------
.. currentmodule:: nifty5.domains.domain
One of the fundamental building blocks of the NIFTy5 framework is the *domain*. One of the fundamental building blocks of the NIFTy5 framework is the *domain*.
Its required capabilities are expressed by the abstract :class:`Domain` class. Its required capabilities are expressed by the abstract :py:class:`Domain` class.
A domain must be able to answer the following queries: A domain must be able to answer the following queries:
m
- its total number of data entries (pixels), which is accessible via the - its total number of data entries (pixels), which is accessible via the
:attr:`~Domain.size` property :attr:`~Domain.size` property
...@@ -51,6 +53,8 @@ A domain must be able to answer the following queries: ...@@ -51,6 +53,8 @@ A domain must be able to answer the following queries:
Unstructured domains Unstructured domains
-------------------- --------------------
.. currentmodule:: nifty5.domains.unstructured_domain
Domains can be either *structured* (i.e. there is geometrical information Domains can be either *structured* (i.e. there is geometrical information
associated with them, like position in space and volume factors), associated with them, like position in space and volume factors),
or *unstructured* (meaning that the data points have no associated manifold). or *unstructured* (meaning that the data points have no associated manifold).
...@@ -62,6 +66,8 @@ Unstructured domains can be described by instances of NIFTy's ...@@ -62,6 +66,8 @@ Unstructured domains can be described by instances of NIFTy's
Structured domains Structured domains
------------------ ------------------
.. currentmodule:: nifty5.domains.structured_domain
In contrast to unstructured domains, these domains have an assigned geometry. In contrast to unstructured domains, these domains have an assigned geometry.
NIFTy requires them to provide the volume elements of their grid cells. NIFTy requires them to provide the volume elements of their grid cells.
The additional methods are specified in the abstract class The additional methods are specified in the abstract class
...@@ -81,15 +87,17 @@ The additional methods are specified in the abstract class ...@@ -81,15 +87,17 @@ The additional methods are specified in the abstract class
NIFTy comes with several concrete subclasses of :class:`StructuredDomain`: NIFTy comes with several concrete subclasses of :class:`StructuredDomain`:
- :class:`RGSpace` represents a regular Cartesian grid with an arbitrary .. currentmodule:: nifty5.domains
- :class:`rg_space.RGSpace` represents a regular Cartesian grid with an arbitrary
number of dimensions, which is supposed to be periodic in each dimension. number of dimensions, which is supposed to be periodic in each dimension.
- :class:`HPSpace` and :class:`GLSpace` describe pixelisations of the - :class:`hp_space.HPSpace` and :class:`gl_space.GLSpace` describe pixelisations of the
2-sphere; their counterpart in harmonic space is :class:`LMSpace`, which 2-sphere; their counterpart in harmonic space is :class:`lm_space.LMSpace`, which
contains spherical harmonic coefficients. contains spherical harmonic coefficients.
- :class:`PowerSpace` is used to describe one-dimensional power spectra. - :class:`power_space.PowerSpace` is used to describe one-dimensional power spectra.
Among these, :class:`RGSpace` can be harmonic or not (depending on constructor arguments), :class:`GLSpace`, :class:`HPSpace`, and :class:`PowerSpace` are Among these, :class:`rg_space.RGSpace` can be harmonic or not (depending on constructor arguments), :class:`gl_space.GLSpace`, :class:`hp_space.HPSpace`, and :class:`power_space.PowerSpace` are
pure position domains (i.e. nonharmonic), and :class:`LMSpace` is always pure position domains (i.e. nonharmonic), and :class:`lm_space.LMSpace` is always
harmonic. harmonic.
...@@ -158,7 +166,7 @@ be extracted first, then changed, and a new field has to be created from the ...@@ -158,7 +166,7 @@ be extracted first, then changed, and a new field has to be created from the
result. result.
Fields defined on a MultiDomain Fields defined on a MultiDomain
------------------------------ -------------------------------
The :class:`MultiField` class can be seen as a dictionary of individual The :class:`MultiField` class can be seen as a dictionary of individual
:class:`Field` s, each identified by a name, which is defined on a :class:`Field` s, each identified by a name, which is defined on a
...@@ -300,7 +308,7 @@ As an example one may consider the following combination of ``x``, which is an o ...@@ -300,7 +308,7 @@ As an example one may consider the following combination of ``x``, which is an o
Basic operators Basic operators
------------ ---------------
# FIXME All this is outdated! # FIXME All this is outdated!
Basic operator classes provided by NIFTy are Basic operator classes provided by NIFTy are
......
import nifty5 import nifty5
import sphinx_rtd_theme
napoleon_google_docstring = False
napoleon_numpy_docstring = True
napoleon_use_ivar = True
napoleon_use_param = False
napoleon_use_keyword = False
autodoc_member_order = 'groupwise'
numpydoc_show_inherited_class_members = False
numpydoc_class_members_toctree = False
extensions = [ extensions = [
'sphinx.ext.autodoc', 'numpydoc', 'sphinx.ext.autosummary', 'sphinx.ext.napoleon', # Support for NumPy and Google style docstrings
'sphinx.ext.napoleon', 'sphinx.ext.imgmath', 'sphinx.ext.viewcode' 'sphinx.ext.imgmath', # Render math as images
'sphinx.ext.viewcode' # Add links to highlighted source code
] ]
templates_path = ['_templates']
source_suffix = '.rst'
master_doc = 'index' master_doc = 'index'
napoleon_google_docstring = False
napoleon_numpy_docstring = True
napoleon_use_ivar = True
project = u'NIFTy5' project = u'NIFTy5'
copyright = u'2013-2019, Max-Planck-Society' copyright = u'2013-2019, Max-Planck-Society'
author = u'Martin Reinecke' author = u'Martin Reinecke'
...@@ -29,19 +21,6 @@ version = release[:-2] ...@@ -29,19 +21,6 @@ version = release[:-2]
language = None language = None
exclude_patterns = [] exclude_patterns = []
add_module_names = False add_module_names = False
pygments_style = 'sphinx'
todo_include_todos = True
html_theme = "sphinx_rtd_theme" html_theme = "sphinx_rtd_theme"
html_theme_options = {
'collapse_navigation': False,
'display_version': False,
}
html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
html_logo = 'nifty_logo_black.png' html_logo = 'nifty_logo_black.png'
html_static_path = []
html_last_updated_fmt = '%b %d, %Y'
html_domain_indices = False
html_use_index = False
html_show_sourcelink = False
htmlhelp_basename = 'NIFTydoc'
...@@ -104,6 +104,7 @@ with :math:`{R}` the measurement response, which maps the continous signal field ...@@ -104,6 +104,7 @@ with :math:`{R}` the measurement response, which maps the continous signal field
This is called a free theory, as the information Hamiltonian This is called a free theory, as the information Hamiltonian
associate professor associate professor
.. math:: .. math::
\mathcal{H}(d,s)= -\log \mathcal{P}(d,s)= \frac{1}{2} s^\dagger S^{-1} s + \frac{1}{2} (d-R\,s)^\dagger N^{-1} (d-R\,s) + \mathrm{const} \mathcal{H}(d,s)= -\log \mathcal{P}(d,s)= \frac{1}{2} s^\dagger S^{-1} s + \frac{1}{2} (d-R\,s)^\dagger N^{-1} (d-R\,s) + \mathrm{const}
...@@ -179,23 +180,22 @@ NIFTy takes advantage of this formulation in several ways: ...@@ -179,23 +180,22 @@ NIFTy takes advantage of this formulation in several ways:
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. 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.
+-------------------------------------------------+ +----------------------------------------------------+
| .. image:: images/getting_started_3_setup.png | | **Output of tomography demo getting_started_3.py** |
| :width: 30 % | +----------------------------------------------------+
+-------------------------------------------------+ | .. image:: images/getting_started_3_setup.png |
| .. image:: images/getting_started_3_results.png | | |
| :width: 30 % | +----------------------------------------------------+
+-------------------------------------------------+ | Non-Gaussian signal field, |
| Output of tomography demo getting_started_3.py. | | data backprojected into the image domain, power |
| **Top row:** Non-Gaussian signal field, | | spectrum of underlying Gausssian process. |
| data backprojected into the image domain, power | +----------------------------------------------------+
| spectrum of underlying Gausssian process. | | .. image:: images/getting_started_3_results.png |
| **Bottom row:** Posterior mean field signal | | |
| reconstruction, its uncertainty, and the power | +----------------------------------------------------+
| spectrum of the process for different posterior | | Posterior mean field signal |
| samples in comparison to the correct one (thick | | reconstruction, its uncertainty, and the power |
| orange line). | | spectrum of the process for different posterior |
+-------------------------------------------------+ | samples in comparison to the correct one (thick |
| orange line). |
+----------------------------------------------------+
...@@ -17,6 +17,7 @@ ...@@ -17,6 +17,7 @@
import sys import sys
from functools import reduce from functools import reduce
import numpy as np import numpy as np
from mpi4py import MPI from mpi4py import MPI
......
...@@ -18,11 +18,10 @@ ...@@ -18,11 +18,10 @@
# Data object module for NIFTy that uses simple numpy ndarrays. # Data object module for NIFTy that uses simple numpy ndarrays.
import numpy as np import numpy as np
from numpy import empty, empty_like, exp, full, log from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log
from numpy import ndarray as data_object from numpy import ndarray as data_object
from numpy import ones, sqrt, tanh, vdot, zeros from numpy import ones, sign, sin, sinc, sinh, sqrt, tan, tanh, vdot, zeros
from numpy import sin, cos, tan, sinh, cosh, sinc
from numpy import absolute, sign, clip
from .random import Random from .random import Random
__all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
......
...@@ -111,23 +111,27 @@ def _SlopePowerSpectrum(logk_space, sm, sv, im, iv): ...@@ -111,23 +111,27 @@ def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv, def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True): keys=['tau', 'phi'], zero_mode=True):
''' ''' Operator for parametrizing smooth power spectra.
Computes a smooth power spectrum. Computes a smooth power spectrum.
Output is defined on a PowerSpace. Output is defined on a PowerSpace.
Parameters Parameters
---------- ----------
Npixdof : int
Npixdof : #pix in dof_space #pix in dof_space
ceps_a : float
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a, k0 = ceps_k0
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 ceps_k0 : float
a = ceps_a, k0 = ceps_k0 Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a, k0 = ceps_k0
sm : float
sm, sv : slope_mean = expected exponent of power law (e.g. -4), slope_mean = expected exponent of power law (e.g. -4)
slope_variance (default=1) sv : float
slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_variance of power_slope im : float
y-intercept_mean
iv : float
y-intercept_variance of power_slope
''' '''
from ..operators.exp_transform import ExpTransform from ..operators.exp_transform import ExpTransform
......
...@@ -28,14 +28,14 @@ from ..field import Field ...@@ -28,14 +28,14 @@ from ..field import Field
from ..operators.linear_operator import LinearOperator from ..operators.linear_operator import LinearOperator
def _gaussian_error_function(x): 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) inc = np.full(len(shp), 1, dtype=np.int64)
for i in range(-2, -len(shp)-1, -1): for i in range(-2, -len(shp)-1, -1):
inc[i] = inc[i+1]*shp[i+1] inc[i] = inc[i+1]*shp[i+1]
...@@ -59,7 +59,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf): ...@@ -59,7 +59,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
dmin += 1e-7 dmin += 1e-7
dmax -= 1e-7 dmax -= 1e-7
if dmin >= dmax: # no intersection 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 continue
# determine coordinates of first cell crossing # determine coordinates of first cell crossing
c_first = np.ceil(start[:, i]+dir*dmin) c_first = np.ceil(start[:, i]+dir*dmin)
...@@ -75,7 +75,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf): ...@@ -75,7 +75,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
tmp = np.arange(start=c_first[j], stop=dmax, tmp = np.arange(start=c_first[j], stop=dmax,
step=abs(1./dir[j])) step=abs(1./dir[j]))
cdist = np.append(cdist, tmp) 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) idx = np.argsort(cdist)
cdist = cdist[idx] cdist = cdist[idx]
add = add[idx] add = add[idx]
...@@ -85,21 +85,19 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf): ...@@ -85,21 +85,19 @@ 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 > mid) & (dist <= hi) mask = (dist > lo) & (dist <= hi)
wgt[mask] *= erf((dist[mask]-mid)/(hi-mid)) wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig)
mask = (dist <= mid) & (dist > lo)
wgt[mask] *= erf((dist[mask]-mid)/(mid-lo))
return wgt return wgt
...@@ -119,17 +117,32 @@ class LOSResponse(LinearOperator): ...@@ -119,17 +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)
For expert use. If unsure, leave blank. 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 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
...@@ -141,20 +154,15 @@ class LOSResponse(LinearOperator): ...@@ -141,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(
...@@ -164,7 +172,11 @@ class LOSResponse(LinearOperator): ...@@ -164,7 +172,11 @@ 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
...@@ -175,8 +187,11 @@ class LOSResponse(LinearOperator): ...@@ -175,8 +187,11 @@ 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),
_gaussian_error_function) difflen,
1./(1./difflen-truncation*sigmas),
sigmas,
_gaussian_sf)
boxsz = 16 boxsz = 16
nlos = len(w_i) nlos = len(w_i)
...@@ -202,7 +217,7 @@ class LOSResponse(LinearOperator): ...@@ -202,7 +217,7 @@ class LOSResponse(LinearOperator):
tmp += (fullidx[j]//boxsz)*fct tmp += (fullidx[j]//boxsz)*fct
fct *= self._local_shape[j] fct *= self._local_shape[j]
tmp += cnt/float(nlos) tmp += cnt/float(nlos)
tmp += iarr[ofs:ofs+nval]/float(nlos*npix) tmp += iarr[ofs:ofs+nval]/(float(nlos)*float(npix))
pri[ofs:ofs+nval] = tmp pri[ofs:ofs+nval] = tmp
ofs += nval ofs += nval
cnt += 1 cnt += 1
......
...@@ -110,11 +110,11 @@ class Linearization(object): ...@@ -110,11 +110,11 @@ class Linearization(object):
def __truediv__(self, other): def __truediv__(self, other):
if isinstance(other, Linearization): if isinstance(other, Linearization):
return self.__mul__(other.inverse()) return self.__mul__(other.one_over())
return self.__mul__(1./other) return self.__mul__(1./other)
def __rtruediv__(self, other): def __rtruediv__(self, other):