ProductSpaces drawing random fields and checking their spectra
Hi, currently I am testing ProductSpace reconstructions (PLN over RGSpace (x) RGSpace). I am facing a problem concerning the spectra. For a simple test scenario I would like to check the power of field drawn from two spectra, as the following
from nifty import *
np.random.seed(42)
#########################################################################
# setting space one
#########################################################################
npix_0 = np.array([200]) # number of pixels
total_volume_0 = 1. # total length
# setting signal parameters
lambda_s_0 = .5 # signal correlation length
sigma_s_0 = 1. # signal variance
#setting length smoothing kernel
kernel_0 = 0.
# calculating parameters
k0_0 = 4./(2*np.pi*lambda_s_0)
a_s_0 = sigma_s_0**2.*lambda_s_0*total_volume_0
# creation of spaces
x_0 = RGSpace(npix_0, distances=total_volume_0/npix_0,
zerocenter=False)
k_0 = RGRGTransformation.get_codomain(x_0)
p_0 = PowerSpace(harmonic_partner=k_0, logarithmic=False)
# creating Power Operator with given spectrum
spec_0 = (lambda k: a_s_0/(1+(k/k0_0)**2)**2)
p_field_0 = Field(p_0, val=spec_0, dtype=np.float64)
# creating FFT-Operator and Response-Operator with Gaussian convolution
FFTOp_0 = FFTOperator(domain=x_0, target=k_0)
#########################################################################
# setting space two
#########################################################################
npix_1 = np.array([100]) # number of pixels
total_volume_1 = 1. # total length
# setting signal parameters
lambda_s_1 = .05 # signal correlation length
sigma_s_1 = 2. # signal variance
#setting length smoothing kernel
kernel_1 = 0.
# calculating parameters
k0_1 = 4./(2*np.pi*lambda_s_1)
a_s_1 = sigma_s_1**2.*lambda_s_1*total_volume_1
# creation of spaces
x_1 = RGSpace(npix_1, distances=total_volume_1/npix_1,
zerocenter=False)
k_1 = RGRGTransformation.get_codomain(x_1)
p_1 = PowerSpace(harmonic_partner=k_1, logarithmic=False)
# creating Power Operator with given spectrum
spec_1 = (lambda k: a_s_0/(1+(k/k0_1)**2)**2)
p_field_1 = Field(p_1, val=spec_1, dtype=np.float64)
# creating FFT-Operator and Response-Operator with Gaussian convolution
FFTOp_1 = FFTOperator(domain=x_1, target=k_1)
#########################################################################
# drawing a random field with above defined correlation structure
#########################################################################
fp_0 = Field(p_0, val=(lambda k: spec_0(k)**1/2.))
fp_1 = Field(p_1, val=(lambda k: spec_1(k)**1/2.))
outer = np.outer(fp_0.val.data, fp_1.val.data)
fp = Field((p_0, p_1), val=outer)
sk = fp.power_synthesize(spaces=(0,1), real_signal=True)
The question is, how do I recover fp_0 and fp_1 from sk (i.e. simply the power of the drawn field).
In my mind this should be:
def test_product_space_power(x):
samples = 100
ps_0 = Field(p_0, val=0.)
ps_1 = Field(p_1, val=1.)
for ii in xrange(samples):
sk = x.power_synthesize(spaces=(0,1), real_signal=True)
p0, p1 = analyze_sk(sk)
ps_1 += p1
ps_0 += p0
p_0_test_result = ps_0/samples
p_1_test_result = ps_1/samples
return p_0_test_result, p_1_test_result
def analyze_sk(x):
sp = x.power_analyze(spaces=0, decompose_power=False)\
.power_analyze(spaces=1, decompose_power=False)
p0= sp.sum(spaces=1)/fp_1.sum()
p1= sp.sum(spaces=0)/fp_0.sum()
return p0, p1
def analyze_fp(x):
p0= x.sum(spaces=1)/fp_1.sum()
p1= x.sum(spaces=0)/fp_0.sum()
return p0, p1
However it seems to be not correct. Analysing the power of sk with test_product_space_power(fp)
in direction of space 0 seems to have a constant offset, while the power in direction of space 1 seems to have the wrong shape. analyze_fp(fp)
only serves as a double check to see if the covariance from which sk is actually drawn is the appropriate one for fp_0 and fp_1.
Am I using .power_analyze
wrong or are we facing an issue when we draw the random fields.
Thanks a lot for your help!