IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26737


Ignore:
Timestamp:
Jan 29, 2010, 11:05:32 AM (16 years ago)
Author:
eugene
Message:

add normalization into the SINGLE as well as the DUAL convolution; catch NANs in the system equation elements for each stamp; use SVD for SINGLE; move a bunch of testing code elsewhere; change penalty function to use the scaled-normalization not the background

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

    r26703 r26737  
    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                                  )
    … …  
    172174    double sumR = 0.0;                  // Sum of the reference
    173175    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
    174177    for (int y = - footprint; y <= footprint; y++) {
    175178        for (int x = - footprint; x <= footprint; x++) {
    … …  
    179182            double rr = PS_SQR(ref);
    180183            double one = 1.0;
     184
     185            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
     186                normI1 += ref;
     187                normI2 += in;
     188            }
     189
    181190            if (weight) {
    182191                float wtVal = weight->kernel[y][x];
    … …  
    202211        }
    203212    }
     213
     214    *norm = normI2 / normI1;
     215    fprintf(stderr, "Sums: %f %f %f, normWindow: %d\n", normI1, normI2, *norm, normWindow);
     216
    204217    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    205218        matrix->data.F64[normIndex][normIndex] = sumRR;
    … …  
    215228        matrix->data.F64[bgIndex][normIndex] = sumR;
    216229    }
     230
     231    // check for any NAN values in the result, skip if found:
     232    for (int iy = 0; iy < matrix->numRows; iy++) {
     233        for (int ix = 0; ix < matrix->numCols; ix++) {
     234            if (!isfinite(matrix->data.F64[iy][ix])) {
     235                fprintf (stderr, "WARNING: NAN in matrix\n");
     236                return false;
     237            }
     238        }
     239    }
     240    for (int ix = 0; ix < vector->n; ix++) {
     241        if (!isfinite(vector->data.F64[ix])) {
     242            fprintf (stderr, "WARNING: NAN in vector\n");
     243            return false;
     244        }
     245    }
     246
    217247    return true;
    218248}
    … …  
    492522
    493523    *norm = normI2 / normI1;
    494 
    495     fprintf(stderr, "Sums: %f %f %f\n", normI1, normI2, *norm);
     524    fprintf(stderr, "Sums: %f %f %f - I1I1: %f\n", normI1, normI2, *norm, sumI1I1);
    496525
    497526    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    … …  
    507536        matrix->data.F64[normIndex][bgIndex] = sumI1;
    508537    }
     538
     539    // check for any NAN values in the result, skip if found:
     540    for (int iy = 0; iy < matrix->numRows; iy++) {
     541        for (int ix = 0; ix < matrix->numCols; ix++) {
     542            if (!isfinite(matrix->data.F64[iy][ix])) {
     543                fprintf (stderr, "WARNING: NAN in matrix\n");
     544                return false;
     545            }
     546        }
     547    }
     548    for (int ix = 0; ix < vector->n; ix++) {
     549        if (!isfinite(vector->data.F64[ix])) {
     550            fprintf (stderr, "WARNING: NAN in vector\n");
     551            return false;
     552        }
     553    }
     554
     555
    509556    return true;
    510557}
    … …  
    512559#if 1
    513560// Add in penalty term to least-squares vector
    514 static bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
     561bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
    515562                             psVector *vector,                    // Vector to which to add in penalty term
    516563                             const pmSubtractionKernels *kernels, // Kernel parameters
    … …  
    527574    int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations
    528575    int numParams = numKernels * numSpatial;                 // Number of kernel parameters
    529     for (int i = 0; i < numParams; i++) {
     576
     577    // order is :
     578    // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0]
     579    // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0]
     580    // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1]
     581    // [norm]
     582    // [bg]
     583    // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0]
     584    // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0]
     585    // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1]
     586
     587    for (int i = 0; i < numKernels; i++) {
    530588        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
    531589            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
    532590                // Contribution to chi^2: a_i^2 P_i
    533591                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
     592                fprintf (stderr, "penalty main: %d %e vs %e x %f : %f\n",
     593                         index,
     594                         matrix->data.F64[index][index],
     595                         norm, penalties->data.F32[i],
     596                         (norm * penalties->data.F32[i]) / matrix->data.F64[index][index] );
    534597                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
    535598                if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     599                    fprintf (stderr, "penalty dual: %d %e vs %e x %f : %f\n",
     600                             index + numParams + 2,
     601                             matrix->data.F64[index + numParams + 2][index + numParams + 2],
     602                             norm, penalties->data.F32[i],
     603                             (norm * penalties->data.F32[i]) / matrix->data.F64[index + numParams + 2][index + numParams + 2] );
    536604                    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
     605                    // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
     606                    // penalties scale with second moments
     607                    //
    537608                }
    538609            }
    … …  
    716787    switch (kernels->mode) {
    717788      case PM_SUBTRACTION_MODE_1:
    718         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
     789        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1,
    719790                                       weight, window, stamp->convolutions1, kernels,
    720                                        polyValues, footprint, mode);
     791                                       polyValues, footprint, stamps->normWindow, mode);
    721792        break;
    722793      case PM_SUBTRACTION_MODE_2:
    723         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
     794        status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2,
    724795                                       weight, window, stamp->convolutions2, kernels,
    725                                        polyValues, footprint, mode);
     796                                       polyValues, footprint, stamps->normWindow, mode);
    726797        break;
    727798      case PM_SUBTRACTION_MODE_DUAL:
    … …  
    893964        psVectorInit(sumVector, 0.0);
    894965        psImageInit(sumMatrix, 0.0);
    895         int numStamps = 0;              // Number of good stamps
    896         for (int i = 0; i < stamps->num; i++) {
    897             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    898 
    899             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
    900 
    901 #ifdef TESTING
    902               // XXX double-check for NAN in data:
    903                 for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
    904                     for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
    905                         if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
    906                             fprintf (stderr, "WARNING: NAN in matrix\n");
    907                         }
    908                     }
    909                 }
    910                 for (int ix = 0; ix < stamp->vector->n; ix++) {
    911                     if (!isfinite(stamp->vector->data.F64[ix])) {
    912                         fprintf (stderr, "WARNING: NAN in vector\n");
    913                     }
    914                 }
    915 #endif
    916 
    917                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    918                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
    919                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    920                 numStamps++;
    921             } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
    922                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
    923             }
    924         }
    925 
    926 #ifdef TESTING
    927         for (int ix = 0; ix < sumVector->n; ix++) {
    928             if (!isfinite(sumVector->data.F64[ix])) {
    929                 fprintf (stderr, "WARNING: NAN in vector\n");
    930             }
    931         }
    932 #endif
    933 
    934 #if 0
    935         int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    936         calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
    937 #endif
    938 
    939 #ifdef TESTING
    940         for (int ix = 0; ix < sumVector->n; ix++) {
    941             if (!isfinite(sumVector->data.F64[ix])) {
    942                 fprintf (stderr, "WARNING: NAN in vector\n");
    943             }
    944         }
    945         {
    946             psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
    947             psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
    948             psFree(inverse);
    949         }
    950         {
    951             psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
    952             psImage *Xt = psMatrixTranspose(NULL, X);
    953             psImage *XtX = psMatrixMultiply(NULL, Xt, X);
    954             psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
    955             psFree(X);
    956             psFree(Xt);
    957             psFree(XtX);
    958         }
    959 #endif
    960 
    961 # define SVD_ANALYSIS 0
    962 # define COEFF_SIG 0.0
    963 # define SVD_TOL 0.0
    964 
    965         psVector *solution = NULL;
    966         psVector *solutionErr = NULL;
    967 
    968 #ifdef TESTING
    969         psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
    970         psVectorWriteFile ("B.dat", sumVector);
    971 #endif
    972 
    973         // Use SVD to determine the kernel coeffs (and validate)
    974         if (SVD_ANALYSIS) {
    975 
    976             // We have sumVector and sumMatrix.  we are trying to solve the following equation:
    977             // sumMatrix * x = sumVector.
    978 
    979             // we can use any standard matrix inversion to solve this.  However, the basis
    980             // functions in general have substantial correlation, so that the solution may be
    981             // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
    982             // system of equations may be statistically ill-conditioned.  Noise in the image
    983             // will drive insignificant, but correlated, terms in the solution.  To avoid these
    984             // problems, we can use SVD to identify numerically unconstrained values and to
    985             // avoid statistically badly determined value.
    986 
    987             // A = sumMatrix, B = sumVector
    988             // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
    989             // x = V (1/w) (U^T B)
    990             // \sigma_x = sqrt(diag(A^{-1}))
    991             // solve for x and A^{-1} to get x & dx
    992             // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
    993             // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
    994 
    995             // If I use the SVD trick to re-condition the matrix, I need to break out the
    996             // kernel and normalization terms from the background term.
    997             // XXX is this true?  or was this due to an error in the analysis?
    998 
    999             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1000 
    1001             // now pull out the kernel elements into their own square matrix
    1002             psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
    1003             psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
    1004 
    1005             for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
    1006                 if (ix == bgIndex) continue;
    1007                 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
    1008                     if (iy == bgIndex) continue;
    1009                     kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
    1010                     ky++;
    1011                 }
    1012                 kernelVector->data.F64[kx] = sumVector->data.F64[ix];
    1013                 kx++;
    1014             }
    1015 
    1016             psImage *U = NULL;
    1017             psImage *V = NULL;
    1018             psVector *w = NULL;
    1019             if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
    1020                 psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
    1021                 return NULL;
    1022             }
    1023 
    1024             // calculate A_inverse:
    1025             // Ainv = V * w * U^T
    1026             psImage *wUt  = p_pmSubSolve_wUt (w, U);
    1027             psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
    1028             psImage *Xvar = NULL;
    1029             psFree (wUt);
    1030 
    1031 # ifdef TESTING
    1032             // kernel terms:
    1033             for (int i = 0; i < w->n; i++) {
    1034                 fprintf (stderr, "w: %f\n", w->data.F64[i]);
    1035             }
    1036 # endif
    1037             // loop over w adding in more and more of the values until chisquare is no longer
    1038             // dropping significantly.
    1039             // XXX this does not seem to work very well: we seem to need all terms even for
    1040             // simple cases...
    1041 
    1042             psVector *Xsvd = NULL;
    1043             {
    1044                 psVector *Ax = NULL;
    1045                 psVector *UtB = NULL;
    1046                 psVector *wUtB = NULL;
    1047 
    1048                 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
    1049                 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
    1050                 psVectorInit (wMask, 1); // start by masking everything
    1051 
    1052                 double chiSquareLast = NAN;
    1053                 int maxWeight = 0;
    1054 
    1055                 double Axx, Bx, y2;
    1056 
    1057                 // XXX this is an attempt to exclude insignificant modes.
    1058                 // it was not successful with the ISIS kernel set: removing even
    1059                 // the least significant mode leaves additional ringing / noise
    1060                 // because the terms are so coupled.
    1061                 for (int k = 0; false && (k < w->n); k++) {
    1062 
    1063                     // unmask the k-th weight
    1064                     wMask->data.U8[k] = 0;
    1065                     p_pmSubSolve_SetWeights(wApply, w, wMask);
    1066 
    1067                     // solve for x:
    1068                     // x = V * w * (U^T * B)
    1069                     p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1070                     p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1071                     p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1072 
    1073                     // chi-square for this system of equations:
    1074                     // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1075                     // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1076                     p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1077                     p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1078                     p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1079                     p_pmSubSolve_y2 (&y2, kernels, stamps);
    1080 
    1081                     // apparently, this works (compare with the brute force value below
    1082                     double chiSquare = Axx - 2.0*Bx + y2;
    1083                     double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
    1084                     chiSquareLast = chiSquare;
    1085 
    1086                     // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
    1087                     if (k && !maxWeight && (deltaChi < 1.0)) {
    1088                         maxWeight = k;
    1089                     }
    1090                 }
    1091 
    1092                 // keep all terms or we get extra ringing
    1093                 maxWeight = w->n;
    1094                 psVectorInit (wMask, 1);
    1095                 for (int k = 0; k < maxWeight; k++) {
    1096                     wMask->data.U8[k] = 0;
    1097                 }
    1098                 p_pmSubSolve_SetWeights(wApply, w, wMask);
    1099 
    1100                 // solve for x:
    1101                 // x = V * w * (U^T * B)
    1102                 p_pmSubSolve_UtB (&UtB, U, kernelVector);
    1103                 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
    1104                 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
    1105 
    1106                 // chi-square for this system of equations:
    1107                 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
    1108                 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
    1109                 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
    1110                 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
    1111                 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
    1112                 p_pmSubSolve_y2 (&y2, kernels, stamps);
    1113 
    1114                 // apparently, this works (compare with the brute force value below
    1115                 double chiSquare = Axx - 2.0*Bx + y2;
    1116                 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
    1117 
    1118                 // re-calculate A^{-1} to get new variances:
    1119                 // Ainv = V * w * U^T
    1120                 // XXX since we keep all terms, this is identical to Ainv
    1121                 psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
    1122                 Xvar = p_pmSubSolve_VwUt (V, wUt);
    1123                 psFree (wUt);
    1124 
    1125                 psFree (Ax);
    1126                 psFree (UtB);
    1127                 psFree (wUtB);
    1128                 psFree (wApply);
    1129                 psFree (wMask);
    1130             }
    1131 
    1132             // copy the kernel solutions to the full solution vector:
    1133             solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1134             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1135 
    1136             for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
    1137                 if (ix == bgIndex) {
    1138                     solution->data.F64[ix] = 0;
    1139                     solutionErr->data.F64[ix] = 0.001;
    1140                     continue;
    1141                 }
    1142                 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
    1143                 solution->data.F64[ix] = Xsvd->data.F64[kx];
    1144                 kx++;
    1145             }
    1146 
    1147             psFree (kernelMatrix);
    1148             psFree (kernelVector);
    1149 
    1150             psFree (U);
    1151             psFree (V);
    1152             psFree (w);
    1153 
    1154             psFree (Ainv);
    1155             psFree (Xsvd);
    1156         } else {
    1157             psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
    1158             psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
    1159             if (!luMatrix) {
    1160                 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
    1161                 psFree(solution);
    1162                 psFree(sumVector);
    1163                 psFree(sumMatrix);
    1164                 psFree(luMatrix);
    1165                 psFree(permutation);
    1166                 return NULL;
    1167             }
    1168 
    1169             solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
    1170             psFree(luMatrix);
    1171             psFree(permutation);
    1172             if (!solution) {
    1173                 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
    1174                 psFree(solution);
    1175                 psFree(sumVector);
    1176                 psFree(sumMatrix);
    1177                 return NULL;
    1178             }
    1179 
    1180             // XXX LUD does not provide A^{-1}?  fake the error for now
    1181             solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
    1182             for (int ix = 0; ix < sumVector->n; ix++) {
    1183                 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
    1184             }
    1185         }
    1186 
    1187         if (!kernels->solution1) {
    1188             kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
    1189             psVectorInit (kernels->solution1, 0.0);
    1190         }
    1191 
    1192         // only update the solutions that we chose to calculate:
    1193         if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    1194             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1195             kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
    1196         }
    1197         if (mode & PM_SUBTRACTION_EQUATION_BG) {
    1198             int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
    1199             kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
    1200         }
    1201         if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
    1202             int numKernels = kernels->num;
    1203             int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
    1204             int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
    1205             for (int i = 0; i < numKernels * numPoly; i++) {
    1206                 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
    1207                 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
    1208                     // XXX fprintf (stderr, "drop\n");
    1209                     kernels->solution1->data.F64[i] = 0.0;
    1210                 } else {
    1211                     // XXX fprintf (stderr, "keep\n");
    1212                     kernels->solution1->data.F64[i] = solution->data.F64[i];
    1213                 }
    1214             }
    1215         }
    1216         // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
    1217         // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
    1218 
    1219         psFree(solution);
    1220         psFree(sumVector);
    1221         psFree(sumMatrix);
    1222 
    1223 #ifdef TESTING
    1224         // XXX double-check for NAN in data:
    1225         for (int ix = 0; ix < kernels->solution1->n; ix++) {
    1226             if (!isfinite(kernels->solution1->data.F64[ix])) {
    1227                 fprintf (stderr, "WARNING: NAN in vector\n");
    1228             }
    1229         }
    1230 #endif
    1231 
    1232     } else {
    1233         // Dual convolution solution
    1234 
    1235         // Accumulation of stamp matrices/vectors
    1236         psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    1237         psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
    1238         psImageInit(sumMatrix, 0.0);
    1239         psVectorInit(sumVector, 0.0);
    1240966
    1241967        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
    … …  
    1247973                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
    1248974                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
    1249 
    1250975                psVectorAppend(norms, stamp->norm);
    1251 
    1252976                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    1253977                numStamps++;
    1254             }
    1255         }
    1256 
    1257 #ifdef TESTING
    1258         psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
    1259         psVectorWriteFile("sumVector.dat", sumVector);
    1260 #endif
    1261 
    1262 #if 1
     978            } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
     979                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
     980            }
     981        }
     982
     983#if 0
    1263984        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
    1264985        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
    … …  
    1276997            // Solve normalisation and background separately
    1277998            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    1278 
    1279 #if 0
    1280             psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
    1281             psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
    1282 
    1283             normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
    1284             normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
    1285             normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
    1286 
    1287             normVector->data.F64[0] = sumVector->data.F64[normIndex];
    1288             normVector->data.F64[1] = sumVector->data.F64[bgIndex];
    1289 
    1290             psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
    1291 
    1292             double normValue = normSolution->data.F64[0];
    1293             double bgValue = normSolution->data.F64[1];
    1294 
    1295             psFree(normMatrix);
    1296             psFree(normVector);
    1297             psFree(normSolution);
    1298 #endif
    1299999
    13001000            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
    … …  
    13091009
    13101010            double normValue = stats->robustMedian;
    1311 //            double bgValue = 0.0;
     1011            // double bgValue = 0.0;
    13121012
    13131013            psFree(stats);
    13141014
    13151015            fprintf(stderr, "Norm: %lf\n", normValue);
    1316 //            fprintf(stderr, "BG: %lf\n", bgValue);
     1016            // fprintf(stderr, "BG: %lf\n", bgValue);
     1017
     1018            // Solve kernel components
     1019            for (int i = 0; i < numSolution1; i++) {
     1020                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
     1021
     1022                sumMatrix->data.F64[i][normIndex] = 0.0;
     1023                sumMatrix->data.F64[normIndex][i] = 0.0;
     1024
     1025                // sumMatrix->data.F64[i][bgIndex] = 0.0;
     1026                // sumMatrix->data.F64[bgIndex][i] = 0.0;
     1027            }
     1028
     1029            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
     1030            // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
     1031            // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1032
     1033            sumVector->data.F64[normIndex] = 0.0;
     1034            // sumVector->data.F64[bgIndex] = 0.0;
     1035
     1036            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1037
     1038            solution->data.F64[normIndex] = normValue;
     1039            // solution->data.F64[bgIndex] = bgValue;
     1040        }
     1041# endif
     1042
     1043        if (!kernels->solution1) {
     1044            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     1045            psVectorInit (kernels->solution1, 0.0);
     1046        }
     1047
     1048        // only update the solutions that we chose to calculate:
     1049        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     1050            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1051            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     1052        }
     1053        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     1054            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1055            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     1056        }
     1057        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     1058            int numKernels = kernels->num;
     1059            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     1060            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     1061            for (int i = 0; i < numKernels * numPoly; i++) {
     1062                kernels->solution1->data.F64[i] = solution->data.F64[i];
     1063            }
     1064        }
     1065
     1066        psFree(solution);
     1067        psFree(sumVector);
     1068        psFree(sumMatrix);
     1069
     1070#ifdef TESTING
     1071        // XXX double-check for NAN in data:
     1072        for (int ix = 0; ix < kernels->solution1->n; ix++) {
     1073            if (!isfinite(kernels->solution1->data.F64[ix])) {
     1074                fprintf (stderr, "WARNING: NAN in vector\n");
     1075            }
     1076        }
     1077#endif
     1078
     1079    } else {
     1080        // Dual convolution solution
     1081
     1082        // Accumulation of stamp matrices/vectors
     1083        psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
     1084        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
     1085        psImageInit(sumMatrix, 0.0);
     1086        psVectorInit(sumVector, 0.0);
     1087
     1088        psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations
     1089
     1090        int numStamps = 0;              // Number of good stamps
     1091        for (int i = 0; i < stamps->num; i++) {
     1092            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     1093            if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
     1094                (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
     1095                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
     1096
     1097                psVectorAppend(norms, stamp->norm);
     1098
     1099                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
     1100                numStamps++;
     1101            }
     1102        }
     1103
     1104#ifdef TESTING
     1105        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
     1106        psVectorWriteFile("sumVector.dat", sumVector);
     1107#endif
     1108
     1109#if 1
     1110        // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
     1111        // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
     1112
     1113        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1114        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0);
     1115#endif
     1116
     1117        psVector *solution = NULL;                       // Solution to equation!
     1118        solution = psVectorAlloc(numParams, PS_TYPE_F64);
     1119        psVectorInit(solution, 0);
     1120
     1121#if 0
     1122        // Regular, straight-forward solution
     1123        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1124#else
     1125        {
     1126            // Solve normalisation and background separately
     1127            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     1128
     1129#if 0
     1130            psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64);
     1131            psVector *normVector = psVectorAlloc(2, PS_TYPE_F64);
     1132
     1133            normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex];
     1134            normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex];
     1135            normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex];
     1136
     1137            normVector->data.F64[0] = sumVector->data.F64[normIndex];
     1138            normVector->data.F64[1] = sumVector->data.F64[bgIndex];
     1139
     1140            psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN);
     1141
     1142            double normValue = normSolution->data.F64[0];
     1143            double bgValue = normSolution->data.F64[1];
     1144
     1145            psFree(normMatrix);
     1146            psFree(normVector);
     1147            psFree(normSolution);
     1148#endif
     1149
     1150            psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
     1151            if (!psVectorStats(stats, norms, NULL, NULL, 0)) {
     1152                psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation");
     1153                psFree(stats);
     1154                psFree(sumMatrix);
     1155                psFree(sumVector);
     1156                psFree(norms);
     1157                return false;
     1158            }
     1159
     1160            double normValue = stats->robustMedian;
     1161            // double bgValue = 0.0;
     1162
     1163            psFree(stats);
     1164
     1165            fprintf(stderr, "Norm: %lf\n", normValue);
     1166            // fprintf(stderr, "BG: %lf\n", bgValue);
    13171167
    13181168            // Solve kernel components
    … …  
    13241174                sumMatrix->data.F64[normIndex][i] = 0.0;
    13251175
    1326 //                sumMatrix->data.F64[i][bgIndex] = 0.0;
    1327 //                sumMatrix->data.F64[bgIndex][i] = 0.0;
     1176                // sumMatrix->data.F64[i][bgIndex] = 0.0;
     1177                // sumMatrix->data.F64[bgIndex][i] = 0.0;
    13281178                sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
    13291179                sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
    1330 //                sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0;
    1331 //                sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0;
     1180                // sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0;
     1181                // sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0;
    13321182            }
    13331183
    13341184            sumMatrix->data.F64[normIndex][normIndex] = 1.0;
    1335 //            sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
    1336 //            sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
     1185            // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
     1186            // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
    13371187
    13381188            sumVector->data.F64[normIndex] = 0.0;
    1339 //            sumVector->data.F64[bgIndex] = 0.0;
     1189            // sumVector->data.F64[bgIndex] = 0.0;
    13401190
    13411191            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    13421192
    13431193            solution->data.F64[normIndex] = normValue;
    1344 //            solution->data.F64[bgIndex] = bgValue;
     1194            // solution->data.F64[bgIndex] = bgValue;
    13451195        }
    13461196#endif
    … …  
    14471297    }
    14481298    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     1299    if (!isfinite(sum))  return false;
     1300    if (!isfinite(dmax)) return false;
     1301    if (!isfinite(dmin)) return false;
     1302    if (!isfinite(peak)) return false;
     1303
    14491304    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
    14501305    psVectorAppend(fSigRes, sigma/sum);
    … …  
    20561911}
    20571912
     1913
     1914# if 0
     1915
     1916#ifdef TESTING
     1917        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
     1918        psVectorWriteFile ("B.dat", sumVector);
     1919#endif
     1920
     1921# define SVD_ANALYSIS 0
     1922# define COEFF_SIG 0.0
     1923# define SVD_TOL 0.0
     1924
     1925        // Use SVD to determine the kernel coeffs (and validate)
     1926        if (SVD_ANALYSIS) {
     1927
     1928            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
     1929            // sumMatrix * x = sumVector.
     1930
     1931            // we can use any standard matrix inversion to solve this.  However, the basis
     1932            // functions in general have substantial correlation, so that the solution may be
     1933            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
     1934            // system of equations may be statistically ill-conditioned.  Noise in the image
     1935            // will drive insignificant, but correlated, terms in the solution.  To avoid these
     1936            // problems, we can use SVD to identify numerically unconstrained values and to
     1937            // avoid statistically badly determined value.
     1938
     1939            // A = sumMatrix, B = sumVector
     1940            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
     1941            // x = V (1/w) (U^T B)
     1942            // \sigma_x = sqrt(diag(A^{-1}))
     1943            // solve for x and A^{-1} to get x & dx
     1944            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
     1945            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
     1946
     1947            // If I use the SVD trick to re-condition the matrix, I need to break out the
     1948            // kernel and normalization terms from the background term.
     1949            // XXX is this true?  or was this due to an error in the analysis?
     1950
     1951            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     1952
     1953            // now pull out the kernel elements into their own square matrix
     1954            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
     1955            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
     1956
     1957            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
     1958                if (ix == bgIndex) continue;
     1959                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
     1960                    if (iy == bgIndex) continue;
     1961                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
     1962                    ky++;
     1963                }
     1964                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
     1965                kx++;
     1966            }
     1967
     1968            psImage *U = NULL;
     1969            psImage *V = NULL;
     1970            psVector *w = NULL;
     1971            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
     1972                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
     1973                return NULL;
     1974            }
     1975
     1976            // calculate A_inverse:
     1977            // Ainv = V * w * U^T
     1978            psImage *wUt  = p_pmSubSolve_wUt (w, U);
     1979            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
     1980            psImage *Xvar = NULL;
     1981            psFree (wUt);
     1982
     1983# ifdef TESTING
     1984            // kernel terms:
     1985            for (int i = 0; i < w->n; i++) {
     1986                fprintf (stderr, "w: %f\n", w->data.F64[i]);
     1987            }
     1988# endif
     1989            // loop over w adding in more and more of the values until chisquare is no longer
     1990            // dropping significantly.
     1991            // XXX this does not seem to work very well: we seem to need all terms even for
     1992            // simple cases...
     1993
     1994            psVector *Xsvd = NULL;
     1995            {
     1996                psVector *Ax = NULL;
     1997                psVector *UtB = NULL;
     1998                psVector *wUtB = NULL;
     1999
     2000                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
     2001                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
     2002                psVectorInit (wMask, 1); // start by masking everything
     2003
     2004                double chiSquareLast = NAN;
     2005                int maxWeight = 0;
     2006
     2007                double Axx, Bx, y2;
     2008
     2009                // XXX this is an attempt to exclude insignificant modes.
     2010                // it was not successful with the ISIS kernel set: removing even
     2011                // the least significant mode leaves additional ringing / noise
     2012                // because the terms are so coupled.
     2013                for (int k = 0; false && (k < w->n); k++) {
     2014
     2015                    // unmask the k-th weight
     2016                    wMask->data.U8[k] = 0;
     2017                    p_pmSubSolve_SetWeights(wApply, w, wMask);
     2018
     2019                    // solve for x:
     2020                    // x = V * w * (U^T * B)
     2021                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
     2022                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     2023                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     2024
     2025                    // chi-square for this system of equations:
     2026                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     2027                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     2028                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     2029                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     2030                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     2031                    p_pmSubSolve_y2 (&y2, kernels, stamps);
     2032
     2033                    // apparently, this works (compare with the brute force value below
     2034                    double chiSquare = Axx - 2.0*Bx + y2;
     2035                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
     2036                    chiSquareLast = chiSquare;
     2037
     2038                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
     2039                    if (k && !maxWeight && (deltaChi < 1.0)) {
     2040                        maxWeight = k;
     2041                    }
     2042                }
     2043
     2044                // keep all terms or we get extra ringing
     2045                maxWeight = w->n;
     2046                psVectorInit (wMask, 1);
     2047                for (int k = 0; k < maxWeight; k++) {
     2048                    wMask->data.U8[k] = 0;
     2049                }
     2050                p_pmSubSolve_SetWeights(wApply, w, wMask);
     2051
     2052                // solve for x:
     2053                // x = V * w * (U^T * B)
     2054                p_pmSubSolve_UtB (&UtB, U, kernelVector);
     2055                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
     2056                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
     2057
     2058                // chi-square for this system of equations:
     2059                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
     2060                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
     2061                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
     2062                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
     2063                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
     2064                p_pmSubSolve_y2 (&y2, kernels, stamps);
     2065
     2066                // apparently, this works (compare with the brute force value below
     2067                double chiSquare = Axx - 2.0*Bx + y2;
     2068                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
     2069
     2070                // re-calculate A^{-1} to get new variances:
     2071                // Ainv = V * w * U^T
     2072                // XXX since we keep all terms, this is identical to Ainv
     2073                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
     2074                Xvar = p_pmSubSolve_VwUt (V, wUt);
     2075                psFree (wUt);
     2076
     2077                psFree (Ax);
     2078                psFree (UtB);
     2079                psFree (wUtB);
     2080                psFree (wApply);
     2081                psFree (wMask);
     2082            }
     2083
     2084            // copy the kernel solutions to the full solution vector:
     2085            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2086            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2087
     2088            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
     2089                if (ix == bgIndex) {
     2090                    solution->data.F64[ix] = 0;
     2091                    solutionErr->data.F64[ix] = 0.001;
     2092                    continue;
     2093                }
     2094                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
     2095                solution->data.F64[ix] = Xsvd->data.F64[kx];
     2096                kx++;
     2097            }
     2098
     2099            psFree (kernelMatrix);
     2100            psFree (kernelVector);
     2101
     2102            psFree (U);
     2103            psFree (V);
     2104            psFree (w);
     2105
     2106            psFree (Ainv);
     2107            psFree (Xsvd);
     2108        } else {
     2109            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
     2110            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
     2111            if (!luMatrix) {
     2112                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
     2113                psFree(solution);
     2114                psFree(sumVector);
     2115                psFree(sumMatrix);
     2116                psFree(luMatrix);
     2117                psFree(permutation);
     2118                return NULL;
     2119            }
     2120
     2121            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
     2122            psFree(luMatrix);
     2123            psFree(permutation);
     2124            if (!solution) {
     2125                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
     2126                psFree(solution);
     2127                psFree(sumVector);
     2128                psFree(sumMatrix);
     2129                return NULL;
     2130            }
     2131
     2132            // XXX LUD does not provide A^{-1}?  fake the error for now
     2133            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
     2134            for (int ix = 0; ix < sumVector->n; ix++) {
     2135                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
     2136            }
     2137        }
     2138
     2139        if (!kernels->solution1) {
     2140            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
     2141            psVectorInit (kernels->solution1, 0.0);
     2142        }
     2143
     2144        // only update the solutions that we chose to calculate:
     2145        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
     2146            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
     2147            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
     2148        }
     2149        if (mode & PM_SUBTRACTION_EQUATION_BG) {
     2150            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
     2151            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
     2152        }
     2153        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
     2154            int numKernels = kernels->num;
     2155            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
     2156            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
     2157            for (int i = 0; i < numKernels * numPoly; i++) {
     2158                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
     2159                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
     2160                    // XXX fprintf (stderr, "drop\n");
     2161                    kernels->solution1->data.F64[i] = 0.0;
     2162                } else {
     2163                    // XXX fprintf (stderr, "keep\n");
     2164                    kernels->solution1->data.F64[i] = solution->data.F64[i];
     2165                }
     2166            }
     2167        }
     2168        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
     2169        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
     2170
     2171        psFree(solution);
     2172        psFree(sumVector);
     2173        psFree(sumMatrix);
     2174# endif
     2175
     2176#ifdef TESTING
     2177              // XXX double-check for NAN in data:
     2178                for (int iy = 0; iy < stamp->matrix->numRows; iy++) {
     2179                    for (int ix = 0; ix < stamp->matrix->numCols; ix++) {
     2180                        if (!isfinite(stamp->matrix->data.F64[iy][ix])) {
     2181                            fprintf (stderr, "WARNING: NAN in matrix\n");
     2182                        }
     2183                    }
     2184                }
     2185                for (int ix = 0; ix < stamp->vector->n; ix++) {
     2186                    if (!isfinite(stamp->vector->data.F64[ix])) {
     2187                        fprintf (stderr, "WARNING: NAN in vector\n");
     2188                    }
     2189                }
     2190#endif
     2191
     2192#ifdef TESTING
     2193        for (int ix = 0; ix < sumVector->n; ix++) {
     2194            if (!isfinite(sumVector->data.F64[ix])) {
     2195                fprintf (stderr, "WARNING: NAN in vector\n");
     2196            }
     2197        }
     2198#endif
     2199
     2200#ifdef TESTING
     2201        for (int ix = 0; ix < sumVector->n; ix++) {
     2202            if (!isfinite(sumVector->data.F64[ix])) {
     2203                fprintf (stderr, "WARNING: NAN in vector\n");
     2204            }
     2205        }
     2206        {
     2207            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
     2208            psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);
     2209            psFree(inverse);
     2210        }
     2211        {
     2212            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
     2213            psImage *Xt = psMatrixTranspose(NULL, X);
     2214            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
     2215            psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);
     2216            psFree(X);
     2217            psFree(Xt);
     2218            psFree(XtX);
     2219        }
     2220#endif
     2221
Note: See TracChangeset for help on using the changeset viewer.