IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36701 for trunk/psphot/src


Ignore:
Timestamp:
Apr 29, 2014, 10:58:27 AM (12 years ago)
Author:
bills
Message:

perform 2D parabolic fit to the chisq distribution around the trial
with the minimum value.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotFullForceSummaryReadout.c

    r36664 r36701  
    44typedef struct {
    55    int     numTrials;
    6     psF32   fRmajorMin;
    7     psF32   fRmajorMax;
    8     psF32   fRmajorDel;
    9     psF32   fRminorMin;
    10     psF32   fRminorMax;
    11     psF32   fRminorDel;
    12     psVector    *rMajor;
    13     psVector    *rMinor;
     6    psF64   fRMajorMin;
     7    psF64   fRMajorMax;
     8    psF64   fRMajorDel;
     9    psF64   fRMinorMin;
     10    psF64   fRMinorMax;
     11    psF64   fRMinorDel;
     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;
    189185
     186    // Loop over sources for this object. Start output source for first input that has
     187    // good pixel frac
    190188    for (int i=0; i < obj->sources->n; i++) {
    191189        pmSource *source = obj->sources->data[i];
    … …  
    220218        }
    221219
     220        // The only parameters that we are summarizing currently are galaxy shapes.
     221        // Continue for inputs that do not have modelEXT which contains the model paramters
     222        if (!source->modelEXT) continue;
     223
    222224        if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) {
    223225            if (numerator == NULL) {
     226                // first source with galaxyFits allocate accumulators
    224227                vectorLength = source->galaxyFits->chisq->n;
    225                 // tmp = psVectorAlloc(vectorLength, PS_TYPE_F32);
    226228                sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
    227229                psVectorInit(sumWeightedFlux, 0.0);
    … …  
    232234            }
    233235
     236            // Die if the lengths of the vectors is not the same in all sources
    234237            psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
    235238                vectorLength, options->numTrials);
    … …  
    237240                    source->galaxyFits->chisq->n, vectorLength);
    238241
    239 
    240242            psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
     243
     244            totalNPix += source->galaxyFits->nPix;
    241245
    242246            for (int k = 0; k < vectorLength; k++) {
    … …  
    245249                psF32 flux  = source->galaxyFits->Flux->data.F32[k]  * scaleFactor;
    246250                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);
    255251
    256252                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 //
     253
    264254                psF32 invSig2 = 1.0 / (dFlux * dFlux);
    265255                sumInvSig2->data.F32[k] += invSig2;
    266256
    267                 // sumWeightedFlux = sum ( Flux / dFlux**2 )
    268 //                tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp);
    269 //                sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp);
    270257                sumWeightedFlux->data.F32[k] += flux * invSig2;
    271258            }
    … …  
    275262    if (outSrc) {
    276263        if (vectorLength) {
     264            // allocate galaxyFits for the output source
    277265            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);
     266            outSrc->galaxyFits->nPix  = totalNPix;
     267            psVector *fluxVec =
     268                outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
     269            psVector *dFluxVec =
     270                outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
     271            psVector *chisqVec =
     272                outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
     273
     274            // Find the nominal radii paramters in the cff table
     275            if (outSrc->seq >= options->cffTable->n) {
     276                psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
     277                    outSrc->seq, options->cffTable->n);
     278                return NULL;
     279            }
     280            psMetadata *cffRow = options->cffTable->data[outSrc->seq];
     281            if (!cffRow) {
     282                psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
     283                psFree(outSrc);
     284                return NULL;
     285            }
     286            bool mdok;
     287            psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
     288            if (!mdok) {
     289                psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
     290                psFree(outSrc);
     291                return NULL;
     292            }
     293            psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
     294            if (!mdok) {
     295                psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
     296                psFree(outSrc);
     297                return NULL;
     298            }
     299
     300            // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq
     301
     302            int min_k = -1;
     303            psF32 minChisq = NAN;
     304
    282305            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
     306                fluxVec->data.F32[k]  = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k];
     307                dFluxVec->data.F32[k] = 1.0 / sumInvSig2->data.F32[k];
     308
     309                psF32 thischisq = chisqVec->data.F32[k] = numerator->data.F32[k] / totalNPix;
     310
     311                if (isfinite(thischisq)  && (!isfinite(minChisq) || thischisq < minChisq)) {
     312                    min_k = k;
     313                    minChisq = thischisq;
     314                }
    297315            }
    298316
    … …  
    300318            psFree(sumInvSig2);
    301319            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;
     320
     321            if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
     322                // copy the best fit params to the model
     323                // fractional radii with the lowest chisq
     324                psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
     325
     326                // examine the params for the trial with minimum chisq.
     327                bool fitMajor = true;
     328                bool fitMinor = true;
     329                psF64 fRMajorBest = options->fRMajor->data.F32[min_k];
     330                if ((fabs(fRMajorBest - options->fRMajorMin) < options->fRMajorDel) ||
     331                    (fabs(fRMajorBest - options->fRMajorMax) < options->fRMajorDel)) {
     332                    fitMajor = false;
    312333                }
    313             }
    314             if (min_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
    315                 // 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;
     334                psF64 fRMinorBest = options->fRMinor->data.F32[min_k];
     335                if ((fabs(fRMinorBest - options->fRMinorMin) < options->fRMinorDel) ||
     336                    (fabs(fRMinorBest - options->fRMinorMax) < options->fRMinorDel)) {
     337                    fitMinor = false;
    321338                }
    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;
     339                // If either major or minor is at one of the limits do not fit report the minimum value
     340                bool useFit = fitMajor && fitMinor;
     341
     342                // save the flux and dFlux values from entry with lowest chisq
     343                psF64 fluxBest = fluxVec->data.F32[min_k];
     344                psF64 dFlux0   = dFluxVec->data.F32[min_k];
     345                psF64 flux0    = NAN;
     346
     347                if (useFit) {
     348                    #define NUM_TRIALS_INIT 9     // this isn't used for anything execpt for the initial vector lengths
     349                    psVector *major = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     350                    psVector *minor = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     351                    psVector *chisq = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     352                    psVector *flux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     353                    psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
     354
     355                    // XXX: use a recipe parameter
     356                    psF64 maxDeltaMaj = 2.2 * options->fRMajorDel;
     357                    psF64 maxDeltaMin = 2.2 * options->fRMinorDel;
     358
     359                    psF64 majorMin = NAN;
     360                    psF64 majorMax = NAN;
     361                    psF64 minorMin = NAN;
     362                    psF64 minorMax = NAN;
     363                    psF64 chisqMin = NAN;
     364                    psF64 chisqMax = NAN;
     365                    for (int k = 0; k < vectorLength; k++) {
     366
     367                        if (fabs(options->fRMajor->data.F32[k] - fRMajorBest) < maxDeltaMaj &&
     368                            fabs(options->fRMinor->data.F32[k] - fRMinorBest) < maxDeltaMin) {
     369
     370                            if (isfinite(chisqVec->data.F32[k]) &&
     371                                isfinite(fluxVec->data.F32[k])  &&
     372                                isfinite(dFluxVec->data.F32[k])) {
     373
     374                                // compute major and minor radius vectors from nominal and trial fractions
     375                                // also find the ranges in the vectors
     376                                psF64 thisMajor = R_MAJOR * options->fRMajor->data.F32[k];
     377                                if (!isfinite(majorMin) || thisMajor < majorMin) {
     378                                    majorMin = thisMajor;
     379                                }
     380                                if (!isfinite(majorMax) || thisMajor > majorMax) {
     381                                    majorMax = thisMajor;
     382                                }
     383                                psVectorAppend(major, thisMajor);
     384
     385                                psF64 thisMinor = R_MINOR * options->fRMinor->data.F32[k];
     386                                if (!isfinite(minorMin) || thisMinor < minorMin) {
     387                                    minorMin = thisMinor;
     388                                }
     389                                if (!isfinite(minorMax) || thisMinor > minorMax) {
     390                                    minorMax = thisMinor;
     391                                }
     392                                psVectorAppend(minor, thisMinor);
     393
     394                                psF64 thisChisq = chisqVec->data.F32[k];
     395                                if (!isfinite(chisqMin) || thisChisq < chisqMin) {
     396                                    chisqMin = thisChisq;
     397                                }
     398                                if (!isfinite(chisqMax) || thisChisq > chisqMax) {
     399                                    chisqMax = thisChisq;
     400                                }
     401                                psVectorAppend(chisq, thisChisq);
     402
     403                                psVectorAppend(flux,  fluxVec->data.F32[k]);
     404                                psVectorAppend(dFlux, dFluxVec->data.F32[k]);
     405                            }
     406                        }
     407                    }
     408
     409                    // see if we ever get too few good values (haven't seen this happen)
     410                    if (major->n < NUM_TRIALS_INIT) {
     411                        fprintf(stderr, "only found %ld good points wanted %d for seq: %d\n",
     412                            major->n, NUM_TRIALS_INIT, outSrc->seq);
     413                    }
     414
     415                    // rescale data - this helps avoid precision errors.
     416                    psF64 majorRange = majorMax - majorMin;
     417                    if (majorRange == 0) {
     418                        majorRange = 1;
     419                    }
     420                    psF64 minorRange = minorMax - minorMin;
     421                    if (minorRange == 0) {
     422                        minorRange = 1;
     423                    }
     424                    psF64 chisqRange = chisqMax - chisqMin;
     425                    if (chisqRange == 0) {
     426                        chisqRange = 1;
     427                    }
     428                    for (int k = 0; k < major->n; k++) {
     429                        major->data.F64[k] = (major->data.F64[k] - majorMin) / majorRange;
     430                        minor->data.F64[k] = (minor->data.F64[k] - minorMin) / minorRange;
     431                        chisq->data.F64[k] = (chisq->data.F64[k] - chisqMin) / chisqRange;
     432                    }
     433
     434                    // Fit chisq versus rMajor and rMinor
     435                    psPolynomial2D *chisqFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     436                    psPolynomial2D *fluxFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
     437                    chisqFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     438                    chisqFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     439                    chisqFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     440                    fluxFit->coeffMask[2][2] = PS_POLY_MASK_SET;
     441                    fluxFit->coeffMask[2][1] = PS_POLY_MASK_SET;
     442                    fluxFit->coeffMask[1][2] = PS_POLY_MASK_SET;
     443
     444                    bool goodFit = psVectorFitPolynomial2D (chisqFit, NULL, 0xff, chisq, NULL, major, minor);
     445
     446                    if (goodFit) {
     447                        // find the axes for the minimum chisq from the formula
     448                        psF64 **coeff = chisqFit->coeff;
     449                        psF64 D = 4*coeff[2][0]*coeff[0][2] - coeff[1][1]*coeff[1][1];
     450                        psF64 major0 = ((coeff[1][1]*coeff[0][1] - 2*coeff[0][2]*coeff[1][0]) / D) * majorRange + majorMin;
     451                        psF64 minor0 = ((coeff[1][1]*coeff[1][0] - 2*coeff[2][0]*coeff[0][1]) / D) * minorRange + minorMin;
     452
     453                        // estimated chisq at the minimum
     454                        psF64 chisq0 = psPolynomial2DEval(chisqFit, major0, minor0) * chisqRange + chisqMin;
     455
     456                        // now fit the flux ...
     457                        bool goodFluxFit = psVectorFitPolynomial2D(fluxFit, NULL, 0xFF, flux, NULL, major, minor);
     458
     459                        // .. and compute flux at the minimum chisq position
     460                        if (goodFluxFit) {
     461                            flux0  = psPolynomial2DEval(fluxFit, major0, minor0);
     462                        } else {
     463                            flux0 = fluxBest;
     464                        }
     465
     466#ifdef PRINTVALS
     467                        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",
     468                            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);
     469#endif
     470
     471
     472                        axes.major = major0;
     473                        axes.minor = minor0;
     474                        outSrc->modelEXT->chisq = chisq0;
     475#ifdef DUMPVECTORS
     476                        char fn[80];
     477
     478                        sprintf(fn, "vectors/%s.%05d.txt", goodFit ? "g" : "b", outSrc->seq);
     479
     480                        FILE *f = fopen(fn, "w");
     481                        fprintf(f, "#major minor chisq i\n");
     482                        for (int k = 0; k < chisq->n; k++) {
     483                            fprintf(f, "%10.6f %10.6f %10.6f %5d\n",
     484                                major->data.F64[k]/R_MAJOR, minor->data.F64[k]/R_MINOR, chisq->data.F64[k], k);
     485                        }
     486                        fclose(f);
     487#endif
     488                    } else {
     489                        outSrc->modelEXT->flags |= PM_MODEL_STATUS_NONCONVERGE;
     490
     491#ifdef PRINTVALS
     492                        fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f bad fit\n",
     493                            outSrc->seq, min_k, major->n, R_MAJOR, fRMajorBest, R_MINOR, fRMinorBest);
     494#endif
     495                        psErrorClear();
     496                        // psFree(outSrc->galaxyFits);
     497
     498                        // XXX drop the model since it is bogus? Perhaps we do this or use the nominal values
     499                        // and add a flag
     500#ifdef notnow
     501                        psFree(outSrc->modelEXT);
     502                        outSrc->modelEXT = 0;
     503                        psFree(outSrc->modelFits);
     504                        outSrc->modelFits = 0;
     505#endif
     506                    }
     507                    psFree(chisqFit);
     508                    psFree(major);
     509                    psFree(minor);
     510                    psFree(chisq);
     511                    psFree(flux);
     512                    psFree(dFlux);
     513                } else {
     514                    // No fit
     515                    // Set the values based on the trial with the best chisq
     516                    axes.major = R_MAJOR * fRMajorBest;
     517                    axes.minor = R_MINOR * fRMinorBest;
     518                    outSrc->modelEXT->chisq = minChisq;
     519                    flux0 = fluxBest;
     520#ifdef PRINTVALS
     521                    fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f skip fit\n",
     522                            outSrc->seq, min_k, 0L, R_MAJOR, fRMajorBest, R_MINOR, fRMinorBest);
     523#endif
    327524                }
    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];
     525                // now save the model parameters in the model structure
     526
    343527                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
     528
     529                outSrc->modelEXT->mag = -2.5 * log10(flux0);
     530                outSrc->modelEXT->magErr = dFlux0 / flux0; // 1 / SN
     531                // XXX: should there be a different flag if we didn't do a fit of the chisq
     532                outSrc->modelEXT->flags |= PM_MODEL_STATUS_FITTED;
    348533            }
    349534        }
    … …  
    367552    psMetadata *md = useAnalysis ? readout->analysis : recipe;
    368553
    369     GETVAL(fRmajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN");
    370     GETVAL(fRmajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX");
    371     GETVAL(fRmajorDel, "GALAXY_SHAPES_FR_MAJOR_DEL");
    372     GETVAL(fRminorMin, "GALAXY_SHAPES_FR_MINOR_MIN");
    373     GETVAL(fRminorMax, "GALAXY_SHAPES_FR_MINOR_MAX");
    374     GETVAL(fRminorDel, "GALAXY_SHAPES_FR_MINOR_DEL");
    375 
    376     opt->numTrials = ceil((opt->fRmajorMax - opt->fRmajorMin + 0.5*opt->fRmajorDel) / opt->fRmajorDel) *
    377                          ceil((opt->fRminorMax - opt->fRminorMin + 0.5*opt->fRminorDel) / opt->fRminorDel) ;
     554    GETVAL(fRMajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN");
     555    GETVAL(fRMajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX");
     556    GETVAL(fRMajorDel, "GALAXY_SHAPES_FR_MAJOR_DEL");
     557    GETVAL(fRMinorMin, "GALAXY_SHAPES_FR_MINOR_MIN");
     558    GETVAL(fRMinorMax, "GALAXY_SHAPES_FR_MINOR_MAX");
     559    GETVAL(fRMinorDel, "GALAXY_SHAPES_FR_MINOR_DEL");
     560
     561    opt->numTrials = ceil((opt->fRMajorMax - opt->fRMajorMin + 0.5*opt->fRMajorDel) / opt->fRMajorDel) *
     562                         ceil((opt->fRMinorMax - opt->fRMinorMin + 0.5*opt->fRMinorDel) / opt->fRMinorDel) ;
    378563
    379564    if (makeVectors) {
    380         opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    381         opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     565        opt->fRMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
     566        opt->fRMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32);
    382567        int i = 0;
    383         for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel;
    384                 fRmajor += opt->fRmajorDel) {
    385             for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel;
    386                     fRminor += opt->fRminorDel) {
    387                 opt->rMinor->data.F32[i] = fRminor;
    388                 opt->rMajor->data.F32[i] = fRmajor;
     568        for (float fRMajor = opt->fRMajorMin; fRMajor < opt->fRMajorMax + 0.5*opt->fRMajorDel;
     569                fRMajor += opt->fRMajorDel) {
     570            for (float fRMinor = opt->fRMinorMin; fRMinor < opt->fRMinorMax + 0.5*opt->fRMinorDel;
     571                    fRMinor += opt->fRMinorDel) {
     572                opt->fRMinor->data.F32[i] = fRMinor;
     573                opt->fRMajor->data.F32[i] = fRMajor;
    389574                i++;
    390575            }
    … …  
    392577        psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials);
    393578    } else {
    394         opt->rMinor = NULL;
    395         opt->rMajor = NULL;
     579        opt->fRMinor = NULL;
     580        opt->fRMajor = NULL;
    396581    }
    397582       
    … …  
    413598    }
    414599    CHECKVAL(options, thisReadoutsOptions, numTrials, "mismatched number of trials")
    415     CHECKVAL(options, thisReadoutsOptions, fRmajorMin, "mismatched fRmajorMin")
    416     CHECKVAL(options, thisReadoutsOptions, fRmajorMax, "mismatched fRmajorMax")
    417     CHECKVAL(options, thisReadoutsOptions, fRmajorDel, "mismatched fRmajorDel")
    418     CHECKVAL(options, thisReadoutsOptions, fRminorMin, "mismatched fRminorMin")
    419     CHECKVAL(options, thisReadoutsOptions, fRminorMax, "mismatched fRminorMax")
    420     CHECKVAL(options, thisReadoutsOptions, fRminorDel, "mismatched fRminorDel")
     600    CHECKVAL(options, thisReadoutsOptions, fRMajorMin, "mismatched fRMajorMin")
     601    CHECKVAL(options, thisReadoutsOptions, fRMajorMax, "mismatched fRMajorMax")
     602    CHECKVAL(options, thisReadoutsOptions, fRMajorDel, "mismatched fRMajorDel")
     603    CHECKVAL(options, thisReadoutsOptions, fRMinorMin, "mismatched fRMinorMin")
     604    CHECKVAL(options, thisReadoutsOptions, fRMinorMax, "mismatched fRMinorMax")
     605    CHECKVAL(options, thisReadoutsOptions, fRMinorDel, "mismatched fRMinorDel")
    421606
    422607    return true;
Note: See TracChangeset for help on using the changeset viewer.