diff --git a/nifty5/__init__.py b/nifty5/__init__.py index b55860393f92a7c86c380f524c4ecc88c00603f5..678dd684965aa8fde0a24dd45e281ecbe6ed4d43 100644 --- a/nifty5/__init__.py +++ b/nifty5/__init__.py @@ -32,8 +32,9 @@ from .operators.domain_distributor import DomainDistributor from .operators.endomorphic_operator import EndomorphicOperator from .operators.exp_transform import ExpTransform from .operators.fft_operator import FFTOperator -from .operators.fft_smoothing_operator import FFTSmoothingOperator from .operators.field_zero_padder import FieldZeroPadder +from .operators.hartley_operator import HartleyOperator +from .operators.harmonic_smoothing_operator import HarmonicSmoothingOperator from .operators.geometry_remover import GeometryRemover from .operators.harmonic_transform_operator import HarmonicTransformOperator from .operators.inversion_enabler import InversionEnabler diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py index d4798e4e5be7612db9e4b40fbf426be32805a7ab..b2c5e2078aee1abbb11f4f896c8c2166ff6bea55 100644 --- a/nifty5/library/correlated_fields.py +++ b/nifty5/library/correlated_fields.py @@ -25,7 +25,7 @@ from ..models.local_nonlinearity import PointwiseExponential from ..models.variable import Variable from ..multi.multi_field import MultiField from ..operators.domain_distributor import DomainDistributor -from ..operators.fft_operator import FFTOperator +from ..operators.hartley_operator import HartleyOperator from ..operators.harmonic_transform_operator import HarmonicTransformOperator from ..operators.power_distributor import PowerDistributor @@ -41,7 +41,7 @@ def make_correlated_field(s_space, amplitude_model): amplitude_model : model for correlation structure ''' h_space = s_space.get_default_codomain() - ht = FFTOperator(h_space, s_space) + ht = HartleyOperator(h_space, s_space) p_space = amplitude_model.value.domain[0] power_distributor = PowerDistributor(h_space, p_space) diff --git a/nifty5/operators/domain_distributor.py b/nifty5/operators/domain_distributor.py index 622214c480cd7f71ca2d6b9775980ceefd6e8c62..eb1899c76f7fb9f34b8388738c93e3d1ff36e043 100644 --- a/nifty5/operators/domain_distributor.py +++ b/nifty5/operators/domain_distributor.py @@ -25,26 +25,16 @@ from ..compat import * from ..domain_tuple import DomainTuple from ..field import Field from .linear_operator import LinearOperator +from .. import utilities -# MR FIXME: this needs to be rewritten in a generic fashion class DomainDistributor(LinearOperator): - def __init__(self, target, axis): - if dobj.ntask > 1: - raise NotImplementedError('UpProj class does not support MPI.') - assert len(target) == 2 - assert axis in [0, 1] - - if axis == 0: - domain = target[1] - self._size = target[0].size - else: - domain = target[0] - self._size = target[1].size - - self._axis = axis - self._domain = DomainTuple.make(domain) + def __init__(self, target, spaces): self._target = DomainTuple.make(target) + self._spaces = utilities.parse_spaces(spaces, len(self._target)) + self._domain = [tgt for i, tgt in enumerate(self._target) + if i in self._spaces] + self._domain = DomainTuple.make(self._domain) @property def domain(self): @@ -57,23 +47,16 @@ class DomainDistributor(LinearOperator): def apply(self, x, mode): self._check_input(x, mode) if mode == self.TIMES: - x = x.local_data - otherDirection = np.ones(self._size) - if self._axis == 0: - res = np.outer(otherDirection, x) - else: - res = np.outer(x, otherDirection) - res = res.reshape(dobj.local_shape(self.target.shape)) - return Field.from_local_data(self.target, res) + ldat = x.local_data if 0 in self._spaces else x.to_global_data() + shp = [] + for i, tgt in enumerate(self._target): + tmp = tgt.shape if i > 0 else tgt.local_shape + shp += tmp if i in self._spaces else(1,)*len(tgt.shape) + ldat = np.broadcast_to(ldat.reshape(shp), self._target.local_shape) + return Field.from_local_data(self._target, ldat) else: - if self._axis == 0: - x = x.local_data.reshape(self._size, -1) - res = np.sum(x, axis=0) - else: - x = x.local_data.reshape(-1, self._size) - res = np.sum(x, axis=1) - res = res.reshape(dobj.local_shape(self.domain.shape)) - return Field.from_local_data(self.domain, res) + return x.sum([s for s in range(len(x.domain)) + if s not in self._spaces]) @property def capability(self): diff --git a/nifty5/operators/exp_transform.py b/nifty5/operators/exp_transform.py index 311c1789fb1c05b5298ff5774a9951bbfee09164..0c39f21592b9bec123d698a446ba32ee5fdd9f0f 100644 --- a/nifty5/operators/exp_transform.py +++ b/nifty5/operators/exp_transform.py @@ -27,12 +27,13 @@ from ..domains.power_space import PowerSpace from ..domains.rg_space import RGSpace from ..field import Field from .linear_operator import LinearOperator +from .. import utilities class ExpTransform(LinearOperator): def __init__(self, target, dof, space=0): self._target = DomainTuple.make(target) - self._space = int(space) + self._space = utilities.infer_space(self._target, space) tgt = self._target[self._space] if not ((isinstance(tgt, RGSpace) and tgt.harmonic) or isinstance(tgt, PowerSpace)): diff --git a/nifty5/operators/fft_operator.py b/nifty5/operators/fft_operator.py index c96687005c1c39f5dcca19066c37bfdd28616b7f..684e55530f62763776cdcad9458393f7b8508374 100644 --- a/nifty5/operators/fft_operator.py +++ b/nifty5/operators/fft_operator.py @@ -43,6 +43,13 @@ class FFTOperator(LinearOperator): The index of the subdomain on which the operator should act If None, it is set to 0 if `domain` contains exactly one space. `domain[space]` must be an RGSpace. + + Notes + ----- + This operator performs full FFTs, which implies that its output field will + always have complex type, regardless of the type of the input field. + If a real field is desired after a forward/backward transform couple, it + must be manually cast to real. """ def __init__(self, domain, target=None, space=None): @@ -67,41 +74,35 @@ class FFTOperator(LinearOperator): utilities.fft_prep() def apply(self, x, mode): + from pyfftw.interfaces.numpy_fft import fftn, ifftn self._check_input(x, mode) - if np.issubdtype(x.dtype, np.complexfloating): - return (self._apply_cartesian(x.real, mode) + - 1j*self._apply_cartesian(x.imag, mode)) + ncells = x.domain[self._space].size + if x.domain[self._space].harmonic: # harmonic -> position + func = fftn + fct = 1. else: - return self._apply_cartesian(x, mode) - - def _apply_cartesian(self, x, mode): + func = ifftn + fct = ncells axes = x.domain.axes[self._space] tdom = self._tgt(mode) oldax = dobj.distaxis(x.val) if oldax not in axes: # straightforward, no redistribution needed ldat = x.local_data - ldat = utilities.hartley(ldat, axes=axes) + ldat = func(ldat, axes=axes) tmp = dobj.from_local_data(x.val.shape, ldat, distaxis=oldax) elif len(axes) < len(x.shape) or len(axes) == 1: - # we can use one Hartley pass in between the redistributions + # we can use one FFT pass in between the redistributions tmp = dobj.redistribute(x.val, nodist=axes) newax = dobj.distaxis(tmp) ldat = dobj.local_data(tmp) - ldat = utilities.hartley(ldat, axes=axes) + ldat = func(ldat, axes=axes) tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=newax) tmp = dobj.redistribute(tmp, dist=oldax) - else: # two separate, full FFTs needed - # ideal strategy for the moment would be: - # - do real-to-complex FFT on all local axes - # - fill up array - # - redistribute array - # - do complex-to-complex FFT on remaining axis - # - add re+im - # - redistribute back + else: # two separate FFTs needed rem_axes = tuple(i for i in axes if i != oldax) tmp = x.val ldat = dobj.local_data(tmp) - ldat = utilities.my_fftn_r2c(ldat, axes=rem_axes) + ldat = func(ldat, axes=rem_axes) if oldax != 0: raise ValueError("bad distribution") ldat2 = ldat.reshape((ldat.shape[0], @@ -110,17 +111,16 @@ class FFTOperator(LinearOperator): tmp = dobj.from_local_data(shp2d, ldat2, distaxis=0) tmp = dobj.transpose(tmp) ldat2 = dobj.local_data(tmp) - ldat2 = utilities.my_fftn(ldat2, axes=(1,)) - ldat2 = ldat2.real+ldat2.imag + ldat2 = func(ldat2, axes=(1,)) tmp = dobj.from_local_data(tmp.shape, ldat2, distaxis=0) tmp = dobj.transpose(tmp) ldat2 = dobj.local_data(tmp).reshape(ldat.shape) tmp = dobj.from_local_data(x.val.shape, ldat2, distaxis=0) Tval = Field(tdom, tmp) if mode & (LinearOperator.TIMES | LinearOperator.ADJOINT_TIMES): - fct = self._domain[self._space].scalar_dvol + fct *= self._domain[self._space].scalar_dvol else: - fct = self._target[self._space].scalar_dvol + fct *= self._target[self._space].scalar_dvol return Tval if fct == 1 else Tval*fct @property diff --git a/nifty5/operators/field_zero_padder.py b/nifty5/operators/field_zero_padder.py index eb2b245631a67a627ef19a986fd828eb1c267b37..5074aade57910913d7584247ed23425ddc348473 100644 --- a/nifty5/operators/field_zero_padder.py +++ b/nifty5/operators/field_zero_padder.py @@ -8,13 +8,14 @@ from ..domain_tuple import DomainTuple from ..domains.rg_space import RGSpace from ..field import Field from .linear_operator import LinearOperator +from .. import utilities class FieldZeroPadder(LinearOperator): def __init__(self, domain, factor, space=0): super(FieldZeroPadder, self).__init__() self._domain = DomainTuple.make(domain) - self._space = int(space) + self._space = utilities.infer_space(self._domain, space) dom = self._domain[self._space] if not isinstance(dom, RGSpace): raise TypeError("RGSpace required") @@ -52,11 +53,11 @@ class FieldZeroPadder(LinearOperator): curax = dobj.distaxis(x) if mode == self.ADJOINT_TIMES: - newarr = np.empty(dobj.local_shape(shp_out), dtype=x.dtype) + newarr = np.empty(dobj.local_shape(shp_out, curax), dtype=x.dtype) newarr[()] = dobj.local_data(x)[(slice(None),)*ax + (slice(0, shp_out[ax]),)] else: - newarr = np.zeros(dobj.local_shape(shp_out), dtype=x.dtype) + newarr = np.zeros(dobj.local_shape(shp_out, curax), dtype=x.dtype) newarr[(slice(None),)*ax + (slice(0, shp_in[ax]),)] = dobj.local_data(x) newarr = dobj.from_local_data(shp_out, newarr, distaxis=curax) diff --git a/nifty5/operators/fft_smoothing_operator.py b/nifty5/operators/harmonic_smoothing_operator.py similarity index 90% rename from nifty5/operators/fft_smoothing_operator.py rename to nifty5/operators/harmonic_smoothing_operator.py index c849a471f441e572e095401676e3b6d43168c942..542738e2bc0c0b4eeb13f514c65a2f90d4039d86 100644 --- a/nifty5/operators/fft_smoothing_operator.py +++ b/nifty5/operators/harmonic_smoothing_operator.py @@ -22,11 +22,11 @@ from ..compat import * from ..domain_tuple import DomainTuple from ..utilities import infer_space from .diagonal_operator import DiagonalOperator -from .fft_operator import FFTOperator +from .hartley_operator import HartleyOperator from .scaling_operator import ScalingOperator -def FFTSmoothingOperator(domain, sigma, space=None): +def HarmonicSmoothingOperator(domain, sigma, space=None): """ This function returns an operator that carries out a smoothing with a Gaussian kernel of width `sigma` on the part of `domain` given by `space`. @@ -59,12 +59,12 @@ def FFTSmoothingOperator(domain, sigma, space=None): space = infer_space(domain, space) if domain[space].harmonic: raise TypeError("domain must not be harmonic") - FFT = FFTOperator(domain, space=space) - codomain = FFT.domain[space].get_default_codomain() + Hartley = HartleyOperator(domain, space=space) + codomain = Hartley.domain[space].get_default_codomain() kernel = codomain.get_k_length_array() smoother = codomain.get_fft_smoothing_kernel_function(sigma) kernel = smoother(kernel) ddom = list(domain) ddom[space] = codomain diag = DiagonalOperator(kernel, ddom, space) - return FFT.inverse*diag*FFT + return Hartley.inverse*diag*Hartley diff --git a/nifty5/operators/harmonic_transform_operator.py b/nifty5/operators/harmonic_transform_operator.py index 9233ac5dd1b325450a437cda1528b1f8fed2e72b..b62d12528f0997c8cc498cb4de2623d6f594306e 100644 --- a/nifty5/operators/harmonic_transform_operator.py +++ b/nifty5/operators/harmonic_transform_operator.py @@ -22,7 +22,7 @@ from .. import utilities from ..compat import * from ..domain_tuple import DomainTuple from ..domains.rg_space import RGSpace -from .fft_operator import FFTOperator +from .hartley_operator import HartleyOperator from .linear_operator import LinearOperator from .sht_operator import SHTOperator @@ -65,7 +65,7 @@ class HarmonicTransformOperator(LinearOperator): raise TypeError( "HarmonicTransformOperator only works on a harmonic space") if isinstance(hspc, RGSpace): - self._op = FFTOperator(domain, target, space) + self._op = HartleyOperator(domain, target, space) else: self._op = SHTOperator(domain, target, space) diff --git a/nifty5/operators/hartley_operator.py b/nifty5/operators/hartley_operator.py new file mode 100644 index 0000000000000000000000000000000000000000..eed423089f111047adf13fc217dfcb201be75952 --- /dev/null +++ b/nifty5/operators/hartley_operator.py @@ -0,0 +1,148 @@ +# 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-2018 Max-Planck-Society +# +# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik +# and financially supported by the Studienstiftung des deutschen Volkes. + +from __future__ import absolute_import, division, print_function + +import numpy as np + +from .. import dobj, utilities +from ..compat import * +from ..domain_tuple import DomainTuple +from ..domains.rg_space import RGSpace +from ..field import Field +from .linear_operator import LinearOperator + + +class HartleyOperator(LinearOperator): + """Transforms between a pair of position and harmonic RGSpaces. + + Parameters + ---------- + domain: Domain, tuple of Domain or DomainTuple + The domain of the data that is input by "times" and output by + "adjoint_times". + target: Domain, optional + The target (sub-)domain of the transform operation. + If omitted, a domain will be chosen automatically. + space: int, optional + The index of the subdomain on which the operator should act + If None, it is set to 0 if `domain` contains exactly one space. + `domain[space]` must be an RGSpace. + + Notes + ----- + This operator always produces output fields with the same data type as + its input. This is achieved by performing so-called Hartley transforms + (https://en.wikipedia.org/wiki/Discrete_Hartley_transform). + For complex input fields, the operator will transform the real and + imaginary parts separately and use the results as real and imaginary parts + of the result field, respectivey. + In many contexts the Hartley transform is a perfect substitute for the + Fourier transform, but in some situations (e.g. convolution with a general, + non-symmetrc kernel, the full FFT must be used instead. + """ + + def __init__(self, domain, target=None, space=None): + super(HartleyOperator, self).__init__() + + # Initialize domain and target + self._domain = DomainTuple.make(domain) + self._space = utilities.infer_space(self._domain, space) + + adom = self._domain[self._space] + if not isinstance(adom, RGSpace): + raise TypeError("HartleyOperator only works on RGSpaces") + if target is None: + target = adom.get_default_codomain() + + self._target = [dom for dom in self._domain] + self._target[self._space] = target + self._target = DomainTuple.make(self._target) + adom.check_codomain(target) + target.check_codomain(adom) + + utilities.fft_prep() + + def apply(self, x, mode): + self._check_input(x, mode) + if np.issubdtype(x.dtype, np.complexfloating): + return (self._apply_cartesian(x.real, mode) + + 1j*self._apply_cartesian(x.imag, mode)) + else: + return self._apply_cartesian(x, mode) + + def _apply_cartesian(self, x, mode): + axes = x.domain.axes[self._space] + tdom = self._tgt(mode) + oldax = dobj.distaxis(x.val) + if oldax not in axes: # straightforward, no redistribution needed + ldat = x.local_data + ldat = utilities.hartley(ldat, axes=axes) + tmp = dobj.from_local_data(x.val.shape, ldat, distaxis=oldax) + elif len(axes) < len(x.shape) or len(axes) == 1: + # we can use one Hartley pass in between the redistributions + tmp = dobj.redistribute(x.val, nodist=axes) + newax = dobj.distaxis(tmp) + ldat = dobj.local_data(tmp) + ldat = utilities.hartley(ldat, axes=axes) + tmp = dobj.from_local_data(tmp.shape, ldat, distaxis=newax) + tmp = dobj.redistribute(tmp, dist=oldax) + else: # two separate, full FFTs needed + # ideal strategy for the moment would be: + # - do real-to-complex FFT on all local axes + # - fill up array + # - redistribute array + # - do complex-to-complex FFT on remaining axis + # - add re+im + # - redistribute back + rem_axes = tuple(i for i in axes if i != oldax) + tmp = x.val + ldat = dobj.local_data(tmp) + ldat = utilities.my_fftn_r2c(ldat, axes=rem_axes) + if oldax != 0: + raise ValueError("bad distribution") + ldat2 = ldat.reshape((ldat.shape[0], + np.prod(ldat.shape[1:]))) + shp2d = (x.val.shape[0], np.prod(x.val.shape[1:])) + tmp = dobj.from_local_data(shp2d, ldat2, distaxis=0) + tmp = dobj.transpose(tmp) + ldat2 = dobj.local_data(tmp) + ldat2 = utilities.my_fftn(ldat2, axes=(1,)) + ldat2 = ldat2.real+ldat2.imag + tmp = dobj.from_local_data(tmp.shape, ldat2, distaxis=0) + tmp = dobj.transpose(tmp) + ldat2 = dobj.local_data(tmp).reshape(ldat.shape) + tmp = dobj.from_local_data(x.val.shape, ldat2, distaxis=0) + Tval = Field(tdom, tmp) + if mode & (LinearOperator.TIMES | LinearOperator.ADJOINT_TIMES): + fct = self._domain[self._space].scalar_dvol + else: + fct = self._target[self._space].scalar_dvol + return Tval if fct == 1 else Tval*fct + + @property + def domain(self): + return self._domain + + @property + def target(self): + return self._target + + @property + def capability(self): + return self._all_ops diff --git a/nifty5/operators/qht_operator.py b/nifty5/operators/qht_operator.py index 13b0607e3b72ee6b92c9c2b723507cebd378fb1c..3eebc338146c167b3ebc40b837291b00accccf0d 100644 --- a/nifty5/operators/qht_operator.py +++ b/nifty5/operators/qht_operator.py @@ -22,7 +22,7 @@ from .. import dobj from ..compat import * from ..domain_tuple import DomainTuple from ..field import Field -from ..utilities import hartley +from ..utilities import hartley, infer_space from .linear_operator import LinearOperator @@ -47,7 +47,7 @@ class QHTOperator(LinearOperator): """ def __init__(self, domain, target, space=0): self._domain = DomainTuple.make(domain) - self._space = int(space) + self._space = infer_space(self._domain, space) from ..domains.log_rg_space import LogRGSpace if not isinstance(self._domain[self._space], LogRGSpace): @@ -57,7 +57,7 @@ class QHTOperator(LinearOperator): if not self._domain[self._space].harmonic: raise TypeError( - "HarmonicTransformOperator only works on a harmonic space") + "QHTOperator only works on a harmonic space") if target.harmonic: raise TypeError("Target is not a codomain of domain") diff --git a/nifty5/operators/symmetrizing_operator.py b/nifty5/operators/symmetrizing_operator.py index d59a6cf66e0c139ee07fe0930a0d9c1ffcb6ddcb..8a2aa881e6c2b3e49cdba19275a2d9b16b54b89a 100644 --- a/nifty5/operators/symmetrizing_operator.py +++ b/nifty5/operators/symmetrizing_operator.py @@ -24,15 +24,16 @@ from ..domain_tuple import DomainTuple from ..domains.log_rg_space import LogRGSpace from ..field import Field from .endomorphic_operator import EndomorphicOperator +from .. import utilities class SymmetrizingOperator(EndomorphicOperator): def __init__(self, domain, space=0): self._domain = DomainTuple.make(domain) - self._space = int(space) + self._space = utilities.infer_space(self._domain, space) dom = self._domain[self._space] if not (isinstance(dom, LogRGSpace) and not dom.harmonic): - raise TypeError + raise TypeError("nonharmonic LogRGSpace needed") @property def domain(self): diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py index 6c93c69fa7cfa577ef809c379a16779861535a3d..f6f9e9542e5cd41d35eb7d9444164e663b2513cd 100644 --- a/test/test_operators/test_adjoint.py +++ b/test/test_operators/test_adjoint.py @@ -56,6 +56,14 @@ class Consistency_Tests(unittest.TestCase): op = ift.FFTOperator(sp.get_default_codomain()) ift.extra.consistency_check(op, dtype, dtype) + @expand(product(_h_RG_spaces+_p_RG_spaces, + [np.float64, np.complex128])) + def testHartley(self, sp, dtype): + op = ift.HartleyOperator(sp) + ift.extra.consistency_check(op, dtype, dtype) + op = ift.HartleyOperator(sp.get_default_codomain()) + ift.extra.consistency_check(op, dtype, dtype) + @expand(product(_h_spaces, [np.float64, np.complex128])) def testHarmonic(self, sp, dtype): op = ift.HarmonicTransformOperator(sp) @@ -93,3 +101,50 @@ class Consistency_Tests(unittest.TestCase): def testGeometryRemover(self, sp, dtype): op = ift.GeometryRemover(sp) ift.extra.consistency_check(op, dtype, dtype) + + @expand(product([0, 1, 2, 3, (0, 1), (0, 2), (0, 1, 2), (0, 2, 3), (1, 3)], + [np.float64, np.complex128])) + def testDomainDistributor(self, spaces, dtype): + dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.GLSpace(5), + ift.HPSpace(4)) + op = ift.DomainDistributor(dom, spaces) + ift.extra.consistency_check(op, dtype, dtype) + + @expand(product([0, 2], [np.float64, np.complex128])) + def testSymmetrizingOperator(self, 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) + + @expand(product([0, 2], [2, 2.7], [np.float64, np.complex128])) + def testZeroPadder(self, space, factor, dtype): + dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.RGSpace(7), + ift.HPSpace(4)) + op = ift.FieldZeroPadder(dom, factor, space) + ift.extra.consistency_check(op, dtype, dtype) + + @expand(product([(ift.RGSpace(10, harmonic=True), 4, 0), + (ift.RGSpace((24, 31), distances=(0.4, 2.34), + harmonic=True), (4, 3), 0), + ((ift.HPSpace(4), ift.RGSpace(27, distances=0.3, + harmonic=True)), (10,), 1), + (ift.PowerSpace(ift.RGSpace(10, distances=0.3, + harmonic=True)), 6, 0)], + [np.float64, np.complex128])) + def testExpTransform(self, args, dtype): + op = ift.ExpTransform(args[0], args[1], args[2]) + ift.extra.consistency_check(op, dtype, dtype) + + @expand(product([(ift.LogRGSpace([10, 17], [2., 3.], [1., 0.]), 0), + ((ift.LogRGSpace(10, [2.], [1.]), + ift.UnstructuredDomain(13)), 0), + ((ift.UnstructuredDomain(13), + ift.LogRGSpace(17, [3.], [.7])), 1)], + [np.float64])) + def testQHTOperator(self, args, dtype): + dom = ift.DomainTuple.make(args[0]) + tgt = list(dom) + tgt[args[1]] = tgt[args[1]].get_default_codomain() + op = ift.QHTOperator(tgt, dom[args[1]], args[1]) + ift.extra.consistency_check(op, dtype, dtype) diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py index 24ec9fe1f74246ed8710388c45d06bf8ed4ec322..783d0a7f9655771eca08d1460ea5970e65547da2 100644 --- a/test/test_operators/test_fft_operator.py +++ b/test/test_operators/test_fft_operator.py @@ -36,14 +36,15 @@ def _get_rtol(tp): class FFTOperatorTests(unittest.TestCase): @expand(product([16, ], [0.1, 1, 3.7], - [np.float64, np.float32, np.complex64, np.complex128])) - def test_fft1D(self, dim1, d, itp): + [np.float64, np.float32, np.complex64, np.complex128], + [ift.HartleyOperator, ift.FFTOperator])) + def test_fft1D(self, dim1, d, itp, op): tol = _get_rtol(itp) a = ift.RGSpace(dim1, distances=d) b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True) np.random.seed(16) - fft = ift.FFTOperator(domain=a, target=b) + fft = op(domain=a, target=b) inp = ift.Field.from_random(domain=a, random_type='normal', std=7, mean=3, dtype=itp) out = fft.inverse_times(fft.times(inp)) @@ -59,14 +60,15 @@ class FFTOperatorTests(unittest.TestCase): @expand(product([12, 15], [9, 12], [0.1, 1, 3.7], [0.4, 1, 2.7], - [np.float64, np.float32, np.complex64, np.complex128])) - def test_fft2D(self, dim1, dim2, d1, d2, itp): + [np.float64, np.float32, np.complex64, np.complex128], + [ift.HartleyOperator, ift.FFTOperator])) + def test_fft2D(self, dim1, dim2, d1, d2, itp, op): tol = _get_rtol(itp) a = ift.RGSpace([dim1, dim2], distances=[d1, d2]) b = ift.RGSpace([dim1, dim2], distances=[1./(dim1*d1), 1./(dim2*d2)], harmonic=True) - fft = ift.FFTOperator(domain=a, target=b) + fft = op(domain=a, target=b) inp = ift.Field.from_random(domain=a, random_type='normal', std=7, mean=3, dtype=itp) out = fft.inverse_times(fft.times(inp)) @@ -81,12 +83,13 @@ class FFTOperatorTests(unittest.TestCase): assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol) @expand(product([0, 1, 2], - [np.float64, np.float32, np.complex64, np.complex128])) - def test_composed_fft(self, index, dtype): + [np.float64, np.float32, np.complex64, np.complex128], + [ift.HartleyOperator, ift.FFTOperator])) + def test_composed_fft(self, index, dtype, op): tol = _get_rtol(dtype) a = [a1, a2, a3] = [ift.RGSpace((32,)), ift.RGSpace((4, 4)), ift.RGSpace((5, 6))] - fft = ift.FFTOperator(domain=a, space=index) + fft = op(domain=a, space=index) inp = ift.Field.from_random(domain=(a1, a2, a3), random_type='normal', std=7, mean=3, dtype=dtype) @@ -96,15 +99,16 @@ class FFTOperatorTests(unittest.TestCase): @expand(product([ift.RGSpace(128, distances=3.76, harmonic=True), ift.RGSpace((15, 27), distances=(.7, .33), harmonic=True), ift.RGSpace(73, distances=0.5643)], - [np.float64, np.float32, np.complex64, np.complex128])) - def test_normalisation(self, space, tp): + [np.float64, np.float32, np.complex64, np.complex128], + [ift.HartleyOperator, ift.FFTOperator])) + def test_normalisation(self, space, tp, op): tol = 10 * _get_rtol(tp) cospace = space.get_default_codomain() - fft = ift.FFTOperator(space, cospace) + fft = op(space, cospace) inp = ift.Field.from_random(domain=space, random_type='normal', std=1, mean=2, dtype=tp) out = fft.times(inp) - fft2 = ift.FFTOperator(cospace, space) + fft2 = op(cospace, space) out2 = fft2.inverse_times(inp) zero_idx = tuple([0]*len(space.shape)) assert_allclose(inp.to_global_data()[zero_idx], out.integrate(), diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py index a07c30088d73114c4922fff0b5f7bed318296183..e5e5d8e7910fcb39946cfd86305273633f821162 100644 --- a/test/test_operators/test_smoothing_operator.py +++ b/test/test_operators/test_smoothing_operator.py @@ -37,13 +37,13 @@ class SmoothingOperator_Tests(unittest.TestCase): @expand(product(spaces, [0., .5, 5.])) def test_property(self, space, sigma): - op = ift.FFTSmoothingOperator(space, sigma=sigma) + op = ift.HarmonicSmoothingOperator(space, sigma=sigma) if op.domain[0] != space: raise TypeError @expand(product(spaces, [0., .5, 5.])) def test_adjoint_times(self, space, sigma): - op = ift.FFTSmoothingOperator(space, sigma=sigma) + op = ift.HarmonicSmoothingOperator(space, sigma=sigma) rand1 = ift.Field.from_random('normal', domain=space) rand2 = ift.Field.from_random('normal', domain=space) tt1 = rand1.vdot(op.times(rand2)) @@ -52,7 +52,7 @@ class SmoothingOperator_Tests(unittest.TestCase): @expand(product(spaces, [0., .5, 5.])) def test_times(self, space, sigma): - op = ift.FFTSmoothingOperator(space, sigma=sigma) + op = ift.HarmonicSmoothingOperator(space, sigma=sigma) fld = np.zeros(space.shape, dtype=np.float64) fld[0] = 1. rand1 = ift.Field.from_global_data(space, fld) @@ -64,7 +64,7 @@ class SmoothingOperator_Tests(unittest.TestCase): def test_smooth_regular1(self, sz, d, sigma, tp): tol = _get_rtol(tp) sp = ift.RGSpace(sz, distances=d) - smo = ift.FFTSmoothingOperator(sp, sigma=sigma) + smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma) inp = ift.Field.from_random(domain=sp, random_type='normal', std=1, mean=4, dtype=tp) out = smo(inp) @@ -75,7 +75,7 @@ class SmoothingOperator_Tests(unittest.TestCase): def test_smooth_regular2(self, sz1, sz2, d1, d2, sigma, tp): tol = _get_rtol(tp) sp = ift.RGSpace([sz1, sz2], distances=[d1, d2]) - smo = ift.FFTSmoothingOperator(sp, sigma=sigma) + smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma) inp = ift.Field.from_random(domain=sp, random_type='normal', std=1, mean=4, dtype=tp) out = smo(inp)