Skip to content

Commit

Permalink
smurf: Allow a diagnostic map to be created from any model
Browse files Browse the repository at this point in the history 
So you can get a 2D map of the COM or FLT model for instance.
  • Loading branch information
David Berry committed Jul 16, 2014
1 parent 1198c43 commit 107a794
Show file tree
Hide file tree
Showing 4 changed files with 146 additions and 13 deletions.
18 changes: 18 additions & 0 deletions applications/smurf/defaults/smurf_makemap.def
View file
Expand Up @@ -685,6 +685,24 @@ diag.bolo = WMEAN
# Only used if a value is supplied for diag.out.
diag.cube = 0

#-----------------------------------------------------------------
# Name: diag.map
# Type: int
# Default: 0
# Purpose: Qualifies the required dignostic information.
# Description: If non-zero, a 2D map containing the binned time-stream
# data for all bolometers is created at each iteration, for
# each required model and set of residuals. These are placed
# in 2D NDFs with names in the following format:
# "<where>_<chunk>_map_<it>", where <chunk> and <where> are
# described above (under "OUT") and <it> is the iteration
# number. If the "diag.map" value is positive, the map will
# contain data for just the subarray specified by "diag.array".
# If "diag.map" is negative, the map will contain data for all
# available subarrays. Only used if a value is supplied for
# "diag.out".
diag.map = 0

#-----------------------------------------------------------------
# Name: diag.mask
# Type: int
Expand Down
2 changes: 1 addition & 1 deletion applications/smurf/libsmf/smf.h.source
View file
Expand Up @@ -962,7 +962,7 @@ void smf_diag( ThrWorkForce *wf, HDSLoc *loc, int *ibolo, int irow,
int power, int time, int isub, smfDIMMData *dat,
smf_modeltype type, smfArray *model, int res,
const char *root, int mask, double mingood, int cube,
int addqual, int *status );
int map, int addqual, int *status );

