Changeset 41454 for trunk

Timestamp:

Dec 7, 2020, 5:45:21 PM (6 years ago)

Author:

eugene

Message:

pre-allocate psf,kron,brightStar arrays; use FitStatData for MosaicOps.c

Location:

trunk/Ohana/src/relphot/src

Files:

• MosaicOps.c (modified) (6 diffs)
• TGroupOps.c (modified) (8 diffs)

trunk/Ohana/src/relphot/src/MosaicOps.c

@@ -925 +925 @@
   off_t Nmax;
   int PoorImages;
-  double *psfMagList;
-  double *psfErrList;
-  double *brightMagList;
-  double *brightErrList;
-  double *kronMagList;
-  double *kronErrList;
+  FitDataSet psfStars;
+  FitDataSet kronStars;
+  FitDataSet brightStars;
 } SetMmosInfo;
 
@@ -961 +958 @@
 
   if (allocLists) {
-    ALLOCATE (info->psfMagList, double, Nmax);
-    ALLOCATE (info->psfErrList, double, Nmax);
-    ALLOCATE (info->kronMagList, double, Nmax);
-    ALLOCATE (info->kronErrList, double, Nmax);
-    ALLOCATE (info->brightMagList, double, Nmax);
-    ALLOCATE (info->brightErrList, double, Nmax);
+    // we can only fit the zero point for individual exposures, not the airmass
+    FitDataSetAlloc (&info->psfStars,    Nmax, 0, 0); 
+    FitDataSetAlloc (&info->kronStars,   Nmax, 0, 0); 
+    FitDataSetAlloc (&info->brightStars, Nmax, 0, 0); 
   }
 }
 
 void SetMmosInfoFree (SetMmosInfo *info) {
-  free (info->psfMagList   );
-  free (info->psfErrList   );
-  free (info->kronMagList  );
-  free (info->kronErrList  );
-  free (info->brightMagList);
-  free (info->brightErrList);
+  FitDataSetFree (&info->brightStars);
+  FitDataSetFree (&info->kronStars);
+  FitDataSetFree (&info->psfStars);
 }
 
@@ -1074 +1066 @@
   off_t j, NimageReal;
 
-  StatType stats;
-  liststats_setmode (&stats, "INNER_WTMEAN");
-
-  double *psfMagList    = info->psfMagList;
-  double *psfErrList    = info->psfErrList;
-  double *kronMagList   = info->kronMagList;
-  double *kronErrList   = info->kronErrList;
-  double *brightMagList = info->brightMagList;
-  double *brightErrList = info->brightErrList;
-
-  assert (Nmos >= 0);
-  assert (Nmos < Nmosaic);
-
-  // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
-  getimages (&NimageReal, NULL);
-
   /* on PoorImages run, skip good images */
   if (info->PoorImages) {
@@ -1095 +1071 @@
     if (!bad) return TRUE;
   }      
+
+  FitDataSet *psfStars = &info->psfStars;
+  FitDataSet *kronStars = &info->kronStars;
+  FitDataSet *brightStars = &info->brightStars;
+
+  assert (Nmos >= 0);
+  assert (Nmos < Nmosaic);
+
+  // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
+  getimages (&NimageReal, NULL);
 
   // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
@@ -1230 +1216 @@
     assert (Nbright >= 0);
 
-    psfMagList[N]  = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
-    psfErrList[N]  = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-    kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
-    kronErrList[N] = psfErrList[N];
+    psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    psfStars->alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+
+    kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
+    kronStars->alldata->dyVector[N] = psfStars->alldata->dyVector[N];
+
     if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-      brightMagList[Nbright] = psfMagList[N];
-      brightErrList[Nbright] = psfErrList[N];
+      brightStars->alldata-> yVector[Nbright] = psfStars->alldata-> yVector[N];
+      brightStars->alldata->dyVector[Nbright] = psfStars->alldata->dyVector[N];
       Nbright ++;
     }
     N++;
+
   }
   /* N_onMosaic[Nmos] is all measurements, N is good measurements */
 
-  if (fout) {
-    fclose (fout);
-  }
+  if (fout) { fclose (fout); }
 
   /* too few good measurements or too many bad measurements (skip in PoorImages run) */
@@ -1262 +1249 @@
   }
 
-  liststats (psfMagList, psfErrList, NULL, N, &stats);
-  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
+  fit1d_irls (psfStars, N);
+  // if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
 
   // for now, I have no reason to measure these separately for camera-level images
