IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42705


Ignore:
Timestamp:
Jun 27, 2024, 10:28:50 AM (2 years ago)
Author:
eugene
Message:

clean method for counting rejections selected in bcatalog, add option to record rejections for test images independently; do not set MaxIteration to 0 for low iterations if we are doing Stack/Warp calibrations; in fit1d_irls/FitDataSetSoften pass dyVector to liststats so both get sorted; only call the FitDataSetSoften once; use 0.5 not 0.1 of altsigma in softening the errors

Location:
branches/eam_branches/ipp-20230313/Ohana/src/relphot/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/ImageOps.c

    r41664 r42705  
    531531
    532532  // until the analysis has converged a bit, do not use the IRLS analysis
    533   // default is MaxIterations = 10
    534   if (UseStandardOLS(ZPT_IMAGES)) {
     533  // default is MaxIterations = 10.
     534
     535  // However, for stack and warp calibration, we are not iteratively solving for the
     536  // zero points.  In this case, use IRLS in a single pass
     537  if (!CALIBRATE_STACKS_AND_WARPS && UseStandardOLS(ZPT_IMAGES)) {
    535538    brightStars.MaxIterations = 0;
    536539    kronStars.MaxIterations = 0;
    … …  
    672675
    673676      if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
    674         fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysKron, MrelKron, Mmos, Mgrid, Mflat, kronStars.alldata->yVector[Nkron]);
     677        fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f %6.3f %6.3f | %5.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, MsysKron, MrelKron, Mmos, Mgrid, Mflat, Moff, psfStars.alldata->dyVector[Nref]);
    675678      }
    676679
    … …  
    702705    if (mark) continue;
    703706
    704     // use liststats to find the 20-pct, median, 80-pct points
    705     StatType stats;
    706     liststats_setmode (&stats, "MEDIAN");
    707     liststats (psfStars.alldata->yVector, NULL, NULL, Nref, &stats);
    708     double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
    709 
    710     // soften the individual errors with 10% of the scatter above
    711     for (j = 0; j < Nref; j++) {
    712       double newSigma = hypot(psfStars.alldata->dyVector[j], 0.1*altSigma);
    713       psfStars.alldata->dyVector[j] = newSigma;
    714     }
    715 
    716     // soften the errors based on the scatter
     707    // soften the errors based on the core scatter
    717708    FitDataSetSoften (&psfStars, Nref);
    718709

  • branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/bcatalog.c

    r42696 r42705  
    11# include "relphot.h"
     2
     3# define BAD_COUNT(FIELD) \
     4  badCount.FIELD ++; \
     5  if (doTest1 && (catalog[0].measure[offset].imageID == TEST_IMAGE1)) { badTest1.FIELD ++; } \
     6  if (doTest2 && (catalog[0].measure[offset].imageID == TEST_IMAGE2)) { badTest2.FIELD ++; }
     7
     8typedef struct {
     9  int Ncode;
     10  int Ntime;
     11  int Ndophot;
     12  int Nmag;
     13  int Nsigma;
     14  int Nimag;
     15  int Nfew;
     16  int Ngalaxy;
     17  int Npsfqf;
     18  int Nnan;
     19  int Nbad;
     20  int Npoor;
     21} MyBadCountType;
    222
    323int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
    … …  
    929  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    1030  float mag;
    11   int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
     31
     32  // int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
    1233
    1334  int Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    2647  int NmMin = catalog[0].Nmeasure;
    2748
    28   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
     49  MyBadCountType badCount = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
     50  MyBadCountType badTest1 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
     51  MyBadCountType badTest2 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
     52
     53  int doTest1 = (TEST_IMAGE1 >= 0);
     54  int doTest2 = (TEST_IMAGE2 >= 0);
     55
     56  FILE *ftest1 = (doTest1) ? fopen ("test1.meas.dat", "a") : NULL;
     57  FILE *ftest2 = (doTest2) ? fopen ("test2.meas.dat", "a") : NULL;
     58
     59  // Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
    2960
    3061  // copy the following fields to the subcatalog:
    … …  
    5788
    5889      offset = catalog[0].average[i].measureOffset + j;
     90
     91      if (doTest1) {
     92        Measure *m = &catalog[0].measure[offset];
     93        if (m->imageID == TEST_IMAGE1) {
     94          fprintf (ftest1, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
     95        }
     96      }
     97      if (doTest2) {
     98        Measure *m = &catalog[0].measure[offset];
     99        if (m->imageID == TEST_IMAGE2) {
     100          fprintf (ftest2, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
     101        }
     102      }
    59103
    60104      /* select measurements by photcode */
    … …  
    69113      }
    70114      if (!found) {
    71         Ncode ++;
     115        BAD_COUNT(Ncode);
    72116        continue;
    73117      }
    … …  
    75119      /* select measurements by time */
    76120      if (TimeSelect) {
    77         if (catalog[0].measure[offset].t < TSTART) { Ntime ++; continue; }
    78         if (catalog[0].measure[offset].t > TSTOP)  { Ntime ++; continue; }
     121        if (catalog[0].measure[offset].t < TSTART) { BAD_COUNT(Ntime); continue; }
     122        if (catalog[0].measure[offset].t > TSTOP)  { BAD_COUNT(Ntime); continue; }
    79123      }
    80124
    81125      /* select measurements by quality */
    82       if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; }
     126      if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { BAD_COUNT(Ndophot); continue; }
    83127
    84128      float Maper = USE_APER_FOR_STARGAL ? catalog[0].measure[offset].Map : catalog[0].measure[offset].Mkron;
    85129
    86130      // skip garbage measurements
    87       if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF     < PSF_QF_MIN)) { Npsfqf ++; continue; }
    88       if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < PSF_QF_MIN)) { Npsfqf ++; continue; }
    89       if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
    90       if (isnan(Maper)) { Nnan ++; continue; }
     131      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF     < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
     132      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
     133      if (isnan(catalog[0].measure[offset].M)) { BAD_COUNT(Nnan); continue; }
     134      if (isnan(Maper)) { BAD_COUNT(Nnan); continue; }
    91135
    92136      // require 0x01 in photFlags (fitted with a PSF) -- add to the test data
    93       if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
     137      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { BAD_COUNT(Nbad); continue; }
    94138
    95139      // very loose cut on PSF - APER mag (Map or Mkron) -- this lightly rejects galaxies & other oddities
    96140      float Mkp = catalog[0].measure[offset].M - Maper;
    97       if (fabs(Mkp) > 1.0) { Nbad ++; continue; }
     141      if (fabs(Mkp) > 1.0) { BAD_COUNT(Nbad); continue; }
    98142
    99143      if (catalog[0].measure[offset].photFlags & code->photomBadMask) {
    100         Nbad++;
     144        BAD_COUNT(Nbad);
    101145        continue;
    102146      }
    103       if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
     147      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { BAD_COUNT(Npoor); continue;}
    104148
    105149      // weak test for galaxies
    … …  
    112156      /* select measurements by mag limit */
    113157      mag = PhotCat (&catalog[0].measure[offset], MAG_CLASS_PSF);
    114       if (mag > MAG_LIM) { Nmag ++; continue; }
    115       if (mag < MAG_MIN) { Nmag ++; continue; }
     158      if (mag > MAG_LIM) { BAD_COUNT(Nmag); continue; }
     159      if (mag < MAG_MIN) { BAD_COUNT(Nmag); continue; }
    116160      // note that these two cuts are applied to all photcodes. these are also applied to
    117161      // the nominal apparent mags : if the expected sensitivity is very different between
    … …  
    120164
    121165      /* select measurements by measurement error */
    122       if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { Nsigma ++; continue; }
     166      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { BAD_COUNT(Nsigma); continue; }
    123167
    124168      /* select measurements by mag limit */
    125169      if (ImagSelect) {
    126170        mag = PhotInst (&catalog[0].measure[offset], MAG_CLASS_PSF);
    127         if (mag < ImagMin) { Nimag ++; continue; }
    128         if (mag > ImagMax) { Nimag ++; continue; }
     171        if (mag < ImagMin) { BAD_COUNT(Nimag); continue; }
     172        if (mag > ImagMax) { BAD_COUNT(Nimag); continue; }
    129173      }
    130174      // note that instrumental mag cuts are reasonable for individual exposures using the same
    … …  
    155199    if (nEXT >= nPSF) {
    156200      Nmeasure -= Nm;
    157       Ngalaxy ++;
     201      BAD_COUNT(Ngalaxy);
    158202      continue;
    159203    }
    … …  
    165209    if (Nm <= STAR_TOOFEW) { /* enough measurements total? */
    166210      Nmeasure -= Nm;
    167       Nfew ++;
     211      BAD_COUNT(Nfew);
    168212      continue;
    169213    }
    … …  
    175219    if (!anySecfiltGood) {
    176220      Nmeasure -= Nm;
    177       Nfew ++;
     221      BAD_COUNT(Nfew);
    178222      continue;
    179223    }
    … …  
    212256             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax);
    213257    fprintf (stderr, "rejections: %d stars have too few measures:\n   %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
    214              Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
     258             badCount.Nfew, badCount.Ncode, badCount.Ntime, badCount.Ndophot, badCount.Nmag, badCount.Nsigma, badCount.Nimag, badCount.Npsfqf, badCount.Nnan, badCount.Ngalaxy, badCount.Nbad, badCount.Npoor);
     259  }
     260
     261  if (doTest1) {
     262    fprintf (stderr, "TEST1 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     263             badTest1.Nfew, badTest1.Ncode, badTest1.Ntime, badTest1.Ndophot, badTest1.Nmag, badTest1.Nsigma, badTest1.Nimag, badTest1.Npsfqf, badTest1.Nnan, badTest1.Ngalaxy, badTest1.Nbad, badTest1.Npoor);
     264  }
     265  if (doTest2) {
     266    fprintf (stderr, "TEST2 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     267             badTest2.Nfew, badTest2.Ncode, badTest2.Ntime, badTest2.Ndophot, badTest2.Nmag, badTest2.Nsigma, badTest2.Nimag, badTest2.Npsfqf, badTest2.Nnan, badTest2.Ngalaxy, badTest2.Nbad, badTest2.Npoor);
    215268  }
    216269
    … …  
    220273  }
    221274
    222   free (Nvalid)
     275  free (Nvalid);
     276
     277  if (ftest1) { fclose (ftest1); }
     278  if (ftest2) { fclose (ftest2); }
    223279
    224280  return (TRUE);

  • branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/fit1d_irls.c

    r41557 r42705  
    374374  StatType stats;
    375375  liststats_setmode (&stats, "MEDIAN");
    376   liststats (dataset->alldata->yVector, NULL, NULL, Nvalues, &stats);
     376  liststats (dataset->alldata->yVector, dataset->alldata->dyVector, NULL, Nvalues, &stats);
     377  // liststats sorts yVector and dyVector together
    377378
    378379  double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
    … …  
    380381  // soften the individual errors with 10% of the scatter above
    381382  for (int j = 0; j < Nvalues; j++) {
    382     double newSigma = hypot(dataset->alldata->dyVector[j], 0.1*altSigma);
     383    double newSigma = hypot(dataset->alldata->dyVector[j], 0.5*altSigma);
    383384    dataset->alldata->dyVector[j] = newSigma;
    384385  }
Note: See TracChangeset for help on using the changeset viewer.