Index: trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.c (revision 35674) +++ trunk/psLib/src/imageops/psImageConvolve.c (revision 35767) @@ -1231,4 +1231,6 @@ } +/*********************** smooth with mask, threaded version ***********************/ + bool psImageSmoothMask_ScanRows_F32(psImage *calculation, psImage *calcMask, const psImage *image, const psImage *mask, psImageMaskType maskVal, psVector *gaussNorm, @@ -1484,4 +1486,224 @@ +/*********************** smooth no-mask ***********************/ + +bool psImageSmoothNoMask_ScanRows_F32(psImage *calculation, const psImage *image, + psVector *gaussNorm, float minGauss, + int size, int rowStart, int rowStop) +{ + const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel + + int numCols = image->numCols; + int xLast = numCols - 1; // Last index + + for (int j = rowStart; j < rowStop; j++) { + for (int i = 0; i < numCols; i++) { + int xMin = PS_MAX(i - size, 0); + int xMax = PS_MIN(i + size, xLast); + const psF32 *imageData = &image->data.F32[j][xMin]; + int uMin = - PS_MIN(i, size); /* Minimum kernel index */ + const psF32 *gaussData = &gauss[uMin]; + double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */ + for (int x = xMin; x <= xMax; x++, imageData++, gaussData++) { + sumIG += *imageData * *gaussData; + sumG += *gaussData; + } + if (sumG > minGauss) { + /* BW */ + calculation->data.F32[i][j] = sumIG / sumG; + } + } + } + return true; +} + +bool psImageSmoothNoMask_ScanCols_F32(psImage *output, psImage *calculation, + psVector *gaussNorm, float minGauss, + int size, int colStart, int colStop) +{ + const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel + + int numRows = output->numRows; + int yLast = numRows - 1; // Last index + + for (int i = colStart; i < colStop; i++) { + for (int j = 0; j < numRows; j++) { + int yMin = PS_MAX(j - size, 0); + int yMax = PS_MIN(j + size, yLast); + const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */ + int vMin = - PS_MIN(j, size); /* Minimum kernel index */ + const psF32 *gaussData = &gauss[vMin]; + double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */ + for (int y = yMin; y <= yMax; y++, imageData++, gaussData++) { + sumIG += *imageData * *gaussData; + sumG += *gaussData; + } + output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN; + } + } + return true; +} + +bool psImageSmoothNoMask_ScanRows_F32_Threaded(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psAssert (job->args, "programming error: job args not alloced?"); + psAssert (job->args->n == 7, "programming error: job Nargs mismatch"); + + psImage *calculation = job->args->data[0]; // calculation image + const psImage *image = job->args->data[1]; // input image + + psVector *gaussNorm = job->args->data[2]; // gauss kernel + float minGauss = PS_SCALAR_VALUE(job->args->data[3],F32); + int size = PS_SCALAR_VALUE(job->args->data[4],S32); + int rowStart = PS_SCALAR_VALUE(job->args->data[5],S32); + int rowStop = PS_SCALAR_VALUE(job->args->data[6],S32); + return psImageSmoothNoMask_ScanRows_F32(calculation, image, + gaussNorm, minGauss, size, + rowStart, rowStop); +} + +bool psImageSmoothNoMask_ScanCols_F32_Threaded(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + psAssert (job->args, "programming error: job args not alloced?"); + psAssert (job->args->n == 7, "programming error: job Nargs mismatch"); + + psImage *output = job->args->data[0]; // input image + psImage *calculation = job->args->data[1]; // calculation image + + psVector *gaussNorm = job->args->data[2]; // gauss kernel + float minGauss = PS_SCALAR_VALUE(job->args->data[3],F32); + int size = PS_SCALAR_VALUE(job->args->data[4],S32); + int colStart = PS_SCALAR_VALUE(job->args->data[5],S32); + int colStop = PS_SCALAR_VALUE(job->args->data[6],S32); + return psImageSmoothNoMask_ScanCols_F32(output, calculation, + gaussNorm, minGauss, size, + colStart, colStop); +} + +psImage *psImageSmoothNoMask_Threaded(psImage *output, const psImage *image, float sigma, float numSigma, float minGauss) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + + int numCols = image->numCols, numRows = image->numRows; // Size of image + + int nThreads = psThreadPoolSize(); + int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols; + int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows; + int size = sigma * numSigma + 0.5; // Half-size of kernel + + // Generate normalized gaussian + psVector *gaussNorm = NULL; + IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32); + + switch (image->type.type) { + // Smooth an image with masked pixels + // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or + // [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW" + + case PS_TYPE_F32: { + psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */ + + if (threaded) { + /** Smooth in X direction **/ + for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) { + int rowStop = PS_MIN (rowStart + scanRows, numRows); + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS"); + psArrayAdd(job->args, 0, calculation); + psArrayAdd(job->args, 0, (psImage *) image); // cast away const + psArrayAdd(job->args, 0, gaussNorm); + PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, rowStop, PS_TYPE_S32); + // -> psImageSmoothNoMask_ScanRows_F32 (calculation, image, gauss, minGauss, size, rowStart, rowStop); + + // if threading is not active, we simply run the job and return + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } + // wait here for the threaded jobs to finish (NOP if threading is not active) + if (!psThreadPoolWait(true, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } else if (!psImageSmoothNoMask_ScanRows_F32(calculation, image, gaussNorm, minGauss, size, 0, numRows)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + + output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32); + + /** Smooth in Y direction **/ + if (threaded) { + for (int colStart = 0; colStart < numCols; colStart+=scanCols) { + int colStop = PS_MIN (colStart + scanCols, numCols); + + // allocate a job, construct the arguments for this job + psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS"); + psArrayAdd(job->args, 0, output); + psArrayAdd(job->args, 0, calculation); + psArrayAdd(job->args, 0, gaussNorm); + PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, size, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, colStop, PS_TYPE_S32); + // -> psImageSmoothMask_ScanCols_F32 (output, calculation, gauss, minGauss, size, colStart, colStop); + + // if threading is not active, we simply run the job and return + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } + + // wait here for the threaded jobs to finish (NOP if threading is not active) + if (!psThreadPoolWait(true, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + } else if (!psImageSmoothNoMask_ScanCols_F32(output, calculation, gaussNorm, minGauss, size, 0, numCols)) { + psError(PS_ERR_UNKNOWN, false, "Unable to smooth image"); + psFree(calculation); + psFree(gaussNorm); + return false; + } + + psFree(calculation); + psFree(gaussNorm); + + return output; + } + default: + psFree(gaussNorm); + char *typeStr; + PS_TYPE_NAME(typeStr,image->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + _("Specified psImage type, %s, is not supported."), + typeStr); + return false; + } + + psAbort("Should never reach here."); + return false; +} + +/****************************************************************************************/ + bool psImageSmoothMaskF32(psImage *image, psImage *mask, psImageMaskType maskVal, double sigma, double Nsigma) @@ -1900,4 +2122,16 @@ } { + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS", 7); + task->function = &psImageSmoothNoMask_ScanRows_F32_Threaded; + psThreadTaskAdd(task); + psFree(task); + } + { + psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS", 7); + task->function = &psImageSmoothNoMask_ScanCols_F32_Threaded; + psThreadTaskAdd(task); + psFree(task); + } + { psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9); task->function = &psImageSmoothMaskPixelsThread; @@ -1909,4 +2143,6 @@ psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANROWS"); psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS"); + psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS"); + psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS"); psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_PIXELS"); } Index: trunk/psLib/src/imageops/psImageConvolve.h =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.h (revision 35674) +++ trunk/psLib/src/imageops/psImageConvolve.h (revision 35767) @@ -242,5 +242,6 @@ ); -/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for masked pixels +/// Smooth an imageby parts using 1D Gaussian independently in x and y, allowing for +/// MASKED PIXELS /// /// Applies a circularly symmetric Gaussian smoothing first in x and then in y @@ -256,5 +257,5 @@ ); -/// Smooth particular pixels on an image, allowing for masked pixels +/// Smooth particular pixels on an image, allowing for MASKED PIXELS /// /// Applies a circularly symmetric Gaussian smoothing first in x and then in y @@ -271,4 +272,9 @@ ); +/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for +/// MASKED PIXELS : THREADED VERSION +/// +/// Applies a circularly symmetric Gaussian smoothing first in x and then in y +/// directions with just a vector. This process is 2N faster than 2D convolutions (in general). psImage *psImageSmoothMask_Threaded(psImage *output, const psImage *image, @@ -279,4 +285,17 @@ float minGauss); +/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for +/// MASKED PIXELS : THREADED VERSION +/// +/// Applies a circularly symmetric Gaussian smoothing first in x and then in y +/// directions with just a vector. This process is 2N faster than 2D convolutions (in general). +psImage *psImageSmoothNoMask_Threaded(psImage *output, + const psImage *image, + float sigma, + float numSigma, + float minGauss); + +/// Smooth an image by parts IN-SITU using 1D Gaussian independently in x and y, allowing +/// for MASKED PIXELS bool psImageSmoothMaskF32( psImage *image, ///< the image to be smoothed Index: trunk/psLib/src/math/psMinimizeLMM.c =================================================================== --- trunk/psLib/src/math/psMinimizeLMM.c (revision 35674) +++ trunk/psLib/src/math/psMinimizeLMM.c (revision 35767) @@ -325,4 +325,8 @@ chisq += PS_SQR(delta) * dy->data.F32[i]; + // XXX remove this later: + // psVector *tmp = x->data[i]; + // fprintf (stderr, "%f %f %f %f %f\n", tmp->data.F32[0], tmp->data.F32[1], y->data.F32[i], dy->data.F32[i], ymodel); + if (isnan(dy->data.F32[i])) goto escape; if (isnan(delta)) goto escape; @@ -360,4 +364,28 @@ } +/****************************************************************************** +psMinimizeLMChi2(): wrapper to call either _Old or _Alt + *****************************************************************************/ +bool psMinimizeLMChi2( + psMinimization *min, + psImage *covar, + psVector *params, + psMinConstraint *constraint, + const psArray *x, + const psVector *y, + const psVector *yWt, + psMinimizeLMChi2Func func) +{ + bool status = psMinimizeLMChi2_Alt( + min, + covar, + params, + constraint, + x, + y, + yWt, + func); + return status; +} /****************************************************************************** @@ -370,5 +398,5 @@ XXX Make an F64 version? *****************************************************************************/ -bool psMinimizeLMChi2( +bool psMinimizeLMChi2_Old( psMinimization *min, psImage *covar, @@ -507,7 +535,7 @@ psF32 rho = (min->value - Chisq) / dLinear; - psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); - - psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); + psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); + + psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); // dump some useful info if trace is defined @@ -580,4 +608,317 @@ } +/****************************************************************************** +psMinimizeLMChi2(): This routine takes a function-pointer (func) which calculates an arbitrary +function and it's derivatives and minimizes the chi-squared match between that function at the +specified points and the specified value at those points. + +The original version of this function used a convergence criterion based on the change in +chisq. this has problems since it depends on the choice of points used to measure the fit. +(consider a gaussian on a background : it 100 pixels are used -- and some or most contribute to +the chisq -- and the delta-chisq is 10%, then the same change in model fit will yield a +delta-chisq of 1% if 1000 pixels are used (all but the 100 measuring the background)). + +This implementation uses changes to the parameters and stops if they are no longer significant. + +This requires F32 input data; all internal calls use F32. + *****************************************************************************/ +bool psMinimizeLMChi2_Alt( + psMinimization *min, + psImage *covar, + psVector *params, + psMinConstraint *constraint, + const psArray *x, + const psVector *y, + const psVector *yWt, + psMinimizeLMChi2Func func) +{ + psTrace("psLib.math", 3, "---- begin ----\n"); + PS_ASSERT_PTR_NON_NULL(min, false); + PS_ASSERT_VECTOR_NON_NULL(params, false); + PS_ASSERT_VECTOR_NON_EMPTY(params, false); + PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); + psVector *paramMask = NULL; + if (constraint != NULL) { + paramMask = constraint->paramMask; + if (paramMask != NULL) { + PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false); + } + } + PS_ASSERT_PTR_NON_NULL(x, false); + for (psS32 i = 0 ; i < x->n ; i++) { + psVector *coord = (psVector *) (x->data[i]); + PS_ASSERT_VECTOR_NON_NULL(coord, false); + PS_ASSERT_VECTOR_TYPE(coord, PS_TYPE_F32, false); + } + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_NON_EMPTY(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, false); + if (yWt != NULL) { + PS_ASSERT_VECTOR_TYPE(yWt, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, yWt, false); + } + PS_ASSERT_PTR_NON_NULL(func, false); + + psMinimizeLMLimitFunc checkLimits = NULL; + if (constraint) { + checkLimits = constraint->checkLimits; + } + + // this function has test and current values for several things + // the current best value is in lower case + // the next guess value is in upper case + + // allocate internal arrays (current vs Guess) + psImage *Alpha = NULL; + psVector *Beta = NULL; + + // Alpha & Beta only contain elements to represent the unmasked parameters + if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) { + psAbort ("programming error: no unmasked parameters to be fit\n"); + } + + int nFitParams = Beta->n; + psImage *alpha = psImageAlloc(nFitParams, nFitParams, PS_TYPE_F32); + psVector *beta = psVectorAlloc(nFitParams, PS_TYPE_F32); + psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); + + psVector *dy = NULL; + psF32 Chisq = 0.0; + psF32 lambda = 0.001; + psF32 dLinear = 0.0; + psF32 nu = 2.0; + + // the user provides the error or NULL. we need to convert + // to appropriate weights + if (yWt != NULL) { + dy = (psVector *) yWt; + } else { + dy = psVectorAlloc(y->n, PS_TYPE_F32); + psVectorInit(dy, 1.0); + } + + // number of degrees of freedom for this fit + int nDOF = dy->n - nFitParams; + + // calculate initial alpha and beta, set chisq (min->value) + min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func); + if (isnan(min->value)) { + min->iter = min->maxIter; + psFree(alpha); + psFree(Alpha); + psFree(beta); + psFree(Beta); + psFree(Params); + return(false); + } + // dump some useful info if trace is defined + if (psTraceGetLevel("psLib.math") >= 6) { + p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)"); + p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)"); + } + if (psTraceGetLevel("psLib.math") >= 5) { + p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)"); + } + + bool done = (min->iter >= min->maxIter); + while (!done) { + psTrace("psLib.math", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); + + if (min->chisqConvergence) { + psTrace("psLib.math", 5, "Last delta is %f. stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol); + } else { + psTrace("psLib.math", 5, "Last delta is %f. stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*nFitParams, min->maxTol*nFitParams); + } + + // set a new guess for Alpha, Beta, Params + if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) { + min->iter ++; + if (min->iter >= min->maxIter) break; + lambda *= 10.0; + // ALT? lambda *= 2.0; + continue; + } + + // dump some useful info if trace is defined + if (psTraceGetLevel("psLib.math") >= 6) { + p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)"); + p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)"); + p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)"); + } + if (psTraceGetLevel("psLib.math") >= 5) { + p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)"); + } + + // calculate the parameter change (rParDelta) and error radius (rParSigma) + // rParDelta : radius of parameter change; + // rParSigma : radius of parameter error + + // note that (before SetABX) Alpha[i][i] is the covariance matrix and + // Beta is the actual parameter change for this pass + + // note that Alpha & Beta only represent unmasked parameters, while params and Params have all + + // dParSigma = Alpha[i][i] : error (squared) on parameter i + // dParDelta = Params->data.F32[i] - params->data.F32[i] : change on parameter i + float rParSigma = 0.0; + for (int j = 0, J = 0; j < Params->n; j++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) { + continue; + } + rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J]; + J++; + } + rParSigma = sqrt(rParSigma); + psTrace("psLib.math", 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter); + // fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter); + + // calculate Chisq for new guess, update Alpha & Beta + Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func); + if (isnan(Chisq)) { + min->iter ++; + if (min->iter >= min->maxIter) break; + lambda *= 10.0; + // ALT lambda *= 2.0; + continue; + } + + // convergence criterion: + // compare the delta (min->value - Chisq) with the + // expected delta from the linear model (dLinear) + // accept new guess if it is an improvement (rho > 0), or else increase lambda + psF32 rho = (min->value - Chisq) / dLinear; + + psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); + + psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); + + // dump some useful info if trace is defined + if (psTraceGetLevel("psLib.math") >= 6) { + p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)"); + p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)"); + } + + // change in chisq/nDOF since last minimum + min->lastDelta = (min->value - Chisq) / nDOF; + + // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) + + // XXX the old version of lambda changes: + // XXX : Madsen gives suggestion for better use of rho + // rho is positive if the new chisq is smaller + if (rho >= -1e-6) { + min->value = Chisq; + alpha = psImageCopy(alpha, Alpha, PS_TYPE_F32); + beta = psVectorCopy(beta, Beta, PS_TYPE_F32); + params = psVectorCopy(params, Params, PS_TYPE_F32); + } + switch (min->gainFactorMode) { + case 0: + if (rho >= -1e-6) { + lambda *= 0.25; + } else { + lambda *= 10.0; + } + break; + + case 1: + // adjust the gain ratio (lambda) based on rho + if (rho < 0.25) { + lambda *= 2.0; + } + if (rho > 0.75) { + lambda *= 0.333; + } + break; + + case 2: + if (rho > 0.0) { + lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0))); + nu = 2.0; + } else { + lambda *= nu; + nu *= 2.0; + } + break; + } + min->iter++; + + // ending conditions: + // 1) hard limit : too many iterations + done = (min->iter >= min->maxIter); + + // 2) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change) + if (min->lastDelta < -1e-6) { + continue; + } + + // save this value in case we stop iterating + min->rParSigma = rParSigma; + + // 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN) + // keep iterating regardless of rParSigma in this case + float chisqDOF = Chisq / nDOF; + if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) { + continue; + } + + // delta-chisq or rParSigma ? + if (min->chisqConvergence) { + done |= (min->lastDelta < min->minTol); + } else { + done |= (rParSigma < min->minTol*nFitParams); + } + } + psTrace("psLib.math", 5, "chisq: %f, last delta: %f, rParSigma: %f, Niter: %d\n", min->value, min->lastDelta, min->rParSigma, min->iter); + + // construct & return the covariance matrix (if requested) + if (covar != NULL) { + if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { + psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n"); + } + // set covar values which are not masked + psImageInit (covar, 0.0); + for (int j = 0, J = 0; j < params->n; j++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) { + covar->data.F32[j][j] = 1.0; + continue; + } + for (int k = 0, K = 0; k < params->n; k++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue; + covar->data.F32[j][k] = Alpha->data.F32[J][K]; + K++; + } + J++; + } + } + + // free the internal temporary data + psFree(alpha); + psFree(Alpha); + psFree(beta); + psFree(Beta); + psFree(Params); + if (yWt == NULL) { + psFree(dy); + } + + // if the last improvement was at least as good as maxTol, accept the fit: + if (min->chisqConvergence) { + if (min->lastDelta <= min->maxTol) { + psTrace("psLib.math", 6, "---- end (true) ----\n"); + return(true); + } + } else { + if (min->rParSigma <= min->maxTol*nFitParams) { + psTrace("psLib.math", 6, "---- end (true) ----\n"); + return(true); + } + } + psTrace("psLib.math", 6, "---- end (false) ----\n"); + return(false); +} + bool psMinLM_AllocAB (psImage **Alpha, psVector **Beta, const psVector *params, const psVector *paramMask) { @@ -625,6 +966,11 @@ min->iter = 0; min->lastDelta = NAN; + min->rParSigma = NAN; min->maxChisqDOF = NAN; + // we default to the old algorithm for convergence + min->chisqConvergence = true; + min->gainFactorMode = 0; + return(min); } Index: trunk/psLib/src/math/psMinimizeLMM.h =================================================================== --- trunk/psLib/src/math/psMinimizeLMM.h (revision 35674) +++ trunk/psLib/src/math/psMinimizeLMM.h (revision 35767) @@ -30,4 +30,7 @@ #include "psConstants.h" +# define PS_MINIMIZE_LMM_GAIN_FACTOR_MODE 0 +# define PS_MINIMIZE_LMM_CHISQ_CONVERGENCE 1 + #define PS_DETERMINE_BRACKET_STEP_SIZE 0.10 #define PS_MAX_LMM_ITERATIONS 100 @@ -82,5 +85,8 @@ int iter; ///< Number of iterations to date float lastDelta; ///< The last difference for the fit + float rParSigma; ///< last fractional change in the parameters float maxChisqDOF; ///< for Chisq minimization, require that we reach here before checking tolerance + int gainFactorMode; + bool chisqConvergence; } psMinimization; @@ -124,4 +130,35 @@ */ bool psMinimizeLMChi2( + psMinimization *min, ///< Minimization specification + psImage *covar, ///< Covariance matrix + psVector *params, ///< "Best Guess" for the parameters that minimize func + psMinConstraint *constraint, ///< Constraints on the parameters + const psArray *x, ///< Measurement ordinates of multiple vectors + const psVector *y, ///< Measurement coordinates + const psVector *yWt, ///< Errors in the measurement coordinates + psMinimizeLMChi2Func func ///< Specified function +); + +/** Minimizes a specified function based on the Levenberg-Marquardt method. + * + * @return bool: True if successful. + */ +bool psMinimizeLMChi2_Old( + psMinimization *min, ///< Minimization specification + psImage *covar, ///< Covariance matrix + psVector *params, ///< "Best Guess" for the parameters that minimize func + psMinConstraint *constraint, ///< Constraints on the parameters + const psArray *x, ///< Measurement ordinates of multiple vectors + const psVector *y, ///< Measurement coordinates + const psVector *yWt, ///< Errors in the measurement coordinates + psMinimizeLMChi2Func func ///< Specified function +); + +/** Minimizes a specified function based on the Levenberg-Marquardt method + Uses alternative convergence criterion. + * + * @return bool: True if successful. + */ +bool psMinimizeLMChi2_Alt( psMinimization *min, ///< Minimization specification psImage *covar, ///< Covariance matrix Index: trunk/psLib/test/imageops/Makefile.am =================================================================== --- trunk/psLib/test/imageops/Makefile.am (revision 35674) +++ trunk/psLib/test/imageops/Makefile.am (revision 35767) @@ -17,4 +17,5 @@ tap_psImagePixelManip \ tap_psImageSmooth \ + tap_psImageSmoothMask_Threaded \ tap_psImageSmooth_PreAlloc \ tap_psImageStructManip \ Index: trunk/psLib/test/imageops/tap_psImageSmoothMask_Threaded.c =================================================================== --- trunk/psLib/test/imageops/tap_psImageSmoothMask_Threaded.c (revision 35767) +++ trunk/psLib/test/imageops/tap_psImageSmoothMask_Threaded.c (revision 35767) @@ -0,0 +1,112 @@ +/** @file tap_psImageSmoothMask_Threaded.c + * + * This code tests the psImageSmoothMask_Threaded() routine. + * + * @author EAM + * @date $Date: 2007-02-27 23:56:12 $ + * + * Copyright 2013 IfA, University of Hawaii + */ +#include +#include +#include +#include "tap.h" +#include "pstap.h" + +# define NX 4800 +# define NY 4800 +# define SIGMA 1.5 +# define NSIGMA 3.0 + +bool saveImage (psImage *image, char *filename) { + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + return true; +} + +int main(int argc, char **argv) +{ + bool status; + + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + plan_tests(14); + + psMemId id = psMemGetId(); + + psTimerStart ("timer"); + psImage *image = psImageAlloc(NX, NY, PS_TYPE_F32); + + // Set a checkboard pattern + for (int iy = 0 ; iy < NY ; iy++) { + for (int ix = 0 ; ix < NX ; ix++) { + if ((ix%2) != (iy%2)) { + image->data.F32[iy][ix] = 1.0; + } else { + image->data.F32[iy][ix] = 0.0; + } + } + } + saveImage (image, "test.im.fits"); + fprintf (stderr, "create image : %f sec\n", psTimerMark("timer")); psTimerClear("timer"); + + psImage *mask = psImageAlloc(NX, NY, PS_TYPE_IMAGE_MASK); + psImageInit (mask, 0); + + // mask a corner + for (int iy = 0 ; iy < 0.5*NY ; iy++) { + for (int ix = 0 ; ix < 0.5*NX ; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 1; + } + } + saveImage (mask, "test.mk.fits"); + fprintf (stderr, "create mask : %f sec\n", psTimerMark("timer")); psTimerClear("timer"); + + psImage *alt = psImageCopy(NULL, image, image->type.type); + + psTimerClear("timer"); + status = psImageSmooth(alt, SIGMA, NSIGMA); + ok(status, "basic image smooth"); + fprintf (stderr, "basic smooth : %f sec\n", psTimerMark("timer")); + + saveImage (alt, "test.sm1.fits"); + psFree(alt); + + psImage *alt2 = psImageCopy(NULL, image, image->type.type); + psImage *mask2 = psImageCopy(NULL, mask, mask->type.type); + + psTimerClear("timer"); + status = psImageSmoothMaskF32(alt2, mask2, 0xff, SIGMA, NSIGMA); + ok(status, "basic image smooth"); + fprintf (stderr, "mask smooth : %f sec\n", psTimerMark("timer")); + + saveImage (alt2, "test.sm2.fits"); + psFree(alt2); + psFree(mask2); + + psThreadPoolInit (4); + psImageConvolveSetThreads(true); + + psTimerClear("timer"); + psImage *alt3 = psImageSmoothMask_Threaded(NULL, image, mask, 0xff, SIGMA, NSIGMA, 0.25); + fprintf (stderr, "mask smooth, threaded : %f sec\n", psTimerMark("timer")); + saveImage (alt3, "test.sm3.fits"); + psFree(alt3); + + psTimerClear("timer"); + psImage *alt4 = psImageSmoothNoMask_Threaded(NULL, image, SIGMA, NSIGMA, 0.25); + fprintf (stderr, "basic smooth, threaded : %f sec\n", psTimerMark("timer")); + saveImage (alt4, "test.sm4.fits"); + psFree(alt4); + + psImageConvolveSetThreads(false); + + psFree(image); + psFree(mask); + psTimerStop(); + + psThreadPoolFinalize(); + + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); +} Index: trunk/psLib/test/math/Makefile.am =================================================================== --- trunk/psLib/test/math/Makefile.am (revision 35674) +++ trunk/psLib/test/math/Makefile.am (revision 35767) @@ -48,4 +48,7 @@ tap_psMatrixVectorArithmetic04 \ tap_psMinimizePowell \ + tap_psMinimizeLMM \ + tap_psMinimizeLMM_Alt \ + tap_psMinimizeLMM_Trail \ tap_psSpline1D \ tap_psPolynomialMD \ Index: trunk/psLib/test/math/tap_psMinimizeLMM_Alt.c =================================================================== --- trunk/psLib/test/math/tap_psMinimizeLMM_Alt.c (revision 35767) +++ trunk/psLib/test/math/tap_psMinimizeLMM_Alt.c (revision 35767) @@ -0,0 +1,299 @@ +/***************************************************************************** + This routine must ensure that psMinimizeLM() works correctly. + + XXX: This code needs a lot of additional test case work. The minimization + currently fails and we don't attempt to check the output values. + XXX: Add tests for + covar arg set to non-NULL + constraint set to non-NULL + Set x->vectors to NULL, or use wrong types + yWt (errors) vector set to incorrect size, type. + *****************************************************************************/ +#include +#include +#include +#include "tap.h" +#include "pstap.h" + +#define NUM_ITER 10 +#define MIN_TOL 0.0001 +#define MAX_TOL 1.0 + +#define NUM_ITERATIONS 100 +#define NUM_X_POINTS 20 +#define NUM_Y_POINTS 20 +#define NUM_PARAMS 3 +#define VERBOSE 0 +float expectedParm[NUM_PARAMS]; + +// test with a simple 2D Gaussian: +// y = p2 + p0 * e^( -0.5 * (x0^2 + x1^2) / p1^2) +psF32 fitFunc(psVector *deriv, + psVector *params, + psVector *x) +{ + if (params == NULL) { + psError(PS_ERR_UNKNOWN, true, "params is NULL.\n"); + } + if (x == NULL) { + psError(PS_ERR_UNKNOWN, true, "x is NULL.\n"); + } + + psF32 X = -0.5*(PS_SQR(x->data.F32[0]) + PS_SQR(x->data.F32[1])); + psF32 Z = X / PS_SQR(params->data.F32[1]); + psF32 dZdP1 = -2.0 * Z / params->data.F32[1]; + + psF32 Q = exp(Z); + psF32 R = params->data.F32[0] * Q; + psF32 F = params->data.F32[2] + R; + + if (deriv) { + deriv->data.F32[0] = Q; + deriv->data.F32[1] = R * dZdP1; // dRdPq = R * dZdP1 since R = exp(Z) + deriv->data.F32[2] = 1.0; + } + + return F; +} + +psS32 main(psS32 argc, char* argv[]) +{ + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + plan_tests(35); + + + // Test psMinimizationAlloc() + { + psMemId id = psMemGetId(); + psMinimization *tmp = psMinimizationAlloc(NUM_ITER, MIN_TOL, MAX_TOL); + ok(tmp != NULL, "psMinimizationAlloc() returned non-NULL"); + skip_start(tmp == NULL, 5, "Skipping tests because psMinimizationAlloc() failed"); + ok(tmp->maxIter == NUM_ITER, "psMinimizationAlloc() properly set ->maxIter"); + ok_float(tmp->minTol, MIN_TOL, "psMinimizationAlloc() properly set ->minTol"); + ok_float(tmp->maxTol, MAX_TOL, "psMinimizationAlloc() properly set ->maxTol"); + ok_float(tmp->value, 0.0, "psMinimizationAlloc() properly set ->value"); + ok(tmp->iter == 0, "psMinimizationAlloc() properly set ->iter (%d)", tmp->iter); + ok(isnan(tmp->lastDelta), "psMinimizationAlloc() properly set ->lastDelta (%f)", tmp->lastDelta); + ok(isnan(tmp->maxChisqDOF), "psMinimizationAlloc() properly set ->maxChisqDOF (%f)", tmp->maxChisqDOF); + skip_end(); + psFree(tmp); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test psMinConstraintAlloc() + { + psMemId id = psMemGetId(); + psMinConstraint *tmp = psMinConstraintAlloc(); + ok(tmp, "psMinConstraintAlloc() returned non-NULL"); + skip_start(!tmp, 2, "Skipping tests because psMinConstraintAlloc() failed"); + ok(!tmp->paramMask, "psMinConstraintAlloc() properly set ->paramMask"); + ok(!tmp->checkLimits, "psMinConstraintAlloc() properly set ->checkLimits"); + psFree(tmp); + skip_end(); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test psMinimizeLMChi2(): unallowed input parameters. + { + int NUM_DATA_POINTS = NUM_X_POINTS * NUM_Y_POINTS; + psMemId id = psMemGetId(); + psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Random number generator; using known seed + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, MIN_TOL, MAX_TOL); + psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS); + psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *coordinatesF64 = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F64); + psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + psVector *paramsF64 = psVectorAlloc(NUM_PARAMS, PS_TYPE_F64); + psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + + // Test with psMinimization set to NULL + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(NULL, NULL, params, NULL, ordinates, + coordinates, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL psMinimization"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with params set to NULL + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(min, NULL, NULL, NULL, ordinates, + coordinates, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL params"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with params wrong type + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(min, NULL, paramsF64, NULL, ordinates, + coordinates, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL params"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with ordinates (x) set to NULL + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, NULL, + coordinates, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL ordinates (x)"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with coordinates (y) set to NULL + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, + NULL, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL coordinates (y)"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with coordinates (y) wrong type (F64) + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, + coordinatesF64, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with coordinates (y) wrong type"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with ordinates and coordinates wrong size + { + psMemId id = psMemGetId(); + coordinates->n--; + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, + coordinates, errors, (psMinimizeLMChi2Func)fitFunc); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL psMinimization"); + coordinates->n++; + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test with function set to NULL + { + psMemId id = psMemGetId(); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, + coordinates, errors, NULL); + ok(!tmpBool, "psMinimizeLMChi2() returned FALSE with NULL fit function"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + psFree(min); + psFree(params); + psFree(paramsF64); + psFree(trueParams); + psFree(ordinates); + psFree(coordinates); + psFree(coordinatesF64); + psFree(errors); + psFree(rng); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // Test psMinimizeLMChi2() with legitimate input values + { + int NUM_DATA_POINTS = NUM_X_POINTS * NUM_Y_POINTS; + psMemId id = psMemGetId(); + psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Random number generator; using known seed + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, MIN_TOL, MAX_TOL); + psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS); + psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + bool tmpBool; + trueParams->data.F32[0] = 100.0; // Normalisation + trueParams->data.F32[1] = 3.0; // Width + trueParams->data.F32[2] = 10.0; // Background + + // Set parameters + for (long i = 0; i < NUM_PARAMS; i++) + { + // Ensure we're not starting right on the true value: + params->data.F32[i] = trueParams->data.F32[i] * + (1.0 - psRandomGaussian(rng) / 10.0); + } + + // create a uniform grid centered at 0,0: + long i = 0; + for (long ix = 0; ix < NUM_X_POINTS; ix++) { + for (long iy = 0; iy < NUM_Y_POINTS; iy++) { + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = ix - 0.5*NUM_X_POINTS; + x->data.F32[1] = iy - 0.5*NUM_Y_POINTS; + ordinates->data[i] = x; + + // generate the data value + float f = fitFunc(NULL, trueParams, x); + float df = sqrt(f); + + // Add some noise + coordinates->data.F32[i] = f + df*psRandomGaussian(rng); + + errors->data.F32[i] = 1.0 / PS_SQR(df); + + printf("%f %f %f %f\n", x->data.F32[0], x->data.F32[1], coordinates->data.F32[i], errors->data.F32[i]); + + if (VERBOSE) { + printf("Data %ld: (%f, %f) --> %f\n", + i, x->data.F32[0], x->data.F32[1], coordinates->data.F32[i]); + } + + i++; + } + } + +# if (1) + psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# else + psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 5); + tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# endif + ok(tmpBool, "psMinimizeLMChi2() suceeded"); + skip_start(!tmpBool, 4, "Skipping tests because psMinimizeLMChi2() failed"); + + printf("Minimisation took %d iterations\n", min->iter); + printf("chi^2 at the minimum is %.3g\n", min->value); + printf("reduced chi^2 at the minimum is %.3g\n", min->value / (float) (NUM_DATA_POINTS - 3)); + for (long i = 0; i < NUM_PARAMS; i++) + { + printf("Parameter %ld at the minimum is %.3f, expected: %f\n", i, + params->data.F32[i], trueParams->data.F32[i]); + } + float diff = 0.0; + for (long i = 0; i < NUM_DATA_POINTS; i++) + { + psVector *x = ordinates->data[i]; + float fitted = fitFunc(NULL, trueParams, x); + float expected = fitFunc(NULL, params, x); + diff += (fitted - expected) / fabsf(expected); + if (VERBOSE) { + printf("Data point %ld: Fitted: %f, expected: %f\n", i, fitted, expected); + } + } + printf("Mean relative difference is %f\n", diff/(float)NUM_DATA_POINTS); + skip_end(); + psFree(min); + psFree(params); + psFree(trueParams); + psFree(ordinates); + psFree(coordinates); + psFree(errors); + psFree(rng); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } +} Index: trunk/psLib/test/math/tap_psMinimizeLMM_Trail.c =================================================================== --- trunk/psLib/test/math/tap_psMinimizeLMM_Trail.c (revision 35767) +++ trunk/psLib/test/math/tap_psMinimizeLMM_Trail.c (revision 35767) @@ -0,0 +1,542 @@ +#include +#include +#include +#include "tap.h" +#include "pstap.h" + +#define NUM_ITER 10 +#define MIN_TOL 0.0001 +#define MAX_TOL 1.0 + +#define NUM_ITERATIONS 100 +#define NUM_PARAMS 5 +#define VERBOSE 0 +float expectedParm[NUM_PARAMS]; + +#define PM_PAR_I0 0 ///< Central intensity +#define PM_PAR_XPOS 1 ///< X center of object +#define PM_PAR_YPOS 2 ///< Y center of object +#define PM_PAR_LENGTH 3 ///< trail length +#define PM_PAR_THETA 4 ///< position angle + +static float SIGMA = 0.0; // cross-width is fixed (not fitted) + +psF32 fitFunc(psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 ST = sin(PAR[PM_PAR_THETA]); + psF32 CT = cos(PAR[PM_PAR_THETA]); + + psF32 S2 = 2.0 * PS_SQR(SIGMA); + + psF32 Zp = (X*CT + Y*ST + 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + psF32 Zm = (X*CT + Y*ST - 0.5*PAR[PM_PAR_LENGTH]) / sqrt(S2); + + psF32 Ep = erf(Zp); + psF32 Em = erf(Zm); + + psF32 Rxy = Y*CT - X*ST; + psF32 Gxy = exp(-Rxy*Rxy/S2); + + psF32 Pxy = Gxy * (Ep - Em); + psF32 f = Pxy * PAR[PM_PAR_I0]; + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + dPAR[PM_PAR_I0] = Pxy; + + float dGdR = -2.0 * Rxy * Gxy / S2; // -R Gxy / (2 Sigma^2) + + // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + float dRdX = +ST; + float dRdY = -CT; + float dRdT = -Y*ST - X*CT; + + float dGdX = dGdR * dRdX; + float dGdY = dGdR * dRdY; + float dGdT = dGdR * dRdT; + + // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo + float dZpdX = -CT / sqrt(S2); + float dZmdX = dZpdX; // float dZmdX = -CT / sqrt(S2); dZmdX = dZpdX + + float dZpdY = -ST / sqrt(S2); // float dZmdY = -ST / sqrt(S2); dZmdY = dZpdY + float dZmdY = dZpdY; + + float dZpdL = +0.5 / sqrt(S2); + float dZmdL = -0.5 / sqrt(S2); + + float dZpdT = (-X*ST + Y*CT) / sqrt(S2); + float dZmdT = dZpdT; // dZpdT = dZmdT + + float dEdZp = exp (-Zp*Zp) * M_2_SQRTPI; + float dEdZm = exp (-Zm*Zm) * M_2_SQRTPI; + + float dEpdX = dEdZp * dZpdX; + float dEmdX = dEdZm * dZmdX; + + float dEpdY = dEdZp * dZpdY; + float dEmdY = dEdZm * dZmdY; + + float dEpdL = dEdZp * dZpdL; + float dEmdL = dEdZm * dZmdL; + + float dEpdT = dEdZp * dZpdT; + float dEmdT = dEdZm * dZmdT; + + float dPdX = dGdX * (Ep - Em) + Gxy * (dEpdX - dEmdX); + float dPdY = dGdY * (Ep - Em) + Gxy * (dEpdY - dEmdY); + + // dGdL is 0.0 because dRdL is 0.0 + float dPdL = Gxy * (dEpdL - dEmdL); + + float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT); + + dPAR[PM_PAR_XPOS] = PAR[PM_PAR_I0] * dPdX; + dPAR[PM_PAR_YPOS] = PAR[PM_PAR_I0] * dPdY; + + dPAR[PM_PAR_LENGTH] = PAR[PM_PAR_I0] * dPdL; + dPAR[PM_PAR_THETA] = PAR[PM_PAR_I0] * dPdT; + } + return(f); +} + +psS32 main(psS32 argc, char* argv[]) +{ + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + plan_tests(35); + + // Test psMinimizeLMChi2() with short-ish trail + if (0) { +# define NUM_X_POINTS 20 +# define NUM_Y_POINTS 20 + int NUM_DATA_POINTS = NUM_X_POINTS * NUM_Y_POINTS; + psMemId id = psMemGetId(); + psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Random number generator; using known seed + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, MIN_TOL, MAX_TOL); + psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS); + psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + + trueParams->data.F32[PM_PAR_I0] = 100.0; + trueParams->data.F32[PM_PAR_XPOS] = 0.0; + trueParams->data.F32[PM_PAR_YPOS] = 0.0; + trueParams->data.F32[PM_PAR_LENGTH] = 10.0; + trueParams->data.F32[PM_PAR_THETA] = 0.0; + + SIGMA = 1.5; + + // Set parameters (Ensure we're not starting right on the true value) + + // scatter flux & length by 10%) + params->data.F32[PM_PAR_I0] = trueParams->data.F32[PM_PAR_I0] * (1.0 - psRandomGaussian(rng) / 10.0); + params->data.F32[PM_PAR_LENGTH] = trueParams->data.F32[PM_PAR_LENGTH] * (1.0 - psRandomGaussian(rng) / 10.0); + + // scatter position by +/- 1 pix + params->data.F32[PM_PAR_XPOS] = trueParams->data.F32[PM_PAR_XPOS] + psRandomGaussian(rng); + params->data.F32[PM_PAR_YPOS] = trueParams->data.F32[PM_PAR_YPOS] + psRandomGaussian(rng); + + // scatter angle by +/- 10 degrees + params->data.F32[PM_PAR_THETA] = trueParams->data.F32[PM_PAR_THETA] + psRandomGaussian(rng) * (10.0 * PS_RAD_DEG); + + // create a uniform grid centered at 0,0: + long i = 0; + for (long ix = 0; ix < NUM_X_POINTS; ix++) { + for (long iy = 0; iy < NUM_Y_POINTS; iy++) { + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = ix - 0.5*NUM_X_POINTS; + x->data.F32[1] = iy - 0.5*NUM_Y_POINTS; + ordinates->data[i] = x; + + // generate the data value (high sky @ 100) + float f = fitFunc(NULL, trueParams, x); + float df = sqrt(f) + 5; + + // Add some noise + coordinates->data.F32[i] = f + df*psRandomGaussian(rng); + + errors->data.F32[i] = 1.0 / PS_SQR(df); + + // printf("%f %f %f %f\n", x->data.F32[0], x->data.F32[1], coordinates->data.F32[i], errors->data.F32[i]); + + if (VERBOSE) { + printf("Data %ld: (%f, %f) --> %f\n", + i, x->data.F32[0], x->data.F32[1], coordinates->data.F32[i]); + } + + i++; + } + } + +# if (1) + psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + bool tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# else + psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 5); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# endif + ok(tmpBool, "psMinimizeLMChi2() suceeded"); + skip_start(!tmpBool, 4, "Skipping tests because psMinimizeLMChi2() failed"); + + printf("Minimisation took %d iterations\n", min->iter); + printf("chi^2 at the minimum is %.3g\n", min->value); + printf("reduced chi^2 at the minimum is %.3g\n", min->value / (float) (NUM_DATA_POINTS - 3)); + for (long i = 0; i < NUM_PARAMS; i++) + { + printf("Parameter %ld at the minimum is %.3f, expected: %f\n", i, + params->data.F32[i], trueParams->data.F32[i]); + } + float diff = 0.0; + for (long i = 0; i < NUM_DATA_POINTS; i++) + { + psVector *x = ordinates->data[i]; + float fitted = fitFunc(NULL, trueParams, x); + float expected = fitFunc(NULL, params, x); + diff += (fitted - expected) / fabsf(expected); + if (VERBOSE) { + printf("Data point %ld: Fitted: %f, expected: %f\n", i, fitted, expected); + } + } + printf("Mean relative difference is %f\n", diff/(float)NUM_DATA_POINTS); + skip_end(); + psFree(min); + psFree(params); + psFree(trueParams); + psFree(ordinates); + psFree(coordinates); + psFree(errors); + psFree(rng); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + # undef NUM_X_POINTS + # undef NUM_Y_POINTS + } + + // Test psMinimizeLMChi2() with long trail & marginal guess + if (0) { +# define NUM_X_POINTS 100 +# define NUM_Y_POINTS 10 + int NUM_DATA_POINTS = NUM_X_POINTS * NUM_Y_POINTS; + psMemId id = psMemGetId(); + // psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Random number generator; using known seed + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, MIN_TOL, MAX_TOL); + psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS); + psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + + trueParams->data.F32[PM_PAR_I0] = 100.0; + trueParams->data.F32[PM_PAR_XPOS] = 0.0; + trueParams->data.F32[PM_PAR_YPOS] = 0.0; + trueParams->data.F32[PM_PAR_LENGTH] = 75.0; + trueParams->data.F32[PM_PAR_THETA] = 0.0; + + SIGMA = 1.5; + + // Set parameters (Ensure we're not starting right on the true value) + + // scatter flux & length by 10%) + params->data.F32[PM_PAR_I0] = trueParams->data.F32[PM_PAR_I0] * (1.0 - psRandomGaussian(rng) / 10.0); + params->data.F32[PM_PAR_LENGTH] = trueParams->data.F32[PM_PAR_LENGTH] * (1.0 - psRandomGaussian(rng) / 3.0); + + // scatter position by +/- 1 pix + params->data.F32[PM_PAR_XPOS] = trueParams->data.F32[PM_PAR_XPOS] + psRandomGaussian(rng); + params->data.F32[PM_PAR_YPOS] = trueParams->data.F32[PM_PAR_YPOS] + psRandomGaussian(rng); + + // scatter angle by +/- 10 degrees + params->data.F32[PM_PAR_THETA] = trueParams->data.F32[PM_PAR_THETA] + psRandomGaussian(rng) * (10.0 * PS_RAD_DEG); + + // create a uniform grid centered at 0,0: + long i = 0; + for (long ix = 0; ix < NUM_X_POINTS; ix++) { + for (long iy = 0; iy < NUM_Y_POINTS; iy++) { + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = ix - 0.5*NUM_X_POINTS; + x->data.F32[1] = iy - 0.5*NUM_Y_POINTS; + ordinates->data[i] = x; + + // generate the data value (high sky @ 100) + float f = fitFunc(NULL, trueParams, x); + float df = sqrt(f) + 5; + + // Add some noise + coordinates->data.F32[i] = f + df*psRandomGaussian(rng); + + errors->data.F32[i] = 1.0 / PS_SQR(df); + + // printf("%f %f %f %f\n", x->data.F32[0], x->data.F32[1], coordinates->data.F32[i], errors->data.F32[i]); + + if (VERBOSE) { + printf("Data %ld: (%f, %f) --> %f\n", + i, x->data.F32[0], x->data.F32[1], coordinates->data.F32[i]); + } + + i++; + } + } + +# if (1) + psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + bool tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# else + psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 5); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# endif + ok(tmpBool, "psMinimizeLMChi2() suceeded"); + skip_start(!tmpBool, 4, "Skipping tests because psMinimizeLMChi2() failed"); + + printf("Minimisation took %d iterations\n", min->iter); + printf("chi^2 at the minimum is %.3g\n", min->value); + printf("reduced chi^2 at the minimum is %.3g\n", min->value / (float) (NUM_DATA_POINTS - 3)); + for (long i = 0; i < NUM_PARAMS; i++) + { + printf("Parameter %ld at the minimum is %.3f, expected: %f\n", i, + params->data.F32[i], trueParams->data.F32[i]); + } + float diff = 0.0; + for (long i = 0; i < NUM_DATA_POINTS; i++) + { + psVector *x = ordinates->data[i]; + float fitted = fitFunc(NULL, trueParams, x); + float expected = fitFunc(NULL, params, x); + diff += (fitted - expected) / fabsf(expected); + if (VERBOSE) { + printf("Data point %ld: Fitted: %f, expected: %f\n", i, fitted, expected); + } + } + printf("Mean relative difference is %f\n", diff/(float)NUM_DATA_POINTS); + skip_end(); + psFree(min); + psFree(params); + psFree(trueParams); + psFree(ordinates); + psFree(coordinates); + psFree(errors); + psFree(rng); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + + // Test psMinimizeLMChi2() with long trail & bad length guess + if (0) { +# define NUM_X_POINTS 100 +# define NUM_Y_POINTS 10 + int NUM_DATA_POINTS = NUM_X_POINTS * NUM_Y_POINTS; + psMemId id = psMemGetId(); + // psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Random number generator; using known seed + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, MIN_TOL, MAX_TOL); + psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS); + psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + + trueParams->data.F32[PM_PAR_I0] = 100.0; + trueParams->data.F32[PM_PAR_XPOS] = 0.0; + trueParams->data.F32[PM_PAR_YPOS] = 0.0; + trueParams->data.F32[PM_PAR_LENGTH] = 75.0; + trueParams->data.F32[PM_PAR_THETA] = 0.0; + + SIGMA = 1.5; + + // Set parameters (Ensure we're not starting right on the true value) + + // scatter flux & length by 10%) + params->data.F32[PM_PAR_I0] = trueParams->data.F32[PM_PAR_I0] * (1.0 - psRandomGaussian(rng) / 10.0); + params->data.F32[PM_PAR_LENGTH] = 10.0; // way too short + + // scatter position by +/- 1 pix + params->data.F32[PM_PAR_XPOS] = trueParams->data.F32[PM_PAR_XPOS] + psRandomGaussian(rng); + params->data.F32[PM_PAR_YPOS] = trueParams->data.F32[PM_PAR_YPOS] + psRandomGaussian(rng); + + // scatter angle by +/- 10 degrees + params->data.F32[PM_PAR_THETA] = trueParams->data.F32[PM_PAR_THETA] + 0.01*psRandomGaussian(rng); + + // create a uniform grid centered at 0,0: + long i = 0; + for (long ix = 0; ix < NUM_X_POINTS; ix++) { + for (long iy = 0; iy < NUM_Y_POINTS; iy++) { + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = ix - 0.5*NUM_X_POINTS; + x->data.F32[1] = iy - 0.5*NUM_Y_POINTS; + ordinates->data[i] = x; + + // generate the data value (high sky @ 100) + float f = fitFunc(NULL, trueParams, x); + float df = sqrt(f) + 5; + + // Add some noise + coordinates->data.F32[i] = f + df*psRandomGaussian(rng); + + errors->data.F32[i] = 1.0 / PS_SQR(df); + + // printf("%f %f %f %f\n", x->data.F32[0], x->data.F32[1], coordinates->data.F32[i], errors->data.F32[i]); + + if (VERBOSE) { + printf("Data %ld: (%f, %f) --> %f\n", + i, x->data.F32[0], x->data.F32[1], coordinates->data.F32[i]); + } + + i++; + } + } + +# if (1) + psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + bool tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# else + psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 5); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# endif + ok(tmpBool, "psMinimizeLMChi2() suceeded"); + skip_start(!tmpBool, 4, "Skipping tests because psMinimizeLMChi2() failed"); + + printf("Minimisation took %d iterations\n", min->iter); + printf("chi^2 at the minimum is %.3g\n", min->value); + printf("reduced chi^2 at the minimum is %.3g\n", min->value / (float) (NUM_DATA_POINTS - 3)); + for (long i = 0; i < NUM_PARAMS; i++) + { + printf("Parameter %ld at the minimum is %.3f, expected: %f\n", i, + params->data.F32[i], trueParams->data.F32[i]); + } + float diff = 0.0; + for (long i = 0; i < NUM_DATA_POINTS; i++) + { + psVector *x = ordinates->data[i]; + float fitted = fitFunc(NULL, trueParams, x); + float expected = fitFunc(NULL, params, x); + diff += (fitted - expected) / fabsf(expected); + if (VERBOSE) { + printf("Data point %ld: Fitted: %f, expected: %f\n", i, fitted, expected); + } + } + printf("Mean relative difference is %f\n", diff/(float)NUM_DATA_POINTS); + skip_end(); + psFree(min); + psFree(params); + psFree(trueParams); + psFree(ordinates); + psFree(coordinates); + psFree(errors); + psFree(rng); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // Test psMinimizeLMChi2() with long trail & bad length guess + if (1) { +# define NUM_X_POINTS 100 +# define NUM_Y_POINTS 10 + int NUM_DATA_POINTS = NUM_X_POINTS * NUM_Y_POINTS; + psMemId id = psMemGetId(); + // psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 1); // Random number generator; using known seed + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psMinimization *min = psMinimizationAlloc(NUM_ITERATIONS, MIN_TOL, MAX_TOL); + psArray *ordinates = psArrayAlloc(NUM_DATA_POINTS); + psVector *coordinates = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *errors = psVectorAlloc(NUM_DATA_POINTS, PS_TYPE_F32); + psVector *params = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + psVector *trueParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + + trueParams->data.F32[PM_PAR_I0] = 100.0; + trueParams->data.F32[PM_PAR_XPOS] = 0.0; + trueParams->data.F32[PM_PAR_YPOS] = 0.0; + trueParams->data.F32[PM_PAR_LENGTH] = 75.0; + trueParams->data.F32[PM_PAR_THETA] = 0.0; + + SIGMA = 1.5; + + // Set parameters (Ensure we're not starting right on the true value) + + // scatter flux & length by 10%) + params->data.F32[PM_PAR_I0] = trueParams->data.F32[PM_PAR_I0] * (1.0 - psRandomGaussian(rng) / 10.0); + params->data.F32[PM_PAR_LENGTH] = 10.0; // way too short + + // scatter position by +/- 1 pix + // params->data.F32[PM_PAR_XPOS] = trueParams->data.F32[PM_PAR_XPOS] + psRandomGaussian(rng); + params->data.F32[PM_PAR_XPOS] = 15.0; + params->data.F32[PM_PAR_YPOS] = trueParams->data.F32[PM_PAR_YPOS] + psRandomGaussian(rng); + + // scatter angle by +/- 10 degrees + params->data.F32[PM_PAR_THETA] = trueParams->data.F32[PM_PAR_THETA] + 0.01*psRandomGaussian(rng); + + // create a uniform grid centered at 0,0: + long i = 0; + for (long ix = 0; ix < NUM_X_POINTS; ix++) { + for (long iy = 0; iy < NUM_Y_POINTS; iy++) { + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = ix - 0.5*NUM_X_POINTS; + x->data.F32[1] = iy - 0.5*NUM_Y_POINTS; + ordinates->data[i] = x; + + // generate the data value (high sky @ 100) + float f = fitFunc(NULL, trueParams, x); + float df = sqrt(f) + 5; + + // Add some noise + coordinates->data.F32[i] = f + df*psRandomGaussian(rng); + + errors->data.F32[i] = 1.0 / PS_SQR(df); + + // printf("%f %f %f %f\n", x->data.F32[0], x->data.F32[1], coordinates->data.F32[i], errors->data.F32[i]); + + if (VERBOSE) { + printf("Data %ld: (%f, %f) --> %f\n", + i, x->data.F32[0], x->data.F32[1], coordinates->data.F32[i]); + } + + i++; + } + } + +# if (1) + // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + bool tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# else + psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 5); + bool tmpBool = psMinimizeLMChi2(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); +# endif + ok(tmpBool, "psMinimizeLMChi2() suceeded"); + skip_start(!tmpBool, 4, "Skipping tests because psMinimizeLMChi2() failed"); + + printf("Minimisation took %d iterations\n", min->iter); + printf("chi^2 at the minimum is %.3g\n", min->value); + printf("reduced chi^2 at the minimum is %.3g\n", min->value / (float) (NUM_DATA_POINTS - 3)); + for (long i = 0; i < NUM_PARAMS; i++) + { + printf("Parameter %ld at the minimum is %.3f, expected: %f\n", i, + params->data.F32[i], trueParams->data.F32[i]); + } + float diff = 0.0; + for (long i = 0; i < NUM_DATA_POINTS; i++) + { + psVector *x = ordinates->data[i]; + float fitted = fitFunc(NULL, trueParams, x); + float expected = fitFunc(NULL, params, x); + diff += (fitted - expected) / fabsf(expected); + if (VERBOSE) { + printf("Data point %ld: Fitted: %f, expected: %f\n", i, fitted, expected); + } + } + printf("Mean relative difference is %f\n", diff/(float)NUM_DATA_POINTS); + skip_end(); + psFree(min); + psFree(params); + psFree(trueParams); + psFree(ordinates); + psFree(coordinates); + psFree(errors); + psFree(rng); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } +} Index: trunk/psLib/test/pstap/src/pstap.h =================================================================== --- trunk/psLib/test/pstap/src/pstap.h (revision 35674) +++ trunk/psLib/test/pstap/src/pstap.h (revision 35767) @@ -144,2 +144,12 @@ } \ } + +# define ok_float is_float +# define ok_float_tol is_float_tol +# define ok_double is_double +# define ok_double_tol is_double_tol +# define ok_str is_str +# define ok_strn is_strn +# define ok_int is_int +# define ok_long is_long +# define ok_bool is_bool