How to loose fewer samples with NOI.BOX_SIZE

At the end of the first iteration, when the common-mode signal (COM model), low frequency drift (FLT model), and the astronomical signal (AST model) have all been removed from the data, the remaining residuals are used to estimate the variance of the high frequency noise in each bolometer (the NOI model). The reciprocal of these variances are used as weights on the second and subsequent iterations, and are applied to the time-series values that are binned to form the map.

By default, the NOI model contains a single variance value for each bolometer , and so all samples from a single bolometer get the same weight in the map. However, the NOI.BOX_SIZE parameter can be used to force each bolometer to have multiple variance estimates, with one value for each box of samples. See this post. 

Whilst this is often in an improvement over the default situation in that it allows for varying noise levels, it has two potential issues:

1. All samples in a box are given the same variance value, thus introducing the potential for steps in the weights and corresponding artifacts in the map.
2. The noise in each box is found by taking the FFT of the data and finding the average power in the 2 to 10 Hz band. Using an FFT is problematic on short time streams, as any missing values (such as caused by steps, jumps, common-mode flagging, etc) have a proportionally larger effect, causing the FFT to be less representative of the data. A requirement that the time stream should have a majority of unflagged values in order to be usable means that a lot of data is often rejected because a reliable FFT cannot be obtained.

An alternative scheme that avoids these two issues has been added to makemap, which is enabled by setting "NOI.BOX_TYPE =1" in your configuration. This option causes the noise at a sample to be set equal to the variance of the neighbouring sample values - those that are less than half a box size (i.e. NOI_BOX_SIZE / 2) away from the sample. This scheme does not involve taking an FFT and produces a continuously varying estimate of the variance and so avoids steps in the weights used to make the map. It also is minimally affected by flagged samples and so results in fewer data samples being flagged as unusable (as shown by the "NOISE" entry in the samples statistics at the end of each iteration).

Weighting different observations in makemake and skyloop

Both makemap and skyloop can be used to produce maps from multiple overlapping observations. Up to now, it has not been possible to control the weights used for each observation when forming the total co-added map. Instead, fixed weights equal to the reciprocal of each pixel variance were used. The new CHUNKWEIGHT configuration parameter can now be used to specify a weight for each observation. The weights determined from the pixel variances are multiplied by these user-specified weights, which default to unity.

The weights can be specified in two ways:

1. The CHUNKWEIGHT configuration parameter can be set to an arbitrary algebraic expression in which the variable names are the names of FITS keywords. The keywords must exist in the input data and have numerical values. The expression should be enclosed in double quotes. As an example, you could include this line in your configuration file:
   
      chunkweight = "AMSTART * WVMTAUST"

The values of the FITS keywords are read from the observation header and substituted into the expression to get the weight to use for the observation.


2. A comma-separated list of explicit numerical weights can be supplied, enclosed in parentheses. If insufficient values are given in the list, the last value is replicated as often as is necessary.
   
      chunkweight = (1.0,2.0,3.0)

Holes in heterodyne maps

HARP has only twelve operational receptors since May 20.  If one or two also suffer from interference and fail quality assurance, then there can be locations in the reduced spectral cubes where no valid data are present.  Such missing spectra (shown in blue) will often be arranged in a grid pattern as in the example below.

This can create issues if you wish to integrate flux or simply improve the cosmetic appearance.

One approach is to smooth the data but this degrades all the spectra.  Instead use KAPPA::FILLBAD. This replaces just the bad values with a smooth function derived from neighbouring values, derived by iterating to a solution of Laplace's equation.

Just using FILLBAD's default values will duly fill in the holes.  You might prefer not to smooth in the spectral axis by setting Parameter SIZE to [s,s,0] where s is the initial scale length in pixels.  The zero means do not smooth along the third axis. The scale length should have a value about half the size of the largest region of bad data to be replaced.  Since the largest bad regions apart from the cube peripheries are two pixels across, a command like this

fillbad in=holey out=filled  size="[1,1,0]"

is appropriate.  The graphic below shows the same region as before, but with the holes filled in.

FILLBAD does what it says on the tin.

There will be some smoothing, but it's not obvious.  Below is a KAPPA::CLINPLOT graphic for a part of the filled cube.  Each spectrum is a pixel in the image above. Can you identify the three spectra which were previously bad?

Inclusion of CSO Tau Fits for Determining FCFs

During the first few months of 2013 the WVM at JCMT had several periods of instability where it was unreliable for determining the value of tau. Because of this, members of the science team collaborated to produce a collection of smooth polynomial fits to the tau values from the 225-GHz tau meter at the CSO for the affected nights, which can be used to perform extinction correction in place of the WVM tau values.

In the latest version of starlink (which you can get by rsyncing from /stardev at the JAC) makemap will now first check to see if the date of the observation is one of the (somewhat sporadic) dates when the WVM was unstable. This occurs as long as the ext.tausrc parameter is set to "auto," which it is by default. If the date is one of the affected dates, makemap will look for an available fit to the CSO data in the file specified by the parameter ext.csofit, which is set up by default to refer to the collection of CSO fits produced at the JAC. If makemap cannot find a fit for an observation in the specified path it will print a warning that no fit or WVM data is available and refuse to reduce the observation, though this shouldn't ever happen in normal operation.

