Preventing self-masks from diverging

Masks can be used in various ways in MAKEMAP. For instance, the dimmconfig_bright_extended.lis file specifies:

ast.zero_snr = 5

which causes the areas of the map that have a signal-to-noise ratio of less than 5 to be masked on each iteration. This is a dynamic mask, in that a new and potentially different mask is generated and used on each iteration. However, the added flexibility introduced by allowing such a mask to change on each iteration can result in eventual divergence of the iterative  map-making algorithm. To prevent this, a set of new configuration parameters has been added that result in the mask being frozen after a specified number of iterations. These parameters are AST.ZERO_FREEZE, COM.ZERO_FREEZE and FLT.ZERO_FREEZE. Each of these parameters may be set to a positive integer value, indicating that he corresponding mask should be frozen after that many iterations. They all default to zero, which reproduces previous behaviour (i.e. the mask is never frozen).

Better core utilisation

Over the past couple of weeks I've modified several parts of the makemap command to split up their calculations between all available threads rather than doing them all in a single thread. This means that makemap should now run significantly faster on multi-core machines. The speed up will be greater for larger data sets than for smaller data sets. As an example, producing a map using a single subarray from a  40 minute observation took  11 mins 49 secs with the new scheme on my 12 core machine, and achieved an average CPU load of 791%. With the old scheme, the same map took 19 minutes 14 secs with an average CPU load of 461 %.

Changes to the common-mode estimation

For the impatient: the new improved algorithm described below is enabled by default and requires no changes to scripts or configuration files. However, the old algorithm is still available, and can be used by setting COM.OLDALG=1 in the configuration file.

For extended sources, MAKEMAP usually uses the mean normalised change in map pixel values between iterations as a measure of convergence. This value drops rapidly over the first few iterations, and eventually levels out, typically somewhere between 0.01 and 0.05. However, if MAKEMAP is allowed to continue iterating beyond this point, the mean normalised change between iterations often starts to increase and decrease by large amounts in an erratic manner. Examining the maps shows that these changes reflect the creation of nasty artifacts in the map, which often take the form of individual bolometer tracks. As an example, Dave Nutter saw the following variations in mean normalised change with iterations for one particular observation:

A relatively smooth decrease until about iteration 50 and then madness. For illustration the following shows the map after 50 iterations (left) and 200 iterations (right). Notice the visible streaks that have appeared by iteration 200.

The cause of this wild behaviour seems to be the way the common-mode signal is estimated and used. On each iteration, the common mode value at each time slice is taken to be the mean of the bolometer values at that time slice, averaged over all remaining bolometers  (i.e. bolometers that have not previously been flagged as unusable for one reason or another). Each bolometer time series is then split into blocks of 30 seconds and compared to the corresponding block of the common-mode signal. This comparison takes the form of a least squares linear fit between the bolometer value and the common-mode value, and generates a gain, offset and correlation coefficient for each bolometer block. Bolometer blocks that look "unusual" - either because they have a low correlation coefficient (i.e. do not look like the common mode) or an unusually low or high gain - are flagged, and all the bolometer values within such blocks are rejected from all further processing. So on the next iteration, the common mode within each block is estimated from only the subset of bolometers that passed this check on all previous iterations. With the passing of iterations, more and more bolometer blocks are rejected, often resulting in the common mode in a block being estimated from a completely different subset of bolometers to those of the neighbouring blocks. This produces discontinuities in the common-mode signal at the 30 second block boundaries, which get larger as more iterations are performed. When the common-mode is subtracted from the data, these discontinuities are transferred into the residuals, which are then filtered within the FLT model. The sharp edges at the discontinuities cause the FFT filtering to introduce ringing (strong oscillations on the scale length of the filter that extend over long times). The following plot shows an example of these discontinuities in the common-mode estimate:

The horizontal axis is sample index. The red curve is the common mode at iteration 85 and the white curve is the common mode at iteration 99. Due to down-sampling, 30 seconds corresponds to 2400 samples, and discontinuities at intervals of 2400 samples are clearly visible. In addition, extensive oscillations are visible, particularly between samples 48000 and 50000.

