Changeset 7130

Timestamp:

May 17, 2006, 2:24:13 PM (20 years ago)

Author:

Paul Price

Message:

Better test for psMinimize --- fitting a Gaussian, and comparing parameters.

File:

• trunk/psLib/test/math/tst_psMinimizeLMM.c (modified) (4 diffs)

trunk/psLib/test/math/tst_psMinimizeLMM.c

@@ -12 +12 @@
 #include <math.h>
 #define NUM_ITERATIONS 100
-#define ERR_TOL 0.1
-#define NUM_DATA_POINTS 20
-#define NUM_PARAMS 5
+#define ERR_TOL 1e-6
+#define NUM_DATA_POINTS 100
+#define NUM_PARAMS 3
 float expectedParm[NUM_PARAMS];
+
+
+float function(const psVector *params,      // Paramters
+               const psVector *x            // Ordinate
+              )
+{
+    return params->data.F32[0] *
+           expf(- params->data.F32[1] * (PS_SQR(x->data.F32[0]) + PS_SQR(x->data.F32[1]))) +
+           params->data.F32[2];
+}
+
+void derivatives(psVector *deriv,       // Derivatives
+                 const psVector *params,      // Paramters
+                 const psVector *x            // Ordinate
+                )
+{
+    deriv->data.F32[0] = expf(- params->data.F32[1] * (PS_SQR(x->data.F32[0]) + PS_SQR(x->data.F32[1])));
+    deriv->data.F32[1] = - params->data.F32[0] * (PS_SQR(x->data.F32[0]) + PS_SQR(x->data.F32[1])) *
+                         expf(- params->data.F32[1] * (PS_SQR(x->data.F32[0]) + PS_SQR(x->data.F32[1])));
+    deriv->data.F32[2] = 1.0;
+}
+
+
 
 /*****************************************************************************
 myFunc():
-    sum = param[0] + x[0] * x[1] +
+    sum = param[0] * x[0] * x[1] +
           param[1] * x[0] +
           param[2] * x[0]^2 +
@@ -26 +49 @@
  
  *****************************************************************************/
-psF32 myFunc(psVector *deriv,
-             psVector *params,
-             psVector *x)
+psF32 fitFunc(psVector *deriv,
+              psVector *params,
+              psVector *x)
 {
     if ((deriv == NULL) || (params == NULL) || (x == NULL)) {
@@ -34 +57 @@
     }
 
-    psF32 sum = (params->data.F32[0] * x->data.F32[0] * x->data.F32[1]) +
-                (params->data.F32[1] * x->data.F32[0]) +
-                (params->data.F32[2] * PS_SQR(x->data.F32[0])) +
-                (params->data.F32[3] * x->data.F32[1]) +
-                (params->data.F32[4] * PS_SQR(x->data.F32[1]));
-
-    deriv->data.F32[0] = x->data.F32[0] * x->data.F32[1];
-    deriv->data.F32[1] = x->data.F32[0];
-    deriv->data.F32[2] = PS_SQR(x->data.F32[0]);
-    deriv->data.F32[3] = x->data.F32[1];
-    deriv->data.F32[4] = PS_SQR(x->data.F32[1]);
-
-    return(sum);
+    derivatives(deriv, params, x);
+    return function(params, x);
 }
 
@@ -54 +66 @@
     psBool testStatus = true;
 
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 1); // Random number generator; using known seed
     psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL);
