Drawing gaussian random fields over multiple domains
Hi,
currently I am facing some inconsistencies for drawing gaussian random fields over multiple domains. Yet I did not find the root of these inconsistencies. However in the minimal test example below, s_1
and s_2
do not match, neither mean nor variance even though they should be drawn from the same gaussian just in two different ways. Further none of the drawn fields, s_1
s_2
, obey the desired spectra.
What am I doing wrong or missing?
from nifty import *
x_0 = RGSpace(50, zerocenter=False, distances=.5)
k_0 = RGRGTransformation.get_codomain(x_0)
p_0 = PowerSpace(harmonic_partner=k_0)
x_1 = RGSpace(75, zerocenter=False, distances=.25)
k_1 = RGRGTransformation.get_codomain(x_1)
p_1 = PowerSpace(harmonic_partner=k_1)
fft_0 = FFTOperator(domain=x_0, target=k_0)
fft_1 = FFTOperator(domain=x_1, target=k_1)
p_spec_0 = Field(p_0, val=lambda l: 1/(1+l)**3)
p_spec_1 = Field(p_1, val=lambda k: 10/(1+k)**5)
def random_field_method_1():
w = Field((x_0, x_1), val=0.)
for ii in xrange(100):
outer = Field((p_0, p_1), val=np.outer(p_spec_0.val, p_spec_1.val))
sk = outer.power_synthesize()
w += fft_0.inverse_times(fft_1.inverse_times(sk, spaces=1), spaces=0)
return w/100.
def random_field_method_2():
w = Field((x_0, x_1), val=0.)
S_0 = create_power_operator(k_0, sqrt(p_spec_0))
S_1 = create_power_operator(k_1, sqrt(p_spec_1))
for ii in xrange(100):
rand = Field.from_random('normal', domain=(x_0, x_1))
rand_k = fft_0.times(fft_1.times(rand, spaces=1), spaces=0)
sk = S_0.times(S_1.times(rand_k, spaces=1), spaces=0)
w += fft_0.inverse_times(fft_1.inverse_times(sk, spaces=1), spaces=0)
return w/100.
def analyze_power(signal):
sk = fft_0.times(fft_1.times(signal,spaces=1), spaces=0)
power = sk.power_analyze()
p_drawn_0 = power.sum(spaces=1)/p_spec_1.sum()
p_drawn_1 = power.sum(spaces=0)/p_spec_0.sum()
return p_drawn_0, p_drawn_1
s_1 = random_field_method_1()
s_2 = random_field_method_2()
method_1_power_0, method_1_power_1 = analyze_power(s_1)
methdo_2_power_0, method_2_power_2 = analyze_power(s_2)