In order to avoid this problem, changes to the way the common-mode is estimated and used have been made. The heart of the change is that bolometer blocks that are flagged as unusual are no longer excluded from the estimation of the common-mode on the next iteration - all such flags are cleared before the new COM estimate is made. So the common mode is always based on data from all bolometers, thus avoiding the discontinuities at the block boundaries. Other more minor changes in the code include:

1. The common-mode value at each time slice can be estimated using a sigma-clipped mean, rather than a simple mean. The new configuration parameters COM.NITER (the number of iterations) and COM.NSIGMA (the number of standard deviations at which to clip on each iteration) control this algorithm. In practice, sigma-clipping does not seem to add much, and may slow down convergence. At the moment, both of these parameters default to 3. Changing COM.NITER to 1 results in a single simple mean being calculated with no clipping. This may become the default in the near future.
2. No attempts are now made to refine the common mode. Previously, once the original COM signal was estimated and unusual bolo-blocks were flagged, a refined COM signal was formed by excluding the flagged bolo-blocks. Further unusual bolo-blocks were then flagged by comparing them with the refined COM signal. This refinement process was repeated until no further bolo-blocks were rejected. In the new code, no attempts are made to refine the original COM signal resulting in the new algorithm being significantly faster than the old algorithm.

Note, the new algorithm still uses COM.GAIN_BOX to define the size of the blocks in which the gain and offset of each bolometer is estimated. Also, "unusual" blocks are still flagged in the same way and excluded from the subsequent FLT model and map .

The new algorithm seems effective at overcoming the problem of divergence illustrated at the start of this post. For instance, with the new code the variation of mean normalised change between iterations for the same data looks as follows:

 Note that:

• The wild variations above iteration 80 have gone.
• The curve drops somewhat faster at the start (i.e. fewer iterations are needed to achieve a given maptol).
• Each estimation of the common-mode signal is faster than before due to the lack of the refining process.
• Far fewer samples are rejected from the final map ( 2.38 % of the data is now flagged by the COM model, as opposed to 17.37 % previously).
• The resulting maps at iteration 50 and 200 are shown below.  There is very little difference between them, and they both look visually like the iteration 50 map from the old algorithm shown above.

Interrupting makemap using control-C

So you're sitting there waiting for a long run of makemap to finish thinking "It's going to take ages! I asked for far more iterations than I need, but I can't restart it because it's taken so long to get to this point". As of git commit efde7077500 you can simply hit control-C and all will be well - makemap will continue running until the end of the next iteration. It will then create a map from the information it has, and terminate normally. Note, it needs to perform one more whole iteration before creating the map in order to get proper pixel values for the areas masked by ast.zeromask. One thing to beware of though - if the input data is split into more than one chunk, subsequent chunks that have not yet been started will not be included in the map.

If you want to abort NOW! without a map,  rather than waiting for the next iteration to complete, then press control-C a second time.

Note, if you are running makemap within a shell script, then the shell may handle the control-C signal itself, leading to some potentially odd behaviour. For instance,  the script may appear to terminate immediately, but in fact may leave the makemap process running in the background until the next iteration has completed. Most shells have ways of controlling what happens when an interrupt signal is detected. For instance, the (t)csh has the "onintr" command.

Masking of FLT and COM models

In previous blogs (e.g. Suppressing large-scale noise with zero-masks), Ed has described how masking the AST model - the astronomical signal - between iterations can reduce the amplitude of low frequency structures in the final map. Over the last month or two, I have extended the uses of masks to allow the COM and FLT models to be masked in a similar way to the AST model.

