IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20864


Ignore:
Timestamp:
Dec 1, 2008, 2:38:47 PM (17 years ago)
Author:
Paul Price
Message:

Removing debugging statement from production.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotApResid.c

    r20581 r20864  
    2222    bool measureAptrend = psMetadataLookupBool (&status, recipe, "MEASURE.APTREND");
    2323    if (!measureAptrend) {
    24         // save nan values since these were not calculated
    25         psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
    26         psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
    27         psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    28         psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
    29         psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    30         psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    31         return true;
     24        // save nan values since these were not calculated
     25        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
     26        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
     27        psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
     28        psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     29        psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
     30        psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
     31        return true;
    3232    }
    3333
    … …  
    9393        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
    9494            Nskip ++;
    95             psTrace ("psphot", 3, "skip : bad source mag");
     95            psTrace ("psphot", 3, "skip : bad source mag");
    9696            continue;
    9797        }
    … …  
    9999        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    100100            Nfail ++;
    101             psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
     101            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
    102102            continue;
    103103        }
    104104
    105         // aperture residual for this source
    106         float dap = source->apMag - source->psfMag;
     105        // aperture residual for this source
     106        float dap = source->apMag - source->psfMag;
    107107
    108108        // sanity clipping : if the model is very discrepant, put your expectation in the recipe
    109109        if ((MAX_AP_OFFSET > 0) && (fabs(dap) > MAX_AP_OFFSET)) {
    110110            Nfail ++;
    111             psTrace ("psphot", 3, "fail : bad dap %f %f", dap, MAX_AP_OFFSET);
     111            psTrace ("psphot", 3, "fail : bad dap %f %f", dap, MAX_AP_OFFSET);
    112112            continue;
    113113        }
    114114
    115         mag->data.F32[Npsf]     = source->psfMag;
     115        mag->data.F32[Npsf]     = source->psfMag;
    116116        apResid->data.F32[Npsf] = dap;
    117117        xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
    … …  
    130130        Npsf ++;
    131131    }
    132    
     132
    133133    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n",
    134134              Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
    … …  
    137137    if (Npsf < APTREND_NSTAR_MIN) {
    138138        psWarning("Only %d valid aperture residual sources (need %d), giving up", Npsf, APTREND_NSTAR_MIN);
    139         // save nan values since these were not calculated
    140         psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
    141         psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
    142         psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    143         psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
    144         psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    145         psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    146 
    147         psFree (mag);
    148         psFree (mask);
    149         psFree (xPos);
    150         psFree (yPos);
    151         psFree (apResid);
    152         psFree (dMag);
     139        // save nan values since these were not calculated
     140        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
     141        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
     142        psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
     143        psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     144        psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
     145        psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
     146
     147        psFree (mag);
     148        psFree (mask);
     149        psFree (xPos);
     150        psFree (yPos);
     151        psFree (apResid);
     152        psFree (dMag);
    153153        return false;
    154154    }
    … …  
    168168    for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
    169169
    170         if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
    171             break;
    172         }
    173 
    174         // store the resulting errorFloor values and the scales, redo the best
    175         if (errorFloor < errorFloorMin) {
    176             errorFloorMin = errorFloor;
    177             entryMin = i;
    178         }
     170        if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
     171            break;
     172        }
     173
     174        // store the resulting errorFloor values and the scales, redo the best
     175        if (errorFloor < errorFloorMin) {
     176            errorFloorMin = errorFloor;
     177            entryMin = i;
     178        }
    179179    }
    180180    if (entryMin == -1) {
    181         psAbort ("failed on ApResid Trend");
    182     }
    183 
    184     // XXX catch error condition 
     181        psAbort ("failed on ApResid Trend");
     182    }
     183
     184    // XXX catch error condition
    185185    psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
    186186
    … …  
    191191
    192192    if (psTraceGetLevel("psphot") >= 2) {
    193         FILE *dumpFile = fopen ("apresid.dat", "w");
    194         for (int i = 0; i < xPos->n; i++) {
    195             fprintf (dumpFile, "%f %f  %f %f %f  %f %f  %x\n",
    196                      xPos->data.F32[i], yPos->data.F32[i],
    197                      mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
    198                      apResid->data.F32[i], apResidRes->data.F32[i],
    199                      mask->data.U8[i]);
    200         }
    201         fclose (dumpFile);
     193        FILE *dumpFile = fopen ("apresid.dat", "w");
     194        for (int i = 0; i < xPos->n; i++) {
     195            fprintf (dumpFile, "%f %f  %f %f %f  %f %f  %x\n",
     196                     xPos->data.F32[i], yPos->data.F32[i],
     197                     mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
     198                     apResid->data.F32[i], apResidRes->data.F32[i],
     199                     mask->data.U8[i]);
     200        }
     201        fclose (dumpFile);
    202202    }
    203203
    … …  
    252252    int nBin = PS_MAX (dMag->n / nGroup, 1);
    253253
    254     // output vectors for ApResid trend 
     254    // output vectors for ApResid trend
    255255    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
    256256    psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
    … …  
    267267    int n = 0;
    268268    for (int i = 0; i < dMo->n; i++) {
    269         int j;
    270         for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
    271             int N = index->data.U32[n];
    272             dMSubset->data.F32[j] = dMag->data.F32[N];
    273             dASubset->data.F32[j] = dap->data.F32[N];
    274             mkSubset->data.U8[j]  = mask->data.U8[N];
    275         }
    276         dMSubset->n = j;
    277         dASubset->n = j;
    278         mkSubset->n = j;
    279 
    280         psStatsInit (statsS);
    281         psStatsInit (statsM);
    282 
    283         if (j > 2) {
    284             bool status = true;
    285             status &= psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff);
    286             status &= psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff);
    287             if (!status) { psErrorClear (); }
    288             dSo->data.F32[i] = statsS->robustStdev;
    289             dMo->data.F32[i] = statsM->sampleMean;
    290             dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
    291             fprintf (stderr, "%d (%d) : sys: %f, phot: %f, rat: %f\n", i, j, dSo->data.F32[i], dMo->data.F32[i], dRo->data.F32[i]);
    292         } else {
    293             dSo->data.F32[i] = NAN;
    294             dMo->data.F32[i] = NAN;
    295             dRo->data.F32[i] = NAN;
    296         }
     269        int j;
     270        for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
     271            int N = index->data.U32[n];
     272            dMSubset->data.F32[j] = dMag->data.F32[N];
     273            dASubset->data.F32[j] = dap->data.F32[N];
     274            mkSubset->data.U8[j]  = mask->data.U8[N];
     275        }
     276        dMSubset->n = j;
     277        dASubset->n = j;
     278        mkSubset->n = j;
     279
     280        psStatsInit (statsS);
     281        psStatsInit (statsM);
     282
     283        if (j > 2) {
     284            bool status = true;
     285            status &= psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff);
     286            status &= psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff);
     287            if (!status) { psErrorClear (); }
     288            dSo->data.F32[i] = statsS->robustStdev;
     289            dMo->data.F32[i] = statsM->sampleMean;
     290            dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
     291            //      fprintf (stderr, "%d (%d) : sys: %f, phot: %f, rat: %f\n", i, j, dSo->data.F32[i], dMo->data.F32[i], dRo->data.F32[i]);
     292        } else {
     293            dSo->data.F32[i] = NAN;
     294            dMo->data.F32[i] = NAN;
     295            dRo->data.F32[i] = NAN;
     296        }
    297297    }
    298298    psFree (dMSubset);
    299299    psFree (dASubset);
    300300    psFree (mkSubset);
    301    
     301
    302302    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    303303    if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
    304         // XXX better testing of raised error
    305         psErrorClear();
     304        // XXX better testing of raised error
     305        psErrorClear();
    306306    }
    307307
    308308    *errorScale = stats->sampleMedian;
    309309    for (int i = 0; i < dSo->n; i++) {
    310         *errorFloor = dSo->data.F32[i];
    311         if (isfinite(*errorFloor)) break;
     310        *errorFloor = dSo->data.F32[i];
     311        if (isfinite(*errorFloor)) break;
    312312    }
    313313
    … …  
    318318    psFree (dMo);
    319319    psFree (dSo);
    320    
     320
    321321    psFree (statsS);
    322322    psFree (statsM);
    … …  
    328328
    329329    int Nx, Ny;
    330        
     330
    331331    if (readout->image->numCols > readout->image->numRows) {
    332332        Nx = scale;
    333         float AR = readout->image->numRows / (float) readout->image->numCols;
    334         Ny = (int) (Nx * AR + 0.5);
     333        float AR = readout->image->numRows / (float) readout->image->numCols;
     334        Ny = (int) (Nx * AR + 0.5);
    335335        Ny = PS_MAX (1, Ny);
    336336    } else {
    337         Ny = scale;
    338         float AR = readout->image->numRows / (float) readout->image->numCols;
    339         Nx = (int) (Ny * AR + 0.5);
    340         Nx = PS_MAX (1, Nx);
    341     }
    342    
     337        Ny = scale;
     338        float AR = readout->image->numRows / (float) readout->image->numCols;
     339        Nx = (int) (Ny * AR + 0.5);
     340        Nx = PS_MAX (1, Nx);
     341    }
     342
    343343    // require at least 10 stars per spatial bin
    344344    if (Npsf < 10*Nx*Ny) {
    345         return false;
     345        return false;
    346346    }
    347347
    … …  
    361361    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
    362362    for (int i = 0; i < dMag->n; i++) {
    363         dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
     363        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
    364364    }
    365365
    366366    // XXX test for errors here
    367367    pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
    368    
     368
    369369    // construct the fitted values and the residuals
    370370    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
    … …  
    387387    return true;
    388388}
    389    
     389
Note: See TracChangeset for help on using the changeset viewer.