-  myMosaic[0].McalPSF    = stats.mean;
-  myMosaic[0].dMcal      = stats.error;
-  myMosaic[0].McalChiSq  = stats.chisq;
-  myMosaic[0].nFitPhotom = N;
-
-  liststats (kronMagList, kronErrList, NULL, N, &stats);
-  myMosaic[0].McalAPER   = stats.mean;
+  myMosaic[0].McalPSF    = psfStars->bSaveArray[0][0];
+  myMosaic[0].dMcal      = psfStars->bSigma[0];     
+  myMosaic[0].McalChiSq  = psfStars->chisq;         
+  myMosaic[0].nFitPhotom = psfStars->Nmeas;         
+
+  // fit1d_irls (kronStars, N);    // does it make sense to calculate a zero point based kron mags?
+  // myMosaic->McalAPER   = ???
 
   if (testImage) {
-    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic[0].start, stats.mean, stats.error, myMosaic[0].nFitPhotom);
+    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic->start, myMosaic->McalPSF, myMosaic->dMcal, myMosaic->nFitPhotom);
   }
 
   if (PLOTSTUFF) {
-    fprintf (stderr, "Mmos: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[MosaicToImage[Nmos][0]].name);
-    plot_setMcal (psfMagList, N);
+    fprintf (stderr, "Mmos: %6.3f +/- %6.3f %5d %5d | %s\n", myMosaic->McalPSF, myMosaic->dMcal, myMosaic->nFitPhotom, N, image[MosaicToImage[Nmos][0]].name);
+    plot_setMcal (psfStars->alldata->yVector, N);
   }
 
   // bright end scatter
-  liststats (brightMagList, brightErrList, NULL, Nbright, &stats);
-  myMosaic[0].dMsys = stats.sigma;
+  // fit1d_irls (brightStars, N);  // does it make sense to calculate a separate bright star zero point?
+  // myMosaic->dMsys = ??? 
 
   if (myMosaic[0].McalPSF < -CLOUD_TOLERANCE) {

trunk/Ohana/src/relphot/src/TGroupOps.c

@@ -527 +527 @@
   off_t Nmax;
   int PoorImages;
-  double *psfMagList;
-  double *psfAirList;
-  double *psfErrList;
-  double *brightMagList;
-  double *brightErrList;
-  double *kronMagList;
-  double *kronErrList;
+  FitDataSet psfStars;
+  FitDataSet kronStars;
+  FitDataSet brightStars;
 } SetMgrpInfo;
 
@@ -565 +561 @@
 
   if (allocLists) {
-    ALLOCATE (info->psfMagList, double, Nmax);
-    ALLOCATE (info->psfAirList, double, Nmax);
-    ALLOCATE (info->psfErrList, double, Nmax);
-    ALLOCATE (info->kronMagList, double, Nmax);
-    ALLOCATE (info->kronErrList, double, Nmax);
-    ALLOCATE (info->brightMagList, double, Nmax);
-    ALLOCATE (info->brightErrList, double, Nmax);
+    int fitOrder = TGROUP_FIT_AIRMASS ? 1 : 0;
+    
+    FitDataSetAlloc (&info->psfStars,    Nmax, fitOrder, 0); 
+    FitDataSetAlloc (&info->kronStars,   Nmax, fitOrder, 0); 
+    FitDataSetAlloc (&info->brightStars, Nmax, fitOrder, 0); 
+
+    if (1) {
+      FitDataSetAddPriors (&info->psfStars);
+      info->psfStars.bPriorValue[1] = 0.0;  // note that we are fitting relative to the nominal slope
+      info->psfStars.bPriorSigma[1] = 0.04; // XXX this prior sigma needs to be user-configured
+    }
   }
 }
 
 void SetMgrpInfoFree (SetMgrpInfo *info) {
-  free (info->psfMagList   );
-  free (info->psfAirList   );
-  free (info->psfErrList   );
-  free (info->kronMagList  );
-  free (info->kronErrList  );
-  free (info->brightMagList);
-  free (info->brightErrList);
+  FitDataSetFree (&info->brightStars);
+  FitDataSetFree (&info->kronStars);
+  FitDataSetFree (&info->psfStars);
 }
 
@@ -694 +690 @@
   }      
 
