1d_separation.py 2.94 KB
Newer Older
1 2 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
from point_separation import build_problem, problem_iteration
from nifty2go import *
import numpy as np
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
from matplotlib import pyplot as plt

np.random.seed(42)
if __name__ == '__main__':
    s_space = RGSpace([1024])
    FFT = FFTOperator(s_space)
    h_space = FFT.target[0]
    p_space = PowerSpace(h_space)
    sp = Field(p_space, val=1./(1+p_space.k_lengths)**2.5 )
    sh = power_synthesize(sp)
    s = FFT.adjoint_times(sh)

    u = Field(s_space, val = -12)
    u.val[200] = 1
    u.val[300] = 3
    u.val[500] = 4
    u.val[700] = 5
    u.val[900] = 2
    u.val[154] = 0.5
    u.val[421] = 0.25
    u.val[652] = 1
    u.val[1002] = 2.5

    d = exp(s) + exp(u)
    data = d.val

    energy1 = build_problem(data,1.25)
    energy2 = build_problem(data,1.5)
    energy3 = build_problem(data,1.75)

    for i in range(20):
        energy1 = problem_iteration(energy1)
        energy2 = problem_iteration(energy2)
        energy3 = problem_iteration(energy3)

Jakob Knollmueller's avatar
Jakob Knollmueller committed
42
    size = 15
43 44 45
    plt.figure()
    # plt.plot(data, 'k-')
    f, (ax0, ax1,ax2) = plt.subplots(3, sharex=True, sharey=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
46
    plt.suptitle('diffuse components', size=size)
47 48 49

    ax0.plot(exp(energy1.s).val, 'k-')
    ax0.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
50
    ax0.set_ylabel(r'$\alpha = 1.25$', size=size)
51 52 53 54 55
    ax0.set_ylim(1e-1,1e3)
    ax0.set_yscale("log")

    ax1.plot(exp(energy2.s).val, 'k-')
    ax1.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
56
    ax1.set_ylabel(r'$\alpha = 1.5$', size=size)
57 58 59

    ax2.plot(exp(energy3.s).val, 'k-')
    ax2.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
60
    ax2.set_ylabel(r'$\alpha = 1.75$', size=size)
61 62 63 64 65 66

    plt.savefig('1d_diffuse.pdf')

    plt.figure()
    f, (ax0, ax1,ax2) = plt.subplots(3, sharex=True, sharey=True)

Jakob Knollmueller's avatar
Jakob Knollmueller committed
67
    plt.suptitle('point-like components', size=size)
68 69 70

    ax0.plot(exp(energy1.u).val, 'k-')
    ax0.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
71
    ax0.set_ylabel(r'$\alpha = 1.25$', size=size)
72 73 74 75 76
    ax0.set_ylim(1e-1,1e3)
    ax0.set_yscale("log")

    ax1.plot(exp(energy2.u).val, 'k-')
    ax1.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
77
    ax1.set_ylabel(r'$\alpha = 1.5$', size=size)
78 79 80

    ax2.plot(exp(energy3.u).val, 'k-')
    ax2.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
81
    ax2.set_ylabel(r'$\alpha = 1.75$', size=size)
82 83 84 85 86 87 88 89 90

    ax0.set_yscale("log")

    ax0.set_ylim(1e-1,1e3)

    # plt.ylim(1e-0)
    plt.savefig('1d_points.pdf')
    plt.figure()
    f, (ax0, ax1,ax2) = plt.subplots(3, sharex=True, sharey=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
91
    plt.suptitle('data and true components', size=size)
92 93 94 95 96

    ax0.plot(data, 'k-')
    ax0.set_yscale("log")
    ax0.set_ylim(1e-1,1e3)
    ax0.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
97
    ax0.set_ylabel(r'data', size=size)
98 99 100 101


    ax1.plot(exp(s).val, 'k-')
    ax1.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
102
    ax1.set_ylabel(r'diffuse', size=size)
103 104
    ax2.plot(exp(u).val, 'k-')
    ax2.yaxis.set_label_position("right")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
105
    ax2.set_ylabel(r'point-like', size=size)
106 107 108 109

    # plt.ylim(1e-0)
    plt.savefig('1d_data.pdf')