IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 1, 2012, 3:00:19 PM (14 years ago)
Author:
eugene
Message:

updates to relastro to support distributed dvo; some optmimization (eg, image matching to area); cleanup names of static helper arrays in ImageOps, MosaicOps; cleanup rules for -reset options (relastro does not reset photom bits); skip objects with NAN error bars (can be used to exclude some photcodes from astrometry); enum for mask bits in FitChip and related; choose max order using same rules as psastro; plug leaks; reduce verbosity; big re-org of main (for parallel code); only use accepted measurements for Tmin,Tmax,Trange,Tmean; require a min parallax factor to fit parallaxes; use MeasureTiny where possible; use new ParsePhotcodeList function; testparallax function; threaded version of UpdateChips; bcatalog limit density selecting objects by number of measurements; high-speed now copies support tables to output database; high-speed may write to a parallel database; test for unsorted database (& fail if true)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/src/high_speed_objects.c

    r31635 r33652  
    66int high_speed_objects (SkyRegion *region, Catalog *catalog) {
    77
    8   // XXX seems silly to require region here just to find the catalog bounds / center
    9 
    10   off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad,l,i1,Noff,Nmeasure,Naverage, NAVERAGE, NMEASURE;
    11   int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj;
    12   int foundA, foundB;
     8  off_t i, j, m, J, ni, nj, *N1;
     9  off_t Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad;
     10  off_t l, i1, NAVERAGE, NMEASURE;
     11  int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas, Nepoch, nv[2], Nmatchmeasobj;
     12  int foundA, foundB, XVERB;
    1313  double *X1, *Y1;
    1414  double dX, dY, dR, RADIUS2;
    1515  Coords tcoords;
    16   Catalog catalog1;
     16  Catalog catalogOut;
    1717
    1818  int Nsecfilt;
    1919  char filename[1024];
    20   Noff = strlen(CATDIR);
    21   sprintf (filename, "%s%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]);
    22   fprintf(stderr, "%s\n",filename);
    23   dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/
    24   catalog1.filename = strcreate(filename);
     20
     21  snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name);
     22  fprintf (stderr, "%s\n",filename);
     23
     24  dvo_catalog_init (&catalogOut, TRUE); /*initialise new catalogue*/
     25  catalogOut.filename = strcreate(filename);
    2526
    2627  SIGMA_LIM = 0.0;
    2728
    2829  Nsecfilt = GetPhotcodeNsecfilt();
    29   Naverage=catalog[0].Naverage;
    30   Nmeasure=catalog[0].Nmeasure;
    31   catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion
    32   catalog1.Nsecfilt  = Nsecfilt;
    33   // always load all of the data (if any exists)
    34   catalog1.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     30  // off_t Naverage = catalog[0].Naverage;
     31  // off_t Nmeasure = catalog[0].Nmeasure;
     32  catalogOut.filename = filename; // based on the input name, need to keep everything below the catdir portion
     33  catalogOut.Nsecfilt = Nsecfilt;
     34  catalogOut.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; // load all data
    3535 
    36   catalog1.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    37   catalog1.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    38 
    39   if (!dvo_catalog_open (&catalog1, region, VERBOSE,"w")) {
     36  catalogOut.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
     37  catalogOut.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
     38
     39  if (!dvo_catalog_open (&catalogOut, region, VERBOSE2,"w")) {
    4040    fprintf (stderr, "ERROR: failure to open catalog file %s\n",
    4141             filename);
    … …  
    4545  NAVERAGE = 1000;
    4646  NMEASURE = 10000;
    47   REALLOCATE (catalog1.average, Average, NAVERAGE);
    48   REALLOCATE (catalog1.measure, Measure, NMEASURE);
    49   REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
     47  REALLOCATE (catalogOut.average, Average, NAVERAGE);
     48  REALLOCATE (catalogOut.measure, Measure, NMEASURE);
     49  REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
    5050
    5151  // high-speed between different surveys (easier case):
    … …  
    7777  // mark (exclude) objects with both sets of target photcodes
    7878  for (i = 0; i < catalog[0].Naverage; i++) {
     79
     80    XVERB  = (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC);
     81    XVERB |= (catalog[0].average[i].objID == OBJ_ID_DST) && (catalog[0].average[i].catID == CAT_ID_DST);
     82    if (XVERB) {
     83      fprintf (stderr, "test object\n");
     84    }
     85
    7986    // do any of the measures for this object match group A?
    8087    m = catalog[0].average[i].measureOffset;
    … …  
    114121            // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n",  i);
    115122        }
    116 
    117     } else if (foundB && !foundA) {
    118 
     123        continue;
     124    }
     125    if (foundB && !foundA) {
    119126        if (applyConstraintsB(catalog, i)) {
    120127            groupB[i] = TRUE;
    … …  
    124131            NgroupBbad++;
    125132        }
     133        continue;
    126134    }
    127135    // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow
    … …  
    173181    nj = N1[j];
    174182
     183    XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
     184    XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
     185    if (XVERB) {
     186      fprintf (stderr, "test object %d or %d\n", (int) ni, (int) nj);
     187    }
     188
    175189    if (slowMoving[ni]) NEXT_I;
    176190    if (slowMoving[nj]) NEXT_J;
    … …  
    185199    // look for pairs that are within the maximum separation
    186200    dX = X1[i] - X1[j];
     201    if (XVERB) {
     202      fprintf (stderr, "%d %d : %f\n", (int) ni, (int) nj, dX);
     203    }
    187204
    188205    if (dX <= -1.02*RADIUS) NEXT_I; // negative dX: i is too small
    … …  
    197214      nj = N1[J];
    198215      if (!groupB[nj]) continue;
     216
     217      XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
     218      XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
     219      if (XVERB) {
     220        fprintf (stderr, "test object pt2 %d or %d : %f\n", (int) ni, (int) nj, dX);
     221        fprintf (stderr, ".\n");
     222      }
    199223
    200224      dY = Y1[i] - Y1[J];
    … …  
    209233      nv[1]=nj;
    210234
    211       catalog1.average[Nmatch]=catalog[0].average[nv[0]];
     235      catalogOut.average[Nmatch]=catalog[0].average[nv[0]];
    212236      /*Loop over index values and set measurements in new catalogue*/
    213237      Nmatchmeasobj=0;
    214       catalog1.average[Nmatch].measureOffset=Nmatchmeas;
     238      catalogOut.average[Nmatch].measureOffset=Nmatchmeas;
    215239      for(l=0;l<Nepoch;l++) {
    216240          m = catalog[0].average[nv[l]].measureOffset;
    217241          for(i1=0;i1<catalog[0].average[nv[l]].Nmeasure;i1++) {
    218               catalog1.measure[Nmatchmeas]=catalog[0].measure[m+i1];
     242              catalogOut.measure[Nmatchmeas]=catalog[0].measure[m+i1];
    219243              /*Set offset RA and Dec wrt correct average value*/
    220               catalog1.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
    221               catalog1.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
    222               catalog1.measure[Nmatchmeas].averef = Nmatch;
     244              catalogOut.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
     245              catalogOut.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
     246              catalogOut.measure[Nmatchmeas].averef = Nmatch;
    223247              Nmatchmeasobj++;
    224248              Nmatchmeas++;
    … …  
    226250              if (Nmatchmeas == NMEASURE - 1) {
    227251                NMEASURE += 10000;
    228                 REALLOCATE (catalog1.measure, Measure, NMEASURE);
     252                REALLOCATE (catalogOut.measure, Measure, NMEASURE);
    229253              }
    230254          }
    231255      }
    232       catalog1.average[Nmatch].Nmeasure=Nmatchmeasobj;
     256      catalogOut.average[Nmatch].Nmeasure = Nmatchmeasobj;
    233257      Nmatch ++;
    234258
    235259      if (Nmatch == NAVERAGE - 1) {
    236260        NAVERAGE += 1000;
    237         REALLOCATE (catalog1.average, Average, NAVERAGE);
    238         REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
     261        REALLOCATE (catalogOut.average, Average, NAVERAGE);
     262        REALLOCATE(catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
    239263      }
    240264    }
    241265    i++;
    242266  }
    243   catalog1.Naverage=Nmatch;
    244   catalog1.Nmeasure=Nmatchmeas;
    245   catalog1.Nsecfilt=Nsecfilt;
    246   catalog1.Nsecf_mem=Nmatch*Nsecfilt;
    247 
    248   UpdateObjects(&catalog1, 1);
     267  catalogOut.Naverage=Nmatch;
     268  catalogOut.Nmeasure=Nmatchmeas;
     269  catalogOut.Nsecfilt=Nsecfilt;
     270  catalogOut.Nsecf_mem=Nmatch*Nsecfilt;
     271
     272  populate_tiny_values (&catalogOut, DVO_TV_MEASURE);
     273  UpdateObjects (&catalogOut, 1);
    249274
    250275  fprintf (stderr, "found %d matches\n", Nmatch);
    251   dvo_catalog_save (&catalog1, VERBOSE);
    252   dvo_catalog_unlock (&catalog1);
    253   dvo_catalog_free (&catalog1);
     276  dvo_catalog_save (&catalogOut, VERBOSE2);
     277  dvo_catalog_unlock (&catalogOut);
     278  dvo_catalog_free (&catalogOut);
    254279  free (slowMoving);
    255280  free (groupA);
    … …  
    262287}
    263288
    264 // return TRUE if any measure->photcodes match the photcodeSet
     289// return TRUE if measure->photcode match any in the photcodeSet
     290// (but, return FALSE if the measurement is bad)
    265291int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset) {
    266292
    267293  int found, k;
    268294
    269   assert (Nset > 0);
    270 
    271   if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
    272   if (!finite(measure[0].M)) return FALSE;
     295  if (!Nset) return TRUE;
     296  // XXX require a set or not?  assert (Nset > 0);
     297
     298  if (!finite(measure[0].dR)) return FALSE;
     299  if (!finite(measure[0].dD)) return FALSE;
     300  if (!finite(measure[0].M))  return FALSE;
    273301 
     302  float dX = GetAstromError (measure, ERROR_MODE_RA);
     303  if (isnan(dX)) return FALSE;
     304
     305  float dY = GetAstromError (measure, ERROR_MODE_DEC);
     306  if (isnan(dY)) return FALSE;
     307
    274308  /* select measurements by photcode, or equiv photcode, if specified */
    275   if (Nset > 0) {
    276     found = FALSE;
    277     for (k = 0; (k < Nset) && !found; k++) {
    278       if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
    279       if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    280     }
    281     if (!found) return FALSE;
    282   }
    283  
    284   return TRUE;
     309  found = FALSE;
     310  for (k = 0; (k < Nset) && !found; k++) {
     311    if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
     312    if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
     313  }
     314
     315  return found;
    285316}
Note: See TracChangeset for help on using the changeset viewer.