Changeset 30266 for branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Jan 12, 2011, 9:32:08 PM (16 years ago)

Author:

eugene

Message:

add errors to kernel coeffs; calculate and plot chisq and moments

File:

• branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (modified) (18 diffs)

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c

@@ -8 +8 @@
 
 #include "pmErrorCodes.h"
+#include "pmVisual.h"
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
@@ -16 +17 @@
 #include "pmSubtractionVisual.h"
 
-//#define TESTING                         // TESTING output for debugging; may not work with threads!
-
-//#define USE_WEIGHT                      // Include weight (1/variance) in equation?
-
-// XXX TEST:
-//# define USE_WINDOW                      // window to avoid neighbor contamination
+//# define TESTING                         // TESTING output for debugging; may not work with threads!
+//# define USE_WEIGHT                      // Include weight (1/variance) in equation?
+# define USE_WINDOW                      // window to avoid neighbor contamination
+// XXX if we want to have a weight and window, we'll need to pass through to the chisq function
 
 # define PENALTY false
@@ -983 +982 @@
     }
 
-    pmSubtractionVisualPlotLeastSquares(stamps);
     pmSubtractionVisualShowKernels((pmSubtractionKernels  *)kernels);
     pmSubtractionVisualShowBasis(stamps);
+    pmSubtractionVisualPlotLeastSquares(stamps);
 
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
@@ -1080 +1079 @@
         }
 
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
+
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
@@ -1087 +1088 @@
 #endif
 
+        psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
         solution = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -1096 +1098 @@
         if (1) {
             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+            invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
         } else {
             psVector *PERM = NULL;
@@ -1104 +1107 @@
         }
 
-# if (0)
+# if (1)
         for (int i = 0; i < solution->n; i++) {
-            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
         }
 # endif
@@ -1113 +1116 @@
             kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
             psVectorInit(kernels->solution1, 0.0);
+            kernels->solution1err = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+            psVectorInit(kernels->solution1err, 0.0);
         }
 
@@ -1120 +1125 @@
         for (int i = 0; i < numKernels * numPoly; i++) {
             kernels->solution1->data.F64[i] = solution->data.F64[i];
+            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]);
         }
 
@@ -1131 +1137 @@
         psFree(sumVector);
         psFree(sumMatrix);
+        psFree(invMatrix);
 
     } else {
@@ -1160 +1167 @@
         }
 
+        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
+
 #if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
@@ -1171 +1180 @@
         }
 
+        psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
         solution = psVectorAlloc(numParams, PS_TYPE_F64);
@@ -1176 +1186 @@
 
         // DUAL solution
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-
-#if (0)
+# if (1)
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7);
+        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
+# endif
+# if (0)
+        psMatrixLUSolve(sumMatrix, sumVector);
+        solution = sumVector;
+        invMatrix = sumMatrix;
+# endif
+
+#if (1)
         for (int i = 0; i < solution->n; i++) {
-            fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
+            fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
         }
 #endif
@@ -1187 +1205 @@
             kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64);
             psVectorInit (kernels->solution1, 0.0);
+            kernels->solution1err = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+            psVectorInit(kernels->solution1err, 0.0);
         }
         if (!kernels->solution2) {
             kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
             psVectorInit (kernels->solution2, 0.0);
+            kernels->solution2err = psVectorAlloc(numSolution2, PS_TYPE_F64);
+            psVectorInit(kernels->solution2err, 0.0);
         }
 
@@ -1205 +1227 @@
             kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue;
             kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
+
+            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]) * stamps->normValue;
+            int i2 = i + numSolution1;
+            kernels->solution2err->data.F64[i] = sqrt(invMatrix->data.F64[i2][i2]);
         }
 
@@ -1213 +1239 @@
         kernels->solution1->data.F64[bgIndex] = 0.0;
 
+        psFree(solution);
+        psFree(sumVector);
         psFree(sumMatrix);
-        psFree(sumVector);
-        psFree(solution);
+        psFree(invMatrix);
     }
 
@@ -1224 +1251 @@
     if (psTraceGetLevel("psModules.imcombine") >= 7) {
         for (int i = 0; i < kernels->solution1->n; i++) {
-            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]);
+            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f +/- %f\n", i, kernels->solution1->data.F64[i], kernels->solution1err->data.F64[i]);
         }
         if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
             for (int i = 0; i < kernels->solution2->n; i++) {
-                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]);
+                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f +/- %f\n", i, kernels->solution2->data.F64[i], kernels->solution2err->data.F64[i]);
             }
         }
@@ -1280 +1307 @@
     psVectorAppend(fResOuter, dOuter/sum);
     psVectorAppend(fResTotal, dTotal/sum);
