Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26735) +++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (revision 26737) @@ -28,4 +28,5 @@ static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated psVector *vector, // Least-squares vector, updated + double *norm, // Normalisation, updated const psKernel *input, // Input image (target) const psKernel *reference, // Reference image (convolution source) @@ -36,4 +37,5 @@ const psImage *polyValues, // Spatial polynomial values int footprint, // (Half-)Size of stamp + int normWindow, // Window (half-)size for normalisation measurement const pmSubtractionEquationCalculationMode mode ) @@ -172,4 +174,5 @@ double sumR = 0.0; // Sum of the reference double sumI = 0.0; // Sum of the input + double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -179,4 +182,10 @@ double rr = PS_SQR(ref); double one = 1.0; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { + normI1 += ref; + normI2 += in; + } + if (weight) { float wtVal = weight->kernel[y][x]; @@ -202,4 +211,8 @@ } } + + *norm = normI2 / normI1; + fprintf(stderr, "Sums: %f %f %f, normWindow: %d\n", normI1, normI2, *norm, normWindow); + if (mode & PM_SUBTRACTION_EQUATION_NORM) { matrix->data.F64[normIndex][normIndex] = sumRR; @@ -215,4 +228,21 @@ matrix->data.F64[bgIndex][normIndex] = sumR; } + + // check for any NAN values in the result, skip if found: + for (int iy = 0; iy < matrix->numRows; iy++) { + for (int ix = 0; ix < matrix->numCols; ix++) { + if (!isfinite(matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + return false; + } + } + } + for (int ix = 0; ix < vector->n; ix++) { + if (!isfinite(vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + return false; + } + } + return true; } @@ -492,6 +522,5 @@ *norm = normI2 / normI1; - - fprintf(stderr, "Sums: %f %f %f\n", normI1, normI2, *norm); + fprintf(stderr, "Sums: %f %f %f - I1I1: %f\n", normI1, normI2, *norm, sumI1I1); if (mode & PM_SUBTRACTION_EQUATION_NORM) { @@ -507,4 +536,22 @@ matrix->data.F64[normIndex][bgIndex] = sumI1; } + + // check for any NAN values in the result, skip if found: + for (int iy = 0; iy < matrix->numRows; iy++) { + for (int ix = 0; ix < matrix->numCols; ix++) { + if (!isfinite(matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + return false; + } + } + } + for (int ix = 0; ix < vector->n; ix++) { + if (!isfinite(vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + return false; + } + } + + return true; } @@ -512,5 +559,5 @@ #if 1 // Add in penalty term to least-squares vector -static bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term +bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term psVector *vector, // Vector to which to add in penalty term const pmSubtractionKernels *kernels, // Kernel parameters @@ -527,12 +574,36 @@ int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations int numParams = numKernels * numSpatial; // Number of kernel parameters - for (int i = 0; i < numParams; i++) { + + // order is : + // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0] + // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] + // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] + // [norm] + // [bg] + // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] + // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] + // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1] + + for (int i = 0; i < numKernels; i++) { for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { // Contribution to chi^2: a_i^2 P_i psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty"); + fprintf (stderr, "penalty main: %d %e vs %e x %f : %f\n", + index, + matrix->data.F64[index][index], + norm, penalties->data.F32[i], + (norm * penalties->data.F32[i]) / matrix->data.F64[index][index] ); matrix->data.F64[index][index] += norm * penalties->data.F32[i]; if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + fprintf (stderr, "penalty dual: %d %e vs %e x %f : %f\n", + index + numParams + 2, + matrix->data.F64[index + numParams + 2][index + numParams + 2], + norm, penalties->data.F32[i], + (norm * penalties->data.F32[i]) / matrix->data.F64[index + numParams + 2][index + numParams + 2] ); matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i]; + // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) + // penalties scale with second moments + // } } @@ -716,12 +787,12 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1, + status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, weight, window, stamp->convolutions1, kernels, - polyValues, footprint, mode); + polyValues, footprint, stamps->normWindow, mode); break; case PM_SUBTRACTION_MODE_2: - status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2, + status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, weight, window, stamp->convolutions2, kernels, - polyValues, footprint, mode); + polyValues, footprint, stamps->normWindow, mode); break; case PM_SUBTRACTION_MODE_DUAL: @@ -893,349 +964,4 @@ psVectorInit(sumVector, 0.0); psImageInit(sumMatrix, 0.0); - int numStamps = 0; // Number of good stamps - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - - if (stamp->status == PM_SUBTRACTION_STAMP_USED) { - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int iy = 0; iy < stamp->matrix->numRows; iy++) { - for (int ix = 0; ix < stamp->matrix->numCols; ix++) { - if (!isfinite(stamp->matrix->data.F64[iy][ix])) { - fprintf (stderr, "WARNING: NAN in matrix\n"); - } - } - } - for (int ix = 0; ix < stamp->vector->n; ix++) { - if (!isfinite(stamp->vector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - - (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); - numStamps++; - } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); - } - } - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - -#if 0 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); -#endif - -#ifdef TESTING - for (int ix = 0; ix < sumVector->n; ix++) { - if (!isfinite(sumVector->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } - { - psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); - psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); - psFree(inverse); - } - { - psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); - psImage *Xt = psMatrixTranspose(NULL, X); - psImage *XtX = psMatrixMultiply(NULL, Xt, X); - psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); - psFree(X); - psFree(Xt); - psFree(XtX); - } -#endif - -# define SVD_ANALYSIS 0 -# define COEFF_SIG 0.0 -# define SVD_TOL 0.0 - - psVector *solution = NULL; - psVector *solutionErr = NULL; - -#ifdef TESTING - psFitsWriteImageSimple("A.fits", sumMatrix, NULL); - psVectorWriteFile ("B.dat", sumVector); -#endif - - // Use SVD to determine the kernel coeffs (and validate) - if (SVD_ANALYSIS) { - - // We have sumVector and sumMatrix. we are trying to solve the following equation: - // sumMatrix * x = sumVector. - - // we can use any standard matrix inversion to solve this. However, the basis - // functions in general have substantial correlation, so that the solution may be - // somewhat poorly determined or unstable. If not numerically ill-conditioned, the - // system of equations may be statistically ill-conditioned. Noise in the image - // will drive insignificant, but correlated, terms in the solution. To avoid these - // problems, we can use SVD to identify numerically unconstrained values and to - // avoid statistically badly determined value. - - // A = sumMatrix, B = sumVector - // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T - // x = V (1/w) (U^T B) - // \sigma_x = sqrt(diag(A^{-1})) - // solve for x and A^{-1} to get x & dx - // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 - // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 - - // If I use the SVD trick to re-condition the matrix, I need to break out the - // kernel and normalization terms from the background term. - // XXX is this true? or was this due to an error in the analysis? - - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - - // now pull out the kernel elements into their own square matrix - psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); - psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { - if (ix == bgIndex) continue; - for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { - if (iy == bgIndex) continue; - kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; - ky++; - } - kernelVector->data.F64[kx] = sumVector->data.F64[ix]; - kx++; - } - - psImage *U = NULL; - psImage *V = NULL; - psVector *w = NULL; - if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { - psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); - return NULL; - } - - // calculate A_inverse: - // Ainv = V * w * U^T - psImage *wUt = p_pmSubSolve_wUt (w, U); - psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); - psImage *Xvar = NULL; - psFree (wUt); - -# ifdef TESTING - // kernel terms: - for (int i = 0; i < w->n; i++) { - fprintf (stderr, "w: %f\n", w->data.F64[i]); - } -# endif - // loop over w adding in more and more of the values until chisquare is no longer - // dropping significantly. - // XXX this does not seem to work very well: we seem to need all terms even for - // simple cases... - - psVector *Xsvd = NULL; - { - psVector *Ax = NULL; - psVector *UtB = NULL; - psVector *wUtB = NULL; - - psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); - psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); - psVectorInit (wMask, 1); // start by masking everything - - double chiSquareLast = NAN; - int maxWeight = 0; - - double Axx, Bx, y2; - - // XXX this is an attempt to exclude insignificant modes. - // it was not successful with the ISIS kernel set: removing even - // the least significant mode leaves additional ringing / noise - // because the terms are so coupled. - for (int k = 0; false && (k < w->n); k++) { - - // unmask the k-th weight - wMask->data.U8[k] = 0; - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; - chiSquareLast = chiSquare; - - // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); - if (k && !maxWeight && (deltaChi < 1.0)) { - maxWeight = k; - } - } - - // keep all terms or we get extra ringing - maxWeight = w->n; - psVectorInit (wMask, 1); - for (int k = 0; k < maxWeight; k++) { - wMask->data.U8[k] = 0; - } - p_pmSubSolve_SetWeights(wApply, w, wMask); - - // solve for x: - // x = V * w * (U^T * B) - p_pmSubSolve_UtB (&UtB, U, kernelVector); - p_pmSubSolve_wUtB (&wUtB, wApply, UtB); - p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); - - // chi-square for this system of equations: - // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 - // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 - p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); - p_pmSubSolve_VdV (&Axx, Ax, Xsvd); - p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); - p_pmSubSolve_y2 (&y2, kernels, stamps); - - // apparently, this works (compare with the brute force value below - double chiSquare = Axx - 2.0*Bx + y2; - psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); - - // re-calculate A^{-1} to get new variances: - // Ainv = V * w * U^T - // XXX since we keep all terms, this is identical to Ainv - psImage *wUt = p_pmSubSolve_wUt (wApply, U); - Xvar = p_pmSubSolve_VwUt (V, wUt); - psFree (wUt); - - psFree (Ax); - psFree (UtB); - psFree (wUtB); - psFree (wApply); - psFree (wMask); - } - - // copy the kernel solutions to the full solution vector: - solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - - for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { - if (ix == bgIndex) { - solution->data.F64[ix] = 0; - solutionErr->data.F64[ix] = 0.001; - continue; - } - solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); - solution->data.F64[ix] = Xsvd->data.F64[kx]; - kx++; - } - - psFree (kernelMatrix); - psFree (kernelVector); - - psFree (U); - psFree (V); - psFree (w); - - psFree (Ainv); - psFree (Xsvd); - } else { - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); - if (!luMatrix) { - psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - return NULL; - } - - // XXX LUD does not provide A^{-1}? fake the error for now - solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); - for (int ix = 0; ix < sumVector->n; ix++) { - solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; - } - } - - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - } - - // only update the solutions that we chose to calculate: - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - int numKernels = kernels->num; - int spatialOrder = kernels->spatialOrder; // Order of spatial variation - int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms - for (int i = 0; i < numKernels * numPoly; i++) { - // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); - if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { - // XXX fprintf (stderr, "drop\n"); - kernels->solution1->data.F64[i] = 0.0; - } else { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - } - } - } - // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); - // fprintf (stderr, "chi square Brute = %f\n", chiSquare); - - psFree(solution); - psFree(sumVector); - psFree(sumMatrix); - -#ifdef TESTING - // XXX double-check for NAN in data: - for (int ix = 0; ix < kernels->solution1->n; ix++) { - if (!isfinite(kernels->solution1->data.F64[ix])) { - fprintf (stderr, "WARNING: NAN in vector\n"); - } - } -#endif - - } else { - // Dual convolution solution - - // Accumulation of stamp matrices/vectors - psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); - psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); - psImageInit(sumMatrix, 0.0); - psVectorInit(sumVector, 0.0); psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations @@ -1247,18 +973,13 @@ (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - psVectorAppend(norms, stamp->norm); - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; - } - } - -#ifdef TESTING - psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); - psVectorWriteFile("sumVector.dat", sumVector); -#endif - -#if 1 + } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); + } + } + +#if 0 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); @@ -1276,25 +997,4 @@ // Solve normalisation and background separately int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - -#if 0 - psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); - psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); - - normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; - normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; - normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; - - normVector->data.F64[0] = sumVector->data.F64[normIndex]; - normVector->data.F64[1] = sumVector->data.F64[bgIndex]; - - psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); - - double normValue = normSolution->data.F64[0]; - double bgValue = normSolution->data.F64[1]; - - psFree(normMatrix); - psFree(normVector); - psFree(normSolution); -#endif psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm @@ -1309,10 +1009,160 @@ double normValue = stats->robustMedian; -// double bgValue = 0.0; + // double bgValue = 0.0; psFree(stats); fprintf(stderr, "Norm: %lf\n", normValue); -// fprintf(stderr, "BG: %lf\n", bgValue); + // fprintf(stderr, "BG: %lf\n", bgValue); + + // Solve kernel components + for (int i = 0; i < numSolution1; i++) { + sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i]; + + sumMatrix->data.F64[i][normIndex] = 0.0; + sumMatrix->data.F64[normIndex][i] = 0.0; + + // sumMatrix->data.F64[i][bgIndex] = 0.0; + // sumMatrix->data.F64[bgIndex][i] = 0.0; + } + + sumMatrix->data.F64[normIndex][normIndex] = 1.0; + // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0; + // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + + sumVector->data.F64[normIndex] = 0.0; + // sumVector->data.F64[bgIndex] = 0.0; + + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + + solution->data.F64[normIndex] = normValue; + // solution->data.F64[bgIndex] = bgValue; + } +# endif + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int ix = 0; ix < kernels->solution1->n; ix++) { + if (!isfinite(kernels->solution1->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + + } else { + // Dual convolution solution + + // Accumulation of stamp matrices/vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix, 0.0); + psVectorInit(sumVector, 0.0); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int numStamps = 0; // Number of good stamps + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + + psVectorAppend(norms, stamp->norm); + + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); + numStamps++; + } + } + +#ifdef TESTING + psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); + psVectorWriteFile("sumVector.dat", sumVector); +#endif + +#if 1 + // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0); +#endif + + psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams, PS_TYPE_F64); + psVectorInit(solution, 0); + +#if 0 + // Regular, straight-forward solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +#else + { + // Solve normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + +#if 0 + psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); + psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); + + normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; + normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; + normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; + + normVector->data.F64[0] = sumVector->data.F64[normIndex]; + normVector->data.F64[1] = sumVector->data.F64[bgIndex]; + + psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); + + double normValue = normSolution->data.F64[0]; + double bgValue = normSolution->data.F64[1]; + + psFree(normMatrix); + psFree(normVector); + psFree(normSolution); +#endif + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(sumMatrix); + psFree(sumVector); + psFree(norms); + return false; + } + + double normValue = stats->robustMedian; + // double bgValue = 0.0; + + psFree(stats); + + fprintf(stderr, "Norm: %lf\n", normValue); + // fprintf(stderr, "BG: %lf\n", bgValue); // Solve kernel components @@ -1324,23 +1174,23 @@ sumMatrix->data.F64[normIndex][i] = 0.0; -// sumMatrix->data.F64[i][bgIndex] = 0.0; -// sumMatrix->data.F64[bgIndex][i] = 0.0; + // sumMatrix->data.F64[i][bgIndex] = 0.0; + // sumMatrix->data.F64[bgIndex][i] = 0.0; sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; -// sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0; -// sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0; + // sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0; + // sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0; } sumMatrix->data.F64[normIndex][normIndex] = 1.0; -// sumMatrix->data.F64[bgIndex][bgIndex] = 1.0; -// sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0; + // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0; + // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0; sumVector->data.F64[normIndex] = 0.0; -// sumVector->data.F64[bgIndex] = 0.0; + // sumVector->data.F64[bgIndex] = 0.0; solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); solution->data.F64[normIndex] = normValue; -// solution->data.F64[bgIndex] = bgValue; + // solution->data.F64[bgIndex] = bgValue; } #endif @@ -1447,4 +1297,9 @@ } float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); + if (!isfinite(sum)) return false; + if (!isfinite(dmax)) return false; + if (!isfinite(dmin)) return false; + if (!isfinite(peak)) return false; + // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak); psVectorAppend(fSigRes, sigma/sum); @@ -2056,2 +1911,311 @@ } + +# if 0 + +#ifdef TESTING + psFitsWriteImageSimple("A.fits", sumMatrix, NULL); + psVectorWriteFile ("B.dat", sumVector); +#endif + +# define SVD_ANALYSIS 0 +# define COEFF_SIG 0.0 +# define SVD_TOL 0.0 + + // Use SVD to determine the kernel coeffs (and validate) + if (SVD_ANALYSIS) { + + // We have sumVector and sumMatrix. we are trying to solve the following equation: + // sumMatrix * x = sumVector. + + // we can use any standard matrix inversion to solve this. However, the basis + // functions in general have substantial correlation, so that the solution may be + // somewhat poorly determined or unstable. If not numerically ill-conditioned, the + // system of equations may be statistically ill-conditioned. Noise in the image + // will drive insignificant, but correlated, terms in the solution. To avoid these + // problems, we can use SVD to identify numerically unconstrained values and to + // avoid statistically badly determined value. + + // A = sumMatrix, B = sumVector + // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T + // x = V (1/w) (U^T B) + // \sigma_x = sqrt(diag(A^{-1})) + // solve for x and A^{-1} to get x & dx + // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 + // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 + + // If I use the SVD trick to re-condition the matrix, I need to break out the + // kernel and normalization terms from the background term. + // XXX is this true? or was this due to an error in the analysis? + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + + // now pull out the kernel elements into their own square matrix + psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); + psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { + if (ix == bgIndex) continue; + for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { + if (iy == bgIndex) continue; + kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; + ky++; + } + kernelVector->data.F64[kx] = sumVector->data.F64[ix]; + kx++; + } + + psImage *U = NULL; + psImage *V = NULL; + psVector *w = NULL; + if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { + psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); + return NULL; + } + + // calculate A_inverse: + // Ainv = V * w * U^T + psImage *wUt = p_pmSubSolve_wUt (w, U); + psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); + psImage *Xvar = NULL; + psFree (wUt); + +# ifdef TESTING + // kernel terms: + for (int i = 0; i < w->n; i++) { + fprintf (stderr, "w: %f\n", w->data.F64[i]); + } +# endif + // loop over w adding in more and more of the values until chisquare is no longer + // dropping significantly. + // XXX this does not seem to work very well: we seem to need all terms even for + // simple cases... + + psVector *Xsvd = NULL; + { + psVector *Ax = NULL; + psVector *UtB = NULL; + psVector *wUtB = NULL; + + psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); + psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); + psVectorInit (wMask, 1); // start by masking everything + + double chiSquareLast = NAN; + int maxWeight = 0; + + double Axx, Bx, y2; + + // XXX this is an attempt to exclude insignificant modes. + // it was not successful with the ISIS kernel set: removing even + // the least significant mode leaves additional ringing / noise + // because the terms are so coupled. + for (int k = 0; false && (k < w->n); k++) { + + // unmask the k-th weight + wMask->data.U8[k] = 0; + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; + chiSquareLast = chiSquare; + + // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); + if (k && !maxWeight && (deltaChi < 1.0)) { + maxWeight = k; + } + } + + // keep all terms or we get extra ringing + maxWeight = w->n; + psVectorInit (wMask, 1); + for (int k = 0; k < maxWeight; k++) { + wMask->data.U8[k] = 0; + } + p_pmSubSolve_SetWeights(wApply, w, wMask); + + // solve for x: + // x = V * w * (U^T * B) + p_pmSubSolve_UtB (&UtB, U, kernelVector); + p_pmSubSolve_wUtB (&wUtB, wApply, UtB); + p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); + + // chi-square for this system of equations: + // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 + // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 + p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); + p_pmSubSolve_VdV (&Axx, Ax, Xsvd); + p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); + p_pmSubSolve_y2 (&y2, kernels, stamps); + + // apparently, this works (compare with the brute force value below + double chiSquare = Axx - 2.0*Bx + y2; + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); + + // re-calculate A^{-1} to get new variances: + // Ainv = V * w * U^T + // XXX since we keep all terms, this is identical to Ainv + psImage *wUt = p_pmSubSolve_wUt (wApply, U); + Xvar = p_pmSubSolve_VwUt (V, wUt); + psFree (wUt); + + psFree (Ax); + psFree (UtB); + psFree (wUtB); + psFree (wApply); + psFree (wMask); + } + + // copy the kernel solutions to the full solution vector: + solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + + for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { + if (ix == bgIndex) { + solution->data.F64[ix] = 0; + solutionErr->data.F64[ix] = 0.001; + continue; + } + solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); + solution->data.F64[ix] = Xsvd->data.F64[kx]; + kx++; + } + + psFree (kernelMatrix); + psFree (kernelVector); + + psFree (U); + psFree (V); + psFree (w); + + psFree (Ainv); + psFree (Xsvd); + } else { + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + + solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); + psFree(luMatrix); + psFree(permutation); + if (!solution) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + return NULL; + } + + // XXX LUD does not provide A^{-1}? fake the error for now + solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); + for (int ix = 0; ix < sumVector->n; ix++) { + solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; + } + } + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); + psVectorInit (kernels->solution1, 0.0); + } + + // only update the solutions that we chose to calculate: + if (mode & PM_SUBTRACTION_EQUATION_NORM) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_BG) { + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; + } + if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { + int numKernels = kernels->num; + int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { + // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); + if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { + // XXX fprintf (stderr, "drop\n"); + kernels->solution1->data.F64[i] = 0.0; + } else { + // XXX fprintf (stderr, "keep\n"); + kernels->solution1->data.F64[i] = solution->data.F64[i]; + } + } + } + // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); + // fprintf (stderr, "chi square Brute = %f\n", chiSquare); + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); +# endif + +#ifdef TESTING + // XXX double-check for NAN in data: + for (int iy = 0; iy < stamp->matrix->numRows; iy++) { + for (int ix = 0; ix < stamp->matrix->numCols; ix++) { + if (!isfinite(stamp->matrix->data.F64[iy][ix])) { + fprintf (stderr, "WARNING: NAN in matrix\n"); + } + } + } + for (int ix = 0; ix < stamp->vector->n; ix++) { + if (!isfinite(stamp->vector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } +#endif + +#ifdef TESTING + for (int ix = 0; ix < sumVector->n; ix++) { + if (!isfinite(sumVector->data.F64[ix])) { + fprintf (stderr, "WARNING: NAN in vector\n"); + } + } + { + psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); + psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); + psFree(inverse); + } + { + psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); + psImage *Xt = psMatrixTranspose(NULL, X); + psImage *XtX = psMatrixMultiply(NULL, Xt, X); + psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); + psFree(X); + psFree(Xt); + psFree(XtX); + } +#endif +