Defining a parameter space and running simulations
An investigation may require a large number of simulations to be carried out in order to explore a parameter space, which typically means that a set of nested loops must be written in order to run all the simulations.
In many cases, only the simulation parameters in the settings tree need to be changed to run a new simulation, but sometimes the sky model and/or telescope model also needs to be changed within a loop. Defining the parameters that need to vary is the first thing to do when writing a new simulation script.
Example: A four-dimensional parameter space
A simple example script which iterates over a 4-dimensional parameter space is shown below. In this case, the observation length (3 values: short, medium, long), the target field (3 values: EoR0, EoR1, EoR2), the ionospheric phase screen (2 values: on, off) and the sky model (2 values: GLEAM and A-team only) were all varied for a total of 36 simulations, and a CASA Measurement Set was written for each case.
Note how nested Python dictionaries are used to define groups of parameters
that need to change on each iteration, and the update() method is used to
merge one dictionary into another.
Each dimension is iterated using the general form:
# Iterate over a dimension.
for key_name, params_to_update in dictionary.items():
# Update the current settings dictionary.
current_settings.update(params_to_update)
# Iterate over the next dimension...
These are all standard Python constructs.
After all the parameters have been set up in the settings tree, an instance
of oskar.Interferometer is created using it. Finally,
calling oskar.Interferometer.run() will run each simulation.
1#!/usr/bin/env python3
2"""
3Run simulations for SKA1-LOW direction-dependent effects.
4https://confluence.skatelescope.org/display/SE/Simulations+with+Direction-Dependent+Effects
5https://jira.skatelescope.org/browse/SIM-489
6"""
7
8import copy
9import os.path
10
11from astropy.io import fits
12import numpy
13import oskar
14
15from .utils import get_start_time
16
17
18def bright_sources():
19 """
20 Returns a list of bright A-team sources.
21 Does not include the Galactic Centre!
22 """
23 # For A: data from the Molonglo Southern 4 Jy sample (VizieR).
24 # Others from GLEAM reference paper, Hurley-Walker et al. (2017), Table 2.
25 return numpy.array(
26 (
27 [
28 50.67375,
29 -37.20833,
30 528,
31 0,
32 0,
33 0,
34 178e6,
35 -0.51,
36 0,
37 0,
38 0,
39 0,
40 ], # For
41 [
42 201.36667,
43 -43.01917,
44 1370,
45 0,
46 0,
47 0,
48 200e6,
49 -0.50,
50 0,
51 0,
52 0,
53 0,
54 ], # Cen
55 [
56 139.52500,
57 -12.09556,
58 280,
59 0,
60 0,
61 0,
62 200e6,
63 -0.96,
64 0,
65 0,
66 0,
67 0,
68 ], # Hyd
69 [
70 79.95833,
71 -45.77889,
72 390,
73 0,
74 0,
75 0,
76 200e6,
77 -0.99,
78 0,
79 0,
80 0,
81 0,
82 ], # Pic
83 [
84 252.78333,
85 4.99250,
86 377,
87 0,
88 0,
89 0,
90 200e6,
91 -1.07,
92 0,
93 0,
94 0,
95 0,
96 ], # Her
97 [
98 187.70417,
99 12.39111,
100 861,
101 0,
102 0,
103 0,
104 200e6,
105 -0.86,
106 0,
107 0,
108 0,
109 0,
110 ], # Vir
111 [
112 83.63333,
113 22.01444,
114 1340,
115 0,
116 0,
117 0,
118 200e6,
119 -0.22,
120 0,
121 0,
122 0,
123 0,
124 ], # Tau
125 [
126 299.86667,
127 40.73389,
128 7920,
129 0,
130 0,
131 0,
132 200e6,
133 -0.78,
134 0,
135 0,
136 0,
137 0,
138 ], # Cyg
139 [
140 350.86667,
141 58.81167,
142 11900,
143 0,
144 0,
145 0,
146 200e6,
147 -0.41,
148 0,
149 0,
150 0,
151 0,
152 ], # Cas
153 )
154 )
155
156
157def main():
158 """Main function."""
159 # Load GLEAM catalogue data as a sky model.
160 sky_dir = "./"
161 gleam = fits.getdata(sky_dir + "GLEAM_EGC.fits", 1)
162 gleam_sky_array = numpy.column_stack(
163 (gleam["RAJ2000"], gleam["DEJ2000"], gleam["peak_flux_wide"])
164 )
165
166 # Define common base settings.
167 tel_dir = "./"
168 tel_model = "SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm"
169 common_settings = {
170 "simulator/max_sources_per_chunk": 65536,
171 "simulator/write_status_to_log_file": True,
172 "observation/start_frequency_hz": 125e6, # First channel at 125 MHz.
173 "observation/frequency_inc_hz": 5e6, # Channels spaced every 5 MHz.
174 "observation/num_channels": 11,
175 "telescope/input_directory": tel_dir + tel_model,
176 "interferometer/channel_bandwidth_hz": 100e3, # 100 kHz-wide channels.
177 "interferometer/time_average_sec": 1.0,
178 "interferometer/max_time_samples_per_block": 4,
179 }
180
181 # Define observations.
182 observations = {
183 "short": {
184 "observation/length": 5 * 60,
185 "observation/num_time_steps": 300,
186 "telescope/external_tec_screen/input_fits_file": "screen_short_300_1.0.fits",
187 },
188 "medium": {
189 "observation/length": 30 * 60,
190 "observation/num_time_steps": 300,
191 "telescope/external_tec_screen/input_fits_file": "screen_medium_300_6.0.fits",
192 },
193 "long": {
194 "observation/length": 4 * 60 * 60,
195 "observation/num_time_steps": 240,
196 "telescope/external_tec_screen/input_fits_file": "screen_long_240_60.0.fits",
197 },
198 }
199
200 # Define fields.
201 fields = {
202 "EoR0": {
203 "observation/phase_centre_ra_deg": 0.0,
204 "observation/phase_centre_dec_deg": -27.0,
205 },
206 "EoR1": {
207 "observation/phase_centre_ra_deg": 60.0,
208 "observation/phase_centre_dec_deg": -30.0,
209 },
210 "EoR2": {
211 "observation/phase_centre_ra_deg": 170.0,
212 "observation/phase_centre_dec_deg": -10.0,
213 },
214 }
215
216 # Define ionosphere settings.
217 ionosphere = {
218 "ionosphere_on": {"telescope/ionosphere_screen_type": "External"},
219 "ionosphere_off": {"telescope/ionosphere_screen_type": "None"},
220 }
221
222 # Define sky model components.
223 sky_models = {
224 "A-team": oskar.Sky.from_array(bright_sources()),
225 "GLEAM": oskar.Sky.from_array(gleam_sky_array),
226 }
227
228 # Loop over observations.
229 for obs_name, obs_params in observations.items():
230 # Copy the base settings dictionary.
231 current_settings = copy.deepcopy(common_settings)
232
233 # Update current settings with observation parameters.
234 current_settings.update(obs_params)
235
236 # Loop over fields.
237 for field_name, field_params in fields.items():
238 # Update current settings with field parameters.
239 current_settings.update(field_params)
240
241 # Update current settings with start time.
242 ra0_deg = current_settings["observation/phase_centre_ra_deg"]
243 length_sec = current_settings["observation/length"]
244 start_time = get_start_time(ra0_deg, length_sec)
245 current_settings["observation/start_time_utc"] = start_time
246
247 # Loop over ionospheric screen on/off.
248 for iono_name, iono_params in ionosphere.items():
249 # Update current settings with ionosphere parameters.
250 current_settings.update(iono_params)
251
252 # Loop over sky model components.
253 for sky_name, sky_model in sky_models.items():
254 # Update output MS filename based on current parameters.
255 ms_name = "SKA_LOW_SIM"
256 ms_name += "_" + obs_name
257 ms_name += "_" + field_name
258 ms_name += "_" + iono_name
259 ms_name += "_" + sky_name
260 ms_name += ".MS"
261
262 # Check if the MS already exists (if so, skip).
263 if os.path.isdir(ms_name):
264 continue
265
266 # Create the settings tree.
267 settings = oskar.SettingsTree("oskar_sim_interferometer")
268 settings.from_dict(current_settings)
269 settings["interferometer/ms_filename"] = ms_name
270
271 # Set up the simulator and run it.
272 sim = oskar.Interferometer(settings=settings)
273 sim.set_sky_model(sky_model)
274 sim.run()
275
276
277if __name__ == "__main__":
278 main()
Example: An irregular frequency axis
Frequency channels which are regularly spaced can be run in one go (and written to a single Measurement Set if required) by specifying multiple channels in the settings. However, when running simulations at spot frequencies across a band, these will need to be run separately by explicitly looping over each one. All that is required is to define a list of frequencies and then loop over them, for example:
axis_freq_MHz = [50, 70, 110, 137, 160, 230, 320]
for freq_MHz in axis_freq_MHz:
settings['observation/start_frequency_hz'] = freq_MHz * 1e6
...