IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 31, 2010, 5:00:42 PM (16 years ago)
Author:
eugene
Message:

updates relative to 20091201, fixes for all psphot variants

Location:
branches/eam_branches/psModules.stack.20100120
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/psModules.stack.20100120

  • branches/eam_branches/psModules.stack.20100120/src/imcombine/pmSubtractionEquation.c

    r26667 r26747  
    2828static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
    2929                                  psVector *vector, // Least-squares vector, updated
     30                                  double *norm,     // Normalisation, updated
    3031                                  const psKernel *input, // Input image (target)
    3132                                  const psKernel *reference, // Reference image (convolution source)
    … …  
    3637                                  const psImage *polyValues, // Spatial polynomial values
    3738                                  int footprint, // (Half-)Size of stamp
     39                                  int normWindow, // Window (half-)size for normalisation measurement
    3840                                  const pmSubtractionEquationCalculationMode mode
    3941                                  )
    … …  
    156158            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    157159                vector->data.F64[iIndex] = sumIC * poly[iTerm];
    158                 // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
    159                 // instead, calculate A - \sum(k x B), with full hermitians
    160                 if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
     160                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
    161161                    // subtract norm * sumRC * poly[iTerm]
    162162                    psAssert (kernels->solution1, "programming error: define solution first!");
    … …  
    174174    double sumR = 0.0;                  // Sum of the reference
    175175    double sumI = 0.0;                  // Sum of the input
     176    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
    176177    for (int y = - footprint; y <= footprint; y++) {
    177178        for (int x = - footprint; x <= footprint; x++) {
    … …  
    181182            double rr = PS_SQR(ref);
    182183            double one = 1.0;
     184
     185            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     186                normI1 += ref;
     187                normI2 += in;
     188            }
     189
    183190            if (weight) {
    184191                float wtVal = weight->kernel[y][x];
    … …  
    204211        }
    205212    }
     213
     214    *norm = normI2 / normI1;
     215
    206216    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    207217        matrix->data.F64[normIndex][normIndex] = sumRR;
    … …  
    217227        matrix->data.F64[bgIndex][normIndex] = sumR;
    218228    }
     229
     230    // check for any NAN values in the result, skip if found:
     231    for (int iy = 0; iy < matrix->numRows; iy++) {
     232        for (int ix = 0; ix < matrix->numCols; ix++) {
     233            if (!isfinite(matrix->data.F64[iy][ix])) {
     234                fprintf (stderr, "WARNING: NAN in matrix\n");
     235                return false;
     236            }
     237        }
     238    }
     239    for (int ix = 0; ix < vector->n; ix++) {
     240        if (!isfinite(vector->data.F64[ix])) {
     241            fprintf (stderr, "WARNING: NAN in vector\n");
     242            return false;
     243        }
     244    }
     245
    219246    return true;
    220247}
    … …  
    224251static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated
    225252                                      psVector *vector, // Least-squares vector, updated
     253                                      double *norm,     // Normalisation, updated
    226254                                      const psKernel *image1, // Image 1
    227255                                      const psKernel *image2, // Image 2
    … …  
    232260                                      const pmSubtractionKernels *kernels, // Kernels
    233261                                      const psImage *polyValues, // Spatial polynomial values
    234                                       int footprint // (Half-)Size of stamp
     262                                      int footprint, // (Half-)Size of stamp
     263                                      int normWindow, // Window (half-)size for normalisation measurement
     264                                      const pmSubtractionEquationCalculationMode mode
    235265                                      )
    236266{
    … …  
    272302    }
    273303
     304
     305    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
     306    // choose to calculate
     307    psImageInit(matrix, 0.0);
     308    psVectorInit(vector, 1.0);
     309    for (int i = 0; i < matrix->numCols; i++) {
     310        matrix->data.F64[i][i] = 1.0;
     311    }
    274312
    275313    for (int i = 0; i < numKernels; i++) {
    … …  
    306344            }
    307345
    308             // Spatial variation
    309             for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    310                 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    311                     double aa = sumAA * poly2[iTerm][jTerm];
    312                     double bb = sumBB * poly2[iTerm][jTerm];
    313                     double ab = sumAB * poly2[iTerm][jTerm];
    314 
    315                     matrix->data.F64[iIndex][jIndex] = aa;
    316                     matrix->data.F64[jIndex][iIndex] = aa;
    317 
    318                     matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
    319                     matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
    320 
    321                     matrix->data.F64[iIndex][jIndex + numParams] = ab;
    322                     matrix->data.F64[jIndex + numParams][iIndex] = ab;
     346            // Spatial variation of kernel coeffs
     347            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     348                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     349                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     350                        double aa = sumAA * poly2[iTerm][jTerm];
     351                        double bb = sumBB * poly2[iTerm][jTerm];
     352                        double ab = sumAB * poly2[iTerm][jTerm];
     353
     354                        matrix->data.F64[iIndex][jIndex] = aa;
     355                        matrix->data.F64[jIndex][iIndex] = aa;
     356
     357                        matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
     358                        matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
     359
     360                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
     361                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
     362                    }
    323363                }
    324364            }
    … …  
    340380            }
    341381
    342             // Spatial variation
    343             for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
    344                 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
    345                     double ab = sumAB * poly2[iTerm][jTerm];
    346                     matrix->data.F64[iIndex][jIndex + numParams] = ab;
    347                     matrix->data.F64[jIndex + numParams][iIndex] = ab;
     382            // Spatial variation of kernel coeffs
     383            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     384                for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
     385                    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
     386                        double ab = sumAB * poly2[iTerm][jTerm];
     387                        matrix->data.F64[iIndex][jIndex + numParams] = ab;
     388                        matrix->data.F64[jIndex + numParams][iIndex] = ab;
     389                    }
    348390                }
    349391            }
    … …  
    403445            double bi2 = sumBI2 * poly[iTerm];
    404446            double ai1 = sumAI1 * poly[iTerm];
    405             double a = sumA * poly[iTerm];
     447            double a   = sumA * poly[iTerm];
    406448            double bi1 = sumBI1 * poly[iTerm];
    407             double b = sumB * poly[iTerm];
    408 
    409             matrix->data.F64[iIndex][normIndex] = ai1;
    410             matrix->data.F64[normIndex][iIndex] = ai1;
    411             matrix->data.F64[iIndex][bgIndex] = a;
    412             matrix->data.F64[bgIndex][iIndex] = a;
    413             matrix->data.F64[iIndex + numParams][normIndex] = bi1;
    414             matrix->data.F64[normIndex][iIndex + numParams] = bi1;
    415             matrix->data.F64[iIndex + numParams][bgIndex] = b;
    416             matrix->data.F64[bgIndex][iIndex + numParams] = b;
    417             vector->data.F64[iIndex] = ai2;
    418             vector->data.F64[iIndex + numParams] = bi2;
     449            double b   = sumB * poly[iTerm];
     450
     451            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     452                matrix->data.F64[iIndex][normIndex] = ai1;
     453                matrix->data.F64[normIndex][iIndex] = ai1;
     454                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
     455                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
     456            }
     457            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
     458                matrix->data.F64[iIndex][bgIndex] = a;
     459                matrix->data.F64[bgIndex][iIndex] = a;
     460                matrix->data.F64[iIndex + numParams][bgIndex] = b;
     461                matrix->data.F64[bgIndex][iIndex + numParams] = b;
     462            }
     463            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     464                vector->data.F64[iIndex] = ai2;
     465                vector->data.F64[iIndex + numParams] = bi2;
     466                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
     467                    // subtract norm * sumRC * poly[iTerm]
     468                    psAssert (kernels->solution1, "programming error: define solution first!");
     469                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     470                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
     471                    vector->data.F64[iIndex] -= norm * ai1;
     472                    vector->data.F64[iIndex + numParams] -= norm * bi1;
     473                }
     474            }
    419475        }
    420476    }
    … …  
    425481    double sumI2 = 0.0;                 // Sum of I_2 (for vector, background)
    426482    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
     483    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
    427484    for (int y = - footprint; y <= footprint; y++) {
    428485        for (int x = - footprint; x <= footprint; x++) {
    … …  
    433490            double one = 1.0;
    434491            double i1i2 = i1 * i2;
     492
     493            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     494                normI1 += i1;
     495                normI2 += i2;
     496            }
    435497
    436498            if (weight) {
    … …  
    457519        }
    458520    }
    459     matrix->data.F64[bgIndex][normIndex] = sumI1;
    460     matrix->data.F64[normIndex][bgIndex] = sumI1;
    461     matrix->data.F64[normIndex][normIndex] = sumI1I1;
    462     matrix->data.F64[bgIndex][bgIndex] = sum1;
    463     vector->data.F64[bgIndex] = sumI2;
    464     vector->data.F64[normIndex] = sumI1I2;
     521
     522    *norm = normI2 / normI1;
     523
     524    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     525        matrix->data.F64[normIndex][normIndex] = sumI1I1;
     526        vector->data.F64[normIndex] = sumI1I2;
     527    }
     528    if (mode & PM_SUBTRACTION_EQUATION_BG) {
     529        matrix->data.F64[bgIndex][bgIndex] = sum1;
     530        vector->data.F64[bgIndex] = sumI2;
     531    }
     532    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
     533        matrix->data.F64[bgIndex][normIndex] = sumI1;
     534        matrix->data.F64[normIndex][bgIndex] = sumI1;
     535    }
     536
     537    // check for any NAN values in the result, skip if found:
     538    for (int iy = 0; iy < matrix->numRows; iy++) {
     539        for (int ix = 0; ix < matrix->numCols; ix++) {
     540            if (!isfinite(matrix->data.F64[iy][ix])) {
     541                fprintf (stderr, "WARNING: NAN in matrix\n");
     542                return false;
     543            }
     544        }
     545    }
     546    for (int ix = 0; ix < vector->n; ix++) {
     547        if (!isfinite(vector->data.F64[ix])) {
     548            fprintf (stderr, "WARNING: NAN in vector\n");
     549            return false;
     550        }
     551    }
     552
    465553
    466554    return true;
    … …  
    469557#if 1
    470558// Add in penalty term to least-squares vector
    471 static bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
     559bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
    472560                             psVector *vector,                    // Vector to which to add in penalty term
    473561                             const pmSubtractionKernels *kernels, // Kernel parameters
    … …  
    482570    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
    483571    int numKernels = kernels->num; // Number of kernel components
     572    int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations
     573    int numParams = numKernels * numSpatial;                 // Number of kernel parameters
     574
     575    // order is :
     576    // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0]
     577    // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0]
     578    // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1]
     579    // [norm]
     580    // [bg]
     581    // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0]
     582    // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0]
     583    // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1]
     584
    484585    for (int i = 0; i < numKernels; i++) {
    485586        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    … …  
    488589                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
    489590                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
     591                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     592                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
     593                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
     594                    // penalties scale with second moments
     595                    //
     596                }
    490597            }
    491598        }
    … …  
    668775    switch (kernels->mode) {
    669776      case PM_SUBTRACTION_MODE_1:
    670         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
     777        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
    671778                                       weight, window, stamp->convolutions1, kernels,
    672                                        polyValues, footprint, mode);
     779                                       polyValues, footprint, stamps->normWindow, mode);
    673780        break;
    674781      case PM_SUBTRACTION_MODE_2:
    675         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
     782        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
    676783                                       weight, window, stamp->convolutions2, kernels,
    677                                        polyValues, footprint, mode);
     784                                       polyValues, footprint, stamps->normWindow, mode);
    678785        break;
    679786      case PM_SUBTRACTION_MODE_DUAL:
    680         status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2,
     787        status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm,
     788                                           stamp->image1, stamp->image2,
    681789                                           weight, window, stamp->convolutions1, stamp->convolutions2,
    682                                            kernels, polyValues, footprint);
     790                                           kernels, polyValues, footprint, stamps->normWindow, mode);
    683791        break;
    684792      default:
    … …  
    721829}
    722830
    723 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode)
     831bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
     832                                    const pmSubtractionEquationCalculationMode mode)
    724833{
    725834    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
    … …  
    843952        psVectorInit(sumVector, 0.0);
    844953        psImageInit(sumMatrix, 0.0);
    845         int numStamps = 0;              // Number of good stamps
    846         for (int i = 0; i < stamps->num; i++) {
    847             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    848 
    849             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
    850 
    851 #ifdef TESTING
    852               // XXX double-check for NAN in data:
    853                 for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
    854                     for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
    855                         if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
    856                             fprintf (stderr, "WARNING: NAN in matrix\n");
    857                         }
    858                     }
    859                 }
    860                 for (int ix = 0; ix < stamp->vector->n; ix++) {
    861                     if (!isfinite(stamp->vector->data.F64[ix])) {
    862                         fprintf (stderr, "WARNING: NAN in vector\n");
    863                     }
    864                 }
    865 #endif
    866 
    867                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    868                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
    869                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    870                 numStamps++;
    871             } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
    872                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
    873             }
    874         }
    875 
    876 #ifdef TESTING
    877         for (int ix = 0; ix < sumVector->n; ix++) {
    878             if (!isfinite(sumVector->data.F64[ix])) {
    879                 fprintf (stderr, "WARNING: NAN in vector\n");
    880             }
    881         }
    882 #endif
    883 
    884 #if 0
    885         int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    886         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
    887 #endif
    888 
    889 #ifdef TESTING
    890         for (int ix = 0; ix < sumVector->n; ix++) {
    891             if (!isfinite(sumVector->data.F64[ix])) {
    892                 fprintf (stderr, "WARNING: NAN in vector\n");
    893             }
    894         }
    895         {
    896             psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
    897             psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
    898             psFree(inverse);
    899         }
    900         {
    901             psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
    902             psImage *Xt = psMatrixTranspose(NULL, X);
    903             psImage *XtX = psMatrixMultiply(NULL, Xt, X);
    904             psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
    905             psFree(X);
    906             psFree(Xt);
    907             psFree(XtX);
    908         }
    909 #endif
    910 
    911 # define SVD_ANALYSIS 0
    912 # define COEFF_SIG 0.0
    913 # define SVD_TOL 0.0
    914 
    915         psVector *solution = NULL;
    916         psVector *solutionErr = NULL;
    917 
    918 #ifdef TESTING
    919         psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
    920         psVectorWriteFile ("B.dat", sumVector);
    921 #endif
    922 
    923         // Use SVD to determine the kernel coeffs (and validate)
    924         if (SVD_ANALYSIS) {
    925 
    926             // We have sumVector and sumMatrix.  we are trying to solve the following equation:
    927             // sumMatrix * x = sumVector.
    928 
    929             // we can use any standard matrix inversion to solve this.  However, the basis
    930             // functions in general have substantial correlation, so that the solution may be
    931             // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
    932             // system of equations may be statistically ill-conditioned.  Noise in the image
    933             // will drive insignificant, but correlated, terms in the solution.  To avoid these
    934             // problems, we can use SVD to identify numerically unconstrained values and to
    935             // avoid statistically badly determined value.
    936 
    937             // A = sumMatrix, B = sumVector
    938             // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
    939             // x = V (1/w) (U^T B)
    940             // \sigma_x = sqrt(diag(A^{-1}))
    941             // solve for x and A^{-1} to get x & dx
    942             // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
    943             // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    944 
    945             // If I use the SVD trick to re-condition the matrix, I need to break out the
    946             // kernel and normalization terms from the background term.
    947             // XXX is this true?  or was this due to an error in the analysis?
    948 
    949             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    950 
    951             // now pull out the kernel elements into their own square matrix
    952             psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
    953             psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
    954 
    955             for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
    956                 if (ix == bgIndex) continue;
    957                 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
    958                     if (iy == bgIndex) continue;
    959                     kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
    960                     ky++;
    961                 }
    962                 kernelVector->data.F64[kx] = sumVector->data.F64[ix];
    963                 kx++;
    964             }
    965 
    966             psImage *U = NULL;
    967             psImage *V = NULL;
    968             psVector *w = NULL;
    969             if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
    970                 psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
    971                 return NULL;
    972             }
    973 
    974             // calculate A_inverse:
    975             // Ainv = V * w * U^T
    976             psImage *wUt  = p_pmSubSolve_wUt (w, U);
    977             psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
    978             psImage *Xvar = NULL;
    979             psFree (wUt);
    980 
    981 # ifdef TESTING
    982             // kernel terms:
    983             for (int i = 0; i < w->n; i++) {
    984                 fprintf (stderr, "w: %f\n", w->data.F64[i]);
    985             }
    986 # endif
    987             // loop over w adding in more and more of the values until chisquare is no longer
    988             // dropping significantly.
    989             // XXX this does not seem to work very well: we seem to need all terms even for
    990             // simple cases...
    991 
    992             psVector *Xsvd = NULL;
    993             {
    994                 psVector *Ax = NULL;
    995                 psVector *UtB = NULL;
    996                 psVector *wUtB = NULL;
    997 
    998                 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
    999                 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
    1000                 psVectorInit (wMask, 1); // start by masking everything
    1001 
    1002                 double chiSquareLast = NAN;
    1003                 int maxWeight = 0;
    1004 
    1005                 double Axx, Bx, y2;
    1006 
    1007                 // XXX this is an attempt to exclude insignificant modes.
    1008                 // it was not successful with the ISIS kernel set: removing even
    1009                 // the least significant mode leaves additional ringing / noise
    1010                 // because the terms are so coupled.
    1011                 for (int k = 0; false && (k < w->n); k++) {
    1012 
    1013                     // unmask the k-th weight
    1014                     wMask->data.U8[k] = 0;
    1015                     p_pmSubSolve_SetWeights(wApply, w, wMask);
    1016 
    1017                     // solve for x:
    1018                     // x = V * w * (U^T * B)
    1019                     p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1020                     p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1021                     p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1022 
    1023                     // chi-square for this system of equations:
    1024                     // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1025                     // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1026                     p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1027                     p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1028                     p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1029                     p_pmSubSolve_y2 (&y2, kernels, stamps);
    1030 
    1031                     // apparently, this works (compare with the brute force value below
    1032                     double chiSquare = Axx - 2.0*Bx + y2;
    1033                     double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
    1034                     chiSquareLast = chiSquare;
    1035 
    1036                     // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
    1037                     if (k && !maxWeight && (deltaChi < 1.0)) {
    1038                         maxWeight = k;
    1039                     }
    1040                 }
    1041 
    1042                 // keep all terms or we get extra ringing
    1043                 maxWeight = w->n;
    1044                 psVectorInit (wMask, 1);
    1045                 for (int k = 0; k < maxWeight; k++) {
    1046                     wMask->data.U8[k] = 0;
    1047                 }
    1048                 p_pmSubSolve_SetWeights(wApply, w, wMask);
    1049 
    1050                 // solve for x:
    1051                 // x = V * w * (U^T * B)
    1052                 p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1053                 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1054                 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1055 
    1056                 // chi-square for this system of equations:
    1057                 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1058                 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1059                 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1060                 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1061                 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1062                 p_pmSubSolve_y2 (&y2, kernels, stamps);
    1063 
    1064                 // apparently, this works (compare with the brute force value below
    1065                 double chiSquare = Axx - 2.0*Bx + y2;
    1066                 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
    1067 
    1068                 // re-calculate A^{-1} to get new variances:
    1069                 // Ainv = V * w * U^T
    1070                 // XXX since we keep all terms, this is identical to Ainv
    1071                 psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
    1072                 Xvar = p_pmSubSolve_VwUt (V, wUt);
    1073                 psFree (wUt);
    1074 
    1075                 psFree (Ax);
    1076                 psFree (UtB);
    1077                 psFree (wUtB);
    1078                 psFree (wApply);
    1079                 psFree (wMask);
    1080             }
    1081 
    1082             // copy the kernel solutions to the full solution vector:
    1083             solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1084             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1085 
    1086             for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
    1087                 if (ix == bgIndex) {
    1088                     solution->data.F64[ix] = 0;
    1089                     solutionErr->data.F64[ix] = 0.001;
    1090                     continue;
    1091                 }
    1092                 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
    1093                 solution->data.F64[ix] = Xsvd->data.F64[kx];
    1094                 kx++;
    1095             }
    1096 
    1097             psFree (kernelMatrix);
    1098             psFree (kernelVector);
    1099 
    1100             psFree (U);
    1101             psFree (V);
    1102             psFree (w);
    1103 
    1104             psFree (Ainv);
    1105             psFree (Xsvd);
    1106         } else {
    1107             psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    1108             psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1109             if (!luMatrix) {
    1110                 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
    1111                 psFree(solution);
    1112                 psFree(sumVector);
    1113                 psFree(sumMatrix);
    1114                 psFree(luMatrix);
    1115                 psFree(permutation);
    1116                 return NULL;
    1117             }
    1118 
    1119             solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1120             psFree(luMatrix);
    1121             psFree(permutation);
    1122             if (!solution) {
    1123                 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    1124                 psFree(solution);
    1125                 psFree(sumVector);
    1126                 psFree(sumMatrix);
    1127                 return NULL;
    1128             }
    1129 
    1130             // XXX LUD does not provide A^{-1}?  fake the error for now
    1131             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1132             for (int ix = 0; ix < sumVector->n; ix++) {
    1133                 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
    1134             }
    1135         }
    1136 
    1137         if (!kernels->solution1) {
    1138             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1139             psVectorInit (kernels->solution1, 0.0);
    1140         }
    1141 
    1142         // only update the solutions that we chose to calculate:
    1143         if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    1144             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1145             kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
    1146         }
    1147         if (mode & PM_SUBTRACTION_EQUATION_BG) {
    1148             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1149             kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
    1150         }
    1151         if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    1152             int numKernels = kernels->num;
    1153             int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
    1154             int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    1155             for (int i = 0; i < numKernels * numPoly; i++) {
    1156                 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
    1157                 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
    1158                     // XXX fprintf (stderr, "drop\n");
    1159                     kernels->solution1->data.F64[i] = 0.0;
    1160                 } else {
    1161                     // XXX fprintf (stderr, "keep\n");
    1162                     kernels->solution1->data.F64[i] = solution->data.F64[i];
    1163                 }
    1164             }
    1165         }
    1166         // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
    1167         // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
    1168 
    1169         psFree(solution);
    1170         psFree(sumVector);
    1171         psFree(sumMatrix);
    1172 
    1173 #ifdef TESTING
    1174         // XXX double-check for NAN in data:
    1175         for (int ix = 0; ix < kernels->solution1->n; ix++) {
    1176             if (!isfinite(kernels->solution1->data.F64[ix])) {
    1177                 fprintf (stderr, "WARNING: NAN in vector\n");
    1178             }
    1179         }
    1180 #endif
    1181 
    1182     } else {
    1183         // Dual convolution solution
    1184 
    1185         // Accumulation of stamp matrices/vectors
    1186         psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    1187         psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
    1188         psImageInit(sumMatrix, 0.0);
    1189         psVectorInit(sumVector, 0.0);
     954
     955        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
    1190956
    1191957        int numStamps = 0;              // Number of good stamps
    … …  
    1195961                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    1196962                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     963                psVectorAppend(norms, stamp->norm);
     964                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
     965                numStamps++;
     966            } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
     967                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
     968            }
     969        }
     970
     971#if 0
     972        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     973        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     974#endif
     975
     976        psVector *solution = NULL;                       // Solution to equation!
     977        solution = psVectorAlloc(numParams, PS_TYPE_F64);
     978        psVectorInit(solution, 0);
     979
     980#if 0
     981        // Regular, straight-forward solution
     982        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     983#else
     984        {
     985            // Solve normalisation and background separately
     986            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     987            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     988
     989            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
     990            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
     991                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
     992                psFree(stats);
     993                psFree(sumMatrix);
     994                psFree(sumVector);
     995                psFree(norms);
     996                return false;
     997            }
     998
     999            double normValue = stats->robustMedian;
     1000            // double bgValue = 0.0;
     1001
     1002            psFree(stats);
     1003
     1004            fprintf(stderr, "Norm: %lf\n", normValue);
     1005
     1006            // Solve kernel components
     1007            for (int i = 0; i < numSolution1; i++) {
     1008                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
     1009
     1010                sumMatrix->data.F64[i][normIndex] = 0.0;
     1011                sumMatrix->data.F64[normIndex][i] = 0.0;
     1012            }
     1013            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
     1014            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1015            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
     1016
     1017            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
     1018            sumVector->data.F64[normIndex] = 0.0;
     1019
     1020            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1021
     1022            solution->data.F64[normIndex] = normValue;
     1023        }
     1024# endif
     1025
     1026        if (!kernels->solution1) {
     1027            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     1028            psVectorInit(kernels->solution1, 0.0);
     1029        }
     1030
     1031        // only update the solutions that we chose to calculate:
     1032        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1033            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1034            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1035        }
     1036        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1037            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1038            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1039        }
     1040        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1041            int numKernels = kernels->num;
     1042            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     1043            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     1044            for (int i = 0; i < numKernels * numPoly; i++) {
     1045                kernels->solution1->data.F64[i] = solution->data.F64[i];
     1046            }
     1047        }
     1048
     1049        psFree(solution);
     1050        psFree(sumVector);
     1051        psFree(sumMatrix);
     1052
     1053#ifdef TESTING
     1054        // XXX double-check for NAN in data:
     1055        for (int ix = 0; ix < kernels->solution1->n; ix++) {
     1056            if (!isfinite(kernels->solution1->data.F64[ix])) {
     1057                fprintf (stderr, "WARNING: NAN in vector\n");
     1058            }
     1059        }
     1060#endif
     1061
     1062    } else {
     1063        // Dual convolution solution
     1064
     1065        // Accumulation of stamp matrices/vectors
     1066        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1067        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
     1068        psImageInit(sumMatrix, 0.0);
     1069        psVectorInit(sumVector, 0.0);
     1070
     1071        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
     1072
     1073        int numStamps = 0;              // Number of good stamps
     1074        for (int i = 0; i < stamps->num; i++) {
     1075            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     1076            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
     1077                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
     1078                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     1079
     1080                psVectorAppend(norms, stamp->norm);
    11971081
    11981082                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    … …  
    12071091
    12081092#if 1
    1209         int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    1210         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     1093        // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     1094        // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     1095
     1096        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1097        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
    12111098#endif
    12121099
    … …  
    12161103
    12171104#if 0
     1105        // Regular, straight-forward solution
     1106        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1107#else
    12181108        {
    1219             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    1220             if (!solution) {
    1221                 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
     1109            // Solve normalisation and background separately
     1110            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1111            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     1112
     1113#if 0
     1114            psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
     1115            psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
     1116
     1117            normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
     1118            normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
     1119            normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
     1120
     1121            normVector->data.F64[0] = sumVector->data.F64[normIndex];
     1122            normVector->data.F64[1] = sumVector->data.F64[bgIndex];
     1123
     1124            psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
     1125
     1126            double normValue = normSolution->data.F64[0];
     1127            double bgValue = normSolution->data.F64[1];
     1128
     1129            psFree(normMatrix);
     1130            psFree(normVector);
     1131            psFree(normSolution);
     1132#endif
     1133
     1134            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
     1135            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
     1136                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
     1137                psFree(stats);
    12221138                psFree(sumMatrix);
    12231139                psFree(sumVector);
    1224                 return NULL;
    1225             }
    1226 
    1227 #ifdef TESTING
    1228             for (int i = 0; i < numParams; i++) {
    1229                 fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
    1230             }
    1231 #endif
    1232 
    1233             int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
    1234             int numKernels = kernels->num; // Number of kernel basis functions
    1235 
    1236 #define MASK_BASIS(INDEX)                                               \
    1237             {                                                           \
    1238                 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \
    1239                     for (int k = 0; k < numParams; k++) {               \
    1240                         sumMatrix->data.F64[index][k] = sumMatrix->data.F64[k][index] = 0.0; \
    1241                     }                                                   \
    1242                     sumMatrix->data.F64[index][index] = 1.0;            \
    1243                     sumVector->data.F64[index] = 0.0;                   \
    1244                 }                                                       \
    1245             }
    1246 
    1247             #define TOL 1.0e-3
    1248             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1249             double norm = fabs(solution->data.F64[normIndex]);  // Normalisation
    1250             double thresh = norm * TOL;                         // Threshold for low parameters
    1251             for (int i = 0; i < numKernels; i++) {
    1252                 // Getting 0th order parameter value.  In the presence of spatial variation, the actual value
    1253                 // of the parameter will vary over the image.  We are in effect getting the value in the
    1254                 // centre of the image.  If we use different polynomial functions (e.g., Chebyshev), we may
    1255                 // have to change this to properly determine the value of the parameter at the centre.
    1256                 double param1 = solution->data.F64[i],
    1257                     param2 = solution->data.F64[numSolution1 + i]; // Parameters of interest
    1258                 bool mask1 = false, mask2 = false;              // Masked the parameter?
    1259                 if (fabs(param1) < thresh) {
    1260                     psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i);
    1261                     MASK_BASIS(i);
    1262                     mask1 = true;
    1263                 }
    1264                 if (fabs(param2) < thresh) {
    1265                     psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i);
    1266                     MASK_BASIS(numSolution1 + i);
    1267                     mask2 = true;
    1268                 }
    1269 
    1270                 if (!mask1 && !mask2) {
    1271                     if (fabs(param1) > fabs(param2)) {
    1272                         psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i);
    1273                         MASK_BASIS(numSolution1 + i);
    1274                     } else {
    1275                         psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i);
    1276                         MASK_BASIS(i);
    1277                     }
    1278                 }
    1279             }
    1280 
    1281 #ifdef TESTING
    1282             psFitsWriteImageSimple ("sumMatrixFix.fits", sumMatrix, NULL);
    1283             psVectorWriteFile("sumVectorFix.dat", sumVector);
    1284 #endif
    1285         }
    1286 #endif /*** kernel-masking block ***/
    1287         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1140                psFree(norms);
     1141                return false;
     1142            }
     1143
     1144            double normValue = stats->robustMedian;
     1145
     1146            psFree(stats);
     1147
     1148            fprintf(stderr, "Norm: %lf\n", normValue);
     1149
     1150            // Solve kernel components
     1151            for (int i = 0; i < numSolution2; i++) {
     1152                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
     1153                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
     1154
     1155                sumMatrix->data.F64[i][normIndex] = 0.0;
     1156                sumMatrix->data.F64[normIndex][i] = 0.0;
     1157
     1158                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
     1159                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
     1160            }
     1161            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
     1162            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1163            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
     1164
     1165            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
     1166
     1167            sumVector->data.F64[normIndex] = 0.0;
     1168
     1169            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1170
     1171            solution->data.F64[normIndex] = normValue;
     1172        }
     1173#endif
     1174
    12881175
    12891176#ifdef TESTING
    … …  
    12961183        psFree(sumVector);
    12971184
     1185        psFree(norms);
     1186
    12981187        if (!kernels->solution1) {
    12991188            kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64);
     1189            psVectorInit (kernels->solution1, 0.0);
    13001190        }
    13011191        if (!kernels->solution2) {
    13021192            kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
    1303         }
     1193            psVectorInit (kernels->solution2, 0.0);
     1194        }
     1195
     1196        // only update the solutions that we chose to calculate:
     1197        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1198            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1199            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1200        }
     1201        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1202            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1203            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1204        }
     1205        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1206            int numKernels = kernels->num;
     1207            for (int i = 0; i < numKernels * numSpatial; i++) {
     1208                // XXX fprintf (stderr, "keep\n");
     1209                kernels->solution1->data.F64[i] = solution->data.F64[i];
     1210                kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
     1211            }
     1212        }
     1213
    13041214
    13051215        memcpy(kernels->solution1->data.F64, solution->data.F64,
    … …  
    13641274    }
    13651275    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     1276    if (!isfinite(sum))  return false;
     1277    if (!isfinite(dmax)) return false;
     1278    if (!isfinite(dmin)) return false;
     1279    if (!isfinite(peak)) return false;
     1280
    13661281    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
    13671282    psVectorAppend(fSigRes, sigma/sum);
    … …  
    19731888}
    19741889
     1890
     1891# if 0
     1892
     1893#ifdef TESTING
     1894        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     1895        psVectorWriteFile ("B.dat", sumVector);
     1896#endif
     1897
     1898# define SVD_ANALYSIS 0
     1899# define COEFF_SIG 0.0
     1900# define SVD_TOL 0.0
     1901
     1902        // Use SVD to determine the kernel coeffs (and validate)
     1903        if (SVD_ANALYSIS) {
     1904
     1905            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
     1906            // sumMatrix * x = sumVector.
     1907
     1908            // we can use any standard matrix inversion to solve this.  However, the basis
     1909            // functions in general have substantial correlation, so that the solution may be
     1910            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
     1911            // system of equations may be statistically ill-conditioned.  Noise in the image
     1912            // will drive insignificant, but correlated, terms in the solution.  To avoid these
     1913            // problems, we can use SVD to identify numerically unconstrained values and to
     1914            // avoid statistically badly determined value.
     1915
     1916            // A = sumMatrix, B = sumVector
     1917            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
     1918            // x = V (1/w) (U^T B)
     1919            // \sigma_x = sqrt(diag(A^{-1}))
     1920            // solve for x and A^{-1} to get x & dx
     1921            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
     1922            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
     1923
     1924            // If I use the SVD trick to re-condition the matrix, I need to break out the
     1925            // kernel and normalization terms from the background term.
     1926            // XXX is this true?  or was this due to an error in the analysis?
     1927
     1928            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1929
     1930            // now pull out the kernel elements into their own square matrix
     1931            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
     1932            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
     1933
     1934            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
     1935                if (ix == bgIndex) continue;
     1936                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
     1937                    if (iy == bgIndex) continue;
     1938                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
     1939                    ky++;
     1940                }
     1941                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
     1942                kx++;
     1943            }
     1944
     1945            psImage *U = NULL;
     1946            psImage *V = NULL;
     1947            psVector *w = NULL;
     1948            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
     1949                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
     1950                return NULL;
     1951            }
     1952
     1953            // calculate A_inverse:
     1954            // Ainv = V * w * U^T
     1955            psImage *wUt  = p_pmSubSolve_wUt (w, U);
     1956            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
     1957            psImage *Xvar = NULL;
     1958            psFree (wUt);
     1959
     1960# ifdef TESTING
     1961            // kernel terms:
     1962            for (int i = 0; i < w->n; i++) {
     1963                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     1964            }
     1965# endif
     1966            // loop over w adding in more and more of the values until chisquare is no longer
     1967            // dropping significantly.
     1968            // XXX this does not seem to work very well: we seem to need all terms even for
     1969            // simple cases...
     1970
     1971            psVector *Xsvd = NULL;
     1972            {
     1973                psVector *Ax = NULL;
     1974                psVector *UtB = NULL;
     1975                psVector *wUtB = NULL;
     1976
     1977                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
     1978                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
     1979                psVectorInit (wMask, 1); // start by masking everything
     1980
     1981                double chiSquareLast = NAN;
     1982                int maxWeight = 0;
     1983
     1984                double Axx, Bx, y2;
     1985
     1986                // XXX this is an attempt to exclude insignificant modes.
     1987                // it was not successful with the ISIS kernel set: removing even
     1988                // the least significant mode leaves additional ringing / noise
     1989                // because the terms are so coupled.
     1990                for (int k = 0; false && (k < w->n); k++) {
     1991
     1992                    // unmask the k-th weight
     1993                    wMask->data.U8[k] = 0;
     1994                    p_pmSubSolve_SetWeights(wApply, w, wMask);
     1995
     1996                    // solve for x:
     1997                    // x = V * w * (U^T * B)
     1998                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
     1999                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     2000                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     2001
     2002                    // chi-square for this system of equations:
     2003                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     2004                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     2005                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     2006                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     2007                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     2008                    p_pmSubSolve_y2 (&y2, kernels, stamps);
     2009
     2010                    // apparently, this works (compare with the brute force value below
     2011                    double chiSquare = Axx - 2.0*Bx + y2;
     2012                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
     2013                    chiSquareLast = chiSquare;
     2014
     2015                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
     2016                    if (k && !maxWeight && (deltaChi < 1.0)) {
     2017                        maxWeight = k;
     2018                    }
     2019                }
     2020
     2021                // keep all terms or we get extra ringing
     2022                maxWeight = w->n;
     2023                psVectorInit (wMask, 1);
     2024                for (int k = 0; k < maxWeight; k++) {
     2025                    wMask->data.U8[k] = 0;
     2026                }
     2027                p_pmSubSolve_SetWeights(wApply, w, wMask);
     2028
     2029                // solve for x:
     2030                // x = V * w * (U^T * B)
     2031                p_pmSubSolve_UtB (&UtB, U, kernelVector);
     2032                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     2033                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     2034
     2035                // chi-square for this system of equations:
     2036                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     2037                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     2038                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     2039                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     2040                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     2041                p_pmSubSolve_y2 (&y2, kernels, stamps);
     2042
     2043                // apparently, this works (compare with the brute force value below
     2044                double chiSquare = Axx - 2.0*Bx + y2;
     2045                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
     2046
     2047                // re-calculate A^{-1} to get new variances:
     2048                // Ainv = V * w * U^T
     2049                // XXX since we keep all terms, this is identical to Ainv
     2050                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
     2051                Xvar = p_pmSubSolve_VwUt (V, wUt);
     2052                psFree (wUt);
     2053
     2054                psFree (Ax);
     2055                psFree (UtB);
     2056                psFree (wUtB);
     2057                psFree (wApply);
     2058                psFree (wMask);
     2059            }
     2060
     2061            // copy the kernel solutions to the full solution vector:
     2062            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2063            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2064
     2065            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
     2066                if (ix == bgIndex) {
     2067                    solution->data.F64[ix] = 0;
     2068                    solutionErr->data.F64[ix] = 0.001;
     2069                    continue;
     2070                }
     2071                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
     2072                solution->data.F64[ix] = Xsvd->data.F64[kx];
     2073                kx++;
     2074            }
     2075
     2076            psFree (kernelMatrix);
     2077            psFree (kernelVector);
     2078
     2079            psFree (U);
     2080            psFree (V);
     2081            psFree (w);
     2082
     2083            psFree (Ainv);
     2084            psFree (Xsvd);
     2085        } else {
     2086            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     2087            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
     2088            if (!luMatrix) {
     2089                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     2090                psFree(solution);
     2091                psFree(sumVector);
     2092                psFree(sumMatrix);
     2093                psFree(luMatrix);
     2094                psFree(permutation);
     2095                return NULL;
     2096            }
     2097
     2098            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     2099            psFree(luMatrix);
     2100            psFree(permutation);
     2101            if (!solution) {
     2102                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     2103                psFree(solution);
     2104                psFree(sumVector);
     2105                psFree(sumMatrix);
     2106                return NULL;
     2107            }
     2108
     2109            // XXX LUD does not provide A^{-1}?  fake the error for now
     2110            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2111            for (int ix = 0; ix < sumVector->n; ix++) {
     2112                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
     2113            }
     2114        }
     2115
     2116        if (!kernels->solution1) {
     2117            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     2118            psVectorInit (kernels->solution1, 0.0);
     2119        }
     2120
     2121        // only update the solutions that we chose to calculate:
     2122        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     2123            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     2124            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     2125        }
     2126        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     2127            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     2128            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     2129        }
     2130        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     2131            int numKernels = kernels->num;
     2132            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     2133            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     2134            for (int i = 0; i < numKernels * numPoly; i++) {
     2135                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     2136                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
     2137                    // XXX fprintf (stderr, "drop\n");
     2138                    kernels->solution1->data.F64[i] = 0.0;
     2139                } else {
     2140                    // XXX fprintf (stderr, "keep\n");
     2141                    kernels->solution1->data.F64[i] = solution->data.F64[i];
     2142                }
     2143            }
     2144        }
     2145        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
     2146        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
     2147
     2148        psFree(solution);
     2149        psFree(sumVector);
     2150        psFree(sumMatrix);
     2151# endif
     2152
     2153#ifdef TESTING
     2154              // XXX double-check for NAN in data:
     2155                for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
     2156                    for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
     2157                        if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
     2158                            fprintf (stderr, "WARNING: NAN in matrix\n");
     2159                        }
     2160                    }
     2161                }
     2162                for (int ix = 0; ix < stamp->vector->n; ix++) {
     2163                    if (!isfinite(stamp->vector->data.F64[ix])) {
     2164                        fprintf (stderr, "WARNING: NAN in vector\n");
     2165                    }
     2166                }
     2167#endif
     2168
     2169#ifdef TESTING
     2170        for (int ix = 0; ix < sumVector->n; ix++) {
     2171            if (!isfinite(sumVector->data.F64[ix])) {
     2172                fprintf (stderr, "WARNING: NAN in vector\n");
     2173            }
     2174        }
     2175#endif
     2176
     2177#ifdef TESTING
     2178        for (int ix = 0; ix < sumVector->n; ix++) {
     2179            if (!isfinite(sumVector->data.F64[ix])) {
     2180                fprintf (stderr, "WARNING: NAN in vector\n");
     2181            }
     2182        }
     2183        {
     2184            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
     2185            psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
     2186            psFree(inverse);
     2187        }
     2188        {
     2189            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
     2190            psImage *Xt = psMatrixTranspose(NULL, X);
     2191            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
     2192            psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
     2193            psFree(X);
     2194            psFree(Xt);
     2195            psFree(XtX);
     2196        }
     2197#endif
     2198
Note: See TracChangeset for help on using the changeset viewer.