-
Notifications
You must be signed in to change notification settings - Fork 0
/
sim_uniform_pareto.py
134 lines (118 loc) · 4.28 KB
/
sim_uniform_pareto.py
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
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
import numpy as np
import typer
from lib.pareto_dist_bounded import sample
from lib.psd import get_snorp_psd
from lib.theory_psd import get_uniform_bounded_pareto_psd
def main(
repeats: int = 1,
n_events: int = 10000,
pulse_magnitude: float = 1,
min_pulse: float = 0,
max_pulse: float = 1e3,
min_gap: float = 1,
max_gap: float = 1e4,
power_gap: float = 1,
min_freq: float = -1,
max_freq: float = -1,
n_freq: int = 100,
archive_dir: str = "data",
seed: int = -1,
) -> None:
"""Simulate rectangular SNORP with uniform pulses and bounded Pareto gaps.
Input:
repeats: (default: 1)
Number of times to generate SNORP. Resulting PSD
will be averaged over all runs.
n_events: (default: 10000)
Number of pulses to simulate.
pulse_magnitude: (default: 1)
Fixed magnitude of the pulses in the signal.
min_pulse: (default: 0)
Minimum pulse duration.
max_pulse: (default: 1e3)
Maximum pulse duration.
min_gap: (default: 1)
Minimum gap duration.
max_gap: (default: 1e4)
Maximum gap duration.
power_gap: (default: 1)
Power law exponent of gap duration distribution.
min_freq: (default: -1)
Minimum frequency to observe. If negative value is
passed (which is the default), then the minimum
frequency is calculated automatically based on
simulation parameters.
max_freq: (default: -1)
Maximum frequency to observe. If negative value is
passed (which is the default), then the maximum
frequency is calculated automatically based on
simulation parameters.
n_freq: (default: 100)
Number of frequencies to consider within the given
(or automatically selected) range. Includes end
points.
archive_dir (default: "data")
Folder in which to save output files.
seed: (default: -1)
RNG seed. If negative value is passed (which is the
default), then seed will be randomly generated by
`np.random.randint(0, int(2**20))`
Output:
Function returns nothing, but saves one file, which
contains the numerically calculated PSD and its
theoretical estimate.
"""
# auto-generate seed
if seed < 0:
np.random.seed()
seed = np.random.randint(0, int(2**20))
# RNG setup
rng = np.random.default_rng(seed)
# sample distributions
def sample_pulse(low: float = 0, high: float = 1000, size: int = 1) -> np.ndarray:
return rng.uniform(low=low, high=high, size=size)
def sample_gap(
power: float, low: float = 1, high: float = 1000, size: int = 1
) -> np.ndarray:
return sample(power, low=low, high=high, size=size, rng=rng)
# simulation archival setup
model_info = f"unif{min_pulse*10000:.0f}_{max_pulse*10000:.0f}.pareto{power_gap*100:.0f}_{min_gap*1000:.0f}_{max_gap:.0f}"
simulation_filename = f"{model_info}.seed{seed:d}"
psd_path = f"{archive_dir}/{simulation_filename}.psd.csv"
# set frequency range
if max_freq < 0:
max_freq = (1 / max_gap) * 0.1 / (2 * np.pi)
if min_freq < 0:
min_freq = (1 / min_gap) * 10 / (2 * np.pi)
freqs = np.logspace(np.log10(min_freq), np.log10(max_freq), n_freq)
# main simulation loop
sim_psds = np.zeros((repeats, n_freq))
for sim_idx in range(repeats):
pulse_duration = sample_pulse(low=min_pulse, high=max_pulse, size=n_events)
gap_duration = sample_gap(power_gap, low=min_gap, high=max_gap, size=n_events)
sim_psds[sim_idx, :] = get_snorp_psd(
freqs,
pulse_duration,
gap_duration,
pulse_magnitude=pulse_magnitude,
)
# numerical PSD
sim_psd = np.mean(sim_psds, axis=0)
# theoretical PSD
theory_psd = get_uniform_bounded_pareto_psd(
freqs,
pulse_magnitude,
min_pulse,
max_pulse,
min_gap,
max_gap,
power_gap,
)
np.savetxt(
psd_path,
np.log10(np.vstack((freqs, sim_psd, theory_psd)).T),
delimiter=",",
fmt="%.4f",
)
if __name__ == "__main__":
typer.run(main)