multi_amplitudes_consistency.py 2.21 KB
Newer Older
1
2
3
4
5
6
7
8
import nifty5 as ift
import numpy as np
#np.random.seed(42)

offset_std = 40
intergated_fluct_std0 = 10.
intergated_fluct_std1 = 2.

Philipp Frank's avatar
Philipp Frank committed
9
10
11
#sspace = ift.RGSpace((32,64),(1.1,0.3))
sspace = ift.HPSpace(64)
sspace = ift.GLSpace(64)
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
hspace = sspace.get_default_codomain()
target0 = ift.PowerSpace(hspace)

fsspace = ift.RGSpace((12,),(0.4,))
fhspace = fsspace.get_default_codomain()
target1 = ift.PowerSpace(fhspace)

fa = ift.CorrelatedFieldMaker()
fa.add_fluctuations(target0, intergated_fluct_std0, 1E-8,
                    1.1, 2., 2.1, .5, -2, 1., 'spatial')
fa.add_fluctuations(target1, intergated_fluct_std1, 1E-8,
                    3.1, 1., .5, .1, -4, 1., 'freq')
op = fa.finalize(offset_std, 1E-8, '')

flucts = [intergated_fluct_std0, intergated_fluct_std1]
tot_flm,totflsig = fa.effective_total_fluctuation(flucts,[1E-8,1E-8])

space = op.target
totaltoalvol = 1.
for s in space[:]:
    totaltoalvol *= s.total_volume

nsam = 1000

zm_std_mean = 0.
fluct_space = 0.
fluct_freq = 0.
fluct_total = 0.

for i in range(nsam):
    x = ift.from_random('normal',op.domain)
    res = op(x)

    zm = res.integrate()/totaltoalvol
    zm2 = res.mean()
    
    fl = ((res-zm)**2).integrate() / totaltoalvol
    
    zm_std_mean += zm**2
    fluct_total += fl

    r = res.integrate(1)/fsspace.total_volume
    r0 = r.integrate()/sspace.total_volume
    tm = ((r-r0)**2).integrate() / sspace.total_volume
    fluct_space += tm
    
    fr = res.integrate(0)/sspace.total_volume
    fr0 = fr.integrate()/fsspace.total_volume
    ftm = ((fr-fr0)**2).integrate() / fsspace.total_volume
    fluct_freq += ftm
    

fluct_total = np.sqrt(fluct_total/nsam)
fluct_space = np.sqrt(fluct_space/nsam)
fluct_freq = np.sqrt(fluct_freq/nsam)
zm_std_mean = np.sqrt(zm_std_mean/nsam)

print("Expected  offset Std: "+str(offset_std))
print("Estimated offset Std: "+str(zm_std_mean))

print("Expected  integrated fluct. space Std: "+str(intergated_fluct_std0))
print("Estimated integrated fluct. space Std: "+str(fluct_space))

print("Expected  integrated fluct. frequency Std: "+str(intergated_fluct_std1))
print("Estimated integrated fluct. frequency Std: "+str(fluct_freq))

print("Expected  total fluct. Std: "+str(tot_flm))
print("Estimated total fluct. Std: "+str(fluct_total))