COM masking: The common-mode signal is the average of all available bolometer time streams. Once found, a least-squares linear fit is performed between the average and each individual bolometer time stream. This fit generates a scale factor and offset, together with a correlation coefficient. Bolometers that have a low correlation with the common-mode are rejected as bad. At any one time, the astronomical source is seen only by a subset of the bolometers, and so will be largely absent in the common-mode estimate. This means that bolometers that see the source will have a lower correlation to the common mode than those that do not see the source, which can result in bolometers being rejected unnecessarily. To avoid this, it is now possible to exclude samples that see the source from each time stream when finding the common-mode signal, and also when performing the least squares linear fit. This can produce a better estimate of the common-mode, unpolluted by bright sources samples, and also better estimates of the correlation coefficients, resulting in fewer bolometers being rejected.

The masking of the COM model is controlled by the new configuration parameters COM.ZERO_MASK, COM.ZERO_CIRCLE, COM.ZERO_LOWHITS and COM.ZERO_SNR. These are all analogous to the corresponding set of AST.xxx parameters used for masking the AST model. As an example, a configuration of:

^/star/share/smurf/dimmconfig_bright_compact.lis

com.zero_circle = 0.02

was used to produce an 850 um map of Mars from observation 40 on 20120202. This excludes samples from the COM model that are within a circle of radius 0.02 degrees, centred on the source. The resulting final map contained 6744407 samples (or 2877.31 bolos). Without the COM masking, the final map contained 4676093 samples (1994.92 bolos).

Sadly, not all observations show such an improvement using COM masking. Indeed, some show little or no improvement, and some actually show a decrease in the number of samples in the final map. For this reason, COM masking is not switched on by default in any of the dimmconfig files distributed with SMURF.

FLT masking: A similar set of configuration parameters - FLT.ZERO_MASK, FLT.ZERO_CIRCLE, FLT.ZERO_LOWHITS, etc., are now available to control masking of the FLT model. The FLT model holds the low frequency background for each bolometer time stream, and is estimated by smoothing the time stream using a filter defined by FLT.FILT_EDGE_LARGESCALE, etc. On the first iteration, any bright astronomical source in the data will result in the smoothed background being over-estimated in the vicinity of the source. This in turn, produces a negative bowl around the source when the background is removed. It can take many iterations for the map-maker to reduce this bowling to acceptable levels. Using the new FLT masking parameters allow the source to be omitted before estimating the background in each bolometer time stream (the source samples are replaced by a values interpolated from the neighbouring data, with added noise). Masking the FLT model in this way on the first two iterations can have a dramatic effect on the number of iterations need to reach convergence. For compact sources, it can typically reduce the number of iterations required to 50% of the original number, and sometimes fewer.

For this reason, FLT masking is switched on by default in the dimmconfig_bright_compact.lis file, using the settings:

flt.zero_circle = 0.0167

flt.zero_niter = 2

The effect of FLT masking on extended sources is marginal, and so is not switched on by default in dimmconfig_bright_extended.

SCUBA-2 Calibration: REDUX.

The short story:

• The heater coupling factors have been adjusted to more realistic values. In practice, this does not change the performance of the instrument - however it does change the absolute value of the FCFs. These values were adjusted in the software in mid-December.
  

• The WVM tau algorithm has been fixed and improved. This will not affect you directly: though the nightly plots now look extremely good and are officially used for weather band determination.

• This has allowed new, and better calculation of the relation between the 225GHz tau derived from the WVM and the opacities at the two SCUBA-2 filter-bands. They are now as follows:

TAU_[850] = 4.6 * (TAU_[225] - 0.0043)
TAU_[450] = 26.0 * (TAU_[225] - 0.019)

• The FCFs (flux conversion factors) have been derived for both wavelengths from an extensive reduction of calibrator sources observed over eight months of SCUBA-2 commissioning and science verification observations. They are as follows:

850um:
FCF_[arcsec] = 2.42 +/- 0.15 Jy/pW/arcsec**2
FCF_[peak] = 556 +/- 45 Jy/pW/beam
Beam area = 229 arcsec**2

