Skip to content

Commit

Permalink
smurf: Do not update quality on each iteration of skyloop
Browse files Browse the repository at this point in the history 
Because we may end up with more than 8 bits, in which case they get
mangled together so that they cannot later be reconstructed.
  • Loading branch information
David Berry committed Mar 5, 2014
1 parent 3cd8b24 commit 90d7ba5
Show file tree
Hide file tree
Showing 3 changed files with 17 additions and 168 deletions.
15 changes: 8 additions & 7 deletions applications/smurf/libsmf/smf_calcmodel_ast.c
View file
Expand Up @@ -137,12 +137,14 @@
* first iteration (we will have noise values on the first iteration
* when running from SKYLOOP).
* 2014-1-29 (DSB):
* Use the map quality array rather than the raw mask array to define
* the areas to mask. The quality array contains the raw mask but also
* masks out all map pixels with bad data values or variances. It is
* the quality array, not the raw mask, that is used in smf_iteratemap
* Use the map quality array rather than the raw mask array to define
* the areas to mask. The quality array contains the raw mask but also
* masks out all map pixels with bad data values or variances. It is
* the quality array, not the raw mask, that is used in smf_iteratemap
* when adding on the previous AST model prior to forming a new map,
* so we really must be consistent and use the same thing here.
* 2014-3-4 (DSB):
* Do despiking even if an initial sky is being subtracted.
* {enter_further_changes_here}
* Copyright:
Expand Down Expand Up @@ -294,9 +296,8 @@ void smf_calcmodel_ast( ThrWorkForce *wf __attribute__((unused)),
have_noi = ( ((double *) noi->sdata[idx]->pntr[0])[ 0 ] != 1.0 );
}

