Error found by dimensional analysis
The test case below was distilled from
wiener_filter_via_curvature.py by Torsten and me.
It indicates that there is a unit error in the way we compute the mock signal.
import numpy as np import nifty as ift import numericalunits as nu np.random.seed(43) correlation_length = 0.05*nu.m field_sigma = 2.* nu.K response_sigma = 0.01*nu.m # stupid power spectrum, but we only care about units now def power_spectrum(k): a = 4 * correlation_length * field_sigma**2 return a / (1 + k * correlation_length) ** 4 # Total side-length of the domain L = 2.*nu.m # Grid resolution (pixels per axis) N_pixels = 512 signal_space = ift.RGSpace([N_pixels], distances=L/N_pixels) harmonic_space = ift.FFTOperator.get_default_codomain(signal_space) fft = ift.FFTOperator(harmonic_space, target=signal_space) power_space = ift.PowerSpace(harmonic_space) mock_power = ift.Field(power_space, val=power_spectrum(power_space.kindex)) mock_harmonic = mock_power.power_synthesize(real_signal=True) mock_signal = fft(mock_harmonic) print "signal mean (measured in K):", mock_signal.mean()/nu.K print "signal mean (measured in K/sqrt(m)): ", mock_signal.mean()*np.sqrt(nu.m)/nu.K
If this code is run several times (I tested with the
nightly branch), the first printed number changes, but the second one does not, which means that the unit for the mock_signal field is actually sqrt(m)/K, which is not the intended one.
My guess is that we either have a dimensional error in
power_spectrum() or in