450um:
FCF_[arcsec] = 6.06 +/- 0.32 Jy/pW/arcsec**2
FCF_[peak] = 606 +/- 55 Jy/pW/beam
Beam area = 97 arcsec**2

• Reminder on how to calibrate your data:

Other posts discuss how best to reduce your data (and what recipes are needed). The latest software releases (since January 2012) all include extinction correction (with the relations above) and the changed coupling factors. If you reduced your data prior to this, you will need to reduce them again to account for these changes. Applying the FCFs reported here to old reductions of your data will be wrong.

• The arcsec FCF: (when you want integrated fluxes)

The arcsec FCF is the factor by which you should multiply your map if you wish to use the calibrated map to do aperture photometry.

• The peak FCF: (the FCF-formerly-known-as-beam):

This FCF is the number by which to multiply your map when you wish to measure absolute peak fluxes of discrete sources.

The whys and the wherefores (the gory details):

• Reductions of the calibration observations:

Uranus and Mars were used as the primary calibrators for these results. In addition, CRL 618, CRL 2688 were also predominant secondary calibrators. All of the calibrators were reduced in January 2012 using the updated dimmconfig_bright_compact.lis. The improvements in the recipe (mostly in how it chooses to stop iterating) have resulted in extremely flat maps with nearly no evidence of the 'bowling' seen around strong sources during S2SRO reductions. To improve the accuracy of the peak-fitting and aperture photometry, the maps were reduced with 1 arcsecond pixels at both wavelengths.

Once the maps were reduced they were then analysed using the PICARD script SCUBA2_FCFNEFD. There have been some changes to this script: some few bugs have been fixed, the average FCF's have been adjusted, as well as the (now emprically derived) beam area, and a few reference fluxes have been adjusted. FCF_beamequiv has been removed entirely, and all calculations of integrated values are now done using AUTOPHOTOM, with a defined aperture and annulus for background subtraction.

Following extensive analysis of the most optimal parameters, all calibrations were reduced using a 60" diameter aperture (at both wavelengths) with an annulus between 90" and 120" from the source position.

• Questions? Let's provide answers to a few we've already seen:
  

• "These FCFs are very different to the old numbers quoted!"
  

Yes they are. The heater coupling factor change and the new tau relations play a significant part in this. But in addition, the large sample of observations has allowed for a much more accurate determination of the beam area. At 450um in particular, optical effects of the telescope show that the error beam is large and the beam is not gaussian. This results in an effective FWHM that is much broader than the 7.5" quoted previously (though that is the approximate FWHM of the fit to the centre of the beam) - it is more like 9.5", taking into account the error beam. Therefore, the measured (and fitted) peak is relatively lower, requiring a higher FCF to calibrate the peak flux in your data.

• "There is a lot of scatter when I calculate the 'beam' or peak FCFs for my calibrators (particularly at 450um)"

No kidding. Peak values are obtained either from reading off the peak of the map (in gaia or by another method) or by fitting to the peak using beamfit (as is done in PICARD). The beam shape (particularly at 450um) can be extremely susceptible to changes in focus and atmospheric instability, amongst other things. The integrated value (FCF_arcsec) is more robust against such changes. If you are measuring a peak fit from a calibrator and see a strong deviation from the expected value things to check are:

- how 'focussed' does the image look? If you see distortion in the shape of a source that should be point-like, or distortion or 'shoulders' in the beam then it is likely that the peak value will be unreliable.

- was the observation taken early in the evening? Focus and atmospheric effects are known to be worst in the early evening hours and sometimes in the morning after sunrise. If you are looking at calibrators, try and look at ones taken later in the night and see if there is improvement.

A 'trap' has been set in PICARD to warn you if the attempted fit to the peak misses the actual peak value by more than 10%. Looking at the fit to the shape also helps in this instance. In any case, the quoted peak FCF value at the top of the post is derived from the arcsec FCF and the empirical beam area derived from nearly 500 observations at both wavelengths and this number has been shown to be robust.

• "How stable are these FCFs? (read: do I need to reduce my own calibrators?)"
  

