Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 29004) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 29543) @@ -19,6 +19,10 @@ //#define USE_WEIGHT // Include weight (1/variance) in equation? -//#define USE_WINDOW // Include weight (1/variance) in equation? - + +// XXX TEST: +//# define USE_WINDOW // window to avoid neighbor contamination + +# define PENALTY false +# define MOMENTS (!PENALTY) ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -27,7 +31,7 @@ // Calculate the least-squares matrix and vector -static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated - psVector *vector, // Least-squares vector, updated - double *norm, // Normalisation, updated +static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated + psVector *vector, // Least-squares vector, updated + double normValue, // Normalisation, supplied const psKernel *input, // Input image (target) const psKernel *reference, // Reference image (convolution source) @@ -37,8 +41,5 @@ const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - int normWindow1, // Window (half-)size for normalisation measurement - int normWindow2, // Window (half-)size for normalisation measurement - const pmSubtractionEquationCalculationMode mode + int footprint // (Half-)Size of stamp ) { @@ -50,13 +51,21 @@ int numKernels = kernels->num; // Number of kernels - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms double poly[numPoly]; // Polynomial terms double poly2[numPoly][numPoly]; // Polynomial-polynomial values + int numParams = numKernels * numPoly; + + psAssert(matrix && + matrix->type.type == PS_TYPE_F64 && + matrix->numCols == numParams && + matrix->numRows == numParams, + "Least-squares matrix is bad."); + psAssert(vector && + vector->type.type == PS_TYPE_F64 && + vector->n == numParams, + "Least-squares vector is bad."); // Evaluate polynomial-polynomial terms - // XXX we can skip this if we are not calculating kernel coeffs for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { @@ -84,6 +93,4 @@ // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] - // normalization - // bg 0, bg 1, bg 2 (only 0 is currently used?) for (int i = 0; i < numKernels; i++) { @@ -107,18 +114,15 @@ // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double value = sumCC * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex] = value; - matrix->data.F64[jIndex][iIndex] = value; - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double value = sumCC * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex] = value; + matrix->data.F64[jIndex][iIndex] = value; + } + } } double sumRC = 0.0; // Sum of the reference-convolution products double sumIC = 0.0; // Sum of the input-convolution products - double sumC = 0.0; // Sum of the convolution for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -128,10 +132,8 @@ double ic = in * conv; double rc = ref * conv; - double c = conv; if (weight) { float wtVal = weight->kernel[y][x]; ic *= wtVal; rc *= wtVal; - c *= wtVal; } if (window) { @@ -139,97 +141,14 @@ ic *= winVal; rc *= winVal; - c *= winVal; } sumIC += ic; sumRC += rc; - sumC += c; - } - } + } + } + // Spatial variation for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - double normTerm = sumRC * poly[iTerm]; - double bgTerm = sumC * poly[iTerm]; - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = normTerm; - matrix->data.F64[normIndex][iIndex] = normTerm; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = bgTerm; - matrix->data.F64[bgIndex][iIndex] = bgTerm; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = sumIC * poly[iTerm]; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * normTerm; - } - } - } - } - - double sumRR = 0.0; // Sum of the reference product - double sumIR = 0.0; // Sum of the input-reference product - double sum1 = 0.0; // Sum of the background - 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++) { - double in = input->kernel[y][x]; - double ref = reference->kernel[y][x]; - double ir = in * ref; - double rr = PS_SQR(ref); - double one = 1.0; - - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { - normI1 += ref; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { - normI2 += in; - } - - if (weight) { - float wtVal = weight->kernel[y][x]; - rr *= wtVal; - ir *= wtVal; - in *= wtVal; - ref *= wtVal; - one *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - rr *= winVal; - ir *= winVal; - in *= winVal; - ref *= winVal; - one *= winVal; - } - sumRR += rr; - sumIR += ir; - sumR += ref; - sumI += in; - sum1 += one; - } - } - - *norm = normI2 / normI1; - - fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); - - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumRR; - vector->data.F64[normIndex] = sumIR; - // subtract sum over kernels * kernel solution - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI; - } - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[normIndex][bgIndex] = sumR; - matrix->data.F64[bgIndex][normIndex] = sumR; + vector->data.F64[iIndex] = (sumIC - normValue*sumRC) * poly[iTerm]; + } } @@ -255,33 +174,31 @@ // Calculate the least-squares matrix and vector for dual convolution -static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated - psVector *vector, // Least-squares vector, updated - double *norm, // Normalisation, updated - const psKernel *image1, // Image 1 - const psKernel *image2, // Image 2 +// XXX we could avoid calculating these values on successive passes *if* the stamp has not changed. +static bool calculateDualMatrixVector(pmSubtractionStamp *stamp, // stamp of interest + double normValue, // Normalisation, updated + double normValue2, // Normalisation, updated const psKernel *weight, // Weight image const psKernel *window, // Window image - const psArray *convolutions1, // Convolutions of image 1 for each kernel - const psArray *convolutions2, // Convolutions of image 2 for each kernel const pmSubtractionKernels *kernels, // Kernels const psImage *polyValues, // Spatial polynomial values - int footprint, // (Half-)Size of stamp - int normWindow1, // Window (half-)size for normalisation measurement - int normWindow2, // Window (half-)size for normalisation measurement - const pmSubtractionEquationCalculationMode mode + int footprint // (Half-)Size of stamp ) { - int numKernels = kernels->num; // Number of kernels - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background - int spatialOrder = kernels->spatialOrder; // Order of spatial variation + int numKernels = kernels->num; // Number of kernels + int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms double poly[numPoly]; // Polynomial terms double poly2[numPoly][numPoly]; // Polynomial-polynomial values - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters - int numParams2 = numKernels * numPoly; // Number of additional parameters for dual - int numDual = numParams + numParams2; // Total number of parameters for dual + int numParams = numKernels * numPoly; // Number of regular parameters + int numParams2 = numKernels * numPoly; // Number of additional parameters for dual + int numDual = numParams + numParams2; // Total number of parameters for dual + + psImage *matrix = stamp->matrix; // Least-squares matrix, updated + psVector *vector = stamp->vector; // Least-squares vector, updated + psKernel *image1 = stamp->image1; // Image 1 + psKernel *image2 = stamp->image2; // Image 2 + psArray *convolutions1 = stamp->convolutions1; // Convolutions of image 1 for each kernel + psArray *convolutions2 = stamp->convolutions2; // Convolutions of image 2 for each kernel psAssert(matrix && @@ -290,4 +207,5 @@ matrix->numRows == numDual, "Least-squares matrix is bad."); + psAssert(vector && vector->type.type == PS_TYPE_F64 && @@ -318,4 +236,18 @@ } + // the order of the elements in the matrix and vector is: + // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] + // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] + // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] + + // for DUAL convolution analysis, we apply the normalization to I1 as follows: + // norm = I2 / I1 + // + // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i + // I2c = I2 + \Sum_i b_i I2 \cross k_i + + // we cannot absorb the normalization into a_i until the analysis is complete, or the + // second moment terms are incorrectly calculated. + for (int i = 0; i < numKernels; i++) { psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i @@ -328,9 +260,14 @@ double sumBB = 0.0; // Sum of convolution products between image 2 double sumAB = 0.0; // Sum of convolution products across images 1 and 2 + + double MxxAA = 0.0; + double MyyAA = 0.0; + double MxxBB = 0.0; + double MyyBB = 0.0; for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x]; + double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue); double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; - double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; if (weight) { float wtVal = weight->kernel[y][x]; @@ -348,27 +285,76 @@ sumBB += bb; sumAB += ab; - } - } + + if (MOMENTS) { + MxxAA += x*x*aa; + MyyAA += y*y*aa; + MxxBB += x*x*bb; + MyyBB += y*y*bb; + } + } + } + + if (MOMENTS) { + MxxAA /= stamp->normSquare1 * PS_SQR(normValue); + MyyAA /= stamp->normSquare1 * PS_SQR(normValue); + MxxBB /= stamp->normSquare2; + MyyBB /= stamp->normSquare2; + } + + // XXX this makes the Chisq portion independent of the normalization and star flux + // but may be mis-scaling between stars of different fluxes + sumAA /= PS_SQR(stamp->normI2); + sumAB /= PS_SQR(stamp->normI2); + sumBB /= PS_SQR(stamp->normI2); + + // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB); // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double aa = sumAA * poly2[iTerm][jTerm]; - double bb = sumBB * poly2[iTerm][jTerm]; - double ab = sumAB * poly2[iTerm][jTerm]; - - matrix->data.F64[iIndex][jIndex] = aa; - matrix->data.F64[jIndex][iIndex] = aa; - - matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; - matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; - - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double aa = sumAA * poly2[iTerm][jTerm]; + double bb = sumBB * poly2[iTerm][jTerm]; + double ab = sumAB * poly2[iTerm][jTerm]; + + matrix->data.F64[iIndex][jIndex] = aa; + matrix->data.F64[jIndex][iIndex] = aa; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; + matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; + + // add in second moments + if (MOMENTS) { + matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA; + matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA; + + matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA; + matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA; + + matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB; + matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB; + + matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB; + matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB; + + // XXX this uses the single moments, first try + // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j]; + // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j]; + // + // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j]; + // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j]; + // + // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j]; + // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j]; + // + // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j]; + // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j]; + } + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } + } + } + + // we need to calculate the lower-diagonal AB elements since they are not symmetric for A <-> B for (int j = 0; j < i; j++) { psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j @@ -376,5 +362,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; + double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; if (weight) { ab *= weight->kernel[y][x]; @@ -387,13 +373,15 @@ } + // XXX this makes the Chisq portion independent of the normalization and star flux + // but may be mis-scaling between stars of different fluxes + sumAB /= PS_SQR(stamp->normI2); + // Spatial variation of kernel coeffs - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { - for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { - double ab = sumAB * poly2[iTerm][jTerm]; - matrix->data.F64[iIndex][jIndex + numParams] = ab; - matrix->data.F64[jIndex + numParams][iIndex] = ab; - } - } + for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { + for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { + double ab = sumAB * poly2[iTerm][jTerm]; + matrix->data.F64[iIndex][jIndex + numParams] = ab; + matrix->data.F64[jIndex + numParams][iIndex] = ab; + } } } @@ -402,8 +390,10 @@ double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) - double sumA = 0.0; // Sum of A (for matrix, background) double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) - double sumB = 0.0; // Sum of B products (for matrix, background) - double sumI2 = 0.0; // Sum of I_2 (for vector, background) + + double MxxAI1 = 0.0; + double MyyAI1 = 0.0; + double MxxBI2 = 0.0; + double MyyBI2 = 0.0; for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -413,8 +403,8 @@ float i2 = image2->kernel[y][x]; - double ai2 = a * i2; + double ai2 = a * i2 * normValue; double bi2 = b * i2; - double ai1 = a * i1; - double bi1 = b * i1; + double ai1 = a * i1 * PS_SQR(normValue); + double bi1 = b * i1 * normValue; if (weight) { @@ -424,7 +414,4 @@ ai1 *= wtVal; bi1 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; } if (window) { @@ -434,17 +421,35 @@ ai1 *= wtVal; bi1 *= wtVal; - a *= wtVal; - b *= wtVal; - i2 *= wtVal; } sumAI2 += ai2; sumBI2 += bi2; sumAI1 += ai1; - sumA += a; sumBI1 += bi1; - sumB += b; - sumI2 += i2; - } - } + + if (MOMENTS) { + MxxAI1 += x*x*ai1; + MyyAI1 += y*y*ai1; + MxxBI2 += x*x*bi2; + MyyBI2 += y*y*bi2; + } + } + } + + if (MOMENTS) { + MxxAI1 /= stamp->normSquare1 * PS_SQR(normValue); + MyyAI1 /= stamp->normSquare1 * PS_SQR(normValue); + MxxBI2 /= stamp->normSquare2; + MyyBI2 /= stamp->normSquare2; + } + + // fprintf (stderr, "i : %d : M(xx,yy)(AI1,BI2) : %f %f %f %f\n", i, MxxAI1, MyyAI1, MxxBI2, MyyBI2); + + // XXX this makes the Chisq portion independent of the normalization and star flux + // but may be mis-scaling between stars of different fluxes + sumAI1 /= PS_SQR(stamp->normI2); + sumBI1 /= PS_SQR(stamp->normI2); + sumAI2 /= PS_SQR(stamp->normI2); + sumBI2 /= PS_SQR(stamp->normI2); + // Spatial variation for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { @@ -452,95 +457,19 @@ double bi2 = sumBI2 * poly[iTerm]; double ai1 = sumAI1 * poly[iTerm]; - double a = sumA * poly[iTerm]; double bi1 = sumBI1 * poly[iTerm]; - double b = sumB * poly[iTerm]; - - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][normIndex] = ai1; - matrix->data.F64[normIndex][iIndex] = ai1; - matrix->data.F64[iIndex + numParams][normIndex] = bi1; - matrix->data.F64[normIndex][iIndex + numParams] = bi1; - } - if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { - matrix->data.F64[iIndex][bgIndex] = a; - matrix->data.F64[bgIndex][iIndex] = a; - matrix->data.F64[iIndex + numParams][bgIndex] = b; - matrix->data.F64[bgIndex][iIndex + numParams] = b; - } - if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { - vector->data.F64[iIndex] = ai2; - vector->data.F64[iIndex + numParams] = bi2; - if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { - // subtract norm * sumRC * poly[iTerm] - psAssert (kernels->solution1, "programming error: define solution first!"); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation - vector->data.F64[iIndex] -= norm * ai1; - vector->data.F64[iIndex + numParams] -= norm * bi1; - } - } - } - } - - double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) - double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) - double sum1 = 0.0; // Sum of 1 (for matrix, background-background) - double sumI2 = 0.0; // Sum of I_2 (for vector, background) - double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) - 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++) { - double i1 = image1->kernel[y][x]; - double i2 = image2->kernel[y][x]; - - double i1i1 = i1 * i1; - double one = 1.0; - double i1i2 = i1 * i2; - - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { - normI1 += i1; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { - normI2 += i2; - } - - if (weight) { - float wtVal = weight->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - i1 *= wtVal; - i1i1 *= wtVal; - one *= wtVal; - i2 *= wtVal; - i1i2 *= wtVal; - } - sumI1 += i1; - sumI1I1 += i1i1; - sum1 += one; - sumI2 += i2; - sumI1I2 += i1i2; - } - } - - *norm = normI2 / normI1; - fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); - - if (mode & PM_SUBTRACTION_EQUATION_NORM) { - matrix->data.F64[normIndex][normIndex] = sumI1I1; - vector->data.F64[normIndex] = sumI1I2; - } - if (mode & PM_SUBTRACTION_EQUATION_BG) { - matrix->data.F64[bgIndex][bgIndex] = sum1; - vector->data.F64[bgIndex] = sumI2; - } - if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { - matrix->data.F64[bgIndex][normIndex] = sumI1; - matrix->data.F64[normIndex][bgIndex] = sumI1; + vector->data.F64[iIndex] = ai2 - ai1; + vector->data.F64[iIndex + numParams] = bi2 - bi1; + + // fprintf (stderr, "i : %d : V(I1,I2) : %f %f\n", i, vector->data.F64[iIndex], vector->data.F64[iIndex + numParams]); + + // add in second moments + if (MOMENTS) { + vector->data.F64[iIndex] -= kernels->penalty * MxxAI1; + vector->data.F64[iIndex] -= kernels->penalty * MyyAI1; + + vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2; + vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2; + } + } } @@ -565,12 +494,14 @@ } -#if 1 // Add in penalty term to least-squares vector 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 - float norm // Normalisation + float normSquare1, // Normalisation for image 1 + float normSquare2 // Normalisation for image 2 ) { + psAssert (kernels->mode == PM_SUBTRACTION_MODE_DUAL, "only use penalties for dual convolution"); + if (kernels->penalty == 0.0) { return true; @@ -589,6 +520,5 @@ // [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] @@ -600,21 +530,19 @@ // Contribution to chi^2: a_i^2 P_i psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty"); - fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]); - matrix->data.F64[index][index] += norm * penalties1->data.F32[i]; - if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], - matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]); - matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i]; - // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) - // penalties scale with second moments - // - } - } - } - } - - return true; -} -# endif + // fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]); + matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i]; + + // fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], + // matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]); + matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i]; + // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) + // penalties scale with second moments + // + } + } + } + + return true; +} ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -647,12 +575,5 @@ int index, bool wantDual) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); - PS_ASSERT_INT_POSITIVE(index, NAN); -#endif - psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, @@ -662,11 +583,5 @@ double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); -#endif - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation return kernels->solution1->data.F64[normIndex]; @@ -676,11 +591,5 @@ const psImage *polyValues) { -#if 0 // This is probably in a tight loop, so don't check inputs - PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); - PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); - PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); -#endif - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); @@ -690,4 +599,208 @@ // Public functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateMoments( + pmSubtractionKernels *kernels, // Kernels + pmSubtractionStampList *stamps) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + // XXX skip this, right? + return true; + + // these are only used by DUAL mode + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true; + + psTimerStart("pmSubtractionCalculateMoments"); + + int footprint = stamps->footprint; // Half-size of stamps + + // Loop over each stamp and calculate its normalization factor + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; + if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; + + pmSubtractionCalculateMomentsStamp(kernels, stamp, footprint, stamps->normWindow2, stamps->normWindow1); + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate moments: %f sec", psTimerClear("pmSubtractionCalculateMoments")); + + return true; +} + +bool pmSubtractionCalculateMomentsStamp( + pmSubtractionKernels *kernels, // Kernels + pmSubtractionStamp *stamp, // stamp on which to save normalization) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ) +{ + double Mxx, Myy; + + int numKernels = kernels->num; + + // Generate convolutions: these are generated once and saved + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + psError(psErrorCodeLast(), false, "Unable to convolve stamp"); + return false; + } + + if (!stamp->MxxI1) { + stamp->MxxI1 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + if (!stamp->MyyI1) { + stamp->MyyI1 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + if (!stamp->MxxI2) { + stamp->MxxI2 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + if (!stamp->MyyI2) { + stamp->MyyI2 = psVectorAlloc (numKernels, PS_TYPE_F32); + } + + for (int i = 0; i < numKernels; i++) { + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions1->data[i], footprint, normWindow1); + stamp->MxxI1->data.F32[i] = Mxx / stamp->normI1; + stamp->MyyI1->data.F32[i] = Myy / stamp->normI1; + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions2->data[i], footprint, normWindow2); + stamp->MxxI2->data.F32[i] = Mxx / stamp->normI2; + stamp->MyyI2->data.F32[i] = Myy / stamp->normI2; + } + + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image1, footprint, normWindow1); + stamp->MxxI1raw = Mxx / stamp->normI1; + stamp->MyyI1raw = Myy / stamp->normI1; + + pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image2, footprint, normWindow2); + stamp->MxxI2raw = Mxx / stamp->normI2; + stamp->MyyI2raw = Myy / stamp->normI2; + + // fprintf (stderr, "Mxx I1: %f, Myy I1: %f, Mxx I2: %f, Myy I2: %f\n", stamp->MxxI1raw, stamp->MyyI1raw, stamp->MxxI2raw, stamp->MyyI2raw); + + return true; +} + +bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window) { + + double Sxx = 0.0; + double Syy = 0.0; + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + if (PS_SQR(x) + PS_SQR(y) > PS_SQR(window)) continue; + + double flux = image->kernel[y][x]; + + Sxx += PS_SQR(x) * flux; + Syy += PS_SQR(y) * flux; + } + } + *Mxx = Sxx; + *Myy = Syy; + return true; +} + +///--------- + +bool pmSubtractionCalculateNormalization( + pmSubtractionStampList *stamps, + const pmSubtractionMode mode) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + psTimerStart("pmSubtractionCalculateNormalization"); + + psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + psVector *norm2 = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations + + int footprint = stamps->footprint; // Half-size of stamps + + // Loop over each stamp and calculate its normalization factor + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; + if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; + + // XXX skip this if we have already calculated it? (stamp->norm does not change, just the median statistic) + // XXX maybe not: the star may have changed for a given stamp -- only if norm is reset to NAN can we do this + if (mode == PM_SUBTRACTION_MODE_2) { + pmSubtractionCalculateNormalizationStamp(stamp, stamp->image2, stamp->image1, footprint, stamps->normWindow2, stamps->normWindow1); + } else { + pmSubtractionCalculateNormalizationStamp(stamp, stamp->image1, stamp->image2, footprint, stamps->normWindow1, stamps->normWindow2); + } + psVectorAppend(norms, stamp->norm); + psVectorAppend(norm2, stamp->normSquare2 / stamp->normSquare1); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm + if (!psVectorStats(stats, norms, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(norms); + psFree(norm2); + return false; + } + stamps->normValue = stats->robustMedian; + + psStatsInit(stats); + if (!psVectorStats(stats, norm2, NULL, NULL, 0)) { + psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); + psFree(stats); + psFree(norms); + psFree(norm2); + return false; + } + stamps->normValue2 = stats->robustMedian; + + psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2); + + psFree(stats); + psFree(norms); + psFree(norm2); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization")); + + return true; +} + +bool pmSubtractionCalculateNormalizationStamp( + pmSubtractionStamp *stamp, // stamp on which to save normalization) + const psKernel *image1, // Input image (target) + const psKernel *image2, // Reference image (convolution source) + int footprint, // (Half-)Size of stamp + int normWindow1, // Window (half-)size for normalisation measurement + int normWindow2 // Window (half-)size for normalisation measurement + ) +{ + double normI1 = 0.0; // Sum of I_1 within the normalisation window (aperture) + double normI2 = 0.0; // Sum of I_2 within the normalisation window (aperture) + double normSquare1 = 0.0; // Sum of (I_1)^2 within the normalisation window (aperture) + double normSquare2 = 0.0; // Sum of (I_2)^2 within the normalisation window (aperture) + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double im1 = image1->kernel[y][x]; + double im2 = image2->kernel[y][x]; + + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { + normI1 += im1; + normSquare1 += PS_SQR(im1); + } + if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { + normI2 += im2; + normSquare2 += PS_SQR(im2); + } + } + } + stamp->norm = normI2 / normI1; + stamp->normI1 = normI1; + stamp->normI2 = normI2; + stamp->normSquare1 = normSquare1; + stamp->normSquare2 = normSquare2; + + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); + + return true; +} bool pmSubtractionCalculateEquationThread(psThreadJob *job) @@ -715,5 +828,4 @@ int numKernels = kernels->num; // Number of kernel basis functions int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). @@ -722,5 +834,6 @@ // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the // number of coefficients for the spatial polynomial, normalisation and a constant background offset. - int numParams = numKernels * numSpatial + 1 + numBackground; + // XXX we no longer solve for the normalization and background in the matrix inversion + int numParams = numKernels * numSpatial; if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { // An additional image is convolved @@ -790,18 +903,13 @@ switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: - status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, - weight, window, stamp->convolutions1, kernels, - polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image2, stamp->image1, + weight, window, stamp->convolutions1, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_2: - status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, - weight, window, stamp->convolutions2, kernels, - polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); + status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, + weight, window, stamp->convolutions2, kernels, polyValues, footprint); break; case PM_SUBTRACTION_MODE_DUAL: - status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, - stamp->image1, stamp->image2, - weight, window, stamp->convolutions1, stamp->convolutions2, - kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); + status = calculateDualMatrixVector(stamp, stamps->normValue, stamps->normValue2, weight, window, kernels, polyValues, footprint); break; default: @@ -928,7 +1036,8 @@ int numKernels = kernels->num; // Number of kernel basis functions int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations - int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms - int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for - int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution + // XXX int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + // XXX int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for + int numParams = numKernels * numSpatial; // Number of parameters being solved for + int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { // An additional image is convolved @@ -960,5 +1069,6 @@ if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { - // Accumulate the least-squares matricies and vectors + // Accumulate the least-squares matrices and vectors. These are generated for the + // kernel elements, excluding the background and normalization. psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector @@ -966,7 +1076,83 @@ psImageInit(sumMatrix, 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); + + 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"); + } + } + +#if 0 + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); +#endif + + psVector *solution = NULL; // Solution to equation! + solution = psVectorAlloc(numParams, PS_TYPE_F64); + psVectorInit(solution, 0); + + // XXX TEST: try some constraint on the svd solution + // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + // SINGLE solution + if (1) { + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + } else { + psVector *PERM = NULL; + psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix); + solution = psMatrixLUSolution(solution, LU, sumVector, PERM); + psFree (LU); + psFree (PERM); + } + +# if (0) + for (int i = 0; i < solution->n; i++) { + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]); + } +# endif + + if (!kernels->solution1) { + kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + psVectorInit(kernels->solution1, 0.0); + } + + 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]; + } + + // Apply the normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + kernels->solution1->data.F64[normIndex] = stamps->normValue; + kernels->solution1->data.F64[bgIndex] = 0.0; + + psFree(solution); + psFree(sumVector); + psFree(sumMatrix); + + } else { + // Dual convolution solution + // Accumulate the least-squares matrices and vectors. These are generated for the + // kernel elements, excluding the background and normalization. + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix, 0.0); + psVectorInit(sumVector, 0.0); + + int numStamps = 0; // Number of good stamps + double normSquare1 = 0.0; // Sum of (I_1)^2 over stamps + double normSquare2 = 0.0; // Sum of (I_2)^2 over stamps for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest @@ -974,5 +1160,8 @@ (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - psVectorAppend(norms, stamp->norm); + + normSquare1 += stamp->normSquare1; + normSquare2 += stamp->normSquare2; + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); numStamps++; @@ -983,7 +1172,13 @@ #if 0 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); + psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); + psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); + psVectorWriteFile("sumVector.dat", sumVector); + psFree (save); #endif + + if (PENALTY) { + calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2); + } psVector *solution = NULL; // Solution to equation! @@ -991,212 +1186,8 @@ 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 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm - if (!psVectorStats(stats, norms, NULL, NULL, 0)) { - psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); - psFree(stats); - psFree(sumMatrix); - psFree(sumVector); - psFree(norms); - return false; - } - - // double normValue = 1.0; - double normValue = stats->robustMedian; - // double bgValue = 0.0; - - psFree(stats); - -#ifdef TESTING - fprintf(stderr, "Norm: %lf\n", normValue); -#endif - // Solve kernel components - for (int i = 0; i < numSolution1; i++) { - sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; - - sumMatrix->data.F64[i][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i] = 0.0; - } - sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; - sumMatrix->data.F64[bgIndex][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][bgIndex] = 0.0; - - sumMatrix->data.F64[normIndex][normIndex] = 1.0; - sumVector->data.F64[normIndex] = 0.0; - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - - solution->data.F64[normIndex] = normValue; - } -# endif - -#if (1) - for (int i = 0; i < solution->n; i++) { - fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); - } -#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(norms); - 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] / 100.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 - int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background - -#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(PM_ERR_DATA, false, "Unable to determine median normalisation"); - psFree(stats); - psFree(sumMatrix); - psFree(sumVector); - psFree(norms); - return false; - } - - double normValue = stats->robustMedian; - - psFree(stats); - -#ifdef TESTING - fprintf(stderr, "Norm: %lf\n", normValue); -#endif - - // Solve kernel components - for (int i = 0; i < numSolution2; i++) { - sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; - sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; - - sumMatrix->data.F64[i][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i] = 0.0; - - sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; - } - sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; - sumMatrix->data.F64[bgIndex][normIndex] = 0.0; - sumMatrix->data.F64[normIndex][bgIndex] = 0.0; - - sumMatrix->data.F64[normIndex][normIndex] = 1.0; - - sumVector->data.F64[normIndex] = 0.0; - - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - - solution->data.F64[normIndex] = normValue; - } -#endif - - -#if (1) + // DUAL solution + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); + +#if (0) for (int i = 0; i < solution->n; i++) { fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); @@ -1204,11 +1195,6 @@ #endif - psFree(sumMatrix); - psFree(sumVector); - - psFree(norms); - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64); + kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); psVectorInit (kernels->solution1, 0.0); } @@ -1218,30 +1204,27 @@ } - // 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; - for (int i = 0; i < numKernels * numSpatial; i++) { - // XXX fprintf (stderr, "keep\n"); - kernels->solution1->data.F64[i] = solution->data.F64[i]; - kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; - } - } - - - memcpy(kernels->solution1->data.F64, solution->data.F64, - numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], - numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); - + // for DUAL convolution analysis, we apply the normalization to I1 as follows: + // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i + // I2c = I2 + \Sum_i b_i I2 \cross k_i + + // We absorb the normalization into a_i after the analysis is complete to be consistent + // with the SINGLE definitions of the convolutions + + int numKernels = kernels->num; + for (int i = 0; i < numKernels * numSpatial; i++) { + // we solve for coefficients + kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue; + kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + } + + // Apply the normalisation and background separately + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + kernels->solution1->data.F64[normIndex] = stamps->normValue; + kernels->solution1->data.F64[bgIndex] = 0.0; + + psFree(sumMatrix); + psFree(sumVector); psFree(solution); - } @@ -1265,7 +1248,10 @@ } -bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { - - // XXX measure some useful stats on the residuals +// measure some useful stats on the stamp residuals: +// fResSigma : the residual stdev / total flux +// fResOuter : the residual fabs / total flux for R > 2 pix +// fResTotal : the residual fabs / total flux for R > 0 pix +bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { + float sum = 0.0; float peak = 0.0; @@ -1277,21 +1263,19 @@ } - // only count pixels with more than X% of the source flux - // calculate stdev(dflux) + // init counters + int npix = 0; float dflux1 = 0.0; float dflux2 = 0.0; - int npix = 0; - - float dmax = 0.0; - float dmin = 0.0; + float dOuter = 0.0; + float dTotal = 0.0; for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm); - if (dflux < 0.02*sum) continue; dflux1 += residual->kernel[y][x]; dflux2 += PS_SQR(residual->kernel[y][x]); - dmax = PS_MAX(residual->kernel[y][x], dmax); - dmin = PS_MIN(residual->kernel[y][x], dmin); + dTotal += fabs(residual->kernel[y][x]); + if (hypot(x,y) > 2.0) { + dOuter += fabs(residual->kernel[y][x]); + } npix ++; } @@ -1299,12 +1283,12 @@ 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); - psVectorAppend(fMaxRes, dmax/peak); - psVectorAppend(fMinRes, dmin/peak); + if (!isfinite(dOuter)) return false; + if (!isfinite(dTotal)) return false; + + // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); + psVectorAppend(fResSigma, sigma/sum); + psVectorAppend(fResOuter, dOuter/sum); + psVectorAppend(fResTotal, dTotal/sum); return true; } @@ -1329,7 +1313,7 @@ pmSubtractionVisualShowFitInit (stamps); - psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); - psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); // we want to save the residual images for the 9 brightest stamps. @@ -1443,5 +1427,5 @@ } - pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint); + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint); } else { @@ -1479,5 +1463,5 @@ } - pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint); + pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint); } @@ -1558,27 +1542,27 @@ psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psVectorStats (stats, fSigRes, NULL, NULL, 0); - kernels->fSigResMean = stats->robustMedian; - kernels->fSigResStdev = stats->robustStdev; + psVectorStats (stats, fResSigma, NULL, NULL, 0); + kernels->fResSigmaMean = stats->robustMedian; + kernels->fResSigmaStdev = stats->robustStdev; psStatsInit (stats); - psVectorStats (stats, fMaxRes, NULL, NULL, 0); - kernels->fMaxResMean = stats->robustMedian; - kernels->fMaxResStdev = stats->robustStdev; + psVectorStats (stats, fResOuter, NULL, NULL, 0); + kernels->fResOuterMean = stats->robustMedian; + kernels->fResOuterStdev = stats->robustStdev; psStatsInit (stats); - psVectorStats (stats, fMinRes, NULL, NULL, 0); - kernels->fMinResMean = stats->robustMedian; - kernels->fMinResStdev = stats->robustStdev; + psVectorStats (stats, fResTotal, NULL, NULL, 0); + kernels->fResTotalMean = stats->robustMedian; + kernels->fResTotalStdev = stats->robustStdev; // XXX save these values somewhere - psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f", - kernels->fSigResMean, kernels->fSigResStdev, - kernels->fMaxResMean, kernels->fMaxResStdev, - kernels->fMinResMean, kernels->fMinResStdev); - - psFree (fSigRes); - psFree (fMaxRes); - psFree (fMinRes); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f", + kernels->fResSigmaMean, kernels->fResSigmaStdev, + kernels->fResOuterMean, kernels->fResOuterStdev, + kernels->fResTotalMean, kernels->fResTotalStdev); + + psFree (fResSigma); + psFree (fResOuter); + psFree (fResTotal); psFree (stats); } @@ -1586,5 +1570,4 @@ psFree(residual); psFree(polyValues); - return deviations; @@ -1912,312 +1895,2 @@ return true; } - - -# 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(psErrorCodeLast(), 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(PM_ERR_DATA, 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(PM_ERR_DATA, 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 -