multi_amplitudes_consistency.py 4 KB
Newer Older
Martin Reinecke's avatar
5->6    
Martin Reinecke committed
1
import nifty6 as ift
2
import numpy as np
Philipp Arras's avatar
Philipp Arras committed
3

4

Philipp Arras's avatar
Philipp Arras committed
5
6
7


def testAmplitudesConsistency(seed, sspace):
Philipp Frank's avatar
fixes    
Philipp Frank committed
8
9
10
11
12
    def stats(op,samples):
        sc = ift.StatCalculator()
        for s in samples:
            sc.add(op(s.extract(op.domain)))
        return sc.mean.to_global_data(), sc.var.sqrt().to_global_data()
13
    np.random.seed(seed)
14
    offset_std = .1
15
16
17
    intergated_fluct_std0 = .003
    intergated_fluct_std1 = 0.1
    
18
    nsam = 1000
Philipp Arras's avatar
Philipp Arras committed
19
20
21
22


    fsspace = ift.RGSpace((12,), (0.4,))

23
    fa = ift.CorrelatedFieldMaker.make(offset_std, 1E-8, '')
24
    fa.add_fluctuations(sspace, intergated_fluct_std0, 1E-8, 1.1, 2., 2.1, .5,
Philipp Arras's avatar
Philipp Arras committed
25
                        -2, 1., 'spatial')
26
    fa.add_fluctuations(fsspace, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1,
Philipp Arras's avatar
Philipp Arras committed
27
                        -4, 1., 'freq')
28
    op = fa.finalize()
Philipp Arras's avatar
Philipp Arras committed
29

30
    samples = [ift.from_random('normal',op.domain) for _ in range(nsam)]
Philipp Frank's avatar
fixes    
Philipp Frank committed
31
32
33
34
    tot_flm, _ = stats(fa.total_fluctuation,samples)
    offset_std,_ = stats(fa.amplitude_total_offset,samples)
    intergated_fluct_std0,_ = stats(fa.average_fluctuation(0),samples)
    intergated_fluct_std1,_ = stats(fa.average_fluctuation(1),samples)
35
    
Philipp Frank's avatar
fixes    
Philipp Frank committed
36
37
    slice_fluct_std0,_ = stats(fa.slice_fluctuation(0),samples)
    slice_fluct_std1,_ = stats(fa.slice_fluctuation(1),samples)
38
39
40
41
42
43
44
45

    sams = [op(s) for s in samples]
    fluct_total = fa.total_fluctuation_realized(sams)
    fluct_space = fa.average_fluctuation_realized(sams,0)
    fluct_freq = fa.average_fluctuation_realized(sams,1)
    zm_std_mean = fa.offset_amplitude_realized(sams)
    sl_fluct_space = fa.slice_fluctuation_realized(sams,0)
    sl_fluct_freq = fa.slice_fluctuation_realized(sams,1)
Philipp Arras's avatar
Philipp Arras committed
46
47
48
49
50
51
52
53
54
55
56

    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))
57
58
59
60
61
62
63
64
    
    print("Expected  slice fluct. space Std: " +
          str(slice_fluct_std0))
    print("Estimated slice fluct. space Std: " + str(sl_fluct_space))

    print("Expected  slice fluct. frequency Std: " +
          str(slice_fluct_std1))
    print("Estimated slice fluct. frequency Std: " + str(sl_fluct_freq))
Philipp Arras's avatar
Philipp Arras committed
65
66
67

    print("Expected  total fluct. Std: " + str(tot_flm))
    print("Estimated total fluct. Std: " + str(fluct_total))
68
    
69
70
71
72
73
74
75
76
77
78
    
    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(intergated_fluct_std1, fluct_freq, rtol=0.5)
    np.testing.assert_allclose(tot_flm, fluct_total, rtol=0.5)
    np.testing.assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5)
    np.testing.assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5)

    
    fa = ift.CorrelatedFieldMaker.make(offset_std, .1, '')
79
    fa.add_fluctuations(fsspace, intergated_fluct_std1, 1., 3.1, 1., .5, .1,
80
81
82
                        -4, 1., 'freq')
    m = 3.
    x = fa.moment_slice_to_average(m)
83
    fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5,
84
                        -2, 1., 'spatial', 0)
85
    op = fa.finalize()
Philipp Frank's avatar
fixes    
Philipp Frank committed
86
    em, estd = stats(fa.slice_fluctuation(0),samples)
87
88
89
    print("Forced   slice fluct. space Std: "+str(m))
    print("Expected slice fluct. Std: " + str(em))
    np.testing.assert_allclose(m, em, rtol=0.5)
90
91
92
    
    assert op.target[0] == sspace
    assert op.target[1] == fsspace
Philipp Arras's avatar
Philipp Arras committed
93
94
95
96
97
98
99
100
101
102
103
104
105
106


# 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)),
Philipp Frank's avatar
fixup    
Philipp Frank committed
107
108
            ift.HPSpace(32),
            ift.GLSpace(32)
Philipp Arras's avatar
Philipp Arras committed
109
110
    ]:
        testAmplitudesConsistency(seed, sp)