Commit c6e6487d authored by Philipp Arras's avatar Philipp Arras
Browse files

Remove everything related to SLAmplitudes

parent 4351c81a
Pipeline #63201 failed with stages
in 4 minutes and 10 seconds
...@@ -11,7 +11,6 @@ from .domains.gl_space import GLSpace ...@@ -11,7 +11,6 @@ from .domains.gl_space import GLSpace
from .domains.hp_space import HPSpace from .domains.hp_space import HPSpace
from .domains.power_space import PowerSpace from .domains.power_space import PowerSpace
from .domains.dof_space import DOFSpace from .domains.dof_space import DOFSpace
from .domains.log_rg_space import LogRGSpace
from .domain_tuple import DomainTuple from .domain_tuple import DomainTuple
from .multi_domain import MultiDomain from .multi_domain import MultiDomain
...@@ -26,7 +25,6 @@ from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter ...@@ -26,7 +25,6 @@ from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
from .operators.contraction_operator import ContractionOperator from .operators.contraction_operator import ContractionOperator
from .operators.linear_interpolation import LinearInterpolator from .operators.linear_interpolation import LinearInterpolator
from .operators.endomorphic_operator import EndomorphicOperator from .operators.endomorphic_operator import EndomorphicOperator
from .operators.exp_transform import ExpTransform
from .operators.harmonic_operators import ( from .operators.harmonic_operators import (
FFTOperator, HartleyOperator, SHTOperator, HarmonicTransformOperator, FFTOperator, HartleyOperator, SHTOperator, HarmonicTransformOperator,
HarmonicSmoothingOperator) HarmonicSmoothingOperator)
...@@ -34,13 +32,10 @@ from .operators.field_zero_padder import FieldZeroPadder ...@@ -34,13 +32,10 @@ from .operators.field_zero_padder import FieldZeroPadder
from .operators.inversion_enabler import InversionEnabler from .operators.inversion_enabler import InversionEnabler
from .operators.linear_operator import LinearOperator from .operators.linear_operator import LinearOperator
from .operators.mask_operator import MaskOperator from .operators.mask_operator import MaskOperator
from .operators.qht_operator import QHTOperator
from .operators.regridding_operator import RegriddingOperator from .operators.regridding_operator import RegriddingOperator
from .operators.sampling_enabler import SamplingEnabler from .operators.sampling_enabler import SamplingEnabler
from .operators.sandwich_operator import SandwichOperator from .operators.sandwich_operator import SandwichOperator
from .operators.scaling_operator import ScalingOperator from .operators.scaling_operator import ScalingOperator
from .operators.slope_operator import SlopeOperator
from .operators.symmetrizing_operator import SymmetrizingOperator
from .operators.block_diagonal_operator import BlockDiagonalOperator from .operators.block_diagonal_operator import BlockDiagonalOperator
from .operators.outer_product_operator import OuterProduct from .operators.outer_product_operator import OuterProduct
from .operators.simple_linear_operators import ( from .operators.simple_linear_operators import (
......
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
import numpy as np
from ..field import Field
from .structured_domain import StructuredDomain
class LogRGSpace(StructuredDomain):
'''Represents a logarithmic Cartesian grid.
Parameters
----------
shape : int or tuple of int
Number of grid points or numbers of gridpoints along each axis.
bindistances : float or tuple of float
Logarithmic distance between two grid points along each axis.
Equidistant spacing of bins on logarithmic scale is assumed.
t_0 : float or tuple of float
Coordinate of pixel ndim*(1,).
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
Default: False.
'''
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
def __init__(self, shape, bindistances, t_0, harmonic=False):
self._harmonic = bool(harmonic)
if np.isscalar(shape):
shape = (shape,)
self._shape = tuple(int(i) for i in shape)
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))
@property
def harmonic(self):
return self._harmonic
@property
def shape(self):
return self._shape
@property
def scalar_dvol(self):
return self._dvol
@property
def bindistances(self):
return np.array(self._bindistances)
@property
def size(self):
return np.prod(self._shape)
@property
def t_0(self):
"""np.ndarray : array of coordinates of pixel ndim*(1,)."""
return np.array(self._t_0)
def __repr__(self):
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
harmonic) partner domain of `self`, depending on `self.harmonic`. The
`bindistances` are transformed and `t_0` stays the same.
Returns
-------
LogRGSpace
The partner domain
"""
codomain_bindistances = 1./(self.bindistances*self.shape)
return LogRGSpace(self.shape, codomain_bindistances, self._t_0, not self.harmonic)
def get_k_length_array(self):
"""Generates array of distances to origin of the space.
Returns
-------
numpy.ndarray
Distances to origin of the space. If any index of the array is
zero then the distance is np.nan if self.harmonic True.
The dtype is float64, the shape is `self.shape`.
Raises
------
NotImplementedError
If `self.harmonic` is False.
"""
if not self.harmonic:
raise NotImplementedError
ks = self.get_k_array()
return Field.from_global_data(self, np.linalg.norm(ks, axis=0))
def get_k_array(self):
"""Generates coordinates of the space.
Returns
-------
numpy.ndarray
Coordinates of the space. If one index of the array is zero the
corresponding coordinate is -np.inf (np.nan) if self.harmonic is
False (True).
The dtype is float64 and shape: `(len(self.shape),) + self.shape`.
"""
ndim = len(self.shape)
k_array = np.zeros((ndim,) + self.shape)
dist = self.bindistances
for i in range(ndim):
ks = np.zeros(self.shape[i])
ks[1:] = np.minimum(self.shape[i] - 1 - np.arange(self.shape[i]-1),
np.arange(self.shape[i]-1)) * dist[i]
if self.harmonic:
ks[0] = np.nan
else:
ks[0] = -np.inf
ks[1:] += self.t_0[i]
k_array[i] += ks.reshape((1,)*i + (self.shape[i],)
+ (1,)*(ndim-i-1))
return k_array
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from .. import dobj
from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace
from ..domains.rg_space import RGSpace
from ..field import Field
from ..utilities import infer_space, special_add_at
from .linear_operator import LinearOperator
class ExpTransform(LinearOperator):
"""Transforms log-space to target
This operator creates a log-space subject to the degrees of freedom and
and its target-domain.
Then it transforms between this log-space and its target, which is defined
in normal units.
FIXME Write something on t_0 of domain space
E.g: A field in log-log-space can be transformed into log-norm-space,
that is the y-axis stays logarithmic, but the x-axis is transformed.
Parameters
----------
target : domain, tuple of domains or DomainTuple
The full output domain
dof : int
The degrees of freedom of the log-domain, i.e. the number of bins.
"""
def __init__(self, target, dof, space=0):
self._target = DomainTuple.make(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._space = infer_space(self._target, space)
tgt = self._target[self._space]
if not ((isinstance(tgt, RGSpace) and tgt.harmonic) or
isinstance(tgt, PowerSpace)):
raise ValueError(
"Target must be a harmonic RGSpace or a power space.")
ndim = len(tgt.shape)
if np.isscalar(dof):
dof = np.full(ndim, int(dof), dtype=np.int)
dof = np.array(dof)
t_mins = np.empty(ndim)
bindistances = np.empty(ndim)
self._bindex = [None] * ndim
self._frac = [None] * ndim
for i in range(ndim):
if isinstance(tgt, RGSpace):
rng = np.arange(tgt.shape[i])
tmp = np.minimum(rng, tgt.shape[i]+1-rng)
k_array = tmp * tgt.distances[i]
else:
k_array = tgt.k_lengths
# avoid taking log of first entry
log_k_array = np.log(k_array[1:])
# Interpolate log_k_array linearly
t_max = np.max(log_k_array)
t_min = np.min(log_k_array)
# Save t_min for later
t_mins[i] = t_min
bindistances[i] = (t_max-t_min) / (dof[i]-1)
coord = np.append(0., 1. + (log_k_array-t_min) / bindistances[i])
self._bindex[i] = np.floor(coord).astype(int)
# Interpolated value is computed via
# (1.-frac)*<value from this bin> + frac*<value from next bin>
# 0 <= frac < 1.
self._frac[i] = coord - self._bindex[i]
from ..domains.log_rg_space import LogRGSpace
log_space = LogRGSpace(2*dof+1, bindistances,
t_mins, harmonic=False)
self._domain = [dom for dom in self._target]
self._domain[self._space] = log_space
self._domain = DomainTuple.make(self._domain)
def apply(self, x, mode):
self._check_input(x, mode)
v = x.val
ndim = len(self.target.shape)
curshp = list(self._dom(mode).shape)
tgtshp = self._tgt(mode).shape
d0 = self._target.axes[self._space][0]
for d in self._target.axes[self._space]:
idx = (slice(None),) * d
wgt = self._frac[d-d0].reshape((1,)*d + (-1,) + (1,)*(ndim-d-1))
v, x = dobj.ensure_not_distributed(v, (d,))
if mode == self.ADJOINT_TIMES:
shp = list(x.shape)
shp[d] = tgtshp[d]
xnew = np.zeros(shp, dtype=x.dtype)
xnew = special_add_at(xnew, d, self._bindex[d-d0], x*(1.-wgt))
xnew = special_add_at(xnew, d, self._bindex[d-d0]+1, x*wgt)
else: # TIMES
xnew = x[idx + (self._bindex[d-d0],)] * (1.-wgt)
xnew += x[idx + (self._bindex[d-d0]+1,)] * wgt
curshp[d] = xnew.shape[d]
v = dobj.from_local_data(curshp, xnew, distaxis=dobj.distaxis(v))
return Field(self._tgt(mode), dobj.ensure_default_distributed(v))
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from .. import dobj
from .. import fft
from ..domain_tuple import DomainTuple
from ..field import Field
from ..utilities import infer_space
from .linear_operator import LinearOperator
class QHTOperator(LinearOperator):
"""Does a Hartley transform on LogRGSpace
This operator takes a field on a LogRGSpace and performs a Hartley
transform. The zero modes are not transformed because they are infinitely
far away in LogRGSpaces.
Parameters
----------
target : domain, tuple of domains or DomainTuple
The full output domain
space : int
The index of the domain on which the operator acts.
target[space] must be a non-harmonic LogRGSpace.
codomain : Domain
The codomain for target[space]. If not supplied, it is inferred.
"""
def __init__(self, target, space=0, codomain=None):
self._target = DomainTuple.make(target)
self._space = infer_space(self._target, space)
from ..domains.log_rg_space import LogRGSpace
if not isinstance(self._target[self._space], LogRGSpace):
raise ValueError("target[space] has to be a LogRGSpace!")
if self._target[self._space].harmonic:
raise TypeError("target[space] must be a nonharmonic space")
self._domain = [dom for dom in self._target]
if codomain is None:
codomain = self._target[self._space].get_default_codomain()
self._domain[self._space] = codomain
self._domain = DomainTuple.make(self._domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
dom = self._domain[self._space]
v = x.val * dom.scalar_dvol
n = self._domain.axes[self._space]
rng = n if mode == self.TIMES else reversed(n)
for i in rng:
sl = (slice(None),)*i + (slice(1, None),)
v, tmp = dobj.ensure_not_distributed(v, (i,))
tmp[sl] = fft.hartley(tmp[sl], axes=(i,))
return Field(self._tgt(mode), dobj.ensure_default_distributed(v))
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.log_rg_space import LogRGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from .linear_operator import LinearOperator
class SlopeOperator(LinearOperator):
"""Evaluates a line on a LogRGSpace given slope and y-intercept
Slope and y-intercept of this line are the two parameters which are
defined on an UnstructeredDomain (in this order) which is the domain of
the operator. Being a LogRGSpace instance each pixel has a well-defined
coordinate value.
The y-intercept is defined to be the value at t_0 of the target.
Parameters
----------
target : LogRGSpace
The target of the operator which needs to be one-dimensional.
"""
def __init__(self, target):
self._target = DomainTuple.make(target)
if len(self._target) > 1:
raise TypeError
if len(self._target[0].shape) > 1:
raise TypeError
if not isinstance(self._target[0], LogRGSpace):
raise TypeError
self._domain = DomainTuple.make(UnstructuredDomain((2,)))
self._capability = self.TIMES | self.ADJOINT_TIMES
pos = self.target[0].get_k_array() - self.target[0].t_0[0]
self._pos = pos[0, 1:]
def apply(self, x, mode):
self._check_input(x, mode)
inp = x.to_global_data()
if mode == self.TIMES:
res = np.empty(self.target.shape, dtype=x.dtype)
res[0] = 0
res[1:] = inp[1] + inp[0]*self._pos
else:
res = np.array(
[np.sum(self._pos*inp[1:]),
np.sum(inp[1:])], dtype=x.dtype)
return Field.from_global_data(self._tgt(mode), res)
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# 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 <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from .. import dobj, utilities
from ..domain_tuple import DomainTuple
from ..domains.log_rg_space import LogRGSpace
from ..field import Field
from .endomorphic_operator import EndomorphicOperator
class SymmetrizingOperator(EndomorphicOperator):
"""Subtracts the field axes-wise in reverse order from itself. The slice of
all elements with at least one index being zero is not touched.
Parameters
----------
domain : Domain, DomainTuple or tuple of Domain
Domain of the operator. `domain[space]` needs to be a non-harmonic
:class:`LogRGSpace`.
space : int
Index of space in domain on which the operator shall act. Default is 0.
"""
def __init__(self, domain, space=0):
self._domain = DomainTuple.make(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._space = utilities.infer_space(self._domain, space)
dom = self._domain[self._space]
if not (isinstance(dom, LogRGSpace) and not dom.harmonic):
raise TypeError("nonharmonic LogRGSpace needed")
def apply(self, x, mode):
self._check_input(x, mode)
v = x.val.copy()
for i in self._domain.axes[self._space]:
lead = (slice(None),)*i
v, loc = dobj.ensure_not_distributed(v, (i,))
loc[lead+(slice(1, None),)] -= loc[lead+(slice(None, 0, -1),)]
loc /= 2
return Field(self.target, dobj.ensure_default_distributed(v))
...@@ -22,7 +22,6 @@ import numpy as np ...@@ -22,7 +22,6 @@ import numpy as np
from . import dobj from . import dobj
from .domains.gl_space import GLSpace from .domains.gl_space import GLSpace
from .domains.hp_space import HPSpace from .domains.hp_space import HPSpace
from .domains.log_rg_space import LogRGSpace
from .domains.power_space import PowerSpace from .domains.power_space import PowerSpace
from .domains.rg_space import RGSpace from .domains.rg_space import RGSpace
from .field import Field from .field import Field
...@@ -306,18 +305,6 @@ def _plot1D(f, ax, **kwargs): ...@@ -306,18 +305,6 @@ def _plot1D(f, ax, **kwargs):
if label != ([None]*len(f)): if label != ([None]*len(f)):
plt.legend() plt.legend()
return return
elif isinstance(dom, LogRGSpace):
plt.yscale(kwargs.pop("yscale", "log"))
npoints = dom.shape[0]
xcoord = dom.t_0 + np.arange(npoints-1)*dom.bindistances[0]
for i, fld in enumerate(f):
ycoord = fld.to_global_data()[1:]
plt.plot(xcoord, ycoord, label=label[i],
linewidth=linewidth[i], alpha=alpha[i])
_limit_xy(**kwargs)
if label != ([None]*len(f)):
plt.legend()
return
elif isinstance(dom, PowerSpace): elif isinstance(dom, PowerSpace):
plt.xscale(kwargs.pop("xscale", "log")) plt.xscale(kwargs.pop("xscale", "log"))
plt.yscale(kwargs.pop("yscale", "log")) plt.yscale(kwargs.pop("yscale", "log"))
...@@ -432,8 +419,8 @@ def _plot(f, ax, **kwargs): ...@@ -432,8 +419,8 @@ def _plot(f, ax, **kwargs):
dom1 = f[0].domain dom1 = f[0].domain
if (len(dom1) == 1 and if (len(dom1) == 1 and
(isinstance(dom1[0], PowerSpace) or (isinstance(dom1[0], PowerSpace) or
(isinstance(dom1[0], (RGSpace, LogRGSpace)) and (isinstance(dom1[0], RGSpace) and
len(dom1[0].shape) == 1))): len(dom1[0].shape) == 1))):
_plot1D(f, ax, **kwargs) _plot1D(f, ax, **kwargs)
return return
else: else:
......
...@@ -91,16 +91,6 @@ def testConjugationOperator(sp): ...@@ -91,16 +91,6 @@ def testConjugationOperator(sp):
only_r_linear=True) only_r_linear=True)
@pmp('args', [(ift.RGSpace(10, harmonic=True), 4, 0), (ift.RGSpace(
(24, 31), distances=(0.4, 2.34), harmonic=True), 3, 0),
(ift.LMSpace(4), 10, 0)])
def testSlopeOperator(args, dtype):
tmp = ift.ExpTransform(ift.PowerSpace(args[0]), args[1], args[2])
tgt = tmp.domain[0]
op = ift.SlopeOperator(tgt)
ift.extra.consistency_check(op, dtype, dtype)
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces) @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
def testOperatorAdaptor(sp, dtype): def testOperatorAdaptor(sp, dtype):
op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype)) op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
...@@ -218,14 +208,6 @@ def testDomainTupleFieldInserter(): ...@@ -218,14 +208,6 @@ def testDomainTupleFieldInserter():
ift.extra.consistency_check(op) ift.extra.consistency_check(op)
@pmp('space', [0, 2])
def testSymmetrizingOperator(space, dtype):
dom = (ift.LogRGSpace(10, [2.], [1.]), ift.UnstructuredDomain(13),
ift.LogRGSpace((5, 27), [1., 2.7], [0., 4.]), ift.HPSpace(4))
op = ift.SymmetrizingOperator(dom, space)
ift.extra.consistency_check(op, dtype, dtype)