/* Despike if we have usable NOI values and this function was not
called as part of subtracting off an initial sky. */
if( (mapspike > 0) && have_noi && !(flags&SMF__DIMM_PREITER) ) {
/* Despike if we have usable NOI values. */
if( (mapspike > 0) && have_noi ) {
size_t nflagged;
smf_map_spikes( wf, res->sdata[idx], noi->sdata[idx], lut->sdata[idx]->pntr[0],
SMF__Q_GOOD, map, mapweight, hitsmap, mapvar, mapspike,
Expand Down
140 changes: 3 additions & 137 deletions applications/smurf/libsmf/smf_initial_sky.c
View file
Expand Up @@ -63,6 +63,9 @@
* - Clear and RING or COM flags in the quality array of the time series data.
* 23-JAN-2014 (DSB):
* Mask the AST map in smf_calcmodel_ast as normal, rather than masking it here.
* 4-MAR-2014 (DSB):
* Do not clear and RING or COM flags in the quality array of the time
* series data (it will never be present anyway).
* {enter_further_changes_here}
* Copyright:
Expand Down Expand Up @@ -100,49 +103,22 @@
#include "libsmf/smf.h"
#include "libsmf/smf_typ.h"

/* Local data types */
typedef struct smfInitialSkyData {
dim_t b1;
dim_t b2;
dim_t ntslice;
size_t bstride;
size_t tstride;
smf_qual_t *qua_data;
} SmfInitialSkyData;

/* Prototypes for local static functions. */
static void smf1_initial_sky( void *job_data_ptr, int *status );

void xxdump( const char *text, smfData *data, int *status );

int smf_initial_sky( ThrWorkForce *wf, AstKeyMap *keymap, smfDIMMData *dat,
int *iters, int *status ) {

/* Local Variables: */
SmfInitialSkyData *job_data = NULL; /* Data describing worker jobs */
SmfInitialSkyData *pdata = NULL; /* Data describing worker jobs */
char refparam[ DAT__SZNAM ];/* Name for reference NDF parameter */
const char *cval; /* The IMPORTSKY string value */
dim_t bolostep; /* Number of bolos per thread */
dim_t idx; /* Index within subgroup */
dim_t nbolo; /* Number of bolometers */
dim_t ntslice; /* Number of time slices */
double *ptr; /* Pointer to NDF Data array */
double *vptr; /* Pointer to NDF Variance array */
int indf1; /* Id. for supplied reference NDF */
int indf2; /* Id. for used section of reference NDF */
int iw; /* Thread index */
int nel; /* Number of mapped NDF pixels */
int nw; /* Number of worker threads */
int result; /* Returned flag */
int there; /* Is there a smurf extension in the NDF? */
size_t bstride; /* bolo stride */
size_t size; /* Size of mapped array */
size_t tstride; /* Time slice stride in data array */
smfArray *qua=NULL; /* Pointer to QUA at chunk */
smf_qfam_t family; /* The family of quality flags found in the NDF */
smf_qual_t *qptr; /* Pointer to mapped quality values */
smf_qual_t *qua_data=NULL; /* Pointer to quality data */

/* Initialise the returned value to indicate no sky has been subtractred. */
result = 0;
Expand Down Expand Up @@ -263,54 +239,6 @@ int smf_initial_sky( ThrWorkForce *wf, AstKeyMap *keymap, smfDIMMData *dat,
/* Remove any existinction correction to the modifed bolometer data. */
if( dat->ext ) smf_calcmodel_ext( wf, dat, 0, keymap, dat->ext,
SMF__DIMM_INVERT, status);

/* Now clear quality flags SMF__Q_COM and SMF__Q_RING in the bolometer data.
This is usually done when undoing the COM model. But the COM model is not
undone at the start of the first iteration, so when running skyloop (which
only performs one iteration on each invocation of makemap) we need
to clear the flags here instead. How many threads do we get to play with */
nw = wf ? wf->nworker : 1;

/* Allocate job data for threads. */
job_data = astCalloc( nw, sizeof(*job_data) );

/* Process each sub-array in turn. */
qua = dat->qua[ 0 ];
for( idx = 0; idx < qua->ndat; idx++ ) {

/* Obtain dimensions of the data */
smf_get_dims( qua->sdata[ idx ], NULL, NULL, &nbolo, &ntslice,
NULL, &bstride, &tstride, status);

/* Get pointers to quality values. */
qua_data = ( qua->sdata[ idx ]->pntr )[0];

/* Find how many bolometers to process in each worker thread. */
bolostep = nbolo/nw;
if( bolostep == 0 ) bolostep = 1;

/* Store information for use by the worker threads. */
for( iw = 0; iw < nw; iw++ ) {
pdata = job_data + iw;
pdata->b1 = iw*bolostep;
if( iw < nw - 1 ) {
pdata->b2 = pdata->b1 + bolostep - 1;
} else {
pdata->b2 = nbolo - 1 ;
}

pdata->ntslice = ntslice;
pdata->qua_data = qua_data;
pdata->bstride = bstride;
pdata->tstride = tstride;
thrAddJob( wf, 0, pdata, smf1_initial_sky, 0, NULL, status );
}

thrWait( wf, status );
}

job_data = astFree( job_data );

}

/* End the AST context. */
Expand All @@ -320,65 +248,3 @@ int smf_initial_sky( ThrWorkForce *wf, AstKeyMap *keymap, smfDIMMData *dat,
return result;
}

static void smf1_initial_sky( void *job_data_ptr, int *status ) {
/*
* Name:
* smf1_initial_sky
* Purpose:
* Executed in a worker thread to do various calculations for
* smf_initial_sky.
* Invocation:
* smf1_initial_sky( void *job_data_ptr, int *status )
* Arguments:
* job_data_ptr = SmfInitialSkyData * (Given)
* Data structure describing the job to be performed by the worker
* thread.
* status = int * (Given and Returned)
* Inherited status.
*/

/* Local Variables: */
SmfInitialSkyData *pdata;
dim_t ibolo;
dim_t itime;
size_t ibase;
smf_qual_t *pq;

/* Check inherited status */
if( *status != SAI__OK ) return;

/* Get a pointer that can be used for accessing the required items in the
supplied structure. */
pdata = (SmfInitialSkyData *) job_data_ptr;

/* Loop round all bolos to be processed by this thread, maintaining the
index of the first time slice for the current bolo. */
ibase = pdata->b1*pdata->bstride;
for( ibolo = pdata->b1; ibolo <= pdata->b2; ibolo++ ) {

/* Check that the whole bolometer has not been flagged as bad. */
pq = pdata->qua_data + ibase;
if( !( *pq & SMF__Q_BADB ) ) {

/* Loop round all time slices. */
for( itime = 0; itime < pdata->ntslice; itime++ ) {

/* Clear any SMF__Q_RING or SMF__Q_COM flags. */
*pq &= ( ~SMF__Q_RING & ~SMF__Q_COM );

/* Move on to the next time slice. */
pq += pdata->tstride;
}
}

/* Increment the index of the first value associated with the next
bolometer. */
ibase += pdata->bstride;
}

}

30 changes: 6 additions & 24 deletions applications/smurf/scripts/skyloop.py
View file
Expand Up @@ -53,7 +53,7 @@
* - First iteration:
* numiter=1
* noi.export=1
* exportNDF=(lut,ext,res,qua)
* exportNDF=(lut,ext)
* noexportsetbad=1
* exportclean=1
* ast.zero_notlast = 0
Expand All @@ -69,7 +69,6 @@
* - Subsequent iterations:
* numiter=1
* noi.import=1
* exportNDF=(res,qua)
* doclean=0
* importsky=ref
* importlut=1
Expand Down Expand Up @@ -280,6 +279,8 @@
* Ensure same map bounds are used on every invocation of makemap.
* 14-FEB-2014 (DSB):
* Ensure downsampling occurs only on the first invocation of makemap.
* 4-MAR-2014 (DSB):
* Do not update quality flags at the end of each iteration.
*-
'''

Expand Down Expand Up @@ -327,8 +328,6 @@ def cleanup():
os.remove( lut )
for noi in new_noi_ndfs:
os.remove( noi )
for res in qua:
os.remove( res )
except:
pass

Expand Down Expand Up @@ -554,9 +553,8 @@ def cleanup():
fd.write("noi.export=1\n") # Export the NOI model. This forces the
# NOI model to be created and exported after
# the first iteration has completed.
fd.write("exportNDF=(lut,ext,res,qua)\n")# Save the EXT, LUT model values to avoid
fd.write("exportNDF=(lut,ext)\n")# Save the EXT, LUT model values to avoid
# re-calculation on each invocation of makemap.
# Also need QUA to update time-series flags
fd.write("noexportsetbad=1\n")# Export good EXT values for bad bolometers
if not precleaned:
fd.write("exportclean=1\n") # Likewise save the cleaned time-series data.
Expand Down Expand Up @@ -670,24 +668,20 @@ def cleanup():
elif os.stat(ndf).st_atime > orig_noi_ndfs[ndf]:
new_noi_ndfs.append(ndf)

# Get the paths to the the moved cleaned files. Also get the paths to the
# files holding the quality flags at the end of each invocation of
# makemap.
# Get the paths to the the moved cleaned files.
if niter > 1:
if not precleaned:
cleaned = NDG( os.path.join( NDG.tempdir,"s*_con_res_cln.sdf"))
qua = NDG( cleaned, "./*|_cln||" )
else:
cleaned = indata
qua = None

# Now do the second and subsequent iterations. These use the cleaned
# time-series data created by the first iteration as their time-series
# input, and use the output map from the previous iteration as their
# initial guess at the sky. First create a map holding things to add
# to the config for subsequent invocations.
add = {}
add["exportNDF"] = "(res,qua)" # Prevent EXT or LUT model being exported.
add["exportNDF"] = 0 # Prevent EXT or LUT model being exported.
add["exportclean"] = 0 # Prevent cleaned time-series data being exported.
add["doclean"] = 0 # Do not clean the supplied data (it has
add["importsky"] = "ref" # Get the initial sky estimate from the REF parameter.
Expand Down Expand Up @@ -793,9 +787,6 @@ def cleanup():
add["com.perarray_last"] = com_perarray_last
newcon = 1

# No need to export quality flags on the last iteration.
add["exportNDF"] = 0

# If this is not the last iteration, get the name of a temporary NDF that
# can be used to store the current iteration's map. This NDF is put in
# the NDG temp directory.
Expand All @@ -813,11 +804,6 @@ def cleanup():
fd.write("{0}={1}\n".format( key, add[key] ))
fd.close()

# Update the quality flags in the cleaned time-series data to be the
# same as the flags exported at the end of the previous iteration.
if qua != None:
invoke( "$KAPPA_DIR/setqual ndf={0} like={1}".format(cleaned,qua) )

# See if the output NDF already exists.
gotit = False
if restart != None:
Expand Down Expand Up @@ -867,10 +853,6 @@ def cleanup():
# itermap cube.
maps.append(newmap)

# Update the NDF from which new quality info is to be read.
if qua:
qua = NDG( cleaned, "./*_con_res" )

# Increment the iteration number
iter += 1

Expand Down

2 comments on commit 90d7ba5

@timj
Copy link
Member

@timj timj commented on 90d7ba5 Mar 5, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@dsberry So your problem is that Skyloop may need to keep more than 8 bits of quality around but that's not possible? I must admit to that not occurring to me. Obviously the fix is to update NDF to support arbitrary bit depth for QUALITY...

@dsberry
Copy link
Member

@dsberry dsberry commented on 90d7ba5 Mar 5, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change was only partly forced by the limit of 8 quality bits. I also realised there was no need to feed the output time-series qualities back into the input on the next iteration anyway. So at the moment, the restriction to 8 bits is not an issue.

Please sign in to comment.