Index: trunk/Ohana/src/addstar/src/ReadImageHeader.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 41801) +++ trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 41891) @@ -70,12 +70,15 @@ } + // 2021.11.03 : consider replacing or supplementing CERROR with CERSTD + // if (!gfits_scan (header, "CERROR", "%lf", 1, &tmp)) tmp = 1.0; + /* require Nastro > 0 unless or ACCEPT_ASTROM */ - Nastro = 0; - gfits_scan (header, "NASTRO", "%d", 1, &Nastro); - Cerror = 0; + FWHM_X = FWHM_Y = Nastro = Cerror = 0; + + gfits_scan (header, "NASTRO", "%d", 1, &Nastro); gfits_scan (header, "CERROR", "%lf", 1, &Cerror); - FWHM_X = FWHM_Y = 0; gfits_scan (header, "FWHM_X", "%lf", 1, &FWHM_X); gfits_scan (header, "FWHM_Y", "%lf", 1, &FWHM_Y); + if (((Nastro == 0) || (Cerror > MAX_CERROR)) && !ACCEPT_ASTROM) { fprintf (stderr, "bad astrometric solution in header\n"); @@ -111,7 +114,6 @@ } - /* CERROR in data file is in arcsec */ - if (!gfits_scan (header, "CERROR", "%lf", 1, &tmp)) tmp = 1.0; - image[0].cerror = tmp * 50.0; + /* CERROR / CERSTD in data file in arcsec; db field is units of 20 mas */ + image[0].cerror = 50.0 * Cerror; /* get photcode from header */ Index: trunk/Ohana/src/getstar/src/ConfigInit_overlaps.c =================================================================== --- trunk/Ohana/src/getstar/src/ConfigInit_overlaps.c (revision 41801) +++ trunk/Ohana/src/getstar/src/ConfigInit_overlaps.c (revision 41891) @@ -44,5 +44,5 @@ if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); - MAX_CERROR = 0.5; // arcseconds + MAX_CERROR = NAN; // NAN means do not cut on this value ScanConfig (config, "OVERLAPS_MAX_CERROR", "%lf", 0, &MAX_CERROR); Index: trunk/Ohana/src/getstar/src/MatchImages.c =================================================================== --- trunk/Ohana/src/getstar/src/MatchImages.c (revision 41801) +++ trunk/Ohana/src/getstar/src/MatchImages.c (revision 41891) @@ -21,5 +21,6 @@ // exclude chips / images with errors larger than OVERLAP_MAX_CERROR: - if (fabs(MAX_CERROR) > 1e-6) { + // 2021.11.03 : MAX_CERROR is now compared to the header value CERSTD + if (isfinite(MAX_CERROR)) { if (image[0].cerror * 0.02 > MAX_CERROR) return NULL; } Index: trunk/Ohana/src/getstar/src/ReadImageHeader.c =================================================================== --- trunk/Ohana/src/getstar/src/ReadImageHeader.c (revision 41801) +++ trunk/Ohana/src/getstar/src/ReadImageHeader.c (revision 41891) @@ -93,6 +93,8 @@ } - /* CERROR in data file is in arcsec, image structure uses units of 20 mas */ - if (!gfits_scan (header, "CERROR", "%lf", 1, &tmp)) tmp = 1.0; + // CERROR & CERSTD in data file are in arcsec, image structure uses units of 20 mas + // if (!gfits_scan (header, "CERROR", "%lf", 1, &tmp)) tmp = 1.0; + // 2021.11.03 : getstar now checks CERSTD, not CERROR + if (!gfits_scan (header, "CERSTD", "%lf", 1, &tmp)) tmp = 1.0; image[0].cerror = tmp * 50.0; Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 41891) @@ -28,4 +28,6 @@ $(SRC)/bisection.$(ARCH).o \ $(SRC)/cast.$(ARCH).o \ +$(SRC)/chebyshev.$(ARCH).o \ +$(SRC)/chebyshev_commands.$(ARCH).o \ $(SRC)/center.$(ARCH).o \ $(SRC)/clear.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.data/applyfit2d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/applyfit2d.c (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/applyfit2d.c (revision 41891) @@ -116,2 +116,117 @@ +int applyfit2d_full (int argc, char **argv) { + + int i, j, n, order; + char *c, name[64]; + double **C, X, Y; + opihi_flt *z; + Vector *xvec, *yvec, *zvec; + + if (argc != 4) { + gprint (GP_ERR, "USAGE: applyfit2d_full x y z\n"); + return (FALSE); + } + + c = get_variable ("Cnn"); + if (c == NULL) { + gprint (GP_ERR, "no fit available\n"); + return (FALSE); + } + order = atof (c); + free (c); + + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((zvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + ALLOCATE (C, double *, order+1); + for (i = 0; i <= order; i++) { + ALLOCATE (C[i], double, order+1); + for (j = 0; j <= order; j++) { + sprintf (name, "CX%dY%d", i, j); + c = get_variable (name); + if (c == NULL) { + gprint (GP_ERR, "missing fit term %d,%d\n", i, j); + for (j = 0; j < i; j++) free (C[j]); + free (C); + return (FALSE); + } + C[i][j] = atof (c); + free (c); + } + } + + ResetVector (zvec, OPIHI_FLT, xvec[0].Nelements); + bzero (zvec[0].elements.Flt, sizeof(opihi_flt)*zvec[0].Nelements); + z = zvec[0].elements.Flt; + + // we have four cases (x == flt or int) and (y == flt or int) + if ((xvec[0].type == OPIHI_FLT) && (yvec[0].type == OPIHI_FLT)) { + opihi_flt *x = xvec[0].elements.Flt; + opihi_flt *y = yvec[0].elements.Flt; + for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) { + Y = X = 1; + for (j = 0; j <= order; j++) { + X = Y; + for (i = 0; i <= order; i++) { + *z += C[i][j]*X; + X = X * (*x); + } + Y = Y * (*y); + } + } + } + if ((xvec[0].type == OPIHI_FLT) && (yvec[0].type == OPIHI_INT)) { + opihi_flt *x = xvec[0].elements.Flt; + opihi_int *y = yvec[0].elements.Int; + for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) { + Y = X = 1; + for (j = 0; j <= order; j++) { + X = Y; + for (i = 0; i <= order; i++) { + *z += C[i][j]*X; + X = X * (*x); + } + Y = Y * (*y); + } + } + } + if ((xvec[0].type == OPIHI_INT) && (yvec[0].type == OPIHI_FLT)) { + opihi_int *x = xvec[0].elements.Int; + opihi_flt *y = yvec[0].elements.Flt; + for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) { + Y = X = 1; + for (j = 0; j <= order; j++) { + X = Y; + for (i = 0; i <= order; i++) { + *z += C[i][j]*X; + X = X * (*x); + } + Y = Y * (*y); + } + } + } + if ((xvec[0].type == OPIHI_INT) && (yvec[0].type == OPIHI_INT)) { + opihi_int *x = xvec[0].elements.Int; + opihi_int *y = yvec[0].elements.Int; + for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++) { + Y = X = 1; + for (j = 0; j <= order; j++) { + X = Y; + for (i = 0; i <= order; i++) { + *z += C[i][j]*X; + X = X * (*x); + } + Y = Y * (*y); + } + } + } + + for (i = 0; i < order + 1; i++) free (C[i]); + free (C); + + return (TRUE); +} + + Index: trunk/Ohana/src/opihi/cmd.data/chebyshev.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/chebyshev.c (revision 41891) +++ trunk/Ohana/src/opihi/cmd.data/chebyshev.c (revision 41891) @@ -0,0 +1,167 @@ +# include "data.h" + +int chebyshev_help (int argc, char **argv); +int chebyshev_list (int argc, char **argv); +int chebyshev_delete (int argc, char **argv); +int chebyshev_rename (int argc, char **argv); + +int chebyshev_poly (int argc, char **argv); +int chebyshev_scale (int argc, char **argv); +int chebyshev_applyscale (int argc, char **argv); + +int chebyshev_fit1d (int argc, char **argv); +int chebyshev_applyfit1d (int argc, char **argv); + +int chebyshev_fit2d (int argc, char **argv); +int chebyshev_applyfit2d (int argc, char **argv); + +// collection of chebyshev commands: +static Command chebyshev_commands[] = { + {1, "help", chebyshev_help, "chebyshev help info"}, + {1, "list", chebyshev_list, "list define chebyshev polynomial sets"}, + {1, "delete", chebyshev_delete, "delete a chebyshev"}, + {1, "rename", chebyshev_rename, "rename a chebyshev"}, + + {1, "poly", chebyshev_poly, "return evaluated chebyshev polynomial of given order for supplied coordinate vector"}, + + {1, "scale", chebyshev_scale, "set domain scale based on supplied coordinate vector"}, + {1, "applyscale", chebyshev_applyscale, "apply domain scale to supplied coordinate vector"}, + + {1, "fit1d", chebyshev_fit1d, "fit chebyshev polynomial to supplied data vector"}, + {1, "applyfit1d", chebyshev_applyfit1d, "apply chebyshev polynomial fit to supplied coordinate vector"}, + + {1, "fit2d", chebyshev_fit2d, "fit chebyshev polynomial to supplied data vector"}, + {1, "applyfit2d", chebyshev_applyfit2d, "apply chebyshev polynomial fit to supplied coordinate vector"}, +}; + +int chebyshev_command (int argc, char **argv) { + + int N; + + if (argc < 2) { + chebyshev_help(0,NULL); + return (FALSE); + } + + N = sizeof (chebyshev_commands) / sizeof (Command); + + /* find the chebyshev sub-command which matches */ + for (int i = 0; i < N; i++) { + if (!strcmp (chebyshev_commands[i].name, argv[1])) { + int status = (*chebyshev_commands[i].func) (argc - 1, argv + 1); + return (status); + } + } + + gprint (GP_ERR, "unknown chebyshev command %s\n", argv[1]); + return (FALSE); +} + +/*** below are the utility commands which just manage the collection of chebyshevs ***/ + +int chebyshev_help (int argc, char **argv) { + OHANA_UNUSED_PARAM(argc); + OHANA_UNUSED_PARAM(argv); + + gprint (GP_ERR, "USAGE: chebyshev (command)\n"); + gprint (GP_ERR, " chebyshev help : this listing\n"); + gprint (GP_ERR, " chebyshev list : list chebyshevs\n"); + gprint (GP_ERR, " chebyshev delete (chebyshev) : delete named chebyshev\n"); + gprint (GP_ERR, " chebyshev rename (chebyshev) (new) : change chebyshev name to new name\n"); + + gprint (GP_ERR, " chebyshev poly (chebyshev) (new) : change chebyshev name to new name\n"); + gprint (GP_ERR, " chebyshev scale (chebyshev) (new) : change chebyshev name to new name\n"); + gprint (GP_ERR, " chebyshev applyscale (chebyshev) (new) : change chebyshev name to new name\n"); + + gprint (GP_ERR, " chebyshev fit1d (chebyshev) (new) : change chebyshev name to new name\n"); + gprint (GP_ERR, " chebyshev applyfit1d (chebyshev) (new) : change chebyshev name to new name\n"); + + gprint (GP_ERR, " chebyshev fit2d (chebyshev) (new) : change chebyshev name to new name\n"); + gprint (GP_ERR, " chebyshev applyfit2d (chebyshev) (new) : change chebyshev name to new name\n"); + + return FALSE; +} + +int chebyshev_list (int argc, char **argv) { + OHANA_UNUSED_PARAM(argv); + + if (argc != 1) { + gprint (GP_ERR, "USAGE: chebyshev list\n"); + return FALSE; + } + + ListChebyshevs(); + return TRUE; +} + +int chebyshev_delete (int argc, char **argv) { + + int N; + + int QUIET = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + QUIET = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 2) { + gprint (GP_ERR, "USAGE: chebyshev delete [-q] (chebyshev)\n"); + return FALSE; + } + + ChebyshevType *chebyshev = FindChebyshev (argv[1]); + if (chebyshev == NULL) { + if (QUIET) return TRUE; + gprint (GP_ERR, "chebyshev %s not found\n", argv[1]); + return FALSE; + } + + int status = DeleteChebyshev (chebyshev); + if (!status) abort (); // status = FALSE means the cheb selected above was not found + return TRUE; +} + +int chebyshev_rename (int argc, char **argv) { + + ChebyshevType *chebyshev; + + if (argc != 3) { + gprint (GP_ERR, "USAGE: chebyshev rename (chebyshev) (newname)\n"); + return FALSE; + } + + chebyshev = FindChebyshev (argv[1]); + if (chebyshev == NULL) { + gprint (GP_ERR, "chebyshev %s not found\n", argv[1]); + return FALSE; + } + + free (chebyshev->name); + chebyshev->name = strcreate (argv[2]); + return TRUE; +} + +/* + + typedef struct { + char *name; + opihi_flt scale; + opihi_flt zero; + int order; + opihi_flt *A; + } ChebyshevType; + + + examples of using the + + cheb fit1d t1 x y 4 -- fit the vector y = f(x) using 4th order + (saves the scaling for x, order, coeffs) + cheb applyfit1d t1 x yfit -- apply the fit defined in t1 to x + cheb applyfit1d t1 X Yfit -- apply the fit defined in t1 to a different vector X + + cheb fit t2 x y 2 + + + + */ + Index: trunk/Ohana/src/opihi/cmd.data/chebyshev_commands.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/chebyshev_commands.c (revision 41891) +++ trunk/Ohana/src/opihi/cmd.data/chebyshev_commands.c (revision 41891) @@ -0,0 +1,344 @@ +# include "data.h" + +int chebyshev_poly (int argc, char **argv) { + + if (argc != 4) { + gprint (GP_ERR, "USAGE: chebyshev poly (xvec) (yvec) (order)\n"); + gprint (GP_ERR, " return requested chebyshev polynomial for given order\n"); + gprint (GP_ERR, " xvec should be normalized (range of -1 to +1)\n"); + return FALSE; + } + + Vector *xvec = NULL, *yvec = NULL; + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + int order = atoi(argv[3]); + + ResetVector (yvec, OPIHI_FLT, xvec->Nelements); + + ChebyshevPolyVector (yvec, xvec, order); + + return TRUE; +} + +int chebyshev_fit1d (int argc, char **argv) { + + int N; + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 5) { + gprint (GP_ERR, "USAGE: chebyshev fit1d (name) (xv) (yv) (order)\n"); + gprint (GP_ERR, " 1D fit of xv to yv using given order\n"); + return FALSE; + } + + ChebyshevType *cheb = FindChebyshev (argv[1]); + if (!cheb) { + if (VERBOSE) gprint (GP_ERR, "chebyshev %s not found, creating it\n", argv[1]); + cheb = CreateChebyshev (argv[1]); + } + + // select the x,y vectors + Vector *xvec = NULL, *yvec = NULL; + if ((xvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + cheb->dimen = 1; + cheb->order = atoi (argv[4]); + int nterm = cheb->order + 1; + + // check type and size: + if (xvec->Nelements < nterm) { + gprint (GP_ERR, "insufficient data to support the requested order\n"); + return FALSE; + } + if (xvec->Nelements != yvec->Nelements) { + gprint (GP_ERR, "vector lengths %s (%d) and %s (%d) do not match\n", xvec->name, xvec->Nelements, yvec->name, yvec->Nelements); + return FALSE; + } + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + + ChebyshevSetScale (cheb, xvec, 0); + + // generate the normalized vector xNorm + // use a local vector? + Vector *xNorm = ChebyshevNormVector (cheb, xvec, 0); + + // generate the N cheb polynomials based on xNorm + ALLOCATE_PTR (poly, Vector *, nterm); + for (int i = 0; i < nterm; i++) { + poly[i] = ChebyshevPolyVector (NULL, xNorm, i); + } + + // fit the data to the cheb polynomials + ChebyshevPolyFit1D (cheb, yvec, poly); + + for (int i = 0; i < nterm; i++) { + fprintf (stderr, "%d : %e\n", i, cheb->A[i]); + } + + // free local temporary variables + for (int i = 0; i < nterm; i++) { + FREE (poly[i]); + } + FREE (poly); + FREE (xNorm); + + return TRUE; +} + +int chebyshev_applyfit1d (int argc, char **argv) { + + if (argc != 4) { + gprint (GP_ERR, "USAGE: chebyshev applyfit1d (name) (xvec) (yfit)\n"); + gprint (GP_ERR, " apply 1D fit to xvec position to yield yfit\n"); + return FALSE; + } + + ChebyshevType *cheb = FindChebyshev (argv[1]); + if (!cheb) { + gprint (GP_ERR, "chebyshev %s not found\n", argv[1]); + return FALSE; + } + int nterm = cheb->order + 1; + + Vector *xvec = NULL, *yvec = NULL; + if ((xvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + ResetVector (yvec, OPIHI_FLT, xvec->Nelements); + + // generate the normalized vector xNorm + Vector *xNorm = ChebyshevNormVector (cheb, xvec, 0); + + // generate the N cheb polynomials based on xNorm + ALLOCATE_PTR (poly, Vector *, nterm); + for (int i = 0; i < nterm; i++) { + poly[i] = ChebyshevPolyVector (NULL, xNorm, i); + } + + ChebyshevPolyApplyFit1D (cheb, yvec, poly); + + // free local temporary variables + for (int i = 0; i < nterm; i++) { + FREE (poly[i]); + } + FREE (poly); + FREE (xNorm); + + return TRUE; +} + +int chebyshev_fit2d (int argc, char **argv) { + + int N; + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 6) { + gprint (GP_ERR, "USAGE: chebyshev fit1d (name) (xv) (yv) (zv) (order)\n"); + gprint (GP_ERR, " 2D fit of (xv,yv) to zv using given order\n"); + return FALSE; + } + + ChebyshevType *cheb = FindChebyshev (argv[1]); + if (!cheb) { + if (VERBOSE) gprint (GP_ERR, "chebyshev %s not found, creating it\n", argv[1]); + cheb = CreateChebyshev (argv[1]); + } + + // select the x,y vectors + Vector *xvec = NULL, *yvec = NULL, *zvec = NULL; + if ((xvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((zvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + cheb->dimen = 2; + cheb->order = atoi (argv[5]); + int nterm = cheb->order + 1; + + // check type and size: + if (xvec->Nelements < nterm) { + gprint (GP_ERR, "insufficient data to support the requested order\n"); + return FALSE; + } + if (xvec->Nelements != yvec->Nelements) { + gprint (GP_ERR, "vector lengths %s (%d) and %s (%d) do not match\n", xvec->name, xvec->Nelements, yvec->name, yvec->Nelements); + return FALSE; + } + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + CastVector (zvec, OPIHI_FLT); + + ChebyshevSetScale (cheb, xvec, 0); + ChebyshevSetScale (cheb, yvec, 1); + + // generate the normalized vector xNorm + // use a local vector? + Vector *xNorm = ChebyshevNormVector (cheb, xvec, 0); + Vector *yNorm = ChebyshevNormVector (cheb, yvec, 1); + + // generate the N cheb polynomials based on xNorm, yNorm + ALLOCATE_PTR (xPoly, Vector *, nterm); + ALLOCATE_PTR (yPoly, Vector *, nterm); + for (int i = 0; i < nterm; i++) { + xPoly[i] = ChebyshevPolyVector (NULL, xNorm, i); + yPoly[i] = ChebyshevPolyVector (NULL, yNorm, i); + } + + // fit the data to the cheb polynomials + ChebyshevPolyFit2D (cheb, zvec, xPoly, yPoly); + + int n = 0; + for (int ix = 0; ix < nterm; ix++) { + for (int iy = 0; iy < nterm; iy++) { + fprintf (stderr, "%d %d : %e\n", ix, iy, cheb->A[n]); + n ++; + } + } + + // free local temporary variables + for (int i = 0; i < nterm; i++) { + FREE (xPoly[i]); + FREE (yPoly[i]); + } + FREE (xPoly); + FREE (yPoly); + FREE (xNorm); + FREE (yNorm); + + return TRUE; +} + +int chebyshev_applyfit2d (int argc, char **argv) { + + if (argc != 5) { + gprint (GP_ERR, "USAGE: chebyshev applyfit1d (name) (xvec) (yvec) (xfit)\n"); + gprint (GP_ERR, " apply 1D fit to xvec,yvec positions to yield zfit\n"); + return FALSE; + } + + ChebyshevType *cheb = FindChebyshev (argv[1]); + if (!cheb) { + gprint (GP_ERR, "chebyshev %s not found\n", argv[1]); + return FALSE; + } + int nterm = cheb->order + 1; + + Vector *xvec = NULL, *yvec = NULL, *zvec = NULL; + if ((xvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((zvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + if (xvec->Nelements != yvec->Nelements) { + gprint (GP_ERR, "vector lengths %s (%d) and %s (%d) do not match\n", xvec->name, xvec->Nelements, yvec->name, yvec->Nelements); + return FALSE; + } + + ResetVector (zvec, OPIHI_FLT, xvec->Nelements); + + // generate the normalized vector xNorm + // use a local vector? + Vector *xNorm = ChebyshevNormVector (cheb, xvec, 0); + Vector *yNorm = ChebyshevNormVector (cheb, yvec, 1); + + // generate the N cheb polynomials based on xNorm, yNorm + ALLOCATE_PTR (xPoly, Vector *, nterm); + ALLOCATE_PTR (yPoly, Vector *, nterm); + for (int i = 0; i < nterm; i++) { + xPoly[i] = ChebyshevPolyVector (NULL, xNorm, i); + yPoly[i] = ChebyshevPolyVector (NULL, yNorm, i); + } + + + ChebyshevPolyApplyFit2D (cheb, zvec, xPoly, yPoly); + + // free local temporary variables + for (int i = 0; i < nterm; i++) { + FREE (xPoly[i]); + FREE (yPoly[i]); + } + FREE (xPoly); + FREE (yPoly); + FREE (xNorm); + FREE (yNorm); + + return TRUE; +} + +int chebyshev_scale (int argc, char **argv) { + + if (argc != 3) { + gprint (GP_ERR, "USAGE: chebyshev scale (name) (vector)\n"); + gprint (GP_ERR, " determine scale factors for the names chebyshev based on the vector\n"); + return FALSE; + } + + ChebyshevType *cheb = FindChebyshev (argv[1]); + if (!cheb) { + gprint (GP_ERR, "chebyshev %s not found, creating it\n", argv[1]); + cheb = CreateChebyshev (argv[1]); + } else { + gprint (GP_ERR, "resetting scale factors for %s\n", argv[1]); + } + + Vector *vec = NULL; + if ((vec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + // check type and size: + if (vec->type != OPIHI_FLT) { + gprint (GP_ERR, "cannot generate chebyshev scale factors for non-float vector\n"); + return FALSE; + } + if (!vec->Nelements) { + gprint (GP_ERR, "cannot generate chebyshev scale factors based on empty vector\n"); + return FALSE; + } + + int dir = 0; + ChebyshevSetScale (cheb, vec, dir); + + return TRUE; +} + +int chebyshev_applyscale (int argc, char **argv) { + + if (argc != 4) { + gprint (GP_ERR, "USAGE: chebyshev applyscale (name) (vector) (output)\n"); + gprint (GP_ERR, " re-scale vector for chebyshev\n"); + return FALSE; + } + + ChebyshevType *cheb = FindChebyshev (argv[1]); + if (!cheb) { + gprint (GP_ERR, "chebyshev %s not found\n", argv[1]); + return FALSE; + } + + Vector *vec = NULL, *out = NULL; + if ((vec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((out = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE); + CastVector (vec, OPIHI_FLT); + + ResetVector (out, OPIHI_FLT, vec->Nelements); + + int dir = 0; + + for (int i = 0; i < vec->Nelements; i++) { + if (!isfinite(vec->elements.Flt[i])) { out->elements.Flt[i] = NAN; continue; } + out->elements.Flt[i] = cheb->scale[dir]*vec->elements.Flt[i] + cheb->zero[dir]; + } + + return TRUE; +} + Index: trunk/Ohana/src/opihi/cmd.data/fit2d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/fit2d.c (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/fit2d.c (revision 41891) @@ -294,2 +294,272 @@ return (TRUE); } + +int fit2d_full (int argc, char **argv) { + + int N; + char name[64]; + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + int Quiet = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + Quiet = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-quiet"))) { + Quiet = TRUE; + remove_argument (N, &argc, argv); + } + + float ClipNSigma = 0; + int ClipNiter = 1; + if ((N = get_argument (argc, argv, "-clip"))) { + remove_argument (N, &argc, argv); + ClipNSigma = atof(argv[N]); + remove_argument (N, &argc, argv); + ClipNiter = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + opihi_flt *dz = NULL; + Vector *dzvec = NULL; + int Weight = FALSE; + if ((N = get_argument (argc, argv, "-dz"))) { + remove_argument (N, &argc, argv); + if ((dzvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + Weight = TRUE; + } + + if (argc != 5) { + gprint (GP_ERR, "USAGE: fit2d_full x y z order [-dz wt] [-quiet/-q] [-clip Nsigma Niter]\n"); + return (FALSE); + } + + Vector *xvec, *yvec, *zvec; + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((zvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + if (xvec[0].Nelements != yvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + if (xvec[0].Nelements != zvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + CastVector (zvec, OPIHI_FLT); + + if (Weight) { + CastVector (dzvec, OPIHI_FLT); + if (xvec[0].Nelements != dzvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + } + + int order = atof (argv[4]); + int nterm = SQ(order + 1); + + ALLOCATE_PTR (zfit, opihi_flt, xvec[0].Nelements); + ALLOCATE_PTR (mask, char, xvec[0].Nelements); + memset (mask, 0, xvec[0].Nelements); + + /* allocate the summation matrices */ + ALLOCATE_PTR (b, double *, nterm); + ALLOCATE_PTR (c, double *, nterm); + for (int i = 0; i < nterm; i++) { + ALLOCATE (c[i], double, nterm); + ALLOCATE (b[i], double, 1); + } + + for (N = 0; N < ClipNiter; N++) { + + /* init registers for current pass */ + for (int i = 0; i < nterm; i++) { + memset (c[i], 0, nterm*sizeof(double)); + memset (b[i], 0, sizeof(double)); + } + + opihi_flt *x = xvec[0].elements.Flt; + opihi_flt *y = yvec[0].elements.Flt; + opihi_flt *z = zvec[0].elements.Flt; + if (Weight) dz = dzvec[0].elements.Flt; + + /* add up the x,y values */ + for (int i = 0; i < xvec[0].Nelements; i++, x++, y++, z++) { + if (mask[i]) continue; + if (!finite(*x) || !finite(*y)) continue; + double dZ = 1.0; + if (Weight) { + dZ = 1.0 / SQ(*dz); + dz ++; + } + + int n = 0; + double iY = 1*dZ; + double iX = iY; + + for (int iy = 0; iy <= order; iy++) { + iX = iY; + for (int ix = 0; ix <= order; ix++) { + b[n][0] += iX * (*z); + + int m = 0; + double jY = 1; + double jX = jY; + + for (int jy = 0; jy <= order; jy++) { + jX = jY; + for (int jx = 0; jx <= order; jx++) { + c[n][m] += iX*jX; + jX = jX * (*x); + m ++; + } + jY = jY * (*y); + } + iX = iX * (*x); + n ++; + } + iY = iY * (*y); + } + } + + /** test print **/ + if (VERBOSE) { + for (int i = 0; i < nterm; i++) { + for (int j = 0; j < nterm; j++) { + gprint (GP_ERR, "%g ", c[i][j]); + } + gprint (GP_ERR, "\n"); + } + gprint (GP_ERR, "-----\n"); + } + + dgaussjordan (c, b, nterm, 1); + + /** test print **/ + if (VERBOSE) { + int n = 0; + for (int ny = 0; ny <= order; ny++) { + for (int nx = 0; nx <= order; nx++) { + gprint (GP_ERR, "x^%d y^%d: %g\n", nx, ny, b[n][0]); + n ++; + } + } + } + + /** test print **/ + if (VERBOSE) { + for (int i = 0; i < nterm; i++) { + for (int j = 0; j < nterm; j++) { + gprint (GP_ERR, "%g ", c[i][j]); + } + gprint (GP_ERR, "\n"); + } + gprint (GP_ERR, "-----\n"); + } + + /* the b[][0] terms are in the following order: + y^0 x^0, y^0 x^1, ... y^0 x^N + y^1 x^0, y^1 x^1, ... y^1 x^N + ... + y^N x^0, y^N x^1, ... y^N x^N + */ + + /* generate fitted values */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + opihi_flt *zf = zfit; + for (int n = 0; n < xvec[0].Nelements; n++, x++, y++, zf++) { + if (!finite(*x) || !finite(*y)) continue; + *zf = 0; + double Y = 1; + double X = Y; + + int n = 0; + for (int ny = 0; ny <= order; ny++) { + X = Y; + for (int nx = 0; nx <= order; nx++) { + *zf += b[n][0]*X; + X = X * (*x); + n ++; + } + Y = Y * (*y); + } + } + + /* measure fit residual scatter */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + zf = zfit; + double dZ = 0.0; + double dZ2 = 0.0; + int Npt = 0; + for (int i = 0, Npt = 0; i < xvec[0].Nelements; i++, x++, y++, z++, zf++) { + if (mask[i]) continue; + if (!finite(*x) || !finite(*y)) continue; + dZ += (*z - *zf); + dZ2 += SQ(*z - *zf); + Npt ++; + } + double mean = dZ / Npt; + double sigma = sqrt (fabs(dZ2/Npt - SQ(mean))); + double maxsigma = ClipNSigma * sigma; + + if (VERBOSE) gprint (GP_ERR, "mean: %g, sigma: %g, maxsigma: %g\n", mean, sigma, maxsigma); + + /* mask outlier points */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + zf = zfit; + int Nmask = 0; + for (int i = 0; ClipNSigma && (i < xvec[0].Nelements); i++, x++, y++, z++, zf++) { + dZ = (*z - *zf); + if (fabs(dZ) > maxsigma) { + mask[i] = TRUE; + Nmask ++; + } else { + mask[i] = FALSE; + } + } + if (VERBOSE) gprint (GP_ERR, "pass: %d, Nmask: %d\n", N, Nmask); + } + + if (!Quiet) gprint (GP_ERR, "z = "); + int n = 0; + for (int ny = 0; ny <= order; ny++) { + for (int nx = 0; nx <= order; nx++) { + sprintf (name, "CX%dY%d", nx, ny); + set_variable (name, b[n][0]); + if (!Quiet) gprint (GP_ERR, "%f x^%d y^%d ", b[n][0], nx, ny); + n++; + } + } + sprintf (name, "Cnn"); + set_variable (name, (double) order); + + if (!Quiet) gprint (GP_ERR, "\n"); + + /* free internal data */ + free (zfit); + free (mask); + + for (int i = 0; i < nterm; i++) { + free (c[i]); + free (b[i]); + } + free (b); + free (c); + + return (TRUE); +} Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 41891) @@ -7,4 +7,5 @@ int applyfit1d PROTO((int, char **)); int applyfit2d PROTO((int, char **)); +int applyfit2d_full PROTO((int, char **)); int applyfit3d PROTO((int, char **)); int box PROTO((int, char **)); @@ -13,4 +14,5 @@ int center PROTO((int, char **)); int cast PROTO((int, char **)); +int chebyshev_command PROTO((int, char **)); int circstats PROTO((int, char **)); int clear PROTO((int, char **)); @@ -45,4 +47,5 @@ int fit1d_irls PROTO((int, char **)); int fit2d PROTO((int, char **)); +int fit2d_full PROTO((int, char **)); int fit3d PROTO((int, char **)); int gaussjordan PROTO((int, char **)); @@ -204,4 +207,5 @@ {1, "applyfit1d", applyfit1d, "apply 1-d fit to new vector"}, {1, "applyfit2d", applyfit2d, "apply 2-d fit to new vector"}, + {1, "applyfit2d_full", applyfit2d_full, "apply 2-d fit to new vector"}, {1, "applyfit3d", applyfit3d, "apply 3-d fit to new vector"}, {1, "bisection", bisection, "use bisection to find threshold in vector"}, @@ -211,4 +215,5 @@ {1, "center", center, "center image on coords"}, {1, "cast", cast, "cast input vector to specified type"}, + {1, "chebyshev", chebyshev_command, "Chebyshev polynomial commands"}, {1, "circstats", circstats, "circular statistics"}, {1, "clear", clear, "erase plot"}, @@ -244,4 +249,5 @@ {1, "fit1d_irls", fit1d_irls, "fit polynomial to vector pair"}, {1, "fit2d", fit2d, "fit 2-d polynomial to vector triplet"}, + {1, "fit2d_full", fit2d_full, "fit 2-d polynomial to vector triplet (all X^n Y^m coeffs)"}, {1, "fit3d", fit3d, "fit 3-d polynomial to vector quad"}, {1, "gaussdev", gaussdeviate, "generate a gaussian deviate vector"}, Index: trunk/Ohana/src/opihi/cmd.data/medimage_calc.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/medimage_calc.c (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/medimage_calc.c (revision 41891) @@ -19,5 +19,5 @@ static float *irls_testval = NULL; -enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN}; +enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN, CALC_INNER_FRACTION}; int medimage_calc (int argc, char **argv) { @@ -39,4 +39,7 @@ } + float minRange = NAN; + float maxRange = NAN; + int mode = CALC_MEDIAN; if ((N = get_argument (argc, argv, "-mean"))) { @@ -45,15 +48,28 @@ } if ((N = get_argument (argc, argv, "-irls"))) { - if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; } + if (mode != CALC_MEDIAN) goto choose_one; mode = CALC_IRLS; remove_argument (N, &argc, argv); } if ((N = get_argument (argc, argv, "-wtmean"))) { - if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; } + if (mode != CALC_MEDIAN) goto choose_one; mode = CALC_WTMEAN; remove_argument (N, &argc, argv); } + if ((N = get_argument (argc, argv, "-inner-fraction"))) { + if (mode != CALC_MEDIAN) goto choose_one; + if (argc < N + 3) goto fraction_options; + mode = CALC_INNER_FRACTION; + remove_argument (N, &argc, argv); + minRange = atof(argv[N]); + if ((minRange < 0.0) || (minRange > 1.0)) goto fraction_options; + remove_argument (N, &argc, argv); + maxRange = atof(argv[N]); + if ((maxRange < 0.0) || (maxRange > 1.0)) goto fraction_options; + if (maxRange <= minRange) goto fraction_options; + remove_argument (N, &argc, argv); + } if ((N = get_argument (argc, argv, "-median"))) { - if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean\n"); return FALSE; } + if (mode != CALC_MEDIAN) goto choose_one; mode = CALC_MEDIAN; remove_argument (N, &argc, argv); @@ -157,4 +173,25 @@ } break; + case CALC_INNER_FRACTION: + fsort (val, N); + int Ns = MIN(MAX(0, minRange * N),N); // e.g., 0.1 * 50 = 5 + int Ne = MIN(MAX(0, maxRange * N),N); // e.g., 0.9 * 50 = 45 + int Npt = Ne - Ns; + + float sum = 0.0; + for (n = Ns; n < Ne; n++) { + sum += val[n]; + } + outvalue[Npix] = sum / (float) Npt; + + if (varvalue) { + float sum = 0.0; + for (n = Ns; n < Ne; n++) { + sum += SQ(val[n] - outvalue[Npix]); + } + // variance on the mean (stdev / sqrt(N))^2 + varvalue[Npix] = sum / (Npt - 1) / Npt; + } + break; case CALC_MEAN: { float sum = 0.0; @@ -195,4 +232,12 @@ irls_free(); return TRUE; + + choose_one: + gprint (GP_ERR, "supply only one of -mean, -median, -irls, -wtmean, -inner-fraction\n"); + return FALSE; + + fraction_options: + gprint (GP_ERR, "USE: -inner-fraction (min) (max) : min & max in range 0.0 to 1.0\n"); + return FALSE; } Index: trunk/Ohana/src/opihi/cmd.data/tdhistogram.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/tdhistogram.c (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/tdhistogram.c (revision 41891) @@ -45,6 +45,8 @@ int valid = (!reuse && (argc == 11)) || (reuse && (argc == 5)); if (!valid) { - gprint (GP_ERR, "USAGE: tdhistogram buffer x y z (Xmin) (Xmax) (Ymin) (Ymax) (Zmin) (Zmax) [-dx dx] [-dy dy] [-dz dz]\n"); + gprint (GP_ERR, "USAGE: tdhistogram buffer x y z (Xmin) (Xmax) (Ymin) (Ymax) (Zmin) (Zmax) [-dx dx] [-dy dy] [-dz dz] [-range range]\n"); gprint (GP_ERR, " OR: tdhistogram buffer x y z -reuse\n"); + gprint (GP_ERR, " -dx, -dy, -dy specify the bin size in these directions\n"); + gprint (GP_ERR, " -range : generate an output vector corresponding to the elements in the z-direction\n"); gprint (GP_ERR, " output buffer is 3D\n"); return (FALSE); Index: trunk/Ohana/src/opihi/cmd.data/vstats.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vstats.c (revision 41801) +++ trunk/Ohana/src/opihi/cmd.data/vstats.c (revision 41891) @@ -48,6 +48,8 @@ } - if (argc != 2) { - gprint (GP_ERR, "USAGE: vstat (vector) [-ignore value] [-q] [-quiet] [-iter Niter] [-sigma-clip Nsigma]\n"); + int valid = (argc == 2); + valid |= (argc > 3) && !strcmp (argv[2], "where"); + if (!valid) { + gprint (GP_ERR, "USAGE: vstat (vector) [-ignore value] [-q] [-quiet] [-iter Niter] [-sigma-clip Nsigma] [where (logical expression)]\n"); gprint (GP_ERR, " default is 1 iteration without sigma clipping or 3 with sigma clipping\n"); return (FALSE); @@ -58,21 +60,48 @@ /* we need two passes, one for max, min, mean, sum, one for median, stdev, etc */ - // set a good / bad mask + // tvec is used for logical test (truth vector) + Vector *tvec = NULL; + if (argc > 3) { + int size; + char *out = dvomath (argc - 3, &argv[3], &size, 1); + if (out == NULL) { + print_error (); + return FALSE; + } + if ((tvec = SelectVector (out, OLDVECTOR, TRUE)) == NULL) { + gprint (GP_ERR, " invalid logic result\n"); + DeleteNamedVector (out); + free (out); + return (FALSE); + } + if (tvec->Nelements != vec->Nelements) { + gprint (GP_ERR, "logical vector does not match in length\n"); + return FALSE; + } + } + + // set a good / bad mask (mask == 1 means 'bad') ALLOCATE (mask, char, vec[0].Nelements); if (vec[0].type == OPIHI_FLT) { opihi_flt *X = vec[0].elements.Flt; for (i = 0; i < vec[0].Nelements; i++, X++) { - mask[i] = 1; - if (!finite (*X)) continue; - if (Ignore && (*X == IgnoreValue)) continue; - mask[i] = 0; + mask[i] = 0; // do not mask unless we have a reason below + if (tvec) { + // note the logical vector is 0 == bad (ignore) while mask[i] has the opposite sense (0 is good) + mask[i] = (tvec->type == OPIHI_FLT) ? (tvec->elements.Flt[i] == 0.0) : (tvec->elements.Int[i] == 0); + } + if (!finite (*X)) { mask[i] = 1; } + if (Ignore && (*X == IgnoreValue)) { mask[i] = 1; } } } else { opihi_int *X = vec[0].elements.Int; for (i = 0; i < vec[0].Nelements; i++, X++) { - mask[i] = 1; - if (!finite (*X)) continue; - if (Ignore && (*X == IgnoreValue)) continue; - mask[i] = 0; + mask[i] = 0; // do not mask unless we have a reason below + if (tvec) { + // note the logical vector is 0 == bad (ignore) while mask[i] has the opposite sense (0 is good) + mask[i] = (tvec->type == OPIHI_FLT) ? (tvec->elements.Flt[i] == 0.0) : (tvec->elements.Int[i] == 0); + } + if (!finite (*X)) { mask[i] = 1; } + if (Ignore && (*X == IgnoreValue)) { mask[i] = 1; } } } Index: trunk/Ohana/src/opihi/include/data.h =================================================================== --- trunk/Ohana/src/opihi/include/data.h (revision 41801) +++ trunk/Ohana/src/opihi/include/data.h (revision 41891) @@ -71,4 +71,15 @@ float **dNabla_w; // a matrix of Nabla_w values for each layer } Nnet; + +// a single chebyshev associates a name with the scaling +// relationship and maybe other stuff +typedef struct { + char *name; + opihi_flt scale[2]; + opihi_flt zero[2]; + int order; + int dimen; + opihi_flt *A; +} ChebyshevType; void InitData (void); @@ -294,3 +305,20 @@ int nnet_apply (int argc, char **argv); +/*** Chebyshev functions (chebyshev.c, chebyshev_commands.c) ***/ + +void InitChebyshevs (); +void FreeChebyshevs (); +void FreeChebyshev (ChebyshevType *chebyshev); +ChebyshevType *FindChebyshev (char *name); +ChebyshevType *CreateChebyshev (char *name); +int DeleteChebyshev (ChebyshevType *chebyshev); +void ListChebyshevs (); +int ChebyshevSetScale (ChebyshevType *cheb, Vector *vec, int dir); +Vector *ChebyshevNormVector (ChebyshevType *cheb, Vector *vec, int dir); +Vector *ChebyshevPolyVector (Vector *Ti, Vector *vec, int order); +int ChebyshevPolyFit1D (ChebyshevType *cheb, Vector *vec, Vector **poly); +int ChebyshevPolyApplyFit1D (ChebyshevType *cheb, Vector *vec, Vector **poly); +int ChebyshevPolyFit2D (ChebyshevType *cheb, Vector *vec, Vector **xPoly, Vector **yPoly); +int ChebyshevPolyApplyFit2D (ChebyshevType *cheb, Vector *vec, Vector **xPoly, Vector **yPoly); + # endif Index: trunk/Ohana/src/opihi/lib.data/ChebyshevOps.c =================================================================== --- trunk/Ohana/src/opihi/lib.data/ChebyshevOps.c (revision 41891) +++ trunk/Ohana/src/opihi/lib.data/ChebyshevOps.c (revision 41891) @@ -0,0 +1,366 @@ +# include "data.h" + +/* this file contains functions to manage the median image data + */ + +ChebyshevType **chebyshevs = NULL; /* book to store the list of all splines */ +int Nchebyshevs; /* number of currently defined chebyshevs */ +int NCHEBYSHEVS; /* number of currently allocated chebyshevs */ + +void InitChebyshevs () { + Nchebyshevs = 0; + NCHEBYSHEVS = 16; + ALLOCATE (chebyshevs, ChebyshevType *, NCHEBYSHEVS); +} + +void FreeChebyshevs () { + + for (int i = 0; i < Nchebyshevs; i++) { + FreeChebyshev (chebyshevs[i]); + } + FREE (chebyshevs); +} + +void FreeChebyshev (ChebyshevType *chebyshev) { + + if (!chebyshev) return; + + FREE (chebyshev[0].name); + FREE (chebyshev[0].A); + free (chebyshev); +} + +/* return the given chebyshev */ +ChebyshevType *FindChebyshev (char *name) { + + if (!chebyshevs) return NULL; + + for (int i = 0; i < Nchebyshevs; i++) { + if (!strcmp (chebyshevs[i][0].name, name)) { + return (chebyshevs[i]); + } + } + return (NULL); +} + +/* make a new named chebyshev */ +ChebyshevType *CreateChebyshev (char *name) { + + if (!chebyshevs) InitChebyshevs (); + + // do not create a chebyshev if one exists + ChebyshevType *chebyshev = FindChebyshev (name); + if (chebyshev != NULL) { + return NULL; + } + + int N = Nchebyshevs; + Nchebyshevs ++; + CHECK_REALLOCATE (chebyshevs, ChebyshevType *, NCHEBYSHEVS, Nchebyshevs, 16); + ALLOCATE (chebyshev, ChebyshevType, 1); + + chebyshev->name = strcreate (name); + chebyshev->zero[0] = chebyshev->zero[1] = NAN; + chebyshev->scale[0] = chebyshev->scale[1] = NAN; + chebyshev->order = -1; + chebyshev->dimen = -1; + chebyshev->A = NULL; + + chebyshevs[N] = chebyshev; + return (chebyshev); +} + +/* delete a chebyshev */ +int DeleteChebyshev (ChebyshevType *chebyshev) { + + if (!chebyshevs) return TRUE; // 'true' otherwise delete will abort if none yet defined + + /* find chebyshev in chebyshev list */ + int N = -1; + for (int i = 0; i < Nchebyshevs; i++) { + if (chebyshevs[i] == chebyshev) { + N = i; + break; + } + } + if (N == -1) return (FALSE); + + for (int i = N; i < Nchebyshevs - 1; i++) { + chebyshevs[i] = chebyshevs[i + 1]; + } + Nchebyshevs --; + int NCHEBYSHEVS_2 = MAX (16, NCHEBYSHEVS / 2); + if (Nchebyshevs < NCHEBYSHEVS_2) { + NCHEBYSHEVS = NCHEBYSHEVS_2; + REALLOCATE (chebyshevs, ChebyshevType *, NCHEBYSHEVS); + } + + FreeChebyshev (chebyshev); + return (TRUE); +} + +/* list known chebyshevs */ +void ListChebyshevs () { + + if (!chebyshevs) return; + + for (int i = 0; i < Nchebyshevs; i++) { + gprint (GP_ERR, "%-15s %2d %2d (%15.8e %15.8e)", chebyshevs[i][0].name, chebyshevs[i][0].order, chebyshevs[i][0].dimen, chebyshevs[i][0].zero[0], chebyshevs[i][0].scale[0]); + if (chebyshevs[i][0].dimen == 2) { gprint (GP_ERR, " (%15.8e %15.8e)", chebyshevs[i][0].zero[1], chebyshevs[i][0].scale[1]); } + gprint (GP_ERR, "\n"); + } + return; +} + +/**** below are technical operations ****/ + +int ChebyshevSetScale (ChebyshevType *cheb, Vector *vec, int dir) { + + if (dir >= 2) { + gprint (GP_ERR, "invalid direction %d\n", dir); + return FALSE; + } + + // find the min and max of the vector + opihi_flt minValue = NAN; + opihi_flt maxValue = NAN; + + for (int i = 0; i < vec->Nelements; i++) { + if (isnan(vec->elements.Flt[i])) continue; + if (isnan(minValue)) { minValue = vec->elements.Flt[i]; } + if (isnan(maxValue)) { maxValue = vec->elements.Flt[i]; } + minValue = MIN(minValue, vec->elements.Flt[i]); + maxValue = MAX(maxValue, vec->elements.Flt[i]); + } + if (minValue == maxValue) { + gprint (GP_ERR, "insufficient data range to determine scale factors\n"); + return FALSE; + } + + cheb->scale[dir] = 2.0 / (maxValue - minValue); + cheb->zero[dir] = 1 - cheb->scale[dir] * maxValue; + return TRUE; +} + +Vector *ChebyshevNormVector (ChebyshevType *cheb, Vector *xvec, int dir) { + + Vector *xNorm = InitVector(); + ResetVector (xNorm, OPIHI_FLT, xvec->Nelements); + + for (int i = 0; i < xvec->Nelements; i++) { + xNorm->elements.Flt[i] = xvec->elements.Flt[i]*cheb->scale[dir] + cheb->zero[dir]; + } + return xNorm; +} + +Vector *ChebyshevPolyVector (Vector *Ti, Vector *vec, int order) { + + if (Ti == NULL) { Ti = InitVector(); } + ResetVector (Ti, OPIHI_FLT, vec->Nelements); + + // quick but inefficient implementation + switch (order) { + case 0: + for (int i = 0; i < vec->Nelements; i++) { Ti->elements.Flt[i] = 1.0; } break; + case 1: + for (int i = 0; i < vec->Nelements; i++) { Ti->elements.Flt[i] = vec->elements.Flt[i]; } break; + case 2: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = 2.0*x2 - 1.0; } break; + case 3: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = vec->elements.Flt[i]*(4.0*x2 - 3.0); } break; + case 4: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = x2*(8.0*x2 - 8.0) + 1.0; } break; + case 5: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = vec->elements.Flt[i] * (x2*(16.0*x2 - 20.0) + 5.0); } break; + case 6: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = x2*(x2*(32.0*x2 - 48.0) + 18.0) - 1.0; } break; + case 7: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = vec->elements.Flt[i] * (x2*(x2*(64.0*x2 - 112.0) + 56.0) - 7.0); } break; + case 8: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = x2*(x2*(x2*(128.0*x2 - 256.0) + 160.0) - 32.0) + 1.0; } break; + case 9: + for (int i = 0; i < vec->Nelements; i++) { opihi_flt x2 = SQ(vec->elements.Flt[i]); Ti->elements.Flt[i] = vec->elements.Flt[i] *(x2*(x2*(x2*(256.0*x2 - 576.0) + 432.0) - 129.0) + 9.0); } break; + default: + gprint (GP_ERR, "Chebyshev orders higher than 9 are not yet coded\n"); + FREE (Ti); + return NULL; + } + + return Ti; +} + + +// we are supplied a vector of length vec->Nelements +// poly is an array of nterm (= order + 1) chebyshev polynomials +int ChebyshevPolyFit1D (ChebyshevType *cheb, Vector *vec, Vector **poly) { + + /* nterm is number of polynomial terms, starting at x^0 */ + int nterm = cheb->order + 1; + + ALLOCATE_PTR (b, double *, nterm); + ALLOCATE_PTR (c, double *, nterm); + for (int i = 0; i < nterm; i++) { + ALLOCATE (c[i], double, nterm); + ALLOCATE (b[i], double, 1); + memset (c[i], 0, nterm*sizeof(double)); + memset (b[i], 0, sizeof(double)); + } + + opihi_flt *y = vec[0].elements.Flt; + + for (int i = 0; i < vec->Nelements; i++, y++) { + if (!finite(*y)) continue; + + // only calculate the upper diagonal + for (int j = 0; j < nterm; j++) { + b[j][0] += *y * poly[j]->elements.Flt[i]; + for (int k = j; k < nterm; k++) { + c[j][k] += poly[j]->elements.Flt[i]*poly[k]->elements.Flt[i]; + } + } + } + + // fill in the lower diagonal + for (int j = 0; j < nterm; j++) { + for (int k = 0; k < j; k++) { + c[j][k] = c[k][j]; + } + } + + if (!dgaussjordan (c, b, nterm, 1)) { + gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n"); + return FALSE; + } + + if (cheb->A) FREE (cheb->A); + ALLOCATE (cheb->A, double, nterm); + for (int j = 0; j < nterm; j++) { + cheb->A[j] = b[j][0]; + } + + for (int i = 0; i < nterm; i++) { + FREE (c[i]); + FREE (b[i]); + } + FREE (b); + FREE (c); + + return TRUE; +} + +// we are supplied a vector of length vec->Nelements +// poly is an array of nterm (= order + 1) chebyshev polynomials +int ChebyshevPolyApplyFit1D (ChebyshevType *cheb, Vector *vec, Vector **poly) { + + int nterm = cheb->order + 1; + + for (int i = 0; i < vec->Nelements; i++) { + opihi_flt value = 0.0; + for (int j = 0; j < nterm; j++) { + value += cheb->A[j] * poly[j]->elements.Flt[i]; + } + vec->elements.Flt[i] = value; + } + + return TRUE; +} + +// we are supplied a vector of length vec->Nelements +// poly is an array of nterm (= order + 1) chebyshev polynomials +int ChebyshevPolyFit2D (ChebyshevType *cheb, Vector *vec, Vector **xPoly, Vector **yPoly) { + + /* nterm is number of polynomial terms, starting at x^0 */ + int nterm = SQ(cheb->order + 1); + + ALLOCATE_PTR (b, double *, nterm); + ALLOCATE_PTR (c, double *, nterm); + for (int i = 0; i < nterm; i++) { + ALLOCATE (c[i], double, nterm); + ALLOCATE (b[i], double, 1); + memset (c[i], 0, nterm*sizeof(double)); + memset (b[i], 0, sizeof(double)); + } + + opihi_flt *y = vec[0].elements.Flt; + + for (int i = 0; i < vec->Nelements; i++, y++) { + if (!finite(*y)) continue; + + // XXX can we only calculate the upper diagonal? + int n = 0; + for (int jx = 0; jx <= cheb->order; jx++) { + for (int jy = 0; jy <= cheb->order; jy++) { + + opihi_flt vj = xPoly[jx]->elements.Flt[i] * yPoly[jy]->elements.Flt[i]; + b[n][0] += *y * vj; + + int m = 0; + for (int kx = 0; kx <= cheb->order; kx++) { + for (int ky = 0; ky <= cheb->order; ky++) { + c[n][m] += vj * xPoly[kx]->elements.Flt[i]*yPoly[ky]->elements.Flt[i]; + m++; + } + } + n++; + } + } + } + + for (int i = 0; i < nterm; i++) { + fprintf (stderr, "c: "); + for (int j = 0; j < nterm; j++) { + fprintf (stderr, "%e ", c[i][j]); + } + fprintf (stderr, " b : %e\n", b[i][0]); + } + + if (!dgaussjordan (c, b, nterm, 1)) { + gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n"); + return FALSE; + } + + for (int i = 0; i < nterm; i++) { + fprintf (stderr, "c: "); + for (int j = 0; j < nterm; j++) { + fprintf (stderr, "%e ", c[i][j]); + } + fprintf (stderr, " b : %e\n", b[i][0]); + } + + if (cheb->A) FREE (cheb->A); + ALLOCATE (cheb->A, double, nterm); + for (int j = 0; j < nterm; j++) { + cheb->A[j] = b[j][0]; + } + + for (int i = 0; i < nterm; i++) { + FREE (c[i]); + FREE (b[i]); + } + FREE (b); + FREE (c); + + return TRUE; +} + +// we are supplied a vector of length vec->Nelements +// poly is an array of nterm (= order + 1) chebyshev polynomials +int ChebyshevPolyApplyFit2D (ChebyshevType *cheb, Vector *vec, Vector **xPoly, Vector **yPoly) { + + for (int i = 0; i < vec->Nelements; i++) { + opihi_flt value = 0.0; + + int n = 0; + for (int jx = 0; jx <= cheb->order; jx++) { + for (int jy = 0; jy <= cheb->order; jy++) { + opihi_flt vj = xPoly[jx]->elements.Flt[i] * yPoly[jy]->elements.Flt[i]; + value += cheb->A[n] * vj; + n ++; + } + } + vec->elements.Flt[i] = value; + } + + return TRUE; +} + Index: trunk/Ohana/src/opihi/lib.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/lib.data/Makefile (revision 41801) +++ trunk/Ohana/src/opihi/lib.data/Makefile (revision 41891) @@ -25,4 +25,5 @@ $(SDIR)/SplineOps.$(ARCH).o \ $(SDIR)/MedImageOps.$(ARCH).o \ +$(SDIR)/ChebyshevOps.$(ARCH).o \ $(SDIR)/mrqmin.$(ARCH).o \ $(SDIR)/mrq2dmin.$(ARCH).o \ Index: trunk/Ohana/src/photdbc/src/args.c =================================================================== --- trunk/Ohana/src/photdbc/src/args.c (revision 41801) +++ trunk/Ohana/src/photdbc/src/args.c (revision 41891) @@ -88,4 +88,5 @@ MAX_MIN_MAG = atof(argv[N]); remove_argument (N, &argc, argv); + fprintf (stderr, "warning -maxminmag may not do what you think: check in make_subcatalog.c\n"); } Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 41801) +++ trunk/Ohana/src/relphot/src/args.c (revision 41891) @@ -525,7 +525,6 @@ STAGES &= ~STAGE_STACK; } - if (!STAGES) { - fprintf (stderr, "ERROR: no valid stages selected\n"); - exit (3); + if (!STAGES && NLOOP > 1) { + fprintf (stderr, "WARNING: no averages requested (all stages skipped), but nloop > 1. is this wasted effort?\n"); } Index: trunk/Ohana/src/uniphot/Makefile =================================================================== --- trunk/Ohana/src/uniphot/Makefile (revision 41801) +++ trunk/Ohana/src/uniphot/Makefile (revision 41891) @@ -1,3 +1,3 @@ -default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client +default: uniphot setphot setphot_client setfwhm setposangle setposangle_client setastrom setastrom_client setgalmodel setgalmodel_client fixhsc fixhsc_client fiximids fiximids_client fixstkids fixstkids_client ckids ckids_client help: @@ -28,4 +28,6 @@ setgalmodel: $(BIN)/setgalmodel.$(ARCH) setgalmodel_client: $(BIN)/setgalmodel_client.$(ARCH) +fixhsc: $(BIN)/fixhsc.$(ARCH) +fixhsc_client: $(BIN)/fixhsc_client.$(ARCH) fiximids: $(BIN)/fiximids.$(ARCH) fiximids_client: $(BIN)/fiximids_client.$(ARCH) @@ -35,5 +37,5 @@ ckids_client: $(BIN)/ckids_client.$(ARCH) -install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client $(DESTBIN)/setgalmodel $(DESTBIN)/setgalmodel_client $(DESTBIN)/fiximids $(DESTBIN)/fiximids_client $(DESTBIN)/fixstkids $(DESTBIN)/fixstkids_client $(DESTBIN)/ckids $(DESTBIN)/ckids_client +install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client $(DESTBIN)/setposangle $(DESTBIN)/setposangle_client $(DESTBIN)/setastrom $(DESTBIN)/setastrom_client $(DESTBIN)/setgalmodel $(DESTBIN)/setgalmodel_client $(DESTBIN)/fixhsc $(DESTBIN)/fixhsc_client $(DESTBIN)/fiximids $(DESTBIN)/fiximids_client $(DESTBIN)/fixstkids $(DESTBIN)/fixstkids_client $(DESTBIN)/ckids $(DESTBIN)/ckids_client UNIPHOT = \ @@ -202,4 +204,30 @@ $(BIN)/setfwhm.$(ARCH): $(SETFWHM) +FIXHSC = \ +$(SRC)/fixhsc.$(ARCH).o \ +$(SRC)/initialize_fixhsc.$(ARCH).o \ +$(SRC)/load_images_fixhsc.$(ARCH).o \ +$(SRC)/load_rules_fixhsc.$(ARCH).o \ +$(SRC)/update_dvo_fixhsc.$(ARCH).o \ +$(SRC)/update_catalog_fixhsc.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o + +$(FIXHSC): $(INC)/fixhsc.h +$(BIN)/fixhsc.$(ARCH): $(FIXHSC) + +FIXHSC_CLIENT = \ +$(SRC)/fixhsc_client.$(ARCH).o \ +$(SRC)/initialize_fixhsc.$(ARCH).o \ +$(SRC)/load_images_fixhsc.$(ARCH).o \ +$(SRC)/load_rules_fixhsc.$(ARCH).o \ +$(SRC)/update_dvo_fixhsc.$(ARCH).o \ +$(SRC)/update_catalog_fixhsc.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o + +$(FIXHSC_CLIENT): $(INC)/fixhsc.h +$(BIN)/fixhsc_client.$(ARCH): $(FIXHSC_CLIENT) + FIXIMIDS = \ $(SRC)/fiximids.$(ARCH).o \ Index: trunk/Ohana/src/uniphot/include/fixhsc.h =================================================================== --- trunk/Ohana/src/uniphot/include/fixhsc.h (revision 41891) +++ trunk/Ohana/src/uniphot/include/fixhsc.h (revision 41891) @@ -0,0 +1,42 @@ +# include +# include +# include + +/* global variables set in parameter file */ +# define DVO_MAX_PATH 1024 +char ImageCat[DVO_MAX_PATH]; +char *CATDIR; +int HOST_ID; +char *HOSTDIR; +char *SINGLE_CPT; +int VERBOSE; +int UPDATE; +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; + +SkyRegion UserPatch; + +/***** prototypes ****/ +int main PROTO((int argc, char **argv)); +void ConfigInit PROTO((int *argc, char **argv)); + +void initialize_fixhsc PROTO((int argc, char **argv)); +void initialize_fixhsc_client PROTO((int argc, char **argv)); + +int load_rules_fixhsc PROTO((char *filename)); +int get_rules_fixhsc PROTO((e_time time)); + +void set_db PROTO((FITS_DB *in)); +int Shutdown PROTO((char *format, ...) OHANA_FORMAT(printf, 1, 2) ); +void TrapSignal PROTO((int sig)); +void SetProtect PROTO((int mode)); +int SetSignals PROTO((void)); + +Image *load_images_fixhsc PROTO((FITS_DB *db, off_t *Nimage)); +void update_images_fixhsc PROTO((Image *image, off_t Nimage)); + +int update_dvo_fixhsc PROTO((void)); +int update_dvo_fixhsc_parallel PROTO((SkyTable *sky)); + +void update_catalog_fixhsc PROTO((Catalog *catalog)); Index: trunk/Ohana/src/uniphot/src/fixhsc.c =================================================================== --- trunk/Ohana/src/uniphot/src/fixhsc.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/fixhsc.c (revision 41891) @@ -0,0 +1,64 @@ +# include "fixhsc.h" + +int main (int argc, char **argv) { + + off_t Nimage; + int status; + FITS_DB db; + + /* get configuration info, args, lockfile */ + SetSignals (); + initialize_fixhsc (argc, argv); + + set_db (&db); + gfits_db_init (&db); + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + + // load images + Image *image = load_images_fixhsc (&db, &Nimage); + + // need to update image table here: + update_images_fixhsc (image, Nimage); + + // write image table (even if some remote clients failed) + if (UPDATE) { + SetProtect (TRUE); + dvo_image_save (&db, VERBOSE); + dvo_image_unlock (&db); + SetProtect (FALSE); + } + + status = update_dvo_fixhsc (); + + if (!status) exit (1); + exit (0); +} + + +/* fixhsc : modify the photcode values based on the time perido + + current HSC photcodes: + g = 20000 - 20111 + r = 21000 - 21111 + + goal HSC photcodes: + g ,North = 20000 - 20111 + g ,South = 20200 - 20311 + r ,North = 21000 - 21111 + r ,South = 21200 - 21311 + r2,North = 21400 - 21511 + r2,South = 21600 - 21711 + + I can define rules which specify the goal photcode for a time period, + then the code can use the value % 1000 to get the chip ID + + rule file: + start stop photcode (start & stop in MJD, photcode = value for chip 00, e.g., 21200) + + ** load images (save subset with imageID, time, photcode, seq) + ** load catalogs + ** update detection + ** count detections / image ID? + */ Index: trunk/Ohana/src/uniphot/src/fixhsc_client.c =================================================================== --- trunk/Ohana/src/uniphot/src/fixhsc_client.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/fixhsc_client.c (revision 41891) @@ -0,0 +1,21 @@ +# include "fixhsc.h" + +// fixhsc_client is run on a remote host and is responsible for updating the catalogs +// owned by that host. + +// fixhsc_client -hostname (hostname) -catdir (catdir) -hostdir (hostdir) -images (images) +// load the SkyTable and HostTable from (catdir) [contains the host responsibilities] +// load the image data from (images) [contains the zp infomation to apply] +// loop over tables from SkyTable with my host ID +// the calling program tells the client which host they are (or host ID?); we do not actually have to run this on the real host + +int main (int argc, char **argv) { + + // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, IMAGES) + SetSignals (); + initialize_fixhsc_client (argc, argv); + + update_dvo_fixhsc (); + + exit (0); +} Index: trunk/Ohana/src/uniphot/src/initialize_fixhsc.c =================================================================== --- trunk/Ohana/src/uniphot/src/initialize_fixhsc.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/initialize_fixhsc.c (revision 41891) @@ -0,0 +1,272 @@ +# include "fixhsc.h" +int args_fixhsc (int argc, char **argv); +int args_fixhsc_client (int argc, char **argv); + +void usage_fixhsc () { + fprintf (stderr, "USAGE: fixhsc rules.txt [options]\n"); + fprintf (stderr, " rules.txt : file with time-based photcode rules:\n\n"); + fprintf (stderr, " options:\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -cpt n0000/0000.00.cpt : only repair this cpt file \n"); + fprintf (stderr, " -region Rmin Rmax Dmin Dmax : only repair this region\n"); + fprintf (stderr, " -update : actually write results to detections tables\n"); + fprintf (stderr, " -parallel : run in parallel mode\n\n"); + + fprintf (stderr, " -h : this help list\n"); + fprintf (stderr, " -help : this help list\n"); + fprintf (stderr, " --h : this help list\n"); + fprintf (stderr, " --help : this help list\n"); + + fprintf (stderr, " Note that the dvo db can be specified by -D CATDIR (directory)\n"); + exit (2); +} + +void initialize_fixhsc (int argc, char **argv) { + + int N; + struct stat statbuffer; + char CatdirPhotcodeFile[256]; + + if ((N = get_argument (argc, argv, "-h"))) usage_fixhsc(); + if ((N = get_argument (argc, argv, "-help"))) usage_fixhsc(); + if ((N = get_argument (argc, argv, "--h"))) usage_fixhsc(); + if ((N = get_argument (argc, argv, "--help"))) usage_fixhsc(); + + /* are these set correctly? */ + ConfigInit (&argc, argv); + args_fixhsc (argc, argv); + + if (stat (CATDIR, &statbuffer)) { + fprintf (stderr, "error accessing dvo database directory '%s'\n", CATDIR); + exit (1); + } + + // load the photcode table + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, NULL, FALSE)) { + fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); + exit (1); + } +} + +int args_fixhsc (int argc, char **argv) { + + int N; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + SINGLE_CPT = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + SINGLE_CPT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + // region of interest + 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); + } + + 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, + // fixhsc 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 fixhsc_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); + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + if (argc != 2) usage_fixhsc(); + + load_rules_fixhsc (argv[1]); + + return (TRUE); +} + +void usage_fixhsc_client () { + fprintf (stderr, "USAGE: fixhsc_client -hostID (hostID) -catdir (catdir) -hostdir (hostdir) -images (images) [options]\n"); + fprintf (stderr, " options:\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -update : actually write results to detections tables\n"); + exit (2); +} + +void initialize_fixhsc_client (int argc, char **argv) { + + int N; + struct stat statbuffer; + char CatdirPhotcodeFile[256]; + + if ((N = get_argument (argc, argv, "-h"))) usage_fixhsc_client(); + if ((N = get_argument (argc, argv, "-help"))) usage_fixhsc_client(); + if ((N = get_argument (argc, argv, "--h"))) usage_fixhsc_client(); + if ((N = get_argument (argc, argv, "--help"))) usage_fixhsc_client(); + + SetZeroPoint (25.0); // XXX is this needed? + args_fixhsc_client (argc, argv); + + if (stat (CATDIR, &statbuffer)) { + fprintf (stderr, "error accessing dvo database directory '%s'\n", CATDIR); + exit (1); + } + + // load the photcode table : XXX needed? + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, NULL, FALSE)) { + fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); + exit (1); + } +} + +int args_fixhsc_client (int argc, char **argv) { + + int N; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + SINGLE_CPT = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + SINGLE_CPT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + // region of interest + 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); + } + + 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) usage_fixhsc_client(); + + 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) usage_fixhsc_client(); + + if ((N = get_argument (argc, argv, "-catdir"))) { + remove_argument (N, &argc, argv); + CATDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!CATDIR) usage_fixhsc_client(); + + if (argc != 1) usage_fixhsc_client(); + return (TRUE); +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr); + +void ConfigInit (int *argc, char **argv) { + + double ZERO_POINT; + char *config, *file; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + // force CATDIR to be absolute (so parallel mode will work) + char *tmpcatdir = NULL; + ALLOCATE (tmpcatdir, char, DVO_MAX_PATH); + GetConfig (config, "CATDIR", "%s", 0, tmpcatdir); + CATDIR = abspath (tmpcatdir, DVO_MAX_PATH); + free (tmpcatdir); + + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + ScanConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); + SetZeroPoint (ZERO_POINT); + + free (config); + free (file); +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: trunk/Ohana/src/uniphot/src/load_images_fixhsc.c =================================================================== --- trunk/Ohana/src/uniphot/src/load_images_fixhsc.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/load_images_fixhsc.c (revision 41891) @@ -0,0 +1,46 @@ +# include "fixhsc.h" + +Image *load_images_fixhsc (FITS_DB *db, off_t *Nimage_load) { + + Image *image; + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + /* read entire db table */ + if (!dvo_image_load (db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db[0].filename); + + /* use a vtable to keep the images to be calibrated */ + image = gfits_table_get_Image (&db[0].ftable, Nimage_load, &db[0].scaledValue, &db[0].nativeOrder); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + + fprintf (stderr, "loaded "OFF_T_FMT" images\n", *Nimage_load); + + return (image); +} + +void update_images_fixhsc (Image *image, off_t Nimage) { + + for (off_t i = 0; i < Nimage; i++) { + + // skip photcode == 0 images (PHU) + short photcode = image[i].photcode; + if (!photcode) continue; + + // we have a measure with a given photcode and time: + e_time time = image[i].tzero; + + int newbase = get_rules_fixhsc (time); + if (!newbase) continue; + + // the rule is the new base photcode + // this is robust against multiple passes + short chipcode = photcode % 200; + + short newcode = newbase + chipcode; + image[i].photcode = newcode; + } +} + Index: trunk/Ohana/src/uniphot/src/load_rules_fixhsc.c =================================================================== --- trunk/Ohana/src/uniphot/src/load_rules_fixhsc.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/load_rules_fixhsc.c (revision 41891) @@ -0,0 +1,58 @@ +# include "fixhsc.h" + +static int Nrule = 0; +static e_time *startRule = NULL; +static e_time *stopRule = NULL; +static int *photbaseRule = NULL; + +int load_rules_fixhsc (char *filename) { + + // the rule file contains + // Start Stop Photbase + + FILE *f = fopen (filename, "r"); + if (!f) { fprintf (stderr, "file %s not found", filename); abort(); } + + int NRULE = 10000; + + ALLOCATE (startRule, e_time, NRULE); + ALLOCATE (stopRule, e_time, NRULE); + ALLOCATE (photbaseRule, int, NRULE); + + double startMJD, stopMJD; + + // rules are in MJD, convert to e_time values + while (fscanf (f, "%lf %lf %d", &startMJD, &stopMJD, &photbaseRule[Nrule]) != EOF) { + + startRule[Nrule] = ohana_mjd_to_sec (startMJD); + stopRule[Nrule] = ohana_mjd_to_sec (stopMJD); + Nrule ++; + + if (Nrule >= NRULE) { + NRULE += 10000; + REALLOCATE (startRule, e_time, NRULE); + REALLOCATE (stopRule, e_time, NRULE); + REALLOCATE (photbaseRule, int, NRULE); + } + } + + // dsortpair (ruleR, ruleD, Nrule); + + return TRUE; +} + +int get_rules_fixhsc (e_time time) { + + // find the rule that encompases this time, return the rule + + for (int i = 0; i < Nrule; i++) { + if (time < startRule[i]) continue; + if (time > stopRule[i]) continue; + + return photbaseRule[i]; + } + + // if we return false, no match was found, keep the old photcode + return FALSE; + +} Index: trunk/Ohana/src/uniphot/src/update_catalog_fixhsc.c =================================================================== --- trunk/Ohana/src/uniphot/src/update_catalog_fixhsc.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/update_catalog_fixhsc.c (revision 41891) @@ -0,0 +1,22 @@ +# include "fixhsc.h" + +void update_catalog_fixhsc (Catalog *catalog) { + + for (int i = 0; i < catalog[0].Nmeasure; i++) { + + // we have a measure with a given photcode and time: + e_time time = catalog[0].measure[i].t; + + int newbase = get_rules_fixhsc (time); + if (!newbase) continue; + + unsigned short photcode = catalog[0].measure[i].photcode; + + // the rule is the new base photcode + // this is robust against multiple passes + short chipcode = photcode % 200; + + short newcode = newbase + chipcode; + catalog[0].measure[i].photcode = newcode; + } +} Index: trunk/Ohana/src/uniphot/src/update_dvo_fixhsc.c =================================================================== --- trunk/Ohana/src/uniphot/src/update_dvo_fixhsc.c (revision 41891) +++ trunk/Ohana/src/uniphot/src/update_dvo_fixhsc.c (revision 41891) @@ -0,0 +1,145 @@ +# include "fixhsc.h" + +int update_dvo_fixhsc (void) { + + SkyTable *sky = NULL; + SkyList *skylist = NULL; + Catalog catalog; + off_t i; + + // load the current sky table (layout of all SkyRegions) + sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, -1, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + if (PARALLEL && !HOST_ID) { + update_dvo_fixhsc_parallel (sky); + return TRUE; + } + + // determine the populated SkyRegions overlapping the requested area (default depth) + if (SINGLE_CPT) { + skylist = SkyRegionByCPT (sky, SINGLE_CPT); + } else { + skylist = SkyListByPatch (sky, -1, &UserPatch); + } + myAssert (skylist, "ooops!"); + + // update measurements for each populated catalog + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + char localFilename[1024]; + snprintf (localFilename, 1024, "%s/%s.cpt", HOSTDIR, skylist->regions[i]->name); + + // set up the basic catalog info + dvo_catalog_init (&catalog, TRUE); + catalog.filename = HOST_ID ? localFilename : skylist[0].filename[i]; + catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); + } + if (!catalog.Naverage_disk) { + if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + update_catalog_fixhsc (&catalog); + + if (!UPDATE) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename); + + SetProtect (TRUE); + if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); } + if (!dvo_catalog_unlock (&catalog)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog.filename); exit (1); } + SetProtect (FALSE); + dvo_catalog_free (&catalog); + } + + return (TRUE); +} + +# define DEBUG 1 + +int update_dvo_fixhsc_parallel (SkyTable *sky) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf_nowarn (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // now launch the fixhsc_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + + // fix this if you ever try on a parallel database + fprintf (stderr, "update code to supply rule file to clients\n"); + exit (2); + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // 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; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/fixhsc.results.%s.dat", table->hosts[i].pathname, uniquer); + + char *command = NULL; + strextend (&command, "fixhsc_client -hostID %d -catdir %s -hostdir %s -region %f %f %f %f", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + if (VERBOSE) { strextend (&command, "-v"); } + if (SINGLE_CPT) { strextend (&command, "-cpt %s", SINGLE_CPT); } + if (UPDATE) { strextend (&command, "-update"); } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running fixhsc_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 + } + } + + // wait for the remote jobs to be completed + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the fixhsc_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + if (!status) return FALSE; + } + + return (TRUE); +} +