Index: trunk/psModules/src/detrend/pmDark.c =================================================================== --- trunk/psModules/src/detrend/pmDark.c (revision 24802) +++ trunk/psModules/src/detrend/pmDark.c (revision 24803) @@ -2,4 +2,5 @@ #include #include +#include #include "psPolynomialMD.h" @@ -12,4 +13,5 @@ #include "pmReadoutStack.h" #include "pmDetrendThreads.h" +#include "pmErrorCodes.h" #include "pmDark.h" @@ -17,4 +19,5 @@ #define PM_DARK_FITS_EXTNAME "PS_DARK" // FITS extension name for ordinates table #define PM_DARK_FITS_NAME "NAME" // Column name for concept name in ordinates table +#define PM_DARK_FITS_RULE "RULE" // Column name for concept rule in ordinates table #define PM_DARK_FITS_ORDER "ORDER" // Column name for polynomial order in ordinates table #define PM_DARK_FITS_SCALE "SCALE" // Column name for scaling option in ordinates table @@ -22,40 +25,24 @@ #define PM_DARK_FITS_MAX "MAX" // Column name for maximum value in ordinates table - - -// Look up the value of an ordinate in a readout -static bool ordinateLookup(float *value, // Value of ordinate, to return - bool *inRange, // is value within min : max range? - const char *name, // Name of ordinate (concept name) - bool scale, // Scale the value? - float min, float max, // Minimum and maximum values for scaling - const pmReadout *ro // Readout of interest - ) -{ - assert(value); - assert(name); - assert(ro); - - *inRange = true; - - pmCell *cell = ro->parent; // Parent cell - if (!cell) { - return false; - } +static bool ordinateParseConcept(double *value, const pmReadout *readout, const char *name) { + + *value = NAN; + + pmCell *cell = readout->parent; // Parent cell + psAssert(cell, "readout is missing cell \n"); + psMetadataItem *item = psMetadataLookup(cell->concepts, name); if (!item) { pmChip *chip = cell->parent; // Parent chip - if (!chip) { - return false; - } + psAssert(chip, "cell is missing chip \n"); + item = psMetadataLookup(chip->concepts, name); if (!item) { pmFPA *fpa = chip->parent; // Parent FPA - if (!fpa) { - return false; - } + psAssert(fpa, "chip is missing fpa \n"); + item = psMetadataLookup(fpa->concepts, name); if (!item) { - psWarning("Unable to find concept %s in readout", name); + psError(PM_ERR_CONFIG, true, "Unable to find concept %s in readout", name); return false; } @@ -63,9 +50,85 @@ } - *value = psMetadataItemParseF32(item); // Value of interest + *value = psMetadataItemParseF64(item); // Value of interest if (!isfinite(*value)) { psWarning("Non-finite value (%f) of concept %s in readout", *value, name); - return false; - } + } + return true; +} + +static bool ordinateParseRule(double *value, const pmReadout *readout, const char *name, const char *rule) { + + psAssert(name, "ordinate name is not defined"); + psAssert(rule, "ordinate rule is not defined"); + + *value = NAN; + + psArray *words = psStringSplitArray(rule, " ", false); + + // we should have a rule of the form (concept) OP (concept) OP (concept) ... + // for now, the only allowed OP is * (eventually, we can steal code from opihi for a better + // RPN parser). + + if (words->n % 2 == 0) { + psError(PM_ERR_CONFIG, true, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + + for (int i = 1; i < words->n; i+=2) { + if (strcmp((char *)words->data[i], "*")) { + psError(PM_ERR_CONFIG, true, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + } + + if (!ordinateParseConcept(value, readout, words->data[0])) { + psError(PM_ERR_CONFIG, false, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + + double value2 = 0.0; + for (int i = 2; i < words->n; i+=2) { + if (!ordinateParseConcept(&value2, readout, words->data[i])) { + psError(PM_ERR_CONFIG, false, "syntax error in DARK.ORDINATE %s rule %s\n", name, rule); + psFree(words); + return false; + } + *value *= value2; + } + psFree(words); + return true; +} + +// Look up the value of an ordinate in a readout +static bool ordinateLookup(double *value, // Value of ordinate, to return + bool *inRange, // is value within min : max range? + const char *name, // Name of ordinate (concept or abstract name) + const char *rule, // Rule for generating the value (if NULL use name as concept) + bool scale, // Scale the value? + float min, float max, // Minimum and maximum values for scaling + const pmReadout *readout // Readout of interest + ) +{ + assert(value); + assert(name); + assert(readout); + + *inRange = true; + + if (rule) { + if (!ordinateParseRule(value, readout, name, rule)) { + psError(PM_ERR_CONFIG, false, "trouble parsing rule %s for DARK.ORDINATE %s", rule, name); + return false; + } + } else { + if (!ordinateParseConcept(value, readout, name)) { + psError(PM_ERR_CONFIG, false, "trouble parsing rule %s for DARK.ORDINATE %s", rule, name); + return false; + } + } + if (scale) { if (*value < min || *value > max) { @@ -82,4 +145,5 @@ { psFree(ord->name); + psFree(ord->rule); return; } @@ -91,4 +155,5 @@ ord->name = psStringCopy(name); + ord->rule = NULL; ord->order = order; ord->scale = false; @@ -118,7 +183,11 @@ pmReadout *readout = inputs->data[i]; // Readout of interest - float normValue; // Normalisation value - if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) { - psWarning("Unable to find value of %s for readout %d", normConcept, i); + double normValue; // Normalisation value + if (!ordinateLookup(&normValue, &inRange, normConcept, NULL, false, NAN, NAN, readout)) { + psError(PM_ERR_CONFIG, false, "problem finding concept %s for DARK.NORM", normConcept); + return false; + } + if (!isfinite(normValue)) { + psWarning("Unable to find acceptable value of %s for readout %d", normConcept, i); roMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; norm->data.F32[i] = NAN; @@ -140,5 +209,5 @@ pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information if (ord->order <= 0) { - psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order); + psError(PS_ERR_UNKNOWN, true, "Bad order for DARK.ORDINATE %s (%d) --- ignored", ord->name, ord->order); psFree(values); psFree(roMask); @@ -157,6 +226,11 @@ pmReadout *readout = inputs->data[j]; // Readout of interest - float value = NAN; // Value of ordinate - if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) { + double value = NAN; // Value of ordinate + if (!ordinateLookup(&value, &inRange, ord->name, ord->rule, ord->scale, ord->min, ord->max, readout)) { + psError(PM_ERR_CONFIG, false, "problem finding rule for DARK.ORDINATE %s", ord->name); + return false; + } + if (!isfinite(value)) { + psWarning("Unable to find acceptable value of DARK.ORDINATE %s for readout %d", ord->name, i); roMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; val->data.F32[i] = NAN; @@ -165,4 +239,5 @@ } if (!inRange) { + psWarning("Value of DARK.ORDINATE %s for readout %d is out of range", ord->name, i); roMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; val->data.F32[i] = NAN; @@ -447,16 +522,21 @@ for (int i = 0; i < numOrdinates; i++) { pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest - float value = NAN; // Value for ordinate - if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) { - psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", ord->name); + double value = NAN; // Value for ordinate + if (!ordinateLookup(&value, &inRange, ord->name, ord->rule, ord->scale, ord->min, ord->max, readout)) { + psError(PS_ERR_UNKNOWN, true, "Unable to find value for DARK.ORDINATE %s", ord->name); psFree(values); return false; } + if (!isfinite(value)) { + psError(PS_ERR_UNKNOWN, true, "Value for DARK.ORDINATE %s is NAN", ord->name); + psFree(values); + return false; + } values->data.F32[i] = value; } - float norm = NAN; // Normalisation value + double norm = NAN; // Normalisation value bool doNorm = false; // Do normalisation? if (normConcept && strlen(normConcept) > 0) { - if (!ordinateLookup(&norm, &inRange, normConcept, false, NAN, NAN, readout)) { + if (!ordinateLookup(&norm, &inRange, normConcept, NULL, false, NAN, NAN, readout)) { psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept); psFree(values); @@ -470,5 +550,5 @@ pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information if (ord->order <= 0) { - psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order); + psError(PS_ERR_UNKNOWN, true, "Bad order for DARK.ORDINATE %s (%d) --- ignored", ord->name, ord->order); psFree(values); psFree(orders); @@ -701,5 +781,13 @@ pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest psMetadata *row = psMetadataAlloc(); // FITS table row - psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name); + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "DARK.ORDINATE name", ord->name); + + // XXX write a dummy value if ord->rule == NULL? (eg, NONE) + if (ord->rule) { + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_RULE, 0, "DARK.ORDINATE rule", ord->rule); + } else { + psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_RULE, 0, "DARK.ORDINATE rule", "NONE"); + } + psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order); psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale); @@ -878,4 +966,12 @@ ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX); + // load the ordinate rule; it is not an error if this field is missing or NULL + // a NULL rule means 'use the name as the single concept' + const char *rule = psMetadataLookupStr(&mdok, row, PM_DARK_FITS_RULE); + if (!strcasecmp(rule, "NONE")) { + ord->rule = NULL; + } else { + ord->rule = psStringCopy(rule); + } ordinates->data[i] = ord; }