-    psImage *myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);
-    psVector *myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
-    psVector *myDerivs = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
-    psArray *myCoords = psArrayAlloc(NUM_DATA_POINTS);
-    myCoords->n = NUM_DATA_POINTS;
-    psVector *y = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32);
-    myParams->n = NUM_PARAMS;
-    myDerivs->n = NUM_PARAMS;
-    y->n = NUM_DATA_POINTS;
-
-    for (psS32 i=0;i<NUM_DATA_POINTS;i++) {
-        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
-        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) (i+10);
-        ((psVector *) (myCoords->data[i]))->data.F32[1] = (float) (i+3);
-        ((psVector *) (myCoords->data[i]))->n = 2;
-        psTrace(__func__, 5, "x[%d] is vector (%f, %f)\n", i, ((psVector *) (myCoords->data[i]))->data.F32[0],
-                ((psVector *) (myCoords->data[i]))->data.F32[1]);
-        y->data.F32[i] = (float) i;
-    }
-    for (psS32 i=0;i<NUM_PARAMS;i++) {
-        expectedParm[i] = 2.42 + (float) (2 * i);
-        myParams->data.F32[i] = (float) (2 * i);
-    }
-
-    psBool rc = psMinimizeLMChi2(min, NULL, myParams, NULL, myCoords, y, NULL,
-                                 (psMinimizeLMChi2Func) myFunc);
-    if (rc == false) {
+    psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS);
+    psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32);
+    psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32);
+    ordinates->n = NUM_DATA_POINTS;
+    coordinates->n = NUM_DATA_POINTS;
+    errors->n = NUM_DATA_POINTS;
+    psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    params->n = NUM_PARAMS;
+    psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    trueParams->n = NUM_PARAMS;
+
+    trueParams->data.F32[0] = 100.0;    // Normalisation
+    trueParams->data.F32[1] = 0.1;      // Shape
+    trueParams->data.F32[2] = 10.0;     // Background
+
+    // Set parameters
+    for (long i = 0; i < NUM_PARAMS; i++) {
+        // So we're not starting right on the true value:
+        params->data.F32[i] = trueParams->data.F32[i] * (1.0 - psRandomGaussian(rng) / 10.0);
+    }
+
+    for (long i = 0; i < NUM_DATA_POINTS; i++) {
+        psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+        x->data.F32[0] = 10.0 * psRandomUniform(rng) - 5.0;
+        x->data.F32[1] = 10.0 * psRandomUniform(rng) - 5.0;
+        ordinates->data[i] = x;
+        coordinates->data.F32[i] = function(params, x) + 0.1 * psRandomGaussian(rng); // Add some noise
+        errors->data.F32[i] = 0.1;
+        printf("Data %ld: (%f, %f) --> %f\n", i, x->data.F32[0], x->data.F32[1], coordinates->data.F32[i]);
+    }
+
+    if (!psMinimizeLMChi2(min, NULL, params, NULL, ordinates, coordinates, errors,
+                          (psMinimizeLMChi2Func)fitFunc)) {
         printf("TEST ERROR: psMinimizeLMChi2() returned FALSE.\n");
         fflush(stdout);
         testStatus = false;
     } else {
-        printf("\nThe value of the function at the minimum is %f\n", min->value);
-        for (psS32 i=0;i<NUM_DATA_POINTS;i++) {
-            float myFuncValue = (float)(myFunc(myDerivs, myParams, (psVector *) myCoords->data[i]));
-            myFuncValue += 0.001;
-            printf("The minimum for data point number %d: f is %.2f\n",
-                   i, myFuncValue);
+        printf("Minimisation took %d iterations\n", min->iter);
+        printf("chi^2 at the minimum is %f\n", min->value);
+        for (long i = 0; i < NUM_PARAMS; i++) {
+            printf("Parameter %ld at the minimum is %f, expected: %f\n", i,
+                   params->data.F32[i], trueParams->data.F32[i]);
             fflush(stdout);
         }
-
-        for (psS32 i=0;i<NUM_PARAMS;i++) {
-            printf("Parameter %d at the minimum is %.1f (expected: %.1f)\n", i,
-                   myParams->data.F32[i], expectedParm[i]);
+        float diff = 0.0;
+        for (long i = 0; i < NUM_DATA_POINTS; i++) {
+            psVector *x = ordinates->data[i];
+            float fitted = function(trueParams, x);
+            float expected = function(params, x);
+            diff += (fitted - expected) / fabsf(expected);
+            printf("Data point %ld: Fitted: %f, expected: %f\n", i, fitted, expected);
             fflush(stdout);
         }
+        printf("Mean relative difference is %f\n", diff/(float)NUM_DATA_POINTS);
+
     }
 
     psFree(min);
-    psFree(myCovar);
-    psFree(myParams);
-    psFree(myDerivs);
-    psFree(myCoords);
-    psFree(y);
+    psFree(params);
+    psFree(trueParams);
+    psFree(ordinates);
+    psFree(coordinates);
+    psFree(errors);
+    psFree(rng);
     psMemCheckCorruption(1);
     psS32 memLeaks = psMemCheckLeaks(currentId, NULL, NULL, false);