IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36664


Ignore:
Timestamp:
Apr 15, 2014, 10:35:18 AM (12 years ago)
Author:
bills
Message:

last revision to this file was made in error

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotFullForceSummaryReadout.c

    r36662 r36664  
    1010    psF32   fRminorMax;
    1111    psF32   fRminorDel;
    12     psVector    *fRMajor;
    13     psVector    *fRMinor;
     12    psVector    *rMajor;
     13    psVector    *rMinor;
    1414    psArray *zeroPt;
    1515    psArray *exptime;
    1616    psArray *cffTable;
    1717} galaxyShapeOptions;
    18 
    19 
    2018
    2119static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options);
    … …  
    4947    psMetadataAddF32(outputCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)",
    5048        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);
    5152
    5253    // Create objects from the various input's sources
    5354    // loop over the available readouts
    54 
     55//    psArray *readouts = psArrayAlloc(num);
    5556    psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32);
    5657    galaxyShapeOptions options;
    … …  
    9091        psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime;
    9192        fluxScaleFactor->data.F32[index] = scaleFactor;
     93
     94//        readouts->data[index] = readout;
    9295
    9396        // find detections
    … …  
    181184    psVector *sumInvSig2 = NULL;
    182185    psVector *numerator   = NULL;
     186//    psVector *tmp = NULL;
    183187    psF32   totalNPix = 0;
    184188    long    vectorLength = 0;
    … …  
    216220        }
    217221
    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);
    225226                sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
    226227                psVectorInit(sumWeightedFlux, 0.0);
    … …  
    236237                    source->galaxyFits->chisq->n, vectorLength);
    237238
     239
    238240            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);
    247255
    248256                numerator->data.F32[k] += chisq * source->galaxyFits->nPix;
    249 
     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//
    250264                psF32 invSig2 = 1.0 / (dFlux * dFlux);
    251265                sumInvSig2->data.F32[k] += invSig2;
    252266
     267                // sumWeightedFlux = sum ( Flux / dFlux**2 )
     268//                tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp);
     269//                sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp);
    253270                sumWeightedFlux->data.F32[k] += flux * invSig2;
    254271            }
    … …  
    259276        if (vectorLength) {
    260277            outSrc->galaxyFits = pmSourceGalaxyFitsAlloc();
    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 
     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);
    300282            for (int k = 0; k < vectorLength; k++) {
    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                 }
     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
    310297            }
    311298
    … …  
    313300            psFree(sumInvSig2);
    314301            psFree(sumWeightedFlux);
    315 
    316             if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
     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) {
    317315                // copy the best fit params to the model
    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                     }
     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;
    378321                }
    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);
     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;
    384327                }
    385 
    386                 // rescale data
    387                 psF64 majorRange = majorMax - majorMin;
    388                 if (majorRange == 0) {
    389                     majorRange = 1;
     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;
    390334                }
    391                 psF64 minorRange = minorMax - minorMin;
    392                 if (minorRange == 0) {
    393                     minorRange = 1;
     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;
    394340                }
    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
     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            }
    505349        }
    506350    }
    … …  
    534378
    535379    if (makeVectors) {
    536         opt->fRMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    537         opt->fRMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     380        opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     381        opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    538382        int i = 0;
    539383        for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel;
    … …  
    541385            for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel;
    542386                    fRminor += opt->fRminorDel) {
    543                 opt->fRMinor->data.F32[i] = fRminor;
    544                 opt->fRMajor->data.F32[i] = fRmajor;
     387                opt->rMinor->data.F32[i] = fRminor;
     388                opt->rMajor->data.F32[i] = fRmajor;
    545389                i++;
    546390            }
    … …  
    548392        psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials);
    549393    } else {
    550         opt->fRMinor = NULL;
    551         opt->fRMajor = NULL;
     394        opt->rMinor = NULL;
     395        opt->rMajor = NULL;
    552396    }
    553397       
Note: See TracChangeset for help on using the changeset viewer.