1.. _example_beam_error_spot_frequencies:
2
3Beam errors at spot frequencies
4===============================
5
6Summary
7-------
8
9This script demonstrates how the OSKAR Python interface can be used to
10run simulations which compare the performance of a set of telescope models.
11
12The idea was to investigate how errors in the station beams would affect the
13quality of images produced by the SKA-Low telescope at a range of spot
14frequencies across the band, at 50, 70, 110, 137, 160, 230 and 320 MHz.
15A number of "corrupted" visibility data sets were generated at each frequency
16by introducing progressively larger errors in the station beams, and these were
17subtracted from a reference "perfect" data set generated without any
18corruptions. Dirty images were then made of the residual visibilities to
19estimate the amount of noise introduced by the beam errors.
20
21The station beam errors were generated by varying the complex beamforming
22weights randomly within each of the 512 stations over a defined range.
23These errors were kept static throughout each observation to simulate the
24effects of residual calibration errors at the station level.
25Element phase errors were calculated as a function of frequency using a
26Gaussian-distributed cable length tolerance of 15 mm.
27Element gain errors were generated from Gaussian distributions with
28widths varying between 0.005 dB and 0.640 dB.
29
30The sky models and other observation parameters were kept the same for each
31of the data sets. Three different fields were used, with each populated by
32sources from the GLEAM extragalactic catalogue, and the 10 brightest "A-team"
33sources in the sky were also added.
34All observations were set up to span 4 hours and to be symmetric about the
35meridian, with station beams and visibilities evaluated every 60 seconds to
36generate a representative aperture-plane coverage.
37
38These results were written up and presented as part of the SKA System CDR
39documentation pack.
40
41The simulation script
42---------------------
43
44The script starts (in ``main``) by declaring the axes of the parameter space
45to iterate over: the seven frequencies of interest, the eight gain error
46values, and the three target fields. The interesting parts of the script are
47mainly in three functions:
48
49- The ``run_set`` function uses three nested loops to map out the parameter
50  space, steering the script by creating the :class:`oskar.Sky` model for
51  each field, updating a :class:`oskar.Telescope` model, and setting the
52  relevant parameters in a :class:`oskar.SettingsTree` before simulating the
53  visibilities using :class:`oskar.Interferometer` (in ``make_vis_data``).
54
55- The ``make_diff_image_stats`` function generates the residual visibilities
56  by taking the difference of the data in two :class:`oskar.VisBlock` classes
57  (the "perfect" and "corrupted" data), which are read from two files.
58  Note that a ``ThreadPoolExecutor`` is used to read data from the two files
59  concurrently to improve performance. The residual visibilities are generated
60  using ``block1.cross_correlations()[...] -= block2.cross_correlations()``
61  and these are passed directly (still inside ``block1``)
62  to :class:`oskar.Imager` via its
63  :meth:`update_from_block() <oskar.Imager.update_from_block()>` method.
64  The images are then returned straight back for analysis using functions
65  from ``numpy``, and the image statistics are stored in a Python dictionary
66  called ``results``.
67
68- Finally, the results dictionary is plotted using ``matplotlib`` in
69  the ``make_plot`` function.
70
71
72.. code-block:: python
73
74    #!/usr/bin/env python3
75
76    """Script to run LOW station beam error simulations at spot frequencies."""
77
78    from __future__ import print_function, division
79    import concurrent.futures
80    import json
81    import logging
82    import os
83    import sys
84
85    from astropy.io import fits
86    from astropy.time import Time, TimeDelta
87    import matplotlib
88    matplotlib.use('Agg')
89    # pylint: disable=wrong-import-position
90    import matplotlib.pyplot as plt
91    import numpy
92    import oskar
93
94
95    LOG = logging.getLogger()
96
97
98    def bright_sources():
99        """Returns a list of bright A-team sources."""
100        # Sgr A: guesstimates only!
101        # For A: data from the Molonglo Southern 4 Jy sample (VizieR).
102        # Others from GLEAM reference paper, Hurley-Walker et al. (2017), Table 2.
103        return numpy.array((
104            [266.41683, -29.00781,  2000,0,0,0,   0,    0,    0, 3600, 3600, 0],
105            [ 50.67375, -37.20833,   528,0,0,0, 178e6, -0.51, 0, 0, 0, 0],  # For
106            [201.36667, -43.01917,  1370,0,0,0, 200e6, -0.50, 0, 0, 0, 0],  # Cen
107            [139.52500, -12.09556,   280,0,0,0, 200e6, -0.96, 0, 0, 0, 0],  # Hyd
108            [ 79.95833, -45.77889,   390,0,0,0, 200e6, -0.99, 0, 0, 0, 0],  # Pic
109            [252.78333,   4.99250,   377,0,0,0, 200e6, -1.07, 0, 0, 0, 0],  # Her
110            [187.70417,  12.39111,   861,0,0,0, 200e6, -0.86, 0, 0, 0, 0],  # Vir
111            [ 83.63333,  22.01444,  1340,0,0,0, 200e6, -0.22, 0, 0, 0, 0],  # Tau
112            [299.86667,  40.73389,  7920,0,0,0, 200e6, -0.78, 0, 0, 0, 0],  # Cyg
113            [350.86667,  58.81167, 11900,0,0,0, 200e6, -0.41, 0, 0, 0, 0]   # Cas
114            ))
115
116
117    def get_start_time(ra0_deg, length_sec):
118        """Returns optimal start time for field RA and observation length."""
119        t = Time('2000-01-01 00:00:00', scale='utc', location=('116.764d', '0d'))
120        dt_hours = 24.0 - t.sidereal_time('apparent').hour + (ra0_deg / 15.0)
121        start = t + TimeDelta(dt_hours * 3600.0 - length_sec / 2.0, format='sec')
122        return start.value
123
124
125    def make_vis_data(settings, sky, tel):
126        """Run simulation using supplied settings."""
127        if os.path.exists(settings['interferometer/oskar_vis_filename']):
128            LOG.info("Skipping simulation, as output data already exist.")
129            return
130        LOG.info("Simulating %s", settings['interferometer/oskar_vis_filename'])
131        sim = oskar.Interferometer(settings=settings)
132        sim.set_sky_model(sky)
133        sim.set_telescope_model(tel)
134        sim.run()
135
136
137    def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy):
138        """Filter sky model.
139
140        Includes all sources within the given radius, and sources above the
141        specified flux outside this radius.
142        """
143        # Get pointing centre.
144        ra0_deg = float(settings['observation/phase_centre_ra_deg'])
145        dec0_deg = float(settings['observation/phase_centre_dec_deg'])
146
147        # Create "inner" and "outer" sky models.
148        sky_inner = sky0.create_copy()
149        sky_outer = sky0.create_copy()
150        sky_inner.filter_by_radius(0.0, radius_deg, ra0_deg, dec0_deg)
151        sky_outer.filter_by_radius(radius_deg, 180.0, ra0_deg, dec0_deg)
152        sky_outer.filter_by_flux(flux_min_outer_jy, 1e9)
153        LOG.info("Number of sources in sky0: %d", sky0.num_sources)
154        LOG.info("Number of sources in inner sky model: %d", sky_inner.num_sources)
155        LOG.info("Number of sources in outer sky model above %.3f Jy: %d",
156                 flux_min_outer_jy, sky_outer.num_sources)
157        sky_outer.append(sky_inner)
158        LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources)
159        return sky_outer
160
161
162    def make_diff_image_stats(filename1, filename2, use_w_projection,
163                              out_image_root=None):
164        """Make an image of the difference between two visibility data sets.
165
166        This function assumes that the observation parameters for both data sets
167        are identical. (It will fail horribly otherwise!)
168        """
169        # Set up an imager.
170        (hdr1, handle1) = oskar.VisHeader.read(filename1)
171        (hdr2, handle2) = oskar.VisHeader.read(filename2)
172        frequency_hz = hdr1.freq_start_hz
173        fov_ref_frequency_hz = 140e6
174        fov_ref_deg = 5.0
175        fov_deg = fov_ref_deg * (fov_ref_frequency_hz / frequency_hz)
176        imager = oskar.Imager(precision='double')
177        imager.set(fov_deg=fov_deg, image_size=8192,
178                   fft_on_gpu=True, grid_on_gpu=True)
179        if out_image_root is not None:
180            imager.output_root = out_image_root
181
182        LOG.info("Imaging differences between '%s' and '%s'", filename1, filename2)
183        block1 = oskar.VisBlock.create_from_header(hdr1)
184        block2 = oskar.VisBlock.create_from_header(hdr2)
185        if hdr1.num_blocks != hdr2.num_blocks:
186            raise RuntimeError("'%s' and '%s' have different dimensions!" %
187                               (filename1, filename2))
188        if use_w_projection:
189            imager.set(algorithm='W-projection')
190            imager.coords_only = True
191            for i_block in range(hdr1.num_blocks):
192                block1.read(hdr1, handle1, i_block)
193                imager.update_from_block(hdr1, block1)
194            imager.coords_only = False
195            imager.check_init()
196            LOG.info("Using %d W-planes", imager.num_w_planes)
197        executor = concurrent.futures.ThreadPoolExecutor(2)
198        for i_block in range(hdr1.num_blocks):
199            tasks_read = []
200            tasks_read.append(executor.submit(block1.read, hdr1, handle1, i_block))
201            tasks_read.append(executor.submit(block2.read, hdr2, handle2, i_block))
202            concurrent.futures.wait(tasks_read)
203            block1.cross_correlations()[...] -= block2.cross_correlations()
204            imager.update_from_block(hdr1, block1)
205        del handle1, handle2, hdr1, hdr2, block1, block2
206
207        # Finalise image and return it to Python.
208        output = imager.finalise(return_images=1)
209        image = output['images'][0]
210
211        LOG.info("Generating image statistics")
212        image_size = imager.image_size
213        box_size = int(0.1 * image_size)
214        centre = image[
215            (image_size - box_size)//2:(image_size + box_size)//2,
216            (image_size - box_size)//2:(image_size + box_size)//2]
217        del imager
218        return {
219            'image_medianabs': numpy.median(numpy.abs(image)),
220            'image_mean': numpy.mean(image),
221            'image_std': numpy.std(image),
222            'image_rms': numpy.sqrt(numpy.mean(image**2)),
223            'image_centre_mean': numpy.mean(centre),
224            'image_centre_std': numpy.std(centre),
225            'image_centre_rms': numpy.sqrt(numpy.mean(centre**2))
226        }
227
228
229    def make_plot(prefix, field_name, metric_key, results,
230                  axis_freq, axis_gain):
231        """Plot selected results."""
232        # Get data for contour plot.
233        X, Y = numpy.meshgrid(axis_freq, axis_gain)
234        Z = numpy.zeros(X.shape)
235        for freq, gain, z in numpy.nditer([X, Y, Z], op_flags=['readwrite']):
236            key = '%s_%s_%d_MHz_%.3f_dB' % (prefix, field_name, freq, gain)
237            if key in results:
238                z[...] = numpy.log10(results[key][metric_key])
239        ax1 = plt.subplot(111)
240        ax1.set_yscale('log')
241
242        # Scatter plot.
243        sp = ax1.scatter(X, Y, c=Z, cmap='plasma')
244
245        # Contour plot.
246        cp = ax1.contour(X, Y, Z, cmap='plasma')
247        levels = cp.levels
248        print(prefix, field_name, len(levels))
249        if len(levels) > 9:
250            levels = levels[::2]
251        clabels = plt.clabel(cp, levels, inline=False, fontsize=10, fmt='%1.1f')
252        for txt in clabels:
253            txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=1))
254
255        # Title and axis labels.
256        metric_name = '[ UNKNOWN ]'
257        if metric_key == 'image_centre_rms':
258            metric_name = 'Central RMS [Jy/beam]'
259        elif metric_key == 'image_medianabs':
260            metric_name = 'MEDIAN(ABS(image)) [Jy/beam]'
261        sky_model = 'GLEAM'
262        if 'A-team' in prefix:
263            sky_model = sky_model + ' + A-team'
264        plt.title('%s for %s field (%s)' % (metric_name, field_name, sky_model))
265        plt.xlabel('Frequency [MHz]')
266        plt.ylabel('Element gain standard deviation [dB]')
267        cbar = plt.colorbar(sp)
268        cbar.set_label('log10(%s)' % metric_name)
269        plt.savefig('%s_%s_%s.png' % (prefix, field_name, metric_key))
270        plt.close('all')
271
272
273    def run_single(prefix_field, settings, sky, tel,
274                   freq_MHz, gain_std_dB, out0_name, results):
275        """Run a single simulation and generate image statistics for it."""
276        out = '%s_%d_MHz_%.3f_dB' % (prefix_field, freq_MHz, gain_std_dB)
277        if out in results:
278            LOG.info("Using cached results for '%s'", out)
279            return
280        out_name = out + '.vis'
281        gain_std = numpy.power(10.0, gain_std_dB / 20.0) - 1.0
282        tel.override_element_gains(1.0, gain_std)
283        tel.override_element_cable_length_errors(0.015)
284        settings['interferometer/oskar_vis_filename'] = out_name
285        make_vis_data(settings, sky, tel)
286        out_image_root = out
287        use_w_projection = True
288        if str(settings['interferometer/ignore_w_components']).lower() == 'true':
289            use_w_projection = False
290        results[out] = make_diff_image_stats(out0_name, out_name, use_w_projection,
291                                             out_image_root)
292
293
294    def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only):
295        """Runs a set of simulations."""
296        if not plot_only:
297            # Load base telescope model.
298            settings = oskar.SettingsTree('oskar_sim_interferometer')
299            settings.from_dict(base_settings)
300            tel = oskar.Telescope(settings=settings)
301
302            # Load base sky model
303            sky0 = oskar.Sky()
304            if 'GLEAM' in prefix:
305                # Load GLEAM catalogue from FITS binary table.
306                hdulist = fits.open('GLEAM_EGC.fits')
307                # pylint: disable=no-member
308                cols = hdulist[1].data[0].array
309                data = numpy.column_stack(
310                    (cols['RAJ2000'], cols['DEJ2000'], cols['peak_flux_wide']))
311                data = data[data[:, 2].argsort()[::-1]]
312                sky_gleam = oskar.Sky.from_array(data)
313                sky0.append(sky_gleam)
314            if 'A-team' in prefix:
315                sky_bright = oskar.Sky.from_array(bright_sources())
316                sky0.append(sky_bright)
317
318        # Iterate over fields.
319        for field_name, field in fields.items():
320            # Load result set, if it exists.
321            prefix_field = prefix + '_' + field_name
322            results = {}
323            json_file = prefix_field + '_results.json'
324            if os.path.exists(json_file):
325                with open(json_file, 'r') as input_file:
326                    results = json.load(input_file)
327
328            # Iterate over frequencies.
329            if not plot_only:
330                for freq_MHz in axis_freq:
331                    # Update settings for field.
332                    settings_dict = base_settings.copy()
333                    settings_dict.update(field)
334                    settings.from_dict(settings_dict)
335                    ra_deg = float(settings['observation/phase_centre_ra_deg'])
336                    dec_deg = float(settings['observation/phase_centre_dec_deg'])
337                    length_sec = float(settings['observation/length'])
338                    settings['observation/start_frequency_hz'] = str(freq_MHz * 1e6)
339                    settings['observation/start_time_utc'] = get_start_time(
340                        ra_deg, length_sec)
341                    tel.set_phase_centre(ra_deg, dec_deg)
342
343                    # Create the sky model.
344                    sky = make_sky_model(sky0, settings, 20.0, 10.0)
345                    settings['interferometer/ignore_w_components'] = 'true'
346                    if 'A-team' in prefix:
347                        settings['interferometer/ignore_w_components'] = 'false'
348
349                    # Simulate the 'perfect' case.
350                    tel.override_element_gains(1.0, 0.0)
351                    tel.override_element_cable_length_errors(0.0)
352                    out0_name = '%s_%d_MHz_no_errors.vis' % (prefix_field, freq_MHz)
353                    settings['interferometer/oskar_vis_filename'] = out0_name
354                    make_vis_data(settings, sky, tel)
355
356                    # Simulate the error cases.
357                    for gain_std_dB in axis_gain:
358                        run_single(prefix_field, settings, sky, tel,
359                                   freq_MHz, gain_std_dB, out0_name, results)
360
361            # Generate plot for the field.
362            make_plot(prefix, field_name, 'image_centre_rms',
363                      results, axis_freq, axis_gain)
364            make_plot(prefix, field_name, 'image_medianabs',
365                      results, axis_freq, axis_gain)
366
367            # Save result set.
368            with open(json_file, 'w') as output_file:
369                json.dump(results, output_file, indent=4)
370
371
372    def main():
373        """Main function."""
374        handler = logging.StreamHandler(sys.stdout)
375        formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
376        handler.setFormatter(formatter)
377        LOG.addHandler(handler)
378        LOG.setLevel(logging.INFO)
379
380        # Define common settings.
381        base_settings = {
382            'simulator': {
383                'double_precision': 'true',
384                'use_gpus': 'true',
385                'max_sources_per_chunk': '23000'
386            },
387            'observation' : {
388                'frequency_inc_hz': '100e3',
389                'length': '14400.0',
390                'num_time_steps': '240'
391            },
392            'telescope': {
393                'input_directory': 'SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm'
394            },
395            'interferometer': {
396                'channel_bandwidth_hz': '100e3',
397                'time_average_sec': '1.0',
398                'max_time_samples_per_block': '4'
399            }
400        }
401
402        # Define axes of parameter space.
403        fields = {
404            'EoR0': {
405                'observation/phase_centre_ra_deg': '0.0',
406                'observation/phase_centre_dec_deg': '-27.0'
407            },
408            'EoR1': {
409                'observation/phase_centre_ra_deg': '60.0',
410                'observation/phase_centre_dec_deg': '-30.0'
411            },
412            'EoR2': {
413                'observation/phase_centre_ra_deg': '170.0',
414                'observation/phase_centre_dec_deg': '-10.0'
415            }
416        }
417        axis_freq = [50, 70, 110, 137, 160, 230, 320]
418        axis_gain = [0.005, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64]
419
420        # GLEAM + A-team sky model simulations.
421        plot_only = False
422        run_set('GLEAM_A-team', base_settings,
423                fields, axis_freq, axis_gain, plot_only)
424
425
426    if __name__ == '__main__':
427        main()
428