From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- pkg/tbtables/cfitsio/eval_f.c | 2293 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2293 insertions(+) create mode 100644 pkg/tbtables/cfitsio/eval_f.c (limited to 'pkg/tbtables/cfitsio/eval_f.c') diff --git a/pkg/tbtables/cfitsio/eval_f.c b/pkg/tbtables/cfitsio/eval_f.c new file mode 100644 index 00000000..20441818 --- /dev/null +++ b/pkg/tbtables/cfitsio/eval_f.c @@ -0,0 +1,2293 @@ +/************************************************************************/ +/* */ +/* CFITSIO Lexical Parser */ +/* */ +/* This file is one of 3 files containing code which parses an */ +/* arithmetic expression and evaluates it in the context of an input */ +/* FITS file table extension. The CFITSIO lexical parser is divided */ +/* into the following 3 parts/files: the CFITSIO "front-end", */ +/* eval_f.c, contains the interface between the user/CFITSIO and the */ +/* real core of the parser; the FLEX interpreter, eval_l.c, takes the */ +/* input string and parses it into tokens and identifies the FITS */ +/* information required to evaluate the expression (ie, keywords and */ +/* columns); and, the BISON grammar and evaluation routines, eval_y.c, */ +/* receives the FLEX output and determines and performs the actual */ +/* operations. The files eval_l.c and eval_y.c are produced from */ +/* running flex and bison on the files eval.l and eval.y, respectively. */ +/* (flex and bison are available from any GNU archive: see www.gnu.org) */ +/* */ +/* The grammar rules, rather than evaluating the expression in situ, */ +/* builds a tree, or Nodal, structure mapping out the order of */ +/* operations and expression dependencies. This "compilation" process */ +/* allows for much faster processing of multiple rows. This technique */ +/* was developed by Uwe Lammers of the XMM Science Analysis System, */ +/* although the CFITSIO implementation is entirely code original. */ +/* */ +/* */ +/* Modification History: */ +/* */ +/* Kent Blackburn c1992 Original parser code developed for the */ +/* FTOOLS software package, in particular, */ +/* the fselect task. */ +/* Kent Blackburn c1995 BIT column support added */ +/* Peter D Wilson Feb 1998 Vector column support added */ +/* Peter D Wilson May 1998 Ported to CFITSIO library. User */ +/* interface routines written, in essence */ +/* making fselect, fcalc, and maketime */ +/* capabilities available to all tools */ +/* via single function calls. */ +/* Peter D Wilson Jun 1998 Major rewrite of parser core, so as to */ +/* create a run-time evaluation tree, */ +/* inspired by the work of Uwe Lammers, */ +/* resulting in a speed increase of */ +/* 10-100 times. */ +/* Peter D Wilson Jul 1998 gtifilter(a,b,c,d) function added */ +/* Peter D Wilson Aug 1998 regfilter(a,b,c,d) function added */ +/* Peter D Wilson Jul 1999 Make parser fitsfile-independent, */ +/* allowing a purely vector-based usage */ +/* Peter D Wilson Aug 1999 Add row-offset capability */ +/* Peter D Wilson Sep 1999 Add row-range capability to ffcalc_rng */ +/* */ +/************************************************************************/ + +#include +#include +#include "eval_defs.h" +#include "region.h" + +typedef struct { + int datatype; /* Data type to cast parse results into for user */ + void *dataPtr; /* Pointer to array of results, NULL if to use iterCol */ + void *nullPtr; /* Pointer to nulval, use zero if NULL */ + long maxRows; /* Max No. of rows to process, -1=all, 0=1 iteration */ + int anyNull; /* Flag indicating at least 1 undef value encountered */ +} parseInfo; + +/* Internal routines needed to allow the evaluator to operate on FITS data */ + +static void Setup_DataArrays( int nCols, iteratorCol *cols, + long fRow, long nRows ); +static int find_column( char *colName, void *itslval ); +static int find_keywd ( char *key, void *itslval ); +static int allocateCol( int nCol, int *status ); +static int load_column( int varNum, long fRow, long nRows, + void *data, char *undef ); + +/*---------------------------------------------------------------------------*/ +int fffrow( fitsfile *fptr, /* I - Input FITS file */ + char *expr, /* I - Boolean expression */ + long firstrow, /* I - First row of table to eval */ + long nrows, /* I - Number of rows to evaluate */ + long *n_good_rows, /* O - Number of rows eval to True */ + char *row_status, /* O - Array of boolean results */ + int *status ) /* O - Error status */ +/* */ +/* Evaluate a boolean expression using the indicated rows, returning an */ +/* array of flags indicating which rows evaluated to TRUE/FALSE */ +/*---------------------------------------------------------------------------*/ +{ + parseInfo Info; + int naxis, constant; + long nelem, naxes[MAXDIMS], elem; + char result; + + if( *status ) return( *status ); + + if( ffiprs( fptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis, + naxes, status ) ) { + ffcprs(); + return( *status ); + } + if( nelem<0 ) { + constant = 1; + nelem = -nelem; + } else + constant = 0; + + if( Info.datatype!=TLOGICAL || nelem!=1 ) { + ffcprs(); + ffpmsg("Expression does not evaluate to a logical scalar."); + return( *status = PARSE_BAD_TYPE ); + } + + if( constant ) { /* No need to call parser... have result from ffiprs */ + result = gParse.Nodes[gParse.resultNode].value.data.log; + *n_good_rows = nrows; + for( elem=0; elem1 ? firstrow : 1); + Info.dataPtr = row_status; + Info.nullPtr = NULL; + Info.maxRows = nrows; + + if( ffiter( gParse.nCols, gParse.colData, firstrow-1, 0, + parse_data, (void*)&Info, status ) == -1 ) + *status = 0; /* -1 indicates exitted without error before end... OK */ + + if( *status ) { + + /***********************/ + /* Error... Do nothing */ + /***********************/ + + } else { + + /***********************************/ + /* Count number of good rows found */ + /***********************************/ + + *n_good_rows = 0L; + for( elem=0; elemHDUposition != (infptr->Fptr)->curhdu ) + ffmahd( infptr, (infptr->HDUposition) + 1, NULL, status ); + if( *status ) { + ffcprs(); + return( *status ); + } + inExt.rowLength = (long) (infptr->Fptr)->rowlength; + inExt.numRows = (infptr->Fptr)->numrows; + inExt.heapSize = (infptr->Fptr)->heapsize; + if( inExt.numRows == 0 ) { /* Nothing to copy */ + ffcprs(); + return( *status ); + } + + if( outfptr->HDUposition != (outfptr->Fptr)->curhdu ) + ffmahd( outfptr, (outfptr->HDUposition) + 1, NULL, status ); + if( (outfptr->Fptr)->datastart < 0 ) + ffrdef( outfptr, status ); + if( *status ) { + ffcprs(); + return( *status ); + } + outExt.rowLength = (long) (outfptr->Fptr)->rowlength; + outExt.numRows = (outfptr->Fptr)->numrows; + if( !outExt.numRows ) + (outfptr->Fptr)->heapsize = 0L; + outExt.heapSize = (outfptr->Fptr)->heapsize; + + if( inExt.rowLength != outExt.rowLength ) { + ffpmsg("Output table has different row length from input"); + ffcprs(); + return( *status = PARSE_BAD_OUTPUT ); + } + + /***********************************/ + /* Fill out Info data for parser */ + /***********************************/ + + Info.dataPtr = (char *)malloc( (inExt.numRows + 1) * sizeof(char) ); + Info.nullPtr = NULL; + Info.maxRows = inExt.numRows; + if( !Info.dataPtr ) { + ffpmsg("Unable to allocate memory for row selection"); + ffcprs(); + return( *status = MEMORY_ALLOCATION ); + } + + /* make sure array is zero terminated */ + ((char*)Info.dataPtr)[inExt.numRows] = 0; + + if( constant ) { /* Set all rows to the same value from constant result */ + + result = gParse.Nodes[gParse.resultNode].value.data.log; + for( ntodo = 0; ntodo 1) + ffirow( outfptr, outExt.numRows, nGood, status ); + } + + do { + if( ((char*)Info.dataPtr)[inloc-1] ) { + ffgtbb( infptr, inloc, 1L, rdlen, buffer+rdlen*nbuff, status ); + nbuff++; + if( nbuff==maxrows ) { + ffptbb( outfptr, outloc, 1L, rdlen*nbuff, buffer, status ); + outloc += nbuff; + nbuff = 0; + } + } + inloc++; + } while( !*status && inloc<=inExt.numRows ); + + if( nbuff ) { + ffptbb( outfptr, outloc, 1L, rdlen*nbuff, buffer, status ); + outloc += nbuff; + } + + if( infptr==outfptr ) { + + if( outloc<=inExt.numRows ) + ffdrow( infptr, outloc, inExt.numRows-outloc+1, status ); + + } else if( inExt.heapSize && nGood ) { + + /* Copy heap, if it exists and at least one row copied */ + + /********************************************************/ + /* Get location information from the output extension */ + /********************************************************/ + + if( outfptr->HDUposition != (outfptr->Fptr)->curhdu ) + ffmahd( outfptr, (outfptr->HDUposition) + 1, NULL, status ); + outExt.dataStart = (outfptr->Fptr)->datastart; + outExt.heapStart = (outfptr->Fptr)->heapstart; + + /*************************************************/ + /* Insert more space into outfptr if necessary */ + /*************************************************/ + + hsize = outExt.heapStart + outExt.heapSize; + freespace = ( ( (hsize + 2879) / 2880) * 2880) - hsize; + ntodo = inExt.heapSize; + + if ( (freespace - ntodo) < 0) { /* not enough existing space? */ + ntodo = (ntodo - freespace + 2879) / 2880; /* number of blocks */ + ffiblk(outfptr, ntodo, 1, status); /* insert the blocks */ + } + ffukyj( outfptr, "PCOUNT", inExt.heapSize+outExt.heapSize, + NULL, status ); + + /*******************************************************/ + /* Get location information from the input extension */ + /*******************************************************/ + + if( infptr->HDUposition != (infptr->Fptr)->curhdu ) + ffmahd( infptr, (infptr->HDUposition) + 1, NULL, status ); + inExt.dataStart = (infptr->Fptr)->datastart; + inExt.heapStart = (infptr->Fptr)->heapstart; + + /**********************************/ + /* Finally copy heap to outfptr */ + /**********************************/ + + ntodo = inExt.heapSize; + inbyteloc = inExt.heapStart + inExt.dataStart; + outbyteloc = outExt.heapStart + outExt.dataStart + outExt.heapSize; + + while ( ntodo && !*status ) { + rdlen = minvalue(ntodo,500000); + ffmbyt( infptr, inbyteloc, REPORT_EOF, status ); + ffgbyt( infptr, rdlen, buffer, status ); + ffmbyt( outfptr, outbyteloc, IGNORE_EOF, status ); + ffpbyt( outfptr, rdlen, buffer, status ); + inbyteloc += rdlen; + outbyteloc += rdlen; + ntodo -= rdlen; + } + + /***********************************************************/ + /* But must update DES if data is being appended to a */ + /* pre-existing heap space. Edit each new entry in file */ + /***********************************************************/ + + if( outExt.heapSize ) { + long repeat, offset, j; + int i; + for( i=1; i<=(outfptr->Fptr)->tfield; i++ ) { + if( (outfptr->Fptr)->tableptr[i-1].tdatatype<0 ) { + for( j=outExt.numRows+1; j<=outExt.numRows+nGood; j++ ) { + ffgdes( outfptr, i, j, &repeat, &offset, status ); + offset += outExt.heapSize; + ffpdes( outfptr, i, j, repeat, offset, status ); + } + } + } + } + + } /* End of HEAP copy */ + + free(buffer); + } + + free(Info.dataPtr); + ffcprs(); + + ffcmph(outfptr, status); /* compress heap, deleting any orphaned data */ + return(*status); +} + +/*---------------------------------------------------------------------------*/ +int ffcrow( fitsfile *fptr, /* I - Input FITS file */ + int datatype, /* I - Datatype to return results as */ + char *expr, /* I - Arithmetic expression */ + long firstrow, /* I - First row to evaluate */ + long nelements, /* I - Number of elements to return */ + void *nulval, /* I - Ptr to value to use as UNDEF */ + void *array, /* O - Array of results */ + int *anynul, /* O - Were any UNDEFs encountered? */ + int *status ) /* O - Error status */ +/* */ +/* Calculate an expression for the indicated rows of a table, returning */ +/* the results, cast as datatype (TSHORT, TDOUBLE, etc), in array. If */ +/* nulval==NULL, UNDEFs will be zeroed out. For vector results, the number */ +/* of elements returned may be less than nelements if nelements is not an */ +/* even multiple of the result dimension. Call fftexp to obtain the */ +/* dimensions of the results. */ +/*---------------------------------------------------------------------------*/ +{ + parseInfo Info; + int naxis; + long nelem1, naxes[MAXDIMS]; + + if( *status ) return( *status ); + + if( ffiprs( fptr, 0, expr, MAXDIMS, &Info.datatype, &nelem1, &naxis, + naxes, status ) ) { + ffcprs(); + return( *status ); + } + if( nelem1<0 ) nelem1 = - nelem1; + + if( nelements1 ? firstrow : 1); + + if( datatype ) Info.datatype = datatype; + + Info.dataPtr = array; + Info.nullPtr = nulval; + Info.maxRows = nelements / nelem1; + + if( ffiter( gParse.nCols, gParse.colData, firstrow-1, 0, + parse_data, (void*)&Info, status ) == -1 ) + *status=0; /* -1 indicates exitted without error before end... OK */ + + *anynul = Info.anyNull; + ffcprs(); + return( *status ); +} + +/*--------------------------------------------------------------------------*/ +int ffcalc( fitsfile *infptr, /* I - Input FITS file */ + char *expr, /* I - Arithmetic expression */ + fitsfile *outfptr, /* I - Output fits file */ + char *parName, /* I - Name of output parameter */ + char *parInfo, /* I - Extra information on parameter */ + int *status ) /* O - Error status */ +/* */ +/* Evaluate an expression for all rows of a table. Call ffcalc_rng with */ +/* a row range of 1-MAX. */ +{ + long start=1, end=LONG_MAX; + + return ffcalc_rng( infptr, expr, outfptr, parName, parInfo, + 1, &start, &end, status ); +} + +/*--------------------------------------------------------------------------*/ +int ffcalc_rng( fitsfile *infptr, /* I - Input FITS file */ + char *expr, /* I - Arithmetic expression */ + fitsfile *outfptr, /* I - Output fits file */ + char *parName, /* I - Name of output parameter */ + char *parInfo, /* I - Extra information on parameter */ + int nRngs, /* I - Row range info */ + long *start, /* I - Row range info */ + long *end, /* I - Row range info */ + int *status ) /* O - Error status */ +/* */ +/* Evaluate an expression using the data in the input FITS file and place */ +/* the results into either a column or keyword in the output fits file, */ +/* depending on the value of parName (keywords normally prefixed with '#') */ +/* and whether the expression evaluates to a constant or a table column. */ +/* The logic is as follows: */ +/* (1) If a column exists with name, parName, put results there. */ +/* (2) If parName starts with '#', as in #NAXIS, put result there, */ +/* with parInfo used as the comment. If expression does not evaluate */ +/* to a constant, flag an error. */ +/* (3) If a keyword exists with name, parName, and expression is a */ +/* constant, put result there, using parInfo as the new comment. */ +/* (4) Else, create a new column with name parName and TFORM parInfo. */ +/* If parInfo is NULL, use a default data type for the column. */ +/*--------------------------------------------------------------------------*/ +{ + parseInfo Info; + int naxis, constant, typecode, newNullKwd=0; + long nelem, naxes[MAXDIMS], repeat, width; + int col_cnt, colNo; + Node *result; + char card[81], tform[16], nullKwd[9], tdimKwd[9]; + int hdutype; + + if( *status ) return( *status ); + + if( ffiprs( infptr, 0, expr, MAXDIMS, &Info.datatype, &nelem, &naxis, + naxes, status ) ) { + ffcprs(); + return( *status ); + } + if( nelem<0 ) { + constant = 1; + nelem = -nelem; + } else + constant = 0; + + /* Case (1): If column exists put it there */ + + colNo = 0; + if( ffgcno( outfptr, CASEINSEN, parName, &colNo, status )==COL_NOT_FOUND ) { + + /* Output column doesn't exist. Test for keyword. */ + + /* Case (2): Does parName indicate result should be put into keyword */ + + *status = 0; + if( parName[0]=='#' ) { + if( ! constant ) { + ffcprs(); + ffpmsg( "Cannot put tabular result into keyword (ffcalc)" ); + return( *status = PARSE_BAD_TYPE ); + } + parName++; + + } else if( constant ) { + + /* Case (3): Does a keyword named parName already exist */ + + if( ffgcrd( outfptr, parName, card, status )==KEY_NO_EXIST ) { + colNo = -1; + } else if( *status ) { + ffcprs(); + return( *status ); + } + + } else + colNo = -1; + + if( colNo<0 ) { + + /* Case (4): Create new column */ + + *status = 0; + ffgncl( outfptr, &colNo, status ); + colNo++; + ffghdt( outfptr, &hdutype, status ); + if( parInfo==NULL || *parInfo=='\0' ) { + /* Figure out best default column type */ + if( hdutype==BINARY_TBL ) { + sprintf(tform,"%ld",nelem); + switch( Info.datatype ) { + case TLOGICAL: strcat(tform,"L"); break; + case TLONG: strcat(tform,"J"); break; + case TDOUBLE: strcat(tform,"D"); break; + case TSTRING: strcat(tform,"A"); break; + case TBIT: strcat(tform,"X"); break; + } + } else { + switch( Info.datatype ) { + case TLOGICAL: + ffcprs(); + ffpmsg("Cannot create LOGICAL column in ASCII table"); + return( *status = NOT_BTABLE ); + break; + case TLONG: strcpy(tform,"I11"); break; + case TDOUBLE: strcpy(tform,"D23.15"); break; + case TSTRING: + case TBIT: sprintf(tform,"A%ld",nelem); break; + } + } + parInfo = tform; + } else if( !(isdigit((int) *parInfo)) && hdutype==BINARY_TBL ) { + if( Info.datatype==TBIT && *parInfo=='B' ) + nelem = (nelem+7)/8; + sprintf(tform,"%ld%s",nelem,parInfo); + parInfo = tform; + } + fficol( outfptr, colNo, parName, parInfo, status ); + if( naxis>1 ) + ffptdm( outfptr, colNo, naxis, naxes, status ); + + /* Setup TNULLn keyword in case NULLs are encountered */ + + ffkeyn("TNULL", colNo, nullKwd, status); + if( ffgcrd( outfptr, nullKwd, card, status )==KEY_NO_EXIST ) { + *status = 0; + if( hdutype==BINARY_TBL ) { + long nullVal=0; + fits_binary_tform( parInfo, &typecode, &repeat, &width, status ); + if( typecode==TBYTE ) + nullVal = UCHAR_MAX; + else if( typecode==TSHORT ) + nullVal = SHRT_MIN; + else if( typecode==TINT ) + nullVal = INT_MIN; + else if( typecode==TLONG ) + nullVal = LONG_MIN; + if( nullVal ) { + ffpkyj( outfptr, nullKwd, nullVal, "Null value", status ); + fits_set_btblnull( outfptr, colNo, nullVal, status ); + newNullKwd = 1; + } + } else if( hdutype==ASCII_TBL ) { + ffpkys( outfptr, nullKwd, "NULL", "Null value string", status ); + fits_set_atblnull( outfptr, colNo, "NULL", status ); + newNullKwd = 1; + } + } + + } + + } else if( *status ) { + ffcprs(); + return( *status ); + } else { + + /********************************************************/ + /* Check if a TDIM keyword should be written/updated. */ + /********************************************************/ + + ffkeyn("TDIM", colNo, tdimKwd, status); + ffgcrd( outfptr, tdimKwd, card, status ); + if( *status==0 ) { + /* TDIM exists, so update it with result's dimension */ + ffptdm( outfptr, colNo, naxis, naxes, status ); + } else if( *status==KEY_NO_EXIST ) { + /* TDIM does not exist, so clear error stack and */ + /* write a TDIM only if result is multi-dimensional */ + *status = 0; + ffcmsg(); + if( naxis>1 ) + ffptdm( outfptr, colNo, naxis, naxes, status ); + } + if( *status ) { + /* Either some other error happened in ffgcrd */ + /* or one happened in ffptdm */ + ffcprs(); + return( *status ); + } + + } + + if( colNo>0 ) { + + /* Output column exists (now)... put results into it */ + + int anyNull = 0; + int nPerLp, i; + long totaln; + + ffgkyj(infptr, "NAXIS2", &totaln, 0, status); + + /*************************************/ + /* Create new iterator Output Column */ + /*************************************/ + + col_cnt = gParse.nCols; + if( allocateCol( col_cnt, status ) ) { + ffcprs(); + return( *status ); + } + + fits_iter_set_by_num( gParse.colData+col_cnt, outfptr, + colNo, 0, OutputCol ); + gParse.nCols++; + + for( i=0; i= 10) && (nRngs == 1) && + (start[0] == 1) && (end[0] == totaln)) + nPerLp = 0; + else + nPerLp = Info.maxRows; + + if( ffiter( gParse.nCols, gParse.colData, start[i]-1, + nPerLp, parse_data, (void*)&Info, status ) == -1 ) + *status = 0; + else if( *status ) { + ffcprs(); + return( *status ); + } + if( Info.anyNull ) anyNull = 1; + } + + if( newNullKwd && !anyNull ) { + ffdkey( outfptr, nullKwd, status ); + } + + } else { + + /* Put constant result into keyword */ + + result = gParse.Nodes + gParse.resultNode; + switch( Info.datatype ) { + case TDOUBLE: + ffukyd( outfptr, parName, result->value.data.dbl, 15, + parInfo, status ); + break; + case TLONG: + ffukyj( outfptr, parName, result->value.data.lng, parInfo, status ); + break; + case TLOGICAL: + ffukyl( outfptr, parName, result->value.data.log, parInfo, status ); + break; + case TBIT: + case TSTRING: + ffukys( outfptr, parName, result->value.data.str, parInfo, status ); + break; + } + } + + ffcprs(); + return( *status ); +} + +/*--------------------------------------------------------------------------*/ +int fftexp( fitsfile *fptr, /* I - Input FITS file */ + char *expr, /* I - Arithmetic expression */ + int maxdim, /* I - Max Dimension of naxes */ + int *datatype, /* O - Data type of result */ + long *nelem, /* O - Vector length of result */ + int *naxis, /* O - # of dimensions of result */ + long *naxes, /* O - Size of each dimension */ + int *status ) /* O - Error status */ +/* */ +/* Evaluate the given expression and return information on the result. */ +/*--------------------------------------------------------------------------*/ +{ + ffiprs( fptr, 0, expr, maxdim, datatype, nelem, naxis, naxes, status ); + ffcprs(); + return( *status ); +} + +/*--------------------------------------------------------------------------*/ +int ffiprs( fitsfile *fptr, /* I - Input FITS file */ + int compressed, /* I - Is FITS file hkunexpanded? */ + char *expr, /* I - Arithmetic expression */ + int maxdim, /* I - Max Dimension of naxes */ + int *datatype, /* O - Data type of result */ + long *nelem, /* O - Vector length of result */ + int *naxis, /* O - # of dimensions of result */ + long *naxes, /* O - Size of each dimension */ + int *status ) /* O - Error status */ +/* */ +/* Initialize the parser and determine what type of result the expression */ +/* produces. */ +/*--------------------------------------------------------------------------*/ +{ + Node *result; + int i,lexpr, tstatus = 0; + static iteratorCol dmyCol; + + if( *status ) return( *status ); + + /* Initialize the Parser structure */ + + gParse.def_fptr = fptr; + gParse.compressed = compressed; + gParse.nCols = 0; + gParse.colData = NULL; + gParse.varData = NULL; + gParse.getData = find_column; + gParse.loadData = load_column; + gParse.Nodes = NULL; + gParse.nNodesAlloc= 0; + gParse.nNodes = 0; + gParse.status = 0; + + if( ffgkyj(fptr, "NAXIS2", &gParse.totalRows, 0, &tstatus) ) + { + /* this might be a 1D or null image with no NAXIS2 keyword */ + gParse.totalRows = 0; + } + /* Copy expression into parser... read from file if necessary */ + + if( expr[0]=='@' ) { + if( ffimport_file( expr+1, &gParse.expr, status ) ) return( *status ); + lexpr = strlen(gParse.expr); + } else { + lexpr = strlen(expr); + gParse.expr = (char*)malloc( (2+lexpr)*sizeof(char)); + strcpy(gParse.expr,expr); + } + strcat(gParse.expr + lexpr,"\n"); + gParse.index = 0; + gParse.is_eobuf = 0; + + /* Parse the expression, building the Nodes and determing */ + /* which columns are neded and what data type is returned */ + + ffrestart(NULL); + if( ffparse() ) { + return( *status = PARSE_SYNTAX_ERR ); + } + /* Check results */ + + *status = gParse.status; + if( *status ) return(*status); + + if( !gParse.nNodes ) { + ffpmsg("Blank expression"); + return( *status = PARSE_SYNTAX_ERR ); + } + if( !gParse.nCols ) { + dmyCol.fptr = fptr; /* This allows iterator to know value of */ + gParse.colData = &dmyCol; /* fptr when no columns are referenced */ + } + + result = gParse.Nodes + gParse.resultNode; + + *naxis = result->value.naxis; + *nelem = result->value.nelem; + for( i=0; i<*naxis && ivalue.naxes[i]; + + switch( result->type ) { + case BOOLEAN: + *datatype = TLOGICAL; + break; + case LONG: + *datatype = TLONG; + break; + case DOUBLE: + *datatype = TDOUBLE; + break; + case BITSTR: + *datatype = TBIT; + break; + case STRING: + *datatype = TSTRING; + break; + default: + *datatype = 0; + ffpmsg("Bad return data type"); + *status = gParse.status = PARSE_BAD_TYPE; + break; + } + gParse.datatype = *datatype; + free(gParse.expr); + + if( result->operation==CONST_OP ) *nelem = - *nelem; + return(*status); +} + +/*--------------------------------------------------------------------------*/ +void ffcprs( void ) /* No parameters */ +/* */ +/* Clear the parser, making it ready to accept a new expression. */ +/*--------------------------------------------------------------------------*/ +{ + int col, node, i; + + if( gParse.nCols > 0 ) { + free( gParse.colData ); + for( col=0; col 0 ) { + node = gParse.nNodes; + while( node-- ) { + i = gParse.Nodes[node].SubNodes[0]; + if( gParse.Nodes[node].operation==gtifilt_fct ) { + free( gParse.Nodes[ i ].value.data.ptr ); + } + else if( gParse.Nodes[node].operation==regfilt_fct ) { + fits_free_region( (SAORegion *)gParse.Nodes[ i ].value.data.ptr ); + } + } + gParse.nNodes = 0; + } + if( gParse.Nodes ) free( gParse.Nodes ); + gParse.Nodes = NULL; +} + +/*---------------------------------------------------------------------------*/ +int parse_data( long totalrows, /* I - Total rows to be processed */ + long offset, /* I - Number of rows skipped at start*/ + long firstrow, /* I - First row of this iteration */ + long nrows, /* I - Number of rows in this iter */ + int nCols, /* I - Number of columns in use */ + iteratorCol *colData, /* IO- Column information/data */ + void *userPtr ) /* I - Data handling instructions */ +/* */ +/* Iterator work function which calls the parser and copies the results */ +/* into either an OutputCol or a data pointer supplied in the userPtr */ +/* structure. */ +/*---------------------------------------------------------------------------*/ +{ + int status, constant=0, anyNullThisTime=0; + long jj, kk, idx, remain, ntodo; + Node *result; + + /* declare variables static to preserve their values between calls */ + static void *Data, *Null; + static int datasize; + static long lastRow, jnull, repeat, resDataSize; + static parseInfo *userInfo; + static long zeros[4] = {0,0,0,0}; + + /*--------------------------------------------------------*/ + /* Initialization procedures: execute on the first call */ + /*--------------------------------------------------------*/ + if (firstrow == offset+1) + { + userInfo = (parseInfo*)userPtr; + userInfo->anyNull = 0; + + if( userInfo->maxRows>0 ) + userInfo->maxRows = minvalue(totalrows,userInfo->maxRows); + else if( userInfo->maxRows<0 ) + userInfo->maxRows = totalrows; + else + userInfo->maxRows = nrows; + + lastRow = firstrow + userInfo->maxRows - 1; + + if( userInfo->dataPtr==NULL ) { + + if( colData[nCols-1].iotype == InputCol ) { + ffpmsg("Output column for parser results not found!"); + return( PARSE_NO_OUTPUT ); + } + /* Data gets set later */ + Null = colData[nCols-1].array; + userInfo->datatype = colData[nCols-1].datatype; + + /* Check for a TNULL keyword for output column */ + + status = 0; + jnull = 0L; + ffgknj( colData[nCols-1].fptr, "TNULL", colData[nCols-1].colnum, + 1, &jnull, (int*)&jj, &status ); + if( status==BAD_INTKEY ) { + /* Probably ASCII table with text TNULL keyword */ + switch( userInfo->datatype ) { + case TSHORT: jnull = SHRT_MIN; break; + case TINT: jnull = INT_MIN; break; + case TLONG: jnull = LONG_MIN; break; + } + } + repeat = colData[nCols-1].repeat; + + } else { + + Data = userInfo->dataPtr; + Null = (userInfo->nullPtr ? userInfo->nullPtr : zeros); + repeat = gParse.Nodes[gParse.resultNode].value.nelem; + + } + + /* Determine the size of each element of the returned result */ + + switch( userInfo->datatype ) { + case TBIT: /* Fall through to TBYTE */ + case TLOGICAL: /* Fall through to TBYTE */ + case TBYTE: datasize = sizeof(char); break; + case TSHORT: datasize = sizeof(short); break; + case TINT: datasize = sizeof(int); break; + case TLONG: datasize = sizeof(long); break; + case TFLOAT: datasize = sizeof(float); break; + case TDOUBLE: datasize = sizeof(double); break; + case TSTRING: datasize = sizeof(char*); break; + } + + /* Determine the size of each element of the calculated result */ + /* (only matters for numeric/logical data) */ + + switch( gParse.Nodes[gParse.resultNode].type ) { + case BOOLEAN: resDataSize = sizeof(char); break; + case LONG: resDataSize = sizeof(long); break; + case DOUBLE: resDataSize = sizeof(double); break; + } + } + + /*-------------------------------------------*/ + /* Main loop: process all the rows of data */ + /*-------------------------------------------*/ + + /* If writing to output column, set first element to appropriate */ + /* null value. If no NULLs encounter, zero out before returning. */ + + if( userInfo->dataPtr == NULL ) { + /* First, reset Data pointer to start of output array */ + Data = (char*)colData[nCols-1].array + datasize; + + switch( userInfo->datatype ) { + case TLOGICAL: *(char *)Null = 'U'; break; + case TBYTE: *(char *)Null = (char )jnull; break; + case TSHORT: *(short *)Null = (short)jnull; break; + case TINT: *(int *)Null = (int )jnull; break; + case TLONG: *(long *)Null = (long )jnull; break; + case TFLOAT: *(float *)Null = FLOATNULLVALUE; break; + case TDOUBLE: *(double*)Null = DOUBLENULLVALUE; break; + case TSTRING: (*(char **)Null)[0] = '\1'; + (*(char **)Null)[1] = '\0'; break; + } + } + + /* Alter nrows in case calling routine didn't want to do all rows */ + + nrows = minvalue(nrows,lastRow-firstrow+1); + + Setup_DataArrays( nCols, colData, firstrow, nrows ); + + /* Parser allocates arrays for each column and calculation it performs. */ + /* Limit number of rows processed during each pass to reduce memory */ + /* requirements... In most cases, iterator will limit rows to less */ + /* than 2500 rows per iteration, so this is really only relevant for */ + /* hk-compressed files which must be decompressed in memory and sent */ + /* whole to parse_data in a single iteration. */ + + remain = nrows; + while( remain ) { + + ntodo = minvalue(remain,2500); + Evaluate_Parser ( firstrow, ntodo ); + if( gParse.status ) break; + + firstrow += ntodo; + remain -= ntodo; + + /* Copy results into data array */ + + result = gParse.Nodes + gParse.resultNode; + if( result->operation==CONST_OP ) constant = 1; + + switch( result->type ) { + + case BOOLEAN: + case LONG: + case DOUBLE: + if( constant ) { + char undef=0; + for( kk=0; kkvalue.data), + &undef, result->value.nelem /* 1 */, + userInfo->datatype, Null, + (char*)Data + (kk*repeat+jj)*datasize, + &anyNullThisTime, &gParse.status ); + } else { + if ( repeat == result->value.nelem ) { + ffcvtn( gParse.datatype, + result->value.data.ptr, + result->value.undef, + result->value.nelem*ntodo, + userInfo->datatype, Null, Data, + &anyNullThisTime, &gParse.status ); + } else if( result->value.nelem == 1 ) { + for( kk=0; kkvalue.data.ptr + kk*resDataSize, + (char*)result->value.undef + kk, + 1, userInfo->datatype, Null, + (char*)Data + (kk*repeat+jj)*datasize, + &anyNullThisTime, &gParse.status ); + } + } else { + int nCopy; + nCopy = minvalue( repeat, result->value.nelem ); + for( kk=0; kkvalue.data.ptr + + kk*result->value.nelem*resDataSize, + (char*)result->value.undef + + kk*result->value.nelem, + nCopy, userInfo->datatype, Null, + (char*)Data + (kk*repeat)*datasize, + &anyNullThisTime, &gParse.status ); + if( nCopy < repeat ) { + memset( (char*)Data + (kk*repeat+nCopy)*datasize, + 0, (repeat-nCopy)*datasize); + } + } + + } + if( result->operation>0 ) { + free( result->value.data.ptr ); + } + } + if( gParse.status==OVERFLOW_ERR ) { + gParse.status = NUM_OVERFLOW; + ffpmsg("Numerical overflow while converting expression to necessary datatype"); + } + break; + + case BITSTR: + switch( userInfo->datatype ) { + case TBYTE: + idx = -1; + for( kk=0; kkvalue.nelem; jj++ ) { + if( jj%8 == 0 ) + ((char*)Data)[++idx] = 0; + if( constant ) { + if( result->value.data.str[jj]=='1' ) + ((char*)Data)[idx] |= 128>>(jj%8); + } else { + if( result->value.data.strptr[kk][jj]=='1' ) + ((char*)Data)[idx] |= 128>>(jj%8); + } + } + } + break; + case TBIT: + case TLOGICAL: + if( constant ) { + for( kk=0; kkvalue.nelem; jj++ ) { + ((char*)Data)[ jj+kk*result->value.nelem ] = + ( result->value.data.str[jj]=='1' ); + } + } else { + for( kk=0; kkvalue.nelem; jj++ ) { + ((char*)Data)[ jj+kk*result->value.nelem ] = + ( result->value.data.strptr[kk][jj]=='1' ); + } + } + break; + case TSTRING: + if( constant ) { + for( jj=0; jjvalue.data.str ); + } + } else { + for( jj=0; jjvalue.data.strptr[jj] ); + } + } + break; + default: + ffpmsg("Cannot convert bit expression to desired type."); + gParse.status = PARSE_BAD_TYPE; + break; + } + if( result->operation>0 ) { + free( result->value.data.strptr[0] ); + free( result->value.data.strptr ); + } + break; + + case STRING: + if( userInfo->datatype==TSTRING ) { + if( constant ) { + for( jj=0; jjvalue.data.str ); + } else { + for( jj=0; jjvalue.undef[jj] ) { + anyNullThisTime = 1; + strcpy( ((char**)Data)[jj], + *(char **)Null ); + } else { + strcpy( ((char**)Data)[jj], + result->value.data.strptr[jj] ); + } + } + } else { + ffpmsg("Cannot convert string expression to desired type."); + gParse.status = PARSE_BAD_TYPE; + } + if( result->operation>0 ) { + free( result->value.data.strptr[0] ); + free( result->value.data.strptr ); + } + break; + } + + if( gParse.status ) break; + + /* Increment Data to point to where the next block should go */ + + if( result->type==BITSTR && userInfo->datatype==TBYTE ) + Data = (char*)Data + + datasize * ( (result->value.nelem+7)/8 ) * ntodo; + else if( result->type==STRING ) + Data = (char*)Data + datasize * ntodo; + else + Data = (char*)Data + datasize * ntodo * repeat; + } + + /* If no NULLs encountered during this pass, set Null value to */ + /* zero to make the writing of the output column data faster */ + + if( anyNullThisTime ) + userInfo->anyNull = 1; + else if( userInfo->dataPtr == NULL ) { + if( userInfo->datatype == TSTRING ) + memcpy( *(char **)Null, zeros, 2 ); + else + memcpy( Null, zeros, datasize ); + } + + /*-------------------------------------------------------*/ + /* Clean up procedures: after processing all the rows */ + /*-------------------------------------------------------*/ + + if( firstrow - 1 == lastRow ) { + if( !gParse.status && userInfo->maxRows UCHAR_MAX ) { + *status = OVERFLOW_ERR; + ((unsigned char*)output)[i] = UCHAR_MAX; + } else + ((unsigned char*)output)[i] = + (unsigned char) ((long*)input)[i]; + } + } + return( *status ); + break; + case TFLOAT: + fffr4i1((float*)input,ntodo,1.,0.,0,0,NULL,NULL, + (unsigned char*)output,status); + break; + case TDOUBLE: + fffr8i1((double*)input,ntodo,1.,0.,0,0,NULL,NULL, + (unsigned char*)output,status); + break; + default: + *status = BAD_DATATYPE; + break; + } + for(i=0;i SHRT_MAX ) { + *status = OVERFLOW_ERR; + ((short*)output)[i] = SHRT_MAX; + } else + ((short*)output)[i] = (short) ((long*)input)[i]; + } + } + return( *status ); + break; + case TFLOAT: + fffr4i2((float*)input,ntodo,1.,0.,0,0,NULL,NULL, + (short*)output,status); + break; + case TDOUBLE: + fffr8i2((double*)input,ntodo,1.,0.,0,0,NULL,NULL, + (short*)output,status); + break; + default: + *status = BAD_DATATYPE; + break; + } + for(i=0;i=0 ) { + found[parNo] = 1; /* Flag this parameter as found */ + switch( gParse.colData[parNo].datatype ) { + case TLONG: + ffgcvj( fptr, gParse.valCol, row, 1L, 1L, + ((long*)gParse.colData[parNo].array)[0], + ((long*)gParse.colData[parNo].array)+currelem, + &anynul, status ); + break; + case TDOUBLE: + ffgcvd( fptr, gParse.valCol, row, 1L, 1L, + ((double*)gParse.colData[parNo].array)[0], + ((double*)gParse.colData[parNo].array)+currelem, + &anynul, status ); + break; + case TSTRING: + ffgcvs( fptr, gParse.valCol, row, 1L, 1L, + ((char**)gParse.colData[parNo].array)[0], + ((char**)gParse.colData[parNo].array)+currelem, + &anynul, status ); + break; + } + if( *status ) return( *status ); + } + } + + if( currelemoperation==CONST_OP ) { + + if( result->value.data.log ) { + *(long*)userPtr = firstrow; + return( -1 ); + } + + } else { + + for( idx=0; idxvalue.data.logptr[idx] && !result->value.undef[idx] ) { + *(long*)userPtr = firstrow + idx; + return( -1 ); + } + } + } + + return( gParse.status ); +} + + +/************************************************************************* + + Functions used by the evaluator to access FITS data + (find_column, find_keywd, allocateCol, load_column) + + *************************************************************************/ + + +static int find_column( char *colName, void *itslval ) +{ + FFSTYPE *thelval = (FFSTYPE*)itslval; + int col_cnt, status; + int colnum, typecode, type, hdutype; + long repeat, width; + fitsfile *fptr; + char temp[80]; + double tzero,tscale; + int istatus; + + if( *colName == '#' ) + return( find_keywd( colName + 1, itslval ) ); + + fptr = gParse.def_fptr; + status = 0; + if( gParse.compressed ) + colnum = gParse.valCol; + else + if( fits_get_colnum( fptr, CASEINSEN, colName, &colnum, &status ) ) { + if( status == COL_NOT_FOUND ) { + type = find_keywd( colName, itslval ); + if( type != pERROR ) ffcmsg(); + return( type ); + } + gParse.status = status; + return pERROR; + } + + if( fits_get_coltype( fptr, colnum, &typecode, + &repeat, &width, &status ) ) { + gParse.status = status; + return pERROR; + } + + col_cnt = gParse.nCols; + if( allocateCol( col_cnt, &gParse.status ) ) return pERROR; + fits_iter_set_by_num( gParse.colData+col_cnt, fptr, colnum, 0, InputCol ); + + /* Make sure we don't overflow variable name array */ + strncpy(gParse.varData[col_cnt].name,colName,MAXVARNAME); + gParse.varData[col_cnt].name[MAXVARNAME] = '\0'; + + switch( typecode ) { + case TBIT: + gParse.varData[col_cnt].type = BITSTR; + gParse.colData[col_cnt].datatype = TBYTE; + type = BITCOL; + break; + case TBYTE: + case TSHORT: + case TLONG: + /* The datatype of column with TZERO and TSCALE keywords might be + float or double. + */ + sprintf(temp,"TZERO%d",colnum); + istatus = 0; + if(fits_read_key(fptr,TDOUBLE,temp,&tzero,NULL,&istatus)) { + tzero = 0.0; + } + sprintf(temp,"TSCAL%d",colnum); + istatus = 0; + if(fits_read_key(fptr,TDOUBLE,temp,&tscale,NULL,&istatus)) { + tscale = 1.0; + } + if (tscale == 1.0 && (tzero == 0.0 || tzero == 32768.0 )) { + gParse.varData[col_cnt].type = LONG; + gParse.colData[col_cnt].datatype = TLONG; + } else if (tscale == 1.0 && tzero == 2147483648.0 ) { + gParse.varData[col_cnt].type = LONG; + gParse.colData[col_cnt].datatype = TULONG; + } else { + gParse.varData[col_cnt].type = DOUBLE; + gParse.colData[col_cnt].datatype = TDOUBLE; + } + type = COLUMN; + break; + case TFLOAT: + case TDOUBLE: + gParse.varData[col_cnt].type = DOUBLE; + gParse.colData[col_cnt].datatype = TDOUBLE; + type = COLUMN; + break; + case TLOGICAL: + gParse.varData[col_cnt].type = BOOLEAN; + gParse.colData[col_cnt].datatype = TLOGICAL; + type = BCOLUMN; + break; + case TSTRING: + gParse.varData[col_cnt].type = STRING; + gParse.colData[col_cnt].datatype = TSTRING; + type = SCOLUMN; + fits_get_hdu_type( fptr, &hdutype, &status ); + if( hdutype == ASCII_TBL ) repeat = width; + break; + default: + gParse.status = PARSE_BAD_TYPE; + return pERROR; + } + gParse.varData[col_cnt].nelem = repeat; + if( repeat>1 && typecode!=TSTRING ) { + if( fits_read_tdim( fptr, colnum, MAXDIMS, + &gParse.varData[col_cnt].naxis, + &gParse.varData[col_cnt].naxes[0], &status ) + ) { + gParse.status = status; + return pERROR; + } + } else { + gParse.varData[col_cnt].naxis = 1; + gParse.varData[col_cnt].naxes[0] = 1; + } + gParse.nCols++; + thelval->lng = col_cnt; + + return( type ); +} + +static int find_keywd(char *keyname, void *itslval ) +{ + FFSTYPE *thelval = (FFSTYPE*)itslval; + int status, type; + char keyvalue[FLEN_VALUE], dtype; + fitsfile *fptr; + double rval; + int bval; + long ival; + + status = 0; + fptr = gParse.def_fptr; + if( fits_read_keyword( fptr, keyname, keyvalue, NULL, &status ) ) { + if( status == KEY_NO_EXIST ) { + /* Do this since ffgkey doesn't put an error message on stack */ + sprintf(keyvalue, "ffgkey could not find keyword: %s",keyname); + ffpmsg(keyvalue); + } + gParse.status = status; + return( pERROR ); + } + + if( fits_get_keytype( keyvalue, &dtype, &status ) ) { + gParse.status = status; + return( pERROR ); + } + + switch( dtype ) { + case 'C': + fits_read_key_str( fptr, keyname, keyvalue, NULL, &status ); + type = STRING; + strcpy( thelval->str , keyvalue ); + break; + case 'L': + fits_read_key_log( fptr, keyname, &bval, NULL, &status ); + type = BOOLEAN; + thelval->log = bval; + break; + case 'I': + fits_read_key_lng( fptr, keyname, &ival, NULL, &status ); + type = LONG; + thelval->lng = ival; + break; + case 'F': + fits_read_key_dbl( fptr, keyname, &rval, NULL, &status ); + type = DOUBLE; + thelval->dbl = rval; + break; + default: + type = pERROR; + break; + } + + if( status ) { + gParse.status=status; + return pERROR; + } + + return( type ); +} + +static int allocateCol( int nCol, int *status ) +{ + if( (nCol%25)==0 ) { + if( nCol ) { + gParse.colData = (iteratorCol*) realloc( gParse.colData, + (nCol+25)*sizeof(iteratorCol) ); + gParse.varData = (DataInfo *) realloc( gParse.varData, + (nCol+25)*sizeof(DataInfo) ); + } else { + gParse.colData = (iteratorCol*) malloc( 25*sizeof(iteratorCol) ); + gParse.varData = (DataInfo *) malloc( 25*sizeof(DataInfo) ); + } + if( gParse.colData == NULL + || gParse.varData == NULL ) { + if( gParse.colData ) free(gParse.colData); + if( gParse.varData ) free(gParse.varData); + gParse.colData = NULL; + gParse.varData = NULL; + return( *status = MEMORY_ALLOCATION ); + } + } + gParse.varData[nCol].data = NULL; + gParse.varData[nCol].undef = NULL; + return 0; +} + +static int load_column( int varNum, long fRow, long nRows, + void *data, char *undef ) +{ + iteratorCol *var = gParse.colData+varNum; + long nelem,nbytes,row,len,idx; + char **bitStrs; + unsigned char *bytes; + int status = 0, anynul; + + nelem = nRows * var->repeat; + + switch( var->datatype ) { + case TBYTE: + nbytes = ((var->repeat+7)/8) * nRows; + bytes = (unsigned char *)malloc( nbytes * sizeof(char) ); + + ffgcvb(var->fptr, var->colnum, fRow, 1L, nbytes, + 0, bytes, &anynul, &status); + + nelem = var->repeat; + bitStrs = (char **)data; + for( row=0; rowfptr, var->colnum, fRow, 1L, nRows, + (char **)data, undef, &anynul, &status); + break; + case TLOGICAL: + ffgcfl(var->fptr, var->colnum, fRow, 1L, nelem, + (char *)data, undef, &anynul, &status); + break; + case TLONG: + ffgcfj(var->fptr, var->colnum, fRow, 1L, nelem, + (long *)data, undef, &anynul, &status); + break; + case TDOUBLE: + ffgcfd(var->fptr, var->colnum, fRow, 1L, nelem, + (double *)data, undef, &anynul, &status); + break; + } + + if( status ) { + gParse.status = status; + return pERROR; + } + + return 0; +} -- cgit