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 Field.power_synthesize()
.