+    return true;
+}
+
+// given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) {
+
+    int npix = 0;
+    float chisq = 0;
+
+    // get the chisq
+    for (int y = residual->yMin; y <= residual->yMax; y++) {
+        for (int x = residual->xMin; x <= residual->xMax; x++) {
+            chisq  += PS_SQR(residual->kernel[y][x]);
+            npix ++;
+        }
+    }
+    psVectorAppend(chisqVector, chisq / npix);
+
+    float value1 = 0;
+    float value2 = 0;
+    float flux2 = 0;
+    float fluxX = 0;
+    float fluxY = 0;
+    float fluxX2 = 0;
+    float fluxY2 = 0;
+
+    float fluxC1 = 0;
+    float fluxC2 = 0;
+
+    float moment = 0;
+
+    // get the moments from convolved1
+    if (convolved1) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved1->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                if (weight) {
+                    value2 *= weight->kernel[y][x];
+                }
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC1 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                fluxX2 += PS_SQR(x) * value2;
+                fluxY2 += PS_SQR(y) * value2;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        float Mxx = fluxX2 / flux2;
+        float Myy = fluxY2 / flux2;
+
+        // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    // get the moments from convolved1
+    if (convolved2) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved2->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                if (weight) {
+                    value2 *= weight->kernel[y][x];
+                }
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC2 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                fluxX2 += PS_SQR(x) * value2;
+                fluxY2 += PS_SQR(y) * value2;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        float Mxx = fluxX2 / flux2;
+        float Myy = fluxY2 / flux2;
+
+        // fprintf (stderr, "conv2, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    float flux = fluxC1 + fluxC2;
+    
+    if (convolved1 && convolved2) {
+        moment *= 0.5;
+        flux *= 0.5;
+    }
+    psVectorAppend(momentVector, moment);
+    psVectorAppend(fluxesVector, flux);
+    return true;
+}
+
+// typedef struct {
+//     float moment;
+//     float chisq;
+// } diffStat;
+
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps,
+                                           pmSubtractionKernels *kernels)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
+
+    psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+
+    int footprint = stamps->footprint; // Half-size of stamps
+    int numKernels = kernels->num;      // Number of kernels
+
+    psImage *polyValues = NULL;         // Polynomial values
+
+    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    // storage for the convolved source image (only convolved1 is used in SINGLE mode)
+    psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            continue;
+        }
+
+        // Calculate coefficients of the kernel basis functions
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+
+        // Calculate residuals
+        psImageInit(residual->image, 0.0);
+
+    psKernel *weight = NULL;
+    psKernel *window = NULL;
+
+#ifdef USE_WEIGHT
+    weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+    window = stamps->window;
+#endif
+
+        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
+
+            // the single-direction psf match code attempts to find the kernel such that:
+            // source * kernel = target.  we need to assign 'source' and 'target' correctly
+            // depending on which of image1 or image2 we asked to be convolved.
+
+            psKernel *target;           // Target postage stamp (convolve source to match the target)
+            psKernel *source;           // Source postage stamp (convolve source to match the target)
+            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
+
+            // init the accumulation image
+            psImageInit(convolved1->image, 0.0);
+
+            switch (kernels->mode) {
+              case PM_SUBTRACTION_MODE_1:
+                target = stamp->image2;
+                source = stamp->image1;
+                convolutions = stamp->convolutions1;
+                break;
+              case PM_SUBTRACTION_MODE_2:
+                target = stamp->image1;
+                source = stamp->image2;
+                convolutions = stamp->convolutions2;
+                break;
+              default:
+                psAbort("Unsupported subtraction mode: %x", kernels->mode);
+            }
+
+            // generate the convolved source image (sum over kernels)
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *convolution = convolutions->data[j]; // Convolution
+                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        convolved1->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+                    }
+                }
+            }
+
+            // generate the residual image
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    convolved1->kernel[y][x] += source->kernel[y][x] * norm;
+                    residual->kernel[y][x] = target->kernel[y][x] - convolved1->kernel[y][x] - background;
+                }
+            }
+            // XXX if we want to have a weight and window, we'll need to pass through to here
+            pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, NULL, residual, weight, window);
+
+        } else {
+
+            // Dual convolution
+            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
+            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
+            psKernel *image1 = stamp->image1; // The first image
+            psKernel *image2 = stamp->image2; // The second image
+
+            // init the accumulation images
+            psImageInit(convolved1->image, 0.0);
+            psImageInit(convolved2->image, 0.0);
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
+                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
+                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
+                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
+
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        // NOTE sign for coeff2
+                        convolved1->kernel[y][x] += +conv1->kernel[y][x] * coeff1;
+                        convolved2->kernel[y][x] += -conv2->kernel[y][x] * coeff2;
+                    }
+                }
+            }
+
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    convolved1->kernel[y][x] += image1->kernel[y][x] * norm;
+                    convolved2->kernel[y][x] += image2->kernel[y][x];
+                    residual->kernel[y][x] = convolved2->kernel[y][x] - convolved1->kernel[y][x] - background;
+                }
+            }
+
+            if (0) {
+                psFitsWriteImageSimple("conv1.fits", convolved1->image, NULL);
+                psFitsWriteImageSimple("conv2.fits", convolved2->image, NULL);
+                psFitsWriteImageSimple("resid.fits", residual->image,   NULL);
+                pmVisualAskUser(NULL);
+            } 
+
+            pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, convolved2, residual, weight, window);
+        }
+    }
+
+    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments);
+
+    // find the mean chisq and mean moment
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+    psVectorStats (stats, chisq, NULL, NULL, 0);
+    float chisqValue = stats->sampleMean;
+
+    psStatsInit(stats);
+    psVectorStats (stats, moments, NULL, NULL, 0);
+    float momentValue = stats->sampleMean;
+
+    fprintf (stderr, "chisq: %f, moment: %f\n", chisqValue, momentValue);
+
+    psFree(stats);
+    psFree(chisq);
+    psFree(fluxes);
+    psFree(moments);
+
+    psFree(residual);
+    psFree(convolved1);
+    psFree(convolved2);
+    psFree(polyValues);
+
     return true;
 }