IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41473


Ignore:
Timestamp:
Jan 24, 2021, 3:31:48 PM (5 years ago)
Author:
eugene
Message:

allow reference stars to have higher weight; push OLS / IRLS choice into ZeroPointModes; update parallel version

Location:
trunk/Ohana/src/relphot
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/include/relphot.h

    r41467 r41473  
    3030
    3131typedef enum {
     32  ZPT_STARS,
    3233  ZPT_TGROUP,
    3334  ZPT_MOSAIC,
    … …  
    458459// int            PhotNsec;
    459460
    460 PhotCode      *refPhotcode;
     461PhotCode *refPhotcode;
     462int       USE_REFERENCE_WEIGHT;
    461463
    462464int MaxDensityUse;
    … …  
    633635int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    634636
    635 void          setMcal             PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr));
    636 int           setMmos             PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr));
    637 int           setMgrp             PROTO((Catalog *catalog, int Nloop, int Poor, FlatCorrectionTable *flatcorr));
     637void          setMcal             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
     638int           setMmos             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
     639int           setMgrp             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
    638640
    639641int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    … …  
    834836void dump_catalog (Catalog *catalog, off_t c, int Npass);
    835837
    836 void SetZeroPointModes (int nloop);
    837 int UseStandardOLS (int nloop, ZptFitModeType mode);
     838void SetZptIteration (int current);
     839int GetZptIteration (void);
     840
     841void SetZeroPointModes (void);
     842int UseStandardOLS (ZptFitModeType mode);

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r41467 r41473  
    581581
    582582/* determine Mcal values for all images */
    583 void setMcal (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
     583void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    584584
    585585  off_t i, j, m, c, n;
    … …  
    611611  // until the analysis has converged a bit, do not use the IRLS analysis
    612612  // default is MaxIterations = 10
    613   if (UseStandardOLS(Nloop, ZPT_IMAGES)) {
     613  if (UseStandardOLS(ZPT_IMAGES)) {
    614614    brightStars.MaxIterations = 0;
    615615    kronStars.MaxIterations = 0;
    … …  
    799799
    800800    // bright end scatter
    801     // fit1d_irls (&brightStars, Nbright);
    802     // image[i].dMagSys = brightStars.sigma;
     801    fit1d_irls (&brightStars, Nbright);
     802    image[i].dMagSys = brightStars.sigma;
    803803
    804804    if (image[i].McalPSF < -CLOUD_TOLERANCE) {

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

    r41467 r41473  
    306306    mosaic[j].nValPhotom = 0;
    307307
    308     mosaic[j].dMsys     = subset[i].flags;
     308    mosaic[j].dMsys     = subset[i].dMagSys;
    309309    mosaic[j].flags     = subset[i].flags;
    310310    mosaic[j].secz      = subset[i].secz;
    … …  
    446446    mosaic[j].dMcal     = 0.0;
    447447    mosaic[j].McalChiSq = 0.0;
    448     mosaic[j].dMsys     = image[i].flags;
     448    mosaic[j].dMsys     = image[i].dMagSys;
    449449    mosaic[j].flags     = image[i].flags;
    450450    mosaic[j].secz      = image[i].secz;
    … …  
    743743      image[im].flags      |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW);
    744744      image[im].flags      |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR);
     745      image[im].flags      |= (mosaic[i].flags & ID_IMAGE_MOSAIC_POOR);
    745746    }
    746747    mosaic[i].McalPSF  = 0.0;
    … …  
    973974  int Nsys;
    974975  int Nskip;
    975   int Nloop;
    976976  off_t Nmax;
    977977  int PoorImages;
    … …  
    994994int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
    995995void *setMmos_worker (void *data);
    996 int setMmos_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr);
    997 
    998 void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages, int Nloop) {
     996int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
     997
     998void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) {
    999999  info->Nfew = 0;
    10001000  info->Nbad = 0;
    … …  
    10061006
    10071007  info->Nmax = Nmax;
    1008   info->Nloop = Nloop;
    10091008  info->PoorImages = PoorImages;
    10101009
    … …  
    10161015
    10171016    // until the analysis has converged a bit, do not use the IRLS analysis
    1018     if (UseStandardOLS(Nloop, ZPT_MOSAIC)) {
     1017    if (UseStandardOLS(ZPT_MOSAIC)) {
    10191018      info->brightStars.MaxIterations = 0;
    10201019      info->kronStars.MaxIterations = 0;
    … …  
    10621061}
    10631062
    1064 int setMmos (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
     1063int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    10651064
    10661065  off_t i, N, Nmax;
    … …  
    10751074  // so do not run setMmos in threaded mode if PLOTSTUFF is set
    10761075  if (NTHREADS && !PLOTSTUFF) {
    1077     int status = setMmos_threaded (catalog, Nloop, PoorImages, flatcorr);
     1076    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
    10781077    return status;
    10791078  }
    … …  
    10951094
    10961095  SetMmosInfo info;
    1097   SetMmosInfoInit (&info, Nmax, TRUE, PoorImages, Nloop);
     1096  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);
    10981097
    10991098  // int savePlotDelay = PLOTDELAY;
    … …  
    13381337
    13391338  fit1d_irls (psfStars, N);
    1340   // if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
     1339  fit1d_irls (brightStars, Nbright);
    13411340
    13421341  // for now, I have no reason to measure these separately for camera-level images
    … …  
    13471346  myMosaic[0].nValPhotom = N;       
    13481347
    1349   // fit1d_irls (kronStars, N);    // does it make sense to calculate a zero point based kron mags?
    1350   // myMosaic->McalAPER   = ???
     1348  fit1d_irls (kronStars, N);
     1349  myMosaic->McalAPER   = kronStars->bSaveArray[0][0];
    13511350
    13521351  if (testImage) {
    … …  
    13601359
    13611360  // bright end scatter
    1362   // fit1d_irls (brightStars, N);  // does it make sense to calculate a separate bright star zero point?
    1363   // myMosaic->dMsys = ???
     1361  fit1d_irls (brightStars, Nbright);
     1362  myMosaic->dMsys = brightStars->sigma; // stdev of bright star mags
    13641363
    13651364  if (myMosaic[0].McalPSF < -CLOUD_TOLERANCE) {
    … …  
    13811380    off_t im = MosaicToImage[Nmos][j];
    13821381    image[im].ubercalDist = myMosaic[0].ubercalDist;
    1383     // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    13841382  }
    13851383 
    … …  
    13901388}
    13911389 
    1392 int setMmos_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
     1390int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    13931391
    13941392  int i;
    … …  
    14111409
    14121410  SetMmosInfo summary;
    1413   SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages, Nloop);
     1411  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);
    14141412
    14151413  pthread_attr_t attr;
    … …  
    14371435    // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop,
    14381436    // PoorImages) used in the threads to allocate the arrays and set the MaxIterations
    1439     SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages, Nloop);
     1437    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
    14401438    pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
    14411439  }
    … …  
    14951493
    14961494  SetMmosInfo results;
    1497   SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages, threadinfo->info.Nloop);
     1495  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages);
    14981496
    14991497  while (1) {

  • trunk/Ohana/src/relphot/src/StarOps.c

    r41462 r41473  
    10841084void dump_catalog (Catalog *catalog, off_t c, int Npass) {
    10851085
     1086  return;
     1087
    10861088  int Nsecfilt = GetPhotcodeNsecfilt ();
    10871089

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

    r41467 r41473  
    427427        image[im].dMagSys     = tgroup[i].dMsys;
    428428        image[im].nFitPhotom  = tgroup[i].nFitPhotom;
    429         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_NIGHT_POOR); // XXX copy this down?
    430         image[im].flags            |= (tgroup[i].flags & ID_IMAGE_MOSAIC_POOR); // XXX copy this down?
     429        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_NIGHT_POOR); // probably not necessary : it is set in markBadTGroup
     430        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_MOSAIC_POOR); // probably not necessary : it is set in markBadTGroup
    431431        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW);
    432432        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR);
    433433        image[im].flags            |= (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL);
    434434       
    435         fprintf (stderr, "TG to IMAGE: %f +/- %f -> %f (%d)\n", tgroup[i].McalPSF, tgroup[i].dMcal, image[im].McalPSF, tgroup[i].nFitPhotom);
     435        // fprintf (stderr, "TG to IMAGE: %f +/- %f -> %f (%d)\n", tgroup[i].McalPSF, tgroup[i].dMcal, image[im].McalPSF, tgroup[i].nFitPhotom);
    436436      }
    437437      tgroup[i].McalPSF  = 0.0;
    … …  
    547547  }
    548548  fprintf (stderr, "Matched %d detections to tgroups\n", Nmatch);
    549   for (i = 0; i < NtgroupTimes; i++) {
     549
     550  // XXX print results for testing
     551  for (i = 0; FALSE && (i < NtgroupTimes); i++) {
    550552    TGroup *tgroup = tgroupTimes[i][0].byCode;
    551553    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
    552       fprintf (stderr, "  TGROUP time %.0f photcode %d, Nmeasure: %d, Nimage: %d\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, (int) tgroup[j].Nmeasure, (int) tgroup[j].Nimage);
     554      if (0) fprintf (stderr, "  TGROUP time %.0f photcode %d, Nmeasure: %d, Nimage: %d\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, (int) tgroup[j].Nmeasure, (int) tgroup[j].Nimage);
    553555    }
    554556  }
    … …  
    643645  int Nsys;
    644646  int Nskip;
    645   int Nloop;
    646647  off_t Nmax;
    647648  int PoorImages;
    … …  
    664665int   setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr);
    665666void *setMgrp_worker (void *data);
    666 int   setMgrp_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr);
    667 
    668 void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages, int Nloop) {
     667int   setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);
     668
     669void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages) {
    669670  info->Nfew = 0;
    670671  info->Nbad = 0;
    … …  
    677678
    678679  info->Nmax = Nmax;
    679   info->Nloop = Nloop;
    680680  info->PoorImages = PoorImages;
    681681
    … …  
    688688
    689689    // until the analysis has converged a bit, do not use the IRLS analysis
    690     if (UseStandardOLS(Nloop, ZPT_TGROUP)) {
     690    if (UseStandardOLS(ZPT_TGROUP)) {
    691691      info->brightStars.MaxIterations = 0;
    692692      info->kronStars.MaxIterations = 0;
    … …  
    740740}
    741741
    742 int setMgrp (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
     742int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    743743
    744744  off_t N;
    … …  
    751751  // so do not run setMgrp in threaded mode if PLOTSTUFF is set
    752752  if (NTHREADS && !PLOTSTUFF) {
    753     int status = setMgrp_threaded (catalog, Nloop, PoorImages, flatcorr);
     753    int status = setMgrp_threaded (catalog, PoorImages, flatcorr);
    754754    return status;
    755755  }
    … …  
    766766
    767767  SetMgrpInfo info;
    768   SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages, Nloop);
     768  SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages);
    769769
    770770  for (off_t i = 0; i < NtgroupTimes; i++) {
    … …  
    772772    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
    773773      setMgrp_tgroup (&tgroup[j], i, image, catalog, &info, flatcorr);
    774       fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE);
     774      // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE);
    775775    }
    776776  }
    … …  
    971971  myTGroup->dKlam      = psfStars->bSaveArray[1][0];
    972972
    973   // XXX drop this for tgroups?
    974   // fit1d_irls (kronStars, N);    // does it make sense to calculate a zero point based kron mags?
    975   // myTGroup->McalAPER   = stats.mean;
     973  // XXX this does not make sense: I need to apply the airmass slope calculated above first
     974  fit1d_irls (kronStars, N);
     975  myTGroup->McalAPER   = kronStars->bSaveArray[0][0];
    976976
    977977  if (testImage) {
    … …  
    979979    fprintf (stderr, "test night %.0f aper: %f %f %d ... ", ohana_sec_to_mjd(parent->start), myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom);
    980980  }
    981 
    982   // fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.2f %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, myTGroup->nValPhotom, myTGroup->McalChiSq, ohana_sec_to_mjd(tgroup->start));
    983 
    984981  if (PLOTSTUFF) {
    985982    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));
    … …  
    990987  // XXX this does not make sense: I need to apply the airmass slope calculated above first
    991988  // redo this by calculating the corrected bright stars mags
    992   // fit1d_irls (brightStars, N);  // does it make sense to calculate a separate bright star zero point?
    993   // myTGroup->dMsys = stats.sigma;
     989  fit1d_irls (brightStars, Nbright);
     990  myTGroup->dMsys = brightStars->sigma;
    994991
    995992  // keep this??
    … …  
    10051002}
    10061003 
    1007 int setMgrp_threaded (Catalog *catalog, int Nloop, int PoorImages, FlatCorrectionTable *flatcorr) {
     1004int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    10081005
    10091006  int i;
    … …  
    10231020
    10241021  SetMgrpInfo summary;
    1025   SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages, Nloop);
     1022  SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages);
    10261023
    10271024  pthread_attr_t attr;
    … …  
    10471044    threadinfo[i].flatcorr = flatcorr;
    10481045
    1049     // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop,
     1046    // we do NOT allocate the arrays here, we only supply basic info (Nmax,
    10501047    // PoorImages) used in the threads to allocate the arrays and set the MaxIterations
    1051     SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages, Nloop);
     1048    SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
    10521049    pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]);
    10531050  }
    … …  
    11041101
    11051102  SetMgrpInfo results;
    1106   SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages, threadinfo->info.Nloop); // allocate list, dlist arrays here
     1103  SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here
    11071104
    11081105  while (1) {
    … …  
    14271424void dump_tgroups (Catalog *catalog, int Npass) {
    14281425
     1426  return;
     1427
    14291428  int Nsecfilt = GetPhotcodeNsecfilt ();
    14301429

  • trunk/Ohana/src/relphot/src/ZeroPointModes.c

    r41467 r41473  
    11# include "relphot.h"
    22
    3 void SetZeroPointModes (int nloop) {
     3static int CurrentLoop = -1; // track the number of iterations
     4
     5void SetZptIteration (int current) {
     6  CurrentLoop = current;
     7}
     8
     9int GetZptIteration (void) {
     10  return CurrentLoop;
     11}
     12
     13void SetZeroPointModes (void) {
    414
    515  if (TGROUP_ZEROPT && MOSAIC_ZEROPT) {
    6     if ((nloop > 4) && (nloop <= 8)) {
     16    if ((CurrentLoop > 4) && (CurrentLoop <= 8)) {
    717      // after iterating a few times on the TGroup zero points:
    818      // * identify the photometric nights
    … …  
    1323      clean_tgroups(); // do this on every pass to update the status of nights
    1424    }
    15     if ((nloop > 8) && (nloop <= 999)) {
     25    if ((CurrentLoop > 8) && (CurrentLoop <= 999)) {
    1626      // after iterating a few times on the TGroup & Mosaic zero points:
    1727      // * identify the good / bad Mosaics
    … …  
    2636
    2737  if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
    28     if ((nloop > 4) && (nloop <= 999)) {
     38    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
    2939      // after iterating a few times on the TGroup zero points:
    3040      // * identify the photometric nights
    … …  
    3848
    3949  if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) {
    40     if ((nloop > 4) && (nloop <= 999)) {
     50    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
    4151      // after iterating a few times on the Mosaic zero points:
    4252      // * identify the good / bad mosaics
    … …  
    5060
    5161  if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
    52     if ((nloop > 4) && (nloop <= 999)) {
     62    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
    5363      // after iterating a few times on the Image zero points:
    5464      // * ???
    … …  
    6171
    6272// mode = TGROUP, MOSAIC, IMAGE
    63 // should we use IRLS or OLS?  depends on nloop and which stage we are in.
    64 int UseStandardOLS (int nloop, ZptFitModeType mode) {
     73// should we use IRLS or OLS?  depends on CurrentLoop and which stage we are in.
     74int UseStandardOLS (ZptFitModeType mode) {
    6575
    6676  if (TGROUP_ZEROPT && MOSAIC_ZEROPT) {
    6777    switch (mode) {
     78      case ZPT_STARS:
     79        if (CurrentLoop < 3) return TRUE;
     80        return FALSE;
    6881      case ZPT_TGROUP:
    69         if (nloop < 3) return TRUE;
     82        if (CurrentLoop < 3) return TRUE;
    7083        return FALSE;
    7184      case ZPT_MOSAIC:
    72         if (nloop < 7) return TRUE;
     85        if (CurrentLoop < 7) return TRUE;
    7386        return FALSE;
    7487      case ZPT_IMAGES:
    7588      default:
    76         if (nloop < 11) return TRUE;
     89        if (CurrentLoop < 11) return TRUE;
    7790        return FALSE;
    7891    }
    … …  
    8194  if (TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
    8295    switch (mode) {
     96      case ZPT_STARS:
     97        if (CurrentLoop < 3) return TRUE;
     98        return FALSE;
    8399      case ZPT_TGROUP:
    84100      case ZPT_MOSAIC:
    85         if (nloop < 3) return TRUE;
     101        if (CurrentLoop < 3) return TRUE;
    86102        return FALSE;
    87103      case ZPT_IMAGES:
    88104      default:
    89         if (nloop < 7) return TRUE;
     105        if (CurrentLoop < 7) return TRUE;
    90106        return FALSE;
    91107    }
    … …  
    94110  if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) {
    95111    switch (mode) {
     112      case ZPT_STARS:
     113        if (CurrentLoop < 3) return TRUE;
     114        return FALSE;
    96115      case ZPT_TGROUP:
    97116      case ZPT_MOSAIC:
    98         if (nloop < 3) return TRUE;
     117        if (CurrentLoop < 3) return TRUE;
    99118        return FALSE;
    100119      case ZPT_IMAGES:
    101120      default:
    102         if (nloop < 7) return TRUE;
     121        if (CurrentLoop < 7) return TRUE;
    103122        return FALSE;
    104123    }
    … …  
    107126  if (!TGROUP_ZEROPT && !MOSAIC_ZEROPT) {
    108127    switch (mode) {
     128      case ZPT_STARS:
     129        if (CurrentLoop < 3) return TRUE;
     130        return FALSE;
    109131      case ZPT_TGROUP:
    110132      case ZPT_MOSAIC:
    111133      case ZPT_IMAGES:
    112134      default:
    113         if (nloop < 3) return TRUE;
     135        if (CurrentLoop < 3) return TRUE;
    114136        return FALSE;
    115137    }

  • trunk/Ohana/src/relphot/src/args.c

    r41462 r41473  
    446446    }
    447447    remove_argument (N, &argc, argv);
     448  }
     449  USE_REFERENCE_WEIGHT = FALSE;
     450  if ((N = get_argument (argc, argv, "-ref-weight"))) {
     451    remove_argument (N, &argc, argv);
     452    USE_REFERENCE_WEIGHT = TRUE;
    448453  }
    449454

  • trunk/Ohana/src/relphot/src/relphot_images.c

    r41467 r41473  
    116116    /* determine fit values */
    117117    for (i = 0; i < NLOOP; i++) {
     118      SetZptIteration (i);
    118119      dump_tgroups (catalog, i);
    119120      dump_catalog (catalog, 0, i); // for a test, just dump a specific catalog
    … …  
    123124      dump_catalog (catalog, 0, i + 20); // for a test, just dump a specific catalog
    124125
    125       setMcal  (catalog, i, FALSE, flatcorr);
    126       setMmos  (catalog, i, FALSE, flatcorr);
    127       setMgrp  (catalog, i, FALSE, flatcorr);
     126      setMcal  (catalog, FALSE, flatcorr);
     127      setMmos  (catalog, FALSE, flatcorr);
     128      setMgrp  (catalog, FALSE, flatcorr);
    128129      dump_tgroups (catalog, i + 40);
    129130      dump_catalog (catalog, 0, i + 40); // for a test, just dump a specific catalog
    … …  
    143144      global_stats (catalog, Ncatalog, flatcorr, i);
    144145
    145       SetZeroPointModes (i);
     146      SetZeroPointModes ();
    146147      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    147148      fprintf (stderr, "----- continue loop %d -----\n", i);
    … …  
    162163    int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
    163164    onlyPoorImages = FALSE; // XXX deactivate for now
    164     setMcal  (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass
    165     setMmos  (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass
    166     setMgrp  (catalog, 999, onlyPoorImages, flatcorr); // max loop since this is a final pass
     165    setMcal  (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass
     166    setMmos  (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass
     167    setMgrp  (catalog, onlyPoorImages, flatcorr); // max loop since this is a final pass
    167168    MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    168169

  • trunk/Ohana/src/relphot/src/relphot_parallel_images.c

    r41462 r41473  
    9595  }
    9696
    97   // XXX : add this in the loop at various points
    98   // if (PLOTSTUFF) plot_scatter (catalog, Ncatalog, flatcorr);
    99 
    10097  /* determine fit values */
    10198  client_logger_message ("starting the loops : %d \n", NLOOP);
    10299  for (i = 0; i < NLOOP; i++) {
     100    SetZptIteration (i);
    103101
    104102    // set the mean stellar mags given the measurements and the image calibrations
    … …  
    114112
    115113    // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements
    116     setMcal  (catalog, i, FALSE, flatcorr);
    117     setMmos  (catalog, i, FALSE, flatcorr);
     114    setMcal  (catalog, FALSE, flatcorr);
     115    setMmos  (catalog, FALSE, flatcorr);
    118116    // setMgrp  (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups
    119117    MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    … …  
    126124    client_logger_message ("slurped image data : loop %d \n", i);
    127125
    128     if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
    129     if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
    130     if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
    131     if ((i > 8) && (i % 8 == 5)) clean_images ();
     126    global_stats (catalog, Ncatalog, flatcorr, i);
    132127
    133     global_stats (catalog, Ncatalog, flatcorr, i);
     128    SetZeroPointModes ();
    134129    MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    135130  }
    … …  
    151146  /* set Mcal & Mmos for bad images (for stack_and_warps, force images to be measured) */
    152147  int onlyPoorImages = !CALIBRATE_STACKS_AND_WARPS;
    153   setMcal  (catalog, 999, onlyPoorImages, flatcorr);
    154   setMmos  (catalog, 999, onlyPoorImages, flatcorr);
     148  setMcal  (catalog, onlyPoorImages, flatcorr);
     149  setMmos  (catalog, onlyPoorImages, flatcorr);
    155150  // setMgrp  (catalog, onlyPoorImages, flatcorr); XXX see above
    156151  MARKTIME("-- finalize Mcal values: %f sec\n", dtime);

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r41462 r41473  
    1919
    2020// # define MAG_STATS_BY_RANKING magStatsByRankingIRLS
    21 # define MAG_STATS_BY_RANKING magStatsByRanking
     21// # define MAG_STATS_BY_RANKING magStatsByRanking
    2222
    2323# define UBERCAL_WEIGHT 100.0
    … …  
    6060FILE *TESTFILE = NULL;
    6161
     62int useOLS = TRUE; // chosen based on iteration below
     63
    6264int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    6365
    … …  
    7375  TESTFILE = fopen ("mytest.dat", "w");
    7476 
     77  useOLS = UseStandardOLS (ZPT_STARS);
    7578  for (j = 0; j < catalog[Nc].Naverage; j++) {
    7679
    … …  
    258261    if (isUbercal) results->haveUbercal[Nsec] = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure XXX define this array
    259262
    260     int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code));
     263    // logic for choosing a modified weight:
     264    float modifiedWeight = 1.0;
     265    if (isUbercal) modifiedWeight = UBERCAL_WEIGHT;
     266    if (refPhotcode && (code->code == refPhotcode->code)) modifiedWeight = UBERCAL_WEIGHT;
     267    if (USE_REFERENCE_WEIGHT && (code->type == PHOT_REF)) modifiedWeight = code->photomErrSys; // we are overloading this field for now
    261268
    262269    // combine the various errors in quadrature here:
    … …  
    277284        results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
    278285        results->aperData[Nsec].errlist[Nap] = dMap;
    279         results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
     286        results->aperData[Nsec].wgtlist[Nap] = modifiedWeight;
    280287        results->aperData[Nsec].ranking[Nap] = measureRank[k];
    281288        results->aperData[Nsec].measSeq[Nap] = k;
    … …  
    297304      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
    298305      results->kronData[Nsec].errlist[Nkron] = dMkron;
    299       results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     306      results->kronData[Nsec].wgtlist[Nkron] = modifiedWeight;
    300307      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
    301308      results->kronData[Nsec].measSeq[Nkron] = k;
    … …  
    310317      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
    311318      results->psfData[Nsec].errlist[Npsf] = dMpsf;
    312       results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
     319      results->psfData[Nsec].wgtlist[Npsf] = modifiedWeight;
    313320      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
    314321      results->psfData[Nsec].measSeq[Npsf] = k;
    … …  
    406413    // if too few valid measurements meet the minimum criteria, go to the next entry
    407414    StatType *psfstats = &results->psfstats;
    408     int NrankingPSF = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
     415
     416    int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
    409417    if (NrankingPSF < Nminmeas) {
    410418      secfilt[Nsec].flags |= ID_OBJ_FEW;
    … …  
    422430
    423431    StatType *kronstats = &results->kronstats;
    424     int NrankingKRON = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
     432    int NrankingKRON = useOLS ? magStatsByRanking (&results->kronData[Nsec], kronstats) : magStatsByRankingIRLS (&results->kronData[Nsec], kronstats);
    425433    if (NrankingKRON) {
    426434      secfilt[Nsec].MkronChp  = kronstats->mean;
    … …  
    455463
    456464      StatType *apstats = &results->apstats;
    457       int NrankingAPER = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
     465      int NrankingAPER = useOLS ? magStatsByRanking (&results->aperData[Nsec], apstats) : magStatsByRankingIRLS (&results->aperData[Nsec], apstats);
    458466      if (NrankingAPER) {
    459467        secfilt[Nsec].MapChp  = apstats->mean;
    … …  
    10321040    // if too few valid measurements meet the minimum criteria, go to the next entry
    10331041    StatType *psfstats = &results->psfstats;
    1034     Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
     1042    // Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
     1043    Nranking = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
    10351044    if (Nranking) {
    10361045      secfilt[Nsec].FpsfWrp  = psfstats->mean;
    … …  
    10451054    // if too few valid measurements meet the minimum criteria, go to the next entry
    10461055    StatType *apstats = &results->apstats;
    1047     Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
     1056    // Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
     1057    Nranking = useOLS ? magStatsByRanking (&results->aperData[Nsec], apstats) : magStatsByRankingIRLS (&results->aperData[Nsec], apstats);
    10481058    if (Nranking) {
    10491059      secfilt[Nsec].FapWrp     = apstats->mean;
    … …  
    10581068    // if too few valid measurements meet the minimum criteria, go to the next entry
    10591069    StatType *kronstats = &results->kronstats;
    1060     Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
     1070    // Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
     1071    Nranking = useOLS ? magStatsByRanking (&results->kronData[Nsec], kronstats) : magStatsByRankingIRLS (&results->kronData[Nsec], kronstats);
    10611072    if (Nranking) {
    10621073      secfilt[Nsec].FkronWrp     = kronstats->mean;
Note: See TracChangeset for help on using the changeset viewer.