Index: trunk/psLib/test/math/tst_psMinimizeLMM.c =================================================================== --- trunk/psLib/test/math/tst_psMinimizeLMM.c (revision 7085) +++ trunk/psLib/test/math/tst_psMinimizeLMM.c (revision 7130) @@ -12,12 +12,35 @@ #include #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,7 +49,7 @@ *****************************************************************************/ -psF32 myFunc(psVector *deriv, - psVector *params, - psVector *x) +psF32 fitFunc(psVector *deriv, + psVector *params, + psVector *x) { if ((deriv == NULL) || (params == NULL) || (x == NULL)) { @@ -34,17 +57,6 @@ } - 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,58 +66,70 @@ 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;idata[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;idata.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;idata[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;idata.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);