Inside-out map-making

If the time-series data supplied to SMURF:MAKEMAP does not consist of a single continuous sequence of time samples, it will divide the data up into two or more "chunks" of continuous samples. A separate map is made from each chunk in turn, using the iterative procedure  described in SUN/258, and these maps are simply co-added at the end to create the final returned map. This procedure can be summarised as follows:

initialise final map to zero

for each chunk:

   initialise chunk map to zero

   for iteration = 1 to numiter

      update chunk map using time-series data in current chunk

   next iteration

   finalise chunk map

   add the chunk map into the total map

next chunk

In general, it seems that the more data that is included in each iteration, the more constrained is the map, and so the fewer the spurious features. So it would be nice if we could turn these loops "inside out" so that the iteration loop occurred outside the chunk loop, causing data from all chunks to be included in the map at the end of each iteration. 

initialise final map to zero

for iteration = 1 to numiter:

   for each chunk:

      update final map using time-series data in current chunk

   next chunk

next iteration

finalise final map

As an aid to experiment, I have added some simple features to MAKEMAP to allow something like this effect to be achieved, at the expense of considerably longer run-times and more effort on the part of the user. If experiments with these facilities show that there is a significant benefit in turning the loops "inside-out", then it may be possible to implement the idea in a more time-efficient manner.

The importsky configuration parameter:

The main feature in the new scheme is the facility to supply an initial estimate of the sky when running MAKEMAP. This is done by specifying 

importsky=ref

in the config file, and then assigning the NDF containing the  initial sky values to the "REF" ADAM parameter on the MAKEMAP command line:

% makemap ref=mysky.sdf

Before the first iteration starts, the specified sky map is sampled at the position of each bolometer value, and each sampled sky value is subtracted from the corresponding cleaned bolometer value. The remaining residuals are then processed as normal by estimating and subtracting (typically) COM and FLT models. The values sampled from the sky map are then added back onto the residuals, which are binned into a new map - the first "itermap". From then on, the algorithm proceeds to do the second and subsequent iterations in the usual manner.

This allows us to produce the "inside-out" looping described above by repeatedly re-running MAKEMAP, doing only one iteration on each run, and  supplying the output from one invocation as the initial sky map for the next invocation. In effect, you are doing the "iteration" loop by hand at the shell level.

Speeding things up:

This is slow. But with patience some experimentation, inside-out looping can be performed. It can be made slightly faster by caching some of the processing performed by the first invocation for use by subsequent invocations. Two things can be done:

1. Caching the cleaned bolometer values. If "exportclean=1" is included in the config file on the first invocation of MAKEMAP, the cleaned time-series data will be written to disk in NDFs that end with the suffix "_cln.sdf". These NDFs can be used in place of the raw data as input for all subsequent invocations of MAKEMAP. In this case you should create a new config file for the second and subsequent invocations and add "doclean=0" to it, so that the cleaned input data will not be re-cleaned.
2. Caching the EXT model values. Calculation of the EXT model values happens only once, before the first iteration starts. So if you add "exportndf=ext"  to the config file for the first invocation, a set of NDFs with suffix "_ext.sdf" will be created holding the EXT model values. These can be supplied as input to the second and subsequent invocations by adding "ext.import=1" to the config file. Note, MAKEMAP uses fixed pre-defined names for the NDFs when writing and reading EXT model values, so the NDFs created on the first invocation should not be moved or re-named. We also add "noexportsetbad=1" to prevent EXT values for flagged bolometers being set bad in the exported NDFs (since a different set of bolometers may potentially be flagged as bad on subsequent invocations). 

There are other possibilities for speeding things up, such as caching the LUT model, but these require changes to the code, which have not yet been implemented.

Other parameter settings:

• Setting "numiter=1" is required to ensure only one iteration is performed by each invocation of MAKEMAP.
• The NOI model, which contains estimates of the noise in each bolometer time stream, is normally calculated at the end of the first iteration, once the final residuals are known. This is no good for us here since each invocation only performs one iteration. The simplest solution is to add "noi.calcfirst=1" to the configuration for every invocation. This forces the noise estimates to be made before the start of the first iteration.
• Care needs to be taken if any AST masking is used. Firstly, since we have set "numiter=1", the first iteration is also the last iteration. So we need to add "ast.zero_notlast=0" to the configuration for every invocation. Without this, no masking would be performed on any of the invocations of MAKEMAP. However, we do want to suppress masking on the final invocation, and so we should revert to the default value of "ast.zero_notlast=1" for the final invocation.
• If an external mask is to be used, it should be supplied as normal as the REF parameter on the first invocation. The output map from the first invocation will have this same mask and so will mask the AST model correctly when supplied as the REF parameter on the second and subsequent invocations. 

An example:

In this example, four invocations are performed. How to choose the best number of invocations is still to be determined. Continuing until the RMS change within the source regions between subsequent maps  is less than some target figure may be the way to go. 

First invocation:

% makemap in=/SCUBA2/s8\*/20111114/00055/\* ref=mymask config=^conf1 out=map1

We make an initial map from the raw 850 um data files for 20111114 obs. 55 using "mymask.sdf" as an external mask. The output map is put into "map1.sdf". The configuration file "conf1" contains:

^/star/share/smurf/dimmconfig_bright_extended.lis

numiter=1          # Only do one iteration

850.flt.filt_edge_largescale=600

450.flt.filt_edge_largescale=600

ast.zero_mask=1

ast.zero_notlast = 0          # Mask the final (i.e. first) iteration

noi.calcfirst=1          # Calculate NOI before iterating

exportNDF=ext          # Export the EXT model values to NDFs
noexportsetbad=1        # Export good EXT values for bad bolometers

exportclean=1          # Export the cleaned  time-series data to NDFs

Second invocation:

% makemap in=./\*_cln.sdf ref=map1 config=^conf2 out=map2

We make a new map from the cleaned data files "*_cln.sdf" created by the first invocation (we assume the current directory contains no other cleaned data, for instance from an earlier run). The map created by the first invocation is supplied as the REF image, and the new map is put in "map2.sdf". The "conf2" file sets the importsky parameter to indicate that the initial sky supplied by the REF image should be subtracted from the time-series before starting the first (and only) iteration. 

^conf1          # Inherit all the settings in "conf1"

exportNDF=0         #  Do not export the EXT model this time

exportclean=0          # Do not export the cleaned data this time

doclean=0          # No need to clean the data this time

importsky=ref          # Subtract off the REF values at the start

ext.import=1          # Import the EXT values created on the first invocation

Third invocation:

% makemap in=./\*_cln.sdf ref=map2 config=^conf2 out=map3

Later invocations use the map created by the previous invocation as REF, and use the same config file as the second invocation.

Last (=fourth) invocation:

% makemap in=./\*_cln.sdf ref=map3 config=^conf3 out=final-map

On the last invocation we need to prevent the mask being applied, so we use a new config:

^conf2          # Inherit all the settings in "conf2"

ast.zero_notlast = 1          # Do not mask the final (i.e. first) iteration