Changeset 24951 for branches/czw_branch/cleanup/psModules/src/detrend/pmDark.c

Timestamp:

Jul 30, 2009, 5:20:29 PM (17 years ago)

Author:

watersc1

Message:

Finished up my edits to the detrend cleanup, and some changes to my
copy of burntool and the pslib astrometry. Detrend cleanup has not
been tested yet. That's up next.

Location:

branches/czw_branch/cleanup

Files:

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

branches/czw_branch/cleanup/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