multi_amplitudes_consistency.py 3.08 KB
Newer Older
1
2
import nifty5 as ift
import numpy as np
Philipp Arras's avatar
Philipp Arras committed
3
4
5
6
7
8
9
10
11
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

np.random.seed(42)


def testAmplitudesConsistency(seed, sspace):
    offset_std = 40
    intergated_fluct_std0 = 10.
    intergated_fluct_std1 = 2.

    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)

Philipp Arras's avatar
Minor    
Philipp Arras committed
68
69
70
71
    np.testing.assert_allclose(offset_std, zm_std_mean, rtol=0.5)
    np.testing.assert_allclose(intergated_fluct_std0, fluct_space, rtol=0.5)
    np.testing.assert_allclose(tot_flm, fluct_total, rtol=0.5)

Philipp Arras's avatar
Philipp Arras committed
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    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))


# Move to tests
# FIXME This test fails but it is not relevant for the final result
# assert_allclose(ampl(from_random('normal', ampl.domain)).val[0], vol) or 1??
# End move to tests

# move to tests
# assert_allclose(
#     smooth(from_random('normal', smooth.domain)).val[0:2], 0)
# end move to tests
for seed in [1, 42]:
    for sp in [
            ift.RGSpace((32, 64), (1.1, 0.3)),
            ift.HPSpace(64),
            ift.GLSpace(64)
    ]:
        testAmplitudesConsistency(seed, sp)