void smf_diagnostics( ThrWorkForce *wf, int where, smfDIMMData *dat,
int chunk, AstKeyMap *keymap, smfArray **allmodel,
Expand Down
116 changes: 108 additions & 8 deletions applications/smurf/libsmf/smf_diag.c
View file
Expand Up @@ -17,7 +17,7 @@
* int power, int time, int isub, smfDIMMData *dat,
* smf_modeltype type, smfArray *model, int res,
* const char *root, int mask, double mingood, int cube,
* int addqual, int *status )
* int map, int addqual, int *status )
* Arguments:
* wf = ThrWorkForce * (Given)
Expand Down Expand Up @@ -67,6 +67,12 @@
* Is a full cube containing time ordered data for all bolometers
* required? If so, it will be stored in an NDF with name
* "<root>_cube_<irow>".
* map = int (Given)
* Is a 2D map containing the binned time-stream data for all
* bolometers required? If so, it will be stored in an NDF with name
* "<root>_map_<irow>". If "map" is positive, the map will contain
* data for just the subarray specified by "isub".If "map" is
* negative, the map will contain data for all available subarrays.
* addqual = int (Given)
* If non-zero, the output NDFs will include a quality array.
* status = int* (Given and Returned)
Expand Down Expand Up @@ -152,7 +158,7 @@ void smf_diag( ThrWorkForce *wf, HDSLoc *loc, int *ibolo, int irow,
int power, int time, int isub, smfDIMMData *dat,
smf_modeltype type, smfArray *model, int res,
const char *root, int mask, double mingood, int cube,
int addqual, int *status ){
int map, int addqual, int *status ){

/* Local Variables: */
AstCmpFrame *totfrm;
Expand All @@ -170,23 +176,33 @@ void smf_diag( ThrWorkForce *wf, HDSLoc *loc, int *ibolo, int irow,
dim_t fdims[2];
dim_t itime;
dim_t jbolo;
dim_t ndata;
dim_t nbolo;
dim_t nbolor;
dim_t ndata;
dim_t ngood;
dim_t ntslice;
double *buffer;
double *ip;
double *ipv;
double *oldcom;
double *oldres = NULL;
double *pd;
double *var;
double *wf_map;
double *wf_mapwgt;
double *wf_mapwgtsq;
double *wf_mapvar;
double scalevar;
int *index;
int *wf_hitsmap;
int bolostep;
int el;
int fax;
int i;
int iax;
int idx;
int idxhi;
int idxlo;
int indf;
int iw;
int lbnd[ 2 ];
Expand All @@ -195,27 +211,28 @@ void smf_diag( ThrWorkForce *wf, HDSLoc *loc, int *ibolo, int irow,
int ndim;
int nw;
int place;
int rebinflags;
int sorted;
int timestep;
int ubnd[ 2 ];
int usebolo;
size_t bstride;
size_t tstride;
size_t bstrider;
size_t tstrider;
size_t nmap;
size_t tstride;
size_t tstrider;
smfArray *array;
smfData *data = NULL;
smfData *data_tmp;
smfData *pow;
smfData *sidequal;
smf_qual_t *oldcomq;
smf_qual_t *pqr;
smf_qual_t *pq;
smf_qual_t *pqr;
smf_qual_t *qbuffer = NULL;
smf_qual_t *qua = NULL;
smf_qual_t *qual;
smf_qual_t *qbuffer = NULL;
smf_qual_t qval;
smfData *sidequal;

/* Check inherited status. */
if( *status != SAI__OK ) return;
Expand Down Expand Up @@ -738,6 +755,89 @@ void smf_diag( ThrWorkForce *wf, HDSLoc *loc, int *ibolo, int irow,

}

/* If required, create a 2D map containing the binned data from all
bolometers. */
if( map && data->ndims == 3 && (data->pntr)[0] ) {

/* Get the name of the NDF to hold the map. */
name = NULL;
nc = 0;
name = astAppendString( name, &nc, root );
name = astAppendString( name, &nc, "_map_" );
sprintf( attr, "%d", irow );
name = astAppendString( name, &nc, attr );

/* Create it and map the Data and Variance components. */
ndfPlace( loc, name, &place, status );
ndfNew( "_DOUBLE", 2, dat->lbnd_out, dat->ubnd_out, &place, &indf,
status );
ndfMap( indf, "DATA", "_DOUBLE", "WRITE", (void **) &ip, &el,
status );
ndfMap( indf, "VARIANCE", "_DOUBLE", "WRITE", (void **) &ipv, &el,
status );

/* If we are dumping the AST model, just copy the existing map. */
if( !array && *status == SAI__OK ) {
memcpy( ip, dat->map, dat->msize*sizeof( *dat->map ) );
memcpy( ipv, dat->mapvar, dat->msize*sizeof( *dat->mapvar ) );

/* Otherwise, bin the time-stream data to form the map. */
} else {

/* Get the indices of the first and last subarray to include in the map. */
if( map > 0 ) {
idxlo = isub;
idxhi = isub;
} else {
idxlo = 0;
idxhi = array->ndat - 1;
}

/* Allocate memory for multithreaded maps. */
wf_map = astMalloc( nw*dat->msize*sizeof( *wf_map ) );
wf_mapwgt = astMalloc( nw*dat->msize*sizeof( *wf_mapwgt ) );
wf_mapwgtsq = astMalloc( nw*dat->msize*sizeof( *wf_mapwgtsq ) );
wf_mapvar = astMalloc( nw*dat->msize*sizeof( *wf_mapvar ) );
wf_hitsmap = astMalloc( nw*dat->msize*sizeof( *wf_hitsmap ) );

/* Loop over subarray. */
for( idx = idxlo; idx <= idxhi; idx++ ) {

/* initialise rebin flags. */
rebinflags = 0;

/* First call to rebin clears the arrays */
if( idx == idxlo ) rebinflags = rebinflags | AST__REBININIT;

/* Final call to rebin re-normalizes */
if( idx == idxhi ) rebinflags = rebinflags | AST__REBINEND;

/* Rebin the residual + astronomical signal into a map */
smf_rebinmap1( wf, array->sdata[ idx ], NULL,
dat->lut[0]->sdata[idx]->pntr[0], 0, 0, 0,
NULL, 0, SMF__Q_GOOD, 1, rebinflags,
wf_map, wf_mapwgt, wf_mapwgtsq, wf_hitsmap,
wf_mapvar, dat->msize, &scalevar, status );
}

/* Copy the data and variance arrays to the NDF. */
if( *status == SAI__OK ) {
memcpy( ip, wf_map, dat->msize*sizeof( *wf_map ) );
memcpy( ipv, wf_mapvar, dat->msize*sizeof( *wf_mapvar ) );
}

/* Free memory. */
wf_map = astFree( wf_map );
wf_mapwgt = astFree( wf_mapwgt );
wf_mapwgtsq = astFree( wf_mapwgtsq );
wf_mapvar = astFree( wf_mapvar );
wf_hitsmap = astFree( wf_hitsmap );
}

/* Annul the NDF identifier. */
ndfAnnul( &indf, status );
}

/* Free the array holding the AST model and re-instate the original RES
values if required. */
if( oldres ) {
Expand Down
23 changes: 19 additions & 4 deletions applications/smurf/libsmf/smf_diagnostics.c
View file
Expand Up @@ -67,6 +67,15 @@
* <chunk> and <where> are described above (under "OUT") and <it>
* is the iteration number.
*
* MAP - If non-zero, a 2D map containing the binned time-stream data
* for all bolometers is created at each iteration, for each required
* model and set of residuals. These are placed in 2D NDFs with
* names in the following format: "<where>_<chunk>_map_<it>", where
* <chunk> and <where> are described above (under "OUT") and <it>
* is the iteration number. If the MAP value is positive, the map will
* contain data for just the subarray specified by "ARRAY".If MAP is
* negative, the map will contain data for all available subarrays.
*
* MODELS - Indicates the models that are to be written out. It
* should be a comma separated list of model names (e.g. COM, FLT,
* AST, RES, etc) contained within parentheses, or a single model
Expand Down Expand Up @@ -143,9 +152,11 @@
* History:
* 25-JAN-2013 (DSB):
* Original version.
* 17-JUL-2014 (DSB):
* Added option to create diagnostic maps from each model.
* Copyright:
* Copyright (C) 2013 Science and Technology Facilities Council.
* Copyright (C) 2013-2014 Science and Technology Facilities Council.
* All Rights Reserved.
* Licence:
Expand Down Expand Up @@ -207,6 +218,7 @@ void smf_diagnostics( ThrWorkForce *wf, int where, smfDIMMData *dat,
int irow;
int isub;
int ivals[ 2 ];
int map;
int mask;
int new;
int nmodel;
Expand Down Expand Up @@ -251,6 +263,9 @@ void smf_diagnostics( ThrWorkForce *wf, int where, smfDIMMData *dat,
each iteration. */
astMapGet0I( kmap, "CUBE", &cube );

/* See if a 2D map of each model is required at each iteration. */
astMapGet0I( kmap, "MAP", &map );

/* See if NDFs are to include Quality arrays. */
astMapGet0I( kmap, "QUAL", &addqual );

Expand Down Expand Up @@ -422,7 +437,7 @@ void smf_diagnostics( ThrWorkForce *wf, int where, smfDIMMData *dat,
sprintf( root, "bef_%d", chunk );
smf_diag( wf, mloc, &ibolo, irow, power, time, isub,
dat, type, NULL, 1, root, 0, mingood, cube,
addqual, status );
map, addqual, status );
}

/* If this function has been called immediately after estimating the new
Expand All @@ -432,14 +447,14 @@ void smf_diagnostics( ThrWorkForce *wf, int where, smfDIMMData *dat,
sprintf( root, "mod_%d", chunk );
smf_diag( wf, mloc, &ibolo, irow, power, time, isub,
dat, type, allmodel ? allmodel[ 0 ] : NULL,
0, root, mask, mingood, cube, addqual, status );
0, root, mask, mingood, cube, map, addqual, status );
if( res_after && type != SMF__RES ) {
msgOutf( "", "Diagnostics: Dumping residuals after subtraction of %s",
status, modname );
sprintf( root, "aft_%d", chunk );
smf_diag( wf, mloc, &ibolo, irow, power, time, isub,
dat, type, NULL, 1, root, 0, mingood, cube,
addqual, status );
map, addqual, status );
}

/* Any other "where" value is currently an error. */
Expand Down

0 comments on commit 107a794

Please sign in to comment.