GitLab

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)