IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 17, 2011, 5:07:46 PM (16 years ago)
Author:
eugene
Message:

test chisq improvements from increasing the spatial order and DUAL vs SINGLE1,2; reject stamps on the basis of the chisq vs flux2 model fit (allowing for systematic errors in the fit); track the best fit with the new pmSubtractionQuality structure; split out pmSubtractionConvolveStamps from pmSubtractionCalculateEquation; replace pmSubtractionCalculateDeviations with pmSubtractionCalculateChisqAndMoments

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

    r30266 r30288  
    2727
    2828static bool useFFT = true;              // Do convolutions using FFT
    29 
    30 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL
    3129
    3230//#define TESTING
    … …  
    473471}
    474472
     473bool pmSubtractionMatchAttempt(pmSubtractionQuality **bestMatch, pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, pmSubtractionMode mode, int spatialOrder, bool final) {
     474
     475    pmSubtractionMode nativeMode = kernels->mode;
     476    pmSubtractionMode nativeOrder = kernels->spatialOrder;
     477
     478    kernels->mode = mode;
     479    kernels->spatialOrder = spatialOrder;
     480
     481    // we always need to recalculate the matrix equation elements...
     482    pmSubtractionStampsResetStatus(stamps);
     483
     484    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
     485    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
     486        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
     487        return false;
     488    }
     489
     490    // step 1: generate the elements of the matrix equation Ax = B
     491    psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
     492    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     493        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     494        return false;
     495    }
     496                   
     497    // step 2: solve the matrix equation Ax = B
     498    psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
     499    if (!pmSubtractionSolveEquation(kernels, stamps)) {
     500        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     501        return false;
     502    }
     503    memCheck("  solve equation");
     504
     505    // calculate the score for this model fit attempt
     506    // XXX store the chisq, flux and moments for stamp rejection
     507    pmSubtractionCalculateChisqAndMoments(bestMatch, stamps, kernels); // Stamp deviations
     508
     509    // display the input and model stamps
     510    pmSubtractionVisualShowFit(stamps, kernels);
     511    pmSubtractionVisualPlotFit(kernels);
     512    pmSubtractionVisualPlotConvKernels(kernels);
     513
     514    // reset the kernel if desired (on final pass, do not reset)
     515    if (!final) {
     516        kernels->mode = nativeMode;
     517        kernels->spatialOrder = nativeOrder;
     518    }
     519    return true;
     520}
    475521
    476522bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
    … …  
    559605    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    560606
     607    pmSubtractionQuality *bestMatch = NULL;
     608
     609    int N_TEST_MODES;
     610    int N_TEST_ORDER = spatialOrder;
     611
     612    pmSubtractionMode TestModes[3];
     613    switch (subMode) {
     614      case PM_SUBTRACTION_MODE_1:
     615        N_TEST_MODES = 1;
     616        TestModes[0] = PM_SUBTRACTION_MODE_1;
     617        break;
     618      case PM_SUBTRACTION_MODE_2:
     619        N_TEST_MODES = 1;
     620        TestModes[0] = PM_SUBTRACTION_MODE_2;
     621        break;
     622      case PM_SUBTRACTION_MODE_DUAL:
     623        N_TEST_MODES = 3;
     624        TestModes[0] = PM_SUBTRACTION_MODE_1;
     625        TestModes[1] = PM_SUBTRACTION_MODE_2;
     626        TestModes[2] = PM_SUBTRACTION_MODE_DUAL;
     627        break;
     628      default:
     629        psError(psErrorCodeLast(), false, "For now, only modes 1, 2, and DUAL are supported.");
     630        goto MATCH_ERROR;
     631    }
     632   
    561633    memCheck("start");
    562634
    … …  
    689761
    690762            int numRejected = -1;       // Number of rejected stamps in each iteration
    691             for (int k = 0; k < iter && numRejected != 0; k++) {
     763            for (int k = 0; (k < iter) && (numRejected != 0); k++) {
    692764                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    693765
    … …  
    711783                }
    712784
    713                 // step 0b : calculate the moments, pass along to the next steps via stamps->normValue
    714                 // XXX this step is now skipped -- delete
    715                 psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
    716                 if (!pmSubtractionCalculateMoments(kernels, stamps)) {
    717                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     785                // on each iteration, we start from scratch
     786                psFree(bestMatch);
     787
     788                // choose the spatial order and subtraction direction (1, 2, dual)
     789                // XXX need to make these respect recipe somewhat
     790                for (int order = 0; order <= N_TEST_ORDER; order++) {
     791                    for (int j = 0; j < N_TEST_MODES; j++) {
     792                        if (!pmSubtractionMatchAttempt(&bestMatch, kernels, stamps, TestModes[j], order, false)) {
     793                            goto MATCH_ERROR;
     794                        }
     795                    }
     796                }
     797               
     798                // reject the deviant stamps based on the stats of the best match
     799                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
     800                numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej);
     801                if (numRejected < 0) {
     802                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    718803                    goto MATCH_ERROR;
    719804                }
    720 
    721                 // step 1: generate the elements of the matrix equation Ax = B
    722                 psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
    723                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    724                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    725                     goto MATCH_ERROR;
    726                 }
    727 
    728                 // step 2: solve the matrix equation Ax = B
    729                 psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
    730                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    731                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    732                     goto MATCH_ERROR;
    733                 }
    734                 memCheck("  solve equation");
    735 
    736                 pmSubtractionVisualPlotConvKernels(kernels);
    737 
    738                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    739                 if (!deviations) {
    740                     psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    741                     goto MATCH_ERROR;
    742                 }
    743                 memCheck("   calculate deviations");
    744 
    745                 pmSubtractionCalculateChisqAndMoments(stamps, kernels); // Stamp deviations
    746 
    747                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    748                 numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    749                 if (numRejected < 0) {
    750                     psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    751                     psFree(deviations);
    752                     goto MATCH_ERROR;
    753                 }
    754                 psFree(deviations);
    755 
    756805                memCheck("  reject stamps");
    757806            }
    758807
    759             // if we hit the max number of iterations and we have rejected stamps, re-solve
    760             if (numRejected > 0) {
    761 
    762                 // step 1: generate the elements of the matrix equation Ax = B
    763                 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
    764                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    765                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    766                     goto MATCH_ERROR;
    767                 }
    768 
    769                 // step 2: solve the matrix equation Ax = B
    770                 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
    771                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    772                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    773                     goto MATCH_ERROR;
    774                 }
    775 
    776                 pmSubtractionVisualPlotConvKernels(kernels);
    777 
    778                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    779                 if (!deviations) {
    780                     psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    781                     goto MATCH_ERROR;
    782                 }
    783                 pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    784                 psFree(deviations);
    785             }
     808            // apply the best fit so we are ready to roll
     809            fprintf (stderr, "applying order: %d, mode: %d\n", bestMatch->spatialOrder, bestMatch->mode);
     810            if (!pmSubtractionMatchAttempt(NULL, kernels, stamps, bestMatch->mode, bestMatch->spatialOrder, true)) {
     811                goto MATCH_ERROR;
     812            }
    786813            psFree(stamps);
    787             stamps = NULL;
    788 
     814            psFree(bestMatch);
    789815            memCheck("solution");
    790816
    … …  
    860886    psFree(variance);
    861887    psFree(rng);
     888    psFree(bestMatch);
    862889    return false;
    863890}
    … …  
    10771104    assert(kernels);
    10781105
     1106    psAbort("this function is not working");
     1107# if (0)
     1108    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
     1109    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
     1110        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
     1111        return false;
     1112    }
     1113
    10791114    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1080     if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
     1115    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    10811116        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    10821117        return false;
    … …  
    10841119
    10851120    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
    1086     if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     1121    if (!pmSubtractionSolveEquation(kernels, stamps)) {
    10871122        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    10881123        return false;
    … …  
    10961131    }
    10971132
     1133    // XXX this needs to be made consistent with the modified 'reject stamps' function
    10981134    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
    10991135    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    … …  
    11081144        // Allow re-fit with reduced stamps set
    11091145        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1110         if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
     1146        if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    11111147            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    11121148            return false;
    … …  
    11141150
    11151151        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    1116         if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
     1152        if (!pmSubtractionSolveEquation(kernels, stamps)) {
    11171153            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    11181154            return false;
    … …  
    11341170        psFree(deviations);
    11351171    }
    1136 
     1172# endif
    11371173    return true;
    11381174}
Note: See TracChangeset for help on using the changeset viewer.