If you have observations taken between January 19 and May 14 of this year, using the latest version of starlink rsynced from /stardev at the JAC will help ensure that you get the best extinction data available.

The graph below shows a comparison between the FCFs derived from WVM data and those derived using the new CSO tau fits over the period of January-May 2013. The blue diamonds represent FCFs from the WVM, and the tan circles are FCFs from the CSO tau fits.

Removing linear striations in maps

You've probably seen several maps that show the faint linear striations visible in the following 850 map.

Well, in at least some cases, these striations seem to be due to bad bolometers that have a very different common-mode signal to the others. The map-maker has an algorithm that looks for and rejects such bolometers, but by default the rejection criteria are fairly  weak, so few samples are rejected. Changing the  com.corr_abstol value from its default of 0.2 to 0.8 produces the following map:

The striations are much less noticeable (both these images use the same grey scale). The down-side to this is that since more bolometer samples are rejected, there are fewer samples to form the map. For the first map, only 1.6 % of the available samples were rejected due to common-mode mis-match, leaving a total of 65.4 % available for the map. For the second map, 11.1 % were rejected due to common-mode mis-match, leaving 55.9 % available for the map. 

Going even further, the following map was created with com.corr_abstol set to 0.97 and  com.gain_fgood set to 0.5:

Here, 38.5% of samples were rejected due to common-mode mismatch, leaving only 28.6 % available for the map. You can see that the striations have been almost completely removed, but at the expense of higher noise.

Automatic FLT masking

The estimation and subtraction of low frequency noise in each bolometer (the FLT model) very often causes ringing around bright sources. This ringing can be very bad for the first few iterations but often improves if sufficient iterations are performed. Masking the FLT model consists of excluding known bright sources from the estimation of the low frequency noise, and thus reducing the ringing. The ability to do this using the "flt.xxx" set of configuration parameters has been around for while, but sadly the ability to do automatic   FLT masking has not been available - an external mask image needs to be supplied. By comparison, the AST model can be masked automatically using the ast.zero_snr and ast.zero_snrlo parameters. The corresponding parameters flt.zero_snr and flt.zero_snrlo do exist, but do not really work since the SNR values are only available once a map has been made. By then it is to late to do the FLT masking - the FLT model has already been applied and the ringing is therefore already present in the residuals.

But a new parameter called ast.skip can get round this. This new parameter defaults to zero, but if set to a positive integer it gives the number of initial iterations for which the AST model (i.e. astronomical signal) is to be skipped. For example, if ast.skip is set to 5 then the first five iterations will not include an AST model. The map will still be created at the end of each iteration, but it willnot be used to find the expected astronomical signal, and the residuals will be left unchanged. This means that all these first five iterations will start from essentially the same time series data - the initial cleaned raw data. However, what this means is that we can now create an FLT mask automatically from the SNR values in the map, and refine that mask five times. So when we get to the sixth iteration, we have a usable FLT mask and we have not  introduced any ringing into the residuals.

It's a bit like doing a separate run of five iterations to determine the FLT mask before starting again to do the main run.

In summary, to use automatic FLT masking, add the following to your config:

ast.skip=5
flt.zero_snr=5
flt.zero_snrlo=3

If you use a fixed number of iterations, you may want to increase your numiter value by the ast.skip value, so that you are doing the same number of real iterations. As usual, you can play around with the precise values of these parameters.

FLT masking speeds up convergence as shown in the following plot of normalised change in the map against iteration number - red uses the new scheme described above, green uses a fixed external FLT mask, and blue is without FLT masking.

Changes to way makemap handles control-C interupts

There have been some changes to the way that makemap behaves when control-C is pressed (or a SIGINT signal is received by the process by any other means). When the first control-C is detected, makemap will continue to complete the current iteration (but a second interupt will cause makemap to terminate immediately). It will then pause and ask the user what to do next. The options are:

1. abort immediately with an error status
2. close the application returning the current output map
3. do one more iteration to finialise the map and then close

Option 3 is important if you are using AST masking, since one final iteration should always be done without masking to ensure correct fluxes in the background regions (see the AST.ZERO_NOTLAST config parameter). If you are not using AST masking, or if you will be re-starting makemap to perform further iterations, then you should select option 2.

At some point I will add a fourth option to ignore the interupt and continue iterating to convergence.

Note if you are running makemap from within a script, then the handling of control-C interupts will probably be quite different, because the shell (or perl, python or whatever) will catch the interupt before it gets to makemap. You may get some sort of useful behaviour by finding the process ID for the makemap process itself, and sending a SIGINT to that process explicitly:

% ps aux | grep makemap

dsb      25407  0.0  0.0  43716  1952 pts/5    S    Jul03   0:00 makemap

% kill -s INT 25407