semicoherent_glitch_robust_directed_grid_search_on_1_glitch.py 2.49 KB
Newer Older
Gregory Ashton's avatar
Gregory Ashton committed
1
2
3
import pyfstat
import numpy as np
import matplotlib.pyplot as plt
4
5
6
7
8
9
10
11
12
13
14
15
16
from make_simulated_data import (
    tstart,
    duration,
    tref,
    F0,
    F1,
    F2,
    Alpha,
    Delta,
    delta_F0,
    outdir,
    dtglitch,
)
Gregory Ashton's avatar
Gregory Ashton committed
17
import time
Gregory Ashton's avatar
Gregory Ashton committed
18
19
20
21
22
23

try:
    from gridcorner import gridcorner
except ImportError:
    raise ImportError(
        "Python module 'gridcorner' not found, please install from "
24
25
        "https://gitlab.aei.uni-hannover.de/GregAshton/gridcorner"
    )
Gregory Ashton's avatar
Gregory Ashton committed
26

27
label = "semicoherent_glitch_robust_directed_grid_search_on_1_glitch"
Gregory Ashton's avatar
Gregory Ashton committed
28

29
plt.style.use("./paper.mplstyle")
Gregory Ashton's avatar
Gregory Ashton committed
30

31
Nstar = 1000
32
33
F0_width = np.sqrt(Nstar) * np.sqrt(12) / (np.pi * duration)
F1_width = np.sqrt(Nstar) * np.sqrt(180) / (np.pi * duration ** 2)
Gregory Ashton's avatar
Gregory Ashton committed
34
N = 20
35
36
F0s = [F0 - F0_width / 2.0, F0 + F0_width / 2.0, F0_width / N]
F1s = [F1 - F1_width / 2.0, F1 + F1_width / 2.0, F1_width / N]
Gregory Ashton's avatar
Gregory Ashton committed
37
38
39
40
41
F2s = [F2]
Alphas = [Alpha]
Deltas = [Delta]

max_delta_F0 = 1e-5
42
43
tglitchs = [tstart + 0.1 * duration, tstart + 0.9 * duration, 0.8 * float(duration) / N]
delta_F0s = [0, max_delta_F0, max_delta_F0 / N]
Gregory Ashton's avatar
Gregory Ashton committed
44
45
delta_F1s = [0]

46

Gregory Ashton's avatar
Gregory Ashton committed
47
t1 = time.time()
Gregory Ashton's avatar
Gregory Ashton committed
48
search = pyfstat.GridGlitchSearch(
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
    label,
    outdir,
    "data/*1_glitch*sft",
    F0s=F0s,
    F1s=F1s,
    F2s=F2s,
    Alphas=Alphas,
    Deltas=Deltas,
    tref=tref,
    minStartTime=tstart,
    maxStartTime=tstart + duration,
    tglitchs=tglitchs,
    delta_F0s=delta_F0s,
    delta_F1s=delta_F1s,
)
Gregory Ashton's avatar
Gregory Ashton committed
64
search.run()
Gregory Ashton's avatar
Gregory Ashton committed
65
dT = time.time() - t1
Gregory Ashton's avatar
Gregory Ashton committed
66
67
68

F0_vals = np.unique(search.data[:, 0]) - F0
F1_vals = np.unique(search.data[:, 1]) - F1
69
delta_F0s_vals = np.unique(search.data[:, 5]) - delta_F0
Gregory Ashton's avatar
Gregory Ashton committed
70
tglitch_vals = np.unique(search.data[:, 7])
71
tglitch_vals_days = (tglitch_vals - tstart) / 86400.0 - dtglitch / 86400.0
Gregory Ashton's avatar
Gregory Ashton committed
72

73
74
75
76
77
78
79
80
81
82
83
twoF = search.data[:, -1].reshape(
    (len(F0_vals), len(F1_vals), len(delta_F0s_vals), len(tglitch_vals))
)
xyz = [F0_vals * 1e6, F1_vals * 1e12, delta_F0s_vals * 1e6, tglitch_vals_days]
labels = [
    "$f - f_\mathrm{s}$\n[$\mu$Hz]",
    "$\dot{f} - \dot{f}_\mathrm{s}$\n[$p$Hz/s]",
    "$\delta f-\delta f_\mathrm{s}$\n[$\mu$Hz]",
    "$t^\mathrm{g} - t^\mathrm{g}_\mathrm{s}$\n[d]",
    "$\widehat{2\mathcal{F}}$",
]
Gregory Ashton's avatar
Gregory Ashton committed
84
fig, axes = gridcorner(
85
86
87
88
89
90
91
92
93
94
    twoF,
    xyz,
    projection="log_mean",
    labels=labels,
    showDvals=False,
    lines=[0, 0, 0, 0],
    label_offset=0.25,
    max_n_ticks=4,
)
fig.savefig("{}/{}_projection_matrix.png".format(outdir, label), bbox_inches="tight")
Gregory Ashton's avatar
Gregory Ashton committed
95
96


97
print(("Prior widths =", F0_width, F1_width))
98
99
print(("Actual run time = {}".format(dT)))
print(("Actual number of grid points = {}".format(search.data.shape[0])))