IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39525 for trunk


Ignore:
Timestamp:
Apr 9, 2016, 12:01:06 PM (10 years ago)
Author:
eugene
Message:

finished setobjflags

Location:
trunk/Ohana/src/addstar
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/include/setobjflags.h

    r39496 r39525  
    22typedef struct {
    33  double R, D;
    4   e_time time;
     4//e_time time;
     5//short filter;
    56  int myBit;
    67  int flag; // in a subset?
    … …  
    2324int args_setobjflags_client (int *argc, char **argv);
    2425
    25 int setobjflags_table (SkyList *skylistInput, HostTable *hosts, char *filename, int myBit);
     26int setobjflags_table (SkyList *skylistInput, HostTable *hosts);
    2627
    2728MyStars *setobjflags_make_subset (MyStars *stars, int Nstars, int start, SkyRegion *region, int *nsubset);
    … …  
    4647MyStars *setobjflags_load_stars (char *filename, int *nstars);
    4748
    48 MyStars *setobjflags_loadfile (char *file, int mybit, int *Nstars);
    49 // MyStars *setobjflags_readstars (char *filename, int *nstars);
     49MyStars *setobjflags_loadfile (int *Nstars);
    5050
    5151int setobjflags_sortStars (MyStars *stars, int Nstars);

  • trunk/Ohana/src/addstar/src/args_setobjflags.c

    r39496 r39525  
    1515  }
    1616
    17   /* we do not allow a subset to be extracted -- all or nothing, babe */
     17  /* specify portion of the sky (useful for testing) */
    1818  UserPatch.Rmin = 0;
    19   UserPatch.Rmax= 360;
     19  UserPatch.Rmax = 360;
    2020  UserPatch.Dmin = -90;
    2121  UserPatch.Dmax = +90;
     22  if ((N = get_argument (*argc, argv, "-region"))) {
     23    remove_argument (N, argc, argv);
     24    UserPatch.Rmin = atof (argv[N]);
     25    remove_argument (N, argc, argv);
     26    UserPatch.Rmax = atof (argv[N]);
     27    remove_argument (N, argc, argv);
     28    UserPatch.Dmin = atof (argv[N]);
     29    remove_argument (N, argc, argv);
     30    UserPatch.Dmax = atof (argv[N]);
     31    remove_argument (N, argc, argv);
     32  }
    2233
    2334  // XXX for the moment, make this selection manual.  it needs to be automatic
    … …  
    5768  }
    5869
    59   /* extra error messages */
    60   SRC_FILES = NULL;
    61   SRC_MYBIT = 0;
    62   if ((N = get_argument (*argc, argv, "-src")) == 0) help ("missing -src argument: required");
    63   if (N > *argc - 3) help("-src is missing (file) or (bits) values");
    64   remove_argument (N, argc, argv);
    65   SRC_FILES = strcreate (argv[N]);
    66   remove_argument (N, argc, argv);
    67   SRC_MYBIT = atoi (argv[N]);
    68   remove_argument (N, argc, argv);
    69 
    70   if (*argc < 3) help ("insufficient arguments");
     70  if (*argc != 1) help ("insufficient arguments");
    7171  return TRUE;
    7272}
    … …  
    143143  if (message) fprintf (stderr, "ERROR: %s\n\n", message);
    144144
    145   fprintf (stderr, "USAGE: setobjflags [options] -src (file) (flag) [..more -src (file) (flag) sets]\n");
     145  fprintf (stderr, "USAGE: setobjflags [options]\n");
    146146  fprintf (stderr, "  add flags for specified detection sets\n\n");
    147147

  • trunk/Ohana/src/addstar/src/find_matches_setobjflags.c

    r39496 r39525  
    2222  ALLOCATE (Y2, double, NAVE);
    2323  ALLOCATE (N2, off_t,  NAVE);
     24
     25  // use this structure to count the number of mops detections
     26  ALLOCATE_ZERO (catalog[0].nOwn_t, int,  NAVE);
    2427
    2528  /* internal counters */
    … …  
    8588    }
    8689
    87     /* this block will match a given detection to the closest object within range of that detection.
    88        XXX note that this matches ALL detections within range of the single object to that same object
    89        this is bad, but I cannot just go in linear order (ie, mark off each object as they are
    90        used).  I should make a list of all Nobj * Ndet pairs in range and choose the matches
    91        based on their separations.  UGH
     90    /* this block will match a given detection to the closest object within range of that
     91       detection.  Note that this matches ALL detections within range of the single object
     92       to that same object.  In this case (setobjflags), this is actually the behavior I
     93       want: I have multiple independent star lists which need to be assigned to the
     94       database.
    9295     */
    9396   
    … …  
    120123    /** set bit for this star **/
    121124    catalog[0].average[n].flags |= stars[N].myBit;
     125    if (stars[N].myBit & ID_OBJ_HAS_SOLSYS_DET) {
     126      catalog[0].nOwn_t[n] ++;
     127      if (catalog[0].nOwn_t[n] >= 0.5*catalog[0].average[n].Nmeasure) {
     128        catalog[0].average[n].flags |= ID_OBJ_MOST_SOLSYS_DET;
     129      }
     130    }
    122131    stars[N].found = TRUE;
    123132    i++;

  • trunk/Ohana/src/addstar/src/setobjflags.c

    r39496 r39525  
    4848  SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch);
    4949
    50   setobjflags_table (skylist, hosts, SRC_FILES, SRC_MYBIT);
     50  setobjflags_table (skylist, hosts);
    5151
    5252  FreeConfig ();

  • trunk/Ohana/src/addstar/src/setobjflags_io.c

    r39496 r39525  
    2929  gfits_define_bintable_column (&theader, "D", "DEC",      "", "degree", 1.0, 0.0);
    3030  gfits_define_bintable_column (&theader, "J", "MYBIT",    "", "", 1.0, FT_BZERO_INT32);
    31 //gfits_define_bintable_column (&theader, "J", "averef",   "", "", 1.0, FT_BZERO_INT32);
    32 //gfits_define_bintable_column (&theader, "J", "objID",    "", "", 1.0, FT_BZERO_INT32);
    33 //gfits_define_bintable_column (&theader, "J", "catID",    "", "", 1.0, FT_BZERO_INT32);
    3431
    3532  int i;
    … …  
    3734  double       *dec      = NULL; ALLOCATE (dec     , double      , Nstars); for (i = 0; i < Nstars; i++) dec     [i] = stars[i].D      ;
    3835  unsigned int *myBit    = NULL; ALLOCATE (myBit   , unsigned int, Nstars); for (i = 0; i < Nstars; i++) myBit   [i] = stars[i].myBit  ;
    39 //unsigned int *averef   = NULL; ALLOCATE (averef  , unsigned int, Nstars); for (i = 0; i < Nstars; i++) averef  [i] = stars[i].starpar.averef  ;
    40 //unsigned int *objID    = NULL; ALLOCATE (objID   , unsigned int, Nstars); for (i = 0; i < Nstars; i++) objID   [i] = stars[i].starpar.objID   ;
    41 //unsigned int *catID    = NULL; ALLOCATE (catID   , unsigned int, Nstars); for (i = 0; i < Nstars; i++) catID   [i] = stars[i].starpar.catID   ;
    4236
    4337  // generate the output array that carries the data
    … …  
    113107  GET_COLUMN(dec     , "DEC"       ,   double      ,   double      );
    114108  GET_COLUMN(myBit   , "MYBIT"     ,   unsigned int,   int         );
    115 //GET_COLUMN(averef  , "averef"    ,   unsigned int,   int         );
    116 //GET_COLUMN(objID   , "objID"     ,   unsigned int,   int         );
    117 //GET_COLUMN(catID   , "catID"     ,   unsigned int,   int         );
    118109
    119110  gfits_free_header (&theader);
    … …  
    128119    stars[i].D                = dec    [i];
    129120    stars[i].myBit            = myBit  [i];
     121    stars[i].flag             = FALSE;
     122    stars[i].found            = FALSE;
    130123  }
    131124

  • trunk/Ohana/src/addstar/src/setobjflags_loadfile.c

    r39496 r39525  
    11# include "addstar.h"
    22# include "setobjflags.h"
     3void sort_exp (int *I, int *T, int *F, int *S, int N);
     4int expname_to_int (char *expname);
     5int value_to_sequence (int *sequence, int Nsequence, int value);
     6
     7MyStars *setobjflags_loadfile_mops (char *metadata, char *detections, int mybit, int *Nstars);
     8MyStars *setobjflags_loadfile_hern (char *filename, int *Nstars);
    39
    410/* read ASCII file with ref star data */
    5 MyStars *setobjflags_loadfile (char *file, int mybit, int *Nstars) {
    6 
    7   FILE *f;
    8   char line[256];
    9 
    10   /* open file */
    11   f = fopen (file, "r");
    12   if (f == NULL) Shutdown ("can't read data from %s", file);
     11MyStars *setobjflags_loadfile (int *Nstars) {
     12
     13  int Nmops;
     14  MyStars *mopsStars = setobjflags_loadfile_mops ("md.fix.dat", "mops_export.txt", ID_OBJ_HAS_SOLSYS_DET, &Nmops);
     15
     16  int Nhern;
     17  MyStars *hernStars = setobjflags_loadfile_hern ("NH_qso_rrlyrae_candidates_catalog.fits", &Nhern);
     18
     19  // merge into a single table
     20  int Ntotal = Nhern + Nmops;
     21  REALLOCATE (mopsStars, MyStars, Ntotal);
     22  for (int i = Nmops, j = 0; i < Ntotal; i++, j++) {
     23    mopsStars[i] = hernStars[j];
     24  }
     25  free (hernStars);
     26
     27  *Nstars = Ntotal;
     28  return mopsStars;
     29}
     30
     31MyStars *setobjflags_loadfile_mops (char *metadata, char *detections, int mybit, int *Nstars) {
     32
     33  // code to read in the MOPS detection dump
     34  int Nexp = 0;
     35  int NEXP = 1000;
     36  ALLOCATE_PTR (expTIM,    int, NEXP);
     37  ALLOCATE_PTR (expSEQ,    int, NEXP);
     38  ALLOCATE_PTR (expFLT,    int, NEXP);
     39  ALLOCATE_PTR (expIDX,    int, NEXP);
     40
     41  // load the metadata file:
     42  FILE *f = fopen (metadata, "r");
     43  if (f == NULL) Shutdown ("can't read data from %s", metadata);
     44   
     45  /* read in stars line-by-line */
     46  char line[1024];
     47  while (scan_line (f, line) != EOF) {
     48    if (Nexp % 1000 == 0) fprintf (stderr, ".");
     49
     50    char filter[64], expname[64];
     51    double expMJD;
     52
     53    int status = sscanf(line, "%lf %s %s", &expMJD, filter, expname);
     54    if (status != 3) {
     55      fprintf (stderr, "error reading line from MD: %s\n", line);
     56      continue;
     57    }
     58
     59    // XXX NOTE this does not handle gpc2 images
     60    expIDX[Nexp] = expname_to_int (expname);
     61
     62    switch (filter[0]) {
     63      case 'g': expFLT[Nexp] = 0; break;
     64      case 'r': expFLT[Nexp] = 1; break;
     65      case 'i': expFLT[Nexp] = 2; break;
     66      case 'z': expFLT[Nexp] = 3; break;
     67      case 'y': expFLT[Nexp] = 4; break;
     68      case 'w': expFLT[Nexp] = 5; break;
     69      default: myAbort ("oops");
     70    }
     71
     72    expSEQ[Nexp] = Nexp;
     73    expTIM[Nexp] = ohana_mjd_to_sec (expMJD);
     74
     75    Nexp++ ;
     76    if (Nexp == NEXP) {
     77      NEXP += 1000;
     78      REALLOCATE (expTIM,    int, NEXP);
     79      REALLOCATE (expSEQ,    int, NEXP);
     80      REALLOCATE (expFLT,    int, NEXP);
     81      REALLOCATE (expIDX,    int, NEXP);
     82    }
     83  }
     84  fclose (f);
     85  sort_exp (expIDX, expTIM, expFLT, expSEQ, Nexp);
     86
     87  fprintf (stderr, "done reading exposure metadata\n");
    1388
    1489  MyStars *stars = NULL;
    … …  
    1691  ALLOCATE (stars, MyStars, NSTARS);
    1792
     93  // load the detection file:
     94  f = fopen (detections, "r");
     95  if (f == NULL) Shutdown ("can't read data from %s", detections);
     96   
     97  scan_line (f, line); // skip the first line
     98
    1899  /* read in stars line-by-line */
    19   int N = 0;
     100  int N = 0, Nbad = 0, Nmiss = 0;
    20101  while (scan_line (f, line) != EOF) {
     102    if (N % 1000 == 0) fprintf (stderr, ".");
    21103    stripwhite (line);
    22104    if (line[0] == 0) continue;
    23105    if (line[0] == '#') continue;
    24 
    25     int status = sscanf (line, "%lf %lf", &stars[N].R, &stars[N].D);
    26     if (status != 2) {
    27       fprintf (stderr, "error reading line: %s\n", line);
     106     
     107    char expname[64];
     108    int status = sscanf (line, "%s %lf %lf", expname, &stars[N].R, &stars[N].D);
     109    if (status != 3) {
     110      if (Nbad < 20) fprintf (stderr, "error reading line: %s\n", line);
     111      Nbad ++;
    28112      continue;
    29113    }
    30114    stars[N].R = ohana_normalize_angle (stars[N].R);
    31    
     115     
     116    int expID = expname_to_int (expname);
     117    int seq = value_to_sequence (expIDX, Nexp, expID);
     118    if (seq < 0) {
     119      if (Nmiss < 20) fprintf (stderr, "error missing expname: %s\n", expname);
     120      Nmiss ++;
     121    }
     122
    32123    stars[N].myBit = mybit;
    33 
    34     stars[N].time = 0;
    35     stars[N].flag = FALSE;
    36     stars[N].found = FALSE;
    37 
     124     
     125    // stars[N].time   = expTIM[seq];
     126    // stars[N].filter = expFLT[seq];
     127    stars[N].flag   = FALSE;
     128    stars[N].found  = FALSE;
     129     
    38130    N++;
    39131    CHECK_REALLOCATE (stars, MyStars, NSTARS, N, 1000);
    40132  }
     133  fclose (f);
    41134  *Nstars = N;
     135
     136  fprintf (stderr, "done reading mops detections\n");
     137
     138  FREE (expTIM);
     139  FREE (expSEQ);
     140  FREE (expFLT);
     141  FREE (expIDX);
     142
    42143  return (stars);
    43144}
     145
     146# define GET_COLUMN(OUT,NAME,TYPE)                                      \
     147  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
     148  myAssert (!strcmp(type, #TYPE), "wrong column type");
     149
     150# define MIN_QSO_P60_VALUE 0.60
     151# define MIN_QSO_P05_VALUE 0.05
     152
     153# define MIN_RRLYRA_P60_VALUE 0.60
     154# define MIN_RRLYRA_P05_VALUE 0.05
     155
     156MyStars *setobjflags_loadfile_hern (char *filename, int *Nstars) {
     157
     158  int i, Ncol;
     159  off_t Nrow;
     160  char type[16];
     161
     162  Header header;
     163  Header theader;
     164  Matrix matrix;
     165  FTable ftable;
     166
     167  header.buffer = NULL;
     168  matrix.buffer = NULL;
     169  ftable.buffer = NULL;
     170  theader.buffer = NULL;
     171 
     172  FILE *f = fopen (filename, "r");
     173  if (!f) {
     174    fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename);
     175    return NULL;
     176  }
     177
     178  /* load in PHU segment (ignore) */
     179  if (!gfits_fread_header (f, &header)) {
     180    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
     181    goto escape;
     182  }
     183  if (!gfits_fread_matrix (f, &matrix, &header)) {
     184    if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n");
     185    goto escape;
     186  }
     187
     188  ftable.header = &theader;
     189
     190  // load data for this header
     191  if (!gfits_load_header (f, &theader)) goto escape;
     192 
     193  // read the fits table bytes
     194  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape;
     195 
     196  // need to create and assign to flat-field correction
     197  GET_COLUMN(R,         "ra",           float);
     198  GET_COLUMN(D,         "dec",          float);
     199  GET_COLUMN(pQSO,      "p_QSO",        float);
     200  GET_COLUMN(pRRLyra,   "p_RRLyrae",    float);
     201  gfits_free_header (&theader);
     202  gfits_free_table  (&ftable);
     203
     204  MyStars *stars = NULL;
     205  ALLOCATE (stars, MyStars, Nrow);
     206
     207  for (i = 0; i < Nrow; i++) {
     208    stars[i].R      = R[i];         
     209    stars[i].D      = D[i];
     210    // stars[i].time   = 0;
     211    // stars[i].filter = 0;
     212    stars[i].myBit  = ID_OBJ_HERN_VARIABLE;
     213
     214    if (pQSO[i] >= MIN_QSO_P60_VALUE) {
     215      stars[i].myBit |= ID_OBJ_HERN_QSO_P60;
     216    }
     217    if (pQSO[i] >= MIN_QSO_P05_VALUE) {
     218      stars[i].myBit |= ID_OBJ_HERN_QSO_P05;
     219    }
     220    if (pRRLyra[i] >= MIN_RRLYRA_P60_VALUE) {
     221      stars[i].myBit |= ID_OBJ_HERN_RRL_P60;
     222    }
     223    if (pRRLyra[i] >= MIN_RRLYRA_P05_VALUE) {
     224      stars[i].myBit |= ID_OBJ_HERN_RRL_P05;
     225    }
     226    stars[i].flag   = FALSE;
     227    stars[i].found  = FALSE;
     228  }
     229
     230  FREE (R);
     231  FREE (D);
     232  FREE (pQSO);
     233  FREE (pRRLyra);
     234
     235  gfits_free_header (&header);
     236  gfits_free_matrix (&matrix);
     237  fclose (f);
     238 
     239  *Nstars = Nrow;
     240  return (stars);
     241
     242 escape:
     243  gfits_free_header (&header);
     244  gfits_free_matrix (&matrix);
     245  gfits_free_header (&theader);
     246  gfits_free_table  (&ftable);
     247
     248  fclose (f);
     249  return NULL;
     250}
     251
     252int expname_to_int (char *expname) {
     253 
     254  // XXX NOTE this does not handle gpc2 images
     255  myAssert (expname[ 0] == 'o', "invalid exposure name");
     256  myAssert (expname[ 5] == 'g', "invalid exposure name");
     257  myAssert (expname[10] == 'o', "invalid exposure name");
     258  int mjdExp = atoi(&expname[1]);
     259  int seqExp = atoi(&expname[6]);
     260  int value  = mjdExp * 10000 + seqExp;
     261  return value;
     262}
     263
     264void sort_exp (int *I, int *T, int *F, int *S, int N) {
     265
     266# define SWAPFUNC(A,B){ int itmp;       \
     267    itmp = I[A]; I[A] = I[B]; I[B] = itmp;      \
     268    itmp = T[A]; T[A] = T[B]; T[B] = itmp;      \
     269    itmp = F[A]; F[A] = F[B]; F[B] = itmp;      \
     270    itmp = S[A]; S[A] = S[B]; S[B] = itmp;      \
     271  }
     272# define COMPARE(A,B)(I[A] < I[B])
     273 
     274  OHANA_SORT (N, COMPARE, SWAPFUNC);
     275 
     276# undef SWAPFUNC
     277# undef COMPARE
     278 
     279}
     280
     281// find the entry by bisection:
     282int value_to_sequence (int *sequence, int Nsequence, int value) {
     283
     284  int Nlo = 0;
     285  int Nhi = Nsequence - 1;
     286
     287  if (value < sequence[Nlo]) return -1;
     288  if (value > sequence[Nhi]) return -1;
     289
     290  int N;
     291  while (Nhi - Nlo > 4) {
     292    N = 0.5*(Nlo + Nhi);
     293    if (sequence[N] < value) {
     294      Nlo = MAX(N, 0);
     295    } else {
     296      Nhi = MIN(N + 1, Nsequence - 1);
     297    }
     298  }
     299  // mySequence->value[Nlo] <  value
     300  // mySequence->value[Nhi] >= value
     301
     302  for (N = Nlo; N <= Nhi; N++) {
     303    if (sequence[N] == value) {
     304      return N;
     305    }
     306  }
     307  return -1;
     308}
     309

  • trunk/Ohana/src/addstar/src/setobjflags_remote_hosts.c

    r39496 r39525  
    2828void free_remote_hosts (void) {
    2929  if (!hosts) return;
     30
     31  for (int i = 0; i < Nhosts; i++) {
     32    if (hosts[i] == NULL) continue;
     33   
     34    HostInfo *host = hosts[i];
     35    free (host->hostname);
     36    free (host->pathname);
     37    FreeIOBuffer (&host->stdout);
     38    FreeIOBuffer (&host->stderr);
     39    free (host);
     40  }
    3041  free (hosts);
    3142}

  • trunk/Ohana/src/addstar/src/setobjflags_save_remote.c

    r39496 r39525  
    5959    fprintf (stderr, "command: %s\n", command);
    6060
     61    if (PARALLEL_MANUAL) {
     62      free (command);
     63      fprintf (stderr, "run the setobjflags_client command above. when done, hit return\n");
     64      getchar ();
     65      save_remote_host (host);
     66      return TRUE;
     67    }
     68
    6169    // launch the job on the remote machine (no handshake)
    6270    int errorInfo = 0;

  • trunk/Ohana/src/addstar/src/setobjflags_table.c

    r39496 r39525  
    22# include "setobjflags.h"
    33
    4 int setobjflags_table (SkyList *skylistInput, HostTable *hosts, char *filename, int mybit) {
     4int setobjflags_table (SkyList *skylistInput, HostTable *hosts) {
    55 
    6   int i, Nstars;
    7 
    8   fprintf (stderr, "loading %s, using bit 0x%08x\n", filename, mybit);
    9   MyStars *stars = setobjflags_loadfile (filename, mybit, &Nstars);
     6  int Nstars;
     7  MyStars *stars = setobjflags_loadfile (&Nstars);
    108  if (!stars) return FALSE;
    119
    … …  
    1513  // scan through the stars, loading the containing catalogs
    1614  // skip through table for unsaved stars
    17   for (i = 0; i < Nstars; i++) {
     15  for (int i = 0; i < Nstars; i++) {
    1816    if (stars[i].flag) continue;
    1917
Note: See TracChangeset for help on using the changeset viewer.