diff --git a/demos/find_amplitude_parameters.py b/demos/find_amplitude_parameters.py new file mode 100644 index 0000000000000000000000000000000000000000..84a68944a1e35d567b96ddd99beb00be096a94e5 --- /dev/null +++ b/demos/find_amplitude_parameters.py @@ -0,0 +1,112 @@ +# 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 . +# +# Copyright(C) 2013-2019 Max-Planck-Society +# Author: Philipp Arras +# +# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. + +import numpy as np + +import nifty5 as ift +import matplotlib.pyplot as plt + + +def _default_pspace(dom): + return ift.PowerSpace(dom.get_default_codomain()) + + +if __name__ == '__main__': + np.random.seed(42) + fa = ift.CorrelatedFieldMaker() + n_samps = 20 + slope_means = [-2, -3] + fa.add_fluctuations(ift.RGSpace(128, 0.1), 10, 2, 1, 1e-6, + 2, 1e-6, slope_means[0], 0.2, 'spatial') + # fa.add_fluctuations(_default_pspace(ift.RGSpace((128, 64))), 10, 2, 1, + # 1e-6, 2, 1e-6, slope_means[0], 0.2, 'spatial') + fa.add_fluctuations(ift.RGSpace(32), 10, 5, 1, 1e-6, 2, + 1e-6, slope_means[1], 1, 'freq') + correlated_field = fa.finalize(10, 0.1, '') + amplitudes = fa.amplitudes + plt.style.use('seaborn-notebook') + + tgt = correlated_field.target + if len(tgt.shape) == 1: + fig, axes = plt.subplots(nrows=1, ncols=2) + fig.set_size_inches(20, 10) + else: + fig, axes = plt.subplots(nrows=3, ncols=3) + fig.set_size_inches(20, 16) + axs = (ax for ax in axes.ravel()) + for ii, aa in enumerate(amplitudes): + ax = next(axs) + pspec = aa**2 + ax.set_xscale('log') + ax.set_yscale('log') + for _ in range(n_samps): + fld = pspec(ift.from_random('normal', pspec.domain)) + klengths = fld.domain[0].k_lengths + ycoord = fld.to_global_data_rw() + ycoord[0] = ycoord[1] + ax.plot(klengths, ycoord, alpha=1) + + ymin, ymax = ax.get_ylim() + color = plt.rcParams['axes.prop_cycle'].by_key()['color'][0] + lbl = 'Mean slope (k^{})'.format(2*slope_means[ii]) + for fac in np.linspace(np.log(ymin), np.log(ymax**2/ymin)): + xs = np.linspace(np.amin(klengths[1:]), np.amax(klengths[1:])) + ys = xs**(2*slope_means[ii])*np.exp(fac) + xs = np.insert(xs, 0, 0) + ys = np.insert(ys, 0, ys[0]) + ax.plot(xs, ys, zorder=1, color=color, linewidth=0.3, label=lbl) + lbl = None + + ax.set_ylim([ymin, ymax]) + ax.set_xlim([None, np.amax(klengths)]) + ax.legend() + + if len(tgt.shape) == 2: + foo = [] + for ax in axs: + pos = ift.from_random('normal', correlated_field.domain) + fld = correlated_field(pos).to_global_data() + foo.append((ax, fld)) + mi, ma = np.inf, -np.inf + for _, fld in foo: + mi = min([mi, np.amin(fld)]) + ma = max([ma, np.amax(fld)]) + nxdx, nydy = tgt.shape + if len(tgt) == 2: + nxdx *= tgt[0].distances[0] + nydy *= tgt[1].distances[0] + else: + nxdx *= tgt[0].distances[0] + nydy *= tgt[0].distances[1] + for ax, fld in foo: + im = ax.imshow(fld.T, + extent=[0, nxdx, 0, nydy], + aspect='auto', + origin='lower', + vmin=mi, + vmax=ma) + fig.colorbar(im, ax=axes.ravel().tolist()) + elif len(tgt.shape) == 1: + ax = next(axs) + flds = [] + for _ in range(n_samps): + pos = ift.from_random('normal', correlated_field.domain) + ax.plot(correlated_field(pos).to_global_data()) + + plt.savefig('correlated_fields.png') + plt.close() diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py index 39db9675997ee87be6962278149abee2adb35bca..528fc19ac0cc9f3d98846a785fc58b7819d0898f 100644 --- a/demos/getting_started_3.py +++ b/demos/getting_started_3.py @@ -56,10 +56,10 @@ if __name__ == '__main__': filename = "getting_started_3_mode_{}_".format(mode) + "{}.png" position_space = ift.RGSpace([128, 128]) - power_space = ift.PowerSpace(position_space.get_default_codomain()) cfmaker = ift.CorrelatedFieldMaker() - cfmaker.add_fluctuations(power_space, 1, 1e-2, 1, .5, .1, .5, -3, 0.5, '') + cfmaker.add_fluctuations(position_space, + 1, 1e-2, 1, .5, .1, .5, -3, 0.5, '') correlated_field = cfmaker.finalize(1e-3, 1e-6, '') A = cfmaker.amplitudes[0] @@ -88,8 +88,8 @@ if __name__ == '__main__': minimizer = ift.NewtonCG(ic_newton) # Set up likelihood and information Hamiltonian - likelihood = ift.GaussianEnergy(mean=data, - inverse_covariance=N.inverse)(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.) diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py new file mode 100644 index 0000000000000000000000000000000000000000..d49872283d183f56d64f4f37356c0bd86d9dbcb6 --- /dev/null +++ b/demos/getting_started_mf.py @@ -0,0 +1,184 @@ +# 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 . +# +# Copyright(C) 2013-2019 Max-Planck-Society +# +# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. + +############################################################ +# Non-linear tomography +# +# The signal is a sigmoid-normal distributed field. +# The data is the field integrated along lines of sight that are +# randomly (set mode=0) or radially (mode=1) distributed +# +# Demo takes a while to compute +############################################################# + +import sys + +import numpy as np + +import nifty5 as ift + + +class SingleDomain(ift.LinearOperator): + def __init__(self, domain, target): + self._domain = ift.makeDomain(domain) + self._target = ift.makeDomain(target) + self._capability = self.TIMES | self.ADJOINT_TIMES + + def apply(self, x, mode): + self._check_input(x, mode) + return ift.from_global_data(self._tgt(mode), x.to_global_data()) + + +def random_los(n_los): + starts = list(np.random.uniform(0, 1, (n_los, 2)).T) + ends = list(np.random.uniform(0, 1, (n_los, 2)).T) + return starts, ends + + +def radial_los(n_los): + starts = list(np.random.uniform(0, 1, (n_los, 2)).T) + ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T) + return starts, ends + + +if __name__ == '__main__': + np.random.seed(45) + + # Choose between random line-of-sight response (mode=0) and radial lines + # of sight (mode=1) + if len(sys.argv) == 2: + mode = int(sys.argv[1]) + else: + mode = 0 + filename = "getting_started_mf_mode_{}_".format(mode) + "{}.png" + + npix1, npix2 = 128, 128 + position_space = ift.RGSpace([npix1, npix2]) + sp1 = ift.RGSpace(npix1) + sp2 = ift.RGSpace(npix2) + + cfmaker = ift.CorrelatedFieldMaker() + amp1 = 0.5 + cfmaker.add_fluctuations(sp1, amp1, 1e-2, 1, .1, .01, .5, -2, 1., 'amp1') + cfmaker.add_fluctuations(sp2, np.sqrt(1. - amp1**2), 1e-2, 1, .1, .01, .5, + -1.5, .5, 'amp2') + correlated_field = cfmaker.finalize(1e-3, 1e-6, '') + + A1 = cfmaker.amplitudes[0] + A2 = cfmaker.amplitudes[1] + DC = SingleDomain(correlated_field.target, position_space) + + # Apply a nonlinearity + signal = DC @ ift.sigmoid(correlated_field) + + # Build the line-of-sight response and define signal response + LOS_starts, LOS_ends = random_los(100) if mode == 0 else radial_los(100) + R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends) + signal_response = R(signal) + + # Specify noise + data_space = R.target + noise = .001 + N = ift.ScalingOperator(noise, data_space) + + # Generate mock signal and data + mock_position = ift.from_random('normal', signal_response.domain) + data = signal_response(mock_position) + N.draw_sample() + + # Minimization parameters + ic_sampling = ift.AbsDeltaEnergyController(name='Sampling', + deltaE=0.01, + iteration_limit=100) + ic_newton = ift.AbsDeltaEnergyController(name='Newton', + deltaE=0.01, + iteration_limit=35) + minimizer = ift.NewtonCG(ic_newton) + + # Set up likelihood and information Hamiltonian + 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.) + mean = initial_mean + + plot = ift.Plot() + plot.add(signal(mock_position), title='Ground Truth') + plot.add(R.adjoint_times(data), title='Data') + plot.add([A1.force(mock_position)], title='Power Spectrum 1') + plot.add([A2.force(mock_position)], title='Power Spectrum 2') + plot.output(ny=2, nx=2, xsize=10, ysize=10, name=filename.format("setup")) + + # number of samples used to estimate the KL + N_samples = 20 + + # Draw new samples to approximate the KL five times + for i in range(10): + # Draw new samples and minimize KL + KL = ift.MetricGaussianKL(mean, H, N_samples) + KL, convergence = minimizer(KL) + mean = KL.position + + # Plot current reconstruction + plot = ift.Plot() + plot.add(signal(mock_position), title="ground truth") + plot.add(signal(KL.position), title="reconstruction") + plot.add([A1.force(KL.position), + A1.force(mock_position)], + title="power1") + plot.add([A2.force(KL.position), + A2.force(mock_position)], + title="power2") + plot.output(nx=2, + ny=2, + ysize=10, + xsize=10, + name=filename.format("loop_{:02d}".format(i))) + + # Draw posterior samples + Nsamples = 20 + KL = ift.MetricGaussianKL(mean, H, N_samples) + sc = ift.StatCalculator() + scA1 = ift.StatCalculator() + scA2 = ift.StatCalculator() + powers1 = [] + powers2 = [] + for sample in KL.samples: + sc.add(signal(sample + KL.position)) + p1 = A1.force(sample + KL.position) + p2 = A2.force(sample + KL.position) + scA1.add(p1) + powers1.append(p1) + scA2.add(p2) + powers2.append(p2) + + # Plotting + filename_res = filename.format("results") + plot = ift.Plot() + plot.add(sc.mean, title="Posterior Mean") + plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation") + + powers1 = [A1.force(s + KL.position) for s in KL.samples] + powers2 = [A2.force(s + KL.position) for s in KL.samples] + plot.add(powers1 + [scA1.mean, A1.force(mock_position)], + title="Sampled Posterior Power Spectrum 1", + linewidth=[1.]*len(powers1) + [3., 3.]) + plot.add(powers2 + [scA2.mean, A2.force(mock_position)], + title="Sampled Posterior Power Spectrum 2", + linewidth=[1.]*len(powers2) + [3., 3.]) + plot.output(ny=2, nx=2, xsize=15, ysize=15, name=filename_res) + print("Saved results as '{}'.".format(filename_res)) diff --git a/demos/multi_amplitudes_consistency.py b/demos/multi_amplitudes_consistency.py new file mode 100644 index 0000000000000000000000000000000000000000..e7b4450819a4e55099aedd6ce183e959ad470a4c --- /dev/null +++ b/demos/multi_amplitudes_consistency.py @@ -0,0 +1,108 @@ +import nifty5 as ift +import numpy as np + + + + +def testAmplitudesConsistency(seed, sspace): + def stats(op,samples): + sc = ift.StatCalculator() + for s in samples: + sc.add(op(s.extract(op.domain))) + return sc.mean.to_global_data(), sc.var.sqrt().to_global_data() + np.random.seed(seed) + offset_std = 30 + intergated_fluct_std0 = .003 + intergated_fluct_std1 = 0.1 + + nsam = 1000 + + + fsspace = ift.RGSpace((12,), (0.4,)) + + fa = ift.CorrelatedFieldMaker() + fa.add_fluctuations(sspace, intergated_fluct_std0, 1E-8, 1.1, 2., 2.1, .5, + -2, 1., 'spatial') + fa.add_fluctuations(fsspace, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1, + -4, 1., 'freq') + op = fa.finalize(offset_std, 1E-8, '') + + samples = [ift.from_random('normal',op.domain) for _ in range(nsam)] + tot_flm, _ = stats(fa.total_fluctuation,samples) + offset_std,_ = stats(fa.amplitude_total_offset,samples) + intergated_fluct_std0,_ = stats(fa.average_fluctuation(0),samples) + intergated_fluct_std1,_ = stats(fa.average_fluctuation(1),samples) + + slice_fluct_std0,_ = stats(fa.slice_fluctuation(0),samples) + slice_fluct_std1,_ = stats(fa.slice_fluctuation(1),samples) + + sams = [op(s) for s in samples] + fluct_total = fa.total_fluctuation_realized(sams) + fluct_space = fa.average_fluctuation_realized(sams,0) + fluct_freq = fa.average_fluctuation_realized(sams,1) + zm_std_mean = fa.offset_amplitude_realized(sams) + sl_fluct_space = fa.slice_fluctuation_realized(sams,0) + sl_fluct_freq = fa.slice_fluctuation_realized(sams,1) + + np.testing.assert_allclose(offset_std, zm_std_mean, rtol=0.5) + np.testing.assert_allclose(intergated_fluct_std0, fluct_space, rtol=0.5) + np.testing.assert_allclose(intergated_fluct_std1, fluct_freq, rtol=0.5) + np.testing.assert_allclose(tot_flm, fluct_total, rtol=0.5) + np.testing.assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5) + np.testing.assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5) + + print("Expected offset Std: " + str(offset_std)) + print("Estimated offset Std: " + str(zm_std_mean)) + + print("Expected integrated fluct. space Std: " + + str(intergated_fluct_std0)) + print("Estimated integrated fluct. space Std: " + str(fluct_space)) + + print("Expected integrated fluct. frequency Std: " + + str(intergated_fluct_std1)) + print("Estimated integrated fluct. frequency Std: " + str(fluct_freq)) + + print("Expected slice fluct. space Std: " + + str(slice_fluct_std0)) + print("Estimated slice fluct. space Std: " + str(sl_fluct_space)) + + print("Expected slice fluct. frequency Std: " + + str(slice_fluct_std1)) + print("Estimated slice fluct. frequency Std: " + str(sl_fluct_freq)) + + print("Expected total fluct. Std: " + str(tot_flm)) + print("Estimated total fluct. Std: " + str(fluct_total)) + + fa = ift.CorrelatedFieldMaker() + fa.add_fluctuations(fsspace, intergated_fluct_std1, 1., 3.1, 1., .5, .1, + -4, 1., 'freq') + m = 3. + x = fa.moment_slice_to_average(m) + fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, + -2, 1., 'spatial', 0) + op = fa.finalize(offset_std, .1, '') + em, estd = stats(fa.slice_fluctuation(0),samples) + print("Forced slice fluct. space Std: "+str(m)) + print("Expected slice fluct. Std: " + str(em)) + np.testing.assert_allclose(m, em, rtol=0.5) + + assert op.target[0] == sspace + assert op.target[1] == fsspace + + +# Move to tests +# FIXME This test fails but it is not relevant for the final result +# assert_allclose(ampl(from_random('normal', ampl.domain)).val[0], vol) or 1?? +# End move to tests + +# move to tests +# assert_allclose( +# smooth(from_random('normal', smooth.domain)).val[0:2], 0) +# end move to tests +for seed in [1, 42]: + for sp in [ + ift.RGSpace((32, 64), (1.1, 0.3)), + ift.HPSpace(32), + ift.GLSpace(32) + ]: + testAmplitudesConsistency(seed, sp) diff --git a/demos/newamplitudes.py b/demos/newamplitudes.py index 5a8e437e91e7ed183f25be4c3007afdd306e6f6d..507cdbb5934ff49c4f6cb32d5fb4f9480eabef75 100644 --- a/demos/newamplitudes.py +++ b/demos/newamplitudes.py @@ -3,11 +3,9 @@ import numpy as np np.random.seed(42) sspace = ift.RGSpace((128,)) -hspace = sspace.get_default_codomain() -target0 = ift.PowerSpace(hspace) fa = ift.CorrelatedFieldMaker() -fa.add_fluctuations(target0, 10, 2, 1, 1e-6, 2, 1e-6, -2, 1e-6, 'spatial') +fa.add_fluctuations(sspace, 10, 2, 1, 1e-6, 2, 1e-6, -2, 1e-6, 'spatial') op = fa.finalize(10, 0.1, '') A = fa.amplitudes[0] diff --git a/docs/source/volume.rst b/docs/source/volume.rst index e12b8b1980ed602ac764d3111a195e9a9d5c8778..89dc9419d28bd59d79bfe0ddd0fa4ecee3a786a1 100644 --- a/docs/source/volume.rst +++ b/docs/source/volume.rst @@ -211,3 +211,18 @@ Consequently, the inverse covariance operator will automatically have lower indi ensuring that the whole log-likelihood expression does not scale with resolution. **This upper-lower index convention is not coded into NIFTy**, in order to not reduce user freedom. One should however have this in mind when constructing log-likelihoods in order to ensure resolution independence. + +Harmonic Transform Convention +............................. + +In NIFTy the convention for the harmonic transformations is set by requiring the zero mode of the transformed field to be the integral over the original field. +This choice is convenient for the Fourier transformation and used throughout the literature. +Note that for the spherical harmonics this convention is only rarely used and might result in unexpected factors in the transformed field. + +To be specific, for the spherical harmonics transformation this means that a monopole of unit amplitude in position-space which is transformed to the spherical harmonics domain and back to the original position space via the adjoint transformation will have a non-unit amplitude afterwards. +The factor between the input field and the twice transformed field is :math:`1/4\pi`. +In comparison to the convention used in HEALPix, this corresponds to dividing the output of the HEALPix transformed field by :math:`\sqrt{4\pi}` in each transformation. + +Depending on the use-case, additional volume factors must be accounted for. +This is especially true if one wants to define the inverse transformation. +Note that the convention for the harmonic transformations used in NIFTy results in non-unitary operators. diff --git a/nifty5/__init__.py b/nifty5/__init__.py index 158419afd43948730d7c639ebdfe2c2f9d60bc88..a5d1953f22d04304aac3d0f9ed427199149d0ce3 100644 --- a/nifty5/__init__.py +++ b/nifty5/__init__.py @@ -11,7 +11,6 @@ from .domains.gl_space import GLSpace from .domains.hp_space import HPSpace from .domains.power_space import PowerSpace from .domains.dof_space import DOFSpace -from .domains.log_rg_space import LogRGSpace from .domain_tuple import DomainTuple from .multi_domain import MultiDomain @@ -20,13 +19,13 @@ from .multi_field import MultiField from .operators.operator import Operator from .operators.adder import Adder +from .operators.log1p import Log1p from .operators.diagonal_operator import DiagonalOperator from .operators.distributors import DOFDistributor, PowerDistributor from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter from .operators.contraction_operator import ContractionOperator from .operators.linear_interpolation import LinearInterpolator from .operators.endomorphic_operator import EndomorphicOperator -from .operators.exp_transform import ExpTransform from .operators.harmonic_operators import ( FFTOperator, HartleyOperator, SHTOperator, HarmonicTransformOperator, HarmonicSmoothingOperator) @@ -34,13 +33,10 @@ from .operators.field_zero_padder import FieldZeroPadder from .operators.inversion_enabler import InversionEnabler from .operators.linear_operator import LinearOperator from .operators.mask_operator import MaskOperator -from .operators.qht_operator import QHTOperator from .operators.regridding_operator import RegriddingOperator from .operators.sampling_enabler import SamplingEnabler from .operators.sandwich_operator import SandwichOperator 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.outer_product_operator import OuterProduct from .operators.simple_linear_operators import ( @@ -51,7 +47,7 @@ from .operators.value_inserter import ValueInserter from .operators.energy_operators import ( EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood, BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator, - Squared2NormOperator) + Squared2NormOperator, StudentTEnergy) from .operators.convolution_operators import FuncConvolutionOperator from .probing import probe_with_posterior_samples, probe_diagonal, \ diff --git a/nifty5/data_objects/distributed_do.py b/nifty5/data_objects/distributed_do.py index 765b77c4fdd615379d8ede0ea448b6c0a18e1b07..004e2949c69c952fb379f05a91bb087a6d912f22 100644 --- a/nifty5/data_objects/distributed_do.py +++ b/nifty5/data_objects/distributed_do.py @@ -32,7 +32,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", "redistribute", "default_distaxis", "is_numpy", "absmax", "norm", "lock", "locked", "uniform_full", "transpose", "to_global_data_rw", "ensure_not_distributed", "ensure_default_distributed", - "tanh", "conjugate", "sin", "cos", "tan", + "tanh", "conjugate", "sin", "cos", "tan", "log10", "sinh", "cosh", "sinc", "absolute", "sign", "clip"] _comm = MPI.COMM_WORLD @@ -297,7 +297,7 @@ def _math_helper(x, function, out): _current_module = sys.modules[__name__] for f in ["sqrt", "exp", "log", "tanh", "conjugate", "sin", "cos", "tan", - "sinh", "cosh", "sinc", "absolute", "sign"]: + "sinh", "cosh", "sinc", "absolute", "sign", "log10"]: def func(f): def func2(x, out=None): return _math_helper(x, f, out) diff --git a/nifty5/data_objects/numpy_do.py b/nifty5/data_objects/numpy_do.py index ed0d96de8f3e855310bf1e5516909b1cb05b7913..02bd03fe819c76751e1495ac6364393e4c263528 100644 --- a/nifty5/data_objects/numpy_do.py +++ b/nifty5/data_objects/numpy_do.py @@ -18,9 +18,11 @@ # Data object module that uses simple numpy ndarrays. import numpy as np -from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log from numpy import ndarray as data_object -from numpy import ones, sign, sin, sinc, sinh, sqrt, tan, tanh, vdot, zeros +from numpy import empty, empty_like, ones, zeros, full +from numpy import absolute, sign, clip, vdot +from numpy import sin, cos, sinh, cosh, tan, tanh +from numpy import exp, log, log10, sqrt, sinc from .random import Random @@ -34,7 +36,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", "lock", "locked", "uniform_full", "to_global_data_rw", "ensure_not_distributed", "ensure_default_distributed", "clip", "sin", "cos", "tan", "sinh", - "cosh", "absolute", "sign", "sinc"] + "cosh", "absolute", "sign", "sinc", "log10"] ntask = 1 rank = 0 diff --git a/nifty5/domain_tuple.py b/nifty5/domain_tuple.py index e6e56da10d1eda8bf01577ce0e46f3944df6a780..9200b85f33851d44ae3a890f7292b2af52f1fda2 100644 --- a/nifty5/domain_tuple.py +++ b/nifty5/domain_tuple.py @@ -19,6 +19,8 @@ from functools import reduce from . import utilities from .domains.domain import Domain +import numpy as np + class DomainTuple(object): """Ordered sequence of Domain objects. @@ -125,6 +127,58 @@ class DomainTuple(object): """ return self._size + def scalar_weight(self, spaces=None): + """Returns the uniform volume element for a sub-domain of `self`. + + Parameters + ---------- + spaces : int, tuple of int or None + Indices of the sub-domains to be considered. + If `None`, the entire domain is used. + + Returns + ------- + float or None + If the requested sub-domain has a uniform volume element, it is + returned. Otherwise, `None` is returned. + """ + if np.isscalar(spaces): + return self._dom[spaces].scalar_dvol + + if spaces is None: + spaces = range(len(self._dom)) + res = 1. + for i in spaces: + tmp = self._dom[i].scalar_dvol + if tmp is None: + return None + res *= tmp + return res + + def total_volume(self, spaces=None): + """Returns the total volume of `self` or of a subspace of it. + + Parameters + ---------- + spaces : int, tuple of int or None + Indices of the sub-domains of the domain to be considered. + If `None`, the total volume of the whole domain is returned. + + Returns + ------- + float + the total volume of the requested (sub-)domain. + """ + if np.isscalar(spaces): + return self._dom[spaces].total_volume + + if spaces is None: + spaces = range(len(self._dom)) + res = 1. + for i in spaces: + res *= self._dom[i].total_volume + return res + @property def axes(self): """tuple of tuple of int : axis indices of the underlying domains""" diff --git a/nifty5/domains/log_rg_space.py b/nifty5/domains/log_rg_space.py deleted file mode 100644 index 8b9c6095c0a503206bd80caa62aa9ff8f59b7f2f..0000000000000000000000000000000000000000 --- a/nifty5/domains/log_rg_space.py +++ /dev/null @@ -1,146 +0,0 @@ -# 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 . -# -# 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 diff --git a/nifty5/extra.py b/nifty5/extra.py index 97146e3778e9613e27a5148d52fcbfc0f2ee04ca..c39ad98b4c7fc092307125bf1717e9514e197d09 100644 --- a/nifty5/extra.py +++ b/nifty5/extra.py @@ -16,11 +16,13 @@ # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. import numpy as np +from numpy.testing import assert_ from .domain_tuple import DomainTuple from .field import Field from .linearization import Linearization from .multi_domain import MultiDomain +from .multi_field import MultiField from .operators.linear_operator import LinearOperator from .sugar import from_random @@ -81,6 +83,38 @@ def _check_linearity(op, domain_dtype, atol, rtol): _assert_allclose(val1, val2, atol=atol, rtol=rtol) +def _actual_domain_check(op, domain_dtype=None, inp=None): + needed_cap = op.TIMES + if (op.capability & needed_cap) != needed_cap: + return + if domain_dtype is not None: + inp = from_random("normal", op.domain, dtype=domain_dtype) + elif inp is None: + raise ValueError('Need to specify either dtype or inp') + assert_(inp.domain is op.domain) + assert_(op(inp).domain is op.target) + + +def _actual_domain_check_nonlinear(op, loc): + assert isinstance(loc, (Field, MultiField)) + assert_(loc.domain is op.domain) + lin = Linearization.make_var(loc, False) + reslin = op(lin) + assert_(lin.domain is op.domain) + assert_(lin.target is op.domain) + assert_(lin.val.domain is lin.domain) + + assert_(reslin.domain is op.domain) + assert_(reslin.target is op.target) + assert_(reslin.val.domain is reslin.target) + + assert_(reslin.target is op.target) + assert_(reslin.jac.domain is reslin.domain) + assert_(reslin.jac.target is reslin.target) + _actual_domain_check(reslin.jac, inp=loc) + _actual_domain_check(reslin.jac.adjoint, inp=reslin.jac(loc)) + + def _domain_check(op): for dd in [op.domain, op.target]: if not isinstance(dd, (DomainTuple, MultiDomain)): @@ -123,6 +157,10 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64, if not isinstance(op, LinearOperator): raise TypeError('This test tests only linear operators.') _domain_check(op) + _actual_domain_check(op, domain_dtype) + _actual_domain_check(op.adjoint, target_dtype) + _actual_domain_check(op.inverse, target_dtype) + _actual_domain_check(op.adjoint.inverse, domain_dtype) _check_linearity(op, domain_dtype, atol, rtol) _check_linearity(op.adjoint, target_dtype, atol, rtol) _check_linearity(op.inverse, target_dtype, atol, rtol) @@ -180,6 +218,7 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100): Tolerance for the check. """ _domain_check(op) + _actual_domain_check_nonlinear(op, loc) for _ in range(ntries): lin = op(Linearization.make_var(loc)) loc2, lin2 = _get_acceptable_location(op, loc, lin) diff --git a/nifty5/field.py b/nifty5/field.py index 5ba5fba3ae79c920fd503ea4b819429e0f184580..685cb34e0f22de6e62f68eece829961b41949f52 100644 --- a/nifty5/field.py +++ b/nifty5/field.py @@ -246,42 +246,23 @@ class Field(object): If the requested sub-domain has a uniform volume element, it is returned. Otherwise, `None` is returned. """ - if np.isscalar(spaces): - return self._domain[spaces].scalar_dvol - - if spaces is None: - spaces = range(len(self._domain)) - res = 1. - for i in spaces: - tmp = self._domain[i].scalar_dvol - if tmp is None: - return None - res *= tmp - return res + return self._domain.scalar_weight(spaces) def total_volume(self, spaces=None): - """Returns the total volume of a sub-domain of `self`. + """Returns the total volume of the field's domain or of a subspace of it. Parameters ---------- spaces : int, tuple of int or None Indices of the sub-domains of the field's domain to be considered. - If `None`, the entire domain is used. + If `None`, the total volume of the whole domain is returned. Returns ------- float - the total volume of the requested sub-domain. + the total volume of the requested (sub-)domain. """ - if np.isscalar(spaces): - return self._domain[spaces].total_volume - - if spaces is None: - spaces = range(len(self._domain)) - res = 1. - for i in spaces: - res *= self._domain[i].total_volume - return res + return self._domain.total_volume(spaces) def weight(self, power=1, spaces=None): """Weights the pixels of `self` with their invidual pixel volumes. @@ -682,7 +663,7 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__", return func2 setattr(Field, op, func(op)) -for f in ["sqrt", "exp", "log", "tanh", +for f in ["sqrt", "exp", "log", "log10", "tanh", "sin", "cos", "tan", "cosh", "sinh", "absolute", "sinc", "sign"]: def func(f): diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py index b218c997a8d1b7286c122450e3ae79fe84971d60..6ea4853ec9960940caa835149209d2165da94d6c 100644 --- a/nifty5/library/correlated_fields.py +++ b/nifty5/library/correlated_fields.py @@ -17,15 +17,11 @@ # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. import numpy as np -from functools import reduce -from numpy.testing import assert_allclose from ..domain_tuple import DomainTuple from ..domains.power_space import PowerSpace from ..domains.unstructured_domain import UnstructuredDomain -from ..extra import check_jacobian_consistency, consistency_check from ..field import Field -from ..multi_domain import MultiDomain from ..operators.adder import Adder from ..operators.contraction_operator import ContractionOperator from ..operators.distributors import PowerDistributor @@ -36,11 +32,12 @@ from ..operators.diagonal_operator import DiagonalOperator from ..operators.operator import Operator from ..operators.simple_linear_operators import VdotOperator, ducktape from ..operators.value_inserter import ValueInserter -from ..sugar import from_global_data, from_random, full, makeDomain, get_default_codomain +from ..probing import StatCalculator +from ..sugar import from_global_data, full, makeDomain, get_default_codomain + def _reshaper(x, shape): x = np.array(x) - if x.shape == shape: return np.asfarray(x) elif x.shape in [(), (1,)]: @@ -48,31 +45,89 @@ def _reshaper(x, shape): else: raise TypeError("Shape of parameters cannot be interpreted") -def _lognormal_moment_matching(mean, sig, key, - domain = DomainTuple.scalar_domain(), space = 0): - domain = makeDomain(domain) - mean, sig = (_reshaper(param, domain.shape) for param in (mean, sig)) - key = str(key) - assert np.all(mean > 0) + +def _lognormal_moments(mean, sig, shape = ()): + mean, sig = (_reshaper(param, shape) for param in (mean, sig)) + assert np.all(mean > 0 ) assert np.all(sig > 0) logsig = np.sqrt(np.log((sig/mean)**2 + 1)) logmean = np.log(mean) - logsig**2/2 - return _normal(logmean, logsig, key, domain).exp() + return logmean, logsig -def _normal(mean, sig, key, - domain = DomainTuple.scalar_domain(), space = 0): +def _normal(mean, sig, key, domain = DomainTuple.scalar_domain()): domain = makeDomain(domain) mean, sig = (_reshaper(param, domain.shape) for param in (mean, sig)) assert np.all(sig > 0) return Adder(from_global_data(domain, mean)) @ ( - DiagonalOperator(from_global_data(domain,sig)) @ ducktape(domain, None, key)) + DiagonalOperator(from_global_data(domain,sig)) + @ ducktape(domain, None, key)) + + +def _log_k_lengths(pspace): + """Log(k_lengths) without zeromode""" + return np.log(pspace.k_lengths[1:]) + + +def _relative_log_k_lengths(power_space): + """Log-distance to first bin + logkl.shape==power_space.shape, logkl[0]=logkl[1]=0""" + power_space = DomainTuple.make(power_space) + assert isinstance(power_space[0], PowerSpace) + assert len(power_space) == 1 + logkl = _log_k_lengths(power_space[0]) + assert logkl.shape[0] == power_space[0].shape[0] - 1 + logkl -= logkl[0] + return np.insert(logkl, 0, 0) + + +def _log_vol(power_space): + power_space = makeDomain(power_space) + assert isinstance(power_space[0], PowerSpace) + logk_lengths = _log_k_lengths(power_space[0]) + return logk_lengths[1:] - logk_lengths[:-1] + + +def _total_fluctuation_realized(samples, space = 0): + res = 0. + for s in samples: + res = res + (s - s.mean(space, keepdims = True))**2 + return np.sqrt((res/len(samples)).mean(space)) + + +def _stats(op, samples): + sc = StatCalculator() + for s in samples: + sc.add(op(s.extract(op.domain))) + return sc.mean.to_global_data(), sc.var.sqrt().to_global_data() + + +class _LognormalMomentMatching(Operator): + def __init__(self, mean, sig, key, + domain = DomainTuple.scalar_domain()): + key = str(key) + logmean, logsig = _lognormal_moments(mean, sig, domain.shape) + self._mean = mean + self._sig = sig + op = _normal(logmean, logsig, key, domain).exp() + self._domain, self._target = op.domain, op.target + self.apply = op.apply + + @property + def mean(self): + return self._mean + + @property + def std(self): + return self._sig class _SlopeRemover(EndomorphicOperator): - def __init__(self, domain, cooridinates, space = 0): + def __init__(self, domain, space = 0): self._domain = makeDomain(domain) - self._sc = cooridinates / float(cooridinates[-1]) + assert isinstance(self._domain[space], PowerSpace) + logkl = _relative_log_k_lengths(self._domain[space]) + self._sc = logkl/float(logkl[-1]) self._space = space axis = self._domain.axes[space][0] @@ -80,46 +135,29 @@ class _SlopeRemover(EndomorphicOperator): self._extender = (None,)*(axis) + (slice(None),) + (None,)*(self._domain.axes[-1][-1]-axis) self._capability = self.TIMES | self.ADJOINT_TIMES - def apply(self,x,mode): - self._check_input(x,mode) + def apply(self, x, mode): + self._check_input(x, mode) x = x.to_global_data() if mode == self.TIMES: res = x - x[self._last]*self._sc[self._extender] else: - #NOTE Why not x.copy()? - res = np.zeros(x.shape,dtype=x.dtype) - res += x - res[self._last] -= np.sum(x*self._sc[self._extender], axis = self._space, keepdims = True) - return from_global_data(self._tgt(mode),res) - -def _make_slope_Operator(smooth,loglogavgslope, space = 0): - tg = smooth.target - logkl = _log_k_lengths(tg[space]) - logkl -= logkl[0] - logkl = np.insert(logkl, 0, 0) - noslope = smooth - noslope = _SlopeRemover(tg,logkl, space) @ smooth - # FIXME Move to tests - consistency_check(_SlopeRemover(tg,logkl, space)) + res = x.copy() + res[self._last] -= (x*self._sc[self._extender]).sum( + axis = self._space, keepdims = True) + return from_global_data(self._tgt(mode), res) - expander = ContractionOperator(tg, spaces = space).adjoint - _t = DiagonalOperator(from_global_data(tg[space], logkl), tg, spaces = space) - return _t @ expander @ loglogavgslope + noslope - -def _log_k_lengths(pspace): - return np.log(pspace.k_lengths[1:]) class _TwoLogIntegrations(LinearOperator): - def __init__(self, target, space = None): + def __init__(self, target, space = 0): self._target = makeDomain(target) assert isinstance(self.target[space], PowerSpace) dom = list(self._target) dom[space] = UnstructuredDomain((2, self.target[space].shape[0]-2)) self._domain = makeDomain(dom) self._space = space - self._capability = self.TIMES | self.ADJOINT_TIMES logk_lengths = _log_k_lengths(self._target[space]) - self._logvol = logk_lengths[1:] - logk_lengths[:-1] + self._log_vol = _log_vol(self._target[space]) + self._capability = self.TIMES | self.ADJOINT_TIMES def apply(self, x, mode): self._check_input(x, mode) @@ -133,20 +171,19 @@ class _TwoLogIntegrations(LinearOperator): from_third = sl + (slice(2,None),) no_border = sl + (slice(1,-1),) reverse = sl + (slice(None,None,-1),) + + x = x.to_global_data_rw() if mode == self.TIMES: - x = x.to_global_data() res = np.empty(self._target.shape) - res[first] = 0 - res[second] = 0 + res[first] = res[second] = 0 res[from_third] = np.cumsum(x[second], axis = axis) - res[from_third] = (res[from_third] + res[no_border])/2*self._logvol[extender_sl] + x[first] + res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first] res[from_third] = np.cumsum(res[from_third], axis = axis) else: - x = x.to_global_data_rw() res = np.zeros(self._domain.shape) x[from_third] = np.cumsum(x[from_third][reverse], axis = axis)[reverse] res[first] += x[from_third] - x[from_third] *= (self._logvol/2.)[extender_sl] + x[from_third] *= (self._log_vol/2.)[extender_sl] x[no_border] += x[from_third] res[second] += np.cumsum(x[from_third][reverse], axis = axis)[reverse] return from_global_data(self._tgt(mode), res) @@ -155,18 +192,17 @@ class _TwoLogIntegrations(LinearOperator): class _Normalization(Operator): def __init__(self, domain, space = 0): self._domain = self._target = makeDomain(domain) + assert isinstance(self._domain[space], PowerSpace) hspace = list(self._domain) hspace[space] = hspace[space].harmonic_partner hspace = makeDomain(hspace) pd = PowerDistributor(hspace, power_space=self._domain[space], space = space) # TODO Does not work on sphere yet mode_multiplicity = pd.adjoint(full(pd.target, 1.)).to_global_data_rw() - zero_mode = (slice(None),)*domain.axes[space][0] + (0,) + zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,) mode_multiplicity[zero_mode] = 0 self._mode_multiplicity = from_global_data(self._domain, mode_multiplicity) self._specsum = _SpecialSum(self._domain, space) - # FIXME Move to tests - consistency_check(self._specsum) def apply(self, x): self._check_input(x) @@ -208,13 +244,9 @@ class _slice_extractor(LinearOperator): return from_global_data(self._tgt(mode), res) -class CorrelatedFieldMaker: - def __init__(self): - self._amplitudes = [] - self._spaces = [] - - def add_fluctuations_from_ops(self, target, fluctuations, flexibility, - asperity, loglogavgslope, key, space = 0): +class _Amplitude(Operator): + def __init__(self, target, fluctuations, flexibility, asperity, + loglogavgslope, key, space = 0): """ fluctuations > 0 flexibility > 0 @@ -228,154 +260,305 @@ class CorrelatedFieldMaker: target = makeDomain(target) assert isinstance(target[space], PowerSpace) + target = makeDomain(target) twolog = _TwoLogIntegrations(target, space) - dt = twolog._logvol - axis = target.axes[space][0] - first = (slice(None),)*axis + (0,) - extender_sl = (None,)*axis + (slice(None),) + (None,)*(target.axes[-1][-1] - axis) - expander = ContractionOperator(twolog.domain, spaces = space).adjoint - - sqrt_t = np.zeros(twolog.domain[space].shape) - sqrt_t[0] = sqrt_t[1] = np.sqrt(dt) - sqrt_t = from_global_data(twolog.domain[space], sqrt_t) - sqrt_t = DiagonalOperator(sqrt_t, twolog.domain, spaces = space) - sigmasq = sqrt_t @ expander @ flexibility - - dist = np.zeros(twolog.domain[space].shape) - dist[0] += 1. - dist = from_global_data(twolog.domain[space], dist) - dist = DiagonalOperator(dist, twolog.domain, spaces = space) - - shift = np.ones(twolog.domain.shape) - shift[first] = (dt**2/12.)[extender_sl] - shift = from_global_data(twolog.domain, shift) - scale = sigmasq*(Adder(shift) @ dist @ expander @ asperity).sqrt() - - smooth = twolog @ (scale*ducktape(scale.target, None, key)) - smoothslope = _make_slope_Operator(smooth,loglogavgslope, space) + dom = twolog.domain + shp = dom[space].shape + totvol = target[space].harmonic_partner.get_default_codomain().total_volume + expander = ContractionOperator(dom, spaces = space).adjoint + ps_expander = ContractionOperator(twolog.target, spaces = space).adjoint + + # Prepare constant fields + foo = np.zeros(shp) + foo[0] = foo[1] = np.sqrt(_log_vol(target[space])) + vflex = DiagonalOperator(from_global_data(dom[space], foo), dom, space) + + foo = np.zeros(shp, dtype=np.float64) + foo[0] += 1 + vasp = DiagonalOperator(from_global_data(dom[space], foo), dom, space) + + foo = np.ones(shp) + foo[0] = _log_vol(target[space])**2/12. + shift = DiagonalOperator(from_global_data(dom[space], foo), dom, space) - # move to tests - assert_allclose( - smooth(from_random('normal', smooth.domain)).val[( - slice(None),)*axis + (slice(0,2),)], 0) - consistency_check(twolog) - check_jacobian_consistency(smooth, from_random('normal', - smooth.domain)) - check_jacobian_consistency(smoothslope, - from_random('normal', smoothslope.domain)) - # end move to tests - - normal_ampl = _Normalization(target, space) @ smoothslope - vol = target[space].harmonic_partner.get_default_codomain().total_volume - arr = np.zeros(target[space].shape) - arr[1:] = vol - expander = ContractionOperator(target, spaces = space).adjoint - expander = DiagonalOperator(from_global_data(target[space], arr) - , target, spaces = space) @ expander - mask = np.zeros(target.shape) - mask[first] = vol - adder = Adder(from_global_data(target, mask)) - ampl = adder @ ((expander @ fluctuations)*normal_ampl) - - # Move to tests - # FIXME This test fails but it is not relevant for the final result - # assert_allclose( - # normal_ampl(from_random('normal', normal_ampl.domain)).val[0], 1) - #assert_allclose(ampl(from_random('normal', ampl.domain)).val[space][0], vol) - op = _Normalization(target, space) - check_jacobian_consistency(op, from_random('normal', op.domain)) - # End move to tests - - self._amplitudes.append(ampl) - self._spaces.append(space) - - def add_fluctuations(self, target, fluctuations_mean, fluctuations_stddev, - flexibility_mean, flexibility_stddev, asperity_mean, - asperity_stddev, loglogavgslope_mean, - loglogavgslope_stddev, prefix, space = 0): - prefix = str(prefix) + vslope = DiagonalOperator( + from_global_data(target[space], _relative_log_k_lengths(target[space])), + target, space) + + foo, bar = [np.zeros(target[space].shape) for _ in range(2)] + bar[1:] = foo[0] = totvol + vol0, vol1 = [DiagonalOperator(from_global_data(target[space], aa), + target, space) for aa in (foo, bar)] + + #Prepare fields for Adder + #NOTE alternative would be adjoint contraction_operator acting + #on every space except the specified on + shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)] + # End prepare constant fields + + slope = vslope @ ps_expander @ loglogavgslope + sig_flex = vflex @ expander @ flexibility + sig_asp = vasp @ expander @ asperity + sig_fluc = vol1 @ ps_expander @ fluctuations + + xi = ducktape(dom, None, key) + sigma = sig_flex*(Adder(shift) @ sig_asp).sqrt() + smooth = _SlopeRemover(target, space) @ twolog @ (sigma*xi) + op = _Normalization(target, space) @ (slope + smooth) + op = Adder(vol0) @ (sig_fluc*op) + + self.apply = op.apply + self._fluc = fluctuations + self._domain, self._target = op.domain, op.target + self._space = space - parameter_domain = list(makeDomain(target)) - del parameter_domain[space] - if parameter_domain != []: - parameter_domain = makeDomain(parameter_domain) - else: - parameter_domain = DomainTuple.scalar_domain() - - fluct = _lognormal_moment_matching(fluctuations_mean, fluctuations_stddev, - prefix + 'fluctuations', parameter_domain, space = space) - flex = _lognormal_moment_matching(flexibility_mean, flexibility_stddev, - prefix + 'flexibility', parameter_domain, space = space) - asp = _lognormal_moment_matching(asperity_mean, asperity_stddev, - prefix + 'asperity', parameter_domain, space = space) + @property + def fluctuation_amplitude(self): + return self._fluc + + +class CorrelatedFieldMaker: + def __init__(self): + self._a = [] + self._spaces = [] + self._azm = None + self._position_spaces = [] + + def add_fluctuations(self, + position_space, + fluctuations_mean, + fluctuations_stddev, + flexibility_mean, + flexibility_stddev, + asperity_mean, + asperity_stddev, + loglogavgslope_mean, + loglogavgslope_stddev, + prefix='', + index=None, + space=0): + position_space = makeDomain(position_space) + power_space = list(position_space) + power_space[space] = PowerSpace(position_space[space].get_default_codomain()) + power_space = makeDomain(power_space) + prefix = str(prefix) + #assert isinstance(position_space[space], (RGSpace, HPSpace, GLSpace) + auxdom = makeDomain(tuple(dom for i, dom in enumerate(position_space) + if i != space)) + + fluct = _LognormalMomentMatching(fluctuations_mean, + fluctuations_stddev, + prefix + 'fluctuations', + auxdom) + flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev, + prefix + 'flexibility', + auxdom) + asp = _LognormalMomentMatching(asperity_mean, asperity_stddev, + prefix + 'asperity', + auxdom) avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev, - prefix + 'loglogavgslope', parameter_domain, space = space) + prefix + 'loglogavgslope', auxdom) + amp = _Amplitude(power_space, + fluct, flex, asp, avgsl, prefix + 'spectrum', space) + if index is not None: + self._a.insert(index, amp) + self._position_spaces.insert(index, position_space) + self._spaces.insert(index, space) + else: + self._a.append(amp) + self._position_spaces.append(position_space) + self._spaces.append(spaces) + + def finalize_from_op(self, zeromode, prefix='', space = 0): + assert isinstance(zeromode, Operator) + self._azm = zeromode + hspace = [] + tuple(hspace.extend(tuple(get_default_codomain(dd, space))) + for dd, space in zip(self._position_spaces, self._spaces)) + hspace = makeDomain(hspace) + zeroind = () + for i, dd in enumerate(self._position_spaces): + zeroind += (slice(None),)*self._spaces[i] + zeroind += (0,)*len(dd[self._spaces[i]].shape) + #tuple(zeroind = zeroind + (slice(None),)*space + len(dd.shape)*(0,) + foo = np.ones(hspace.shape) + foo[zeroind] = 0 - return self.add_fluctuations_from_ops(target, fluct, flex, asp, avgsl, - prefix + 'spectrum', space) + ZeroModeInserter = _slice_extractor(hspace, + self._azm.target, zeroind).adjoint + azm = Adder(from_global_data(hspace, foo)) @ ZeroModeInserter @ zeromode + + n_amplitudes = len(self._a) + ht = HarmonicTransformOperator(hspace, + self._position_spaces[0][self._spaces[0]], + space=self._spaces[0]) + for i in range(1, n_amplitudes): + ht = (HarmonicTransformOperator(ht.target, + self._position_spaces[i][self._spaces[i]], + space=self._spaces[i]) @ ht) + + pd = PowerDistributor(hspace, self._a[0].target[self._spaces[0]], self._spaces[0]) + for i in range(1, n_amplitudes): + pd = (PowerDistributor(pd.domain, + self._a[i].target[self._spaces[i]], + space=self._spaces[i]) @ pd) - def finalize_from_op(self, zeromode): - raise NotImplementedError + spaces = np.cumsum(self._spaces) + np.arange(len(self._spaces)) + a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ self._a[0] + for i in range(1, n_amplitudes): + co = ContractionOperator(pd.domain, spaces[:i] + spaces[(i + 1):]) + a = a*(co.adjoint @ self._a[i]) + + return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi')) def finalize(self, offset_amplitude_mean, offset_amplitude_stddev, - prefix, - offset=None): + prefix='', + offset=None, + prior_info=100): """ offset vs zeromode: volume factor """ prefix = str(prefix) if offset is not None: - offset = float(offset) - hspace = [] - zeroind = () - for amp, space in zip(self._amplitudes, self._spaces): - dd = list(amp.target) - dd[space] = dd[space].harmonic_partner - hspace.extend(dd) - zeroind += (slice(None),)*space + (0,)*len(dd[space].shape) - hspace = makeDomain(hspace) - spaces = np.cumsum(self._spaces) + np.arange(len(self._spaces)) - - parameter_domain = list(makeDomain(hspace)) - for space in self._spaces: - del parameter_domain[space] - if parameter_domain != []: - parameter_domain = makeDomain(parameter_domain) - else: - parameter_domain = DomainTuple.scalar_domain() + raise NotImplementedError + offset = float(offset) + + auxdom = makeDomain((dom for i, domT in enumerate(self._position_spaces) + for j, dom in enumerate(domT) + if j != self._spaces[i])) + azm = _LognormalMomentMatching(offset_amplitude_mean, + offset_amplitude_stddev, + prefix + 'zeromode', + auxdom) + op = self.finalize_from_op(azm, prefix) + if prior_info > 0: + from ..sugar import from_random + samps = [ + from_random('normal', op.domain) for _ in range(prior_info) + ] + self.statistics_summary(samps) + return op + + def statistics_summary(self, samples): + lst = [('Offset amplitude', self.amplitude_total_offset), + ('Total fluctuation amplitude', self.total_fluctuation)] + + namps = len(self.amplitudes) + if namps > 1: + for ii in range(namps): + lst.append(('Slice fluctuation (space {})'.format(ii), + self.slice_fluctuation(ii))) + lst.append(('Average fluctuation (space {})'.format(ii), + self.average_fluctuation(ii))) + + for kk, op in lst: + mean, stddev = _stats(op, samples) + for m, s in zip(mean.flatten(), stddev.flatten()): + print('{}: {:.02E} ± {:.02E}'.format(kk, m, s)) + + def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000): + fluctuations_slice_mean = float(fluctuations_slice_mean) + assert fluctuations_slice_mean > 0 + scm = 1. + for a in self._a: + m, std = a.fluctuation_amplitude.mean, a.fluctuation_amplitude.std + mu, sig = _lognormal_moments(m, std) + flm = np.exp(mu + sig*np.random.normal(size=nsamples)) + scm *= flm**2 + 1. + return fluctuations_slice_mean/np.mean(np.sqrt(scm)) - azm = _lognormal_moment_matching(offset_amplitude_mean, - offset_amplitude_stddev, - prefix + 'zeromode', parameter_domain, - space = tuple(self._spaces)) - - foo = np.ones(hspace.shape) - foo[zeroind] = 0 - - ZeroModeInserter = _slice_extractor(hspace, azm.target, zeroind).adjoint - - azm = Adder(from_global_data(hspace, foo)) @ ZeroModeInserter @ azm - - #NOTE ht and pd operator able to act on several spaces might be nice - ht = HarmonicTransformOperator(hspace, space = spaces[0]) - pd = PowerDistributor(hspace, - self._amplitudes[0].target[spaces[0]], spaces[0]) - for i in range(1, len(self._amplitudes)): - ht = HarmonicTransformOperator(ht.target, space = spaces[i]) @ ht - pd = pd @ PowerDistributor( pd.domain, - self._amplitudes[i].target[spaces[i]], space = spaces[i]) - - a = ContractionOperator(pd.domain, - spaces[1:]).adjoint(self._amplitudes[0]) - for i in range(1, len(self._amplitudes)): - a = a*(ContractionOperator(pd.domain, spaces[:i] + spaces[ - (i + 1):]).adjoint(self._amplitudes[i])) + @property + def amplitudes(self): + return self._a - A = pd @ a - return ht(azm*A*ducktape(hspace, None, prefix + 'xi')) + @property + def amplitude_total_offset(self): + return self._azm @property - def amplitudes(self): - return self._amplitudes + def total_fluctuation(self): + """Returns operator which acts on prior or posterior samples""" + if len(self._a) == 0: + raise NotImplementedError + if len(self._a) == 1: + return self._a[0].fluctuation_amplitude + q = 1. + for a in self._a: + fl = a.fluctuation_amplitude + q = q*(Adder(full(fl.target, 1.)) @ fl**2) + return (Adder(full(q.target, -1.)) @ q).sqrt() + + def slice_fluctuation(self, space): + """Returns operator which acts on prior or posterior samples""" + if len(self._a) == 0: + raise NotImplementedError + assert space < len(self._a) + if len(self._a) == 1: + return self._a[0].fluctuation_amplitude + q = 1. + for j in range(len(self._a)): + fl = self._a[j].fluctuation_amplitude + if j == space: + q = q*fl**2 + else: + q = q*(Adder(full(fl.target, 1.)) @ fl**2) + return q.sqrt() + + def average_fluctuation(self, space): + """Returns operator which acts on prior or posterior samples""" + if len(self._a) == 0: + raise NotImplementedError + assert space < len(self._a) + if len(self._a) == 1: + return self._a[0].fluctuation_amplitude + return self._a[space].fluctuation_amplitude + + @staticmethod + def offset_amplitude_realized(samples): + res = 0. + for s in samples: + res = res + s.mean()**2 + return np.sqrt(res/len(samples)) + + @staticmethod + def total_fluctuation_realized(samples): + return _total_fluctuation_realized(samples) + + @staticmethod + def slice_fluctuation_realized(samples, space): + """Computes slice fluctuations from collection of field (defined in signal + space) realizations.""" + ldom = len(samples[0].domain) + assert space < ldom + if ldom == 1: + return _total_fluctuation_realized(samples) + res1, res2 = 0., 0. + for s in samples: + res1 = res1 + s**2 + res2 = res2 + s.mean(space)**2 + res1 = res1/len(samples) + res2 = res2/len(samples) + res = res1.mean() - res2.mean() + return np.sqrt(res) + + + @staticmethod + def average_fluctuation_realized(samples, space): + """Computes average fluctuations from collection of field (defined in signal + space) realizations.""" + ldom = len(samples[0].domain) + assert space < ldom + if ldom == 1: + return _total_fluctuation_realized(samples) + spaces = () + for i in range(ldom): + if i != space: + spaces += (i,) + res = 0. + for s in samples: + r = s.mean(spaces) + res = res + (r - r.mean())**2 + res = res/len(samples) + return np.sqrt(res.mean()) diff --git a/nifty5/library/light_cone_operator.py b/nifty5/library/light_cone_operator.py index c680d5f0c14086e549636a6732d1bc441b0839fe..f02129d27f00962c127a44ff52c7ef91e6d29885 100644 --- a/nifty5/library/light_cone_operator.py +++ b/nifty5/library/light_cone_operator.py @@ -63,7 +63,7 @@ class _LightConeDerivative(LinearOperator): if mode == self.TIMES: res += self._derivatives[i]*x[i] else: - res[i] = np.sum(self._derivatives[i]*x) + res[i] = np.sum(self._derivatives[i]*x.real) return Field.from_global_data(self._tgt(mode), res) diff --git a/nifty5/linearization.py b/nifty5/linearization.py index 6c5b8cb20dcccd9a042dc4aa3ed36841c41fac2f..816c74d321433c0e9d9ef62f6c76637a68005459 100644 --- a/nifty5/linearization.py +++ b/nifty5/linearization.py @@ -330,6 +330,11 @@ class Linearization(object): tmp = self._val.log() return self.new(tmp, makeOp(1./self._val)(self._jac)) + def log10(self): + tmp = self._val.log10() + tmp2 = 1. / (self._val * np.log(10)) + return self.new(tmp, makeOp(tmp2)(self._jac)) + def sinh(self): tmp = self._val.sinh() tmp2 = self._val.cosh() diff --git a/nifty5/minimization/metric_gaussian_kl_mpi.py b/nifty5/minimization/metric_gaussian_kl_mpi.py index f3c280d2c9724a37ccdfda9f98e3a8c6f034743b..5c20f603ebbd7f5afa1908fd6cf6c3052218d086 100644 --- a/nifty5/minimization/metric_gaussian_kl_mpi.py +++ b/nifty5/minimization/metric_gaussian_kl_mpi.py @@ -73,7 +73,6 @@ class KLMetric(EndomorphicOperator): return self._KL.metric_sample(from_inverse, dtype) - class MetricGaussianKL_MPI(Energy): """Provides the sampled Kullback-Leibler divergence between a distribution and a Metric Gaussian. @@ -116,7 +115,7 @@ class MetricGaussianKL_MPI(Energy): _samples : None Only a parameter for internal uses. Typically not to be set by users. seed_offset : int - A parameter with which one can controll from which seed the samples + A parameter with which one can controll from which seed the samples are drawn. Per default, the seed is different for MPI tasks, but the same every time this class is initialized. @@ -132,7 +131,6 @@ class MetricGaussianKL_MPI(Energy): Torsten A. Enßlin, ``_ """ - def __init__(self, mean, hamiltonian, n_samples, constants=[], point_estimates=[], mirror_samples=False, napprox=0, _samples=None, seed_offset=0): @@ -161,6 +159,7 @@ class MetricGaussianKL_MPI(Energy): if napprox > 1: met._approximation = makeOp(approximation2endo(met, napprox)) _samples = [] + rand_state = np.random.get_state() for i in range(lo, hi): if mirror_samples: np.random.seed(i//2+seed_offset) @@ -171,8 +170,9 @@ class MetricGaussianKL_MPI(Energy): _samples.append(((i % 2)*2-1) * met.draw_sample(from_inverse=True)) else: - np.random.seed(i) + np.random.seed(i+seed_offset) _samples.append(met.draw_sample(from_inverse=True)) + np.random.set_state(rand_state) _samples = tuple(_samples) if mirror_samples: n_samples *= 2 @@ -242,8 +242,11 @@ class MetricGaussianKL_MPI(Energy): raise NotImplementedError() lin = self._lin.with_want_metric() samp = full(self._hamiltonian.domain, 0.) + rand_state = np.random.get_state() + np.random.seed(rank+np.random.randint(99999)) for v in self._samples: samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False, dtype=dtype) + np.random.set_state(rand_state) return allreduce_sum_field(samp) def metric_sample(self, from_inverse=False, dtype=np.float64): diff --git a/nifty5/operators/energy_operators.py b/nifty5/operators/energy_operators.py index 2f03e4d2a1d66bae1e3c5e97b871516f407aa6de..d6d3170ca5cc6494e7f4a2d80c27452c65bf497d 100644 --- a/nifty5/operators/energy_operators.py +++ b/nifty5/operators/energy_operators.py @@ -27,6 +27,7 @@ from .linear_operator import LinearOperator from .operator import Operator from .sampling_enabler import SamplingEnabler from .sandwich_operator import SandwichOperator +from .scaling_operator import ScalingOperator from .simple_linear_operators import VdotOperator @@ -248,6 +249,43 @@ class InverseGammaLikelihood(EnergyOperator): return res.add_metric(metric) +class StudentTEnergy(EnergyOperator): + """Computes likelihood energy of expected event frequency constrained by + event data. + + .. math :: + E(f) = -\\log \\text{Bernoulli}(d|f) + = -d^\\dagger \\log f - (1-d)^\\dagger \\log(1-f), + + where f is a field defined on `d.domain` with the expected + frequencies of events. + + Parameters + ---------- + d : Field + Data field with events (1) or non-events (0). + theta : Scalar + Degree of freedom parameter for the student t distribution + """ + + def __init__(self, domain, theta): + self._domain = DomainTuple.make(domain) + self._theta = theta + from .log1p import Log1p + self._l1p = Log1p(domain) + + def apply(self, x): + self._check_input(x) + v = ((self._theta+1)/2)*self._l1p(x**2/self._theta).sum() + if not isinstance(x, Linearization): + return Field.scalar(v) + if not x.want_metric: + return v + met = ScalingOperator((self._theta+1)/(self._theta+3), self.domain) + met = SandwichOperator.make(x.jac, met) + return v.add_metric(met) + + class BernoulliEnergy(EnergyOperator): """Computes likelihood energy of expected event frequency constrained by event data. diff --git a/nifty5/operators/exp_transform.py b/nifty5/operators/exp_transform.py deleted file mode 100644 index 47e0c7c97a6d629f7cd6cd05e9bb6303dd098200..0000000000000000000000000000000000000000 --- a/nifty5/operators/exp_transform.py +++ /dev/null @@ -1,128 +0,0 @@ -# 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 . -# -# 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)* + frac* - # 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)) diff --git a/nifty5/operators/log1p.py b/nifty5/operators/log1p.py new file mode 100644 index 0000000000000000000000000000000000000000..2dca0c3144b6ab11a9b731acc51f75aa83ee4e12 --- /dev/null +++ b/nifty5/operators/log1p.py @@ -0,0 +1,42 @@ +# 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 . +# +# Copyright(C) 2013-2019 Max-Planck-Society +# +# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. + +from ..field import Field +from ..multi_field import MultiField +from .operator import Operator +from .diagonal_operator import DiagonalOperator +from ..linearization import Linearization +from ..sugar import from_local_data +from numpy import log1p + + +class Log1p(Operator): + """computes x -> log(1+x) + """ + def __init__(self, dom): + self._domain = dom + self._target = dom + + def apply(self, x): + lin = isinstance(x, Linearization) + xval = x.val if lin else x + xlval = xval.local_data + res = from_local_data(xval.domain, log1p(xlval)) + if not lin: + return res + jac = DiagonalOperator(1/(1+xval)) + return x.new(res, jac@x.jac) diff --git a/nifty5/operators/operator.py b/nifty5/operators/operator.py index 21a8c6a891ef84e9742f1ae15e2231c78cd28fc3..90b88f784704ed9cf57c3bc2cd97d34a1785761b 100644 --- a/nifty5/operators/operator.py +++ b/nifty5/operators/operator.py @@ -159,7 +159,7 @@ class Operator(metaclass=NiftyMeta): for f in ["sqrt", "exp", "log", "tanh", "sigmoid", 'sin', 'cos', 'tan', - 'sinh', 'cosh', 'absolute', 'sinc', 'one_over']: + 'sinh', 'cosh', 'absolute', 'sinc', 'one_over', 'log10']: def func(f): def func2(self): fa = _FunctionApplier(self.target, f) diff --git a/nifty5/operators/qht_operator.py b/nifty5/operators/qht_operator.py deleted file mode 100644 index ac8f5b52c2371b471da96479a8c9835f46342d91..0000000000000000000000000000000000000000 --- a/nifty5/operators/qht_operator.py +++ /dev/null @@ -1,72 +0,0 @@ -# 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 . -# -# 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)) diff --git a/nifty5/operators/simple_linear_operators.py b/nifty5/operators/simple_linear_operators.py index 9fd8efab3fbc08f1295bfa80812ee2b949b651a0..f8c7c75d99ab69eb312c3bb22958301c66a8fb93 100644 --- a/nifty5/operators/simple_linear_operators.py +++ b/nifty5/operators/simple_linear_operators.py @@ -187,9 +187,7 @@ class _SlowFieldAdapter(LinearOperator): self._check_input(x, mode) if isinstance(x, MultiField): return x[self._name] - else: - return MultiField.from_dict({self._name: x}, - domain=self._tgt(mode)) + return MultiField.from_dict({self._name: x}, domain=self._tgt(mode)) def __repr__(self): return '_SlowFieldAdapter' @@ -338,12 +336,17 @@ class PartialExtractor(LinearOperator): if self._domain[key] is not self._target[key]: raise ValueError("domain mismatch") self._capability = self.TIMES | self.ADJOINT_TIMES + self._compldomain = MultiDomain.make({kk: self._domain[kk] + for kk in self._domain.keys() + if kk not in self._target.keys()}) def apply(self, x, mode): self._check_input(x, mode) if mode == self.TIMES: return x.extract(self._target) - return MultiField.from_dict({key: x[key] for key in x.domain.keys()}) + res0 = MultiField.from_dict({key: x[key] for key in x.domain.keys()}) + res1 = MultiField.full(self._compldomain, 0.) + return res0.unite(res1) class MatrixProductOperator(EndomorphicOperator): @@ -359,20 +362,19 @@ class MatrixProductOperator(EndomorphicOperator): `dot()` and `transpose()` in the style of numpy arrays. """ def __init__(self, domain, matrix): - self._domain = domain + self._domain = DomainTuple.make(domain) + shp = self._domain.shape + if len(shp) > 1: + raise TypeError('Only 1D-domain supported yet.') + if matrix.shape != (*shp, *shp): + raise ValueError self._capability = self.TIMES | self.ADJOINT_TIMES - self._mat = matrix - self._mat_tr = matrix.transpose() + self._mat_tr = matrix.transpose().conjugate() def apply(self, x, mode): self._check_input(x, mode) res = x.to_global_data() - if mode == self.TIMES: - res = self._mat.dot(res) - if mode == self.ADJOINT_TIMES: - res = self._mat_tr.dot(res) + f = self._mat.dot if mode == self.TIMES else self._mat_tr.dot + res = f(res) return Field.from_global_data(self._domain, res) - - def __repr__(self): - return "MatrixProductOperator" diff --git a/nifty5/operators/slope_operator.py b/nifty5/operators/slope_operator.py deleted file mode 100644 index c9486e488f18ecf085878338b42d34753c8bdf89..0000000000000000000000000000000000000000 --- a/nifty5/operators/slope_operator.py +++ /dev/null @@ -1,67 +0,0 @@ -# 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 . -# -# 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) diff --git a/nifty5/operators/symmetrizing_operator.py b/nifty5/operators/symmetrizing_operator.py deleted file mode 100644 index c0abb67b4a656f8a57176907ed699087d2602d47..0000000000000000000000000000000000000000 --- a/nifty5/operators/symmetrizing_operator.py +++ /dev/null @@ -1,53 +0,0 @@ -# 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 . -# -# 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)) diff --git a/nifty5/plot.py b/nifty5/plot.py index 7a2d1600860c52d9735af66cf5dc25cd864aa206..b53f232fbb982c8e15a93e126230284ab9638809 100644 --- a/nifty5/plot.py +++ b/nifty5/plot.py @@ -22,7 +22,6 @@ import numpy as np from . import dobj from .domains.gl_space import GLSpace from .domains.hp_space import HPSpace -from .domains.log_rg_space import LogRGSpace from .domains.power_space import PowerSpace from .domains.rg_space import RGSpace from .field import Field @@ -171,14 +170,15 @@ def _find_closest(A, target): return idx -def _makeplot(name): +def _makeplot(name, block=True): import matplotlib.pyplot as plt if dobj.rank != 0: plt.close() return if name is None: - plt.show() - plt.close() + plt.show(block=block) + if block: + plt.close() return extension = os.path.splitext(name)[1] if extension in (".pdf", ".png", ".svg"): @@ -306,18 +306,6 @@ def _plot1D(f, ax, **kwargs): if label != ([None]*len(f)): plt.legend() 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): plt.xscale(kwargs.pop("xscale", "log")) plt.yscale(kwargs.pop("yscale", "log")) @@ -432,8 +420,8 @@ def _plot(f, ax, **kwargs): dom1 = f[0].domain if (len(dom1) == 1 and (isinstance(dom1[0], PowerSpace) or - (isinstance(dom1[0], (RGSpace, LogRGSpace)) and - len(dom1[0].shape) == 1))): + (isinstance(dom1[0], RGSpace) and + len(dom1[0].shape) == 1))): _plot1D(f, ax, **kwargs) return else: @@ -511,6 +499,9 @@ class Plot(object): If left empty, the plot will be shown on the screen, otherwise it will be written to a file with the given name. Supported extensions: .png and .pdf. Default: None. + block: bool + Override the blocking behavior of the non-interactive plotting + mode. The plot will not be closed in this case but is left open! """ import matplotlib.pyplot as plt nplot = len(self._plots) @@ -537,4 +528,4 @@ class Plot(object): ax = fig.add_subplot(ny, nx, i+1) _plot(self._plots[i], ax, **self._kwargs[i]) fig.tight_layout() - _makeplot(kwargs.pop("name", None)) + _makeplot(kwargs.pop("name", None), block=kwargs.pop("block", True)) diff --git a/nifty5/sugar.py b/nifty5/sugar.py index 1f9257a964eb18c8c481360c7bcfa78ccdd0eea3..c2d5142710b34437ba69e696e3038fed33e16127 100644 --- a/nifty5/sugar.py +++ b/nifty5/sugar.py @@ -35,7 +35,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator', 'create_harmonic_smoothing_operator', 'from_random', 'full', 'from_global_data', 'from_local_data', 'makeDomain', 'sqrt', 'exp', 'log', 'tanh', 'sigmoid', - 'sin', 'cos', 'tan', 'sinh', 'cosh', + 'sin', 'cos', 'tan', 'sinh', 'cosh', 'log10', 'absolute', 'one_over', 'clip', 'sinc', 'conjugate', 'get_signal_variance', 'makeOp', 'domain_union', 'get_default_codomain', 'single_plot'] @@ -389,7 +389,7 @@ def domain_union(domains): _current_module = sys.modules[__name__] -for f in ["sqrt", "exp", "log", "tanh", "sigmoid", +for f in ["sqrt", "exp", "log", "log10", "tanh", "sigmoid", "conjugate", 'sin', 'cos', 'tan', 'sinh', 'cosh', 'absolute', 'one_over', 'sinc']: def func(f): diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py index cd4d1b50e43cc90a5f40b6e19d84a605a3f5b5a7..0e6a9371f0a40ea69be5816d64a02e6fc0ec5f49 100644 --- a/test/test_energy_gradients.py +++ b/test/test_energy_gradients.py @@ -47,13 +47,18 @@ def test_gaussian(field): ift.extra.check_jacobian_consistency(energy, field) +def test_studentt(field): + energy = ift.StudentTEnergy(domain=field.domain, theta=.5) + ift.extra.check_jacobian_consistency(energy, field, tol=1e-6) + + def test_inverse_gamma(field): field = field.exp() space = field.domain d = np.random.normal(10, size=space.shape)**2 d = ift.Field.from_global_data(space, d) energy = ift.InverseGammaLikelihood(d) - ift.extra.check_jacobian_consistency(energy, field, tol=1e-7) + ift.extra.check_jacobian_consistency(energy, field, tol=1e-5) def testPoissonian(field): @@ -83,4 +88,4 @@ def test_bernoulli(field): d = np.random.binomial(1, 0.1, size=space.shape) d = ift.Field.from_global_data(space, d) energy = ift.BernoulliEnergy(d) - ift.extra.check_jacobian_consistency(energy, field, tol=1e-6) + ift.extra.check_jacobian_consistency(energy, field, tol=1e-5) diff --git a/test/test_field.py b/test/test_field.py index 325cae05ffac98c810988315935d98c2c577868c..a39d214da0f95317b2bb606652ad9f3687ec36fd 100644 --- a/test/test_field.py +++ b/test/test_field.py @@ -219,17 +219,23 @@ def test_weight(): f = ift.Field.full(s1, 10.) f2 = f.weight(1) assert_equal(f.weight(1).local_data, f2.local_data) + assert_equal(f.domain.total_volume(), 1) + assert_equal(f.domain.total_volume(0), 1) + assert_equal(f.domain.total_volume((0,)), 1) assert_equal(f.total_volume(), 1) assert_equal(f.total_volume(0), 1) assert_equal(f.total_volume((0,)), 1) + assert_equal(f.domain.scalar_weight(), 0.1) + assert_equal(f.domain.scalar_weight(0), 0.1) + assert_equal(f.domain.scalar_weight((0,)), 0.1) assert_equal(f.scalar_weight(), 0.1) assert_equal(f.scalar_weight(0), 0.1) assert_equal(f.scalar_weight((0,)), 0.1) s1 = ift.GLSpace(10) f = ift.Field.full(s1, 10.) - assert_equal(f.scalar_weight(), None) - assert_equal(f.scalar_weight(0), None) - assert_equal(f.scalar_weight((0,)), None) + assert_equal(f.domain.scalar_weight(), None) + assert_equal(f.domain.scalar_weight(0), None) + assert_equal(f.domain.scalar_weight((0,)), None) @pmp('dom', [ift.RGSpace(10), ift.GLSpace(10)]) @@ -329,7 +335,7 @@ def test_emptydomain(): @pmp('dom', [ift.RGSpace((8,), harmonic=True), ()]) @pmp('func', [ "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "sinc", "absolute", - "sign" + "sign", "log10" ]) def test_funcs(num, dom, func): num = 5 diff --git a/test/test_linearization.py b/test/test_linearization.py index f9973900ec9055db59fcfeb20b3935fbbeaf1890..fa0c49e68daa88d4170fcf13b52b8bf1dcd45651 100644 --- a/test/test_linearization.py +++ b/test/test_linearization.py @@ -54,7 +54,7 @@ def test_special_gradients(): @pmp('f', [ 'log', 'exp', 'sqrt', 'sin', 'cos', 'tan', 'sinc', 'sinh', 'cosh', 'tanh', - 'absolute', 'one_over', 'sigmoid' + 'absolute', 'one_over', 'sigmoid', 'log10' ]) def test_actual_gradients(f): dom = ift.UnstructuredDomain((1,)) diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py index a6782a6c0d96d81498c0263724692af4de1873d3..3cfddd04d894003d06bb6dbf5e840f06f7a1a133 100644 --- a/test/test_operators/test_adjoint.py +++ b/test/test_operators/test_adjoint.py @@ -32,7 +32,11 @@ _p_RG_spaces = [ ift.RGSpace((1, 2, 3, 6), distances=(0.2, 0.25, 0.34, .8)) ] _p_spaces = _p_RG_spaces + [ift.HPSpace(17), ift.GLSpace(8, 13)] -_pow_spaces = [ift.PowerSpace(ift.RGSpace((17, 38), harmonic=True))] +_pow_spaces = [ + ift.PowerSpace(ift.RGSpace((17, 38), (0.99, 1340), harmonic=True)), + ift.PowerSpace(ift.LMSpace(18), + ift.PowerSpace.useful_binbounds(ift.LMSpace(18), False)) +] pmp = pytest.mark.parametrize dtype = list2fixture([np.float64, np.complex128]) @@ -87,16 +91,6 @@ def testConjugationOperator(sp): 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) def testOperatorAdaptor(sp, dtype): op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype)) @@ -214,14 +208,6 @@ def testDomainTupleFieldInserter(): 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) - - @pmp('space', [0, 2]) @pmp('factor', [1, 2, 2.7]) @pmp('central', [False, True]) @@ -233,28 +219,6 @@ def testZeroPadder(space, factor, dtype, central): ift.extra.consistency_check(op, dtype, dtype) -@pmp('args', - [(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)]) -def testExpTransform(args, dtype): - op = ift.ExpTransform(args[0], args[1], args[2]) - ift.extra.consistency_check(op, dtype, dtype) - - -@pmp('args', - [(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)]) -def testQHTOperator(args): - dtype = np.float64 - tgt = ift.DomainTuple.make(args[0]) - op = ift.QHTOperator(tgt, args[1]) - ift.extra.consistency_check(op, dtype, dtype) - - @pmp('args', [[ift.RGSpace( (13, 52, 40)), (4, 6, 25), None], [ift.RGSpace( (128, 128)), (45, 48), 0], [ift.RGSpace(13), (7,), None], [ @@ -295,3 +259,52 @@ def testValueInserter(sp, seed): ind.append(np.random.randint(0, ss-1)) op = ift.ValueInserter(sp, ind) ift.extra.consistency_check(op) + + +@pmp('sp', _pow_spaces) +def testSlopeRemover(sp): + op = ift.library.correlated_fields._SlopeRemover(sp) + ift.extra.consistency_check(op) + + +@pmp('sp', _pow_spaces) +def testTwoLogIntegrations(sp): + op = ift.library.correlated_fields._TwoLogIntegrations(sp) + ift.extra.consistency_check(op) + + +@pmp('sp', _h_spaces + _p_spaces + _pow_spaces) +def testSpecialSum(sp): + op = ift.library.correlated_fields._SpecialSum(sp) + ift.extra.consistency_check(op) + + +@pmp('sp', [ift.RGSpace(10)]) +@pmp('seed', [12, 3]) +def testMatrixProductOperator(sp, seed): + np.random.seed(seed) + mat = np.random.randn(*sp.shape, *sp.shape) + op = ift.MatrixProductOperator(sp, mat) + ift.extra.consistency_check(op) + mat = mat + 1j*np.random.randn(*sp.shape, *sp.shape) + op = ift.MatrixProductOperator(sp, mat) + ift.extra.consistency_check(op) + + +@pmp('seed', [12, 3]) +def testPartialExtractor(seed): + np.random.seed(seed) + tgt = {'a': ift.RGSpace(1), 'b': ift.RGSpace(2)} + dom = tgt.copy() + dom['c'] = ift.RGSpace(3) + dom = ift.MultiDomain.make(dom) + tgt = ift.MultiDomain.make(tgt) + op = ift.PartialExtractor(dom, tgt) + ift.extra.consistency_check(op) + + +@pmp('seed', [12, 3]) +def testSlowFieldAdapter(seed): + dom = {'a': ift.RGSpace(1), 'b': ift.RGSpace(2)} + op = ift.operators.simple_linear_operators._SlowFieldAdapter(dom, 'a') + ift.extra.consistency_check(op) diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py new file mode 100644 index 0000000000000000000000000000000000000000..4beb608244b1791296611f6b8b19a4334177b8ab --- /dev/null +++ b/test/test_operators/test_correlated_fields.py @@ -0,0 +1,89 @@ +# 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 . +# +# Copyright(C) 2013-2019 Max-Planck-Society +# +# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. + +import pytest +from numpy.testing import assert_allclose +from numpy.random import seed + +import nifty5 as ift + + +@pytest.mark.parametrize('sspace', [ + ift.RGSpace(4), + ift.RGSpace((4, 4), (0.123,0.4)), + ift.HPSpace(8), + ift.GLSpace(4) +]) +@pytest.mark.parametrize('rseed', [13, 2]) +@pytest.mark.parametrize('Astds', [[1.,3.],[0.2,1.4]]) +@pytest.mark.parametrize('offset_std', [1.,10.]) +def testAmplitudesConsistency(rseed, sspace, Astds, offset_std): + def stats(op,samples): + sc = ift.StatCalculator() + for s in samples: + sc.add(op(s.extract(op.domain))) + return sc.mean.to_global_data(), sc.var.sqrt().to_global_data() + seed(rseed) + nsam = 100 + + fsspace = ift.RGSpace((12,), (0.4,)) + + fa = ift.CorrelatedFieldMaker() + fa.add_fluctuations(sspace, Astds[0], 1E-8, 1.1, 2., 2.1, .5, + -2, 1., 'spatial') + fa.add_fluctuations(fsspace, Astds[1], 1E-8, 3.1, 1., .5, .1, + -4, 1., 'freq') + op = fa.finalize(offset_std, 1E-8, '') + + samples = [ift.from_random('normal',op.domain) for _ in range(nsam)] + tot_flm, _ = stats(fa.total_fluctuation,samples) + offset_std,_ = stats(fa.amplitude_total_offset,samples) + intergated_fluct_std0,_ = stats(fa.average_fluctuation(0),samples) + intergated_fluct_std1,_ = stats(fa.average_fluctuation(1),samples) + + slice_fluct_std0,_ = stats(fa.slice_fluctuation(0),samples) + slice_fluct_std1,_ = stats(fa.slice_fluctuation(1),samples) + + sams = [op(s) for s in samples] + fluct_total = fa.total_fluctuation_realized(sams) + fluct_space = fa.average_fluctuation_realized(sams,0) + fluct_freq = fa.average_fluctuation_realized(sams,1) + zm_std_mean = fa.offset_amplitude_realized(sams) + sl_fluct_space = fa.slice_fluctuation_realized(sams,0) + sl_fluct_freq = fa.slice_fluctuation_realized(sams,1) + + assert_allclose(offset_std, zm_std_mean, rtol=0.5) + assert_allclose(intergated_fluct_std0, fluct_space, rtol=0.5) + assert_allclose(intergated_fluct_std1, fluct_freq, rtol=0.5) + assert_allclose(tot_flm, fluct_total, rtol=0.5) + assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5) + assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5) + + fa = ift.CorrelatedFieldMaker() + fa.add_fluctuations(fsspace, Astds[1], 1., 3.1, 1., .5, .1, + -4, 1., 'freq') + m = 3. + x = fa.moment_slice_to_average(m) + fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, + -2, 1., 'spatial', 0) + op = fa.finalize(offset_std, .1, '') + em, estd = stats(fa.slice_fluctuation(0),samples) + + assert_allclose(m, em, rtol=0.5) + + assert op.target[0] == sspace + assert op.target[1] == fsspace \ No newline at end of file