Commit 6e22052f authored by Lukas Platz's avatar Lukas Platz

Merge branch 'NIFTy_5' into mf_new

parents bacaa93f ad4211af
Pipeline #53279 passed with stages
in 6 minutes and 54 seconds
......@@ -85,7 +85,7 @@ Support for spherical harmonic transforms is added via:
Support for the radio interferometry gridder is added via:
pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
MPI support is added via:
......
......@@ -9,7 +9,7 @@ np.random.seed(40)
N0s, a0s, b0s, c0s = [], [], [], []
for ii in range(10, 23):
for ii in range(10, 26):
nu = 1024
nv = 1024
N = int(2**ii)
......@@ -27,17 +27,15 @@ for ii in range(10, 23):
img = ift.from_global_data(uvspace, img)
t0 = time()
GM = ift.GridderMaker(uvspace, eps=1e-7)
idx = GM.getReordering(uv)
uv = uv[idx]
vis = vis[idx]
GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
vis = ift.from_global_data(visspace, vis)
op = GM.getFull(uv).adjoint
op = GM.getFull().adjoint
t1 = time()
op(img).to_global_data()
t2 = time()
op.adjoint(vis).to_global_data()
t3 = time()
print(t2-t1, t3-t2)
N0s.append(N)
a0s.append(t1 - t0)
b0s.append(t2 - t1)
......
......@@ -109,7 +109,8 @@ if __name__ == '__main__':
minimizer = ift.NewtonCG(ic_newton)
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
likelihood = ift.GaussianEnergy(mean=data,
inverse_covariance=N.inverse)(signal_response)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
initial_mean = ift.MultiField.full(H.domain, 0.)
......
......@@ -19,9 +19,9 @@ Support for spherical harmonic transforms is added via::
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
Support for the radio interferometry gridder is added via:
Support for the radio interferometry gridder is added via::
pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
MPI support is added via::
......
......@@ -50,6 +50,8 @@ class LogRGSpace(StructuredDomain):
self._bindistances = tuple(bindistances)
self._t_0 = tuple(t_0)
if min(self._bindistances) <= 0:
raise ValueError('Non-positive bindistances encountered')
self._dim = int(reduce(lambda x, y: x*y, self._shape))
self._dvol = float(reduce(lambda x, y: x*y, self._bindistances))
......@@ -80,8 +82,8 @@ class LogRGSpace(StructuredDomain):
return np.array(self._t_0)
def __repr__(self):
return ("LogRGSpace(shape={}, harmonic={})".format(
self.shape, self.harmonic))
return ("LogRGSpace(shape={}, bindistances={}, t_0={}, harmonic={})".format(
self.shape, self.bindistances, self.t_0, self.harmonic))
def get_default_codomain(self):
"""Returns a :class:`LogRGSpace` object representing the (position or
......
......@@ -165,6 +165,8 @@ class PowerSpace(StructuredDomain):
if binbounds is not None:
binbounds = tuple(binbounds)
if min(binbounds) < 0:
raise ValueError('Negative binbounds encountered')
key = (harmonic_partner, binbounds)
if self._powerIndexCache.get(key) is None:
......
......@@ -54,6 +54,8 @@ class RGSpace(StructuredDomain):
if np.isscalar(shape):
shape = (shape,)
self._shape = tuple(int(i) for i in shape)
if min(self._shape) < 0:
raise ValueError('Negative number of pixels encountered')
if distances is None:
if self.harmonic:
......@@ -66,6 +68,8 @@ class RGSpace(StructuredDomain):
temp = np.empty(len(self.shape), dtype=np.float64)
temp[:] = distances
self._distances = tuple(temp)
if min(self._distances) <= 0:
raise ValueError('Non-positive distances encountered')
self._dvol = float(reduce(lambda x, y: x*y, self._distances))
self._size = int(reduce(lambda x, y: x*y, self._shape))
......
......@@ -15,8 +15,6 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from .utilities import iscomplextype
import numpy as np
import pypocketfft
_nthreads = 1
......@@ -31,81 +29,14 @@ def set_nthreads(nthr):
_nthreads = nthr
# FIXME this should not be necessary ... no one should call a complex FFT
# with a float array.
def _make_complex(a):
if a.dtype in (np.complex64, np.complex128):
return a
if a.dtype == np.float64:
return a.astype(np.complex128)
if a.dtype == np.float32:
return a.astype(np.complex64)
raise NotImplementedError
def fftn(a, axes=None):
return pypocketfft.fftn(_make_complex(a), axes=axes, nthreads=_nthreads)
def rfftn(a, axes=None):
return pypocketfft.rfftn(a, axes=axes, nthreads=_nthreads)
return pypocketfft.c2c(a, axes=axes, nthreads=_nthreads)
def ifftn(a, axes=None):
# FIXME this is a temporary fix and can be done more elegantly
if axes is None:
fct = 1./a.size
else:
fct = 1./np.prod(np.take(a.shape, axes))
return pypocketfft.ifftn(_make_complex(a), axes=axes, fct=fct,
nthreads=_nthreads)
return pypocketfft.c2c(a, axes=axes, inorm=2, forward=False,
nthreads=_nthreads)
def hartley(a, axes=None):
return pypocketfft.hartley2(a, axes=axes, nthreads=_nthreads)
# Do a real-to-complex forward FFT and return the _full_ output array
def my_fftn_r2c(a, axes=None):
# Check if the axes provided are valid given the shape
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes do not match array shape")
if iscomplextype(a.dtype):
raise TypeError("Transform requires real-valued input arrays.")
tmp = rfftn(a, axes=axes)
def _fill_complex_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
res[tuple(slice1)] = tmp
def _fill_upper_half_complex(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
# np.conjugate(tmp[slice2], out=res[slice1])
res[tuple(slice1)] = np.conjugate(tmp[tuple(slice2)])
for i, ax in enumerate(axes[:-1]):
dim1 = tuple([slice(None)]*ax + [slice(0, 1)])
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half_complex(tmp[dim1], res[dim1], axes2)
_fill_upper_half_complex(tmp, res, axes)
return res
return _fill_complex_array(tmp, np.empty_like(a, dtype=tmp.dtype), axes)
def my_fftn(a, axes=None):
return fftn(a, axes=axes)
return pypocketfft.genuine_hartley(a, axes=axes, nthreads=_nthreads)
......@@ -15,108 +15,96 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..fft import hartley
from ..operators.linear_operator import LinearOperator
from ..sugar import from_global_data, makeDomain
import numpy as np
class GridderMaker(object):
def __init__(self, domain, eps=2e-13):
from nifty_gridder import get_w
domain = makeDomain(domain)
if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or
not len(domain.shape) == 2):
raise ValueError("need domain with exactly one 2D RGSpace")
nu, nv = domain.shape
if nu % 2 != 0 or nv % 2 != 0:
raise ValueError("dimensions must be even")
nu2, nv2 = 2*nu, 2*nv
w = get_w(eps)
nsafe = (w+1)//2
nu2 = max([nu2, 2*nsafe])
nv2 = max([nv2, 2*nsafe])
oversampled_domain = RGSpace(
[nu2, nv2], distances=[1, 1], harmonic=False)
self._eps = eps
self._rest = _RestOperator(domain, oversampled_domain, eps)
def getReordering(self, uv):
from nifty_gridder import peanoindex
nu2, nv2 = self._rest._domain.shape
return peanoindex(uv, nu2, nv2)
def getGridder(self, uv):
return RadioGridder(self._rest.domain, self._eps, uv)
def __init__(self, dirty_domain, uv, eps=2e-13):
import nifty_gridder
dirty_domain = makeDomain(dirty_domain)
if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace)
or not len(dirty_domain.shape) == 2):
raise ValueError("need dirty_domain with exactly one 2D RGSpace")
if uv.ndim != 2:
raise ValueError("uv must be a 2D array")
if uv.shape[1] != 2:
raise ValueError("second dimension of uv must have length 2")
dstx, dsty = dirty_domain[0].distances
# wasteful hack to adjust to shape required by nifty_gridder
uvw = np.empty((uv.shape[0], 3), dtype=np.float64)
uvw[:, 0:2] = uv
uvw[:, 2] = 0.
# Scale uv such that 0<uv<=1 which is assumed by nifty_gridder
uvw[:, 0] = uvw[:, 0]*dstx
uvw[:, 1] = uvw[:, 1]*dsty
speedOfLight = 299792458.
bl = nifty_gridder.Baselines(uvw, np.array([speedOfLight]))
nxdirty, nydirty = dirty_domain.shape
gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.)
nu, nv = gconf.Nu(), gconf.Nv()
self._idx = nifty_gridder.getIndices(
bl, gconf, np.zeros((uv.shape[0], 1), dtype=np.bool))
self._bl = bl
du, dv = 1./(nu*dstx), 1./(nv*dsty)
grid_domain = RGSpace([nu, nv], distances=[du, dv], harmonic=True)
self._rest = _RestOperator(dirty_domain, grid_domain, gconf)
self._gridder = RadioGridder(grid_domain, bl, gconf, self._idx)
def getGridder(self):
return self._gridder
def getRest(self):
return self._rest
def getFull(self, uv):
return self.getRest() @ self.getGridder(uv)
def getFull(self):
return self.getRest() @ self._gridder
def ms2vis(self, x):
return self._bl.ms2vis(x, self._idx)
class _RestOperator(LinearOperator):
def __init__(self, domain, oversampled_domain, eps):
from nifty_gridder import correction_factors
self._domain = makeDomain(oversampled_domain)
self._target = domain
nu, nv = domain.shape
nu2, nv2 = oversampled_domain.shape
fu = correction_factors(nu2, nu//2+1, eps)
fv = correction_factors(nv2, nv//2+1, eps)
# compute deconvolution operator
rng = np.arange(nu)
k = np.minimum(rng, nu-rng)
self._deconv_u = np.roll(fu[k], -nu//2).reshape((-1, 1))
rng = np.arange(nv)
k = np.minimum(rng, nv-rng)
self._deconv_v = np.roll(fv[k], -nv//2).reshape((1, -1))
class _RestOperator(LinearOperator):
def __init__(self, dirty_domain, grid_domain, gconf):
self._domain = makeDomain(grid_domain)
self._target = makeDomain(dirty_domain)
self._gconf = gconf
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
nu, nv = self._target.shape
res = x.to_global_data()
if mode == self.TIMES:
res = hartley(res)
res = np.roll(res, (nu//2, nv//2), axis=(0, 1))
res = res[:nu, :nv]
res *= self._deconv_u
res *= self._deconv_v
res = self._gconf.grid2dirty(res)
else:
res = res*self._deconv_u
res *= self._deconv_v
nu2, nv2 = self._domain.shape
res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant',
constant_values=0)
res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1))
res = hartley(res)
res = self._gconf.dirty2grid(res)
return from_global_data(self._tgt(mode), res)
class RadioGridder(LinearOperator):
def __init__(self, target, eps, uv):
def __init__(self, grid_domain, bl, gconf, idx):
self._domain = DomainTuple.make(
UnstructuredDomain((uv.shape[0],)))
self._target = DomainTuple.make(target)
UnstructuredDomain((bl.Nrows())))
self._target = DomainTuple.make(grid_domain)
self._bl = bl
self._gconf = gconf
self._idx = idx
self._capability = self.TIMES | self.ADJOINT_TIMES
self._eps = float(eps)
self._uv = uv # FIXME: should we write-protect this?
def apply(self, x, mode):
from nifty_gridder import to_grid, from_grid
import nifty_gridder
self._check_input(x, mode)
if mode == self.TIMES:
nu2, nv2 = self._target.shape
res = to_grid(self._uv, x.to_global_data(), nu2, nv2, self._eps)
x = self._bl.ms2vis(x.to_global_data().reshape((-1, 1)), self._idx)
res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x)
else:
res = from_grid(self._uv, x.to_global_data(), self._eps)
res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx,
x.to_global_data())
res = self._bl.vis2ms(res, self._idx).reshape((-1,))
return from_global_data(self._tgt(mode), res)
......@@ -20,6 +20,7 @@ import numpy as np
from .field import Field
from .multi_field import MultiField
from .sugar import makeOp
from . import utilities
class Linearization(object):
......@@ -108,7 +109,6 @@ class Linearization(object):
return self._metric
def __getitem__(self, name):
from .operators.simple_linear_operators import ducktape
return self.new(self._val[name], self._jac.ducktape_left(name))
def __neg__(self):
......@@ -298,6 +298,10 @@ class Linearization(object):
tmp2 = makeOp(1. - (tmp == min) - (tmp == max))
return self.new(tmp, tmp2(self._jac))
def sqrt(self):
tmp = self._val.sqrt()
return self.new(tmp, makeOp(0.5/tmp)(self._jac))
def sin(self):
tmp = self._val.sin()
tmp2 = self._val.cos()
......@@ -315,7 +319,11 @@ class Linearization(object):
def sinc(self):
tmp = self._val.sinc()
tmp2 = (self._val.cos()-tmp)/self._val
tmp2 = ((np.pi*self._val).cos()-tmp)/self._val
ind = self._val.local_data == 0
loc = tmp2.local_data.copy()
loc[ind] = 0
tmp2 = Field.from_local_data(tmp.domain, loc)
return self.new(tmp, makeOp(tmp2)(self._jac))
def log(self):
......@@ -342,8 +350,16 @@ class Linearization(object):
return self.new(tmp2, makeOp(0.5*(1.-tmp**2))(self._jac))
def absolute(self):
if utilities.iscomplextype(self._val.dtype):
raise TypeError("Argument must not be complex")
tmp = self._val.absolute()
tmp2 = self._val.sign()
ind = self._val.local_data == 0
loc = tmp2.local_data.copy().astype(float)
loc[ind] = np.nan
tmp2 = Field.from_local_data(tmp.domain, loc)
return self.new(tmp, makeOp(tmp2)(self._jac))
def one_over(self):
......
......@@ -28,12 +28,15 @@ class Adder(Operator):
field : Field or MultiField
The field by which the input is shifted.
"""
def __init__(self, field):
def __init__(self, field, neg=False):
if not isinstance(field, (Field, MultiField)):
raise TypeError
self._field = field
self._domain = self._target = field.domain
self._neg = bool(neg)
def apply(self, x):
self._check_input(x)
if self._neg:
return x - self._field
return x + self._field
......@@ -15,19 +15,16 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domains.rg_space import RGSpace
from ..domains.lm_space import LMSpace
from ..domains.hp_space import HPSpace
from .. import utilities
from ..domain_tuple import DomainTuple
from ..domains.gl_space import GLSpace
from ..domains.hp_space import HPSpace
from ..domains.rg_space import RGSpace
from .diagonal_operator import DiagonalOperator
from .endomorphic_operator import EndomorphicOperator
from .harmonic_operators import HarmonicTransformOperator
from .diagonal_operator import DiagonalOperator
from .simple_linear_operators import WeightApplier
from ..domain_tuple import DomainTuple
from ..field import Field
from .. import utilities
def FuncConvolutionOperator(domain, func, space=None):
......@@ -77,4 +74,25 @@ def _ConvolutionOperator(domain, kernel, space=None):
HT = HarmonicTransformOperator(lm, domain[space], space)
diag = DiagonalOperator(kernel*domain[space].total_volume, lm, (space,))
wgt = WeightApplier(domain, space, 1)
return HT(diag(HT.adjoint(wgt)))
op = HT(diag(HT.adjoint(wgt)))
return _ApplicationWithoutMeanOperator(op)
class _ApplicationWithoutMeanOperator(EndomorphicOperator):
def __init__(self, op):
self._capability = self.TIMES | self.ADJOINT_TIMES
if op.domain != op.target:
raise TypeError("Operator needs to be endomorphic")
self._domain = op.domain
self._op = op
def apply(self, x, mode):
self._check_input(x, mode)
mean = x.mean()
return mean + self._op.apply(x - mean, mode)
def __repr__(self):
from ..utilities import indent
return "\n".join((
"_ApplicationWithoutMeanOperator:",
indent(self._op.__repr__())))
......@@ -30,7 +30,7 @@ class DOFDistributor(LinearOperator):
"""Operator which distributes actual degrees of freedom (dof) according to
some distribution scheme into a higher dimensional space. This distribution
scheme is defined by the dofdex, a degree of freedom index, which
associates the entries within the operators domain to locations in its
associates the entries within the operator's domain to locations in its
target. This operator's domain is a DOFSpace, which is defined according to
the distribution scheme.
......
......@@ -110,8 +110,8 @@ class GaussianEnergy(EnergyOperator):
----------
mean : Field
Mean of the Gaussian. Default is 0.
covariance : LinearOperator
Covariance of the Gaussian. Default is the identity operator.
inverse_covariance : LinearOperator
Inverse covariance of the Gaussian. Default is the identity operator.
domain : Domain, DomainTuple, tuple of Domain or MultiDomain
Operator domain. By default it is inferred from `mean` or
`covariance` if specified
......@@ -121,28 +121,27 @@ class GaussianEnergy(EnergyOperator):
At least one of the arguments has to be provided.
"""
def __init__(self, mean=None, covariance=None, domain=None):
def __init__(self, mean=None, inverse_covariance=None, domain=None):
if mean is not None and not isinstance(mean, (Field, MultiField)):
raise TypeError
if covariance is not None and not isinstance(covariance,
LinearOperator):
if inverse_covariance is not None and not isinstance(inverse_covariance, LinearOperator):
raise TypeError
self._domain = None
if mean is not None:
self._checkEquivalence(mean.domain)
if covariance is not None:
self._checkEquivalence(covariance.domain)
if inverse_covariance is not None:
self._checkEquivalence(inverse_covariance.domain)
if domain is not None:
self._checkEquivalence(domain)
if self._domain is None:
raise ValueError("no domain given")
self._mean = mean
if covariance is None:
if inverse_covariance is None:
self._op = SquaredNormOperator(self._domain).scale(0.5)
else:
self._op = QuadraticFormOperator(covariance.inverse)
self._icov = None if covariance is None else covariance.inverse
self._op = QuadraticFormOperator(inverse_covariance)
self._icov = None if inverse_covariance is None else inverse_covariance
def _checkEquivalence(self, newdom):
newdom = makeDomain(newdom)
......
......@@ -202,7 +202,7 @@ class HartleyOperator(LinearOperator):
rem_axes = tuple(i for i in axes if i != oldax)
tmp = x.val
ldat = dobj.local_data(tmp)
ldat = fft.my_fftn_r2c(ldat, axes=rem_axes)
ldat = fft.fftn(ldat, axes=rem_axes)
if oldax != 0:
raise ValueError("bad distribution")
ldat2 = ldat.reshape((ldat.shape[0],
......@@ -211,7 +211,7 @@ class HartleyOperator(LinearOperator):
tmp = dobj.from_local_data(shp2d, ldat2, distaxis=0)
tmp = dobj.transpose(tmp)
ldat2 = dobj.local_data(tmp)
ldat2 = fft.my_fftn(ldat2, axes=(1,))
ldat2 = fft.fftn(ldat2, axes=(1,))
ldat2 = ldat2.real+ldat2.imag
tmp = dobj.from_local_data(tmp.shape, ldat2, distaxis=0)
tmp = dobj.transpose(tmp)
......
......@@ -23,7 +23,6 @@ from ..sugar import domain_union
from ..utilities import indent
from .block_diagonal_operator import BlockDiagonalOperator
from .linear_operator import LinearOperator
from .simple_linear_operators import NullOperator
class SumOperator(LinearOperator):
......
......@@ -157,7 +157,6 @@ def _rgb_data(spectral_cube):
xyz_data /= xyz_data.max()
xyz_data = to_logscale(xyz_data, max(1e-3, xyz_data.min()), 1.)
rgb_data = xyz_data.copy()
it = np.nditer(xyz_data[:, 0], flags=['multi_index'])
for x in range(xyz_data.shape[0]):
rgb_data[x] = _gammacorr(np.matmul(MATRIX_SRGB_D65, xyz_data[x]))
rgb_data = rgb_data.clip(0., 1.)
......@@ -323,7 +322,8 @@ def _plot1D(f, ax, **kwargs):
plt.yscale(kwargs.pop("yscale", "log"))
xcoord = dom.k_lengths
for i, fld in enumerate(f):
ycoord = fld.to_global_data()
ycoord = fld.to_global_data_rw()
ycoord[0] = ycoord[1]
plt.plot(xcoord, ycoord, label=label[i],
linewidth=linewidth[i], alpha=alpha[i])
_limit_xy(**kwargs)
......
......@@ -15,9 +15,9 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from .field import Field
from .operators.endomorphic_operator import EndomorphicOperator
from .operators.operator import Operator
from .sugar import from_random
class StatCalculator(object):
......@@ -131,7 +131,6 @@ def probe_diagonal(op, nprobes, random_type="pm1"):
'''
sc = StatCalculator()
for i in range(nprobes):
input = Field.from_random(random_type, op.domain)
output = op(input)
sc.add(output.conjugate()*input)
x = from_random(random_type, op.domain)
sc.add(op(x).conjugate()*x)
return sc.mean
......@@ -240,6 +240,8 @@ def special_add_at(a, axis, index, b):
_iscomplex_tpl = (np.complex64, np.complex128)
def iscomplextype(dtype):
return dtype.type in _iscomplex_tpl
......
......@@ -57,8 +57,6 @@ def test_power_synthesize_analyze(space1, space2):
fp1 = ift.PS_field(p1, _spec1)
p2 = ift.PowerSpace(space2)
fp2 = ift.PS_field(p2, _spec2)
outer = np.outer(fp1.to_global_data(), fp2.to_global_data())
fp = ift.Field.from_global_data((p1, p2), outer)
op1 = ift.create_power_operator((space1, space2), _spec1, 0)
op2 = ift.create_power_operator((space1, space2), _spec2, 1)
......@@ -345,11 +343,11 @@ def test_funcs(num, dom, func):
@pmp('dtype', [np.float64, np.complex128])
def test_from_random(rtype, dtype):
sp = ift.RGSpace(3)
f = ift.Field.from_random(rtype, sp, dtype=dtype)
ift.Field.from_random(rtype, sp, dtype=dtype)
def test_field_of_objects():
arr = np.array(['x', 'y', 'z'])
sp = ift.RGSpace(3)