Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h (revision 37699) @@ -229,4 +229,5 @@ int RESET; +int RESET_IMAGES; int RESET_BAD_IMAGES; int NLOOP; Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FitChip.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FitChip.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FitChip.c (revision 37699) @@ -75,5 +75,10 @@ return FALSE; } - if (!CHIPMAP) order_use = MIN(order_use, 3); // can only go up to 3rd order + if (CHIPMAP) { + order_use = MIN(MIN(order_use, CHIPMAP), 6); // can only go up to 6th order map (can be user limited) + } else { + order_use = MIN(order_use, 3); // can only go up to 3rd order for polynomials + } + if (VERBOSE2) fprintf (stderr, "using %d for %s\n", order_use, image[0].name); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImageOps.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImageOps.c (revision 37699) @@ -125,6 +125,6 @@ if (imageSeq) { - myAssert (ID >= minImageID, "oops"); - myAssert (ID <= maxImageID, "oops"); + if (ID < minImageID) return (-1); + if (ID > maxImageID) return (-1); off_t N = imageSeq[ID]; return N; @@ -276,4 +276,6 @@ } } + + if (!ID) return; // detection not associated with an image idx = getImageByID (ID); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImagePosIO.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImagePosIO.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ImagePosIO.c (revision 37699) @@ -85,4 +85,5 @@ ALLOCATE (image_pos, ImagePos, Nrow); for (i = 0; i < Nrow; i++) { + InitCoords(&image_pos[i].coords, NULL); image_pos[i].coords.crval1 = crval1 [i]; image_pos[i].coords.crval2 = crval2 [i]; Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateChips.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateChips.c (revision 37699) @@ -564,5 +564,5 @@ AstromOffsetTableNewMap(table, tmpMap->Nx, tmpMap->Ny, image); // registers the map with the image unlockUpdateChips (); - AstromOffsetMapCopyData (truMap, tmpMap); + AstromOffsetMapCopyData (image->coords.offsetMap, tmpMap); return; } Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateMeasures.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateMeasures.c (revision 37699) @@ -1,9 +1,11 @@ # include "relastro.h" +void dump_measures (Average *average, Measure *measure); // in ImageOps.c -// this function operates on Measure, not MeasureTiny +// this function is called by UpdateObjectOffsets or relastro_images.c for reset int UpdateMeasures (Catalog *catalog, int Ncatalog) { - off_t i, Nimage; + off_t i, j, Nimage; Image *image; + double DPOS_MAX_ASEC; int badImage = @@ -18,20 +20,115 @@ image = getimages (&Nimage, NULL); + // track measurements which are far from their original position (poor at pole) + int NoffRAave = 0; int NoffRAori = 0; + int NoffDECave = 0; int NoffDECori = 0; + for (i = 0; i < Ncatalog; i++) { + myAssert (!catalog[i].Nmeasure || catalog[i].measureT, "programming error"); for (j = 0; j < catalog[i].Nmeasure; j++) { MeasureTiny *measureT = &catalog[i].measureT[j]; - Measure *measure = catalog[i].measure ? &catalog[i].measure[j] : NULL; + Measure *measureB = catalog[i].measure ? &catalog[i].measure[j] : NULL; off_t im = getImageByID (measureT->imageID); + if (im < 0) continue; // detections without imageIDs are not associated with a known image (eg, refs) - // check that we have the + // check that we have the right image: + myAssert (measureT->photcode == image[im].photcode, "image photcode mismatch"); + myAssert (measureT->t >= image[im].tzero, "image time mismatch (1)"); + myAssert (measureT->t <= image[im].tzero + image[im].NX*image[im].trate/5000.0, "image time mismatch (2)"); + // I'm allowing the trate*NX to be a factor of 2x too small - // skip images that have failed solutions (divergent or otherwise) - // XXX apply this or not?? - if (FALSE && (image[i].flags & badImage)) continue; + // skip measurements on images that have failed solutions (divergent or otherwise) + // XXX apply this or not?? + if (FALSE && (image[im].flags & badImage)) continue; + Coords *moscoords = image[im].coords.mosaic; + Coords *imcoords = &image[im].coords; + + if (moscoords) { + DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(moscoords->cdelt1, moscoords->cdelt2); + } else { + DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(imcoords->cdelt1, imcoords->cdelt2); + } + + double X = measureT->Xccd; + double Y = measureT->Yccd; + if (USE_FIXED_PIXCOORDS) { + if (isfinite(measureT->Xfix) && isfinite(measureT->Yfix)) { + float dX = measureT->Xfix - measureT->Xccd; + float dY = measureT->Yfix - measureT->Yccd; + if (hypot(dX,dY) < 2.0) { + X = measureT->Xfix; + Y = measureT->Yfix; + } + } + } + + double R, D; + double dR = NAN, dD = NAN; + XY_to_RD (&R, &D, X, Y, imcoords); + + int n = measureT[0].averef; + Average *average = &catalog[i].average[n]; + + // make sure detection is on the same side of the 0,360 boundary as the average + // this will give some funny results withing ~1 arcsec of the pol + if (dR > +180.0*3600.0) { + // average on high end of boundary, move star up + R += 360.0; + dR = 3600.0*(average[0].R - R); + } + if (dR < -180.0*3600.0) { + // average on low end of boundary, move star down + R -= 360.0; + dR = 3600.0*(average[0].R - R); + } + + float csdec = cos(average[0].D * RAD_DEG); + + // complain if the new location is far from the average location + // NOTE: This should never happen, or our StarMap tests are not working + if (fabs(dR*csdec) > 3.0*ADDSTAR_RADIUS) { + NoffRAave ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, dR*csdec, dD); + dump_measures (&average[0], catalog[i].measure); + } + } + if (fabs(dD) > 3.0*ADDSTAR_RADIUS) { + NoffDECave ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", average[0].R, average[0].D, dR, dD); + dump_measures (&average[0], catalog[i].measure); + } + } + + // complain if the new location is far from the old location + if (fabs(csdec*(measureT[0].R - R)) > DPOS_MAX_ASEC) { + NoffRAori ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement moves far from original location (R): %f %f (%f %f %f %f)\n", average[0].R, average[0].D, measureT[0].R, dR, csdec*(measureT[0].R - R), dD); + dump_measures (&average[0], catalog[i].measure); + } + } + if (fabs(measureT[0].D - D) > DPOS_MAX_ASEC) { + NoffDECori ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement moves far from original location (D): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, measureT[0].D, dD); + dump_measures (&average[0], catalog[i].measure); + } + } + measureT->R = R; + measureT->D = D; + if (measureB) { + measureB->R = R; + measureB->D = D; + } + } } + // printNcatTotal (); - // printNcatTotal (); + int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori; + if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori); return (TRUE); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 37699) @@ -80,4 +80,5 @@ int UpdateObjectOffsets_parallel_group (HostTableGroup *group, SkyList *sky); +int UpdateObjectOffsets_parallel_table (HostTable *table, SkyList *sky); // CATDIR is supplied globally @@ -94,4 +95,8 @@ } + + UpdateObjectOffsets_parallel_table (table, sky); + +# if (0) int Ngroups; HostTableGroup *groups = HostTableGroups (table, &Ngroups); @@ -103,4 +108,5 @@ UpdateObjectOffsets_parallel_group (&groups[i], sky); } +# endif return TRUE; @@ -224,2 +230,120 @@ return TRUE; } + +// CATDIR is supplied globally +# define DEBUG 1 +int UpdateObjectOffsets_parallel_table (HostTable *table, SkyList *sky) { + + // launch the relastro_client jobs to the parallel hosts + + int i, j; + for (i = 0; i < table->Nhosts; i++) { + + if (sky->Nregions < table->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + // options / arguments that can affect relastro_client -load: + // VERBOSE, VERBOSE2 + // RESET (-reset) + // TimeSelect -time + // (note that psfQF is applied rigidly at 0.85, as is the galaxy test) + // ImagSelect, ImagMin, ImagMax + // MaxDensityUse, MaxDensityValue + + char command[1024]; + snprintf (command, 1024, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR); + + if (FIT_MODE == FIT_PM_ONLY) strextend (command, "-pm"); + if (FIT_MODE == FIT_PAR_ONLY) strextend (command, "-par"); + if (FIT_MODE == FIT_PM_AND_PAR) strextend (command, "-pmpar"); + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (UPDATE) strextend (command, "-update"); + + if (RESET_BAD_IMAGES) strextend (command, "-reset-bad-images"); + + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + + if (USE_BASIC_CHECK) strextend (command, "-basic-image-search"); + if (FlagOutlier) strextend (command, "-clip %d", CLIP_THRESH); + if (ExcludeBogus) strextend (command, "-exclude-bogus %f", ExcludeBogusRadius); + + if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); + + if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + if (PhotFlagSelect) strextend (command, "+photflags"); + if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); + // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + + if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { + strextend (command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); + } + if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { + strextend (command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); + } + + 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_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running relastro_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the relastro_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/args.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/args.c (revision 37699) @@ -372,8 +372,10 @@ } - CHIPMAP = FALSE; + CHIPMAP = 0; if ((N = get_argument (argc, argv, "-chipmap"))) { remove_argument (N, &argc, argv); - CHIPMAP = TRUE; + CHIPMAP = atoi(argv[N]); + remove_argument (N, &argc, argv); + } @@ -405,9 +407,13 @@ RESET = TRUE; } - RESET_BAD_IMAGES = FALSE; if ((N = get_argument (argc, argv, "-reset-bad-images"))) { remove_argument (N, &argc, argv); RESET_BAD_IMAGES = TRUE; + } + RESET_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-reset-images"))) { + remove_argument (N, &argc, argv); + RESET_IMAGES = TRUE; } @@ -874,4 +880,5 @@ fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n"); fprintf (stderr, " -hpm (radius) (output catdir)\n"); + fprintf (stderr, " -chipmap (MaxOrder)\n"); fprintf (stderr, " -merge-source (objID) (catID) into (objID) (catID)\n\n"); @@ -886,4 +893,5 @@ fprintf (stderr, " -statmode (mode)\n"); fprintf (stderr, " -reset\n"); + fprintf (stderr, " -reset-images\n"); fprintf (stderr, " -nloop (N) : number of image-fit iterations\n"); fprintf (stderr, " -update : apply new fit to database\n"); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c (revision 37699) @@ -128,4 +128,7 @@ if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); + if (CHIPMAP) strextend (command, "-chipmap %d", CHIPMAP); + if (RESET_IMAGES) strextend (command, "-reset-images"); + if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c (revision 37699) @@ -70,4 +70,5 @@ case TARGET_CHIPS: + if (RESET_IMAGES) UpdateMeasures (catalog, Ncatalog); for (i = 0; i < NLOOP; i++) { UpdateObjects (catalog, Ncatalog, (i > 0)); // calculate ,; if (i > 0), apply Galaxy Motion Model (if desired) Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c (revision 37699) @@ -103,4 +103,5 @@ case TARGET_CHIPS: for (i = 0; i < NLOOP; i++) { + if (RESET_IMAGES) UpdateMeasures (catalog, Ncatalog); UpdateObjects (catalog, Ncatalog, (i > 0)); if ((i > 1) || !USE_GALAXY_MODEL) { Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 37699) @@ -63,5 +63,7 @@ if (PARALLEL) { - // save the updated image parameters + // save the updated image parameters. if we are NOT PARALLEL, we have to wait + // until after UpdateObjectOffsets, or the image parameters are swapped and/or freed + save_astrom_table (); dvo_image_save (&db, VERBOSE); dvo_image_unlock (&db); @@ -73,4 +75,5 @@ if (!PARALLEL) { // save the changes to the image parameters + save_astrom_table (); dvo_image_save (&db, VERBOSE); dvo_image_unlock (&db); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_images_pos.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_images_pos.c (revision 37698) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_images_pos.c (revision 37699) @@ -84,5 +84,8 @@ continue; } + Coords *moscoords = images[seq].coords.mosaic; images[seq].coords = image_pos[i].coords ; + images[seq].coords.mosaic = moscoords; + images[seq].dXpixSys = image_pos[i].dXpixSys ; images[seq].dYpixSys = image_pos[i].dYpixSys ; @@ -108,7 +111,19 @@ if (seq < 0) continue; - images[seq].coords.offsetMap = table->map[j]; - // is this sufficient? + AstromOffsetMap *oldMap = images[seq].coords.offsetMap; + AstromOffsetMap *newMap = table->map[j]; + + if (oldMap) { + AstromOffsetMapSetOrder (oldMap, newMap->Nx, newMap->Ny, &images[seq]); + AstromOffsetMapCopyData (oldMap, newMap); + } else { + lockUpdateChips (); + AstromOffsetTable *FullTable = get_astrom_table (); + AstromOffsetTableNewMap(FullTable, newMap->Nx, newMap->Ny, &images[seq]); // registers the map with the image + unlockUpdateChips (); + AstromOffsetMapCopyData (images[seq].coords.offsetMap, newMap); // oldMap is now not valid + } } + AstromOffsetTableFree(table); }