Index: /branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit1D.c =================================================================== --- /branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit1D.c (revision 42496) +++ /branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit1D.c (revision 42497) @@ -8,4 +8,6 @@ XXX: Try null stats. + +XXX: this function is confused about POLY_ORDER (5 should have x^5, but setData stops at x^4) *****************************************************************************/ #include @@ -350,4 +352,5 @@ plan_tests(104); + // psTraceSetLevel ("psLib.math.psMatrixGJSolve", 4); // psVectorFitPolynomial1D() Index: /branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit_IRLS.c =================================================================== --- /branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit_IRLS.c (revision 42496) +++ /branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit_IRLS.c (revision 42497) @@ -5,5 +5,5 @@ #include "pstap.h" -#define POLY_ORDER 5 +// #define POLY_ORDER 5 // #define A -560.0 // #define B +1116.0 @@ -20,68 +20,130 @@ // #define F +2.0 -#define A -300.0 -#define B +1200.0 -#define C -400.0 -#define D +500.0 -#define E -100.0 -#define F +3.0 +// #define A -300.0 +// #define B +1200.0 +// #define C -400.0 +// #define D +500.0 +// #define E -100.0 +// #define F +3.0 + +#define POLY_ORDER 1 +#define A -1.0 +#define B +2.0 #define MASK_VALUE 1 +static psRandom *rng = NULL; + psF32 setData(psF32 x) { - return(A + (B * x) + (C * x * x) + (D * x * x * x) + (E * x * x * x * x) + (F * x * x * x * x * x)); +// return(A + (B * x) + (C * x * x) + (D * x * x * x) + (E * x * x * x * x) + (F * x * x * x * x * x)); + return(A + (B * x)); +} + +# define DUMP_FIT \ + /* compare the fitted values to the input coefficients */ \ + fprintf (stderr, "%5.3f vs %5.3f : %5.2f\n", myPoly->coeff[0], A, myPoly->coeff[0] - A); \ + fprintf (stderr, "%5.3f vs %5.3f : %5.2f\n", myPoly->coeff[1], B, myPoly->coeff[1] - B); + +// fprintf (stderr, "%5.2f vs %5.2f : %5.2f = %5.2f\n", myPoly->coeff[2], C, myPoly->coeff[2] - C, (myPoly->coeff[2] - C)/myPoly->coeff[2]); +// fprintf (stderr, "%5.2f vs %5.2f : %5.2f = %5.2f\n", myPoly->coeff[3], D, myPoly->coeff[3] - D, (myPoly->coeff[3] - D)/myPoly->coeff[3]); +// fprintf (stderr, "%5.2f vs %5.2f : %5.2f = %5.2f\n", myPoly->coeff[4], E, myPoly->coeff[4] - E, (myPoly->coeff[4] - E)/myPoly->coeff[4]); +// fprintf (stderr, "%5.2f vs %5.2f : %5.2f = %5.2f\n", myPoly->coeff[5], F, myPoly->coeff[5] - F, (myPoly->coeff[5] - F)/myPoly->coeff[5]); + +# define IS_IRLS 1 +int genericFit (psVector *xTruth, psVector *fTruth, float sigma, float outfrac, int isIRLS) { + + psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, POLY_ORDER); + psVector *x = psVectorAlloc(xTruth->n, PS_TYPE_F64); + psVector *f = psVectorAlloc(xTruth->n, PS_TYPE_F64); + psVector *fErr = psVectorAlloc(xTruth->n, PS_TYPE_F64); + // psVector *mask = psVectorAlloc(numData, PS_TYPE_U8); + + // set values with an error of 0.5 + for (int i = 0; i < xTruth->n; i++) { + x->data.F64[i] = xTruth->data.F64[i]; + f->data.F64[i] = fTruth->data.F64[i] + sigma*psRandomGaussian(rng); + fErr->data.F64[i] = sigma; + } + + // add in 5% random outliers + for (int i = 0; i < outfrac*xTruth->n; i++) { + int n = (xTruth->n - 1)*psRandomUniform(rng); // n is in range 0 - (numData-1) + float nsig = 20.0*psRandomUniform(rng) +10.0; + // symmetric outliers do not change the fit parameters + // nsig = (psRandomUniform(rng) > 0.5) ? nsig : -1.0 * nsig; + f->data.F64[n] += nsig*sigma; + } + psMemId id = psMemGetId(); + psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); + stats->clipIter = 10; // max number of iterations + + if (isIRLS) { + diag ("IRLS fits with %f outliers: sigma = %f", outfrac, sigma); + bool rc = psVectorIRLSFitPolynomial1D(myPoly, stats, NULL, MASK_VALUE, f, fErr, x); + ok(rc == true, "fit succeeded mechanically"); + + // XXX test: + FILE *ftest = fopen ("irls.ft.dat", "w"); + for (int i = 0; i < f->n; i++) { + fprintf (ftest, "%d %f %f %f\n", i, x->data.F64[i], f->data.F64[i], fErr->data.F64[i]); + } + fclose (ftest); + } else { + diag ("ORD fits with %f outliers: sigma = %f", outfrac, sigma); + bool rc = psVectorFitPolynomial1D(myPoly, NULL, MASK_VALUE, f, fErr, x); + ok(rc == true, "fit succeeded mechanically"); + } + + // compare the fitted values to the input coefficients + DUMP_FIT; + + // is_float_tol (myPoly->coeff[0], A, 1e-5, "A coeffs match"); + // is_float_tol (myPoly->coeff[1], B, 1e-5, "B coeffs match"); + // is_float_tol (myPoly->coeff[2], C, 1e-5, "C coeffs match"); + // is_float_tol (myPoly->coeff[3], D, 1e-5, "D coeffs match"); + // is_float_tol (myPoly->coeff[4], E, 1e-5, "E coeffs match"); + // is_float_tol (myPoly->coeff[5], F, 1e-5, "F coeffs match"); + + psFree (stats); + psFree (x); + psFree (f); + psFree (fErr); + psFree (myPoly); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + + return true; } int main() { - psLogSetFormat("HLNM"); - psLogSetLevel(PS_LOG_INFO); - plan_tests(104); + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + plan_tests(104); - // create x vector running from -1 to 11 in steps of 0.1 - int numData = 120; - psVector *xTruth = psVectorAlloc(numData, PS_TYPE_F64); - psVector *fTruth = psVectorAlloc(numData, PS_TYPE_F64); - psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Using a known seed + // create x vector running from -1 to 11 in steps of 0.1 + int numData = 120; + psVector *xTruth = psVectorAlloc(numData, PS_TYPE_F64); + psVector *fTruth = psVectorAlloc(numData, PS_TYPE_F64); + rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Using a known seed + + // create reference data + for (int i = 0; i < numData; i++) { + xTruth->data.F64[i] = 0.1*i - 1.0; + fTruth->data.F64[i] = setData(xTruth->data.F64[i]); + } - psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, POLY_ORDER); - psVector *x = psVectorAlloc(numData, PS_TYPE_F64); - psVector *f = psVectorAlloc(numData, PS_TYPE_F64); - psVector *fErr = psVectorAlloc(numData, PS_TYPE_F64); - // psVector *mask = psVectorAlloc(numData, PS_TYPE_U8); - - // XXX F64 or F32?? - for (int i = 0; i < numData; i++) { - xTruth->data.F64[i] = 0.1*i - 1.0; - fTruth->data.F64[i] = setData(xTruth->data.F64[i]); - } + genericFit (xTruth, fTruth, 0.1, 0.20, !IS_IRLS); + genericFit (xTruth, fTruth, 0.1, 0.20, IS_IRLS); - // set values with an error of sigma - double sigma = 0.01; - for (int i = 0; i < numData; i++) { - x->data.F64[i] = xTruth->data.F64[i]; - f->data.F64[i] = fTruth->data.F64[i] + sigma*psRandomGaussian(rng); - fErr->data.F64[i] = sigma; - // fprintf (stdout, "%d %f %f %f\n", i, x->data.F64[i], f->data.F64[i], fErr->data.F64[i]); - } - - { - psMemId id = psMemGetId(); - bool rc = psVectorFitPolynomial1D(myPoly, NULL, MASK_VALUE, f, fErr, x); - ok(rc == true, "fit succeeded mechanically"); - - // compare the fitted values to the input coefficients - - fprintf (stderr, "%5.2f vs %5.2f\n", myPoly->coeff[0], A); - fprintf (stderr, "%5.2f vs %5.2f\n", myPoly->coeff[1], B); - fprintf (stderr, "%5.2f vs %5.2f\n", myPoly->coeff[2], C); - fprintf (stderr, "%5.2f vs %5.2f\n", myPoly->coeff[3], D); - fprintf (stderr, "%5.2f vs %5.2f\n", myPoly->coeff[4], E); - fprintf (stderr, "%5.2f vs %5.2f\n", myPoly->coeff[5], F); - - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - + // genericFit (xTruth, fTruth, 0.01, 0.00, !IS_IRLS); + // genericFit (xTruth, fTruth, 0.50, 0.00, !IS_IRLS); + // genericFit (xTruth, fTruth, 0.01, 0.05, !IS_IRLS); + // genericFit (xTruth, fTruth, 0.50, 0.05, !IS_IRLS); + // + // genericFit (xTruth, fTruth, 0.01, 0.00, IS_IRLS); + // genericFit (xTruth, fTruth, 0.50, 0.00, IS_IRLS); + // genericFit (xTruth, fTruth, 0.01, 0.05, IS_IRLS); + // genericFit (xTruth, fTruth, 0.50, 0.05, IS_IRLS); }