-  int fitOrder = TGROUP_FIT_AIRMASS ? 1 : 0;
-
-  // allocate the arrays used by the fit1d_irls analysis below.
-  // the last 0 means we are NOT using bootstrap analysis here
-
-  FitDataSet psfStars, kronStars, brightStars;
-  FitDataSetAlloc (&psfStars,    myTGroup->Nmeasure, fitOrder, 0); 
-  FitDataSetAlloc (&kronStars,   myTGroup->Nmeasure, fitOrder, 0); 
-  FitDataSetAlloc (&brightStars, myTGroup->Nmeasure, fitOrder, 0); 
-
-  // add prior constraint on airmass:
-  if (1) {
-    FitDataSetAddPriors (&psfStars);
-    psfStars.bPriorValue[1] = 0.0;  // note that we are fitting relative to the nominal slope
-    psfStars.bPriorSigma[1] = 0.04; // XXX this prior sigma needs to be user-configured
-  }
+  FitDataSet *psfStars = &info->psfStars;
+  FitDataSet *kronStars = &info->kronStars;
+  FitDataSet *brightStars = &info->brightStars;
+
+  assert (Ngrp >= 0);
+  assert (Ngrp < NtgroupTimes);
 
   // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
@@ -820 +806 @@
     assert (Nbright >= 0);
 
-    psfStars.alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
-    psfStars.alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0);
-    psfStars.alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-
-    kronStars.alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
-    kronStars.alldata-> xVector[N] = psfStars.alldata-> xVector[N];
-    kronStars.alldata->dyVector[N] = psfStars.alldata->dyVector[N];
+    psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    psfStars->alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0);
+    psfStars->alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+
+    kronStars->alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
+    kronStars->alldata-> xVector[N] = psfStars->alldata-> xVector[N];
+    kronStars->alldata->dyVector[N] = psfStars->alldata->dyVector[N];
 
     if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-      brightStars.alldata-> yVector[Nbright] = psfStars.alldata-> yVector[N];
-      brightStars.alldata-> xVector[Nbright] = psfStars.alldata-> xVector[N];
-      brightStars.alldata->dyVector[Nbright] = psfStars.alldata->dyVector[N];
+      brightStars->alldata-> yVector[Nbright] = psfStars->alldata-> yVector[N];
+      brightStars->alldata-> xVector[Nbright] = psfStars->alldata-> xVector[N];
+      brightStars->alldata->dyVector[Nbright] = psfStars->alldata->dyVector[N];
       Nbright ++;
     }
@@ -859 +845 @@
   // (myTGroup->flags & ID_IMAGE_PHOTOM_UBERCAL)
 
-  fit1d_irls (&psfStars, N);
+  fit1d_irls (psfStars, N);
 
   // if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
 
   // for now, I have no reason to measure these separately for camera-level images
-  myTGroup->McalPSF    = psfStars.bSaveArray[0][0];
-  myTGroup->dMcal      = psfStars.bSigma[0];
-  myTGroup->McalChiSq  = psfStars.chisq;
-  myTGroup->nFitPhotom = psfStars.Nmeas;
-  myTGroup->dKlam      = psfStars.bSaveArray[1][0];
+  myTGroup->McalPSF    = psfStars->bSaveArray[0][0];
+  myTGroup->dMcal      = psfStars->bSigma[0];
+  myTGroup->McalChiSq  = psfStars->chisq;
+  myTGroup->nFitPhotom = psfStars->Nmeas;
+  myTGroup->dKlam      = psfStars->bSaveArray[1][0];
 
   // XXX drop this for tgroups?
-  // fit1d_irls (&kronStars, N);    // does it make sense to calculate a zero point based kron mags?
+  // fit1d_irls (kronStars, N);    // does it make sense to calculate a zero point based kron mags?
   // myTGroup->McalAPER   = stats.mean;
 
@@ -881 +867 @@
   if (PLOTSTUFF) {
     fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, N, ohana_sec_to_mjd(tgroup->start));
-    plot_setMcal (psfStars.alldata->yVector, N);
+    plot_setMcal (psfStars->alldata->yVector, N);
   }
 
@@ -887 +873 @@
   // XXX this does not make sense: I need to apply the airmass slope calculated above first
   // redo this by calculating the corrected bright stars mags 
-  // fit1d_irls (&brightStars, N);  // does it make sense to calculate a separate bright star zero point?
+  // fit1d_irls (brightStars, N);  // does it make sense to calculate a separate bright star zero point?
   // myTGroup->dMsys = stats.sigma;
 
@@ -899 +885 @@
   }
 
-  FitDataSetFree (&brightStars);
-  FitDataSetFree (&kronStars);
-  FitDataSetFree (&psfStars);
+  // FitDataSetFree (&brightStars);
+  // FitDataSetFree (&kronStars);
+  // FitDataSetFree (&psfStars);
 
   return TRUE;