IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2011, 12:19:53 PM (15 years ago)
Author:
eugene
Message:

some code reorganzation (move all types into pmSubtractionTypes.h); add stage to check on different mode and order options, choosing best based on chisq of subtraction; careful with the window and the kernel sizes: if the measured kron radius is too large, allow the window to grow as needed; use the 1st radial moment measured on the stamps to set the kernel scaling; remove some cruft from the code; calling program now passes in kernel scaling options to be used when 1st radial moment is measured; new function to update the kernel description string used for kernel I/O; update the kernel description after the spatial order and direction is chosen; need to carry the kernel fwhms and spatial order to allow for update; calculate the psf-match solution errors; psf solution now uses the weights to get sensible chisq values, window is used to calculate the moments; penalty scale is now adjusted to be consistent with sensible reduced chisq; improved visualizations; use range-limiter in SVD of 1e-10; calculate chisq and moments for the solution and use in the evaluation; split out the stamp selection / convolution steps; reject sources after fitting a model to the chisq assuming non-zero systematic error

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionEquation.c

    r29598 r30622  
    88
    99#include "pmErrorCodes.h"
     10#include "pmVisual.h"
     11#include "pmFPA.h"
     12#include "pmSubtractionTypes.h"
    1013#include "pmSubtraction.h"
    1114#include "pmSubtractionKernels.h"
    … …  
    1619#include "pmSubtractionVisual.h"
    1720
    18 //#define TESTING                         // TESTING output for debugging; may not work with threads!
    19 
    20 //#define USE_WEIGHT                      // Include weight (1/variance) in equation?
    21 
    22 // XXX TEST:
    23 //# define USE_WINDOW                      // window to avoid neighbor contamination
     21//# define TESTING                         // TESTING output for debugging; may not work with threads!
     22# define USE_WEIGHT                      // Include weight (1/variance) in equation?
     23# define USE_WINDOW                      // window to avoid neighbor contamination
     24
     25/* I believe we want to apply the WEIGHT to the chisq portions of the calculation (but not the WINDOW),
     26 * and the WINDOW to the moments portiosn of the calculations (but not the WEIGHT)
     27 *
     28 */
    2429
    2530# define PENALTY false
    2631# define MOMENTS (!PENALTY)
    27 # define MOMENTS_PENALTY_SCALE 2e-5 // XXX this value is not completely arbitrary, but I don't understand why it needs to be this value...
     32# define MOMENTS_PENALTY_SCALE 20 // up-weight the moments somewhat
    2833
    2934//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    107112                        cc *= weight->kernel[y][x];
    108113                    }
    109                     if (window) {
     114                    // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
     115                    if (false && window) {
    110116                        cc *= window->kernel[y][x];
    111117                    }
    … …  
    138144                    rc *= wtVal;
    139145                }
    140                 if (window) {
     146                // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
     147                if (false && window) {
    141148                    float winVal = window->kernel[y][x];
    142149                    ic *= winVal;
    … …  
    173180}
    174181
     182# define RENORM_BY_FLUX 0
    175183
    176184// Calculate the least-squares matrix and vector for dual convolution
    … …  
    268276            for (int y = - footprint; y <= footprint; y++) {
    269277                for (int x = - footprint; x <= footprint; x++) {
     278
     279                    // XXX NOTE: clipping low S/N pixels does not seem to work very well
     280                    if (false && weight) {
     281                        float i1 = image1->kernel[y][x];
     282                        float i2 = image2->kernel[y][x];
     283                        float sn = (i1 + i2) / sqrt (weight->kernel[y][x]);
     284                        if (sn < 0.5) continue;
     285                    }
     286
    270287                    double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue);
    271288                    double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
    272289                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
    273                     if (weight) {
    274                         float wtVal = weight->kernel[y][x];
    275                         aa *= wtVal;
    276                         bb *= wtVal;
    277                         ab *= wtVal;
    278                     }
    279                     if (window) {
    280                         float wtVal = window->kernel[y][x];
    281                         aa *= wtVal;
    282                         bb *= wtVal;
    283                         ab *= wtVal;
    284                     }
    285                     sumAA += aa;
    286                     sumBB += bb;
    287                     sumAB += ab;
     290
     291                    float wtVal = (weight) ? weight->kernel[y][x] : 1.0;
     292                    sumAA += wtVal*aa;
     293                    sumBB += wtVal*bb;
     294                    sumAB += wtVal*ab;
    288295
    289296                    if (MOMENTS) {
    290                         MxxAA += x*x*aa;
    291                         MyyAA += y*y*aa;
    292                         MxxBB += x*x*bb;
    293                         MyyBB += y*y*bb;
     297                        float winVal = (window) ? window->kernel[y][x] : 1.0;
     298                        MxxAA += winVal*x*x*aa;
     299                        MyyAA += winVal*y*y*aa;
     300                        MxxBB += winVal*x*x*bb;
     301                        MyyBB += winVal*y*y*bb;
    294302                    }
    295303                }
    296304            }
    297305
     306            // XXX does normSquare1,2 mess up the relative scaling?
     307            // XXX no: normSquare1,2 is the sum of the flux^2 for the source
    298308            if (MOMENTS) {
    299309                MxxAA /= stamp->normSquare1 * PS_SQR(normValue);
    … …  
    305315            // XXX this makes the Chisq portion independent of the normalization and star flux
    306316            // but may be mis-scaling between stars of different fluxes
     317# if (RENORM_BY_FLUX)       
    307318            sumAA /= PS_SQR(stamp->normI2);
    308319            sumAB /= PS_SQR(stamp->normI2);
    309320            sumBB /= PS_SQR(stamp->normI2);
     321# endif
    310322
    311323            // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB);
    … …  
    328340                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
    329341                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
    330 
    331342                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
    332343                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
    … …  
    334345                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
    335346                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
    336 
    337347                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
    338348                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
    … …  
    354364                        ab *= weight->kernel[y][x];
    355365                    }
    356                     if (window) {
     366                    // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
     367                    if (false && window) {
    357368                        ab *= window->kernel[y][x];
    358369                    }
    … …  
    363374            // XXX this makes the Chisq portion independent of the normalization and star flux
    364375            // but may be mis-scaling between stars of different fluxes
     376# if (RENORM_BY_FLUX)
    365377            sumAB /= PS_SQR(stamp->normI2);
     378# endif
    366379
    367380            // Spatial variation of kernel coeffs
    … …  
    391404                float i2 = image2->kernel[y][x];
    392405
     406                // XXX NOTE: clipping low S/N pixels does not seem to work very well
     407                if (false && weight) {
     408                    float sn = (i1 + i2) / sqrt (weight->kernel[y][x]);
     409                    if (sn < 0.5) continue;
     410                }
     411
    393412                double ai2 = a * i2 * normValue;
    394413                double bi2 = b * i2;
    … …  
    396415                double bi1 = b * i1 * normValue;
    397416
    398                 if (weight) {
    399                     float wtVal = weight->kernel[y][x];
    400                     ai2 *= wtVal;
    401                     bi2 *= wtVal;
    402                     ai1 *= wtVal;
    403                     bi1 *= wtVal;
    404                 }
    405                 if (window) {
    406                     float wtVal = window->kernel[y][x];
    407                     ai2 *= wtVal;
    408                     bi2 *= wtVal;
    409                     ai1 *= wtVal;
    410                     bi1 *= wtVal;
    411                 }
    412                 sumAI2 += ai2;
    413                 sumBI2 += bi2;
    414                 sumAI1 += ai1;
    415                 sumBI1 += bi1;
     417                float wtVal = (weight) ? weight->kernel[y][x] : 1.0;
     418                sumAI2 += wtVal*ai2;
     419                sumBI2 += wtVal*bi2;
     420                sumAI1 += wtVal*ai1;
     421                sumBI1 += wtVal*bi1;
    416422
    417423                if (MOMENTS) {
    418                     MxxAI1 += x*x*ai1;
    419                     MyyAI1 += y*y*ai1;
    420                     MxxBI2 += x*x*bi2;
    421                     MyyBI2 += y*y*bi2;
     424                    float winVal = (window) ? window->kernel[y][x] : 1.0;
     425                    MxxAI1 += winVal*x*x*ai1;
     426                    MyyAI1 += winVal*y*y*ai1;
     427                    MxxBI2 += winVal*x*x*bi2;
     428                    MyyBI2 += winVal*y*y*bi2;
    422429                }
    423430            }
    … …  
    435442        // XXX this makes the Chisq portion independent of the normalization and star flux
    436443        // but may be mis-scaling between stars of different fluxes
     444# if (RENORM_BY_FLUX)
    437445        sumAI1 /= PS_SQR(stamp->normI2);
    438446        sumBI1 /= PS_SQR(stamp->normI2);
    439447        sumAI2 /= PS_SQR(stamp->normI2);
    440448        sumBI2 /= PS_SQR(stamp->normI2);
     449# endif
    441450
    442451        // Spatial variation
    … …  
    788797    stamp->normSquare2 = normSquare2;
    789798
    790     psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
     799    // psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
    791800
    792801    return true;
    … …  
    800809    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
    801810    int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
    802     pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
    803 
    804     return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode);
    805 }
    806 
    807 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
    808                                          int index, const pmSubtractionEquationCalculationMode mode)
     811
     812    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
     813}
     814
     815bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, int index)
    809816{
    810817    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    833840    psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
    834841
    835     // Generate convolutions: these are generated once and saved
    836     if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
    837         psError(psErrorCodeLast(), false, "Unable to convolve stamp %d.", index);
    838         return NULL;
    839     }
    840 
    841 #ifdef TESTING
    842     for (int j = 0; j < numKernels; j++) {
    843         if (stamp->convolutions1) {
    844             psString convName = NULL;
    845             psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
    846             psFits *fits = psFitsOpen(convName, "w");
    847             psFree(convName);
    848             psKernel *conv = stamp->convolutions1->data[j];
    849             psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
    850             psFitsClose(fits);
    851         }
    852 
    853         if (stamp->convolutions2) {
    854             psString convName = NULL;
    855             psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
    856             psFits *fits = psFitsOpen(convName, "w");
    857             psFree(convName);
    858             psKernel *conv = stamp->convolutions2->data[j];
    859             psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
    860             psFitsClose(fits);
    861         }
    862     }
    863 #endif
    864 
    865     // XXX visualize the set of convolved stamps
    866 
    867     psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
    868                                                     stamp->xNorm, stamp->yNorm); // Polynomial terms
    869 
    870     bool new = stamp->vector ? false : true; // Is this a new run?
     842    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
     843
     844    // Is this a new run? Have we allocated the correct sized vector/matrix?
     845    bool new = stamp->vector ? false : true;
     846    if (!new && (stamp->vector->n != numParams)) {
     847        psFree (stamp->vector);
     848        psFree (stamp->matrix);
     849        new = true;
     850    }
     851
    871852    if (new) {
    872853        stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    … …  
    941922}
    942923
    943 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
    944                                     const pmSubtractionEquationCalculationMode mode)
     924bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
    945925{
    946926    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    953933    for (int i = 0; i < stamps->num; i++) {
    954934        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     935
    955936        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
    956937            continue;
    … …  
    969950            psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
    970951            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    971             PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
    972952            if (!psThreadJobAddPending(job)) {
    973953                return false;
    974954            }
    975955        } else {
    976             pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);
     956            pmSubtractionCalculateEquationStamp(stamps, kernels, i);
    977957        }
    978958    }
    … …  
    983963    }
    984964
    985     pmSubtractionVisualPlotLeastSquares(stamps);
    986965    pmSubtractionVisualShowKernels((pmSubtractionKernels  *)kernels);
    987966    pmSubtractionVisualShowBasis(stamps);
    988 
    989     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
    990              psTimerClear("pmSubtractionCalculateEquation"));
    991 
     967    pmSubtractionVisualPlotLeastSquares(stamps);
     968
     969    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", psTimerClear("pmSubtractionCalculateEquation"));
    992970
    993971    return true;
    … …  
    998976bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header);
    999977
    1000 psImage *p_pmSubSolve_wUt (psVector *w, psImage *U);
    1001 psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt);
    1002 
    1003 bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask);
    1004 
    1005 bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B);
    1006 bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB);
    1007 bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB);
    1008 
    1009 bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x);
    1010 bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y);
    1011 bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    1012 
    1013 psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w);
    1014 
    1015 double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
    1016 
    1017 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
    1018                                 const pmSubtractionStampList *stamps,
    1019                                 const pmSubtractionEquationCalculationMode mode)
     978bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
    1020979{
    1021980    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    1022981    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     982
     983    psTimerStart("pmSubtractionSolveEquation");
    1023984
    1024985    // Check inputs
    … …  
    10421003        }
    10431004
     1005        if (stamp->vector->n != numParams) {
     1006            fprintf (stderr, "mismatch length\n");
     1007        }
    10441008        PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
    10451009        PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
    … …  
    10801044        }
    10811045
     1046        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
     1047
    10821048#if 0
    10831049        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
    … …  
    10871053#endif
    10881054
     1055        // XXX TEST : print the matrix & vector
     1056        if (0) {
     1057            for (int iy = 0; iy < sumMatrix->numRows; iy++) {
     1058                for (int ix = 0; ix < sumMatrix->numCols; ix++) {
     1059                    fprintf (stderr, "%e  ", sumMatrix->data.F64[iy][ix]);
     1060                }
     1061                fprintf (stderr, " : %e\n", sumVector->data.F64[iy]);
     1062            }
     1063        }
     1064
     1065        psImage *invMatrix = NULL;
    10891066        psVector *solution = NULL;                       // Solution to equation!
    10901067        solution = psVectorAlloc(numParams, PS_TYPE_F64);
    … …  
    10941071        // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    10951072        // SINGLE solution
    1096         if (1) {
    1097             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    1098         } else {
    1099             psVector *PERM = NULL;
    1100             psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix);
    1101             solution = psMatrixLUSolution(solution, LU, sumVector, PERM);
    1102             psFree (LU);
    1103             psFree (PERM);
    1104         }
     1073# if (1)
     1074        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
     1075        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
     1076# endif
     1077# if (0)
     1078        psMatrixLUSolve(sumMatrixLU, sumVector);
     1079        solution = psMemIncrRefCounter(sumVector);
     1080        invMatrix = psMemIncrRefCounter(sumMatrix);
     1081# endif
     1082# if (0)
     1083        psMatrixGJSolve(sumMatrix, sumVector);
     1084        invMatrix = psMemIncrRefCounter(sumMatrix);
     1085        solution = psMemIncrRefCounter(sumVector);
     1086# endif
    11051087
    11061088# if (0)
    11071089        for (int i = 0; i < solution->n; i++) {
    1108             psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
     1090            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(fabs(invMatrix->data.F64[i][i])));
    11091091        }
    11101092# endif
    11111093
    1112         if (!kernels->solution1) {
    1113             kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
    1114             psVectorInit(kernels->solution1, 0.0);
    1115         }
     1094        // ensure we have a solution vector of the right size
     1095        kernels->solution1    = psVectorRecycle(kernels->solution1,    sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
     1096        kernels->solution1err = psVectorRecycle(kernels->solution1err, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
     1097        psVectorInit(kernels->solution1, 0.0);
     1098        psVectorInit(kernels->solution1err, 0.0);
    11161099
    11171100        int numKernels = kernels->num;
    11181101        int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
    11191102        int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     1103
    11201104        for (int i = 0; i < numKernels * numPoly; i++) {
    11211105            kernels->solution1->data.F64[i] = solution->data.F64[i];
     1106            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]);
    11221107        }
    11231108
    … …  
    11311116        psFree(sumVector);
    11321117        psFree(sumMatrix);
     1118        psFree(invMatrix);
    11331119
    11341120    } else {
    … …  
    11601146        }
    11611147
     1148        pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps);
     1149
    11621150#if 0
    11631151        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
    … …  
    11711159        }
    11721160
     1161        // XXX TEST : print the matrix & vector
     1162        if (0) {
     1163            for (int iy = 0; iy < sumMatrix->numRows; iy++) {
     1164                for (int ix = 0; ix < sumMatrix->numCols; ix++) {
     1165                    fprintf (stderr, "%e  ", sumMatrix->data.F64[iy][ix]);
     1166                }
     1167                fprintf (stderr, " : %e\n", sumVector->data.F64[iy]);
     1168            }
     1169        }
     1170
     1171        psImage *invMatrix = NULL;
    11731172        psVector *solution = NULL;                       // Solution to equation!
    11741173        solution = psVectorAlloc(numParams, PS_TYPE_F64);
    … …  
    11761175
    11771176        // DUAL solution
    1178         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1177# if (1)
     1178        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
     1179        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
     1180# endif
     1181# if (0)
     1182        psMatrixLUSolve(sumMatrix, sumVector);
     1183        solution = psMemIncrRefCounter(sumVector);
     1184        invMatrix = psMemIncrRefCounter(sumMatrix);
     1185# endif
    11791186
    11801187#if (0)
    11811188        for (int i = 0; i < solution->n; i++) {
    1182             fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
     1189            fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
    11831190        }
    11841191#endif
    11851192
    1186         if (!kernels->solution1) {
    1187             kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64);
    1188             psVectorInit (kernels->solution1, 0.0);
    1189         }
    1190         if (!kernels->solution2) {
    1191             kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
    1192             psVectorInit (kernels->solution2, 0.0);
    1193         }
     1193        // XXX TEST: manually set the coeffs to a desired solution
     1194        // solution->data.F64[0] = +1.826;
     1195        // solution->data.F64[1] = -0.115;
     1196        // solution->data.F64[2] =  0.0;
     1197        // solution->data.F64[3] =  0.0;
     1198
     1199        // ensure we have solution vectors of the right size
     1200        kernels->solution1    = psVectorRecycle(kernels->solution1,    numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
     1201        kernels->solution1err = psVectorRecycle(kernels->solution1err, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
     1202        kernels->solution2    = psVectorRecycle(kernels->solution2,    numSolution2,     PS_TYPE_F64); // 1 for norm, 1 for bg
     1203        kernels->solution2err = psVectorRecycle(kernels->solution2err, numSolution2,     PS_TYPE_F64); // 1 for norm, 1 for bg
     1204
     1205        psVectorInit(kernels->solution1, 0.0);
     1206        psVectorInit(kernels->solution1err, 0.0);
     1207        psVectorInit(kernels->solution2, 0.0);
     1208        psVectorInit(kernels->solution2err, 0.0);
    11941209
    11951210        // for DUAL convolution analysis, we apply the normalization to I1 as follows:
    … …  
    12051220            kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue;
    12061221            kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
     1222
     1223            kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]) * stamps->normValue;
     1224            int i2 = i + numSolution1;
     1225            kernels->solution2err->data.F64[i] = sqrt(invMatrix->data.F64[i2][i2]);
    12071226        }
    12081227
    … …  
    12131232        kernels->solution1->data.F64[bgIndex] = 0.0;
    12141233
     1234        psFree(solution);
     1235        psFree(sumVector);
    12151236        psFree(sumMatrix);
    1216         psFree(sumVector);
    1217         psFree(solution);
     1237        psFree(invMatrix);
    12181238    }
    12191239
    … …  
    12241244    if (psTraceGetLevel("psModules.imcombine") >= 7) {
    12251245        for (int i = 0; i < kernels->solution1->n; i++) {
    1226             psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]);
     1246            psTrace("psModules.imcombine", 7, "Solution 1 %d: %f +/- %f\n", i, kernels->solution1->data.F64[i], kernels->solution1err->data.F64[i]);
    12271247        }
    12281248        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    12291249            for (int i = 0; i < kernels->solution2->n; i++) {
    1230                 psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]);
     1250                psTrace("psModules.imcombine", 7, "Solution 2 %d: %f +/- %f\n", i, kernels->solution2->data.F64[i], kernels->solution2err->data.F64[i]);
    12311251            }
    12321252        }
    12331253     }
     1254
     1255    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Solve equation: %f sec", psTimerClear("pmSubtractionSolveEquation"));
    12341256
    12351257    // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
    … …  
    12831305}
    12841306
    1285 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
    1286                                            pmSubtractionKernels *kernels)
     1307// given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq
     1308bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window) {
     1309
     1310# ifndef USE_WEIGHT
     1311    psAssert(weight == NULL, "impossible!");
     1312# endif
     1313# ifndef USE_WINDOW
     1314    psAssert(window == NULL, "impossible!");
     1315# endif
     1316
     1317    int npix = 0;
     1318    float chisqR = 0;
     1319    float chisqD = 0;
     1320
     1321    // get the chisq
     1322    for (int y = residual->yMin; y <= residual->yMax; y++) {
     1323        for (int x = residual->xMin; x <= residual->xMax; x++) {
     1324            float valueR = PS_SQR(residual->kernel[y][x]);
     1325            if (weight) {
     1326                valueR *= weight->kernel[y][x];
     1327            }
     1328            // XXX NOTE: do NOT apply the window to the chisq portions of the calculation (that would bias the chisq)
     1329            chisqR += valueR;
     1330
     1331            float valueD = PS_SQR(difference->kernel[y][x]);
     1332            if (weight) {
     1333                valueD *= weight->kernel[y][x];
     1334            }
     1335            chisqD += valueD;
     1336            npix ++;
     1337        }
     1338    }
     1339    psVectorAppend(chisqRVector, chisqR / npix);
     1340    psVectorAppend(chisqDVector, chisqD / npix);
     1341
     1342    float value1 = 0;
     1343    float value2 = 0;
     1344    float flux2 = 0;
     1345    float fluxX = 0;
     1346    float fluxY = 0;
     1347    float fluxX2 = 0;
     1348    float fluxY2 = 0;
     1349
     1350    float fluxC1 = 0;
     1351    float fluxC2 = 0;
     1352
     1353    float moment = 0;
     1354
     1355    // get the moments from convolved1
     1356    if (convolved1) {
     1357        for (int y = residual->yMin; y <= residual->yMax; y++) {
     1358            for (int x = residual->xMin; x <= residual->xMax; x++) {
     1359                value1  = convolved1->kernel[y][x];
     1360                value2  = PS_SQR(value1);
     1361
     1362                if (window) {
     1363                    value1 *= window->kernel[y][x];
     1364                    value2 *= window->kernel[y][x];
     1365                }
     1366
     1367                fluxC1 += value1;
     1368                flux2  += value2;
     1369                fluxX  += x * value2;
     1370                fluxY  += y * value2;
     1371                // fluxX2 += PS_SQR(x) * value2;
     1372                // fluxY2 += PS_SQR(y) * value2;
     1373                fluxX2 += PS_SQR(x) * value1;
     1374                fluxY2 += PS_SQR(y) * value1;
     1375            }
     1376        }
     1377        // float Mx = fluxX / flux2;
     1378        // float My = fluxY / flux2;
     1379        // float Mxx = fluxX2 / flux2;
     1380        // float Myy = fluxY2 / flux2;
     1381        float Mxx = fluxX2 / fluxC1;
     1382        float Myy = fluxY2 / fluxC1;
     1383
     1384        // 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);
     1385        moment += Mxx + Myy;
     1386    }
     1387
     1388    // get the moments from convolved1
     1389    if (convolved2) {
     1390        for (int y = residual->yMin; y <= residual->yMax; y++) {
     1391            for (int x = residual->xMin; x <= residual->xMax; x++) {
     1392                value1  = convolved2->kernel[y][x];
     1393                value2  = PS_SQR(value1);
     1394
     1395                // XXX NOTE: do NOT apply the weight to the moments calculation
     1396                if (false && weight) {
     1397                    value2 *= weight->kernel[y][x];
     1398                }
     1399                if (window) {
     1400                    value1 *= window->kernel[y][x];
     1401                    value2 *= window->kernel[y][x];
     1402                }
     1403
     1404                fluxC2 += value1;
     1405                flux2  += value2;
     1406                fluxX  += x * value2;
     1407                fluxY  += y * value2;
     1408                // fluxX2 += PS_SQR(x) * value2;
     1409                // fluxY2 += PS_SQR(y) * value2;
     1410                fluxX2 += PS_SQR(x) * value1;
     1411                fluxY2 += PS_SQR(y) * value1;
     1412            }
     1413        }
     1414        // float Mx = fluxX / flux2;
     1415        // float My = fluxY / flux2;
     1416        // float Mxx = fluxX2 / flux2;
     1417        // float Myy = fluxY2 / flux2;
     1418        float Mxx = fluxX2 / fluxC2;
     1419        float Myy = fluxY2 / fluxC2;
     1420
     1421        // 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);
     1422        moment += Mxx + Myy;
     1423    }
     1424
     1425    float flux = fluxC1 + fluxC2;
     1426   
     1427    if (convolved1 && convolved2) {
     1428        moment *= 0.5;
     1429        flux *= 0.5;
     1430    }
     1431    psVectorAppend(momentVector, moment);
     1432    psVectorAppend(fluxesVector, flux);
     1433
     1434    // check that the last appended values are ok:
     1435    int Nelem = fluxesVector->n - 1;
     1436    bool valid = true;
     1437    valid &= isfinite(chisqRVector->data.F32[Nelem]);
     1438    valid &= isfinite(fluxesVector->data.F32[Nelem]);
     1439    valid &= isfinite(momentVector->data.F32[Nelem]);
     1440    if (valid) {
     1441      psVectorAppend(stampMask, 0);
     1442    } else {
     1443      psVectorAppend(stampMask, 0x02);
     1444    }
     1445    return true;
     1446}
     1447
     1448bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch,
     1449                                           pmSubtractionStampList *stamps,
     1450                                           pmSubtractionKernels *kernels)
    12871451{
    12881452    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
    12891453    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
    12901454    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
     1455
     1456    psTimerStart("pmSubtractionCalculateChisqAndMoments");
     1457
     1458    // XXX need to save these somewhere
     1459    psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1460    psVector *chisqD = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1461    psVector *chisqR = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1462    psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1463    psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK);
     1464
     1465    int footprint = stamps->footprint; // Half-size of stamps
     1466    int numKernels = kernels->num;      // Number of kernels
     1467
     1468    psImage *polyValues = NULL;         // Polynomial values
     1469
     1470    // storage for the image (convolved2 is not used in SINGLE mode)
     1471    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     1472    psKernel *difference = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     1473    psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     1474    psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     1475
     1476    int nGood = 0;
     1477    for (int i = 0; i < stamps->num; i++) {
     1478        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
     1479        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
     1480            // mark this stamp as unused (note that we have to append NANs to the other vectors to keep the lengths in sync)
     1481            psVectorAppend(moments, NAN);
     1482            psVectorAppend(fluxes, NAN);
     1483            psVectorAppend(chisqD, NAN);
     1484            psVectorAppend(chisqR, NAN);
     1485            psVectorAppend(stampMask, 0x01);
     1486            continue;
     1487        }
     1488        nGood ++;
     1489
     1490        // Calculate coefficients of the kernel basis functions
     1491        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
     1492        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
     1493        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
     1494
     1495        // Calculate residuals
     1496        psImageInit(residual->image, 0.0);
     1497        psImageInit(difference->image, 0.0);
     1498
     1499        psKernel *weight = NULL;
     1500        psKernel *window = NULL;
     1501   
     1502#ifdef USE_WEIGHT
     1503    weight = stamp->weight;
     1504#endif
     1505#ifdef USE_WINDOW
     1506    window = stamps->window;
     1507#endif
     1508
     1509        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
     1510
     1511            // the single-direction psf match code attempts to find the kernel such that:
     1512            // source * kernel = target.  we need to assign 'source' and 'target' correctly
     1513            // depending on which of image1 or image2 we asked to be convolved.
     1514
     1515            psKernel *target;           // Target postage stamp (convolve source to match the target)
     1516            psKernel *source;           // Source postage stamp (convolve source to match the target)
     1517            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
     1518
     1519            // init the accumulation image
     1520            psImageInit(convolved1->image, 0.0);
     1521
     1522            switch (kernels->mode) {
     1523              case PM_SUBTRACTION_MODE_1:
     1524                target = stamp->image2;
     1525                source = stamp->image1;
     1526                convolutions = stamp->convolutions1;
     1527                break;
     1528              case PM_SUBTRACTION_MODE_2:
     1529                target = stamp->image1;
     1530                source = stamp->image2;
     1531                convolutions = stamp->convolutions2;
     1532                break;
     1533              default:
     1534                psAbort("Unsupported subtraction mode: %x", kernels->mode);
     1535            }
     1536
     1537            // generate the convolved source image (sum over kernels)
     1538            for (int j = 0; j < numKernels; j++) {
     1539                psKernel *convolution = convolutions->data[j]; // Convolution
     1540                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
     1541                for (int y = - footprint; y <= footprint; y++) {
     1542                    for (int x = - footprint; x <= footprint; x++) {
     1543                        convolved1->kernel[y][x] += convolution->kernel[y][x] * coefficient;
     1544                    }
     1545                }
     1546            }
     1547
     1548            // Generate the difference, residual, and convolved source images.  Note the we
     1549            // accumulate the convolution of (A-B), so we need to replace it to generate the
     1550            // images of the convolved source image.
     1551            for (int y = - footprint; y <= footprint; y++) {
     1552                for (int x = - footprint; x <= footprint; x++) {
     1553                    difference->kernel[y][x] = target->kernel[y][x] - source->kernel[y][x] * norm - background;
     1554                    residual->kernel[y][x] = difference->kernel[y][x] - convolved1->kernel[y][x];
     1555                    convolved1->kernel[y][x] += source->kernel[y][x] * norm;
     1556                }
     1557            }
     1558
     1559            // XXX if we want to have a weight and window, we'll need to pass through to here
     1560            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, NULL, difference, residual, weight, window);
     1561
     1562        } else {
     1563
     1564            // Dual convolution
     1565            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
     1566            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
     1567            psKernel *image1 = stamp->image1; // The first image
     1568            psKernel *image2 = stamp->image2; // The second image
     1569
     1570            // init the accumulation images
     1571            psImageInit(convolved1->image, 0.0);
     1572            psImageInit(convolved2->image, 0.0);
     1573
     1574            for (int j = 0; j < numKernels; j++) {
     1575                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
     1576                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
     1577                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
     1578                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
     1579
     1580                for (int y = - footprint; y <= footprint; y++) {
     1581                    for (int x = - footprint; x <= footprint; x++) {
     1582                        // NOTE sign for coeff2
     1583                        convolved1->kernel[y][x] += +conv1->kernel[y][x] * coeff1;
     1584                        convolved2->kernel[y][x] += -conv2->kernel[y][x] * coeff2;
     1585                    }
     1586                }
     1587            }
     1588
     1589            // Generate the difference, residual, and convolved source images.  Note the we
     1590            // accumulate the convolutions of (A-B), so we need to replace (A or B) to generate
     1591            // the images of the convolved source images.
     1592            for (int y = - footprint; y <= footprint; y++) {
     1593                for (int x = - footprint; x <= footprint; x++) {
     1594                    difference->kernel[y][x] = image2->kernel[y][x] - image1->kernel[y][x] * norm - background;
     1595                    residual->kernel[y][x] = difference->kernel[y][x] + convolved2->kernel[y][x] - convolved1->kernel[y][x];
     1596                    convolved1->kernel[y][x] += image1->kernel[y][x] * norm;
     1597                    convolved2->kernel[y][x] += image2->kernel[y][x];
     1598                }
     1599            }
     1600
     1601            if (0) {
     1602                psFitsWriteImageSimple("conv1.fits", convolved1->image, NULL);
     1603                psFitsWriteImageSimple("conv2.fits", convolved2->image, NULL);
     1604                psFitsWriteImageSimple("resid.fits", residual->image,   NULL);
     1605                pmVisualAskUser(NULL);
     1606            }
     1607
     1608            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, convolved2, difference, residual, weight, window);
     1609        }
     1610    }
     1611
     1612    // find the mean chisq and mean moment
     1613    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     1614    psVectorStats (stats, chisqD, NULL, stampMask, 0xff);
     1615    float chisqDValue = stats->sampleMean;
     1616
     1617    psStatsInit(stats);
     1618    psVectorStats (stats, chisqR, NULL, stampMask, 0xff);
     1619    float chisqRValue = stats->sampleMean;
     1620
     1621    psStatsInit(stats);
     1622    psVectorStats (stats, moments, NULL, stampMask, 0xff);
     1623    float momentValue = stats->sampleMean;
     1624
     1625    double sumKernel1 = 0.0, sumKernel2 = 0.0; // Sum of the kernel
     1626
     1627    // calculate the variance contribution from this smoothing kernel
     1628    psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, false);
     1629    for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) {
     1630        for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) {
     1631            if (!isfinite(modelKernel->kernel[y][x])) {
     1632                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y);
     1633                return NULL;
     1634            }
     1635            sumKernel1 += PS_SQR(modelKernel->kernel[y][x]);
     1636        }
     1637    }
     1638    psFree (modelKernel);
     1639
     1640    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     1641        psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, true);
     1642        for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) {
     1643            for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) {
     1644                if (!isfinite(modelKernel->kernel[y][x])) {
     1645                    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y);
     1646                    return NULL;
     1647                }
     1648                sumKernel2 += PS_SQR(modelKernel->kernel[y][x]);
     1649            }
     1650        }
     1651        psFree (modelKernel);
     1652    } else {
     1653        sumKernel2 = 1.0;
     1654    }
     1655
     1656    // if we modify the chisq value by the (sumKernel1 + sumKernel2), we account for the
     1657    // smoothing coming from larger kernels adding additional spatial fit terms should be
     1658    // penalized by increasing the score somewhat.  the 0.01 value is not well-chosen.
     1659    float orderFactor = 0.01 * kernels->spatialOrder;
     1660    float score = 2.0 * chisqRValue / (sumKernel1 + sumKernel2) + orderFactor;
     1661    psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, chisqD: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqRValue, chisqDValue, momentValue, sumKernel1, sumKernel2, score);
     1662
     1663    // save this result if it is the first or the best (skip if bestMatch is NULL)
     1664    if (bestMatch) {
     1665        pmSubtractionQuality *match = *bestMatch;
     1666        bool keep = false;
     1667        if (match == NULL) {
     1668            *bestMatch = match = pmSubtractionQualityAlloc();
     1669            keep = true;
     1670        } else {
     1671            if (score < match->score) {
     1672                psFree(match->fluxes);
     1673                psFree(match->chisq);
     1674                psFree(match->moments);
     1675                psFree(match->stampMask);
     1676                keep = true;
     1677            }
     1678        }
     1679        if (keep) {
     1680            psLogMsg("psModules.imcombine", PS_LOG_INFO, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score);
     1681            match->score        = score;
     1682            match->spatialOrder = kernels->spatialOrder;
     1683            match->mode         = kernels->mode;
     1684            match->nGood        = nGood;
     1685            match->fluxes       = psMemIncrRefCounter(fluxes);
     1686            match->chisq        = psMemIncrRefCounter(chisqR);
     1687            match->moments      = psMemIncrRefCounter(moments);
     1688            match->stampMask    = psMemIncrRefCounter(stampMask);
     1689        }           
     1690    }
     1691
     1692    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisqR, moments);
     1693
     1694    psFree(stats);
     1695    psFree(chisqR);
     1696    psFree(chisqD);
     1697    psFree(fluxes);
     1698    psFree(moments);
     1699    psFree(stampMask);
     1700
     1701    psFree(residual);
     1702    psFree(difference);
     1703    psFree(convolved1);
     1704    psFree(convolved2);
     1705    psFree(polyValues);
     1706
     1707    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Chisq and Moments: %f sec", psTimerClear("pmSubtractionCalculateChisqAndMoments"));
     1708
     1709    return true;
     1710}
     1711
     1712// XXX for now, let's not use this, and let's instead just use values from pmSubtractionCalculateChisqAndMoments
     1713psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
     1714{
     1715    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
     1716    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
     1717    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
     1718
     1719    psTimerStart("pmSubtractionCalculateDeviations");
    12911720
    12921721    psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps
    … …  
    12981727    psImage *polyValues = NULL;         // Polynomial values
    12991728    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
    1300 
    1301     // set up holding images for the visualization
    1302     pmSubtractionVisualShowFitInit (stamps);
    13031729
    13041730    psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    … …  
    14011827            }
    14021828
    1403             // XXX visualize the target, source, convolution and residual
    1404             pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
    1405 
    14061829            for (int y = - footprint; y <= footprint; y++) {
    14071830                for (int x = - footprint; x <= footprint; x++) {
    … …  
    14381861            }
    14391862
    1440             // XXX visualize the target, source, convolution and residual
    1441             pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
    1442 
    14431863            for (int y = - footprint; y <= footprint; y++) {
    14441864                for (int x = - footprint; x <= footprint; x++) {
    … …  
    14551875        }
    14561876
     1877        double flux = 0.0;
    14571878        double deviation = 0.0;         // Sum of differences
    14581879        for (int y = - footprint; y <= footprint; y++) {
    … …  
    14601881                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
    14611882                deviation += dev;
    1462 #ifdef TESTING
    1463                 residual->kernel[y][x] = dev;
    1464 #endif
     1883                flux += stamp->image1->kernel[y][x] + stamp->image2->kernel[y][x];
    14651884            }
    14661885        }
    14671886        deviations->data.F32[i] = devNorm * deviation;
    1468         psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    1469                 i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
     1887        psTrace("psModules.imcombine", 5, "Deviation and Flux for stamp %d (%d,%d): %f %f\n",
     1888                i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i], flux);
    14701889        psStringAppend(&log, "Stamp %d (%d,%d): %f\n",
    14711890                       i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
    14721891        if (!isfinite(deviations->data.F32[i])) {
    14731892            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    1474             psTrace("psModules.imcombine", 5,
    1475                     "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
    1476                     i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
     1893            psTrace("psModules.imcombine", 5, "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
    14771894            continue;
    14781895        }
    1479 
    1480 #ifdef TESTING
    1481         {
    1482             psString filename = NULL;
    1483             psStringAppend(&filename, "resid_%03d.fits", i);
    1484             psFits *fits = psFitsOpen(filename, "w");
    1485             psFree(filename);
    1486             psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
    1487             psFitsClose(fits);
    1488         }
    1489         if (stamp->image1) {
    1490             psString filename = NULL;
    1491             psStringAppend(&filename, "stamp_image1_%03d.fits", i);
    1492             psFits *fits = psFitsOpen(filename, "w");
    1493             psFree(filename);
    1494             psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
    1495             psFitsClose(fits);
    1496         }
    1497         if (stamp->image2) {
    1498             psString filename = NULL;
    1499             psStringAppend(&filename, "stamp_image2_%03d.fits", i);
    1500             psFits *fits = psFitsOpen(filename, "w");
    1501             psFree(filename);
    1502             psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
    1503             psFitsClose(fits);
    1504         }
    1505         if (stamp->weight) {
    1506             psString filename = NULL;
    1507             psStringAppend(&filename, "stamp_weight_%03d.fits", i);
    1508             psFits *fits = psFitsOpen(filename, "w");
    1509             psFree(filename);
    1510             psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
    1511             psFitsClose(fits);
    1512         }
    1513 #endif
    1514 
    15151896    }
    15161897
    15171898    psFree(keepStamps);
    15181899
    1519     psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
     1900    psLogMsg("psModules.imcombine", PS_LOG_MINUTIA, "%s", log);
    15201901    psFree(log);
    15211902
    … …  
    15271908        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
    15281909
    1529         pmSubtractionVisualShowFit(norm);
    1530         pmSubtractionVisualPlotFit(kernels);
    1531 
    15321910        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    15331911        psVectorStats (stats, fResSigma, NULL, NULL, 0);
    … …  
    15601938    psFree(polyValues);
    15611939
     1940    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Deviations: %f sec", psTimerClear("pmSubtractionCalculateDeviations"));
     1941
    15621942    return deviations;
    1563 }
    1564 
    1565 // we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w)
    1566 psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) {
    1567 
    1568     psAssert (w->n == U->numCols, "w and U dimensions do not match");
    1569 
    1570     // wUt has dimensions transposed relative to Ut.
    1571     psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64);
    1572     psImageInit (wUt, 0.0);
    1573 
    1574     for (int i = 0; i < wUt->numCols; i++) {
    1575         for (int j = 0; j < wUt->numRows; j++) {
    1576             if (!isfinite(w->data.F64[j])) continue;
    1577             if (w->data.F64[j] == 0.0) continue;
    1578             wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
    1579         }
    1580     }
    1581     return wUt;
    1582 }
    1583 
    1584 // XXX this is just standard matrix multiplication: use psMatrixMultiply?
    1585 psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) {
    1586 
    1587     psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match");
    1588 
    1589     psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64);
    1590 
    1591     for (int i = 0; i < Ainv->numCols; i++) {
    1592         for (int j = 0; j < Ainv->numRows; j++) {
    1593             double sum = 0.0;
    1594             for (int k = 0; k < V->numCols; k++) {
    1595                 sum += V->data.F64[j][k] * wUt->data.F64[k][i];
    1596             }
    1597             Ainv->data.F64[j][i] = sum;
    1598         }
    1599     }
    1600     return Ainv;
    1601 }
    1602 
    1603 // we are supplied U, not Ut
    1604 bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) {
    1605 
    1606     psAssert (U->numRows == B->n, "U and B dimensions do not match");
    1607 
    1608     UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64);
    1609 
    1610     for (int i = 0; i < U->numCols; i++) {
    1611         double sum = 0.0;
    1612         for (int j = 0; j < U->numRows; j++) {
    1613             sum += B->data.F64[j] * U->data.F64[j][i];
    1614         }
    1615         UtB[0]->data.F64[i] = sum;
    1616     }
    1617     return true;
    1618 }
    1619 
    1620 // w is diagonal
    1621 bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) {
    1622 
    1623     psAssert (w->n == UtB->n, "w and UtB dimensions do not match");
    1624 
    1625     // wUt has dimensions transposed relative to Ut.
    1626     wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64);
    1627     psVectorInit (wUtB[0], 0.0);
    1628 
    1629     for (int i = 0; i < w->n; i++) {
    1630         if (!isfinite(w->data.F64[i])) continue;
    1631         if (w->data.F64[i] == 0.0) continue;
    1632         wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
    1633     }
    1634     return true;
    1635 }
    1636 
    1637 // this is basically matrix * vector
    1638 bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) {
    1639 
    1640     psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match");
    1641 
    1642     VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64);
    1643 
    1644     for (int j = 0; j < V->numRows; j++) {
    1645         double sum = 0.0;
    1646         for (int i = 0; i < V->numCols; i++) {
    1647             sum += V->data.F64[j][i] * wUtB->data.F64[i];
    1648         }
    1649         VwUtB[0]->data.F64[j] = sum;
    1650     }
    1651     return true;
    1652 }
    1653 
    1654 // this is basically matrix * vector
    1655 bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) {
    1656 
    1657     psAssert (A->numCols == x->n, "A and x dimensions do not match");
    1658 
    1659     B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64);
    1660 
    1661     for (int j = 0; j < A->numRows; j++) {
    1662         double sum = 0.0;
    1663         for (int i = 0; i < A->numCols; i++) {
    1664             sum += A->data.F64[j][i] * x->data.F64[i];
    1665         }
    1666         B[0]->data.F64[j] = sum;
    1667     }
    1668     return true;
    1669 }
    1670 
    1671 // this is basically Vector * vector
    1672 bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) {
    1673 
    1674     psAssert (x->n == y->n, "x and y dimensions do not match");
    1675 
    1676     double sum = 0.0;
    1677     for (int i = 0; i < x->n; i++) {
    1678         sum += x->data.F64[i] * y->data.F64[i];
    1679     }
    1680     *value = sum;
    1681     return true;
    1682 }
    1683 
    1684 bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
    1685 
    1686     int footprint = stamps->footprint; // Half-size of stamps
    1687 
    1688     double sum = 0.0;
    1689     for (int i = 0; i < stamps->num; i++) {
    1690 
    1691         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1692         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1693 
    1694         psKernel *weight = NULL;
    1695         psKernel *window = NULL;
    1696         psKernel *input = NULL;
    1697 
    1698 #ifdef USE_WEIGHT
    1699         weight = stamp->weight;
    1700 #endif
    1701 #ifdef USE_WINDOW
    1702         window = stamps->window;
    1703 #endif
    1704 
    1705         switch (kernels->mode) {
    1706             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1707           case PM_SUBTRACTION_MODE_1:
    1708             input = stamp->image2;
    1709             break;
    1710           case PM_SUBTRACTION_MODE_2:
    1711             input = stamp->image1;
    1712             break;
    1713           default:
    1714             psAbort ("programming error");
    1715         }
    1716 
    1717         for (int y = - footprint; y <= footprint; y++) {
    1718             for (int x = - footprint; x <= footprint; x++) {
    1719                 double in = input->kernel[y][x];
    1720                 double value = in*in;
    1721                 if (weight) {
    1722                     float wtVal = weight->kernel[y][x];
    1723                     value *= wtVal;
    1724                 }
    1725                 if (window) {
    1726                     float  winVal = window->kernel[y][x];
    1727                     value *= winVal;
    1728                 }
    1729                 sum += value;
    1730             }
    1731         }
    1732     }
    1733     *y2 = sum;
    1734     return true;
    1735 }
    1736 
    1737 double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
    1738 
    1739     int footprint = stamps->footprint; // Half-size of stamps
    1740     int numKernels = kernels->num;      // Number of kernels
    1741 
    1742     double sum = 0.0;
    1743 
    1744     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
    1745     psImageInit(residual->image, 0.0);
    1746 
    1747     psImage *polyValues = NULL;         // Polynomial values
    1748 
    1749     for (int i = 0; i < stamps->num; i++) {
    1750 
    1751         pmSubtractionStamp *stamp = stamps->stamps->data[i];
    1752         if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
    1753 
    1754         psKernel *weight = NULL;
    1755         psKernel *window = NULL;
    1756         psKernel *target = NULL;
    1757         psKernel *source = NULL;
    1758 
    1759         psArray *convolutions = NULL;
    1760 
    1761 #ifdef USE_WEIGHT
    1762         weight = stamp->weight;
    1763 #endif
    1764 #ifdef USE_WINDOW
    1765         window = stamps->window;
    1766 #endif
    1767 
    1768         switch (kernels->mode) {
    1769             // MODE_1 : convolve image 1 to match image 2 (and vice versa)
    1770           case PM_SUBTRACTION_MODE_1:
    1771             target = stamp->image2;
    1772             source = stamp->image1;
    1773             convolutions = stamp->convolutions1;
    1774             break;
    1775           case PM_SUBTRACTION_MODE_2:
    1776             target = stamp->image1;
    1777             source = stamp->image2;
    1778             convolutions = stamp->convolutions2;
    1779             break;
    1780           default:
    1781             psAbort ("programming error");
    1782         }
    1783 
    1784         // Calculate coefficients of the kernel basis functions
    1785         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
    1786         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
    1787         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    1788 
    1789         psImageInit(residual->image, 0.0);
    1790         for (int j = 0; j < numKernels; j++) {
    1791             psKernel *convolution = convolutions->data[j]; // Convolution
    1792             double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
    1793             for (int y = - footprint; y <= footprint; y++) {
    1794                 for (int x = - footprint; x <= footprint; x++) {
    1795                     residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
    1796                 }
    1797             }
    1798         }
    1799 
    1800         for (int y = - footprint; y <= footprint; y++) {
    1801             for (int x = - footprint; x <= footprint; x++) {
    1802                 double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
    1803                 double value = PS_SQR(resid);
    1804                 if (weight) {
    1805                     float wtVal = weight->kernel[y][x];
    1806                     value *= wtVal;
    1807                 }
    1808                 if (window) {
    1809                     float  winVal = window->kernel[y][x];
    1810                     value *= winVal;
    1811                 }
    1812                 sum += value;
    1813             }
    1814         }
    1815     }
    1816     psFree (polyValues);
    1817     psFree (residual);
    1818 
    1819     return sum;
    1820 }
    1821 
    1822 bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) {
    1823 
    1824     for (int i = 0; i < w->n; i++) {
    1825         wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
    1826     }
    1827     return true;
    1828 }
    1829 
    1830 // we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w)
    1831 psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) {
    1832 
    1833     psAssert (w->n == V->numCols, "w and U dimensions do not match");
    1834 
    1835     psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
    1836     psImageInit (Vn, 0.0);
    1837 
    1838     // generate Vn = V * w^{-1}
    1839     for (int j = 0; j < Vn->numRows; j++) {
    1840         for (int i = 0; i < Vn->numCols; i++) {
    1841             if (!isfinite(w->data.F64[i])) continue;
    1842             if (w->data.F64[i] == 0.0) continue;
    1843             Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
    1844         }
    1845     }
    1846 
    1847     psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
    1848     psImageInit (Xvar, 0.0);
    1849 
    1850     // generate Xvar = Vn * Vn^T
    1851     for (int j = 0; j < Vn->numRows; j++) {
    1852         for (int i = 0; i < Vn->numCols; i++) {
    1853             double sum = 0.0;
    1854             for (int k = 0; k < Vn->numCols; k++) {
    1855                 sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
    1856             }
    1857             Xvar->data.F64[j][i] = sum;
    1858         }
    1859     }
    1860     return Xvar;
    18611943}
    18621944
    … …  
    18641946// of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix
    18651947// multiplication is: A_k,j * B_i,k = C_i,j
    1866 
    18671948
    18681949bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) {
Note: See TracChangeset for help on using the changeset viewer.