Index: trunk/Ohana/src/libdvo/Makefile =================================================================== --- trunk/Ohana/src/libdvo/Makefile (revision 36623) +++ trunk/Ohana/src/libdvo/Makefile (revision 36630) @@ -103,4 +103,5 @@ $(SRC)/ImageMetadata.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/RegionHostTable.$(ARCH).o \ $(SRC)/match_image.$(ARCH).o \ $(SRC)/db_utils.$(ARCH).o \ Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 36623) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 36630) @@ -285,4 +285,46 @@ } HostTable; +// A RegionHost processes data for some region in parallel with other regions +typedef struct RegionHostInfo { + double Rmin; // (Rmin,Rmax),(Dmin,Dmax) arehard RA,DEC boundaries of the + double Rmax; // region for which each host is responsible. A given host + double Dmin; // calibrates the images for which the fiducial point (center) + double Dmax; // lands in the region, and all objects in the region + + double RminCat; // (RminCat,RmaxCat),(DminCat,DmaxCat) are the region for which + double RmaxCat; // the catalogs need to be loaded : this is the outer bounds + double DminCat; // of the region containing all images completely + double DmaxCat; + + char *hostname; + + int hostID; // remove machine ID in SkyTable + int stdio[3]; // fd's for communication with the remote host + int pid; // remote process ID + int status; + IOBuffer stdout; + IOBuffer stderr; + + off_t Nimage; + off_t NIMAGE; + Image *image; + off_t *imseq; + + int *neighbors; // list of neighbor index values + int Nneighbors; // number of neighbors + char isNeighbor; // TRUE if I am a neighbor to the current region host +} RegionHostInfo; + +typedef struct { + double Rmin; + double Rmax; + double Dmin; + double Dmax; + + int Nhosts; + RegionHostInfo *hosts; + short *index; +} RegionHostTable; + // special-case function: CMF_PS1_V2 *gfits_table_get_CMF_PS1_V1_Alt (FTable *ftable, off_t *Ndata, char *swapped); @@ -385,4 +427,6 @@ uint32_t flags; int catID; + int objID; + int nOwn; } AverageTiny; @@ -410,4 +454,5 @@ short dYccd; short dRsys; + char myDet; } MeasureTiny; @@ -555,4 +600,5 @@ float *X; float *Y; + int *nOwn; // relastro uses this to count owned detections per object } Catalog; @@ -616,6 +662,4 @@ float PhotInst (Measure *measure); float PhotCat (Measure *measure); -float PhotAper (Measure *measure); -float PhotKron (Measure *measure); float PhotSys (Measure *measure, Average *average, SecFilt *secfilt); float PhotRel (Measure *measure, Average *average, SecFilt *secfilt); @@ -623,12 +667,24 @@ float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt); float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +float PhotAveErr (PhotCode *code, Average *average, SecFilt *secfilt); + +float PhotAperInst (Measure *measure); +float PhotAperCat (Measure *measure); +float PhotAperSys (Measure *measure, Average *average, SecFilt *secfilt); +float PhotAperRel (Measure *measure, Average *average, SecFilt *secfilt); +float PhotAperCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotAperRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); + +float PhotKronInst (Measure *measure); +float PhotKronCat (Measure *measure); +float PhotKronSys (Measure *measure, Average *average, SecFilt *secfilt); +float PhotKronRel (Measure *measure, Average *average, SecFilt *secfilt); +float PhotKronCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code); +float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt); +float PhotKronRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure); +float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt); + float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt); -float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt); - -float PhotAperInst (Measure *measure); -float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt); -float PhotKronInst (Measure *measure); -float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt); -float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt); float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt); @@ -855,3 +911,11 @@ void dvo_measureT_init (MeasureTiny *measure); +void InitRegionHosts (RegionHostInfo *hosts, int Nhosts, int NHOSTS); +void FreeRegionHosts (RegionHostInfo *hosts, int Nhosts); +void FreeRegionHostTable (RegionHostTable *table); +RegionHostTable *RegionHostTableLoad (char *catdir, char *rootname); +int RegionHostTableWaitJobs (RegionHostTable *regionHosts, char *file, int lineno); +int RegionHostTableWaitJobsGetIO (RegionHostTable *regionHosts, char *file, int lineno, int VERBOSE); +int RegionHostFindNeighbors (RegionHostTable *table, int Nhost); + # endif // DVO_H Index: trunk/Ohana/src/libdvo/include/dvodb.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvodb.h (revision 36623) +++ trunk/Ohana/src/libdvo/include/dvodb.h (revision 36630) @@ -23,31 +23,60 @@ /* magnitude types */ enum {MAG_NONE, - MAG_INST, - MAG_CAT, - MAG_APER, - MAG_APER_INST, - MAG_KRON, - MAG_KRON_INST, - MAG_KRON_ERR, - MAG_SYS, - MAG_REL, - MAG_CAL, + + // these magnitude types are from the measurement table (or derived therein) + // the following imply PSF magnitudes: + MAG_INST, // -2.5*log(DN) [ie, not DN/sec] + MAG_CAT, // MAG_INST + 2.5*log(exptime) + C_lambda + K_lambda*(airmass - 1) + MAG_SYS, // MAG_CAT + \sum_i A_i * (color_i - color_o) [color correction for measure photcode] + MAG_REL, // MAG_SYS - Mcal [specific zero point for image] + MAG_CAL, // MAG_REL + \sum_i A_i * (color_i - color_o) [color correction for average photcode] + MAG_ERR, // error on PSF MAG + + // the following imply APERTURE magnitudes: + MAG_APER_INST, // -2.5*log(DN) [ie, not DN/sec] + MAG_APER_CAT, // MAG_APER_INST + 2.5*log(exptime) + C_lambda + K_lambda*(airmass - 1) + MAG_APER_SYS, // MAG_APER_CAT + \sum_i A_i * (color_i - color_o) [color correction for measure photcode] + MAG_APER_REL, // MAG_APER_SYS - Mcal [specific zero point for image] + MAG_APER_CAL, // MAG_APER_REL + \sum_i A_i * (color_i - color_o) [color correction for average photcode] + MAG_APER_ERR, // error on APER MAG (if known, else psf error) + MAG_APER, // in mextract -> MAG_APER_REL, in avextract -> MAG_APER_AVE + + // the following imply KRON magnitudes: + MAG_KRON_INST, // -2.5*log(DN) [ie, not DN/sec] + MAG_KRON_CAT, // MAG_KRON_INST + 2.5*log(exptime) + C_lambda + K_lambda*(airmass - 1) + MAG_KRON_SYS, // MAG_KRON_CAT + \sum_i A_i * (color_i - color_o) [color correction for measure photcode] + MAG_KRON_REL, // MAG_KRON_SYS - Mcal [specific zero point for image] + MAG_KRON_CAL, // MAG_KRON_REL + \sum_i A_i * (color_i - color_o) [color correction for average photcode] + MAG_KRON_ERR, // error on KRON MAG + MAG_KRON, // in mextract -> MAG_KRON_REL, in avextract -> MAG_KRON_AVE + + // these magnitude types are from the average (secfilt) table: + // the following imply PSF magnitudes (or fluxes) MAG_AVE, MAG_REF, - MAG_ERR, MAG_STDEV, MAG_AVE_ERR, + MAG_CHISQ, + MAG_20, + MAG_80, + MAG_FLUX_PSF, + MAG_FLUX_PSF_ERR, + + // the following imply APERTURE magnitudes: + MAG_APER_AVE, + MAG_APER_REF, + + // the following imply KRON magnitudes or fluxes: + MAG_KRON_AVE, + MAG_KRON_REF, + MAG_FLUX_KRON, + MAG_FLUX_KRON_ERR, + + // other secfilt fields of interest MAG_PHOT_FLAGS, - MAG_CHISQ, MAG_NCODE, MAG_NPHOT, - MAG_20, - MAG_80, MAG_UC_DIST, MAG_STACK_DET_ID, - MAG_FLUX_PSF, - MAG_FLUX_PSF_ERR, - MAG_FLUX_KRON, - MAG_FLUX_KRON_ERR, }; Index: trunk/Ohana/src/libdvo/src/RegionHostTable.c =================================================================== --- trunk/Ohana/src/libdvo/src/RegionHostTable.c (revision 36630) +++ trunk/Ohana/src/libdvo/src/RegionHostTable.c (revision 36630) @@ -0,0 +1,497 @@ +# include "dvo.h" +# include +# include + +void InitRegionHosts (RegionHostInfo *hosts, int Nhosts, int NHOSTS) { + + int i; + for (i = Nhosts; i < NHOSTS; i++) { + hosts[i].Rmin = NAN; + hosts[i].Rmax = NAN; + hosts[i].Dmin = NAN; + hosts[i].Dmax = NAN; + + hosts[i].RminCat = NAN; + hosts[i].RmaxCat = NAN; + hosts[i].DminCat = NAN; + hosts[i].DmaxCat = NAN; + + hosts[i].hostname = NULL; + hosts[i].stdio[HOST_STDIN] = -1; + hosts[i].stdio[HOST_STDOUT] = -1; + hosts[i].stdio[HOST_STDERR] = -1; + hosts[i].pid = 0; + + hosts[i].image = NULL; + hosts[i].imseq = NULL; + hosts[i].Nimage = 0; + hosts[i].NIMAGE = 0; + + hosts[i].neighbors = NULL; + hosts[i].Nneighbors = 0; + hosts[i].isNeighbor = FALSE; + } + return; +} + +void FreeRegionHosts (RegionHostInfo *hosts, int Nhosts) { + + int i; + for (i = 0; i < Nhosts; i++) { + free (hosts[i].hostname); + } + free (hosts); + return; +} + +void FreeRegionHostTable (RegionHostTable *table) { + + if (!table) return; + if (table->hosts) { + FreeRegionHosts (table->hosts, table->Nhosts); + } + if (table->index) free (table->index); + free (table); + return; +} + +RegionHostTable *RegionHostTableLoad (char *catdir, char *rootname) { + + int i, Nline; + + char *filename = NULL; + + ALLOCATE (filename, char, strlen(catdir) + strlen(rootname) + 2); // one slash and one EOL + sprintf (filename, "%s/%s", catdir, rootname); + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "failed to open host table %s\n", filename); + free (filename); + return NULL; + } + + // simple format: ID hostname pathname + + int NHOSTS = 16; + int Nhosts = 0; + RegionHostInfo *hosts = NULL; + ALLOCATE (hosts, RegionHostInfo, NHOSTS); + InitRegionHosts (hosts, Nhosts, NHOSTS); + + int maxID = 0; + + RegionHostTable *table = NULL; + ALLOCATE (table, RegionHostTable, 1); + table->Rmin = 360; + table->Rmax = 0; + table->Dmin = +90; + table->Dmax = -90; + + for (Nline = 0; TRUE; Nline ++) { + int ID; + char tmphost[1024]; + char line[1024]; + + // XXXX use this for safety: int status = scan_line_maxlen (f, line, 1024); + int status = scan_line (f, line); + if (status == EOF) break; + + // find first non-whitespace char & skip commented lines + for (i = 0; OHANA_WHITESPACE (line[i]); i++); + if (line[i] == '#') continue; + if (line[i] == 0) continue; + + double Rmin, Rmax, Dmin, Dmax; + status = sscanf (line, "%d %1023s %lf %lf %lf %lf", &ID, tmphost, &Rmin, &Rmax, &Dmin, &Dmax); + if (status != 6) { + fprintf (stderr, "error reading line %d of region host table %s\n", Nline, filename); + FreeRegionHosts (hosts, Nhosts); + free (filename); + fclose (f); + return NULL; + } + + // check the validity of ID (0 < ID < MAX_SHORT) + + if (ID < 1) { + fprintf (stderr, "invalid host ID %d\n", ID); + exit (1); + } + if (ID > 255) { + fprintf (stderr, "invalid host ID %d\n", ID); + exit (1); + } + maxID = MAX(maxID, ID); + + hosts[Nhosts].hostID = ID; + hosts[Nhosts].hostname = strcreate(tmphost); + + InitIOBuffer (&hosts[Nhosts].stdout, 1000); + InitIOBuffer (&hosts[Nhosts].stderr, 1000); + + hosts[Nhosts].Rmin = Rmin; + hosts[Nhosts].Rmax = Rmax; + hosts[Nhosts].Dmin = Dmin; + hosts[Nhosts].Dmax = Dmax; + + hosts[Nhosts].RminCat = Rmin; + hosts[Nhosts].RmaxCat = Rmax; + hosts[Nhosts].DminCat = Dmin; + hosts[Nhosts].DmaxCat = Dmax; + + table->Rmin = MIN(Rmin, table->Rmin); + table->Rmax = MAX(Rmax, table->Rmax); + table->Dmin = MIN(Dmin, table->Dmin); + table->Dmax = MAX(Dmax, table->Dmax); + + // InitIOBuffer (&hosts[Nhosts].stdout, 1000); + // InitIOBuffer (&hosts[Nhosts].stderr, 1000); + + Nhosts ++; + if (Nhosts >= NHOSTS) { + NHOSTS += 16; + REALLOCATE (hosts, RegionHostInfo, NHOSTS); + InitRegionHosts (hosts, Nhosts, NHOSTS); + } + } + + table->Nhosts = Nhosts; + table->hosts = hosts; + + ALLOCATE (table->index, short, maxID + 1); + for (i = 0; i <= maxID; i++) table->index[i] = -1; + + for (i = 0; i < table->Nhosts; i++) { + if (table->index[table->hosts[i].hostID] != -1) { + fprintf (stderr, "error: duplicate hostID %d\n", table->hosts[i].hostID); + exit (1); + } + table->index[table->hosts[i].hostID] = i; + } + + free (filename); + fclose (f); + return table; +} + +int RegionHostFindNeighbors (RegionHostTable *table, int Nhost) { + + int i; + + // given a specific host (by table sequence), find all of its neighbors + // (eg, (host->Rmin == althost->Rmax && (host->Dmin <= althost->Dmax) && (host->Dmax >= althost->Dmin) + + RegionHostInfo *myhost = &table->hosts[Nhost]; + myAssert (!myhost->neighbors, "myhost neighbors already allocated"); + myAssert (!myhost->Nneighbors, "myhost Nneighbors not zero?"); + ALLOCATE (myhost->neighbors, int, 1); // always allocate 1 extra + + for (i = 0; i < table->Nhosts; i++) { + + if (i == Nhost) continue; + RegionHostInfo *altHost = &table->hosts[i]; + + int onBorder; + + onBorder = + (myhost->Rmin == altHost->Rmax) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Rmax == altHost->Rmin) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + // handle the 0,360 boundary (Rmin,Rmax are in the range 0,360, but 0 == 360) + onBorder = + (myhost->Rmin == 0.0) && + (altHost->Rmax == 360.0) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Rmax == 360.0) && + (altHost->Rmin == 0.0) && + (myhost->Dmin <= altHost->Dmax) && + (myhost->Dmax >= altHost->Dmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Dmin == altHost->Dmax) && + (myhost->Rmin <= altHost->Rmax) && + (myhost->Rmax >= altHost->Rmin); + if (onBorder) goto add_neighbor; + + onBorder = + (myhost->Dmax == altHost->Dmin) && + (myhost->Rmin <= altHost->Rmax) && + (myhost->Rmax >= altHost->Rmin); + if (onBorder) goto add_neighbor; + + continue; + + add_neighbor: + myhost->neighbors[myhost->Nneighbors] = i; + myhost->Nneighbors ++; + REALLOCATE (myhost->neighbors, int, myhost->Nneighbors + 1); + altHost->isNeighbor = TRUE; + } + return TRUE; +} + +// wait for all children to complete, report output to stdout +int RegionHostTableWaitJobs (RegionHostTable *table, char *file, int lineno) { + + int i; + + // we have launched table->Nhosts jobs; wait for all of them to complete... + // if one (N) failed to launch, we will get an ECHILD error from the last (N) calls + int done = FALSE; + for (i = 0; !done && (i < table->Nhosts); i++) { + int status = 0; + + // XXX we'll need to pass in WNOHANG and keep retrying if we want to have a timeout... + int pid = waitpid (-1, &status, 0); + if (!pid) { + // this should only occur if we called waitpid with the WNOHANG option + fprintf (stderr, "programming error (1)? %s %d", file, lineno); + exit (2); + } + if (pid == -1) { + switch (errno) { + case ECHILD: + done = TRUE; + break; + default: + fprintf (stderr, "programming error (2)? %s %d", file, lineno); + exit (2); + } + } + + // when the host has finished, close the open sockets + + // find the host which has finished + int Nout, j; + int found = FALSE; + for (j = 0; j < table->Nhosts; j++) { + if (table->hosts[j].pid != pid) continue; + found = TRUE; + // check on the status of this and report any output? + fprintf (stderr, "job finished for %s (%d)\n", table->hosts[j].hostname, pid); + // read the stderr and stdout + IOBuffer buffer; + InitIOBuffer (&buffer, 100); + EmptyIOBuffer (&buffer, 100, table->hosts[j].stdio[HOST_STDOUT]); + fprintf (stderr, "--- stdout from %s ---\n", table->hosts[j].hostname); + Nout = write (STDOUT_FILENO, buffer.buffer, buffer.Nbuffer); + if (Nout != buffer.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + fprintf (stderr, "\n"); + + InitIOBuffer (&buffer, 100); + EmptyIOBuffer (&buffer, 100, table->hosts[j].stdio[HOST_STDERR]); + fprintf (stderr, "--- stderr from %s ---\n", table->hosts[j].hostname); + Nout = write (STDOUT_FILENO, buffer.buffer, buffer.Nbuffer); + if (Nout != buffer.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + fprintf (stderr, "\n"); + if (WIFEXITED(status)) { + fprintf (stderr, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + table->hosts[j].status = WEXITSTATUS(status); + if (table->hosts[j].status) { + fprintf (stderr, "job failed on %s\n", table->hosts[j].hostname); + continue; + } + } else { + table->hosts[j].status = -1; + fprintf (stderr, "job exited abnormally on %s\n", table->hosts[j].hostname); + continue; + } + } + if (!found) { + fprintf (stderr, "Programming error: failed to matched finished job to known host!\n"); + exit (2); + } + } + return TRUE; +} + + +// wait for all children to complete, report output to stdout +int RegionHostTableWaitJobsGetIO (RegionHostTable *table, char *file, int lineno, int VERBOSE) { + + // we have launched table->Nhosts jobs; wait for all of them to complete... + // if one (N) failed to launch, we will get an ECHILD error from the last (N) calls + + // we need to read any data waiting on stderr or stdout from these jobs, or the overfull + // buffers can cause a problem. we alternate between 'select' and 'waitpid' calls with + // timeouts for both + + // add all hosts' sockets to the fd_sets + fd_set rdSet, wtSet; + FD_ZERO (&rdSet); + FD_ZERO (&wtSet); + + // XXX can I set the fd_sets once, since I am not actually closing the fd's? + + int globalStatus = TRUE; + + int i; + int Nmax = 0; + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + FD_SET (table->hosts[i].stdio[HOST_STDIN], &wtSet); + Nmax = MAX (Nmax, table->hosts[i].stdio[HOST_STDIN]); + FD_SET (table->hosts[i].stdio[HOST_STDOUT], &rdSet); + Nmax = MAX (Nmax, table->hosts[i].stdio[HOST_STDOUT]); + FD_SET (table->hosts[i].stdio[HOST_STDERR], &rdSet); + Nmax = MAX (Nmax, table->hosts[i].stdio[HOST_STDERR]); + } + Nmax ++; + + // need the list of connected hosts for exit test below + int Nrunning = 0; + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + Nrunning ++; + } + + int Nfound = 0; + + // this loop has 2 chunks: (a) check for I/O + (b) check for jobs done + while (1) { + + // Wait up to 0.5 second for host to provide I/O + // timeout gets mucked: need to reset before each select + struct timeval timeout; + timeout.tv_sec = 10; + timeout.tv_usec = 500000; + + int status = select (Nmax, NULL, &wtSet, NULL, &timeout); + if (status == -1) { + perror("select()"); + exit (2); + } + + // we have some sockets to check, check sockets for all hosts + for (i = 0; (status > 0) && (i < table->Nhosts); i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + + if (FALSE && FD_ISSET (table->hosts[i].stdio[HOST_STDIN], &wtSet)) { + // this host is waiting for input : this is an error, so exit + fprintf (stderr, "host %s is waiting for input\n", table->hosts[i].hostname); + abort(); + } + + if ((table->hosts[i].stdio[HOST_STDOUT] > 0) && FD_ISSET (table->hosts[i].stdio[HOST_STDOUT], &rdSet)) { + // this host has waiting output : read to buffer, and dump if necessary + ReadtoIOBuffer (&table->hosts[i].stdout, table->hosts[i].stdio[HOST_STDOUT]); + // if (table->hosts[i].stdout.Nbuffer > 0x10000) { + if (table->hosts[i].stdout.Nbuffer > 0x1000) { + int printHead = VERBOSE || (table->hosts[i].stdout.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v1)\n", table->hosts[i].hostname, table->hosts[i].stdout.Nbuffer); + int Nout = write (STDOUT_FILENO, table->hosts[i].stdout.buffer, table->hosts[i].stdout.Nbuffer); + if (Nout != table->hosts[i].stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&table->hosts[i].stdout); + if (printHead) fprintf (stdout, "\n"); + } + } + + if ((table->hosts[i].stdio[HOST_STDERR] > 0) && FD_ISSET (table->hosts[i].stdio[HOST_STDERR], &rdSet)) { + // this host has waiting output : read to buffer, and dump if necessary + ReadtoIOBuffer (&table->hosts[i].stderr, table->hosts[i].stdio[HOST_STDERR]); + // if (table->hosts[i].stderr.Nbuffer > 0x10000) { + if (table->hosts[i].stderr.Nbuffer > 0x1000) { + int printHead = VERBOSE || (table->hosts[i].stderr.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v1)\n", table->hosts[i].hostname, table->hosts[i].stderr.Nbuffer); + int Nout = write (STDOUT_FILENO, table->hosts[i].stderr.buffer, table->hosts[i].stderr.Nbuffer); + if (Nout != table->hosts[i].stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&table->hosts[i].stderr); + if (printHead) fprintf (stdout, "\n"); + } + } + } + + // now check if any children have finished... + while (TRUE) { + int status = 0; + int pid = waitpid (-1, &status, WNOHANG); + if (!pid) { + // fprintf (stderr, "no hosts to harvest\n"); + usleep (500000); + break; // no outstanding jobs have finished + } + if ((pid == -1) && (errno == ECHILD)) goto escape; // no more jobs on which to wait + if ((pid == -1) && (errno != ECHILD)) { + fprintf (stderr, "programming error (2)? %s %d", file, lineno); + exit (2); + } + + // find the host which has finished + int found = FALSE; + for (i = 0; (i < table->Nhosts) && !found; i++) { + if (table->hosts[i].pid != pid) continue; + found = TRUE; + + RegionHostInfo *host = &table->hosts[i]; + + // check on the status of this and report any output? + if (VERBOSE) fprintf (stdout, "job finished for %s (%d)\n", host->hostname, pid); + + // read stdout + int printHead; + printHead = VERBOSE || (host->stdout.Nbuffer > 0); + EmptyIOBuffer (&host->stdout, 100, host->stdio[HOST_STDOUT]); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v2)\n", host->hostname, host->stdout.Nbuffer); + int Nout = write (STDOUT_FILENO, host->stdout.buffer, host->stdout.Nbuffer); + if (Nout != host->stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stdout); + if (printHead) fprintf (stdout, "\n"); + + // read stderr + printHead = VERBOSE || (host->stderr.Nbuffer > 0); + EmptyIOBuffer (&host->stderr, 100, host->stdio[HOST_STDERR]); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v2)\n", host->hostname, host->stderr.Nbuffer); + Nout = write (STDOUT_FILENO, host->stderr.buffer, host->stderr.Nbuffer); + if (Nout != host->stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stderr); + if (printHead) fprintf (stdout, "\n"); + + if (WIFEXITED(status)) { + if (VERBOSE) fprintf (stdout, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + host->status = WEXITSTATUS(status); + if (host->status) { + fprintf (stdout, "job failed on %s\n", host->hostname); + globalStatus = FALSE; + } + } else { + host->status = -1; + fprintf (stdout, "job exited abnormally on %s\n", host->hostname); + globalStatus = FALSE; + continue; + } + } + if (!found) { + fprintf (stderr, "Programming error: failed to matched finished job to known host!\n"); + exit (2); + } + Nfound ++; + if (Nfound == Nrunning) goto escape; // we've harvested all jobs + } + } + +escape: + + // close all opened connections + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i].pid) continue; // any unconnected hosts should be skipped + close (table->hosts[i].stdio[HOST_STDIN]); + close (table->hosts[i].stdio[HOST_STDOUT]); + close (table->hosts[i].stdio[HOST_STDERR]); + } + + return globalStatus; +} Index: trunk/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 36630) @@ -211,5 +211,5 @@ } } else { - value.Flt = PhotdM (field->photcode, average, secfilt); + value.Flt = PhotAveErr (field->photcode, average, secfilt); } break; @@ -244,9 +244,19 @@ case MAG_APER: + case MAG_APER_AVE: value.Flt = PhotAperAve (field->photcode, average, secfilt); break; + case MAG_APER_REF: + value.Flt = PhotAperRef (field->photcode, average, secfilt, measure); + break; + case MAG_KRON: + case MAG_KRON_AVE: value.Flt = PhotKronAve (field->photcode, average, secfilt); + break; + + case MAG_KRON_REF: + value.Flt = PhotKronRef (field->photcode, average, secfilt, measure); break; @@ -283,5 +293,5 @@ break; case AVE_dMAG: - value.Flt = PhotdM (field->photcode, average, secfilt); + value.Flt = PhotAveErr (field->photcode, average, secfilt); break; case AVE_Xm: Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 36630) @@ -123,10 +123,20 @@ case MAG_REL: case MAG_CAL: - case MAG_APER: + case MAG_ERR: + case MAG_APER_INST: - case MAG_KRON: + case MAG_APER_CAT: + case MAG_APER_SYS: + case MAG_APER_REL: + case MAG_APER_CAL: + case MAG_APER_ERR: + case MAG_KRON_INST: + case MAG_KRON_CAT: + case MAG_KRON_SYS: + case MAG_KRON_REL: + case MAG_KRON_CAL: case MAG_KRON_ERR: - case MAG_ERR: + case MAG_PHOT_FLAGS: equiv = myEquiv; @@ -137,12 +147,23 @@ case MAG_AVE: case MAG_REF: + case MAG_STDEV: + case MAG_AVE_ERR: case MAG_CHISQ: - case MAG_AVE_ERR: + case MAG_FLUX_PSF: + case MAG_FLUX_PSF_ERR: + + case MAG_APER_AVE: + + case MAG_KRON_AVE: + case MAG_FLUX_KRON: + case MAG_FLUX_KRON_ERR: + + case MAG_20: + case MAG_80: + case MAG_UC_DIST: + case MAG_STACK_DET_ID: + case MAG_NCODE: case MAG_NPHOT: - case MAG_FLUX_PSF: - case MAG_FLUX_PSF_ERR: - case MAG_FLUX_KRON: - case MAG_FLUX_KRON_ERR: equiv = field->photcode; goto valid_photcode; @@ -178,40 +199,71 @@ value.Flt = PhotRef (equiv, average, secfilt, measure); break; - case MAG_APER: - value.Flt = PhotAper (measure); - break; + case MAG_ERR: + value.Flt = measure[0].dM; + break; + case MAG_AVE_ERR: + value.Flt = PhotAveErr (equiv, average, secfilt); + break; + case MAG_APER_INST: value.Flt = PhotAperInst (measure); break; - case MAG_KRON: - value.Flt = PhotKron (measure); - break; + case MAG_APER_CAT: + value.Flt = PhotAperCat (measure); + break; + case MAG_APER_SYS: + value.Flt = PhotAperSys (measure, average, secfilt); + break; + case MAG_APER: + case MAG_APER_REL: + value.Flt = PhotAperRel (measure, average, secfilt); + break; + case MAG_APER_CAL: + value.Flt = PhotAperCal (measure, average, secfilt, measure, equiv); + break; + case MAG_APER_AVE: + value.Flt = PhotAperAve (equiv, average, secfilt); + break; + case MAG_APER_REF: + value.Flt = PhotAperRef (equiv, average, secfilt, measure); + break; + case MAG_APER_ERR: + value.Flt = measure[0].dM; + break; + case MAG_KRON_INST: value.Flt = PhotKronInst (measure); break; + case MAG_KRON_CAT: + value.Flt = PhotKronCat (measure); + break; + case MAG_KRON_SYS: + value.Flt = PhotKronSys (measure, average, secfilt); + break; + case MAG_KRON: + case MAG_KRON_REL: + value.Flt = PhotKronRel (measure, average, secfilt); + break; + case MAG_KRON_CAL: + value.Flt = PhotKronCal (measure, average, secfilt, measure, equiv); + break; + case MAG_KRON_AVE: + value.Flt = PhotKronAve (equiv, average, secfilt); + break; + case MAG_KRON_REF: + value.Flt = PhotKronRef (equiv, average, secfilt, measure); + break; case MAG_KRON_ERR: - value.Flt = measure[0].dMkron; - break; - case MAG_ERR: - value.Flt = measure[0].dM; - break; - case MAG_AVE_ERR: - value.Flt = PhotdM (equiv, average, secfilt); - break; + value.Flt = measure[0].dMkron; + break; + case MAG_PHOT_FLAGS: if ((field->photcode->type == PHOT_REF) || (field->photcode->type == PHOT_DEP)) { value.Int = measure[0].photFlags; - } else { - value.Int = 0; } break; - case MAG_CHISQ: - value.Flt = PhotXm (equiv, average, secfilt); - break; + case MAG_NCODE: - if (equiv == NULL) { - value.Int = 0; - break; - } + if (equiv == NULL) break; Nsec = GetPhotcodeNsec (equiv->code); if (Nsec == -1) break; @@ -219,12 +271,26 @@ break; case MAG_NPHOT: - if (equiv == NULL) { - value.Int = 0; - break; - } + if (equiv == NULL) break; Nsec = GetPhotcodeNsec (equiv->code); if (Nsec == -1) break; value.Int = secfilt[Nsec].Nused; break; + + case MAG_CHISQ: + value.Flt = PhotXm (equiv, average, secfilt); + break; + case MAG_STDEV: + value.Flt = PhotMstdev (equiv, average, secfilt); + break; + case MAG_20: + value.Flt = PhotM20 (equiv, average, secfilt); + break; + case MAG_80: + value.Flt = PhotM80 (equiv, average, secfilt); + break; + case MAG_UC_DIST: + value.Flt = PhotUCdist (equiv, average, secfilt); + break; + case MAG_FLUX_PSF: value.Flt = PhotAveFluxPSF (field->photcode, average, secfilt); @@ -531,7 +597,8 @@ fieldc = MatchFieldMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); // fieldc = MatchField (measure[0].t, measure[0].photcode); + fieldc = MatchFieldMetadata (measure[0].imageID); } if (fieldc == NULL) break; @@ -547,7 +614,8 @@ fieldc = MatchFieldMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); // fieldc = MatchField (measure[0].t, measure[0].photcode); + fieldc = MatchFieldMetadata (measure[0].imageID); } if (fieldc == NULL) break; @@ -564,7 +632,8 @@ mosaic = MatchMosaicMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); + // mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + mosaic = MatchMosaicMetadata (measure[0].imageID); } if (mosaic == NULL) break; @@ -580,7 +649,8 @@ mosaic = MatchMosaicMetadata (measure[0].imageID); } else { - fprintf (stderr, "non-parallel Xmos broken\n"); - abort(); - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + // fprintf (stderr, "non-parallel Xmos broken\n"); + // abort(); + // mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + mosaic = MatchMosaicMetadata (measure[0].imageID); } if (mosaic == NULL) break; Index: trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbFields.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dbFields.c (revision 36630) @@ -30,4 +30,5 @@ int GetMagMode (char *string) { + // these all imply PSF mags: if (!strcasecmp (string, "inst")) return (MAG_INST); if (!strcasecmp (string, "cat")) return (MAG_CAT); @@ -37,14 +38,45 @@ if (!strcasecmp (string, "ave")) return (MAG_AVE); if (!strcasecmp (string, "ref")) return (MAG_REF); + if (!strcasecmp (string, "err")) return (MAG_ERR); + if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); + + // these are explicit PSF mags: + if (!strcasecmp (string, "psf_inst")) return (MAG_INST); + if (!strcasecmp (string, "psf_cat")) return (MAG_CAT); + if (!strcasecmp (string, "psf_sys")) return (MAG_SYS); + if (!strcasecmp (string, "psf_rel")) return (MAG_REL); + if (!strcasecmp (string, "psf_cal")) return (MAG_CAL); + if (!strcasecmp (string, "psf_ave")) return (MAG_AVE); + if (!strcasecmp (string, "psf_ref")) return (MAG_REF); + if (!strcasecmp (string, "psf_err")) return (MAG_ERR); + if (!strcasecmp (string, "psf_aveerr")) return (MAG_AVE_ERR); + + // these are aper mags: if (!strcasecmp (string, "ap")) return (MAG_APER); if (!strcasecmp (string, "aper")) return (MAG_APER); if (!strcasecmp (string, "aperinst")) return (MAG_APER_INST); if (!strcasecmp (string, "aper_inst")) return (MAG_APER_INST); + if (!strcasecmp (string, "aper_cat")) return (MAG_APER_CAT); + if (!strcasecmp (string, "aper_sys")) return (MAG_APER_SYS); + if (!strcasecmp (string, "aper_rel")) return (MAG_APER_REL); + if (!strcasecmp (string, "aper_cal")) return (MAG_APER_CAL); + if (!strcasecmp (string, "aper_ave")) return (MAG_APER_AVE); + if (!strcasecmp (string, "aper_ref")) return (MAG_APER_REF); + if (!strcasecmp (string, "aper_err")) return (MAG_APER_ERR); + + // these are kron mags: if (!strcasecmp (string, "kron")) return (MAG_KRON); + if (!strcasecmp (string, "kron_inst")) return (MAG_KRON_INST); + if (!strcasecmp (string, "kron_cat")) return (MAG_KRON_CAT); + if (!strcasecmp (string, "kron_sys")) return (MAG_KRON_SYS); + if (!strcasecmp (string, "kron_rel")) return (MAG_KRON_REL); + if (!strcasecmp (string, "kron_cal")) return (MAG_KRON_CAL); + if (!strcasecmp (string, "kron_ave")) return (MAG_KRON_AVE); + if (!strcasecmp (string, "kron_ref")) return (MAG_KRON_REF); + if (!strcasecmp (string, "kron_err")) return (MAG_KRON_ERR); if (!strcasecmp (string, "kroninst")) return (MAG_KRON_INST); - if (!strcasecmp (string, "kron_inst")) return (MAG_KRON_INST); if (!strcasecmp (string, "kronerr")) return (MAG_KRON_ERR); - if (!strcasecmp (string, "err")) return (MAG_ERR); - if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); + + // other fields: if (!strcasecmp (string, "photflags")) return (MAG_PHOT_FLAGS); if (!strcasecmp (string, "flags")) return (MAG_PHOT_FLAGS); Index: trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 36630) @@ -119,4 +119,5 @@ average->measureOffset = -1; average->catID = 0; + average->nOwn = 0; } @@ -242,4 +243,5 @@ measure->dYccd = 0; measure->dRsys = 0; + measure->myDet = FALSE; } Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 36630) @@ -549,4 +549,5 @@ } +// XXX note that there are 2 bad versions of PS1_V4 : the other one does have have PSF_QF_PERFECT, but has PAD instead int gfits_convert_Measure_PS1_V4alt (Measure_PS1_V4alt *data, off_t size, off_t nitems) { @@ -708,5 +709,5 @@ # ifdef BYTE_SWAP byte = (unsigned char *) data; - for (i = 0; i < nitems; i++, byte += 129) { + for (i = 0; i < nitems; i++, byte += 120) { /** BYTE SWAP **/ SWAP_DBLE (0); // RA @@ -754,5 +755,5 @@ Ncols = ftable[0].header[0].Naxis[0]; if (Ncols != 120) { - fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, 176); + fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, 120); return NULL; } Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 36630) @@ -169,4 +169,5 @@ } +// returns Nsec if code is PRI/SEC, else -1 int GetPhotcodeNsec (int code) { @@ -259,76 +260,4 @@ } -float PhotAper (Measure *measure) { - - int Np; - float Mcat; - PhotCode *code; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Mcat = measure[0].Map; - return (Mcat); - } - code = &photcodes[0].code[Np]; - Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; - - return (Mcat); -} - -float PhotAperInst (Measure *measure) { - - int Np; - float Minst; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Minst = measure[0].Map; - return (Minst); - } - Minst = measure[0].Map - measure[0].dt - ZERO_POINT; - - return (Minst); -} - -float PhotKron (Measure *measure) { - - int Np; - float Mcat; - PhotCode *code; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Mcat = measure[0].Mkron; - return (Mcat); - } - code = &photcodes[0].code[Np]; - Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; - - return (Mcat); -} - -float PhotKronInst (Measure *measure) { - - int Np; - float Minst; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - Minst = measure[0].Mkron; - return (Minst); - } - Minst = measure[0].Mkron - measure[0].dt - ZERO_POINT; - - return (Minst); -} - float PhotSys (Measure *measure, Average *average, SecFilt *secfilt) { @@ -362,23 +291,4 @@ } -float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt) { - - int Np; - float ZP; - PhotCode *code; - - Np = photcodes[0].hashcode[measure[0].photcode]; - if (Np == -1) return (NAN); - - if (photcodes[0].code[Np].type == PHOT_REF) { - ZP = 0.0; - return (ZP); - } - code = &photcodes[0].code[Np]; - ZP = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; - - return (ZP); -} - float PhotRel (Measure *measure, Average *average, SecFilt *secfilt) { @@ -444,4 +354,435 @@ Mcal += Mcol; return (Mcal); +} + +/***/ +float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Mave; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].M; + return (Mave); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].M; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +float PhotAveErr (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dM; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dM = (Ns == -1) ? NAN : secfilt[Ns].dM; + return (dM); +} + +/************************************** APERTURE-related Magnitudes ***********************************************/ + +float PhotAperInst (Measure *measure) { + + int Np; + float Minst; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Minst = measure[0].Map; + return (Minst); + } + Minst = measure[0].Map - measure[0].dt - ZERO_POINT; + + return (Minst); +} + +float PhotAperCat (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +float PhotAperSys (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Msys, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mcat); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Msys = Mcat + Mcol; + return (Msys); +} + +float PhotAperRel (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Mrel, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mrel = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mrel); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mrel += Mcol; + return (Mrel); +} + +/* return calibrated magnitude from measure for given photcode */ +float PhotAperCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + int i, Np; + float Mcal, Mrel, Mcol, mc, Mc; + + if (code == NULL) return NAN; + + /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */ + Np = photcodes[0].hashcode[thisone[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mrel = thisone[0].Map; + return (Mrel); + } + if (code[0].code != photcodes[0].code[Np].equiv) return (NAN); + + Mcal = PhotAperRel (thisone, average, secfilt) + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mcal); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Mcal += Mcol; + return (Mcal); +} + +float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Maper; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Maper = (Ns == -1) ? NAN : secfilt[Ns].Map; + return (Maper); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotAperRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].Map; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +/************************************** KRON-related Magnitudes ***********************************************/ + +float PhotKronInst (Measure *measure) { + + int Np; + float Minst; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Minst = measure[0].Mkron; + return (Minst); + } + Minst = measure[0].Mkron - measure[0].dt - ZERO_POINT; + + return (Minst); +} + +float PhotKronCat (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Mkron; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +float PhotKronSys (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Msys, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Msys = measure[0].Mkron; + return (Msys); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mcat); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Msys = Mcat + Mcol; + return (Msys); +} + +float PhotKronRel (Measure *measure, Average *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Mrel, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Mkron; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mrel = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCode (average, secfilt, NULL, code); + if (isnan(mc)) return (Mrel); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mrel += Mcol; + return (Mrel); +} + +/* return calibrated magnitude from measure for given photcode */ +float PhotKronCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) { + + int i, Np; + float Mcal, Mrel, Mcol, mc, Mc; + + if (code == NULL) return NAN; + + /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */ + Np = photcodes[0].hashcode[thisone[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mrel = thisone[0].Mkron; + return (Mrel); + } + if (code[0].code != photcodes[0].code[Np].equiv) return (NAN); + + Mcal = PhotKronRel (thisone, average, secfilt) + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mcal); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Mcal += Mcol; + return (Mcal); +} + +float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Mkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mkron = (Ns == -1) ? NAN : secfilt[Ns].Mkron; + return (Mkron); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotKronRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].Mkron; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCode (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dMkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dMkron = (Ns == -1) ? NAN : secfilt[Ns].dMkron; + return (dMkron); +} + +/********************* other support functions ********************************/ + +float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt) { + + int Np; + float ZP; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + ZP = 0.0; + return (ZP); + } + code = &photcodes[0].code[Np]; + ZP = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (ZP); } @@ -496,56 +837,4 @@ } -/* return calibrated magnitude from average/secfilt for given photcode */ -float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) { - - int i, Ns; - float Mave, Mref, Mcol, mc; - double Mc; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mave = (Ns == -1) ? NAN : secfilt[Ns].M; - Mref = Mave + SCALE*code[0].C; - - mc = PhotColorForCode (average, secfilt, measure, code); - if (isnan(mc)) return (Mref); - mc = mc - SCALE*code[0].dX; - - Mc = mc; - Mcol = 0; - for (i = 0; i < code[0].Nc; i++) { - Mcol += code[0].X[i]*Mc; - Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ - } - Mref += Mcol; - return (Mref); -} - -/***/ -float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Mave; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mave = (Ns == -1) ? NAN : secfilt[Ns].M; - return (Mave); -} - -float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float dM; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - dM = (Ns == -1) ? NAN : secfilt[Ns].dM; - return (dM); -} - float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt) { @@ -558,40 +847,4 @@ Mstdev = (Ns == -1) ? NAN : 0.001*secfilt[Ns].Mstdev; return (Mstdev); -} - -float PhotAperAve (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Maper; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Maper = (Ns == -1) ? NAN : secfilt[Ns].Map; - return (Maper); -} - -float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Mkron; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mkron = (Ns == -1) ? NAN : secfilt[Ns].Mkron; - return (Mkron); -} - -float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float dMkron; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - dMkron = (Ns == -1) ? NAN : secfilt[Ns].dMkron; - return (dMkron); } Index: trunk/Ohana/src/libdvo/src/dvo_tiny_values.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 36623) +++ trunk/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 36630) @@ -8,4 +8,6 @@ averageT[0].measureOffset = average[0].measureOffset; averageT[0].catID = average[0].catID; + averageT[0].objID = average[0].objID; + averageT[0].nOwn = -1; return; } @@ -33,4 +35,5 @@ measureT[0].dYccd = measure[0].dYccd; measureT[0].dRsys = measure[0].dRsys; + measureT[0].myDet = 0; return ; } Index: trunk/Ohana/src/libfits/table/F_get_column.c =================================================================== --- trunk/Ohana/src/libfits/table/F_get_column.c (revision 36623) +++ trunk/Ohana/src/libfits/table/F_get_column.c (revision 36630) @@ -174,118 +174,13 @@ int gfits_get_bintable_column (Header *header, FTable *table, char *label, void **data) { - off_t Nx, Ny; - int i, N, Nfields, Nval, Nbytes, Nstart, Nv, Nb; - char tlabel[80], field[80], format[80], type[16], tmpline[16]; - char *Pin, *Pout, *array; - double Bscale, Bzero; - - if (label == (char *) NULL) return (FALSE); - if (label[0] == 0) return (FALSE); - - /* find label in header */ - tlabel[0] = 0; - gfits_scan (header, "TFIELDS", "%d", 1, &Nfields); - for (i = 1; strcasecmp (label, tlabel) && (i < Nfields + 1); i++) { - sprintf (field, "TTYPE%d", i); - gfits_scan (header, field, "%s", 1, tlabel); - } - if (strcasecmp (label, tlabel)) return (FALSE); - N = i - 1; - - Bscale = 1; - Bzero = 0; - - /* interpret format */ - sprintf (field, "TSCAL%d", N); - gfits_scan (header, field, "%lf", 1, &Bscale); - sprintf (field, "TZERO%d", N); - gfits_scan (header, field, "%lf", 1, &Bzero); - sprintf (field, "TFORM%d", N); - gfits_scan (header, field, "%s", 1, format); - - if (!gfits_bintable_format (format, type, &Nval, &Nbytes)) return (FALSE); - - /* check existing table dimensions */ - gfits_scan (header, "NAXIS1", OFF_T_FMT, 1, &Nx); - gfits_scan (header, "NAXIS2", OFF_T_FMT, 1, &Ny); - - /* scan columns to find insert point */ - Nstart = 0; - for (i = 1; i < N; i++) { - sprintf (field, "TFORM%d", i); - gfits_scan (header, field, "%s", 1, format); - gfits_bintable_format (format, tmpline, &Nv, &Nb); - Nstart += Nv*Nb; - } - - /* extract bytes from table into temporary array */ - ALLOCATE (array, char, Nbytes*Nval*Ny); - Pin = table[0].buffer + Nstart; - Pout = array; - for (i = 0; i < Ny; i++, Pin += Nx, Pout += Nval*Nbytes) { - memcpy (Pout, Pin, Nval*Nbytes); - } - - /* convert data in-situ with correct type, byte swap and Bzero/Bscale */ - Pin = array; - Pout = array; - if (!strcmp (type, "char")) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { - *(char *)Pout = *(char *)Pin*Bscale + Bzero; - } - } - if (!strcmp (type, "short")) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_BYTE; -# endif - *(short *)Pout = *(short *)Pin*Bscale + Bzero; - } - } - if (!strcmp (type, "int")) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_WORD; -# endif - *(int *)Pout = *(int *)Pin*Bscale + Bzero; - } - } - if (!strcmp (type, "int64_t")) { - // XXX 64 bit int operations with Bzero & Bscale are inaccurate even with doubles - if ((Bzero == 0.0) && (Bscale == 1.0)) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - *(int64_t *)Pout = *(int64_t *)Pin; - } - } else { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero; - } - } - } - if (!strcmp (type, "float")) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_WORD; -# endif - *(float *)Pout = *(float *)Pin*Bscale + Bzero; - } - } - if (!strcmp (type, "double")) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - *(double *)Pout = *(double *)Pin*Bscale + Bzero; - } - } + char type[16]; + off_t Nrow; + int Ncol; + + char *array = gfits_get_bintable_column_data (header, table, label, type, &Nrow, &Ncol); + if (array == NULL) return (FALSE); *data = array; - return (TRUE); + return TRUE; } Index: trunk/Ohana/src/opihi/dvo/photometry.c =================================================================== --- trunk/Ohana/src/opihi/dvo/photometry.c (revision 36623) +++ trunk/Ohana/src/opihi/dvo/photometry.c (revision 36630) @@ -662,5 +662,5 @@ break; case AVE_dMAG: - value = PhotdM (code, average, secfilt); + value = PhotAveErr (code, average, secfilt); break; case AVE_Xm: @@ -865,5 +865,5 @@ /* for ErrSelect, check average errors */ if (ErrSelect) { - dM = PhotdM (code, average, secfilt); + dM = PhotAveErr (code, average, secfilt); if (dM > ErrValue) return (FALSE); } Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 36623) +++ trunk/Ohana/src/relastro/Makefile (revision 36630) @@ -69,4 +69,16 @@ $(SRC)/resort_catalog.$(ARCH).o \ $(SRC)/BrightCatalog.$(ARCH).o \ +$(SRC)/assign_images.$(ARCH).o \ +$(SRC)/launch_region_hosts.$(ARCH).o \ +$(SRC)/relastro_parallel_images.$(ARCH).o \ +$(SRC)/relastro_parallel_regions.$(ARCH).o \ +$(SRC)/MeanPosIO.$(ARCH).o \ +$(SRC)/share_mean_pos.$(ARCH).o \ +$(SRC)/share_images_pos.$(ARCH).o \ +$(SRC)/ImagePosIO.$(ARCH).o \ +$(SRC)/ImageTable.$(ARCH).o \ +$(SRC)/markObjects.$(ARCH).o \ +$(SRC)/indexCatalogs.$(ARCH).o \ +$(SRC)/syncfile.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 36623) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 36630) @@ -6,4 +6,13 @@ # include +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); + // choose off_t or int depending on full-scale relphot analysis resources // # define IDX_T off_t @@ -20,5 +29,7 @@ typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode; -typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS, TARGET_HPM} FitTarget; +typedef enum {OP_NONE, OP_IMAGES, OP_HIGH_SPEED, OP_MERGE_SOURCE, OP_UPDATE_OBJECTS, OP_UPDATE_OFFSETS, OP_LOAD_OBJECTS, OP_HPM, OP_PARALLEL_REGIONS, OP_PARALLEL_IMAGES} RelastroOp; + +typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget; typedef enum { @@ -30,4 +41,20 @@ MARK_BIG_OFFSET = 0x0010, } MeasurementMask; + +typedef struct { + float R; + float D; + unsigned int objID; + unsigned int catID; +} MeanPos; + +typedef struct { + Coords coords; + float dXpixSys; + float dYpixSys; + unsigned int imageID; + int nFitAstrom; + int flags; +} ImagePos; typedef struct { @@ -93,4 +120,5 @@ unsigned int start; unsigned int stop; + off_t myImage; float Mcal; float dMcal; @@ -113,10 +141,4 @@ } StatType; -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } - /* global variables set in parameter file */ # define DVO_MAX_PATH 1024 @@ -129,4 +151,10 @@ char SKY_TABLE[DVO_MAX_PATH]; int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ + +// globals for parallel region operations +char *REGION_FILE; +char *IMAGE_TABLE; +int REGION_HOST_ID; +int PARALLEL_REGIONS_MANUAL; int HOST_ID; @@ -208,4 +236,5 @@ FitMode FIT_MODE; +RelastroOp RELASTRO_OP; FitTarget FIT_TARGET; @@ -450,2 +479,44 @@ int hpm_catalogs_parallel (SkyList *skylist); int hpm_objects (SkyRegion *region, Catalog *catalog); + +int strextend (char *input, char *format,...); +int launch_region_hosts (RegionHostTable *regionHosts); + +int assign_images (FITS_DB *db, RegionHostTable *regionHosts); +int select_images_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage); +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid); +int calculate_host_image_bounds (RegionHostTable *regionHosts); +int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc); + +int relastro_parallel_regions (); +int relastro_parallel_images (); + +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname); +int check_sync_file (char *filename, int nloop); +int clear_sync_file (char *filename); +int update_sync_file (char *filename, int nloop); + +int share_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int slurp_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int set_mean_pos (MeanPos *meanpos, Average *average); +MeanPos *merge_mean_pos (MeanPos *target, int *ntarget, MeanPos *source, int Nsource); + +int MeanPosSave(char *filename, MeanPos *meanpos, off_t Nmeanpos); +MeanPos *MeanPosLoad(char *filename, off_t *nmeanpos); + +int indexCatalogs (Catalog *catalog, int Ncatalog); +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq); + +int markObjects (Catalog *catalog, int Ncatalog); + +int ImagePosSave(char *filename, ImagePos *image_pos, off_t Nimage_pos); +ImagePos *ImagePosLoad(char *filename, off_t *nimage_pos); + +int share_image_pos (RegionHostTable *regionHosts, int nloop); +int slurp_image_pos (RegionHostTable *regionHosts, int nloop); +ImagePos *merge_image_pos (ImagePos *target, int *ntarget, ImagePos *source, int Nsource); +int set_image_pos (ImagePos *image_pos, Image *image); + +Image *ImageTableLoad(char *filename, off_t *nimage); +int ImageTableSave (char *filename, Image *images, off_t Nimages); +int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 36623) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 36630) @@ -18,6 +18,7 @@ if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); - GetConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM); // exclude measurements on this basis - GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW); + // set defaults for all of these if they are not used by parallel / remote clients + if (!ScanConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM)) SIGMA_LIM = 0.01; + if (!ScanConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW)) SRC_MEAS_TOOFEW = 3; if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NITER", "%d", 0, &IMFIT_CLIP_NITER)) IMFIT_CLIP_NITER = 3; @@ -25,18 +26,14 @@ if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01; - GetConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN); - GetConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW); - GetConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW); + if (!ScanConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN)) PM_DT_MIN = 0.25; + if (!ScanConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW)) PM_TOOFEW = 4; + if (!ScanConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW)) POS_TOOFEW = 1; + if (!ScanConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN)) PAR_FACTOR_MIN = 0.2; + if (!ScanConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP)) NX_MAP = 5; + if (!ScanConfig (config, "RELASTRO_MAP_NY", "%d", 0, &NY_MAP)) NY_MAP = 5; + if (!ScanConfig (config, "RELASTRO_DPOS_MAX", "%lf", 0, &DPOS_MAX)) DPOS_MAX = 6.0; + if (!ScanConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS)) ADDSTAR_RADIUS = 1.0; - GetConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN); - - GetConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP); - GetConfig (config, "RELASTRO_MAP_NY", "%d", 0, &NY_MAP); - GetConfig (config, "RELASTRO_DPOS_MAX", "%lf", 0, &DPOS_MAX); - GetConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS); - - if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS)) { - USE_FIXED_PIXCOORDS = FALSE; - } + if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS)) USE_FIXED_PIXCOORDS = FALSE; // force CATDIR to be absolute (so parallel mode will work) @@ -56,10 +53,6 @@ sprintf (ImageCat, "%s/Images.dat", CATDIR); - if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) { - SKY_DEPTH = 2; - } - if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) { - SKY_TABLE[0] = 0; - } + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) SKY_DEPTH = 2; + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) SKY_TABLE[0] = 0; if (*CATMODE == 0) strcpy (CATMODE, "RAW"); Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 36623) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 36630) @@ -226,6 +226,6 @@ Measure *measureBig = &catalog[cat].measure[meas]; int TESTPT = FALSE; - TESTPT |= (measureBig->imageID == CAT_ID_SRC) && (measureBig->detID == OBJ_ID_SRC); - TESTPT |= (measureBig->imageID == CAT_ID_DST) && (measureBig->detID == OBJ_ID_DST); + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measureBig->imageID == CAT_ID_SRC) && (measureBig->detID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measureBig->imageID == CAT_ID_DST) && (measureBig->detID == OBJ_ID_DST); if (TESTPT) { fprintf (stderr, "got test det\n"); @@ -462,6 +462,6 @@ int TESTPT = FALSE; - TESTPT |= (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC); - TESTPT |= (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST); + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST); if (TESTPT) { fprintf (stderr, "got test det\n"); Index: trunk/Ohana/src/relastro/src/ImagePosIO.c =================================================================== --- trunk/Ohana/src/relastro/src/ImagePosIO.c (revision 36630) +++ trunk/Ohana/src/relastro/src/ImagePosIO.c (revision 36630) @@ -0,0 +1,315 @@ +# include "relastro.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +ImagePos *ImagePosLoad(char *filename, off_t *nimage_pos) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage_pos = 0; + ImagePos *image_pos = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // we read the entire block of data, then extract the columns, then set the image structure values. + // I free the FITS table data after extracting the colums to avoid having 3 copies in memory. + + char type[16]; + + GET_COLUMN (crval1, "CRVAL1", double); + GET_COLUMN (crval2, "CRVAL2", double); + GET_COLUMN (crpix1, "CRPIX1", float); + GET_COLUMN (crpix2, "CRPIX2", float); + GET_COLUMN (cdelt1, "CDELT1", float); + GET_COLUMN (cdelt2, "CDELT2", float); + GET_COLUMN (pc1_1, "PC1_1", float); + GET_COLUMN (pc1_2, "PC1_2", float); + GET_COLUMN (pc2_1, "PC2_1", float); + GET_COLUMN (pc2_2, "PC2_2", float); + GET_COLUMN (polyterms, "POLYTERMS", float); // verify that we got 14 columns? + GET_COLUMN (ctype, "CTYPE", char); // verify that we got 15 columns? + GET_COLUMN (Npolyterms, "NPOLYTERMS", char); + GET_COLUMN (dXpixSys, "XPIX_SYS_ERR", float); + GET_COLUMN (dYpixSys, "YPIX_SYS_ERR", float); + GET_COLUMN (imageID, "ID", int); + GET_COLUMN (nFitAstrom, "NFIT", int); + GET_COLUMN (flags, "FLAGS", int); + + // free the memory associated with the FITS files + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + ALLOCATE (image_pos, ImagePos, Nrow); + for (i = 0; i < Nrow; i++) { + image_pos[i].coords.crval1 = crval1 [i]; + image_pos[i].coords.crval2 = crval2 [i]; + image_pos[i].coords.crpix1 = crpix1 [i]; + image_pos[i].coords.crpix2 = crpix2 [i]; + image_pos[i].coords.cdelt1 = cdelt1 [i]; + image_pos[i].coords.cdelt2 = cdelt2 [i]; + image_pos[i].coords.pc1_1 = pc1_1 [i]; + image_pos[i].coords.pc1_2 = pc1_2 [i]; + image_pos[i].coords.pc2_1 = pc2_1 [i]; + image_pos[i].coords.pc2_2 = pc2_2 [i]; + image_pos[i].coords.Npolyterms = Npolyterms[i]; + image_pos[i].dXpixSys = dXpixSys [i]; + image_pos[i].dYpixSys = dYpixSys [i]; + image_pos[i].imageID = imageID [i]; + image_pos[i].nFitAstrom = nFitAstrom[i]; + image_pos[i].flags = flags [i]; + + // polyterms and ctype are a bit different + memcpy (&image_pos[i].coords.polyterms, &polyterms[i*14], 14*sizeof(float)); + memcpy (&image_pos[i].coords.ctype , &ctype [i*15], 15*sizeof(char)); + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (crval1 ); + free (crval2 ); + free (crpix1 ); + free (crpix2 ); + free (cdelt1 ); + free (cdelt2 ); + free (pc1_1 ); + free (pc1_2 ); + free (pc2_1 ); + free (pc2_2 ); + free (polyterms ); + free (ctype ); + free (Npolyterms); + free (dXpixSys ); + free (dYpixSys ); + free (imageID ); + free (nFitAstrom); + free (flags ); + + *nimage_pos = Nrow; + return image_pos; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int ImagePosSave(char *filename, ImagePos *image_pos, off_t Nimage_pos) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "IMAGE_POS"); + + gfits_define_bintable_column (&theader, "D", "CRVAL1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "CRVAL2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CRPIX1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CRPIX2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CDELT1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CDELT2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC1_1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC1_2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC2_1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC2_2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "14E", "POLYTERMS", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "15A", "CTYPE", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "L", "NPOLYTERMS", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "XPIX_SYS_ERR", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "YPIX_SYS_ERR", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + double *crval1 ; + double *crval2 ; + float *crpix1 ; + float *crpix2 ; + float *cdelt1 ; + float *cdelt2 ; + float *pc1_1 ; + float *pc1_2 ; + float *pc2_1 ; + float *pc2_2 ; + float *polyterms ; + char *ctype ; + char *Npolyterms; + float *dXpixSys ; + float *dYpixSys ; + int *imageID ; + int *nFitAstrom; + int *flags ; + + // create intermediate storage arrays + ALLOCATE (crval1 , double, Nimage_pos); + ALLOCATE (crval2 , double, Nimage_pos); + ALLOCATE (crpix1 , float , Nimage_pos); + ALLOCATE (crpix2 , float , Nimage_pos); + ALLOCATE (cdelt1 , float , Nimage_pos); + ALLOCATE (cdelt2 , float , Nimage_pos); + ALLOCATE (pc1_1 , float , Nimage_pos); + ALLOCATE (pc1_2 , float , Nimage_pos); + ALLOCATE (pc2_1 , float , Nimage_pos); + ALLOCATE (pc2_2 , float , Nimage_pos); + ALLOCATE (polyterms , float , 14*Nimage_pos); + ALLOCATE (ctype , char , 15*Nimage_pos); + ALLOCATE (Npolyterms, char , Nimage_pos); + ALLOCATE (dXpixSys , float , Nimage_pos); + ALLOCATE (dYpixSys , float , Nimage_pos); + ALLOCATE (imageID , int , Nimage_pos); + ALLOCATE (nFitAstrom, int , Nimage_pos); + ALLOCATE (flags , int , Nimage_pos); + + // assign the storage arrays + for (i = 0; i < Nimage_pos; i++) { + crval1 [i] = image_pos[i].coords.crval1 ; + crval2 [i] = image_pos[i].coords.crval2 ; + crpix1 [i] = image_pos[i].coords.crpix1 ; + crpix2 [i] = image_pos[i].coords.crpix2 ; + cdelt1 [i] = image_pos[i].coords.cdelt1 ; + cdelt2 [i] = image_pos[i].coords.cdelt2 ; + pc1_1 [i] = image_pos[i].coords.pc1_1 ; + pc1_2 [i] = image_pos[i].coords.pc1_2 ; + pc2_1 [i] = image_pos[i].coords.pc2_1 ; + pc2_2 [i] = image_pos[i].coords.pc2_2 ; + Npolyterms[i] = image_pos[i].coords.Npolyterms; + dXpixSys [i] = image_pos[i].dXpixSys ; + dYpixSys [i] = image_pos[i].dYpixSys ; + imageID [i] = image_pos[i].imageID ; + nFitAstrom[i] = image_pos[i].nFitAstrom ; + flags [i] = image_pos[i].flags ; + + // polyterms and ctype are a bit different + memcpy (&polyterms[i*14], &image_pos[i].coords.polyterms, 14*sizeof(float)); + memcpy (&ctype [i*15], &image_pos[i].coords.ctype , 15*sizeof(char)); + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "CRVAL1", crval1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRVAL2", crval2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRPIX1", crpix1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRPIX2", crpix2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CDELT1", cdelt1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CDELT2", cdelt2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC1_1", pc1_1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC1_2", pc1_2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC2_1", pc2_1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC2_2", pc2_2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "POLYTERMS", polyterms , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CTYPE", ctype , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "NPOLYTERMS", Npolyterms, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "XPIX_SYS_ERR",dXpixSys , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "YPIX_SYS_ERR",dYpixSys , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitAstrom, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_pos); + + free (crval1 ); + free (crval2 ); + free (crpix1 ); + free (crpix2 ); + free (cdelt1 ); + free (cdelt2 ); + free (pc1_1 ); + free (pc1_2 ); + free (pc2_1 ); + free (pc2_2 ); + free (polyterms ); + free (ctype ); + free (Npolyterms); + free (dXpixSys ); + free (dYpixSys ); + free (imageID ); + free (nFitAstrom); + free (flags ); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image_pos file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for image_pos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for image_pos %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for image_pos %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for image_pos %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file image_pos %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file image_pos %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing image_pos file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/ImageTable.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageTable.c (revision 36630) +++ trunk/Ohana/src/relastro/src/ImageTable.c (revision 36630) @@ -0,0 +1,62 @@ +# include "relastro.h" + +Image *ImageTableLoad(char *filename, off_t *nimage) { + + int status; + off_t Nimage; + FITS_DB db; + + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, filename, 60.0, LCK_SOFT); + if (!status) { + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + // if the file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + Shutdown ("ERROR: failed to read images"); + } + + *nimage = Nimage; + return image; +} + +int ImageTableSave (char *filename, Image *images, off_t Nimages) { + + int status; + FITS_DB db; + + // setup image table format and lock + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&db, filename, 60.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + // load or create the image table + if (db.dbstate != LCK_EMPTY) { + Shutdown ("image catalog already exists %s", db.filename); + } + + dvo_image_create (&db, GetZeroPoint()); + + /* add the new images and save */ + dvo_image_addrows (&db, images, Nimages); + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/MeanPosIO.c =================================================================== --- trunk/Ohana/src/relastro/src/MeanPosIO.c (revision 36630) +++ trunk/Ohana/src/relastro/src/MeanPosIO.c (revision 36630) @@ -0,0 +1,183 @@ +# include "relastro.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +MeanPos *MeanPosLoad(char *filename, off_t *nmeanpos) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nmeanpos = 0; + MeanPos *meanpos = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // a bit annoying : we read the entire block of data, then extract the columns, then set the image structure values. + // this means I need 3 copies in memory at some point. ugh. + + char type[16]; + + GET_COLUMN (R, "RA", double); + GET_COLUMN (D, "DEC", double); + GET_COLUMN (objID, "OBJ_ID", int); + GET_COLUMN (catID, "CAT_ID", int); + + ALLOCATE (meanpos, MeanPos, Nrow); + for (i = 0; i < Nrow; i++) { + meanpos[i].R = R [i]; + meanpos[i].D = D [i]; + meanpos[i].objID = objID[i]; + meanpos[i].catID = catID[i]; + } + fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); + + free (R ); + free (D ); + free (objID); + free (catID); + + // free FITS table pieces... + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + *nmeanpos = Nrow; + return meanpos; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int MeanPosSave(char *filename, MeanPos *meanpos, off_t Nmeanpos) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "MEANPOS"); + + gfits_define_bintable_column (&theader, "D", "RA", "mean position, ra", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "mean position, dec", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ID", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ID", NULL, 1.0, 1.0*0x8000); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + double *R, *D; + unsigned int *objID, *catID; + + // create intermediate storage arrays + ALLOCATE (R, double, Nmeanpos); + ALLOCATE (D, double, Nmeanpos); + ALLOCATE (objID, unsigned int, Nmeanpos); + ALLOCATE (catID, unsigned int, Nmeanpos); + + // assign the storage arrays + for (i = 0; i < Nmeanpos; i++) { + R [i] = meanpos[i].R ; + D [i] = meanpos[i].D ; + objID[i] = meanpos[i].objID; + catID[i] = meanpos[i].catID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, Nmeanpos); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, Nmeanpos); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeanpos); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeanpos); + + free (R ); + free (D ); + free (objID); + free (catID); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for meanpos %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for meanpos %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file meanpos %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file meanpos %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing meanpos file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relastro/src/MosaicOps.c (revision 36623) +++ trunk/Ohana/src/relastro/src/MosaicOps.c (revision 36630) @@ -11,5 +11,5 @@ // list of mosaic associated with each image -static off_t Nmosaic_for_images; // number of images (for off_ternal checks) +static off_t Nmosaic_for_images; // number of images (for internal checks) static off_t *mosaic_for_images; // array of: image -> mosaic @@ -105,4 +105,5 @@ mosaic[Nmosaic].secz = image[i].secz; mosaic[Nmosaic].coords = image[i].coords; + mosaic[Nmosaic].myImage = i; // init the mosaic_own_images array data @@ -135,4 +136,8 @@ ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images); + // emit an error if we miss mosaics, but stop if we miss too many + int NmissMosaic = 0; + int NtestMosaic = 0; + /* now assign the WRP images to these mosaics */ for (i = 0; i < Nimage; i++) { @@ -140,4 +145,5 @@ if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue; + NtestMosaic ++; /* set image time range */ @@ -147,5 +153,8 @@ Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos); if (Nmos == -1) { - fprintf (stderr, "cannot match mosaic for %s\n", image[i].name); + if (NmissMosaic < 1000) { + fprintf (stderr, "cannot match mosaic for %s\n", image[i].name); + } + NmissMosaic ++; continue; } @@ -164,4 +173,10 @@ } + fprintf (stderr, "mosaic matching : %d of possible %d failed to match\n", NmissMosaic, NtestMosaic); + if (NmissMosaic > 0.5*NtestMosaic) { + fprintf (stderr, "serious problem with mosaic matching\n"); + exit (5); + } + free (startMos); free (stopMos); @@ -250,5 +265,5 @@ if (im >= Nmosaic_for_images) abort(); - // search for the mosaic that + // search for the mosaic that matches this image mos = mosaic_for_images[im]; if (mos < 0) return NULL; @@ -256,2 +271,60 @@ return &mosaic[mos]; } + +// extend each host image table to include the mosaic 'images' needed by the host +int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage) { + + int i; + off_t j; + char *mosaicUsed; + + ALLOCATE (mosaicUsed, char, Nmosaic); + + // we need to add the mosaics to each of the region hosts lists of images + for (i = 0; i < regionHosts->Nhosts; i++) { + + int Nadd = 0; + int NADD = 100; + off_t *addMosaic = NULL; + ALLOCATE (addMosaic, off_t, NADD); + + // reset the mosaicUsed flags (valid only for this host) + memset (mosaicUsed, 0, Nmosaic * sizeof(char)); + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // find the mosaics associated with a given + for (j = 0; j < host->Nimage; j++) { + + int im = host->imseq[j]; + + if (im < 0) abort(); + if (im >= Nmosaic_for_images) abort(); + + // search for the mosaic that matches this image (skip unmatched images) + off_t mos = mosaic_for_images[im]; + if (mos < 0) continue; + + if (mosaicUsed[mos]) continue; + + mosaicUsed[mos] = TRUE; + addMosaic[Nadd] = mos; + Nadd ++; + + CHECK_REALLOCATE (addMosaic, off_t, NADD, Nadd, 100); + } + + REALLOCATE (host->image, Image, host->Nimage + Nadd); + + for (j = 0; j < Nadd; j++) { + off_t mos = addMosaic[j]; + off_t mos_im = mosaic[mos].myImage; + + host->image[host->Nimage + j] = image[mos_im]; + } + + host->Nimage += Nadd; + } + return TRUE; +} + Index: trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StarMaps.c (revision 36623) +++ trunk/Ohana/src/relastro/src/StarMaps.c (revision 36630) @@ -59,7 +59,6 @@ off_t i, N, Nimages; int xbin, ybin; - struct timeval start, stop; - gettimeofday (&start, (void *) NULL); + INITTIME; // Images *images = getimages(&Nimages, NULL); return value ignored Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 36623) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 36630) @@ -1,5 +1,5 @@ # include "relastro.h" -static off_t Nmax; +static off_t Nmax; static double *X, *dX; static double *Y, *dY; @@ -36,4 +36,22 @@ ALLOCATE (pX, double, MAX (1, Nmax)); ALLOCATE (pY, double, MAX (1, Nmax)); +} + +void freeObjectData () { + + free (R); + free (D); + free (T); + free (X); + free (Y); + + free (dR); + free (dD); + free (dT); + free (dX); + free (dY); + + free (pX); + free (pY); } @@ -205,5 +223,5 @@ if (((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) && (N <= PM_TOOFEW)) mode = FIT_AVERAGE; - if (FIT_TARGET == TARGET_HIGH_SPEED) { + if (RELASTRO_OP == OP_HIGH_SPEED) { Tmean = 0.5*(Tmax - Tmin); } else { @@ -412,4 +430,6 @@ } + freeObjectData (); + if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (NaveSum + NpmSum + NparSum), NaveSum, NpmSum, NparSum, NskipSum, NoffSum); return (TRUE); Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 36623) +++ trunk/Ohana/src/relastro/src/args.c (revision 36630) @@ -13,4 +13,5 @@ /* possible operations */ FIT_TARGET = TARGET_NONE; + RELASTRO_OP = OP_NONE; FIT_MODE = FIT_AVERAGE; @@ -21,5 +22,5 @@ if (N > argc - 6) usage_merge_source(); if (strcmp(argv[N+3], "into")) usage_merge_source(); - FIT_TARGET = TARGET_MERGE_SOURCE; + RELASTRO_OP = OP_MERGE_SOURCE; remove_argument (N, &argc, argv); OBJ_ID_SRC = strtol(argv[N], &endptr, 0); @@ -43,10 +44,51 @@ if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OBJECTS; + RELASTRO_OP = OP_UPDATE_OBJECTS; } if ((N = get_argument (argc, argv, "-update-offsets"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OFFSETS; + RELASTRO_OP = OP_UPDATE_OFFSETS; + } + + // elements needed for parallel regions / parallel images + REGION_FILE = NULL; + if ((N = get_argument (argc, argv, "-region-hosts"))) { + remove_argument (N, &argc, argv); + REGION_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + REGION_HOST_ID = 0; + if ((N = get_argument (argc, argv, "-region-hostID"))) { + remove_argument (N, &argc, argv); + REGION_HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + IMAGE_TABLE = NULL; + if ((N = get_argument (argc, argv, "-parallel-images"))) { + remove_argument (N, &argc, argv); + RELASTRO_OP = OP_PARALLEL_IMAGES; + if (N >= argc) usage(); + IMAGE_TABLE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + if (!REGION_FILE) usage(); + } + + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + RELASTRO_OP = OP_IMAGES; + } + + PARALLEL_REGIONS_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-regions"))) { + remove_argument (N, &argc, argv); + RELASTRO_OP = OP_PARALLEL_REGIONS; + if (!REGION_FILE) usage(); + if ((N = get_argument (argc, argv, "-parallel-regions-manual"))) { + remove_argument (N, &argc, argv); + PARALLEL_REGIONS_MANUAL = TRUE; + } } @@ -90,5 +132,5 @@ // XXX include a parallax / no-parallax option if (N >= argc - 4) usage(); - FIT_TARGET = TARGET_HIGH_SPEED; + RELASTRO_OP = OP_HIGH_SPEED; remove_argument (N, &argc, argv); PHOTCODE_A_LIST = strcreate(argv[N]); @@ -104,5 +146,5 @@ if ((N = get_argument (argc, argv, "-hpm"))) { if (N >= argc - 2) usage(); - FIT_TARGET = TARGET_HPM; + RELASTRO_OP = OP_HPM; remove_argument (N, &argc, argv); RADIUS = atof(argv[N]); @@ -116,5 +158,5 @@ remove_argument (N, &argc, argv); PARALLEL_OUTPUT = TRUE; - if ((FIT_TARGET != TARGET_HIGH_SPEED) && (FIT_TARGET != TARGET_HPM)) { + if ((RELASTRO_OP != OP_HIGH_SPEED) && (RELASTRO_OP != OP_HPM)) { fprintf (stderr, "-parallel-output only valid for -high-speed or -hpm modes\n"); exit (1); @@ -143,5 +185,7 @@ } - if (FIT_TARGET == TARGET_NONE) usage(); + if (RELASTRO_OP == OP_NONE) usage(); + + if (((RELASTRO_OP == OP_IMAGES) || (RELASTRO_OP == OP_PARALLEL_REGIONS) || (RELASTRO_OP == OP_PARALLEL_IMAGES)) && (FIT_TARGET == TARGET_NONE)) usage(); /* specify portion of the sky : allow default of all sky? */ @@ -427,5 +471,5 @@ /* possible operations */ - FIT_TARGET = TARGET_NONE; + RELASTRO_OP = TARGET_NONE; FIT_MODE = FIT_AVERAGE; @@ -437,4 +481,9 @@ BCATALOG = NULL; + REGION_FILE = NULL; + REGION_HOST_ID = 0; + IMAGE_TABLE = NULL; + PARALLEL_REGIONS_MANUAL = FALSE; + HOST_ID = 0; if ((N = get_argument (argc, argv, "-hostID"))) { @@ -457,15 +506,15 @@ BCATALOG = strcreate(argv[N]); remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_LOAD_OBJECTS; + RELASTRO_OP = OP_LOAD_OBJECTS; } if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OBJECTS; + RELASTRO_OP = OP_UPDATE_OBJECTS; } if ((N = get_argument (argc, argv, "-update-offsets"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_UPDATE_OFFSETS; + RELASTRO_OP = OP_UPDATE_OFFSETS; } @@ -487,5 +536,5 @@ // XXX include a parallax / no-parallax option if (N >= argc - 5) usage_client(); - FIT_TARGET = TARGET_HIGH_SPEED; + RELASTRO_OP = OP_HIGH_SPEED; remove_argument (N, &argc, argv); PHOTCODE_A_LIST = strcreate(argv[N]); @@ -501,5 +550,5 @@ if ((N = get_argument (argc, argv, "-hpm"))) { if (N >= argc - 3) usage(); - FIT_TARGET = TARGET_HPM; + RELASTRO_OP = OP_HPM; remove_argument (N, &argc, argv); RADIUS = atof(argv[N]); @@ -539,5 +588,5 @@ } - if (FIT_TARGET == TARGET_NONE) usage_client(); + if (RELASTRO_OP == OP_NONE) usage_client(); /* specify portion of the sky : allow default of all sky? */ @@ -718,7 +767,13 @@ void usage () { - fprintf (stderr, "ERROR: USAGE: relastro -update-simple [options]\n"); - fprintf (stderr, " OR: relastro -update-chips [options]\n"); - fprintf (stderr, " OR: relastro -update-mosaic [options]\n"); + fprintf (stderr, "ERROR: USAGE: relastro -images -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -images -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -images -update-mosaic [options]\n"); + fprintf (stderr, " OR: relastro -parallel-regions -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -parallel-regions -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -parallel-regions -update-mosaic [options]\n"); + fprintf (stderr, " OR: relastro -parallel-images -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -parallel-images -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -parallel-images -update-mosaic [options]\n"); fprintf (stderr, " OR: relastro -update-objects [options]\n"); fprintf (stderr, " OR: relastro -high-speed [options]\n"); Index: trunk/Ohana/src/relastro/src/assign_images.c =================================================================== --- trunk/Ohana/src/relastro/src/assign_images.c (revision 36630) +++ trunk/Ohana/src/relastro/src/assign_images.c (revision 36630) @@ -0,0 +1,280 @@ +# include "relastro.h" + +// This function generates a subset of the images based on selections. Input db has already +// been loaded with the raw fits table data +int assign_images (FITS_DB *db, RegionHostTable *regionHosts) { + + off_t Nimage; + + INITTIME; + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + Image *image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + MARKTIME("convert image table to internal structure: %f sec\n", dtime); + + // *** NOTE : for the moment, regions must be in the range 0 - 360, -90 - +90 + + // generate the chip match here so we can define the mosaic centers (if needed) + BuildChipMatch (image, Nimage); + MARKTIME("build chip match for %d images: %f sec\n", (int) Nimage, dtime); + + initMosaics (image, Nimage); + MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); + + // register the image array with ImageOps.c for later getimageByID calls + initImages (image, NULL, Nimage); + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + // for each regionHost, select images which are contained by the region + // even faster would be to use a tree to get to the real regions... + select_images_hostregion (regionHosts, image, Nimage); + + // supply the mosaics to the image table for the regionHosts : we already have the image + // <-> mosaic relationship, we just need to select these mosaics (once per regionHost) + select_mosaics_hostregion (regionHosts, image, Nimage); + + return TRUE; +} + +# define D_NIMAGE 1000 + +// assign images to the region hosts; at the end, each host will have its list of images +int select_images_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage) { + + off_t i, j; + + // INITTIME; + + for (i = 0; i < regionHosts->Nhosts; i++) { + regionHosts->hosts[i].Nimage = 0; + regionHosts->hosts[i].NIMAGE = D_NIMAGE; + ALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); + ALLOCATE (regionHosts->hosts[i].imseq, off_t, regionHosts->hosts[i].NIMAGE); + } + + for (j = 0; j < Nimage; j++) { + + /* select images by photcode, or equiv photcode, if specified */ + if (NphotcodesKeep > 0) { + int found = FALSE; + // XXX this bit of code excludes DIS mosaics and should be fixed + for (i = 0; (i < NphotcodesKeep) && !found; i++) { + if (photcodesKeep[i][0].code == image[j].photcode) found = TRUE; + if (photcodesKeep[i][0].code == GetPhotcodeEquivCodebyCode(image[j].photcode)) found = TRUE; + } + if (!found) continue; + } + if (NphotcodesSkip > 0) { + int found = FALSE; + for (i = 0; (i < NphotcodesSkip) && !found; i++) { + if (photcodesSkip[i][0].code == image[j].photcode) found = TRUE; + if (photcodesSkip[i][0].code == GetPhotcodeEquivCodebyCode(image[j].photcode)) found = TRUE; + } + if (found) continue; + } + + /* exclude images by time */ + if (TimeSelect) { + if (image[j].tzero < TSTART) continue; + if (image[j].tzero > TSTOP) continue; + } + + // do not include DIS (PHU-level mosaics) in the output list + if (!strcmp(&image[j].coords.ctype[4], "-DIS")) continue; + + // match the image to its corresponding ASTROMETRIC mosaic (not the same as the Mosaic above) + if (!FindMosaicForImage (image, Nimage, j)) { + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + // Exclude images with crazy astrometry + // XXX NOTE : this is gpc1-specific + { + double dP1 = hypot(image[j].coords.pc1_1, image[j].coords.pc1_2); + double dP2 = hypot(image[j].coords.pc2_1, image[j].coords.pc2_2); + if (fabs(dP1 - 1.0) > 0.02) continue; + if (fabs(dP2 - 1.0) > 0.02) continue; + + double X00, Y00, X10, Y10, X01, Y01; + XY_to_LM (&X00, &Y00, 0.0, 0.0, &image[j].coords); + XY_to_LM (&X10, &Y10, image[j].NX, 0.0, &image[j].coords); + XY_to_LM (&X01, &Y01, 0.0, image[j].NY, &image[j].coords); + double dS0 = hypot ((X00 - X10), (Y00 - Y10)); + double dS1 = hypot ((X00 - X01), (Y00 - Y01)); + if (dS0 > 6000) continue; + if (dS1 > 6500) continue; + } + + // use a reference coordinate for each image to assign to hosts + // define image center - note the DIS images (mosaic phu) are special + double Xc, Yc; + double Rc, Dc; + Xc = 0.5*image[j].NX; + Yc = 0.5*image[j].NY; + + XY_to_RD (&Rc, &Dc, Xc, Yc, &image[j].coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + image[j].RAo = Rc; + image[j].DECo = Dc; + + i = find_host_for_coords (regionHosts, Rc, Dc); + + if (i == -1) continue; + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // image bounds are defined for a range centered on the image center thus, an image + // with center 0.5 will have chips bounds ranging from ~ -1.5 to +2.5 or so, while an + // image with center 359.5 will have chips bounds ranging from ~ 357.5 to 361.5 or so + double Rmin, Rmax, Dmin, Dmax; + calculate_image_bounds (&image[j], &Rmin, &Rmax, &Dmin, &Dmax, Rc); + + host->RminCat = MIN(Rmin, host->RminCat); + host->RmaxCat = MAX(Rmax, host->RmaxCat); + host->DminCat = MIN(Dmin, host->DminCat); + host->DmaxCat = MAX(Dmax, host->DmaxCat); + + // regionHosts needs to have the full outer boundary + // (so reload_catalogs covers the correct region) + regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); + regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); + regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); + regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); + + // this is a bit memory expensive : I am making a complete copy of the image table here + off_t Nsubset = regionHosts->hosts[i].Nimage; + regionHosts->hosts[i].image[Nsubset] = image[j]; + regionHosts->hosts[i].imseq[Nsubset] = j; + + regionHosts->hosts[i].Nimage ++; + if (regionHosts->hosts[i].Nimage == regionHosts->hosts[i].NIMAGE) { + regionHosts->hosts[i].NIMAGE += D_NIMAGE; + REALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); + REALLOCATE (regionHosts->hosts[i].imseq, off_t, regionHosts->hosts[i].NIMAGE); + } + } + + return TRUE; +} + +double Xf[] = {0.0, 1.0, 0.0, 1.0}; +double Yf[] = {0.0, 0.0, 1.0, 1.0}; + +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid) { + + int n; + + double Rmin = 360.0; + double Rmax = 0.0; + double Dmin = +90.0; + double Dmax = -90.0; + + // define image corners + for (n = 0; n < 4; n++) { + double Xc, Yc, Rc, Dc; + Xc = Xf[n]*image->NX; + Yc = Yf[n]*image->NY; + XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, Rmid); + + Rmin = MIN (Rmin, Rc); + Rmax = MAX (Rmax, Rc); + Dmin = MIN (Dmin, Dc); + Dmax = MAX (Dmax, Dc); + } + + *rmin = Rmin; + *rmax = Rmax; + *dmin = Dmin; + *dmax = Dmax; + + return TRUE; +} + +int calculate_host_image_bounds (RegionHostTable *regionHosts) { + + int i, n; + off_t j; + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // XXX clear the chip match? + + BuildChipMatch (host->image, host->Nimage); + + double Rmin = 360.0; + double Rmax = 0.0; + double Dmin = +90.0; + double Dmax = -90.0; + + for (j = 0; j < host->Nimage; j++) { + + Image *image = &host->image[j]; + + if (!FindMosaicForImage (host->image, host->Nimage, j)) { + fprintf (stderr, "missing astrometry? programming error?\n"); + abort (); + } + + // define image corners + for (n = 0; n < 4; n++) { + double Xc, Yc, Rc, Dc; + Xc = Xf[n]*image->NX; + Yc = Yf[n]*image->NY; + XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + + Rmin = MIN (Rmin, Rc); + Dmin = MIN (Dmin, Rc); + + Rmax = MAX (Rmax, Rc); + Dmax = MAX (Dmax, Rc); + } + + } + + host->RminCat = Rmin; + host->DminCat = Dmin; + host->RmaxCat = Rmax; + host->DmaxCat = Dmax; + + // regionHosts needs to have the full outer boundary + // (so reload_catalogs covers the correct region) + regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); + regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); + regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); + regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); + } + return TRUE; +} + + +// XXX add a search tree to speed this up? +int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc) { + + int i; + + if (isnan(Rc)) return -1; + if (isnan(Dc)) return -1; + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + if (Rc < host->Rmin) continue; + if (Rc >= host->Rmax) continue; + if (Dc < host->Dmin) continue; + if (Dc >= host->Dmax) continue; + + return i; + } + return -1; +} Index: trunk/Ohana/src/relastro/src/indexCatalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/indexCatalogs.c (revision 36630) +++ trunk/Ohana/src/relastro/src/indexCatalogs.c (revision 36630) @@ -0,0 +1,66 @@ +# include "relastro.h" + +static int catIDmax = 0; +static int *catIDseq = NULL; +static int *objIDmax = NULL; +static int **objIDseq = NULL; + +int indexCatalogs (Catalog *catalog, int Ncatalog) { + + int i; + off_t j; + + if (!Ncatalog) return TRUE; + + // find the max value of catID + for (i = 0; i < Ncatalog; i++) { + catIDmax = MAX (catalog[i].catID, catIDmax); + } + + ALLOCATE (catIDseq, int, catIDmax + 1); + for (i = 0; i < catIDmax + 1; i++) { + catIDseq[i] = -1; + } + + for (i = 0; i < Ncatalog; i++) { + int catID = catalog[i].catID; + catIDseq[catID] = i; + } + + ALLOCATE (objIDmax, int, Ncatalog); + ALLOCATE (objIDseq, int *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + objIDmax[i] = 0; + for (j = 0; j < catalog[i].Naverage; j++) { + objIDmax[i] = MAX (catalog[i].average[j].objID, objIDmax[i]); + } + + ALLOCATE (objIDseq[i], int, objIDmax[i] + 1); + for (j = 0; j < objIDmax[i] + 1; j++) { + objIDseq[i][j] = -1; + } + + for (j = 0; j < catalog[i].Naverage; j++) { + int objID = catalog[i].average[j].objID; + objIDseq[i][objID] = j; + } + } + return TRUE; +} + +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq) { + + if (catID > catIDmax) return FALSE; + + int cat = catIDseq[catID]; + if (cat < 1) return FALSE; + + if (objID > objIDmax[cat]) return FALSE; + + int obj = objIDseq[cat][objID]; + if (obj < 1) return FALSE; + + *catSeq = cat; + *objSeq = obj; + return TRUE; +} Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 36623) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 36630) @@ -6,5 +6,5 @@ args (argc, argv); - if (FIT_TARGET == TARGET_MERGE_SOURCE) return; + if (RELASTRO_OP == OP_MERGE_SOURCE) return; fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 36630) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 36630) @@ -0,0 +1,144 @@ +# include "relastro.h" +# define DEBUG 0 + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + +int launch_region_hosts (RegionHostTable *regionHosts) { + + int i; + + // clear the I/O files + for (i = 0; i < regionHosts->Nhosts; i++) { + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync"); + truncate (syncfile, 0); + free (syncfile); + + char *fitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits"); + truncate (fitsfile, 0); + free (fitsfile); + + char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync"); + truncate (imsyncfile, 0); + free (imsyncfile); + + char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits"); + truncate (imfitsfile, 0); + free (imfitsfile); + + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + truncate (loopsyncfile, 0); + free (loopsyncfile); + } + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // communication files: + // subset images per host : CATDIR/Image.HOSTNAME.fits + char filename[1024]; + snprintf (filename, 1024, "%s/Image.%d.fits", CATDIR, host->hostID); + + // write the image subset for this host + ImageTableSave (filename, host->image, host->Nimage); + + char command[1024]; + snprintf (command, 1024, "relastro -parallel-images %s", filename); + strextend (command, "-region-hosts %s", REGION_FILE); + strextend (command, "-region-hostID %d", host->hostID); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-region %f %f %f %f", host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat); + strextend (command, "-statmode %s", STATMODE); + strextend (command, "-minerror %f", MIN_ERROR); + strextend (command, "-nloop %d", NLOOP); + strextend (command, "-threads %d", NTHREADS); + + switch (FIT_TARGET) { + case TARGET_SIMPLE: + strextend (command, "-update-simple"); + break; + case TARGET_CHIPS: + strextend (command, "-update-chips"); + break; + case TARGET_MOSAICS: + strextend (command, "-update-mosaics"); + break; + case TARGET_NONE: + abort(); + } + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (UPDATE) strextend (command, "-update"); + if (PARALLEL) strextend (command, "-parallel"); + if (PARALLEL_MANUAL) strextend (command, "-parallel-manual"); + if (PARALLEL_SERIAL) strextend (command, "-parallel-serial"); + if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (PhotFlagSelect) strextend (command, "+photflags"); + if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); + + if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); + + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + strextend (command, "-time %s %s", tstart, tstop); + free (tstart); + free (tstop); + } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_REGIONS_MANUAL) continue; + + // launch the job, then wait for it to be done loading catalogs. force the remote + // client to generate the file + char *syncfile = make_filename (CATDIR, host->hostname, host->hostID, "loadcat.sync"); + clear_sync_file (syncfile); + + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", host->hostname, command, host->stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", host->hostname, errorInfo); + exit (1); + } + host->pid = pid; // save for future reference + + // remove client is done, go ahead with next client + check_sync_file (syncfile, 1); + free (syncfile); + } + + if (PARALLEL_REGIONS_MANUAL) { + fprintf (stderr, "run the relastro_client commands above. when these are done, hit return\n"); + getchar(); + } else { + RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE); + } + + int status = TRUE; + for (i = 0; i < regionHosts->Nhosts; i++) { + status = status && (regionHosts->hosts[i].status == 0); + } + + return status; +} Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 36623) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 36630) @@ -115,4 +115,9 @@ Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + // load the list of hosts HostTable *table = HostTableLoad (CATDIR, sky->hosts); @@ -131,5 +136,5 @@ ALLOCATE (table->hosts[i].results, char, 1024); - snprintf (table->hosts[i].results, 1024, "%s/relastro.catalog.subset.dat", table->hosts[i].pathname); + snprintf (table->hosts[i].results, 1024, "%s/relastro.catalog.%s.dat", table->hosts[i].pathname, uniquer); // options / arguments that can affect relastro_client -load: Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 36623) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 36630) @@ -1,9 +1,3 @@ # include "relastro.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) { @@ -12,7 +6,6 @@ off_t Nimage, Nsubset; off_t *LineNumber; - struct timeval start, stop; - gettimeofday (&start, (void *) NULL); + INITTIME; // convert database table to internal structure @@ -35,5 +28,5 @@ /* unlock, if we can (else, unlocked below) */ - int unlockImages = !UPDATE || (FIT_TARGET == TARGET_UPDATE_OFFSETS); + int unlockImages = !UPDATE || (RELASTRO_OP == OP_UPDATE_OFFSETS); if (unlockImages) dvo_image_unlock (db); Index: trunk/Ohana/src/relastro/src/markObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/markObjects.c (revision 36630) +++ trunk/Ohana/src/relastro/src/markObjects.c (revision 36630) @@ -0,0 +1,24 @@ +# include "relastro.h" + +int markObjects (Catalog *catalog, int Ncatalog) { + + int i, n; + off_t j; + + // How strongly do I own this object? + for (i = 0; i < Ncatalog; i++) { + ALLOCATE (catalog[i].nOwn, int, catalog[i].Naverage); + memset (catalog[i].nOwn, 0, catalog[i].Naverage*sizeof(int)); + for (j = 0; j < catalog[i].Naverage; j++) { + int nOwn = 0; + int m = catalog[i].average[j].measureOffset; + for (n = 0; n < catalog[i].average[j].Nmeasure; n++) { + if (!catalog[i].measureT[m+n].myDet) continue; + nOwn ++; + } + catalog[i].nOwn[j] = nOwn; + } + } + return TRUE; +} + Index: trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro.c (revision 36623) +++ trunk/Ohana/src/relastro/src/relastro.c (revision 36630) @@ -10,33 +10,41 @@ SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - switch (FIT_TARGET) { - case TARGET_UPDATE_OBJECTS: + switch (RELASTRO_OP) { + case OP_UPDATE_OBJECTS: /* the object analysis is a separate process iterating over catalogs */ relastro_objects (skylist, 0, NULL); exit (0); - case TARGET_HIGH_SPEED: + case OP_HIGH_SPEED: /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */ high_speed_catalogs (sky, skylist, 0, NULL); exit (0); - case TARGET_HPM: + case OP_HPM: hpm_catalogs (sky, skylist, 0, NULL); exit (0); - case TARGET_MERGE_SOURCE: + case OP_MERGE_SOURCE: /* a special method to manually merge unlinked detections of sources togther (not parallel) */ relastro_merge_source (sky); exit (0); - case TARGET_SIMPLE: - case TARGET_CHIPS: - case TARGET_MOSAICS: + case OP_IMAGES: relastro_images (skylist); exit (0); - case TARGET_UPDATE_OFFSETS: + case OP_UPDATE_OFFSETS: // iterate over catalogs to make detection coordinates consistant UpdateObjectOffsets (skylist, 0, NULL); + exit (0); + + case OP_PARALLEL_REGIONS: + // run image updates in parallel across multiple remote machines + relastro_parallel_regions (); + exit (0); + + case OP_PARALLEL_IMAGES: + // operation on the remote machines in the PARALLEL_REGION mode + relastro_parallel_images (); exit (0); Index: trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_client.c (revision 36623) +++ trunk/Ohana/src/relastro/src/relastro_client.c (revision 36630) @@ -28,23 +28,23 @@ SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - switch (FIT_TARGET) { + switch (RELASTRO_OP) { - case TARGET_LOAD_OBJECTS: { + case OP_LOAD_OBJECTS: { // USAGE: relastro_client -load-objects int Ncatalog; Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR); if (!catalog) { - fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); - exit (2); + fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); + exit (2); } BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog); if (!BrightCatalogSave (BCATALOG, bcatalog)) { - fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); - exit (2); + fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); + exit (2); } break; } - case TARGET_UPDATE_OBJECTS: { + case OP_UPDATE_OBJECTS: { // USAGE: relastro_client -update-objects relastro_objects (skylist, HOST_ID, HOSTDIR); @@ -52,5 +52,5 @@ } - case TARGET_HIGH_SPEED: { + case OP_HIGH_SPEED: { // USAGE: relastro_client -high-speed high_speed_catalogs (sky, skylist, HOST_ID, HOSTDIR); @@ -58,5 +58,5 @@ } - case TARGET_HPM: { + case OP_HPM: { // USAGE: relastro_client -high-speed hpm_catalogs (sky, skylist, HOST_ID, HOSTDIR); @@ -64,6 +64,6 @@ } - // XXX loading the images is fairly costly -- see if we can do an image subset - case TARGET_UPDATE_OFFSETS: { + // XXX loading the images is fairly costly -- see if we can do an image subset? + case OP_UPDATE_OFFSETS: { FITS_DB db; Index: trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_images.c (revision 36623) +++ trunk/Ohana/src/relastro/src/relastro_images.c (revision 36630) @@ -6,7 +6,6 @@ Catalog *catalog; FITS_DB db; - struct timeval start, stop; - gettimeofday (&start, (void *) NULL); + INITTIME; /* register database handle with shutdown procedure */ Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 36623) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 36630) @@ -144,5 +144,4 @@ // PM_TOOFEW // SRC_MEAS_TOOFEW - // FIT_TARGET char command[1024]; Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 36630) +++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 36630) @@ -0,0 +1,131 @@ +# include "relastro.h" + +/* This function is essentially identical to relastro_images, except: + + * load the subset images saved by the master node + * distinguish detections we own (touch our images) and those we don't + * distinguish objects we own (in region) and those we don't + * update the unowned detections for owned objects from neighbor regions + * update the unowned objects for owned detections + */ + +int relastro_parallel_images () { + + int i, Ncatalog; + Catalog *catalog = NULL; + + INITTIME; + + // load the RegionTable (UserRegion should not be used at this level) + RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); + int myHost = regionHosts->index[REGION_HOST_ID]; + RegionHostFindNeighbors (regionHosts, myHost); + + // load the subset images belonging to this host + off_t Nimage; + Image *image = ImageTableLoad (IMAGE_TABLE, &Nimage); + if (!image) { + fprintf (stderr, "ERROR loading image %s\n", IMAGE_TABLE); + exit (2); + } + + // once we have read this table, we should remove it for repeat runs + // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection? + + // XXX need to deal with mosaic vs image... + initMosaics (image, Nimage); + + initImages (image, NULL, Nimage); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByBounds (sky, -1, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + initCoords(); + + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); + + // NOTE: if I let all hosts load blindly, I saturate the data clients with too many + // relastro_client requests. I need to have the master mediate this. the master + // will not launch the next remote job until this one says it is done + char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); + update_sync_file (syncfile, 1); + + // generate tables go from catID,objID -> catSeq,objSeq + indexCatalogs (catalog, Ncatalog); + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + // set test points based on the starmap + createStarMap (catalog, Ncatalog); + + markObjects (catalog, Ncatalog); + + SAVEPLOT = FALSE; + + /* major modes */ + switch (FIT_TARGET) { + case TARGET_SIMPLE: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + share_mean_pos (catalog, Ncatalog, regionHosts, i); + slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateSimple (catalog, Ncatalog); + share_image_pos (regionHosts, i); + slurp_image_pos (regionHosts, i); + } + break; + + case TARGET_CHIPS: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + share_mean_pos (catalog, Ncatalog, regionHosts, i); + slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateChips (catalog, Ncatalog); + share_image_pos (regionHosts, i); + slurp_image_pos (regionHosts, i); + MARKTIME("update chips: %f sec\n", dtime); + } + // create summary plots of the process + // relastroVisualSummaryChips(); + break; + + case TARGET_MOSAICS: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + share_mean_pos (catalog, Ncatalog, regionHosts, i); + slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateMosaic (catalog, Ncatalog); + share_image_pos (regionHosts, i); + slurp_image_pos (regionHosts, i); + } + break; + + default: + fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__); + exit (2); + } + + // this is a checkpoint to make sure all hosts have finished the loop above + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loop.sync"); + update_sync_file (loopsyncfile, 0); + free (loopsyncfile); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + char *loopsync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + check_sync_file (loopsync, 0); + free (loopsync); + } + + share_image_pos (regionHosts, -1); + + exit (0); +} Index: trunk/Ohana/src/relastro/src/relastro_parallel_regions.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 36630) +++ trunk/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 36630) @@ -0,0 +1,79 @@ +# include "relastro.h" + +int relastro_parallel_regions () { + + int status; + FITS_DB db; + + INITTIME; + + // load the RegionTable (UserRegion should not be used at this level) + RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); + if (!regionHosts) { + fprintf (stderr, "ERROR: problem with region host table\n"); + exit (2); + } + + // register database handle with shutdown procedure + set_db (&db); + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + // lock and load the image db table + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status && UPDATE) { + fprintf (stderr, "error\n"); + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + // if the Image.dat file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("-- load image data: %f sec\n", dtime); + + /* assign the images to the different region hosts */ + if (!assign_images (&db, regionHosts)) Shutdown ("error assigning images to region hosts"); + MARKTIME("-- assign images: %f sec\n", dtime); + + /* launch processing on the parallel region hosts */ + if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts"); + + // retrieve updated image parameters from the remote hosts (also set Image.mcal) + if (!slurp_image_pos (regionHosts, -1)) Shutdown ("error loading image updates"); + + if (!UPDATE) { + dvo_image_unlock (&db); + MARKTIME ("finished relastro -parallel-regions: %f sec total\n", dtime); + fprintf (stderr, "NOTE: UPDATE is OFF (results are not saved)\n"); + exit (0); + } + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByBounds (sky, -1, regionHosts->Rmin, regionHosts->Rmax, regionHosts->Dmin, regionHosts->Dmax); + + if (PARALLEL) { + // save the updated image parameters + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + } + + // iterate over catalogs to make detection coordinates consistant + UpdateObjectOffsets (skylist, 0, NULL); + + if (!PARALLEL) { + // save the changes to the image parameters + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + } + + MARKTIME ("finished relastro -parallel-regions: %f sec total\n", dtime); + + exit (0); +} + Index: trunk/Ohana/src/relastro/src/resort_catalog.c =================================================================== --- trunk/Ohana/src/relastro/src/resort_catalog.c (revision 36623) +++ trunk/Ohana/src/relastro/src/resort_catalog.c (revision 36630) @@ -30,6 +30,5 @@ if (catalog[0].sorted == TRUE) return; - // struct timeval start, stop; - // gettimeofday (&start, NULL); + // INITTIME; /* internal counters */ Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 36623) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 36630) @@ -15,10 +15,4 @@ void dsortindex (double *X, off_t *Y, int N); off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap) { @@ -31,5 +25,4 @@ Coords tcoords; SkyRegionCoords *skycoords; - struct timeval start, stop; double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; @@ -50,5 +43,5 @@ } - gettimeofday (&start, (void *) NULL); + INITTIME; // the comparison is made in the catalog local projection. below we set crval1,2 @@ -132,5 +125,6 @@ if (NphotcodesKeep > 0) { found = FALSE; - // XXX this bit of code excludes DIS mosaics and should be fixed + // we have to keep DIS mosaics explicitly (photcode = 0) + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) found = TRUE; for (k = 0; (k < NphotcodesKeep) && !found; k++) { if (photcodesKeep[k][0].code == timage[i].photcode) found = TRUE; Index: trunk/Ohana/src/relastro/src/share_images_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_images_pos.c (revision 36630) +++ trunk/Ohana/src/relastro/src/share_images_pos.c (revision 36630) @@ -0,0 +1,117 @@ +# include "relastro.h" + +// we are sharing image astrometry calibrations for all images which (a) I own and (b) which have unowned detections + +# define D_NIMAGEPOS 1000 +int share_image_pos (RegionHostTable *regionHosts, int nloop) { + + off_t i, Nimages; + Image *images = getimages (&Nimages, NULL); + + off_t Nimage_pos = 0; + off_t NIMAGE_POS = D_NIMAGEPOS; + + ImagePos *image_pos = NULL; + ALLOCATE (image_pos, ImagePos, NIMAGE_POS); + + for (i = 0; i < Nimages; i++) { + // XXX does this image have missing detections (does someone else need it?) + // XXX : NOTE NEED TO FIX THIS: if (imageExtra[i].Nmiss == 0) continue; + + set_image_pos (&image_pos[Nimage_pos], &images[i]); + Nimage_pos ++; + + CHECK_REALLOCATE (image_pos, ImagePos, NIMAGE_POS, Nimage_pos, D_NIMAGEPOS); + } + + // write out the image_mag fits table AND write state in some file + int myHost = regionHosts->index[REGION_HOST_ID]; + char *hostname = regionHosts->hosts[myHost].hostname; + + char *iposfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagepos.fits"); + ImagePosSave (iposfile, image_pos, Nimage_pos); + free (image_pos); + free (iposfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagepos.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_image_pos (RegionHostTable *regionHosts, int nloop) { + + off_t Nimage, i; + Image *images = getimages (&Nimage, NULL); + + int Nimage_pos = 0; + ImagePos *image_pos = NULL; + ALLOCATE (image_pos, ImagePos, 1); + + fprintf (stderr, "grabbing image mags from other hosts...\n"); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset; + char *iposfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.fits"); + ImagePos *image_pos_subset = ImagePosLoad (iposfile, &Nsubset); + free (iposfile); + + image_pos = merge_image_pos (image_pos, &Nimage_pos, image_pos_subset, Nsubset); + } + + for (i = 0; i < Nimage_pos; i++) { + off_t seq = getImageByID (image_pos[i].imageID); + if (seq < 0) { + // XXX is this a problem? (no, other hosts don't know which images I own) + continue; + } + images[seq].coords = image_pos[i].coords ; + images[seq].dXpixSys = image_pos[i].dXpixSys ; + images[seq].dYpixSys = image_pos[i].dYpixSys ; + images[seq].imageID = image_pos[i].imageID ; + images[seq].nFitAstrom = image_pos[i].nFitAstrom; + images[seq].flags = image_pos[i].flags ; + } + free (image_pos); + + fprintf (stderr, "DONE grabbing image mags from other hosts\n"); + + return TRUE; +} + +int set_image_pos (ImagePos *image_pos, Image *image) { + + image_pos->coords = image->coords; + image_pos->dXpixSys = image->dXpixSys; + image_pos->dYpixSys = image->dYpixSys; + image_pos->imageID = image->imageID; + image_pos->nFitAstrom = image->nFitAstrom; + image_pos->flags = image->flags; + + return TRUE; +} + +ImagePos *merge_image_pos (ImagePos *target, int *ntarget, ImagePos *source, int Nsource) { + + off_t i; + + REALLOCATE (target, ImagePos, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relastro/src/share_mean_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_mean_pos.c (revision 36630) +++ trunk/Ohana/src/relastro/src/share_mean_pos.c (revision 36630) @@ -0,0 +1,132 @@ +# include "relastro.h" +// we are sharing mean positions for all objects which (a) I own and (b) which have unowned detections + +# define D_NMEANPOS 10000 +int share_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + int i; + off_t j; + + off_t Nmeanpos = 0; + off_t NMEANPOS = D_NMEANPOS; + + MeanPos *meanpos = NULL; + ALLOCATE (meanpos, MeanPos, NMEANPOS); + + int myHost = regionHosts->index[REGION_HOST_ID]; + double Rmin = regionHosts->hosts[myHost].Rmin; + double Rmax = regionHosts->hosts[myHost].Rmax; + double Dmin = regionHosts->hosts[myHost].Dmin; + double Dmax = regionHosts->hosts[myHost].Dmax; + + // XXX skip some catalogs based on UserPatch? + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + // do I own this object? (in region range?) --- CAREFUL HERE!! + if (catalog[i].average[j].R < Rmin) continue; + if (catalog[i].average[j].R >= Rmax) continue; + if (catalog[i].average[j].D < Dmin) continue; + if (catalog[i].average[j].D >= Dmax) continue; + + // does this object have missing detections (does someone else need it?) + // XXX : sky objects without missing detections + // XXX watch out for detections which are not associated with an image (REF) + if (catalog[i].nOwn[j] == catalog[i].average[j].Nmeasure) continue; + + set_mean_pos (&meanpos[Nmeanpos], &catalog[i].average[j]); + Nmeanpos ++; + CHECK_REALLOCATE (meanpos, MeanPos, NMEANPOS, Nmeanpos, D_NMEANPOS); + } + } + + // write out the meanmag fits table AND write state in some file + char *hostname = regionHosts->hosts[myHost].hostname; + + char *posfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanpos.fits"); + MeanPosSave (posfile, meanpos, Nmeanpos); + free (meanpos); + free (posfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanpos.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_mean_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + off_t i; + + int Nmeanpos = 0; + MeanPos *meanpos = NULL; + ALLOCATE (meanpos, MeanPos, 1); + + fprintf (stderr, "grabbing mean object pos from other hosts...\n"); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset = 0; + char *posfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.fits"); + MeanPos *meanposSubset = MeanPosLoad (posfile, &Nsubset); + free (posfile); + + // merge_mean_pos reallocs meanpos and frees the input meanposSubset + meanpos = merge_mean_pos (meanpos, &Nmeanpos, meanposSubset, Nsubset); + } + + for (i = 0; i < Nmeanpos; i++) { + int objID = meanpos[i].objID; + int catID = meanpos[i].catID; + + // set the mean mag + int catSeq; + off_t objSeq; + if (!catID_and_objID_to_seq (catID, objID, &catSeq, &objSeq)) { + // XXX what should I do if this does not match? + continue; + } + + catalog[catSeq].average[objSeq].R = meanpos[i].R; + catalog[catSeq].average[objSeq].D = meanpos[i].D; + } + free (meanpos); + + fprintf (stderr, "DONE grabbing mean object pos from other hosts...\n"); + + return TRUE; +} + +int set_mean_pos (MeanPos *meanpos, Average *average) { + + meanpos->R = average->R; + meanpos->D = average->D; + meanpos->objID = average->objID; + meanpos->catID = average->catID; + + return TRUE; +} + +MeanPos *merge_mean_pos (MeanPos *target, int *ntarget, MeanPos *source, int Nsource) { + + off_t i; + + REALLOCATE (target, MeanPos, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relastro/src/syncfile.c =================================================================== --- trunk/Ohana/src/relastro/src/syncfile.c (revision 36630) +++ trunk/Ohana/src/relastro/src/syncfile.c (revision 36630) @@ -0,0 +1,77 @@ +# include "relastro.h" +# define MSG_LENGTH 12 + +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname) { + + char tmp[10], *line; + int Nchar = snprintf (tmp, 0, "%s/%s.%03d.%s", dirname, hostname, hostID, tailname); + + ALLOCATE (line, char, Nchar + 1); + snprintf (line, Nchar + 1, "%s/%s.%03d.%s", dirname, hostname, hostID, tailname); + + return line; +} + +int check_sync_file (char *filename, int nloop) { + + char message[MSG_LENGTH]; + + FILE *f = NULL; + + while (TRUE) { + + f = fopen (filename, "r"); + if (!f) { + usleep (2000000); + continue; + } + + // XXX MSG_LENGTH : 0 EOL byte? + int Nread = fread (message, 1, MSG_LENGTH, f); + if (Nread < MSG_LENGTH) { + fclose (f); + usleep (2000000); + continue; + } + fclose (f); + + // message is of the form: NLOOP: %03d + int loop; + sscanf (message, "%*s %d", &loop); + if (loop != nloop) { + usleep (2000000); + continue; + } + return TRUE; + } + return FALSE; +} + +int clear_sync_file (char *filename) { + // delete file contents + truncate (filename, 0); + + return TRUE; +} + +int update_sync_file (char *filename, int nloop) { + + char message[MSG_LENGTH]; + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "failure to open sync file for write\n"); + exit (4); + } + + snprintf (message, MSG_LENGTH, "NLOOP: %03d\n", nloop); + + int Nwrite = fwrite (message, 1, MSG_LENGTH, f); + if (Nwrite != MSG_LENGTH) { + fprintf (stderr, "failure to write sync message\n"); + exit (3); + } + + fclose (f); + return TRUE; +} Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 36623) +++ trunk/Ohana/src/relphot/Makefile (revision 36630) @@ -47,5 +47,16 @@ $(SRC)/relphot_images.$(ARCH).o \ $(SRC)/relphot_objects.$(ARCH).o \ +$(SRC)/relphot_parallel_regions.$(ARCH).o \ +$(SRC)/relphot_parallel_images.$(ARCH).o \ $(SRC)/select_images.$(ARCH).o \ +$(SRC)/assign_images.$(ARCH).o \ +$(SRC)/launch_region_hosts.$(ARCH).o \ +$(SRC)/ImageTable.$(ARCH).o \ +$(SRC)/ImageMagIO.$(ARCH).o \ +$(SRC)/syncfile.$(ARCH).o \ +$(SRC)/share_image_mags.$(ARCH).o \ +$(SRC)/share_mean_mags.$(ARCH).o \ +$(SRC)/MeanMagIO.$(ARCH).o \ +$(SRC)/indexCatalog.$(ARCH).o \ $(SRC)/client_logger.$(ARCH).o \ $(SRC)/setExclusions.$(ARCH).o \ Index: trunk/Ohana/src/relphot/doc/parallel.txt =================================================================== --- trunk/Ohana/src/relphot/doc/parallel.txt (revision 36623) +++ trunk/Ohana/src/relphot/doc/parallel.txt (revision 36630) @@ -1,2 +1,71 @@ + +2014.02.15 + +Nearly done with the relphot mods. some outstanding questions; + + * image vs mosaic? + * how do decide which images need to be shared? + +2014.02.14 + + more relphot notes: + + * I load the catalogs, with a + +2014.02.12 + +I am making progess on the relphot -parallel-regions implementation. some things I need to deal with: + + * assign_images: select mosaics for all images saved for a given host (need to supply them as well) + ** I need to decide on the Image / Mosaic split. If I am calibrating by image, then the image center defines ownership + ** if I am calibrating by exposure, then the mosaic center must define ownership + o determine SkyList covering the images for region host + o tag detections and objects which I own or do not own + o ImageOps.c / matchImage matches detections to images and needs to handle mine / not mine cases + o ditto for matchMosaics + o load_catalogs : I need to make my requested catalogs unique (just add uniquer to load_catalogs.c:130 + o figure out what fields I'm saving in meanmags and how to construct it + -- mag, objID, catID, photcode, (dmag or other stats?) + * from which hosts do I slurp mean mags? + o function to merge a new meanmag array into the existing one + o how to go from (objID,catID) to a given catalog[i].average[j] + +2014.02.06 + +** a single region-level host owns images for which the center lands + in its region. + +** it owns detections which come from images which it owns + +** it owns objects which land in its region + +Extending parallel relphot processing to split the sky (and images) +into regions, each of which runs in parallel at the same time: + +Top Level (relphot -parallel-images): + * define regions of the sky -> hosts + * load images, assign to hosts + * launch region-level jobs on remote hosts + + Region Level (relphot -parallel-images-region) + * load my image subset table + * request objects and detections for my skyregion + * match images & objects, etc + + * update my image parameters + * update my detections + * write out detections + * load detections from my border hosts + * match to my objects + * update my objects + * write out objects + * load objects from my border hosts + * update images + (iterate N times) + * write out image parameters + +Top Level + * read image parameters, update + * update objects 2012.02.13 Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 36623) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 36630) @@ -23,4 +23,13 @@ // # define IDX_T off_t # define IDX_T int + +typedef enum { + MODE_ERROR = 0, + UPDATE_IMAGES, + UPDATE_AVERAGES, + PARALLEL_REGIONS, + PARALLEL_IMAGES, + APPLY_OFFSETS, +} RelphotMode; typedef enum { @@ -108,4 +117,25 @@ typedef struct { + float M; + float dM; + float Xm; + int Nsec; + unsigned int objID; + unsigned int catID; + int photcode; +} MeanMag; + +typedef struct { + float Mcal; + float dMcal; + float dMagSys; + float Xm; + int nFitPhotom; + int flags; + unsigned int imageID; + short ubercalDist; +} ImageMag; + +typedef struct { AverageTiny *average; // array of (minimal) average data MeasureTiny *measure; // array of (minimal) measure data @@ -155,4 +185,9 @@ int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ +// globals for parallel region operations +char *REGION_FILE; +char *IMAGE_TABLE; +int REGION_HOST_ID; + int HOST_ID; char *HOSTDIR; @@ -179,4 +214,6 @@ int PARALLEL_MANUAL; int PARALLEL_SERIAL; + +int PARALLEL_REGIONS_MANUAL; int NTHREADS; @@ -266,5 +303,5 @@ void InterpolateGrid PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords)); off_t *SelectRefMosaic PROTO((Mosaic **refmosaic, off_t *Nimage)); -int args PROTO((int argc, char **argv)); +RelphotMode args PROTO((int argc, char **argv)); int args_client PROTO((int argc, char **argv)); int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog, int Ncat)); @@ -279,4 +316,8 @@ void findImages PROTO((Catalog *catalog, int Ncatalog, int doImageList)); int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); + +void makeMosaics (Image *image, off_t Nimage); +Mosaic *getMosaicForImage (off_t im); +void setMosaicCenters (Image *image, off_t Nimage); void set_db (FITS_DB *in); @@ -315,5 +356,5 @@ void initMosaics PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage)); void initMrel PROTO((Catalog *catalog, int Ncatalog)); -void initialize PROTO((int argc, char **argv)); +RelphotMode initialize PROTO((int argc, char **argv)); void initialize_client PROTO((int argc, char **argv)); void liststats_setmode PROTO((StatType *stats, char *strmode)); @@ -426,2 +467,44 @@ int init_synthetic_mags (); int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm); + +int relphot_parallel_regions (); + +int assign_images (FITS_DB *db, RegionHostTable *regionHosts); +int select_images_hostregion (RegionHostTable *hosts, Image *image, off_t Nimage); +int find_host_for_coords (RegionHostTable *regionHosts, double R, double d); +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid); + +int launch_region_hosts (RegionHostTable *regionHosts); +int strextend (char *input, char *format,...); + +Image *ImageTableLoad(char *filename, off_t *nimage); +int ImageTableSave (char *filename, Image *images, off_t Nimages); + +int indexCatalogs (Catalog *catalog, int Ncatalog); +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq); + +int check_sync_file (char *filename, int nloop); +int clear_sync_file (char *filename); +int update_sync_file (char *filename, int nloop); +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname); + +int share_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int slurp_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int set_mean_mags (MeanMag *meanmags, AverageTiny *average, SecFilt *secfilt, int Nsec); +MeanMag *merge_mean_mags (MeanMag *target, int *ntarget, MeanMag *source, int Nsource); + +MeanMag *MeanMagLoad(char *filename, off_t *nmeanmags); +int MeanMagSave(char *filename, MeanMag *meanmags, off_t Nmeanmags); + +int share_image_mags (RegionHostTable *regionHosts, int nloop); +int slurp_image_mags (RegionHostTable *regionHosts, int nloop); +int set_image_mags (ImageMag *image_mags, Image *image); +ImageMag *merge_image_mags (ImageMag *target, int *ntarget, ImageMag *source, int Nsource); + +ImageMag *ImageMagLoad(char *filename, off_t *nimage_mags); +int ImageMagSave(char *filename, ImageMag *image_mags, off_t Nimage_mags); + +int markObjects (Catalog *catalog, int Ncatalog); + +int relphot_parallel_images (); +int relphot_parallel_regions (); Index: trunk/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 36623) +++ trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 36630) @@ -101,4 +101,5 @@ measure[i].catID = catID[i]; measure[i].photcode = photcode[i]; + measure[i].myDet = FALSE; } fprintf (stderr, "loaded data for %lld measure\n", (long long) Nrow); @@ -141,4 +142,5 @@ GET_COLUMN(flags, "FLAGS", int); GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(objID, "OBJ_ID", int); gfits_free_header (&theader); gfits_free_table (&ftable); @@ -153,4 +155,6 @@ average[i].flags = flags[i]; average[i].catID = catID[i]; + average[i].objID = objID[i]; + average[i].nOwn = 0; } fprintf (stderr, "loaded data for %lld average\n", (long long) Nrow); @@ -162,4 +166,5 @@ free (flags ); free (catID ); + free (objID ); catalog->average = average; @@ -370,4 +375,5 @@ gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ref", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ref", NULL, 1.0, 0.0); // generate the output array that carries the data @@ -381,4 +387,5 @@ int *flags ; ALLOCATE (flags, int, catalog->Naverage); int *catID ; ALLOCATE (catID, int, catalog->Naverage); + int *objID ; ALLOCATE (objID, int, catalog->Naverage); // assign the storage arrays @@ -391,4 +398,5 @@ flags[i] = average[i].flags; catID[i] = average[i].catID; + objID[i] = average[i].objID; } @@ -400,4 +408,5 @@ gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, catalog->Naverage); free (R ); @@ -407,4 +416,5 @@ free (flags ); free (catID ); + free (objID ); gfits_fwrite_Theader (f, &theader); Index: trunk/Ohana/src/relphot/src/ImageMagIO.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 36630) +++ trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 36630) @@ -0,0 +1,221 @@ +# include "relphot.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +ImageMag *ImageMagLoad(char *filename, off_t *nimage_mags) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage_mags = 0; + ImageMag *image_mags = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // we read the entire block of data, then extract the columns, then set the image structure values. + // I free the FITS table data after extracting the colums to avoid having 3 copies in memory. + + char type[16]; + + GET_COLUMN (Mcal, "MCAL", float); + GET_COLUMN (dMcal, "MCAL_ERR", float); + GET_COLUMN (dMagSys, "MCAL_SYSERR", float); + GET_COLUMN (Xm, "CHISQ", float); + GET_COLUMN (nFitPhotom, "NFIT", int); + GET_COLUMN (flags, "FLAGS", int); + GET_COLUMN (ubercalDist, "UDIST", int); + GET_COLUMN (imageID, "ID", int); + + // free the memory associated with the FITS files + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + ALLOCATE (image_mags, ImageMag, Nrow); + for (i = 0; i < Nrow; i++) { + image_mags[i].Mcal = Mcal [i]; + image_mags[i].dMcal = dMcal [i]; + image_mags[i].dMagSys = dMagSys [i]; + image_mags[i].Xm = Xm [i]; + image_mags[i].nFitPhotom = nFitPhotom [i]; + image_mags[i].flags = flags [i]; + image_mags[i].ubercalDist = ubercalDist[i]; + image_mags[i].imageID = imageID [i]; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (Mcal ); + free (dMcal ); + free (dMagSys ); + free (Xm ); + free (nFitPhotom ); + free (flags ); + free (ubercalDist); + free (imageID ); + + *nimage_mags = Nrow; + return image_mags; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int ImageMagSave(char *filename, ImageMag *image_mags, off_t Nimage_mags) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "IMAGE_MAGS"); + + gfits_define_bintable_column (&theader, "E", "MCAL", "cal offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "cal error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR", "systematic error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ", "cal chisq", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "UDIST", "distance to ubercal images", "images", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + float *Mcal ; + float *dMcal ; + float *dMagSys ; + float *Xm ; + int *nFitPhotom ; + int *flags ; + int *ubercalDist ; + int *imageID ; + + // create intermediate storage arrays + ALLOCATE (Mcal , float, Nimage_mags); + ALLOCATE (dMcal , float, Nimage_mags); + ALLOCATE (dMagSys , float, Nimage_mags); + ALLOCATE (Xm , float, Nimage_mags); + ALLOCATE (nFitPhotom , int, Nimage_mags); + ALLOCATE (flags , int, Nimage_mags); + ALLOCATE (ubercalDist , int, Nimage_mags); + ALLOCATE (imageID , int, Nimage_mags); + + // assign the storage arrays + for (i = 0; i < Nimage_mags; i++) { + Mcal [i] = image_mags[i].Mcal ; + dMcal [i] = image_mags[i].dMcal ; + dMagSys [i] = image_mags[i].dMagSys ; + Xm [i] = image_mags[i].Xm ; + nFitPhotom [i] = image_mags[i].nFitPhotom ; + flags [i] = image_mags[i].flags ; + ubercalDist[i] = image_mags[i].ubercalDist; + imageID [i] = image_mags[i].imageID ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MCAL", Mcal , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_ERR", dMcal , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "MCAL_SYSERR", dMagSys , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "CHISQ", Xm , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitPhotom , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "UDIST", ubercalDist, Nimage_mags); + gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_mags); + + free (Mcal ); + free (dMcal ); + free (dMagSys ); + free (Xm ); + free (nFitPhotom ); + free (flags ); + free (ubercalDist); + free (imageID ); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image_mags file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for image_mags %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for image_mags %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for image_mags %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for image_mags %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file image_mags %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file image_mags %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing image_mags file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 36623) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 36630) @@ -244,7 +244,11 @@ int Nmatch = 0; - for (i = 0; i < Ncatalog; i++) { + for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { + catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode); + + // skip measurements which do not match one of the requested photcodes ( + // (do we not already exclude in bcatalog -- maybe needed for reload_objects? found = FALSE; for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { @@ -252,4 +256,6 @@ } if (!found) continue; + + // if we match one of our images, myDet gets set to TRUE matchImage (catalog, j, i, doImageList); Nmatch ++; @@ -322,4 +328,5 @@ return; } + catalog[cat].measureT[meas].myDet = TRUE; if (USE_GRID) { Index: trunk/Ohana/src/relphot/src/ImageTable.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageTable.c (revision 36630) +++ trunk/Ohana/src/relphot/src/ImageTable.c (revision 36630) @@ -0,0 +1,62 @@ +# include "relphot.h" + +Image *ImageTableLoad(char *filename, off_t *nimage) { + + int status; + off_t Nimage; + FITS_DB db; + + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, filename, 60.0, LCK_SOFT); + if (!status) { + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + // if the file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + Shutdown ("ERROR: failed to read images"); + } + + *nimage = Nimage; + return image; +} + +int ImageTableSave (char *filename, Image *images, off_t Nimages) { + + int status; + FITS_DB db; + + // setup image table format and lock + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&db, filename, 60.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + // load or create the image table + if (db.dbstate != LCK_EMPTY) { + Shutdown ("image catalog already exists %s", db.filename); + } + + dvo_image_create (&db, GetZeroPoint()); + + /* add the new images and save */ + dvo_image_addrows (&db, images, Nimages); + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/MeanMagIO.c =================================================================== --- trunk/Ohana/src/relphot/src/MeanMagIO.c (revision 36630) +++ trunk/Ohana/src/relphot/src/MeanMagIO.c (revision 36630) @@ -0,0 +1,200 @@ +# include "relphot.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +MeanMag *MeanMagLoad(char *filename, off_t *nmeanmags) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nmeanmags = 0; + MeanMag *meanmags = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // a bit annoying : we read the entire block of data, then extract the columns, then set the image structure values. + // this means I need 3 copies in memory at some point. ugh. + + char type[16]; + + GET_COLUMN (M, "MAG", float); + GET_COLUMN (dM, "MAG_ERR", float); + GET_COLUMN (Xm, "MAG_CHISQ", float); + GET_COLUMN (Nsec, "NSEC", int); + GET_COLUMN (objID, "OBJ_ID", int); + GET_COLUMN (catID, "CAT_ID", int); + + // free the memory associated with the FITS files + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + ALLOCATE (meanmags, MeanMag, Nrow); + for (i = 0; i < Nrow; i++) { + meanmags[i].M = M [i]; + meanmags[i].dM = dM [i]; + meanmags[i].Xm = Xm [i]; + meanmags[i].Nsec = Nsec [i]; + meanmags[i].objID = objID[i]; + meanmags[i].catID = catID[i]; + } + fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); + + free (M ); + free (dM ); + free (Xm ); + free (Nsec ); + free (objID); + free (catID); + + *nmeanmags = Nrow; + return meanmags; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int MeanMagSave(char *filename, MeanMag *meanmags, off_t Nmeanmags) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "MEANMAGS"); + + gfits_define_bintable_column (&theader, "E", "MAG", "mean magnitude", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "mean magnitude error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "mean magnitude chisq", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NSEC", "secfilt sequence", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ID", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ID", NULL, 1.0, 1.0*0x8000); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + float *M, *dM, *Xm; + int *Nsec; + unsigned int *objID, *catID; + + // create intermediate storage arrays + ALLOCATE (M, float, Nmeanmags); + ALLOCATE (dM, float, Nmeanmags); + ALLOCATE (Xm, float, Nmeanmags); + ALLOCATE (Nsec, int, Nmeanmags); + ALLOCATE (objID, unsigned int, Nmeanmags); + ALLOCATE (catID, unsigned int, Nmeanmags); + + // assign the storage arrays + for (i = 0; i < Nmeanmags; i++) { + M [i] = meanmags[i].M ; + dM [i] = meanmags[i].dM ; + Xm [i] = meanmags[i].Xm ; + Nsec [i] = meanmags[i].Nsec ; + objID[i] = meanmags[i].objID; + catID[i] = meanmags[i].catID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MAG", M, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Xm, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "NSEC", Nsec, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeanmags); + + free (M ); + free (dM ); + free (Xm ); + free (Nsec ); + free (objID); + free (catID); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanmags %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanmags %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for meanmags %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for meanmags %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file meanmags %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file meanmags %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 36623) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 36630) @@ -320,4 +320,157 @@ } +/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes but we + match images by MOSAICNAME. this last point is weak: it forces a single camera at a time. + we can extend the logic to multiple cameras if we make list of MOSAICNAMES (better to assign a camera ID) + */ +void makeMosaics (Image *image, off_t Nimage) { + + off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; + unsigned int start, stop, *startTimes, *startTimesMosaic; + char *pname; + + if (!MOSAIC_ZEROPT) return; + + INITTIME; + + /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no + * entry in the image table that represents this mosaic. Instead, it is an + * internal construct that defines a group of related images + */ + + // generate a list of all image start times + ALLOCATE (startTimes, unsigned int, Nimage); + for (i = 0; i < Nimage; i++) { + startTimes[i] = image[i].tzero; + } + sort_times (startTimes, Nimage); + MARKTIME("create array of all image obstimes: %f sec\n", dtime); + + Nmosaic = 0; + NMOSAIC = 1000; + ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC); + startTimesMosaic[0] = startTimes[0]; + + // generate a list of the unique start times (these define the mosaics) + for (i = 0; i < Nimage; i++) { + if (startTimes[i] < startTimesMosaic[Nmosaic]) { + fprintf (stderr, "error?\n"); + abort(); + } + if (startTimes[i] == startTimesMosaic[Nmosaic]) continue; + Nmosaic ++; + if (Nmosaic >= NMOSAIC) { + NMOSAIC += 1000; + REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC); + } + startTimesMosaic[Nmosaic] = startTimes[i]; + } + Nmosaic ++; + MARKTIME("create array of mosaic obstimes: %f sec\n", dtime); + + // now I have a list of uniq start times, and they are in order + // create the mosaic arrays for these times + ALLOCATE (mosaic, Mosaic, Nmosaic); + + ALLOCATE (MosaicToImage, off_t *, Nmosaic); + ALLOCATE (MosaicN_Image, off_t, Nmosaic); + ALLOCATE (MosaicN_IMAGE, off_t, Nmosaic); + + // init the mosaic array values + for (i = 0; i < Nmosaic; i++) { + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].Xm = 0.0; + mosaic[i].flags = 0; + mosaic[i].secz = NAN; + mosaic[i].photcode = 0; + mosaic[i].skipCal = FALSE; + + memset (&mosaic[i].coords, 0, sizeof(Coords)); + + MosaicN_IMAGE[i] = 10; + MosaicN_Image[i] = 0; + ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]); + MosaicToImage[i][0] = -1; + } + + ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs + + // assign each image to a mosaic + for (i = 0; i < Nimage; i++) { + ImageToMosaic[i] = -1; + + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (image[i].photcode); + if (!pname) continue; + + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + + start = image[i].tzero; + stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + j = findMosaic(startTimesMosaic, Nmosaic, start); + if (j == -1) { + fprintf (stderr, "programming error? all image images should belong to a mosaic\n"); + abort(); + } + + // add reference from image to mosaic + ImageToMosaic[i] = j; + + // have we already found this mosaic? + found = (MosaicN_Image[j] > 0); + + /* add image to mosaic image list */ + MosaicToImage[j][MosaicN_Image[j]] = i; + MosaicN_Image[j] ++; + if (MosaicN_Image[j] == MosaicN_IMAGE[j]) { + MosaicN_IMAGE[j] += 10; + REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]); + } + if (found) continue; + + /* a new mosaic, define ranges */ + if (mosaic[j].start != start) { + fprintf (stderr, "error?\n"); + abort(); + } + mosaic[j].stop = stop; + mosaic[j].Mcal = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].Xm = 0.0; + mosaic[j].dMsys = image[i].flags; + mosaic[j].flags = image[i].flags; + mosaic[j].secz = image[i].secz; + mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + } + MARKTIME("assign images to mosaic: %f sec\n", dtime); + + // free this or not? + free (MosaicN_IMAGE); + free (startTimes); + free (startTimesMosaic); + + fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); + + return; +} + +Mosaic *getMosaicForImage (off_t im) { + + if (im < 0) return NULL; + if (!ImageToMosaic) return NULL; + + off_t m = ImageToMosaic[im]; + if (m < 0) return NULL; + if (m >= Nmosaic) return NULL; + + return (&mosaic[m]); +} + // use bisection to find the overlapping mosaic (returns exact match) // startTimes is a sorted, unique list of times @@ -350,4 +503,123 @@ } return (-1); +} + +void setMosaicCenters (Image *image, off_t Nimage) { + + /* find max dR, dD range for all mosaics */ + /* define mosaic.coords to cover dR, dD */ + /* send results to initGridBins */ + + off_t i, j, m, NX, NY, NC, Nc; + double R, D, Rmid, Dmid; + double Mcal, dMcal, Xm; + double *Rc, *Dc; + + NC = 100; + ALLOCATE (Rc, double, NC); + ALLOCATE (Dc, double, NC); + + for (i = 0; i < Nmosaic; i++) { + Nc = 0; + Rmid = Dmid = NAN; + Mcal = dMcal = Xm = 0; + for (j = 0; j < MosaicN_Image[i]; j++) { + m = MosaicToImage[i][j]; + + if (!FindMosaicForImage (image, Nimage, m)) { + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + NX = image[m].NX; + NY = image[m].NY; + XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &image[m].coords); + R = ohana_normalize_angle_to_midpoint (R, 180.0); + + // Exclude images with crazy astrometry + // XXX NOTE : this is gpc1-specific + { + double dP1 = hypot(image[m].coords.pc1_1, image[m].coords.pc1_2); + double dP2 = hypot(image[m].coords.pc2_1, image[m].coords.pc2_2); + if (fabs(dP1 - 1.0) > 0.02) continue; + if (fabs(dP2 - 1.0) > 0.02) continue; + + double X00, Y00, X10, Y10, X01, Y01; + XY_to_LM (&X00, &Y00, 0.0, 0.0, &image[m].coords); + XY_to_LM (&X10, &Y10, image[m].NX, 0.0, &image[m].coords); + XY_to_LM (&X01, &Y01, 0.0, image[m].NY, &image[m].coords); + double dS0 = hypot ((X00 - X10), (Y00 - Y10)); + double dS1 = hypot ((X00 - X01), (Y00 - Y01)); + if (dS0 > 6000) continue; + if (dS1 > 6500) continue; + } + + Rc[Nc] = R; + Dc[Nc] = D; + Nc ++; + if (Nc >= NC) { + NC += 100; + REALLOCATE (Rc, double, NC); + REALLOCATE (Dc, double, NC); + } + + Mcal += image[m].Mcal; + dMcal += image[m].dMcal; + Xm += image[m].Xm; + + // for ubercal images, we (elsewhere) keep Mcal frozen + + /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ + image[m].Mcal = 0.0; + image[m].dMcal = NAN; + image[m].Xm = NAN_S_SHORT; + } + + if (Nc > 0) { + dsort (Rc, Nc); + if (Rc[Nc-1] - Rc[0] > 180.0) { + // in our list, Rc is in the range 0.0 to 360.0. + // any mosaic which is close to the 0.0, 360.0 boundary may have some on + // one side or the other. count how many have values more than Rc[0] + 180. + // if more than half are at the large end, re-normalize to that range + int Nbig = 0; + for (j = 1; j < Nc; j++) { + if (Rc[j] - Rc[0] > 180.0) Nbig ++; + } + if (Nbig > 0.5*Nc) { + for (j = 0; j < Nc; j++) { + Rc[j] = ohana_normalize_angle_to_midpoint (Rc[j], 360.0); + } + dsort (Rc, Nc); + } else if (Nbig > 0) { + for (j = 0; j < Nc; j++) { + Rc[j] = ohana_normalize_angle_to_midpoint (Rc[j], 0.0); + } + dsort (Rc, Nc); + } + } + dsort (Dc, Nc); + + Rmid = Rc[(int)(0.5*Nc)]; + Dmid = Dc[(int)(0.5*Nc)]; + } + + strcpy (mosaic[i].coords.ctype, "DEC--TAN"); + mosaic[i].coords.crval1 = Rmid; + mosaic[i].coords.crval2 = Dmid; + mosaic[i].coords.crpix1 = 0.0; + mosaic[i].coords.crpix2 = 0.0; + mosaic[i].coords.cdelt1 = 1.0 / 3600.0; + mosaic[i].coords.cdelt2 = 1.0 / 3600.0; + mosaic[i].coords.pc1_1 = 1.0; + mosaic[i].coords.pc2_2 = 1.0; + mosaic[i].coords.pc1_2 = 0.0; + mosaic[i].coords.pc2_1 = 0.0; + + mosaic[i].Mcal = Mcal / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].Xm = Xm / MosaicN_Image[i]; + } + return; } @@ -527,8 +799,11 @@ if (!MOSAIC_ZEROPT) return (FALSE); + // if we are calibrating by mosaic, redefine myDet == on one of my mosaics int Nmatch = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { + catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise + if (TimeSelect) { if (catalog[i].measureT[j].t < TSTART) continue; @@ -582,4 +857,6 @@ } + // this measurement is on one of my mosaics, mark it as mine. + catalog[cat].measureT[meas].myDet = TRUE; MeasureToMosaic[cat][meas] = mosID; Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 36623) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 36630) @@ -908,4 +908,43 @@ } +int markObjects (Catalog *catalog, int Ncatalog) { + + int i, n; + off_t j; + + // How strongly do I own this object? + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + int nOwn = 0; + int m = catalog[i].averageT[j].measureOffset; + for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) { + if (!catalog[i].measureT[m+n].myDet) continue; + nOwn ++; + } + catalog[i].averageT[j].nOwn = nOwn; + } + } + return TRUE; +} + +int dumpObjects (char *filename, Catalog *catalog, int Ncatalog) { + + int i, n; + off_t j; + + FILE *ftest = fopen (filename, "w"); + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + int m = catalog[i].averageT[j].measureOffset; + for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) { + fprintf (ftest, "%08x %08x %10.6f %10.6f %3d %1d\n", catalog[i].averageT[j].catID, catalog[i].averageT[j].objID, catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].averageT[j].Nmeasure, catalog[i].measureT[m+n].myDet); + } + } + } + fclose (ftest); + return TRUE; +} + void clean_stars (Catalog *catalog, int Ncatalog) { Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 36623) +++ trunk/Ohana/src/relphot/src/args.c (revision 36630) @@ -1,8 +1,538 @@ # include "relphot.h" -int args (int argc, char **argv) { +RelphotMode args (int argc, char **argv) { int N; double trange; + + /* define time */ + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &TSTART)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (MODE_ERROR); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + if (!ohana_str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (MODE_ERROR); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* specify region file by name (eg n0000/0000.00) */ + UserCatalog = NULL; + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserCatalog = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + USE_BASIC_CHECK = FALSE; + if ((N = get_argument (argc, argv, "-basic-image-search"))) { + remove_argument (N, &argc, argv); + USE_BASIC_CHECK = TRUE; + } + + USE_FULL_OVERLAP = TRUE; + if ((N = get_argument (argc, argv, "-sloppy-image-overlap"))) { + remove_argument (N, &argc, argv); + USE_FULL_OVERLAP = FALSE; + } + + SET_MREL_VERSION = 1; + if ((N = get_argument (argc, argv, "-set-mrel-version"))) { + remove_argument (N, &argc, argv); + SET_MREL_VERSION = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE2 = VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + NTHREADS = 0; + if ((N = get_argument (argc, argv, "-threads"))) { + remove_argument (N, &argc, argv); + NTHREADS = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + // XXX for the moment, make this selection manual. it needs to be automatic + // based on the state of the SkyTable + HOST_ID = 0; + PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // relphot will simply list the remote command and wait for the user to signal completion + PARALLEL_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-manual"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the relphot_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; + remove_argument (N, &argc, argv); + } + + PLOTSTUFF = FALSE; + if ((N = get_argument (argc, argv, "-plot"))) { + PLOTSTUFF = TRUE; + remove_argument (N, &argc, argv); + } + + PLOTDELAY = 500000; + if ((N = get_argument (argc, argv, "-plotdelay"))) { + remove_argument (N, &argc, argv); + PLOTDELAY = 1e6*atof(argv[N]); + PLOTSTUFF = TRUE; // always turn on plotting if i request a plot delay + remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-outroot"))) { + remove_argument (N, &argc, argv); + OUTROOT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } else { + OUTROOT = strcreate ("relphot"); + } + + strcpy (STATMODE, "WT_MEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + NLOOP = 8; + if ((N = get_argument (argc, argv, "-n"))) { + remove_argument (N, &argc, argv); + NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-nloop"))) { + remove_argument (N, &argc, argv); + NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + NGRID = 8; + if ((N = get_argument (argc, argv, "-ngrid"))) { + remove_argument (N, &argc, argv); + NGRID = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + RESET_ZEROPTS = FALSE; + if ((N = get_argument (argc, argv, "-reset-zpts"))) { + remove_argument (N, &argc, argv); + RESET_ZEROPTS = TRUE; + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + UPDATE_CATFORMAT = NULL; + if ((N = get_argument (argc, argv, "-update-catformat"))) { + remove_argument (N, &argc, argv); + UPDATE_CATFORMAT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + SAVE_IMAGE_UPDATES = TRUE; + if ((N = get_argument (argc, argv, "-skip-image-updates"))) { + remove_argument (N, &argc, argv); + SAVE_IMAGE_UPDATES = FALSE; + } + + MaxDensityUse = FALSE; + if ((N = get_argument (argc, argv, "-max-density"))) { + remove_argument (N, &argc, argv); + MaxDensityValue = atof(argv[N]); + remove_argument (N, &argc, argv); + MaxDensityUse = TRUE; + } + + CLOUD_TOLERANCE = 0.02; + if ((N = get_argument (argc, argv, "-cloud-limit"))) { + remove_argument (N, &argc, argv); + CLOUD_TOLERANCE = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + // XXX should we load a tree from CATDIR by default? + // NOTE: a given catdir needs an appropriate boundary tree + BOUNDARY_TREE = NULL; + if ((N = get_argument (argc, argv, "-boundary-tree"))) { + remove_argument (N, &argc, argv); + BOUNDARY_TREE = strcreate(argv[N]); + load_tess (BOUNDARY_TREE); + remove_argument (N, &argc, argv); + } + + SHOW_PARAMS = FALSE; + if ((N = get_argument (argc, argv, "-params"))) { + remove_argument (N, &argc, argv); + SHOW_PARAMS = TRUE; + } + + PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0; + if ((N = get_argument (argc, argv, "-plrange"))) { + remove_argument (N, &argc, argv); + PlotMmin = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotMmax = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotdMmin = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotdMmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* XXX this argument used to do two things: specify the camera name and tell the analysis to + calculate a common zero point for a single mosaic. I've moved the MOSAICNAME concept into the + config system, but need to use this argument to specify that the mosaic zeropoints should be + calculated. */ + MOSAIC_ZEROPT = FALSE; + if ((N = get_argument (argc, argv, "-mosaic"))) { + remove_argument (N, &argc, argv); + MOSAIC_ZEROPT = TRUE; + if (!strcasecmp (MOSAICNAME, "none")) { + fprintf (stderr, "mosaic astrometry selected by MOSAICNAME not defined\n"); + exit (2); + } + } + + FREEZE_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-imfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_IMAGES = TRUE; + } + + FREEZE_MOSAICS = FALSE; + if ((N = get_argument (argc, argv, "-mosfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_MOSAICS = TRUE; + } + + // USE_GRID is not valid for all cases, probably should be its own mode... + USE_GRID = FALSE; + if ((N = get_argument (argc, argv, "-grid"))) { + remove_argument (N, &argc, argv); + USE_GRID = TRUE; + } + + KEEP_UBERCAL = TRUE; + if ((N = get_argument (argc, argv, "-reset-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = FALSE; + } + + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } + + AreaSelect = FALSE; + if ((N = get_argument (argc, argv, "-area"))) { + remove_argument (N, &argc, argv); + AreaXmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaXmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaSelect = TRUE; + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; + } + + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = TRUE; + init_synthetic_mags(); + } + + refPhotcode = NULL; + if ((N = get_argument (argc, argv, "-refcode"))) { + remove_argument (N, &argc, argv); + refPhotcode = GetPhotcodebyName (argv[N]); + if (!refPhotcode) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); + exit (1); + } + remove_argument (N, &argc, argv); + } + + REGION_FILE = NULL; + if ((N = get_argument (argc, argv, "-region-hosts"))) { + remove_argument (N, &argc, argv); + REGION_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + REGION_HOST_ID = 0; + if ((N = get_argument (argc, argv, "-region-hostID"))) { + remove_argument (N, &argc, argv); + REGION_HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + RelphotMode mode = MODE_ERROR; + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + mode = UPDATE_IMAGES; + } + if ((N = get_argument (argc, argv, "-averages"))) { + remove_argument (N, &argc, argv); + mode = UPDATE_AVERAGES; + } + if ((N = get_argument (argc, argv, "-apply-offsets"))) { + remove_argument (N, &argc, argv); + mode = APPLY_OFFSETS; + } + IMAGE_TABLE = NULL; + if ((N = get_argument (argc, argv, "-parallel-images"))) { + remove_argument (N, &argc, argv); + mode = PARALLEL_IMAGES; + if (N >= argc) relphot_usage(); + IMAGE_TABLE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + if (!REGION_FILE) relphot_usage(); + } + + PARALLEL_REGIONS_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-regions"))) { + remove_argument (N, &argc, argv); + mode = PARALLEL_REGIONS; + if (!REGION_FILE) relphot_usage(); + if ((N = get_argument (argc, argv, "-parallel-regions-manual"))) { + remove_argument (N, &argc, argv); + PARALLEL_REGIONS_MANUAL = TRUE; + } + } + + switch (mode) { + case UPDATE_AVERAGES: + if (argc != 1) relphot_usage(); + break; + + case UPDATE_IMAGES: + case PARALLEL_IMAGES: + case PARALLEL_REGIONS: + PhotcodeList = strcreate (argv[1]); + photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes + remove_argument (1, &argc, argv); + break; + + default: + fprintf (stderr, "no valid mode selected\n"); + relphot_usage(); + break; + } + if (argc != 1) relphot_usage (); + + return mode; +} + +int args_client (int argc, char **argv) { + + int N; + double trange; + + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + HOST_ID = 0; + if ((N = get_argument (argc, argv, "-hostID"))) { + remove_argument (N, &argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOST_ID) relphot_client_usage(); + + HOSTDIR = NULL; + if ((N = get_argument (argc, argv, "-hostdir"))) { + remove_argument (N, &argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOSTDIR) relphot_client_usage(); + + IMAGES = NULL; // used in -update mode + BCATALOG = NULL; // used in -load mode + MODE = MODE_NONE; + if ((N = get_argument (argc, argv, "-load"))) { + MODE = MODE_LOAD; + remove_argument (N, &argc, argv); + BCATALOG = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-update-catalogs"))) { + if (MODE) { + fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); + relphot_client_usage(); + } + MODE = MODE_UPDATE; + remove_argument (N, &argc, argv); + IMAGES = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-update-objects"))) { + if (MODE) { + fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); + relphot_client_usage(); + } + MODE = MODE_UPDATE_OBJECTS; + remove_argument (N, &argc, argv); + } + if (!MODE) relphot_client_usage(); + + strcpy (STATMODE, "WT_MEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + BOUNDARY_TREE = NULL; + if ((N = get_argument (argc, argv, "-boundary-tree"))) { + remove_argument (N, &argc, argv); + BOUNDARY_TREE = strcreate(argv[N]); + load_tess (BOUNDARY_TREE); + remove_argument (N, &argc, argv); + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + SET_MREL_VERSION = 1; + if ((N = get_argument (argc, argv, "-set-mrel-version"))) { + remove_argument (N, &argc, argv); + SET_MREL_VERSION = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE2 = VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + RESET_ZEROPTS = FALSE; + if ((N = get_argument (argc, argv, "-reset-zpts"))) { + remove_argument (N, &argc, argv); + RESET_ZEROPTS = TRUE; + } + + KEEP_UBERCAL = TRUE; + if ((N = get_argument (argc, argv, "-reset-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = FALSE; + } /* define time */ @@ -33,153 +563,12 @@ } - /* specify portion of the sky */ - UserPatch.Rmin = 0; - UserPatch.Rmax = 360; - UserPatch.Dmin = -90; - UserPatch.Dmax = +90; - if ((N = get_argument (argc, argv, "-region"))) { - remove_argument (N, &argc, argv); - UserPatch.Rmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Rmax = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmax = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - /* specify region file by name (eg n0000/0000.00) */ - UserCatalog = NULL; - if ((N = get_argument (argc, argv, "-catalog"))) { - remove_argument (N, &argc, argv); - UserCatalog = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - - USE_BASIC_CHECK = FALSE; - if ((N = get_argument (argc, argv, "-basic-image-search"))) { - remove_argument (N, &argc, argv); - USE_BASIC_CHECK = TRUE; - } - - USE_FULL_OVERLAP = TRUE; - if ((N = get_argument (argc, argv, "-sloppy-image-overlap"))) { - remove_argument (N, &argc, argv); - USE_FULL_OVERLAP = FALSE; - } - - SET_MREL_VERSION = 1; - if ((N = get_argument (argc, argv, "-set-mrel-version"))) { - remove_argument (N, &argc, argv); - SET_MREL_VERSION = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - VERBOSE = VERBOSE2 = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { - VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-vv"))) { - VERBOSE2 = VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - - NTHREADS = 0; - if ((N = get_argument (argc, argv, "-threads"))) { - remove_argument (N, &argc, argv); - NTHREADS = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - // XXX for the moment, make this selection manual. it needs to be automatic - // based on the state of the SkyTable - HOST_ID = 0; - PARALLEL = FALSE; - if ((N = get_argument (argc, argv, "-parallel"))) { - PARALLEL = TRUE; - remove_argument (N, &argc, argv); - } - // this is a test mode : rather than launching the remote jobs and waiting for completion, - // relphot will simply list the remote command and wait for the user to signal completion - PARALLEL_MANUAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-manual"))) { - PARALLEL = TRUE; // -parallel-manual implies -parallel - PARALLEL_MANUAL = TRUE; - remove_argument (N, &argc, argv); - } - // this is a test mode : rather than launching the relphot_client jobs remotely, they are - // run in serial via 'system' - PARALLEL_SERIAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-serial"))) { - if (PARALLEL_MANUAL) { - fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); - exit (1); - } - PARALLEL = TRUE; // -parallel-serial implies -parallel - PARALLEL_SERIAL = TRUE; - remove_argument (N, &argc, argv); - } - - PLOTSTUFF = FALSE; - if ((N = get_argument (argc, argv, "-plot"))) { - PLOTSTUFF = TRUE; - remove_argument (N, &argc, argv); - } - - PLOTDELAY = 500000; - if ((N = get_argument (argc, argv, "-plotdelay"))) { - remove_argument (N, &argc, argv); - PLOTDELAY = 1e6*atof(argv[N]); - PLOTSTUFF = TRUE; // always turn on plotting if i request a plot delay - remove_argument (N, &argc, argv); - } - - if ((N = get_argument (argc, argv, "-outroot"))) { - remove_argument (N, &argc, argv); - OUTROOT = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } else { - OUTROOT = strcreate ("relphot"); - } - - strcpy (STATMODE, "WT_MEAN"); - if ((N = get_argument (argc, argv, "-statmode"))) { - remove_argument (N, &argc, argv); - strcpy (STATMODE, argv[N]); - remove_argument (N, &argc, argv); - } - - NLOOP = 8; - if ((N = get_argument (argc, argv, "-n"))) { - remove_argument (N, &argc, argv); - NLOOP = atof (argv[N]); - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-nloop"))) { - remove_argument (N, &argc, argv); - NLOOP = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - NGRID = 8; - if ((N = get_argument (argc, argv, "-ngrid"))) { - remove_argument (N, &argc, argv); - NGRID = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - RESET = FALSE; - if ((N = get_argument (argc, argv, "-reset"))) { - remove_argument (N, &argc, argv); - RESET = TRUE; - } - - RESET_ZEROPTS = FALSE; - if ((N = get_argument (argc, argv, "-reset-zpts"))) { - remove_argument (N, &argc, argv); - RESET_ZEROPTS = TRUE; - } + + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } UPDATE = FALSE; @@ -196,8 +585,41 @@ } - SAVE_IMAGE_UPDATES = TRUE; - if ((N = get_argument (argc, argv, "-skip-image-updates"))) { - remove_argument (N, &argc, argv); - SAVE_IMAGE_UPDATES = FALSE; + AreaSelect = FALSE; + if ((N = get_argument (argc, argv, "-area"))) { + remove_argument (N, &argc, argv); + AreaXmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaXmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaSelect = TRUE; + } + + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = FALSE; + init_synthetic_mags(); + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; } @@ -210,376 +632,8 @@ } - CLOUD_TOLERANCE = 0.02; - if ((N = get_argument (argc, argv, "-cloud-limit"))) { - remove_argument (N, &argc, argv); - CLOUD_TOLERANCE = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - // XXX should we load a tree from CATDIR by default? - // NOTE: a given catdir needs an appropriate boundary tree - BOUNDARY_TREE = NULL; - if ((N = get_argument (argc, argv, "-boundary-tree"))) { - remove_argument (N, &argc, argv); - BOUNDARY_TREE = strcreate(argv[N]); - load_tess (BOUNDARY_TREE); - remove_argument (N, &argc, argv); - } - - SHOW_PARAMS = FALSE; - if ((N = get_argument (argc, argv, "-params"))) { - remove_argument (N, &argc, argv); - SHOW_PARAMS = TRUE; - } - - PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0; - if ((N = get_argument (argc, argv, "-plrange"))) { - remove_argument (N, &argc, argv); - PlotMmin = atof (argv[N]); - remove_argument (N, &argc, argv); - PlotMmax = atof (argv[N]); - remove_argument (N, &argc, argv); - PlotdMmin = atof (argv[N]); - remove_argument (N, &argc, argv); - PlotdMmax = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - /* XXX this argument used to do two things: specify the camera name and tell the analysis to - calculate a common zero point for a single mosaic. I've moved the MOSAICNAME concept into the - config system, but need to use this argument to specify that the mosaic zeropoints should be - calculated. */ - MOSAIC_ZEROPT = FALSE; - if ((N = get_argument (argc, argv, "-mosaic"))) { - remove_argument (N, &argc, argv); - MOSAIC_ZEROPT = TRUE; - if (!strcasecmp (MOSAICNAME, "none")) { - fprintf (stderr, "mosaic astrometry selected by MOSAICNAME not defined\n"); - exit (2); - } - } - - FREEZE_IMAGES = FALSE; - if ((N = get_argument (argc, argv, "-imfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_IMAGES = TRUE; - } - - FREEZE_MOSAICS = FALSE; - if ((N = get_argument (argc, argv, "-mosfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_MOSAICS = TRUE; - } - - USE_GRID = FALSE; - if ((N = get_argument (argc, argv, "-grid"))) { - remove_argument (N, &argc, argv); - USE_GRID = TRUE; - } - - KEEP_UBERCAL = TRUE; - if ((N = get_argument (argc, argv, "-reset-ubercal"))) { - remove_argument (N, &argc, argv); - KEEP_UBERCAL = FALSE; - } - - MIN_ERROR = 0.001; - if ((N = get_argument (argc, argv, "-minerror"))) { - remove_argument (N, &argc, argv); - MIN_ERROR = atof (argv[N]); - remove_argument (N, &argc, argv); - /* require MIN_ERROR > 0 */ - } - - AreaSelect = FALSE; - if ((N = get_argument (argc, argv, "-area"))) { - remove_argument (N, &argc, argv); - AreaXmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaXmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaSelect = TRUE; - } - - ImagSelect = FALSE; - if ((N = get_argument (argc, argv, "-instmag"))) { - remove_argument (N, &argc, argv); - ImagMin = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagMax = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagSelect = TRUE; - } - - DophotSelect = FALSE; - if ((N = get_argument (argc, argv, "-dophot"))) { - remove_argument (N, &argc, argv); - DophotValue = atof (argv[N]); - remove_argument (N, &argc, argv); - DophotSelect = TRUE; - } - - SyntheticPhotometry = FALSE; - if ((N = get_argument (argc, argv, "-synthphot"))) { - remove_argument (N, &argc, argv); - SyntheticPhotometry = TRUE; - init_synthetic_mags(); - } - - refPhotcode = NULL; - if ((N = get_argument (argc, argv, "-refcode"))) { - remove_argument (N, &argc, argv); - refPhotcode = GetPhotcodebyName (argv[N]); - if (!refPhotcode) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); - exit (1); - } - remove_argument (N, &argc, argv); - } - - UpdateAverages = FALSE; - if ((N = get_argument (argc, argv, "-averages"))) { - remove_argument (N, &argc, argv); - UpdateAverages = TRUE; - } - - ApplyOffsets = FALSE; - if ((N = get_argument (argc, argv, "-apply-offsets"))) { - remove_argument (N, &argc, argv); - ApplyOffsets = TRUE; - } - - if (UpdateAverages && (argc == 1)) return TRUE; - if (argc != 2) relphot_usage (); + if ((MODE == MODE_UPDATE_OBJECTS) && (argc == 1)) return TRUE; + if (argc != 2) relphot_client_usage (); return TRUE; } -int args_client (int argc, char **argv) { - - int N; - double trange; - - // by definition, the client is not parallel - PARALLEL = FALSE; - PARALLEL_MANUAL = FALSE; - PARALLEL_SERIAL = FALSE; - - HOST_ID = 0; - if ((N = get_argument (argc, argv, "-hostID"))) { - remove_argument (N, &argc, argv); - HOST_ID = atoi (argv[N]); - remove_argument (N, &argc, argv); - } - if (!HOST_ID) relphot_client_usage(); - - HOSTDIR = NULL; - if ((N = get_argument (argc, argv, "-hostdir"))) { - remove_argument (N, &argc, argv); - HOSTDIR = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if (!HOSTDIR) relphot_client_usage(); - - IMAGES = NULL; // used in -update mode - BCATALOG = NULL; // used in -load mode - MODE = MODE_NONE; - if ((N = get_argument (argc, argv, "-load"))) { - MODE = MODE_LOAD; - remove_argument (N, &argc, argv); - BCATALOG = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-update-catalogs"))) { - if (MODE) { - fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); - relphot_client_usage(); - } - MODE = MODE_UPDATE; - remove_argument (N, &argc, argv); - IMAGES = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-update-objects"))) { - if (MODE) { - fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n"); - relphot_client_usage(); - } - MODE = MODE_UPDATE_OBJECTS; - remove_argument (N, &argc, argv); - } - if (!MODE) relphot_client_usage(); - - strcpy (STATMODE, "WT_MEAN"); - if ((N = get_argument (argc, argv, "-statmode"))) { - remove_argument (N, &argc, argv); - strcpy (STATMODE, argv[N]); - remove_argument (N, &argc, argv); - } - - BOUNDARY_TREE = NULL; - if ((N = get_argument (argc, argv, "-boundary-tree"))) { - remove_argument (N, &argc, argv); - BOUNDARY_TREE = strcreate(argv[N]); - load_tess (BOUNDARY_TREE); - remove_argument (N, &argc, argv); - } - - /* specify portion of the sky */ - UserPatch.Rmin = 0; - UserPatch.Rmax = 360; - UserPatch.Dmin = -90; - UserPatch.Dmax = +90; - if ((N = get_argument (argc, argv, "-region"))) { - remove_argument (N, &argc, argv); - UserPatch.Rmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Rmax = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmin = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmax = atof (argv[N]); - remove_argument (N, &argc, argv); - } - - SET_MREL_VERSION = 1; - if ((N = get_argument (argc, argv, "-set-mrel-version"))) { - remove_argument (N, &argc, argv); - SET_MREL_VERSION = atof(argv[N]); - remove_argument (N, &argc, argv); - } - - VERBOSE = VERBOSE2 = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { - VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-vv"))) { - VERBOSE2 = VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - - RESET = FALSE; - if ((N = get_argument (argc, argv, "-reset"))) { - remove_argument (N, &argc, argv); - RESET = TRUE; - } - - RESET_ZEROPTS = FALSE; - if ((N = get_argument (argc, argv, "-reset-zpts"))) { - remove_argument (N, &argc, argv); - RESET_ZEROPTS = TRUE; - } - - KEEP_UBERCAL = TRUE; - if ((N = get_argument (argc, argv, "-reset-ubercal"))) { - remove_argument (N, &argc, argv); - KEEP_UBERCAL = FALSE; - } - - /* define time */ - TimeSelect = FALSE; - if ((N = get_argument (argc, argv, "-time"))) { - remove_argument (N, &argc, argv); - if (!ohana_str_to_time (argv[N], &TSTART)) { - fprintf (stderr, "ERROR: syntax error\n"); - return (FALSE); - } - remove_argument (N, &argc, argv); - if (!ohana_str_to_dtime (argv[N], &trange)) { - if (!ohana_str_to_time (argv[N], &TSTOP)) { - fprintf (stderr, "ERROR: syntax error\n"); - return (FALSE); - } - } else { - if (trange < 0) { - trange = fabs (trange); - TSTOP = TSTART; - TSTART -= trange; - } else { - TSTOP = TSTART + trange; - } - } - remove_argument (N, &argc, argv); - TimeSelect = TRUE; - } - - - MIN_ERROR = 0.001; - if ((N = get_argument (argc, argv, "-minerror"))) { - remove_argument (N, &argc, argv); - MIN_ERROR = atof (argv[N]); - remove_argument (N, &argc, argv); - /* require MIN_ERROR > 0 */ - } - - UPDATE = FALSE; - if ((N = get_argument (argc, argv, "-update"))) { - remove_argument (N, &argc, argv); - UPDATE = TRUE; - } - - UPDATE_CATFORMAT = NULL; - if ((N = get_argument (argc, argv, "-update-catformat"))) { - remove_argument (N, &argc, argv); - UPDATE_CATFORMAT = strcreate (argv[N]); - remove_argument (N, &argc, argv); - } - - AreaSelect = FALSE; - if ((N = get_argument (argc, argv, "-area"))) { - remove_argument (N, &argc, argv); - AreaXmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaXmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaSelect = TRUE; - } - - SyntheticPhotometry = FALSE; - if ((N = get_argument (argc, argv, "-synthphot"))) { - remove_argument (N, &argc, argv); - SyntheticPhotometry = FALSE; - init_synthetic_mags(); - } - - ImagSelect = FALSE; - if ((N = get_argument (argc, argv, "-instmag"))) { - remove_argument (N, &argc, argv); - ImagMin = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagMax = atof (argv[N]); - remove_argument (N, &argc, argv); - ImagSelect = TRUE; - } - - DophotSelect = FALSE; - if ((N = get_argument (argc, argv, "-dophot"))) { - remove_argument (N, &argc, argv); - DophotValue = atof (argv[N]); - remove_argument (N, &argc, argv); - DophotSelect = TRUE; - } - - MaxDensityUse = FALSE; - if ((N = get_argument (argc, argv, "-max-density"))) { - remove_argument (N, &argc, argv); - MaxDensityValue = atof(argv[N]); - remove_argument (N, &argc, argv); - MaxDensityUse = TRUE; - } - - if ((MODE == MODE_UPDATE_OBJECTS) && (argc == 1)) return TRUE; - if (argc != 2) relphot_client_usage (); - - return TRUE; -} - Index: trunk/Ohana/src/relphot/src/assign_images.c =================================================================== --- trunk/Ohana/src/relphot/src/assign_images.c (revision 36630) +++ trunk/Ohana/src/relphot/src/assign_images.c (revision 36630) @@ -0,0 +1,289 @@ +# include "relphot.h" + +// This function generates a subset of the images based on selections. Input db has already +// been loaded with the raw fits table data +int assign_images (FITS_DB *db, RegionHostTable *regionHosts) { + + off_t Nimage; + + INITTIME; + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + Image *image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + MARKTIME("convert image table to internal structure: %f sec\n", dtime); + + // *** NOTE : for the moment, regions must be in the range 0 - 360, -90 - +90 + + // generate the chip match here so we can define the mosaic centers (if needed) + BuildChipMatch (image, Nimage); + MARKTIME("build chip match for %d images: %f sec\n", (int) Nimage, dtime); + + if (MOSAIC_ZEROPT) { + makeMosaics (image, Nimage); + + // center coords and Mcal, dMcal, Xm for the mosaics + setMosaicCenters (image, Nimage); + MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); + } + + // register the image array with ImageOps.c for later getimageByID calls + initImages (image, NULL, Nimage); + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + // for each regionHost, select images which are contained by the region + // even faster would be to use a tree to get to the real regions... + select_images_hostregion (regionHosts, image, Nimage); + + return TRUE; +} + +# define D_NIMAGE 1000 + +// assign images to the region hosts; at the end, each host will have its list of images +int select_images_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage) { + + int ecode, found; + off_t i, j; + + // INITTIME; + + for (i = 0; i < regionHosts->Nhosts; i++) { + regionHosts->hosts[i].Nimage = 0; + regionHosts->hosts[i].NIMAGE = D_NIMAGE; + ALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); + } + + for (j = 0; j < Nimage; j++) { + + /* exclude images by photcode */ + ecode = GetPhotcodeEquivCodebyCode (image[j].photcode); + found = FALSE; + int Ns; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) continue; + + /* exclude images by time */ + if (TimeSelect) { + if (image[j].tzero < TSTART) continue; + if (image[j].tzero > TSTOP) continue; + } + + // do not include DIS (PHU-level mosaics) in the output list + if (!strcmp(&image[j].coords.ctype[4], "-DIS")) continue; + + // match the image to its corresponding ASTROMETRIC mosaic (not the same as the Mosaic above) + if (!FindMosaicForImage (image, Nimage, j)) { + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + // Exclude images with crazy astrometry + // XXX NOTE : this is gpc1-specific + { + double dP1 = hypot(image[j].coords.pc1_1, image[j].coords.pc1_2); + double dP2 = hypot(image[j].coords.pc2_1, image[j].coords.pc2_2); + if (fabs(dP1 - 1.0) > 0.02) continue; + if (fabs(dP2 - 1.0) > 0.02) continue; + + double X00, Y00, X10, Y10, X01, Y01; + XY_to_LM (&X00, &Y00, 0.0, 0.0, &image[j].coords); + XY_to_LM (&X10, &Y10, image[j].NX, 0.0, &image[j].coords); + XY_to_LM (&X01, &Y01, 0.0, image[j].NY, &image[j].coords); + double dS0 = hypot ((X00 - X10), (Y00 - Y10)); + double dS1 = hypot ((X00 - X01), (Y00 - Y01)); + if (dS0 > 6000) continue; + if (dS1 > 6500) continue; + } + + // use a reference coordinate for each image to assign to hosts + // define image center - note the DIS images (mosaic phu) are special, but we have + // already excluded them above. we also save the image centers for reference + double Xc, Yc; + double Rc, Dc; + + Xc = 0.5*image[j].NX; + Yc = 0.5*image[j].NY; + + XY_to_RD (&Rc, &Dc, Xc, Yc, &image[j].coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + image[j].RAo = Rc; + image[j].DECo = Dc; + + // NOTE: if we are NOT using mosaic centers, we don't need to keep chips with centers + // on opposite sides of 0,360 together. in which case, it is OK for the range of Rc + // to be 0.0 to 360.0. also NOTE: RAo,DECo are only used for debugging reference. + + if (MOSAIC_ZEROPT) { + // use the coords of the associated mosaic to select + Mosaic *mosaic = getMosaicForImage (j); + Rc = mosaic->coords.crval1; + Dc = mosaic->coords.crval2; + // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the + // chips are located. but, for host assignment, we rationalize to 0.0 - 360.0 + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + } + + i = find_host_for_coords (regionHosts, Rc, Dc); + + if (i == -1) continue; + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // image bounds are defined for a range centered on the mosaic center thus, a mosaic + // with center 0.5 will have chips bounds ranging from ~ -1.5 to +2.5 or so, while a mosaic + // with center 359.5 will have chips bounds ranging from ~ 357.5 to 361.5 or so + double Rmin, Rmax, Dmin, Dmax; + calculate_image_bounds (&image[j], &Rmin, &Rmax, &Dmin, &Dmax, Rc); + + host->RminCat = MIN(Rmin, host->RminCat); + host->RmaxCat = MAX(Rmax, host->RmaxCat); + host->DminCat = MIN(Dmin, host->DminCat); + host->DmaxCat = MAX(Dmax, host->DmaxCat); + + // regionHosts needs to have the full outer boundary + // (so reload_catalogs covers the correct region) + regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); + regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); + regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); + regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); + + // this is a bit memory expensive : I am making a complete copy of the image table here + off_t Nsubset = regionHosts->hosts[i].Nimage; + regionHosts->hosts[i].image[Nsubset] = image[j]; + // regionHosts->hosts[i].line_number[Nsubset] = j; + + regionHosts->hosts[i].Nimage ++; + if (regionHosts->hosts[i].Nimage == regionHosts->hosts[i].NIMAGE) { + regionHosts->hosts[i].NIMAGE += D_NIMAGE; + REALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); + // REALLOCATE (regionHosts->hosts[i].line_number, off_t, regionHosts->hosts[i].NIMAGE); + } + } + + return TRUE; +} + +double Xf[] = {0.0, 1.0, 0.0, 1.0}; +double Yf[] = {0.0, 0.0, 1.0, 1.0}; + +int calculate_image_bounds (Image *image, double *rmin, double *rmax, double *dmin, double *dmax, double Rmid) { + + int n; + + double Rmin = 360.0; + double Rmax = 0.0; + double Dmin = +90.0; + double Dmax = -90.0; + + // define image corners + for (n = 0; n < 4; n++) { + double Xc, Yc, Rc, Dc; + Xc = Xf[n]*image->NX; + Yc = Yf[n]*image->NY; + XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, Rmid); + + Rmin = MIN (Rmin, Rc); + Rmax = MAX (Rmax, Rc); + Dmin = MIN (Dmin, Dc); + Dmax = MAX (Dmax, Dc); + } + + *rmin = Rmin; + *rmax = Rmax; + *dmin = Dmin; + *dmax = Dmax; + + return TRUE; +} + +int calculate_host_image_bounds (RegionHostTable *regionHosts, double Rmid) { + + int i, n; + off_t j; + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // XXX clear the chip match? + + BuildChipMatch (host->image, host->Nimage); + + double Rmin = 720.0; + double Rmax = -360.0; + double Dmin = +90.0; + double Dmax = -90.0; + + for (j = 0; j < host->Nimage; j++) { + + Image *image = &host->image[j]; + + if (!FindMosaicForImage (host->image, host->Nimage, j)) { + fprintf (stderr, "missing astrometry? programming error?\n"); + abort (); + } + + // define image corners + for (n = 0; n < 4; n++) { + double Xc, Yc, Rc, Dc; + Xc = Xf[n]*image->NX; + Yc = Yf[n]*image->NY; + XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); + Rc = ohana_normalize_angle_to_midpoint (Rc, Rmid); + + Rmin = MIN (Rmin, Rc); + Dmin = MIN (Dmin, Rc); + + Rmax = MAX (Rmax, Rc); + Dmax = MAX (Dmax, Rc); + } + + } + + // RminCat, RmaxCat may extended beyond 0.0 - 360.0 + host->RminCat = Rmin; + host->DminCat = Dmin; + host->RmaxCat = Rmax; + host->DmaxCat = Dmax; + + // regionHosts needs to have the full outer boundary + // (so reload_catalogs covers the correct region) + regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); + regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); + regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); + regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); + } + return TRUE; +} + + +// Rc is in range 0.0 - 360.0, hosts Rmin,Rmax in range 0.0 - 360.0 +int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc) { + + int i; + + if (isnan(Rc)) return -1; + if (isnan(Dc)) return -1; + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + if (Rc < host->Rmin) continue; + if (Rc >= host->Rmax) continue; + if (Dc < host->Dmin) continue; + if (Dc >= host->Dmax) continue; + + return i; + } + return -1; +} + Index: trunk/Ohana/src/relphot/src/client_logger.c =================================================================== --- trunk/Ohana/src/relphot/src/client_logger.c (revision 36623) +++ trunk/Ohana/src/relphot/src/client_logger.c (revision 36630) @@ -5,13 +5,13 @@ static FILE *logfile = NULL; -int client_logger_init () { +int client_logger_init (char *dirname) { char filename[DVO_MAX_PATH]; - snprintf (filename, DVO_MAX_PATH, "%s/log.rlpc.XXXXXX", HOSTDIR); + snprintf (filename, DVO_MAX_PATH, "%s/log.rlpc.XXXXXX", dirname); int fd = mkstemp (filename); if (fd == -1) { - fprintf (stderr, "failed to open client logger, exiting\n"); + fprintf (stderr, "failed to open client logger %s, exiting\n", filename); exit (50); } @@ -32,4 +32,6 @@ vfprintf (logfile, format, argp); va_end (argp); + + fflush (logfile); return TRUE; } Index: trunk/Ohana/src/relphot/src/help.c =================================================================== --- trunk/Ohana/src/relphot/src/help.c (revision 36623) +++ trunk/Ohana/src/relphot/src/help.c (revision 36630) @@ -2,7 +2,11 @@ void relphot_usage (void) { - fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -region RA RA DEC DEC\n"); - fprintf (stderr, " or: relphot (photcodes) -catalog (name)\n"); - fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); + fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -images\n"); + fprintf (stderr, " or: relphot -averages\n"); + fprintf (stderr, " or: relphot -apply-offsets\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-regions -region-hosts (RegionFile)\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-images (ImageTable) -region-hosts (RegionFile)\n\n"); + fprintf (stderr, " regions: -region RA RA DEC DEC)\n"); + fprintf (stderr, " or: -catalog (name)\n"); fprintf (stderr, " use -h for more usage information\n"); exit (2); @@ -18,7 +22,11 @@ show_help: - fprintf (stderr, "ERROR: USAGE: relphot (photcode) -region RA RA DEC DEC\n"); - fprintf (stderr, " or: relphot (photcode) -catalog (name)\n"); - fprintf (stderr, " or: relphot -averages -region RA RA DEC DEC\n"); + fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -images\n"); + fprintf (stderr, " or: relphot -averages\n"); + fprintf (stderr, " or: relphot -apply-offsets\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-regions -region-hosts (RegionFile)\n"); + fprintf (stderr, " or: relphot (photcodes) -parallel-images (ImageTable) -region-hosts (RegionFile)\n\n"); + fprintf (stderr, " regions: -region RA RA DEC DEC)\n"); + fprintf (stderr, " or: -catalog (name)\n"); fprintf (stderr, " options: \n"); fprintf (stderr, " -time (start) (stop)\n"); Index: trunk/Ohana/src/relphot/src/indexCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/indexCatalog.c (revision 36630) +++ trunk/Ohana/src/relphot/src/indexCatalog.c (revision 36630) @@ -0,0 +1,64 @@ +# include "relphot.h" + +static int catIDmax = 0; +static int *catIDseq = NULL; +static int *objIDmax = NULL; +static int **objIDseq = NULL; + +int indexCatalogs (Catalog *catalog, int Ncatalog) { + + int i; + off_t j; + + // find the max value of catID + for (i = 0; i < Ncatalog; i++) { + catIDmax = MAX (catalog[i].catID, catIDmax); + } + + ALLOCATE (catIDseq, int, catIDmax + 1); + for (i = 0; i < catIDmax + 1; i++) { + catIDseq[i] = -1; + } + + for (i = 0; i < Ncatalog; i++) { + int catID = catalog[i].catID; + catIDseq[catID] = i; + } + + ALLOCATE (objIDmax, int, Ncatalog); + ALLOCATE (objIDseq, int *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + objIDmax[i] = 0; + for (j = 0; j < catalog[i].Naverage; j++) { + objIDmax[i] = MAX (catalog[i].averageT[j].objID, objIDmax[i]); + } + + ALLOCATE (objIDseq[i], int, objIDmax[i] + 1); + for (j = 0; j < objIDmax[i] + 1; j++) { + objIDseq[i][j] = -1; + } + + for (j = 0; j < catalog[i].Naverage; j++) { + int objID = catalog[i].averageT[j].objID; + objIDseq[i][objID] = j; + } + } + return TRUE; +} + +int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq) { + + if (catID > catIDmax) return FALSE; + + int cat = catIDseq[catID]; + if (cat < 1) return FALSE; + + if (objID > objIDmax[cat]) return FALSE; + + int obj = objIDseq[cat][objID]; + if (obj < 1) return FALSE; + + *catSeq = cat; + *objSeq = obj; + return TRUE; +} Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 36623) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 36630) @@ -1,16 +1,13 @@ # include "relphot.h" -void initialize (int argc, char **argv) { +RelphotMode initialize (int argc, char **argv) { relphot_help (argc, argv); ConfigInit (&argc, argv); - args (argc, argv); + RelphotMode mode = args (argc, argv); + if (!mode) exit (2); - if (!UpdateAverages) { - // load the list of photcodes into the globals (photcodes, Nphotcodes) - // only a single remaining argument in this mode (the list of photcodes, eg g,r,i) - PhotcodeList = strcreate (argv[1]); - photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes - } else { + // UPDATE_AVERAGES always operates on all photcodes? + if (mode == UPDATE_AVERAGES) { char tmpline1[256]; int Ns; @@ -67,4 +64,6 @@ for (B = 0; A == time(NULL); B++); srand48(B); + + return mode; } Index: trunk/Ohana/src/relphot/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 36630) +++ trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 36630) @@ -0,0 +1,117 @@ +# include "relphot.h" +# define DEBUG 0 + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + +int launch_region_hosts (RegionHostTable *regionHosts) { + + int i; + + // clear the I/O files + for (i = 0; i < regionHosts->Nhosts; i++) { + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync"); + truncate (syncfile, 0); + free (syncfile); + + char *fitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits"); + truncate (fitsfile, 0); + free (fitsfile); + + char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync"); + truncate (imsyncfile, 0); + free (imsyncfile); + + char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits"); + truncate (imfitsfile, 0); + free (imfitsfile); + + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + truncate (loopsyncfile, 0); + free (loopsyncfile); + } + + for (i = 0; i < regionHosts->Nhosts; i++) { + + RegionHostInfo *host = ®ionHosts->hosts[i]; + + // communication files: + // subset images per host : CATDIR/Image.HOSTNAME.fits + char filename[1024]; + snprintf (filename, 1024, "%s/Image.%d.fits", CATDIR, host->hostID); + + // write the image subset for this host + ImageTableSave (filename, host->image, host->Nimage); + + char command[1024]; + snprintf (command, 1024, "relphot %s -parallel-images %s -region-hosts %s -region-hostID %d -D CATDIR %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f", + PhotcodeList, filename, REGION_FILE, host->hostID, CATDIR, host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR); + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (RESET_ZEROPTS) strextend (command, "-reset-zpts"); + if (!KEEP_UBERCAL) strextend (command, "-reset-ubercal"); + if (DophotSelect) strextend (command, "-dophot %d", DophotValue); + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + if (SyntheticPhotometry) strextend (command, "-synthphot"); + + if (UPDATE) strextend (command, "-update"); + if (MOSAIC_ZEROPT) strextend (command, "-mosaic"); + if (FREEZE_IMAGES) strextend (command, "-imfreeze"); + if (FREEZE_MOSAICS) strextend (command, "-mosfreeze"); + if (PARALLEL) strextend (command, "-parallel"); + if (PARALLEL_MANUAL) strextend (command, "-parallel-manual"); + if (PARALLEL_SERIAL) strextend (command, "-parallel-serial"); + + // XXX deprecate this if we are happy with the new version + if (SET_MREL_VERSION != 1) strextend (command, "-set-mrel-version %d", SET_MREL_VERSION); + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_REGIONS_MANUAL) continue; + + // launch the job, then wait for it to be done loading catalogs. force the remote + // client to generate the file + char *syncfile = make_filename (CATDIR, host->hostname, host->hostID, "loadcat.sync"); + clear_sync_file (syncfile); + + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", host->hostname, command, host->stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", host->hostname, errorInfo); + exit (1); + } + host->pid = pid; // save for future reference + + // remote client is done, go ahead with next client + check_sync_file (syncfile, 1); + free (syncfile); + } + + if (PARALLEL_REGIONS_MANUAL) { + fprintf (stderr, "run the relphot_client commands above. when these are done, hit return\n"); + getchar(); + } else { + RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE); + } + + int status = TRUE; + for (i = 0; i < regionHosts->Nhosts; i++) { + status = status && (regionHosts->hosts[i].status == 0); + } + + return status; +} Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 36623) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 36630) @@ -89,7 +89,6 @@ fprintf (stderr, "using "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measurements)\n", Nstar, Nstar_total, Nmeas, Nmeas_total); - if (!hostID && (Nstar < 1)) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n"); + if (!hostID && !REGION_HOST_ID && (Nstar < 1)) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n"); // in regular relphot, we shutdown here; in relphot_client, we generate and return an empty table (for consistency) - // XXX consider only returning the populated catalogs @@ -107,4 +106,9 @@ # define DEBUG 1 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); int Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog @@ -128,5 +132,5 @@ ALLOCATE (table->hosts[i].results, char, 1024); - snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname); + snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.%s.dat", table->hosts[i].pathname, uniquer); // options / arguments that can affect relphot_client -load: Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 36623) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 36630) @@ -156,4 +156,8 @@ // load the list of hosts HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: problem with parallel host table\n"); + exit (2); + } if (BOUNDARY_TREE) { Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 36623) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 36630) @@ -3,26 +3,42 @@ int main (int argc, char **argv) { - /* get configuration info, args */ - initialize (argc, argv); + // get configuration info, args + RelphotMode mode = initialize (argc, argv); + if (!mode) exit (2); - /* the object analysis is a separate process iterating over catalogs */ - if (UpdateAverages) { - relphot_objects (0, NULL); - exit (0); + switch (mode) { + case UPDATE_IMAGES: + relphot_images (); + exit (0); + + case UPDATE_AVERAGES: + // take the current set of detections and set the mean magnitudes + relphot_objects (0, NULL); + exit (0); + + case PARALLEL_REGIONS: + // run image updates in parallel across multiple remote machines + relphot_parallel_regions (); + exit (0); + + case PARALLEL_IMAGES: + // operation on the remote machines in the PARALLEL_REGION mode + relphot_parallel_images (); + exit (0); + + case APPLY_OFFSETS: + // re-run this step from a previous attempt (assumes an existing Images.subset.dat file) + if (!PARALLEL) { + fprintf (stderr, "-apply-offsets only makes sense in an parallel context\n"); + exit (2); + } + reload_catalogs (NULL, NULL, 0, NULL); + exit (0); + + default: + fprintf (stderr, "ERROR: no valid relphot mode chosen\n"); + exit (2); } - - if (ApplyOffsets) { - // re-run this step from a previous attempt (assumes an existing Images.subset.dat file) - if (!PARALLEL) { - fprintf (stderr, "-apply-offsets only makes sense in an parallel context\n"); - exit (2); - } - reload_catalogs (NULL, NULL, 0, NULL); - exit (0); - } - - relphot_images (); - - exit (0); + fprintf (stderr, "IMPOSSIBLE: skipped out of switch?\n"); + exit (1); } - Index: trunk/Ohana/src/relphot/src/relphot_client.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_client.c (revision 36623) +++ trunk/Ohana/src/relphot/src/relphot_client.c (revision 36630) @@ -19,5 +19,5 @@ // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) initialize_client (argc, argv); - client_logger_init (); + client_logger_init (HOSTDIR); // load the current sky table (layout of all SkyRegions) Index: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 36623) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 36630) @@ -244,5 +244,7 @@ if (!UPDATE) exit (0); - /* load catalog data from region files, update Mrel include all data */ + /* Load catalog data from region files, update Mrel include all data. In a parallel + context, this function writes the image parameters as a subset table for the remote + clients */ reload_catalogs (skylist, flatcorr, 0, NULL); MARKTIME("-- updated all catalogs: %f sec\n", dtime); Index: trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 36630) +++ trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 36630) @@ -0,0 +1,165 @@ +# include "relphot.h" + +/* This function is essentially identical to relphot_images, except: + + * load the subset images saved by the master node + * distinguish detections we own (touch our images) and those we don't + * distinguish objects we own (in region) and those we don't + * update the unowned detections for owned objects from neighbor regions + * update the unowned objects for owned detections + */ + +int dumpObjects (char *filename, Catalog *catalog, int Ncatalog); + +int relphot_parallel_images () { + + int i, Ncatalog; + Catalog *catalog = NULL; + + INITTIME; + + client_logger_init (CATDIR); + + // load the RegionTable (UserRegion should not be used at this level) + RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); + int myHost = regionHosts->index[REGION_HOST_ID]; + RegionHostFindNeighbors (regionHosts, myHost); + client_logger_message ("started parallel images on %s\n", regionHosts->hosts[myHost].hostname); + + // load the subset images belonging to this host + off_t Nimage; + Image *image = ImageTableLoad (IMAGE_TABLE, &Nimage); + if (!image) { + fprintf (stderr, "ERROR loading image %s\n", IMAGE_TABLE); + exit (2); + } + + // once we have read this table, we should remove it for repeat runs + // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection? + + // XXX need to deal with mosaic vs image... + makeMosaics (image, Nimage); + + initImages (image, NULL, Nimage); + + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // UserPatch.Rmin,Rmax may have range from a few degrees < 0.0 to few degrees > 360.0. + // the following function correctly chooses the sky regions on the 0,360 boundary + SkyList *skylist = SkyListByBounds (sky, -1, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); + client_logger_message ("loaded catalog data\n"); + + // NOTE: if I let all hosts load blindly, I saturate the data clients with too many + // relphot_client requests. I need to have the master mediate this. the master + // will not launch the next remote job until this one says it is done + char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); + update_sync_file (syncfile, 1); + + // generate tables go from catID,objID -> catSeq,objSeq + indexCatalogs (catalog, Ncatalog); + client_logger_message ("indexed Catalogs\n"); + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + initMosaicBins (catalog, Ncatalog, TRUE); + initMrel (catalog, Ncatalog); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + client_logger_message ("done setting up indexes\n"); + + // dumpObjects ("test.obj.dat", catalog, Ncatalog); + markObjects (catalog, Ncatalog); + + SAVEPLOT = FALSE; + + setExclusions (catalog, Ncatalog, TRUE); + + global_stats (catalog, Ncatalog, flatcorr); + + if (PLOTSTUFF) { + plot_star_coords (catalog, Ncatalog); + // plot_mosaic_fields (catalog); + } + + // XXX : add this in the loop at various points + // if (PLOTSTUFF) plot_scatter (catalog, Ncatalog, flatcorr); + + /* determine fit values */ + client_logger_message ("starting the loops : %d \n", NLOOP); + for (i = 0; i < NLOOP; i++) { + + // set the mean stellar mags given the measurements and the image calibrations + setMrel (catalog, Ncatalog, flatcorr); + + // share mean mags for objects at the boundary (number of unowned meas > 0) + share_mean_mags (catalog, Ncatalog, regionHosts, i); + client_logger_message ("shared mean mag data : loop %d \n", i); + + // load mean mags from other region hosts + slurp_mean_mags (catalog, Ncatalog, regionHosts, i); + client_logger_message ("slurped mean mag data : loop %d \n", i); + + // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements + setMcal (catalog, FALSE, flatcorr); + setMmos (catalog, FALSE, flatcorr); + MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); + + // share image mags for images with non-zero unowned detections + share_image_mags (regionHosts, i); + client_logger_message ("shared image data : loop %d \n", i); + + slurp_image_mags (regionHosts, i); + client_logger_message ("slurped image data : loop %d \n", i); + + if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); + if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); + if ((i > 8) && (i % 8 == 5)) clean_mosaics (); + if ((i > 8) && (i % 8 == 5)) clean_images (); + + if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); + } + client_logger_message ("done with loops\n"); + + // this is a checkpoint to make sure all hosts have finished the loop above + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loop.sync"); + update_sync_file (loopsyncfile, 0); + free (loopsyncfile); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + char *loopsync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + check_sync_file (loopsync, 0); + free (loopsync); + } + client_logger_message ("all hosts are done the loops\n"); + + /* set Mcal & Mmos for bad images */ + setMcal (catalog, TRUE, flatcorr); + setMmos (catalog, TRUE, flatcorr); + MARKTIME("-- finalize Mcal values: %f sec\n", dtime); + + setMcalFinal (); // copy per-mosaic calibrations to the images + + share_image_mags (regionHosts, -1); + client_logger_message ("done with parallel images\n"); + + exit (0); +} Index: trunk/Ohana/src/relphot/src/relphot_parallel_regions.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 36630) +++ trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 36630) @@ -0,0 +1,75 @@ +# include "relphot.h" + +int relphot_parallel_regions () { + + int status; + FITS_DB db; + + INITTIME; + + // load the RegionTable (UserRegion should not be used at this level) + RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); + if (!regionHosts) { + fprintf (stderr, "ERROR: problem with region host table\n"); + exit (2); + } + + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + // register database handle with shutdown procedure + set_db (&db); + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + // lock and load the image db table + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status && UPDATE) { + fprintf (stderr, "error\n"); + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + // if the Image.dat file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("-- load image data: %f sec\n", dtime); + + /* assign the images to the different region hosts */ + if (!assign_images (&db, regionHosts)) Shutdown ("error assigning images to region hosts"); + MARKTIME("-- assign images: %f sec\n", dtime); + + /* launch processing on the parallel region hosts */ + if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts"); + + // retrieve updated image parameters from the remote hosts (also set Image.mcal) + if (!slurp_image_mags (regionHosts, -1)) Shutdown ("error loading image updates"); + + if (!UPDATE) { + dvo_image_unlock (&db); + MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime); + fprintf (stderr, "NOTE: UPDATE is OFF (results are not saved)\n"); + exit (0); + } + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByBounds (sky, -1, regionHosts->Rmin, regionHosts->Rmax, regionHosts->Dmin, regionHosts->Dmax); + + /* update catalogs (in parallel) */ + reload_catalogs (skylist, flatcorr, 0, NULL); + + // save the changes to the image parameters + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + MARKTIME ("finished relphot -parallel-regions: %f sec total\n", dtime); + + exit (0); +} + Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 36623) +++ trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 36630) @@ -255,5 +255,5 @@ // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter - float Map = PhotAper (&measure[k]); + float Map = PhotAperCat (&measure[k]); if (!isnan(Map)) { Maplist[Nap] = Map - Mcal - Mmos - Mgrid; @@ -263,5 +263,5 @@ } - float Mkron = PhotKron (&measure[k]); + float Mkron = PhotKronCat (&measure[k]); if (!isnan(Mkron)) { Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; @@ -282,4 +282,8 @@ wpsflist[Npsf] = isUbercal ? 10.0 : 1.0; + // NOTE: + // Msys is measure[i].M + zp corrections + // Mcal is image[j].Mcal + // Mmos and Mgrid are offsets for mosaic and grid // tie down reference photometry if the -refcode (code) option is selected Index: trunk/Ohana/src/relphot/src/share_image_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/share_image_mags.c (revision 36630) +++ trunk/Ohana/src/relphot/src/share_image_mags.c (revision 36630) @@ -0,0 +1,120 @@ +# include "relphot.h" + +// we are sharing image calibrations mags for all images which (a) I own and (b) which have unowned detections + +# define D_NIMAGEMAGS 1000 +int share_image_mags (RegionHostTable *regionHosts, int nloop) { + + off_t i, Nimages; + Image *images = getimages (&Nimages, NULL); + + off_t Nimage_mags = 0; + off_t NIMAGE_MAGS = D_NIMAGEMAGS; + + ImageMag *image_mags = NULL; + ALLOCATE (image_mags, ImageMag, NIMAGE_MAGS); + + for (i = 0; i < Nimages; i++) { + // XXX does this image have missing detections (does someone else need it?) + // XXX : NOTE NEED TO FIX THIS: if (imageExtra[i].Nmiss == 0) continue; + + set_image_mags (&image_mags[Nimage_mags], &images[i]); + Nimage_mags ++; + + CHECK_REALLOCATE (image_mags, ImageMag, NIMAGE_MAGS, Nimage_mags, D_NIMAGEMAGS); + } + + // write out the image_mag fits table AND write state in some file + int myHost = regionHosts->index[REGION_HOST_ID]; + char *hostname = regionHosts->hosts[myHost].hostname; + + char *imagfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagemags.fits"); + ImageMagSave (imagfile, image_mags, Nimage_mags); + free (image_mags); + free (imagfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagemags.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_image_mags (RegionHostTable *regionHosts, int nloop) { + + off_t Nimage, i; + Image *images = getimages (&Nimage, NULL); + + int Nimage_mags = 0; + ImageMag *image_mags = NULL; + ALLOCATE (image_mags, ImageMag, 1); + + fprintf (stderr, "grabbing image mags from other hosts...\n"); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset; + char *imagfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits"); + ImageMag *image_mags_subset = ImageMagLoad (imagfile, &Nsubset); + free (imagfile); + + image_mags = merge_image_mags (image_mags, &Nimage_mags, image_mags_subset, Nsubset); + } + + for (i = 0; i < Nimage_mags; i++) { + off_t seq = getImageByID (image_mags[i].imageID); + if (seq < 0) { + // XXX is this a problem? (no, other hosts don't know which images I own) + continue; + } + images[seq].Mcal = image_mags[i].Mcal; + images[seq].dMcal = image_mags[i].dMcal; + images[seq].dMagSys = image_mags[i].dMagSys; + images[seq].Xm = image_mags[i].Xm; + images[seq].nFitPhotom = image_mags[i].nFitPhotom; + images[seq].flags = image_mags[i].flags; + images[seq].ubercalDist = image_mags[i].ubercalDist; + } + free (image_mags); + + fprintf (stderr, "DONE grabbing image mags from other hosts\n"); + + return TRUE; +} + +int set_image_mags (ImageMag *image_mags, Image *image) { + + image_mags->Mcal = image->Mcal; + image_mags->dMcal = image->dMcal; + image_mags->dMagSys = image->dMagSys; + image_mags->Xm = image->Xm; + image_mags->nFitPhotom = image->nFitPhotom; + image_mags->flags = image->flags; + image_mags->ubercalDist = image->ubercalDist; + image_mags->imageID = image->imageID; + + return TRUE; +} + +ImageMag *merge_image_mags (ImageMag *target, int *ntarget, ImageMag *source, int Nsource) { + + off_t i; + + REALLOCATE (target, ImageMag, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relphot/src/share_mean_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 36630) +++ trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 36630) @@ -0,0 +1,149 @@ +# include "relphot.h" +// we are sharing mean mags for all objects which (a) I own and (b) which have unowned detections + +# define D_NMEANMAGS 10000 +int share_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + int i; + off_t j; + + off_t Nmeanmags = 0; + off_t NMEANMAGS = D_NMEANMAGS; + + MeanMag *meanmags = NULL; + ALLOCATE (meanmags, MeanMag, NMEANMAGS); + + int Ns; + int Nsecfilt = GetPhotcodeNsecfilt(); + + int myHost = regionHosts->index[REGION_HOST_ID]; + double Rmin = regionHosts->hosts[myHost].Rmin; + double Rmax = regionHosts->hosts[myHost].Rmax; + double Dmin = regionHosts->hosts[myHost].Dmin; + double Dmax = regionHosts->hosts[myHost].Dmax; + + // XXX skip some catalogs based on UserPatch? + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + // do I own this object? (in region range?) + if (catalog[i].averageT[j].R < Rmin) continue; + if (catalog[i].averageT[j].R >= Rmax) continue; + if (catalog[i].averageT[j].D < Dmin) continue; + if (catalog[i].averageT[j].D >= Dmax) continue; + + // does this object have missing detections (does someone else need it?) + // XXX : sky objects without missing detections + // XXX watch out for detections which are not associated with an image (REF) + if (catalog[i].averageT[j].nOwn == catalog[i].averageT[j].Nmeasure) continue; + + for (Ns = 0; Ns < Nphotcodes; Ns++) { + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + set_mean_mags (&meanmags[Nmeanmags], &catalog[i].averageT[j], &catalog[i].secfilt[Nsecfilt*j + Nsec], Nsec); + Nmeanmags ++; + CHECK_REALLOCATE (meanmags, MeanMag, NMEANMAGS, Nmeanmags, D_NMEANMAGS); + } + } + } + + // write out the meanmag fits table AND write state in some file + char *hostname = regionHosts->hosts[myHost].hostname; + + char *magsfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanmags.fits"); + MeanMagSave (magsfile, meanmags, Nmeanmags); + free (meanmags); + free (magsfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "meanmags.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_mean_mags (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + off_t i; + + int Nmeanmags = 0; + MeanMag *meanmags = NULL; + ALLOCATE (meanmags, MeanMag, 1); + + fprintf (stderr, "grabbing mean object mags from other hosts...\n"); + + int Nsecfilt = GetPhotcodeNsecfilt(); + + for (i = 0; i < regionHosts->Nhosts; i++) { + // if (not_neighbor(host[i])) continue; + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset = 0; + char *magsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits"); + MeanMag *meanmagsSubset = MeanMagLoad (magsfile, &Nsubset); + free (magsfile); + + // merge_mean_mags reallocs meanmags and frees the input meanmagsSubset + meanmags = merge_mean_mags (meanmags, &Nmeanmags, meanmagsSubset, Nsubset); + } + + for (i = 0; i < Nmeanmags; i++) { + int objID = meanmags[i].objID; + int catID = meanmags[i].catID; + + // set the mean mag + int catSeq; + off_t objSeq; + if (!catID_and_objID_to_seq (catID, objID, &catSeq, &objSeq)) { + // XXX what should I do if this does not match? + continue; + } + + int ecode = GetPhotcodeEquivCodebyCode (meanmags[i].photcode); + if (ecode < 0) continue; + int Nsec = GetPhotcodeNsec (ecode); + if (Nsec < 0) continue; + + catalog[catSeq].secfilt[objSeq*Nsecfilt + Nsec].M = meanmags[i].M; + } + free (meanmags); + + fprintf (stderr, "DONE grabbing mean object mags from other hosts...\n"); + + return TRUE; +} + +int set_mean_mags (MeanMag *meanmags, AverageTiny *average, SecFilt *secfilt, int Nsec) { + + meanmags->M = secfilt->M; + meanmags->dM = secfilt->dM; + meanmags->Xm = secfilt->Xm; + meanmags->Nsec = Nsec; // key to secfilt entry + + meanmags->objID = average->objID; + meanmags->catID = average->catID; + + return TRUE; +} + +MeanMag *merge_mean_mags (MeanMag *target, int *ntarget, MeanMag *source, int Nsource) { + + off_t i; + + REALLOCATE (target, MeanMag, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relphot/src/syncfile.c =================================================================== --- trunk/Ohana/src/relphot/src/syncfile.c (revision 36630) +++ trunk/Ohana/src/relphot/src/syncfile.c (revision 36630) @@ -0,0 +1,77 @@ +# include "relphot.h" +# define MSG_LENGTH 12 + +char *make_filename (char *dirname, char *hostname, int hostID, char *tailname) { + + char tmp[10], *line; + int Nchar = snprintf (tmp, 0, "%s/%s.%03d.%s", dirname, hostname, hostID, tailname); + + ALLOCATE (line, char, Nchar + 1); + snprintf (line, Nchar + 1, "%s/%s.%03d.%s", dirname, hostname, hostID, tailname); + + return line; +} + +int check_sync_file (char *filename, int nloop) { + + char message[MSG_LENGTH]; + + FILE *f = NULL; + + while (TRUE) { + + f = fopen (filename, "r"); + if (!f) { + usleep (2000000); + continue; + } + + // XXX MSG_LENGTH : 0 EOL byte? + int Nread = fread (message, 1, MSG_LENGTH, f); + if (Nread < MSG_LENGTH) { + fclose (f); + usleep (2000000); + continue; + } + fclose (f); + + // message is of the form: NLOOP: %03d + int loop; + sscanf (message, "%*s %d", &loop); + if (loop != nloop) { + usleep (2000000); + continue; + } + return TRUE; + } + return FALSE; +} + +int clear_sync_file (char *filename) { + // delete file contents + truncate (filename, 0); + + return TRUE; +} + +int update_sync_file (char *filename, int nloop) { + + char message[MSG_LENGTH]; + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "failure to open sync file for write\n"); + exit (4); + } + + snprintf (message, MSG_LENGTH, "NLOOP: %03d\n", nloop); + + int Nwrite = fwrite (message, 1, MSG_LENGTH, f); + if (Nwrite != MSG_LENGTH) { + fprintf (stderr, "failure to write sync message\n"); + exit (3); + } + + fclose (f); + return TRUE; +} Index: trunk/Ohana/src/uniphot/src/update_dvo_setphot.c =================================================================== --- trunk/Ohana/src/uniphot/src/update_dvo_setphot.c (revision 36623) +++ trunk/Ohana/src/uniphot/src/update_dvo_setphot.c (revision 36630) @@ -2,4 +2,5 @@ int update_dvo_setphot_parallel (SkyTable *sky, Image *image, off_t Nimage); +// XXX convert this to use the ImageSubset to save space int update_dvo_setphot (Image *image, off_t Nimage, FlatCorrectionTable *flatcorr) { Index: trunk/psphot/src/psphotFullForceReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceReadout.c (revision 36623) +++ trunk/psphot/src/psphotFullForceReadout.c (revision 36630) @@ -100,4 +100,10 @@ psphotMagnitudes(config, view, filerule); + // calculate lensing parameters + if (!psphotLensing(config, view, filerule)) { + psErrorStackPrint(stderr, "Unable to do lensing parameters."); + psErrorClear(); + } + // replace background in residual image psphotSkyReplace (config, view, filerule); Index: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 36623) +++ trunk/psphot/src/psphotStackReadout.c (revision 36630) @@ -418,4 +418,10 @@ // calculate source magnitudes psphotMagnitudes(config, view, STACK_RAW); + + // calculate lensing parameters + if (!psphotLensing(config, view, STACK_RAW)) { + psErrorStackPrint(stderr, "Unable to do lensing parameters."); + psErrorClear(); + } if (!psphotEfficiency(config, view, STACK_RAW)) {