diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py index bf435fec4e9af306ca17bf33e7a038171c700aa5..e088f427ea8fdd6e2d34d815055442311a56f944 100644 --- a/demos/getting_started_3.py +++ b/demos/getting_started_3.py @@ -46,32 +46,32 @@ if __name__ == '__main__': mode = 1 position_space = ift.RGSpace([128, 128]) + harmonic_space = position_space.get_default_codomain() + ht = ift.HarmonicTransformOperator(harmonic_space, position_space) + power_space = ift.PowerSpace(harmonic_space) # Set up an amplitude operator for the field - # The parameters mean: - # 64 spectral bins - # - # Spectral smoothness (affects Gaussian process part) - # 3 = relatively high variance of spectral curbvature - # 0.4 = quefrency mode below which cepstrum flattens - # - # Power-law part of spectrum: - # -5 = preferred power-law slope - # 0.5 = low variance of power-law slope - # 0.4 = y-intercept mean - # 0.3 = relatively high y-intercept variance - A = ift.AmplitudeOperator(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3) + dct = { + 'target': power_space, + 'n_pix': 64, # 64 spectral bins + + # Spectral smoothness (affects Gaussian process part) + 'a': 3, # relatively high variance of spectral curbvature + 'k0': .4, # quefrency mode below which cepstrum flattens + + # Power-law part of spectrum: + 'sm': -5, # preferred power-law slope + 'sv': .5, # low variance of power-law slope + 'im': .4, # y-intercept mean + 'iv': .3 # relatively high y-intercept variance + } + A = ift.AmplitudeOperator(**dct) # Build the operator for a correlated signal - harmonic_space = position_space.get_default_codomain() - ht = ift.HarmonicTransformOperator(harmonic_space, position_space) - power_space = A.target[0] power_distributor = ift.PowerDistributor(harmonic_space, power_space) - - vol = ift.ScalingOperator(harmonic_space.scalar_dvol**(-0.5), - harmonic_space) - correlated_field = ht( - vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi')) + vol = harmonic_space.scalar_dvol**-0.5 + xi = ift.ducktape(harmonic_space, None, 'xi') + correlated_field = ht(vol*power_distributor(A)*xi) # Alternatively, one can use: # correlated_field = ift.CorrelatedField(position_space, A) diff --git a/nifty5/library/amplitude_operator.py b/nifty5/library/amplitude_operator.py index a3b2fbfab5b3646a6bca6fb9c4326fb6dad1b836..c9e6172447709b4d1e5cb6d4816e719ea1c55d76 100644 --- a/nifty5/library/amplitude_operator.py +++ b/nifty5/library/amplitude_operator.py @@ -19,76 +19,51 @@ import numpy as np from ..domains.power_space import PowerSpace from ..field import Field -from ..sugar import makeOp, sqrt +from ..sugar import makeOp def _ceps_kernel(dof_space, k, a, k0): return a**2/(1 + (k/k0)**2)**2 -def create_cepstrum_amplitude_field(domain, cepstrum): - """Creates a ... - Writes the sum of all modes into the zero-mode. - - Parameters - ---------- - domain: ??? - ??? - cepstrum: Callable - ??? - """ - +def _create_cepstrum_amplitude_field(domain, cepstrum): dim = len(domain.shape) shape = domain.shape - q_array = domain.get_k_array() - # Fill cepstrum field (all non-zero modes) - no_zero_modes = (slice(1, None),)*dim - ks = q_array[(slice(None),) + no_zero_modes] + # Fill all non-zero modes + no_zero_modes = (slice(1, None), )*dim + ks = q_array[(slice(None), ) + no_zero_modes] cepstrum_field = np.zeros(shape) cepstrum_field[no_zero_modes] = cepstrum(ks) - # Fill cepstrum field (zero-mode subspaces) + # Fill zero-mode subspaces for i in range(dim): - # Prepare indices - fst_dims = (slice(None),)*i - sl = fst_dims + (slice(1, None),) - sl2 = fst_dims + (0,) - - # Do summation + fst_dims = (slice(None), )*i + sl = fst_dims + (slice(1, None), ) + sl2 = fst_dims + (0, ) cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i) - return Field.from_global_data(domain, cepstrum_field) -def _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True): +def CepstrumOperator(domain, a, k0): ''' - Parameters - ---------- - ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel - eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 - a = ceps_a, k0 = ceps_k0 + .. math:: + C(k) = \\left(\\frac{a}{1+(k/k0)^2}\\right)^2 ''' - from ..operators.qht_operator import QHTOperator from ..operators.symmetrizing_operator import SymmetrizingOperator - qht = QHTOperator(target=logk_space) + + # FIXME a>0 k0>0 + qht = QHTOperator(target=domain) dof_space = qht.domain[0] - sym = SymmetrizingOperator(logk_space) - kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k) - cepstrum = create_cepstrum_amplitude_field(dof_space, kern) - res = sym(qht(makeOp(sqrt(cepstrum)))) - if not zero_mode: - shp = res.target.shape - mask = np.ones(shp) - mask[(0,)*len(shp)] = 0. - mask = makeOp(Field.from_global_data(res.target, mask)) - res = mask(res) - return res - - -def _SlopePowerSpectrum(logk_space, sm, sv, im, iv): + sym = SymmetrizingOperator(domain) + kern = lambda k: _ceps_kernel(dof_space, k, a, k0) + cepstrum = _create_cepstrum_amplitude_field(dof_space, kern) + return sym @ qht @ makeOp(cepstrum.sqrt()) + + +def SlopeOperator(domain, sm, sv, im, iv): ''' Parameters ---------- @@ -97,54 +72,90 @@ def _SlopePowerSpectrum(logk_space, sm, sv, im, iv): slope_variance (default=1) im, iv : y-intercept_mean, y-intercept_std of power_slope - ''' + ''' from ..operators.slope_operator import SlopeOperator from ..operators.offset_operator import OffsetOperator - phi_mean = np.array([sm, im + sm*logk_space.t_0[0]]) + + # sv, iv>0 + + phi_mean = np.array([sm, im + sm*domain.t_0[0]]) phi_sig = np.array([sv, iv]) - slope = SlopeOperator(logk_space) + slope = SlopeOperator(domain) phi_mean = Field.from_global_data(slope.domain, phi_mean) phi_sig = Field.from_global_data(slope.domain, phi_sig) return slope(OffsetOperator(phi_mean)(makeOp(phi_sig))) -def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv, - keys=['tau', 'phi'], zero_mode=True): - ''' Operator for parametrizing smooth power spectra. +def AmplitudeOperator( + target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']): + '''Operator for parametrizing smooth power spectra. + + The general guideline for setting up generative models in IFT is to + transform the problem into the eigenbase of the prior and formulate the + generative model in this base. This is done here for the case of a power + spectrum which is smooth and has a linear component (both on + double-logarithmic scale). + + This function assembles an :class:`Operator` which maps two a-priori white + Gaussian random fields to a smooth power spectrum which is composed out of + a linear and a smooth component. + + On double-logarithmic scale, i.e. both x and y-axis on logarithmic scale, + the output of the generated operator is: + + AmplitudeOperator = 0.5*(smooth_component + linear_component) - Computes a smooth power spectrum. - Output is defined on a PowerSpace. + This is then exponentiated and exponentially binned (via a :class:`ExpTransform`). + + The prior on the linear component is parametrized by four real numbers, + being expected value and prior variance on the slope and the y-intercept + of the linear function. + + The prior on the smooth component is parametrized by two real numbers: the + strength and the cutoff of the smoothness prior (see :class:`CepstrumOperator`). Parameters ---------- - Npixdof : int - #pix in dof_space - ceps_a : float - Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a, k0 = ceps_k0 - ceps_k0 : float - Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a, k0 = ceps_k0 + n_pix : int + Number of pixels of the space in which the . + target : PowerSpace + Target of the Operator. + a : float + Strength of smoothness prior (see :class:`CepstrumOperator`). + k0 : float + Cutoff of smothness prior in quefrency space (see :class:`CepstrumOperator`). sm : float - slope_mean = expected exponent of power law (e.g. -4) + Expected exponent of power law (see :class:`SlopeOperator`). sv : float - slope_variance (default=1) + Prior variance of exponent of power law (see :class:`SlopeOperator`). im : float - y-intercept_mean + Expected y-intercept of power law (see :class:`SlopeOperator`). iv : float - y-intercept_variance of power_slope + Prior variance of y-intercept of power law (see :class:`SlopeOperator`). + + Returns + ------- + Operator + Operator which is defined on the space of white excitations fields and + which returns on its target a power spectrum which consists out of a + smooth and a linear part. ''' - from ..operators.exp_transform import ExpTransform - from ..operators.scaling_operator import ScalingOperator - h_space = s_space.get_default_codomain() - et = ExpTransform(PowerSpace(h_space), Npixdof) - logk_space = et.domain[0] + if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)): + raise TypeError + + a, k0 = float(a), float(k0) + sm, sv, im, iv = float(sm), float(sv), float(im), float(iv) + + et = ExpTransform(target, n_pix) + dom = et.domain[0] + + dct = {'a': a, 'k0': k0} + smooth = CepstrumOperator(dom, **dct).ducktape(keys[0]) - smooth = _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode) - smooth = smooth.ducktape(keys[0]) - linear = _SlopePowerSpectrum(logk_space, sm, sv, im, iv) - linear = linear.ducktape(keys[1]) + dct = {'sm': sm, 'sv': sv, 'im': im, 'iv': iv} + linear = SlopeOperator(dom, **dct).ducktape(keys[1]) - fac = ScalingOperator(0.5, smooth.target) - return et((fac(smooth + linear)).exp()) + return et((0.5*(smooth + linear)).exp()) diff --git a/test/test_model_gradients.py b/test/test_model_gradients.py index 340f2b98b6c0964dfa2291dfd8a22e9ed0383e6c..90fbdd4e5fe2dda2d341f69caaaf8a704cd92eb6 100644 --- a/test/test_model_gradients.py +++ b/test/test_model_gradients.py @@ -17,6 +17,7 @@ import numpy as np import pytest +from numpy.testing import assert_ import nifty5 as ift @@ -53,11 +54,11 @@ def testBasics(space, seed): @pmp('type2', ['Variable']) def testBinary(type1, type2, space, seed): dom1 = ift.MultiDomain.make({'s1': space}) - # FIXME Remove? - lin1 = _make_linearization(type1, dom1, seed) dom2 = ift.MultiDomain.make({'s2': space}) - # FIXME Remove? - lin2 = _make_linearization(type2, dom2, seed) + + # FIXME Remove this? + _make_linearization(type1, dom1, seed) + _make_linearization(type2, dom2, seed) dom = ift.MultiDomain.union((dom1, dom2)) select_s1 = ift.ducktape(None, dom, "s1") @@ -74,8 +75,7 @@ def testBinary(type1, type2, space, seed): model = ift.ScalingOperator(2.456, space)(select_s1*select_s2) pos = ift.from_random("normal", dom) ift.extra.check_value_gradient_consistency(model, pos, ntries=20) - model = ift.sigmoid( - ift.ScalingOperator(2.456, space)(select_s1*select_s2)) + model = ift.sigmoid(2.456*(select_s1*select_s2)) pos = ift.from_random("normal", dom) ift.extra.check_value_gradient_consistency(model, pos, ntries=20) pos = ift.from_random("normal", dom) @@ -91,8 +91,10 @@ def testModelLibrary(space, seed): # Tests amplitude model and coorelated field model Npixdof, ceps_a, ceps_k, sm, sv, im, iv = 4, 0.5, 2., 3., 1.5, 1.75, 1.3 np.random.seed(seed) - model = ift.AmplitudeOperator(space, Npixdof, ceps_a, ceps_k, sm, sv, im, + domain = ift.PowerSpace(space.get_default_codomain()) + model = ift.AmplitudeOperator(domain, Npixdof, ceps_a, ceps_k, sm, sv, im, iv) + assert_(isinstance(model, ift.Operator)) S = ift.ScalingOperator(1., model.domain) pos = S.draw_sample() ift.extra.check_value_gradient_consistency(model, pos, ntries=20) @@ -112,28 +114,38 @@ def testPointModel(space, seed): # FIXME All those cdfs and ppfs are not very accurate ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2, ntries=20) -@pmp('domain', [ift.RGSpace(64, distances=.789), - ift.RGSpace([32, 32], distances=.789), - ift.RGSpace([32, 32, 8], distances=.789)]) + +@pmp('domain', [ + ift.RGSpace(64, distances=.789), + ift.RGSpace([32, 32], distances=.789), + ift.RGSpace([32, 32, 8], distances=.789) +]) @pmp('causal', [True, False]) @pmp('minimum_phase', [True, False]) @pmp('seed', [4, 78, 23]) def testDynamicModel(domain, causal, minimum_phase, seed): - model, _ = ift.dynamic_operator(domain,None,1.,1.,'f', - causal = causal, - minimum_phase = minimum_phase) + model, _ = ift.dynamic_operator( + domain, None, 1., 1., 'f', causal=causal, minimum_phase=minimum_phase) S = ift.ScalingOperator(1., model.domain) pos = S.draw_sample() # FIXME I dont know why smaller tol fails for 3D example - ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, - ntries=20) + ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, ntries=20) if len(domain.shape) > 1: - model, _ = ift.dynamic_lightcone_operator(domain,None,3.,1., - 'f','c',1.,5, - causal = causal, - minimum_phase = minimum_phase) + dct = { + 'domain': domain, + 'harmonic_padding': None, + 'sm_s0': 3., + 'sm_x0': 1., + 'key': 'f', + 'lightcone_key': 'c', + 'sigc': 1., + 'quant': 5, + 'causal': causal, + 'minimum_phase': minimum_phase + } + model, _ = ift.dynamic_lightcone_operator(**dct) S = ift.ScalingOperator(1., model.domain) pos = S.draw_sample() # FIXME I dont know why smaller tol fails for 3D example - ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, - ntries=20) + ift.extra.check_value_gradient_consistency( + model, pos, tol=1e-5, ntries=20)