Changeset 42497

Timestamp:

Aug 15, 2023, 12:16:20 PM (3 years ago)

Author:

eugene

Message:

working on the polynomial-fitting tests

Location:

branches/eam_branches/ipp-20230313/psLib/test/math

Files:

• tap_psPolyFit1D.c (modified) (2 diffs)
• tap_psPolyFit_IRLS.c (modified) (2 diffs)

branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit1D.c

@@ -8 +8 @@
  
 XXX: Try null stats.
+
+XXX: this function is confused about POLY_ORDER (5 should have x^5, but setData stops at x^4)
  *****************************************************************************/
 #include <stdio.h>
@@ -350 +352 @@
     plan_tests(104);
 
+    // psTraceSetLevel ("psLib.math.psMatrixGJSolve", 4);
 
     // psVectorFitPolynomial1D()

branches/eam_branches/ipp-20230313/psLib/test/math/tap_psPolyFit_IRLS.c

@@ -5 +5 @@
 #include "pstap.h"
 
-#define POLY_ORDER 5
+// #define POLY_ORDER 5
 // #define A -560.0
 // #define B +1116.0
@@ -20 +20 @@
 // #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);
 }