Very. The absolute errors at both wavelengths are within 5% and no significant trends have been seen in the last six months. Instrument performance is being monitored very closely and any deviations are likely to be noted specifically. However, we do not discourage you taking calibrators from the nights your data was taken and reducing them yourself - we appreciate the sanity checks! Another handy rule: if you do it to your data, do it to your calibrator. If you have specific methods you plan to use on your data, apply the same methods to your calibrator in order to ensure your calibration is correct. We are now happy to say though, that these FCFs look stable and correct, so using these numbers should provide you with well-calibrated data.

• "What happened to FCF_beamequiv?"

FCF_beamequiv is a seductive, evil little value that tempted us to stray to the dark side, albeit temporarily. In essence it was created to use as a comparison to SCUBA performance, but should never have been used to actively calibrate SCUBA-2 data as it assumed a perfect gaussian beam. The statements above explain that this is patently untrue, especially at 450um. The beamequiv number was quoted previously, and incorrectly, as the true FCF, and it is largely the reason that the new (and correct) numbers seem so much larger. We have banished it from PICARD and it shall now be known as the FCF-that-shall-not-be-named.

Supplying an externally-generated zero mask

In an earlier post I described zero-masking , in which the map is forced to a value of zero in certain locations to help with convergence (i.e. reduce saddles, negative bowls around sources etc.). The two methods that have been available until now are: defining circular regions beyond which the map is set to zero (good for known compact sources); and using a threshold on the local SNR (an automatic way to identify source and blank sky).

Since SMURF still handles non-continuous pieces of data independently (i.e. multiple maps of the same field), the SNR of an individual map may not be sufficiently large to produce a good mask of source / blank sky, whereas the combination of all data sets do. For those cases, and other examples where you may know, a priori, where to expect emission, a new facility has been added so that the user can define a zero mask externally and then provide it to SMURF.

In this illustrative example, I use one of the OMC-1 maps from 2010 with the s4a array: obs 18,22,23 on 2010216 and obs 27,29,33,34 on 20100218. First, we make a map using the bright_extended configuration (which uses SNR-based zero masking; I also use a lot of down-sampling and big map pixels to speed things up for this example):

makemap ^names_450.txt map_450_be pixsize=8 \
method=iterate \
config='"^/stardev/share/smurf/dimmconfig_bright_extended.lis"'

The resulting image looks reasonably good, but there are obvious negative bowls around the bright extended sources. However, the SNR of the final image is quite good, so I create a new mask based by thresholding the SNR map (4-sigma), and then smoothing it with a 5-arcsec FWHM Gaussian to spread out from where the sources are slightly:

makesnr map_450_be map_450_be_snr

thresh map_450_be_snr mask_450 thrlo=4 newlo=0 \
thrhi=4 newhi=1

gausmooth mask_450 mask_450_sm fwhm=5

thresh mask_450_sm zero_mask_450 thrlo=0.05 newlo=bad \
thrhi=0.05 newhi=1

I can then feed this new "zero_mask" back into makemap and re-run the data:

makemap ^names_450.txt map_450_zm pixsize=8 \
method=iterate \
config='"^/stardev/share/smurf/dimmconfig_bright_extended.lis,itermap=1,ast.zero_mask=1,ast.zero_snr=0"' \
ref=zero_mask_450


The key things to note are: (i) I turned off the SNR thresholding, ast.zero_snr, that is the default for the bright_extended configuration, and (ii) the new zero mask is provided as the REF image (also ensuring that the map and mask are on the same pixel grid), and to use it we set ast.zero_mask=1.

In the following images I show (from top-left): (i) an arm sticking out towards the north of OMC-1 using the default bright_extended reduction; (ii) the same region after using the updated zero-mask; (iii) the first map subtracted from the second map to illustrate the improved response to large-scale structure; (iv) the original combined zero mask based on the SNR of map pixels in individual scans; (v) the new zero mask based on thresholding/smoothing the combined SNR map: