Changeset 25022 for branches/eam_branches/20090715/psModules/src/detrend/pmDark.c

Timestamp:

Aug 7, 2009, 12:45:24 PM (17 years ago)

Author:

eugene

Message:

update from mainline

Location:

branches/eam_branches/20090715

Files:

• . (modified) (1 prop)
• psModules/src/detrend/pmDark.c (modified) (17 diffs)

branches/eam_branches/20090715/psModules/src/detrend/pmDark.c

@@ +1 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
 #include <stdio.h>
 #include <pslib.h>
 #include <string.h>
+#include <strings.h>
 
 #include "psPolynomialMD.h"
@@ -12 +17 @@
 #include "pmReadoutStack.h"
 #include "pmDetrendThreads.h"
-
+#include "pmErrorCodes.h"
 #include "pmDark.h"
 
 #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 +28 @@
 #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 +53 @@
     }
 
-    *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 +148 @@
 {
     psFree(ord->name);
+    psFree(ord->rule);
     return;
 }
@@ -91 +158 @@
 
     ord->name = psStringCopy(name);
+    ord->rule = NULL;
     ord->order = order;
     ord->scale = false;
@@ -118 +186 @@
 
         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 +212 @@
         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 +229 @@
 
             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 +242 @@
             }
             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;
@@ -309 +387 @@
         return false;
     }
+
+    pmDarkVisualInit(values);
 
     pmReadout *outReadout = output->readouts->data[0];
@@ -352 +432 @@
                 psVectorInit(poly->coeff, NAN);
             }
+
+            pmDarkVisualPixelFit(pixels, mask);
+            pmDarkVisualPixelModel(poly, values);
+
             for (int k = 0; k < poly->coeff->n; k++) {
                 pmReadout *ro = output->readouts->data[k]; // Readout of interest
@@ -447 +531 @@
     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 +559 @@
         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 +790 @@
             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 +975 @@
               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 (!rule || !strcasecmp(rule, "NONE")) {
+                  ord->rule = NULL;
+              } else {
+                  ord->rule = psStringCopy(rule);
+              }
               ordinates->data[i] = ord;
           }
@@ -892 +997 @@
 }
 
+#if (HAVE_KAPA)
+#include <kapa.h>
+#include "pmKapaPlots.h"
+#include "pmVisual.h"
+
+static int nKapa = 0;
+static int *kapa = NULL;
+static bool plotFlag = true;
+static psArray *xVectors = NULL;
+
+// this init function only gets the ordinates for the first readout...
+bool pmDarkVisualInit(psArray *values) {
+    
+    if (!pmVisualIsVisual()) return true;
+
+    // skip if we have already opened the windows (or if none are requested...)
+    if (nKapa) return true;
+
+    // values has Ninput vectors with Norder elements; we need Norder vectors with Ninput elements...
+    int nOrders = 0;
+    for (int i = 0; i < values->n; i++) {
+        psVector *vect = values->data[i];
+        if (!nOrders) {
+            nOrders = vect->n;
+        } else {
+            psAssert (nOrders == vect->n, "mismatch in order vector lengths");
+        }
+    }
+    xVectors = psArrayAlloc(nOrders);
+    for (int i = 0; i < nOrders; i++) {
+        xVectors->data[i] = psVectorAlloc(values->n, PS_TYPE_F32);
+    }
+
+    for (int i = 0; i < values->n; i++) {
+        psVector *vect = values->data[i];
+        for (int j = 0; j < vect->n; j++) {
+            psVector *xVec = xVectors->data[j];
+            xVec->data.F32[i] = vect->data.F32[j];
+        }
+    }
+
+    nKapa = nOrders;
+
+    kapa = psAlloc(nKapa*sizeof(int));
+    
+    for (int i = 0; i < nKapa; i++) {
+        kapa[i] = -1;
+        pmVisualInitWindow(&kapa[i], "ppmerge");
+    }
+    return true;
+}
+
+bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask) {
+
+    Graphdata graphdata;
+
+    if (!pmVisualIsVisual()) return true;
+
+    KapaInitGraph(&graphdata);
+
+    psAssert(nKapa == xVectors->n, "inconsistent number of orders %d vs %ld\n", nKapa, xVectors->n);
+
+    psVector *xSub = psVectorAlloc(pixels->n, PS_TYPE_F32);
+    psVector *ySub = psVectorAlloc(pixels->n, PS_TYPE_F32);
+
+    for (int i = 0; i < xVectors->n; i++) {
+        psVector *x = xVectors->data[i];
+
+        // generate vectors of the unmasked values
+        int nSub = 0;
+        for (int j = 0; j < pixels->n; j++) {
+            if (mask && mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) continue;
+            xSub->data.F32[nSub] = x->data.F32[j];
+            ySub->data.F32[nSub] = pixels->data.F32[j];
+            nSub ++;
+        }
+        xSub->n = ySub->n = nSub;
+        
+        // plot the unmasked values
+        pmVisualScaleGraphdata (&graphdata, xSub, ySub, false);
+        KapaSetGraphData(kapa[i], &graphdata);
+        KapaSetLimits(kapa[i], &graphdata);
+        KapaClearPlots (kapa[i]);
+
+        KapaSetFont (kapa[i], "courier", 14);
+        KapaBox (kapa[i], &graphdata);
+        KapaSendLabel (kapa[i], "ordinate", KAPA_LABEL_XM);
+        KapaSendLabel (kapa[i], "pixel values", KAPA_LABEL_YM);
+
+        graphdata.color = KapaColorByName("black");
+        graphdata.style = 2;
+        graphdata.ptype = 2;
+        KapaPrepPlot  (kapa[i], xSub->n, &graphdata);
+        KapaPlotVector(kapa[i], xSub->n, xSub->data.F32, "x");
+        KapaPlotVector(kapa[i], xSub->n, ySub->data.F32, "y");
+    }
+    pmVisualAskUser (&plotFlag);
+    return true;
+}
+
+bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values) {
+
+    Graphdata graphdata;
+
+    if (!pmVisualIsVisual()) return true;
+
+    KapaInitGraph(&graphdata);
+
+    psAssert(nKapa == xVectors->n, "inconsistent number of orders %d vs %ld\n", nKapa, xVectors->n);
+
+    psVector *yFit = psVectorAlloc(values->n, PS_TYPE_F32);
+
+    for (int i = 0; i < values->n; i++) {
+        psVector *coord = values->data[i];
+        yFit->data.F32[i] = psPolynomialMDEval (poly, coord);
+    }
+
+    for (int i = 0; i < xVectors->n; i++) {
+        psVector *xFit = xVectors->data[i];
+
+        KapaGetGraphData(kapa[i], &graphdata);
+        graphdata.color = KapaColorByName("red");
+        graphdata.style = 2;
+        graphdata.ptype = 7;
+        KapaPrepPlot  (kapa[i], xFit->n, &graphdata);
+        KapaPlotVector(kapa[i], xFit->n, xFit->data.F32, "x");
+        KapaPlotVector(kapa[i], xFit->n, yFit->data.F32, "y");
+    }
+    pmVisualAskUser (&plotFlag);
+    return true;
+}
+
+bool pmDarkVisualCleanup() {
+
+    for (int i = 0; i < nKapa; i++) {
+        KapaClose(kapa[i]);
+    }
+    psFree (kapa);
+    psFree (xVectors);
+    return true;
+}
+
+# else
+
+bool pmDarkVisualInit(psArray *values) { return true; }
+bool pmDarkVisualPixelFit(psVector *pixels, psVector *mask) { return true; }
+bool pmDarkVisualPixelModel(psPolynomialMD *poly, psArray *values) { return true; }
+bool pmDarkVisualCleanup() { return true; }
+
+# endif