IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 15, 2014, 9:51:56 AM (12 years ago)
Author:
bills
Message:

detect failure to find psf model and exit gracefully
if no sources pass selection cuts in lensing code do not fault

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotFullForceSummaryReadout.c

    r36633 r36662  
    1010    psF32   fRminorMax;
    1111    psF32   fRminorDel;
    12     psVector    *rMajor;
    13     psVector    *rMinor;
     12    psVector    *fRMajor;
     13    psVector    *fRMinor;
    1414    psArray *zeroPt;
    1515    psArray *exptime;
    1616    psArray *cffTable;
    1717} galaxyShapeOptions;
     18
     19
    1820
    1921static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options);
    … …  
    4749    psMetadataAddF32(outputCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)",
    4850        outputExptime);
    49     // don't think this one matters
    50 //    psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)",
    51  //       outputExptime);
    5251
    5352    // Create objects from the various input's sources
    5453    // loop over the available readouts
    55 //    psArray *readouts = psArrayAlloc(num);
     54
    5655    psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32);
    5756    galaxyShapeOptions options;
    … …  
    9190        psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime;
    9291        fluxScaleFactor->data.F32[index] = scaleFactor;
    93 
    94 //        readouts->data[index] = readout;
    9592
    9693        // find detections
    … …  
    184181    psVector *sumInvSig2 = NULL;
    185182    psVector *numerator   = NULL;
    186 //    psVector *tmp = NULL;
    187183    psF32   totalNPix = 0;
    188184    long    vectorLength = 0;
    … …  
    220216        }
    221217
     218        // The only parameters that we are summarizing currently are galaxy shapes.
     219        // Continue for inputs that do not have modelEXT
     220        if (!source->modelEXT) continue;
     221
    222222        if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) {
    223223            if (numerator == NULL) {
    224224                vectorLength = source->galaxyFits->chisq->n;
    225                 // tmp = psVectorAlloc(vectorLength, PS_TYPE_F32);
    226225                sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
    227226                psVectorInit(sumWeightedFlux, 0.0);
    … …  
    237236                    source->galaxyFits->chisq->n, vectorLength);
    238237
    239 
    240238            psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
     239
     240            totalNPix += source->galaxyFits->nPix;
    241241
    242242            for (int k = 0; k < vectorLength; k++) {
    … …  
    245245                psF32 flux  = source->galaxyFits->Flux->data.F32[k]  * scaleFactor;
    246246                psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor;
    247                 // for chisq
    248                 // Numerator = sum( nPix[i] * chisq[i] )
    249                 // denominator will be sum( nPix )
    250 
    251 
    252 //                tmp = (psVector *) psBinaryOp(tmp,
    253  //                   source->galaxyFits->chisq, "*", psScalarAlloc(source->galaxyFits->nPix, PS_TYPE_F32));
    254 //                numerator = (psVector *) psBinaryOp(numerator, numerator, "+", tmp);
    255247
    256248                numerator->data.F32[k] += chisq * source->galaxyFits->nPix;
    257                 totalNPix += source->galaxyFits->nPix;
    258 
    259                 // sumInvSig2 = sum( 1 / dFlux**2 )
    260 //                tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->dFlux, "*", source->galaxyFits->dFlux);
    261 //                tmp = (psVector *) psBinaryOp(tmp, psScalarAlloc(1.0, PS_TYPE_F32), "/", tmp);
    262 //                sumInvSig2 = (psVector *) psBinaryOp(sumInvSig2, sumInvSig2, "+", tmp);
    263 //
     249
    264250                psF32 invSig2 = 1.0 / (dFlux * dFlux);
    265251                sumInvSig2->data.F32[k] += invSig2;
    266252
    267                 // sumWeightedFlux = sum ( Flux / dFlux**2 )
    268 //                tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp);
    269 //                sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp);
    270253                sumWeightedFlux->data.F32[k] += flux * invSig2;
    271254            }
    … …  
    276259        if (vectorLength) {
    277260            outSrc->galaxyFits = pmSourceGalaxyFitsAlloc();
    278             outSrc->galaxyFits->nPix = totalNPix;
    279             outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux, vectorLength, PS_TYPE_F32);
    280             outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
    281             outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
     261            outSrc->galaxyFits->nPix  = totalNPix;
     262            psVector *fluxVec =
     263                outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
     264            psVector *dFluxVec =
     265                outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
     266            psVector *chisqVec =
     267                outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
     268
     269            // Get nominal radii paramters from the cff table
     270            if (outSrc->seq >= options->cffTable->n) {
     271                psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
     272                    outSrc->seq, options->cffTable->n);
     273                return NULL;
     274            }
     275            psMetadata *cffRow = options->cffTable->data[outSrc->seq];
     276            if (!cffRow) {
     277                psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
     278                psFree(outSrc);
     279                return NULL;
     280            }
     281            bool mdok;
     282            psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
     283            if (!mdok) {
     284                psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
     285                psFree(outSrc);
     286                return NULL;
     287            }
     288            psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
     289            if (!mdok) {
     290                psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
     291                psFree(outSrc);
     292                return NULL;
     293            }
     294
     295            // fill the summary galaxyFits vectors
     296
     297            int min_k = -1;
     298            psF32 minChisq = NAN;
     299
    282300            for (int k = 0; k < vectorLength; k++) {
    283                 outSrc->galaxyFits->Flux->data.F32[k]  = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k];
    284                 outSrc->galaxyFits->dFlux->data.F32[k] = 1.0 / sumInvSig2->data.F32[k];
    285                 outSrc->galaxyFits->chisq->data.F32[k] = numerator->data.F32[k] / totalNPix;
    286 
    287 #ifdef nodef
    288                 outSrc->galaxyFits->Flux = (psVector *) psBinaryOp(outSrc->galaxyFits->Flux,
    289                         sumWeightedFlux, "/", sumInvSig2);
    290 
    291                 outSrc->galaxyFits->dFlux = (psVector *) psBinaryOp(outSrc->galaxyFits->dFlux,
    292                         psScalarAlloc(1.0, PS_TYPE_F32), "/", sumInvSig2);
    293 
    294                 outSrc->galaxyFits->chisq = (psVector *) psBinaryOp(outSrc->galaxyFits->chisq,
    295                         numerator, "/", psScalarAlloc(totalNPix, PS_TYPE_F32));
    296 #endif
     301                fluxVec->data.F32[k]  = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k];
     302                dFluxVec->data.F32[k] = 1.0 / sumInvSig2->data.F32[k];
     303                psF32 thischisq =
     304                    chisqVec->data.F32[k] = numerator->data.F32[k] / totalNPix;
     305
     306                if (isfinite(thischisq)  && (!isfinite(minChisq) || thischisq < minChisq)) {
     307                    min_k = k;
     308                    minChisq = thischisq;
     309                }
    297310            }
    298311
    … …  
    300313            psFree(sumInvSig2);
    301314            psFree(sumWeightedFlux);
    302 //            psFree(tmp);
    303 
    304             int min_j = -1;
    305             psF32 minChisq = NAN;
    306             psVector *chisq = outSrc->galaxyFits->chisq;
    307             for (int j=0; j < chisq->n; j++) {
    308                 psF32 thischisq = chisq->data.F32[j];
    309                 if (isfinite(thischisq)  && (!isfinite(minChisq) || thischisq < minChisq)) {
    310                     min_j = j;
    311                     minChisq = thischisq;
    312                 }
    313             }
    314             if (min_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
     315
     316            if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
    315317                // copy the best fit params to the model
    316                 psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
    317                 if (outSrc->seq >= options->cffTable->n) {
    318                     psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
    319                         outSrc->seq, options->cffTable->n);
    320                     return NULL;
    321                 }
    322                 psMetadata *row = options->cffTable->data[outSrc->seq];
    323                 if (!row) {
    324                     psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
    325                     psFree(outSrc);
    326                     return NULL;
    327                 }
    328                 bool mdok;
    329                 psF32 rMajor = psMetadataLookupF32(&mdok, row, "R_MAJOR");
    330                 if (!mdok) {
    331                     psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
    332                     psFree(outSrc);
    333                     return NULL;
    334                 }
    335                 psF32 rMinor = psMetadataLookupF32(&mdok, row, "R_MINOR");
    336                 if (!mdok) {
    337                     psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
    338                     psFree(outSrc);
    339                     return NULL;
    340                 }
    341                 axes.major = rMajor * options->rMajor->data.F32[min_j];
    342                 axes.minor = rMinor * options->rMinor->data.F32[min_j];
    343                 pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
    344                 outSrc->modelEXT->chisq = minChisq;
    345                 outSrc->modelEXT->mag = -2.5 * log10(outSrc->galaxyFits->Flux->data.F32[min_j]);
    346                 outSrc->modelEXT->magErr = 1.0 /
    347                     (outSrc->galaxyFits->Flux->data.F32[min_j]/outSrc->galaxyFits->dFlux->data.F32[min_j]); // 1 / SN
    348             }
     318
     319#define NUM_TRIALS_TO_FIT 9
     320                psVector *major = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64);
     321                psVector *minor = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64);
     322                psVector *chisq = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64);
     323                psVector *flux = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64);
     324                psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64);
     325
     326                // fractional radii with the lowest chisq
     327                psF32 fRMajorBest = options->fRMajor->data.F32[min_k];
     328                psF32 fRMinorBest = options->fRMinor->data.F32[min_k];
     329                // XXX: use a recipe parameter
     330                psF32 maxDelta = 0.11;
     331                psF64 majorMin = NAN;
     332                psF64 majorMax = NAN;
     333                psF64 minorMin = NAN;
     334                psF64 minorMax = NAN;
     335                psF64 chisqMin = NAN;
     336                psF64 chisqMax = NAN;
     337                for (int k = 0; k < vectorLength; k++) {
     338
     339                    if (fabs(options->fRMajor->data.F32[k] - fRMajorBest) < maxDelta &&
     340                        fabs(options->fRMinor->data.F32[k] - fRMinorBest) < maxDelta) {
     341
     342                        if (isfinite(chisqVec->data.F32[k]) &&
     343                            isfinite(fluxVec->data.F32[k])  &&
     344                            isfinite(dFluxVec->data.F32[k])) {
     345
     346                            // compute major and minor radius vectors are in pixel space from nominal and trial fractions
     347                            psF64 thisMajor = R_MAJOR * options->fRMajor->data.F32[k];
     348                            if (!isfinite(majorMin) || thisMajor < majorMin) {
     349                                majorMin = thisMajor;
     350                            }
     351                            if (!isfinite(majorMax) || thisMajor > majorMax) {
     352                                majorMax = thisMajor;
     353                            }
     354                            psVectorAppend(major, thisMajor);
     355
     356                            psF64 thisMinor = R_MINOR * options->fRMinor->data.F32[k];
     357                            if (!isfinite(minorMin) || thisMinor < minorMin) {
     358                                minorMin = thisMinor;
     359                            }
     360                            if (!isfinite(minorMax) || thisMinor > minorMax) {
     361                                minorMax = thisMinor;
     362                            }
     363                            psVectorAppend(minor, thisMinor);
     364
     365                            psF64 thisChisq = chisqVec->data.F32[k];
     366                            if (!isfinite(chisqMin) || thisChisq < chisqMin) {
     367                                chisqMin = thisChisq;
     368                            }
     369                            if (!isfinite(chisqMax) || thisChisq > chisqMax) {
     370                                chisqMax = thisChisq;
     371                            }
     372                            psVectorAppend(chisq, thisChisq);
     373
     374                            psVectorAppend(flux,  fluxVec->data.F32[k]);
     375                            psVectorAppend(dFlux, dFluxVec->data.F32[k]);
     376                        }
     377                    }
     378                }
     379
     380                // see if we ever get too few good values (haven't seen this happen)
     381                if (major->n < NUM_TRIALS_TO_FIT) {
     382                    fprintf(stderr, "only found %ld good points wanted %d for seq: %d\n",
     383                        major->n, NUM_TRIALS_TO_FIT, outSrc->seq);
     384                }
     385
     386                // rescale data
     387                psF64 majorRange = majorMax - majorMin;
     388                if (majorRange == 0) {
     389                    majorRange = 1;
     390                }
     391                psF64 minorRange = minorMax - minorMin;
     392                if (minorRange == 0) {
     393                    minorRange = 1;
     394                }
     395                psF64 chisqRange = chisqMax - chisqMin;
     396                if (chisqRange == 0) {
     397                    chisqRange = 1;
     398                }
     399                for (int k = 0; k < major->n; k++) {
     400                    major->data.F64[k] = (major->data.F64[k] - majorMin) / majorRange;
     401                    minor->data.F64[k] = (minor->data.F64[k] - minorMin) / minorRange;
     402                    chisq->data.F64[k] = (chisq->data.F64[k] - chisqMin) / chisqRange;
     403                }
     404
     405                // Fit chisq versus rMajor and rMinor
     406                psPolynomial2D *chisqFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     407                psPolynomial2D *fluxFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     408                chisqFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     409                chisqFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     410                chisqFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     411                fluxFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     412                fluxFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     413                fluxFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     414
     415                bool goodFit = psVectorFitPolynomial2D (chisqFit, NULL, 0xff, chisq, NULL, major, minor);
     416
     417                if (goodFit) {
     418                    // find the axes for the minimum chisq from the formula
     419                    psF64 **coeff = chisqFit->coeff;
     420                    psF64 D = 4*coeff[2][0]*coeff[0][2] - coeff[1][1]*coeff[1][1];
     421                    psF64 major0 = ((coeff[1][1]*coeff[0][1] - 2*coeff[0][2]*coeff[1][0]) / D) * majorRange + majorMin;
     422                    psF64 minor0 = ((coeff[1][1]*coeff[1][0] - 2*coeff[2][0]*coeff[0][1]) / D) * minorRange + minorMin;
     423
     424                    psF64 chisq0 = psPolynomial2DEval(chisqFit, major0, minor0) * chisqRange + chisqMin;
     425
     426                    // flux and dFlux values from entry with lowest chisq
     427                    psF64 fluxBest = fluxVec->data.F32[min_k];
     428                    psF64 dFlux0 = dFluxVec->data.F32[min_k];
     429
     430                    bool goodFluxFit = psVectorFitPolynomial2D(fluxFit, NULL, 0xFF, flux, NULL, major, minor);
     431                    // flux from fit
     432                    psF64 flux0;
     433                    if (goodFluxFit) {
     434                        flux0  = psPolynomial2DEval(fluxFit, major0, minor0);
     435                    } else {
     436                        flux0 = fluxBest;
     437                    }
     438
     439                    // compare to the minimum value in the vector
     440
     441                    fprintf (stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %4.2f |  %6.3f %6.3f %6.3f %4.2f | %7.4f %7.4f %7.1f %7.1f\n",
     442                        outSrc->seq, min_k, major->n, R_MAJOR, fRMajorBest*R_MAJOR, major0, major0/R_MAJOR, R_MINOR, fRMinorBest*R_MINOR, minor0, minor0/R_MINOR, minChisq, chisq0, fluxBest, flux0);
     443
     444                    psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
     445
     446#ifdef USE_MINIMUM_CHISQ
     447                    axes.major = R_MAJOR * fRMajorBest;
     448                    axes.minor = R_MINOR * fRMinorBest;
     449                    outSrc->modelEXT->chisq = minChisq;
     450                    flux0 = fluxBest;
     451#else
     452                    axes.major = major0;
     453                    axes.minor = minor0;
     454                    outSrc->modelEXT->chisq = chisq0;
     455#endif
     456                    pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
     457
     458                    outSrc->modelEXT->mag = -2.5 * log10(flux0);
     459                    outSrc->modelEXT->magErr = dFlux0 / flux0; // 1 / SN
     460                    outSrc->modelEXT->flags |= PM_MODEL_STATUS_FITTED;
     461                } else {
     462#ifdef notdef
     463                    psErrorStackPrint(stderr, NULL);
     464#endif
     465                    outSrc->modelEXT->flags |= PM_MODEL_STATUS_NONCONVERGE;
     466
     467                    fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f bad fit\n",
     468                        outSrc->seq, min_k, major->n, R_MAJOR, fRMajorBest, R_MINOR, fRMinorBest);
     469                    psErrorClear();
     470                    // psFree(outSrc->galaxyFits);
     471
     472                    // XXX drop the model since it is bogus? Perhaps we do this or use the nominal values
     473                    // and add a flag
     474#ifdef notnow
     475                    psFree(outSrc->modelEXT);
     476                    outSrc->modelEXT = 0;
     477                    psFree(outSrc->modelFits);
     478                    outSrc->modelFits = 0;
     479#endif
     480                }
     481                char fn[80];
     482
     483                sprintf(fn, "vectors/%s.%05d.txt", goodFit ? "g" : "b", outSrc->seq);
     484
     485                FILE *f = fopen(fn, "w");
     486                fprintf(f, "#major minor chisq i\n");
     487                for (int k = 0; k < chisq->n; k++) {
     488                    fprintf(f, "%10.6f %10.6f %10.6f %5d\n",
     489                        major->data.F64[k]/R_MAJOR, minor->data.F64[k]/R_MINOR, chisq->data.F64[k], k);
     490                }
     491                fclose(f);
     492                psFree(chisqFit);
     493                psFree(major);
     494                psFree(minor);
     495                psFree(chisq);
     496                psFree(flux);
     497                psFree(dFlux);
     498            }
     499        } else {
     500#ifdef notnow
     501            // drop the extended source information so that we can tell which objects that we fit
     502            psFree(outSrc->modelEXT);
     503            psFree(outSrc->modelFits);
     504#endif
    349505        }
    350506    }
    … …  
    378534
    379535    if (makeVectors) {
    380         opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    381         opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     536        opt->fRMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     537        opt->fRMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    382538        int i = 0;
    383539        for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel;
    … …  
    385541            for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel;
    386542                    fRminor += opt->fRminorDel) {
    387                 opt->rMinor->data.F32[i] = fRminor;
    388                 opt->rMajor->data.F32[i] = fRmajor;
     543                opt->fRMinor->data.F32[i] = fRminor;
     544                opt->fRMajor->data.F32[i] = fRmajor;
    389545                i++;
    390546            }
    … …  
    392548        psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials);
    393549    } else {
    394         opt->rMinor = NULL;
    395         opt->rMajor = NULL;
     550        opt->fRMinor = NULL;
     551        opt->fRMajor = NULL;
    396552    }
    397553       
Note: See TracChangeset for help on using the changeset viewer.