Commits
Vincent Geers authored 331abc88d30 Merge
350 350 | # Select data for current polarisation. |
351 351 | ampdata = np.squeeze(data['corrected_data'], axis=1)[pol] |
352 352 | flagdata = np.squeeze(data['flag'], axis=1)[pol] |
353 353 | weightdata = data['weight'][pol] |
354 354 | |
355 355 | # Select for non-flagged data and non-NaN data. |
356 356 | id_nonbad = np.where(np.logical_and(np.logical_not(flagdata), np.isfinite(ampdata))) |
357 357 | amplitudes = ampdata[id_nonbad] |
358 358 | weights = weightdata[id_nonbad] |
359 359 | |
360 + | # Determine number of non-flagged antennas, baselines, and |
361 + | # integrations. |
362 + | ant1 = data['antenna1'][id_nonbad] |
363 + | ant2 = data['antenna2'][id_nonbad] |
364 + | n_ants = len(set(ant1) | set(ant2)) |
365 + | n_baselines = n_ants * (n_ants - 1) // 2 |
366 + | n_ints = len(amplitudes) // n_baselines |
367 + | |
368 + | # PIPE-644: Determine scale factor for variance. |
369 + | var_scale = np.mean([len(amplitudes), n_ints * n_ants]) |
370 + | |
360 371 | # Compute mean flux and stdev for current polarisation. |
361 372 | mean_flux = np.abs(np.average(amplitudes, weights=weights)) |
362 - | variance = np.average((np.abs(amplitudes) - mean_flux)**2, weights=weights) |
373 + | variance = np.average((np.abs(amplitudes) - mean_flux)**2, weights=weights) / var_scale |
363 374 | |
364 375 | # Store for this polarisation. |
365 376 | mean_fluxes.append(mean_flux) |
366 377 | variances.append(variance) |
367 378 | |
368 379 | # Compute mean flux and mean stdev for all polarisations. |
369 380 | mean_flux = np.mean(mean_fluxes) |
370 381 | mean_std = np.sqrt(np.mean(variances)) |
371 382 | |
372 383 | # Combine into single flux measurement object. |
584 595 | openms.msselect(data_selection) |
585 596 | |
586 597 | # Read in the weights from MS separately, prior to channel |
587 598 | # averaging. |
588 599 | weights = openms.getdata(['weight']) |
589 600 | |
590 601 | # Set channel selection to take the average of all channels. |
591 602 | openms.selectchannel(1, 0, nchans, 1) |
592 603 | |
593 604 | # Extract data from MS. |
594 - | data = openms.getdata(['corrected_data', 'flag']) |
605 + | data = openms.getdata(['corrected_data', 'flag', 'antenna1', 'antenna2']) |
595 606 | except: |
596 607 | # Log a warning and return without data. |
597 608 | LOG.warning('Unable to compute mean vis for intent(s) {}, field {}, spw {}' |
598 609 | ''.format(intent, fieldid, spwid)) |
599 610 | return |
600 611 | |
601 612 | # Combine the two dictionaries. |
602 613 | data = {**data, **weights} |
603 614 | |
604 615 | return data |