Index: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/Makefile (revision 41827) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/Makefile (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/applyfit2d.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/applyfit2d.c (revision 41827) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/applyfit2d.c (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/chebyshev.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/chebyshev.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/chebyshev.c (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/chebyshev_commands.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/chebyshev_commands.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/chebyshev_commands.c (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/fit2d.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/fit2d.c (revision 41827) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/fit2d.c (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/init.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/init.c (revision 41827) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/cmd.data/init.c (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/include/data.h =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/include/data.h (revision 41827) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/include/data.h (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/lib.data/ChebyshevOps.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/lib.data/ChebyshevOps.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/lib.data/ChebyshevOps.c (revision 41830) @@ -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: branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/lib.data/Makefile =================================================================== --- branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/lib.data/Makefile (revision 41827) +++ branches/eam_branches/ipp-dev-20210817/Ohana/src/opihi/lib.data/Makefile (revision 41830) @@ -25,4 +25,5 @@ $(SDIR)/SplineOps.$(ARCH).o \ $(SDIR)/MedImageOps.$(ARCH).o \ +$(SDIR)/ChebyshevOps.$(ARCH).o \ $(SDIR)/mrqmin.$(ARCH).o \ $(SDIR)/mrq2dmin.$(ARCH).o \