Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1831) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 1836) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-18 01:50:45 $ + * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-20 23:16:10 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -94,29 +94,4 @@ /*****************************************************************************/ -typedef struct -{ - size_t n; // Number of data points points in domain. - int paramCount; // Number of non-masked parameters. - psVector* restrict initialGuess; - const psImage* restrict domain; - const psVector* restrict data; - const psVector* restrict errors; - const psVector* restrict paramMask; - psMinimizeFunction evalModel; - psMinimizeFunctionDeriv d_evalModel; -} -psMinChi2Data; - -typedef struct -{ - int paramCount; // Number of non-masked parameters. - psVector* restrict initialGuess; - const psVector* restrict coord; - const psVector* restrict paramMask; - psMinimizeFunction evalModel; - psMinimizeFunctionDeriv d_evalModel; -} -psMinimizeData; - /*****************************************************************************/ /* GLOBAL VARIABLES */ @@ -137,281 +112,4 @@ /* FUNCTION IMPLEMENTATION - LOCAL */ /*****************************************************************************/ - -/****************************************************************************** -p_psMinFunc(*params, *funcData): We use the GSL procedure -gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied -by the user. That GSL procedure requires the function to be minimized to be -in a different format than the psLib format. The purpose of this procedure -is to serve as a GSL-format wrapper for the user-supplied procedure which is -to be minimized. - - params: The parameters of the function to be minimized. These will be - varied by GSL in order to minimize the function. - - funcData: a private psLib struct which contains the data point to be - minimized, the function and derivative function pointers, an initial - guess at the parameters, an option parameter mask, etc. - *****************************************************************************/ -double p_psMinFunc(const gsl_vector * params, void *funcData) -{ - int i; // Loop index variable. - int j; // Loop index variable. - float tmpf; // Temporary floating point variable. - const psVector* restrict coord = ((psMinimizeData* ) funcData)->coord; - const psVector* restrict mask = ((psMinimizeData* ) funcData)->paramMask; - psVector* restrict initialGuess = ((psMinimizeData* ) funcData)->initialGuess; - psMinimizeFunction evalModel = ((psMinimizeData* ) funcData)->evalModel; - psVector* inputParameterList = NULL; - - // The GSL routines will call this function with the masked parameters - // removed. However, the user-supplied function (to be modified) does not - // have those parameters removed. Here will create a new parameter list - // with the masked parameters added (we expand initialGuess). - inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32); - if (mask != NULL) { - j = 0; - for (i = 0; i < mask->n; i++) { - if (mask->data.U8[i] != 0) { - inputParameterList->data.F32[i] = initialGuess->data.F32[i]; - } else { - inputParameterList->data.F32[i] = gsl_vector_get(params, j++); - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - inputParameterList->data.F32[i] = gsl_vector_get(params, i); - } - } - - // Call the user-supplied function. - tmpf = evalModel(inputParameterList, coord); - - // Free allocated memory and return the value of the function. - psFree(inputParameterList); - return (tmpf); -} - -/****************************************************************************** -p_psMinFuncDeriv(*params, *funcData): a GSL-like wrapper for the -user-supplied procedure which calculates the derviative of the function to be -minimized. - - params: The parameters of the function to be minimized. These will be - varied by GSL in order to minimize the function. - - funcData: a private psLib struct which contains the data point to be - minimized, the function and derivative function pointers, an initial - guess at the parameters, an option parameter mask, etc. - - df: we calculate the derivative of the function w.r.t. to each parameter - in "params" and return those derivatives in this psVector. - *****************************************************************************/ -void p_psMinFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * df) -{ - int i; // Loop index variable. - int j; // Loop index variable. - float tmpf; // Temporary floating point variable. - const psVector* restrict coord = ((psMinimizeData* ) funcData)->coord; - const psVector* restrict mask = ((psMinimizeData* ) funcData)->paramMask; - psVector* restrict initialGuess = ((psMinimizeData* ) funcData)->initialGuess; - psMinimizeFunctionDeriv d_evalModel = ((psMinimizeData* ) funcData)->d_evalModel; - psVector* inputParameterList = NULL; - - // The GSL routines will call this function with the masked parameters - // removed. However, the user-supplied function (to be modified) does not - // have those parameters removed. Here will create a new parameter list - // with the masked parameters added (we expand initialGuess). - inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32); - if (mask != NULL) { - j = 0; - for (i = 0; i < mask->n; i++) { - if (mask->data.U8[i] != 0) { - inputParameterList->data.F32[i] = initialGuess->data.F32[i]; - } else { - inputParameterList->data.F32[i] = gsl_vector_get(params, j++); - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - inputParameterList->data.F32[i] = gsl_vector_get(params, i); - } - } - - // Evaluate the derivative w.r.t. each parameter. - // NOTE: we can probably remove the calls for masked parameters. - for (i = 0; i < initialGuess->n; i++) { - tmpf = d_evalModel(inputParameterList, coord, i); - gsl_vector_set(df, i, tmpf); - } - - // Free allocated memory. - psFree(inputParameterList); -} - -/****************************************************************************** - Compute both p_psMinFunc and p_psMinFuncDeriv together. - *****************************************************************************/ -void p_psMinFuncFuncDeriv(const gsl_vector * params, void *funcData, double *f, gsl_vector * df) -{ - *f = p_psMinFunc(params, funcData); - p_psMinFuncDeriv(params, funcData, df); -} - -// The first argument to evalModel() and d_evalModel() specifies the data -// point. It must have the same size as the second dimension of *domain. -// The second argument must have the same size as *initialGuess and -// *paramMask. - -/****************************************************************************** -p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL procedure -gsl_multifit_fdfsolver_iterate() to fit an arbitrary function, supplied by -the user, to a set of data points. That GSL procedure requires the function -to be fit to be in a different format than the psLib format. The purpose of -this procedure is to serve as a GSL-format wrapper for the user-supplied -procedure which is to be fit to the data. - - params: These are the parameters which are to be varied by GSL in order - to minimize chi2 over the data set. - - funcData: this data structure contains the input values over which the - function will be evaluated, the expected value of the function at - those points, the amount of error tolerable at those points, a mask - vector which specifies which parameters to the function are to be - constant, and an initial guess at the parameters. - - outData: The function is evaluated at each point, then subtract the - expected value and divide by the error. - *****************************************************************************/ -int p_psMinChi2Func(const gsl_vector * params, void *funcData, gsl_vector * outData) -{ - int i; // Loop index variable. - int j; // Loop index variable. - float tmpf; // Temporary floating point variable. - const psImage* restrict domain = ((psMinChi2Data* ) funcData)->domain; - const psVector* restrict data = ((psMinChi2Data* ) funcData)->data; - const psVector* restrict errors = ((psMinChi2Data* ) funcData)->errors; - const psVector* restrict mask = ((psMinChi2Data* ) funcData)->paramMask; - psVector* restrict initialGuess = ((psMinChi2Data* ) funcData)->initialGuess; - psMinimizeFunction evalModel = ((psMinChi2Data* ) funcData)->evalModel; - psVector* inputParameterList = NULL; - psVector* tmpVecPtr = NULL; - - tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); - - // The GSL routines will call this function with the masked parameters - // removed. However, the user-supplied function (to be modified) does not - // have those parameters removed. Here will create a new parameter list - // with the masked parameters added (we expand initialGuess). - - inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32); - if (mask != NULL) { - j = 0; - for (i = 0; i < mask->n; i++) { - if (mask->data.U8[i] != 0) { - inputParameterList->data.F32[i] = initialGuess->data.F32[i]; - } else { - inputParameterList->data.F32[i] = gsl_vector_get(params, j++); - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - inputParameterList->data.F32[i] = gsl_vector_get(params, i); - } - } - - // Evaluate the function at each data point. - for (i = 0; i < domain->numRows; i++) { - for (j = 0; j < domain->numCols; j++) { - tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; - } - tmpf = evalModel(tmpVecPtr, inputParameterList); - - gsl_vector_set(outData, i, (tmpf - data->data.F32[i]) / errors->data.F32[i]); - } - - // Free allocated memory. - psFree(inputParameterList); - psFree(tmpVecPtr); - - return GSL_SUCCESS; -} - -/****************************************************************************** -p_psMinChi2FuncDeriv(*x, *funcData, *outdata): a GSL-like wrapper for the -user-supplied procedure which calculates the derviative of the function to be -minimized. - params: These are the parameters which are to be varied by GSL in order - to minimize chi2 over the data set. - - funcData: this data structure contains the input values over which the - function will be evaluated, the expected value of the function at - those points, the amount of error tolerable at those points, a mask - vector which specifies which parameters to the function are to be - constant, and an initial guess at the parameters. - - J: The derivative is evaluated at each point and w.r.t. each parameter - and returned in this data structure. - *****************************************************************************/ -int p_psMinChi2FuncDeriv(const gsl_vector * params, void *funcData, gsl_matrix * J) -{ - const psImage* restrict domain = ((psMinChi2Data* ) funcData)->domain; - const psVector* restrict errors = ((psMinChi2Data* ) funcData)->errors; - const psVector* restrict mask = ((psMinChi2Data* ) funcData)->paramMask; - psVector* restrict initialGuess = ((psMinChi2Data* ) funcData)->initialGuess; - psVector* inputParameterList = NULL; - psVector* tmpVecPtr = NULL; - psMinimizeFunctionDeriv d_evalModel = ((psMinChi2Data* ) funcData)->d_evalModel; - - size_t i; - int j; - float tmpf; - - tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); - - // The GSL routines will call this functions with the masked parameters - // removed. However, the user-supplied function (to be modified) does not - // have those parameters removed. Here will create a new parameter list - // with the masked parameters added (we expand initialGuess). - - inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32); - if (mask != NULL) { - j = 0; - for (i = 0; i < mask->n; i++) { - if (mask->data.U8[i] != 0) { - inputParameterList->data.F32[i] = initialGuess->data.F32[i]; - } else { - inputParameterList->data.F32[i] = gsl_vector_get(params, j++); - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - inputParameterList->data.F32[i] = gsl_vector_get(params, i); - } - } - - // Evaluate the derivtaive at each data point, and w.r.t. each parameter. - for (i = 0; i < domain->numRows; i++) { - for (j = 0; j < tmpVecPtr->n; j++) { - tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; - } - - for (j = 0; j < inputParameterList->n; j++) { - tmpf = d_evalModel(tmpVecPtr, inputParameterList, j); - gsl_matrix_set(J, i, j, (tmpf / errors->data.F32[i])); - } - } - - psFree(inputParameterList); - psFree(tmpVecPtr); - return GSL_SUCCESS; -} - -int p_psMinChi2FuncFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * f, gsl_matrix * J) -{ - p_psMinChi2Func(params, funcData, f); - p_psMinChi2FuncDeriv(params, funcData, J); - - return GSL_SUCCESS; -} /****************************************************************************** @@ -433,604 +131,30 @@ /****************************************************************************** -p_psBuildSums1D(x): this routine returns a psVector with "x" elements. The +VectorNormalizeGen(): this routine returns a psVector with "x" elements. The values of the vector will be scaled uniformly between -1.0 and 1.0. - *****************************************************************************/ -psVector* psBuildImageScalingFactors(int x) +*****************************************************************************/ +psVector* VectorNormalizeGen(int x) { - int i = 0; // loop index variable. - psVector* imageScalingFactors = NULL; - - imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32); + int i = 0; + psVector *tmp = NULL; + + tmp = psVectorAlloc(x, PS_TYPE_F32); for (i = 0; i < x; i++) { - imageScalingFactors->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0; - } - - return (imageScalingFactors); -} + tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0; + } + + return (tmp); +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ /****************************************************************************** -CURRENTLY NOT IN USE. - -p_psPolyOrderCheck(A, N, *indx, *B, polyOrder,*flag) This routine checks if -all polyOrder-th terms in the polyOrder-th order sky background polynomial -defined by the coefficients in the array B[] are consistent with zero. If -true, then *flag is set to 1. Otherwise, *flag is set to 0. The matrix -inversion code in the middle of this procedure draws from Numerical Recipes -in C page 48. -Input: - A This is the LUD decomposition of the original matrix A. - N The size of the matrix (plus 1, actually, since offset 1). - indx misc Numerical Recipes data structure. - B The coefficients of the sky polynomial. - polyOrder The degree of the sky polynomial. -Output: - *flag Set this to 1 if we must recalculate the coefficients. - *****************************************************************************/ -void p_psPolyOrderCheck(float **A, int N, int *indx, float *B, int polyOrder, int *flag) -{ - float **y = NULL; // This 2-D matrix will hold A^-1 - float *col = NULL; // misc NumerRecipes data structure - float *error = NULL; // will hold the sqrt() of the - - // diagonal of y[][]. - int i = 0; // loop-index variable - int j = 0; // loop-index variable - int numPolyTerms = 0; // The number of terms in the - - // polynomial. - int lastTerm = 0; // The index location of the first - - // n-th order term in array B[]. - int firstTerm = 0; // Index location of last such term. - - // Allocate the necessary data structures for this procedure... - error = (float *)psAlloc((N + 1) * sizeof(float)); - col = (float *)psAlloc((N + 1) * sizeof(float)); - y = (float **)psAlloc((N + 1) * sizeof(float *)); - for (i = 1; i <= N; i++) { - y[i] = (float *)psAlloc((N + 1) * sizeof(float)); - } - - // Invert the matrix A and put the result in y[][]. This code is taken - // from Numerical Recipes in C page 48. - for (j = 1; j <= N; j++) { - for (i = 1; i <= N; i++) { - col[i] = 0.0; - } - col[j] = 1.0; - // NOTE: substitue the LUD rotine - // lubksb(A, N, indx, col); - for (i = 1; i <= N; i++) { - y[i][j] = col[i]; - } - } - - // Determine where the first n-th order (in this comment, n equals - // polyOrder) polynomial term is stored in the matrix B[], and also were - // the last n-order term is stored. Then we loop over all the n-order - // terms and check if they are consistent with zero. - - numPolyTerms = (((polyOrder + 1) * (polyOrder + 2)) / 2); - lastTerm = numPolyTerms + 1; - firstTerm = lastTerm - polyOrder; - *flag = 1; - for (i = firstTerm; i <= lastTerm; i++) { - #ifdef DARWIN - error[i] = (float)sqrt(y[i][i]); - #else - - error[i] = sqrtf(y[i][i]); - #endif - - if (!((B[i] <= (2.0f * error[i])) && ((-2.0f * error[i]) <= B[i]))) { - *flag = 0; - } - } - - // Free all memory allocated in this routine. - psFree(error); - psFree(col); - for (j = 1; j <= N; j++) { - psFree(y[j]); - } - psFree(y); -} - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ - -/****************************************************************************** -psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask): - -This routine must minimize an arbitrary function; it determines the set of -parameters of that function such that the ... - *****************************************************************************/ -psVector* psMinimize(psVector* restrict initialGuess, - psMinimizeFunction myFunction, - psMinimizeFunctionDeriv myFunctionDeriv, - const psVector* restrict coord, - const psVector* restrict paramMask) -{ - int status; - int i = 0; - int j = 0; - int iter = 0; - gsl_multimin_function_fdf f; - const gsl_multimin_fdfminimizer_type *T; - gsl_multimin_fdfminimizer *s; - psMinimizeData inputData; - gsl_vector *x; - - psTrace(".psLib.dataManip.psMinimize", 4, - "---- psMinimize() begin ----\n"); - - PS_CHECK_NULL_VECTOR(initialGuess); - PS_CHECK_EMPTY_VECTOR(initialGuess); - PS_CHECK_NULL_VECTOR(coord); - PS_CHECK_EMPTY_VECTOR(coord); - if (paramMask != NULL) { - PS_CHECK_NULL_VECTOR(paramMask); - PS_CHECK_EMPTY_VECTOR(paramMask); - PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask); - } - - inputData.initialGuess = initialGuess; - inputData.coord = coord; - inputData.paramMask = paramMask; - inputData.evalModel = myFunction; - inputData.d_evalModel = myFunctionDeriv; - inputData.paramCount = 0; - - // If the user supplied a parameter mask, then count the number of - // non-masked elements. This will be used later in allocating a vector - // for the parameters. - if (paramMask != NULL) { - for (i = 0; i < paramMask->n; i++) { - if (paramMask->data.U8[i] != 0) { - inputData.paramCount++; - } - } - } else { - inputData.paramCount = initialGuess->n; - } - - // The initial guess at the parameters for the function are written into - // the vector inputParameterList. If the paramMask is not NULL, then - // masked parameters are masked out. - x = gsl_vector_alloc(inputData.paramCount); - if (paramMask != NULL) { - j = 0; - for (i = 0; i < initialGuess->n; i++) { - if (paramMask->data.U8[i] == 0) { - gsl_vector_set(x, j++, initialGuess->data.F32[i]); - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - gsl_vector_set(x, i, initialGuess->data.F32[i]); - } - } - f.f = &p_psMinFunc; - f.df = &p_psMinFuncDeriv; - f.fdf = &p_psMinFuncFuncDeriv; - f.n = inputData.paramCount; - f.params = &inputData; - - T = gsl_multimin_fdfminimizer_conjugate_fr; - s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount); - gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4); - do { - iter++; - status = gsl_multimin_fdfminimizer_iterate(s); - - if (status) - break; - - status = gsl_multimin_test_gradient(s->gradient, 1e-3); - - if (status == GSL_SUCCESS) - printf("Minimum found at:\n"); - - } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS); - - // In the above steps we had removed the masked elements from the - // the solver. This next code blocks puts those masked elements - // into the solution. - if (paramMask != NULL) { - j = 0; - for (i = 0; i < initialGuess->n; i++) { - if (paramMask->data.U8[i] == 0) { - initialGuess->data.F32[i] = gsl_vector_get(s->x, j++); - } else { - initialGuess->data.F32[i] = initialGuess->data.F32[i]; - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - initialGuess->data.F32[i] = gsl_vector_get(s->x, i); - } - } - psTrace(".psLib.dataManip.psMinimize", 4, - "---- psMinimize() end ----\n"); - return (initialGuess); -} - -/****************************************************************************** - This routine must determine the parameters of an arbitrary function - such that they best fit the supplied data points. - *****************************************************************************/ -psVector* psMinimizeChi2(psMinimizeFunction evalModel, - psMinimizeFunctionDeriv DevalModel, - const psImage* restrict domain, - const psVector* restrict data, - const psVector* restrict errors, - psVector* restrict initialGuess, - const psVector* restrict paramMask, - float *chiSq) -{ - int numData = domain->numRows; // Number of data points - int status; // Return status for the GSL solver. - int i = 0; // Loop index variable. - int j = 0; // Loop index variable. - int iter = 0; // Iteration counter. - gsl_multifit_function_fdf f; // GSL structure that contains the - - // functions/derivative to be solved. - double *xInit = NULL; // The initial guess at the parameters - - // with masked parameters removed. - const gsl_multifit_fdfsolver_type *T; - - // This tells GSL to use the Levenberg- - // Marquardt algorithm for chi2 - // minimization. - gsl_multifit_fdfsolver *s; // GSL data structure. - psMinChi2Data inputData; - float chiSqOld = 0.0; - - psTrace(".psLib.dataManip.psMinimizeChi2", 4, - "---- psMinimizeChi2() begin ----\n"); - - PS_CHECK_NULL_IMAGE(domain); - PS_CHECK_EMPTY_IMAGE(domain); - PS_CHECK_NULL_VECTOR(data); - PS_CHECK_EMPTY_VECTOR(data); - PS_CHECK_NULL_VECTOR(errors); - PS_CHECK_EMPTY_VECTOR(errors); - PS_CHECK_NULL_VECTOR(initialGuess); - PS_CHECK_EMPTY_VECTOR(initialGuess); - PS_CHECK_VECTOR_SIZE_EQUAL(data, errors); - if (domain->numRows != data->n) { - psAbort(__func__, "Number of data points and data values not equal."); - } - if (paramMask != NULL) { - PS_CHECK_NULL_VECTOR(paramMask); - PS_CHECK_EMPTY_VECTOR(paramMask); - PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask); - } - - inputData.n = numData; - inputData.paramCount = 0; - inputData.initialGuess = initialGuess; - inputData.domain = domain; - inputData.data = data; - inputData.errors = errors; - inputData.paramMask = paramMask; - inputData.evalModel = evalModel; - inputData.d_evalModel = DevalModel; - - // If the user supplied a parameter mask, then count the number of - // non-masked elements. This will be used later in allocating a vector - // for the parameters. - if (paramMask != NULL) { - for (i = 0; i < paramMask->n; i++) { - if (paramMask->data.U8[i] != 0) { - inputData.paramCount++; - } - } - } else { - inputData.paramCount = initialGuess->n; - } - - // The initial guess at the parameters for the function are written into - // the vector inputParameterList. If the paramMask is not NULL, then those - // parameters are masked out. - xInit = (double *)psAlloc(inputData.paramCount * sizeof(double)); - if (paramMask != NULL) { - j = 0; - for (i = 0; i < initialGuess->n; i++) { - if (paramMask->data.U8[i] == 0) { - xInit[j++] = initialGuess->data.F32[i]; - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - xInit[i] = initialGuess->data.F32[i]; - } - } - - const gsl_rng_type *type; - gsl_rng *r; - - gsl_rng_env_setup(); - - type = gsl_rng_default; - r = gsl_rng_alloc(type); - - // Initialize the main data structure used by the GSL solver. This will - // contain pointers to the function to be minimized, it's derivative - // function, the number of data points, the number of free parameters, - // and the data structures those functions use. - - f.f = &p_psMinChi2Func; - f.df = &p_psMinChi2FuncDeriv; - f.fdf = &p_psMinChi2FuncFuncDeriv; - f.n = numData; - f.p = inputData.paramCount; - f.params = &inputData; - - gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount); - - T = gsl_multifit_fdfsolver_lmsder; - s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount); - gsl_multifit_fdfsolver_set(s, &f, &x.vector); - *chiSq = 0.0; - chiSqOld = 0.0; - do { - iter++; - for (i = 0; i < initialGuess->n; i++) { - printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i)); - } - // Perform an iteration of the GSL solver. - status = gsl_multifit_fdfsolver_iterate(s); - printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status)); - for (i = 0; i < initialGuess->n; i++) { - printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i)); - } - - // If there was a problem, abort. - if (status) { - psAbort(__func__, "gsl_multifit_fdfsolver_iterate(%s)\n", gsl_strerror(status)); - } - // Test if the parameters changed by a small enough amount. - // NOTE: This wasn't working right when the parameters fit exactly. - // Figure out why. - // status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); - - // We test for convergence if chiSquared changes by less than 1.0 - // as specified in the ADD. - *chiSq = gsl_blas_dnrm2(s->f); - printf("psMinimize.c: chiSq is %.3f\n", *chiSq); - if (fabs(*chiSq - chiSqOld) < 1.0) { - status = GSL_SUCCESS; - } else { - status = GSL_CONTINUE; - } - chiSqOld = *chiSq; - - } while (status == GSL_CONTINUE && iter < MAX_LMM_ITERATIONS); - - // In the above steps we had removed the masked elements from the - // the solver. This next code blocks puts those masked elements - // into the solution. - if (paramMask != NULL) { - j = 0; - for (i = 0; i < initialGuess->n; i++) { - if (paramMask->data.U8[i] == 0) { - initialGuess->data.F32[i] = gsl_vector_get(s->x, j++); - } else { - initialGuess->data.F32[i] = initialGuess->data.F32[i]; - } - } - } else { - for (i = 0; i < initialGuess->n; i++) { - initialGuess->data.F32[i] = gsl_vector_get(s->x, i); - } - } - - // Calculate the chi-squared for the derived solution. - *chiSq = gsl_blas_dnrm2(s->f); - - // Free all allocated memory - // NOTE: Free x. - gsl_multifit_fdfsolver_free(s); - psFree(xInit); - - // Bye bye. - psTrace(".psLib.dataManip.psMinimizeChi2", 4, - "---- psMinimizeChi2() end ----\n"); - return (initialGuess); -} - -/****************************************************************************** -This routine will take an procedure which calculates an arbitrary function -and it's derivative and minimize it. - -XXX: Do this: - After checking that all entries in the paramMask are 1 or 0, when - forming the A matrix from alpha, try this: - - A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; - *****************************************************************************/ -bool psMinimizeLM(psMinimization *min, - psImage *covar, - psVector *params, - const psVector *paramMask, - const psArray *coords, - psMinimizeLMFunc func) -{ - psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64); - psVector *perm = psVectorAlloc(params->n, PS_TYPE_F64); - psVector *newParamsF64 = psVectorAlloc(params->n, PS_TYPE_F64); - psVector *newParamsF32 = psVectorAlloc(params->n, PS_TYPE_F32); - psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32); - psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); - psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F64); - psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F64); - psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32); - psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32); - psVector *NRparams = psVectorAlloc(params->n, PS_TYPE_F32); - int i; - int j; - int k; - float newValue; - float oldValue; - float lamda = 0.00005; - - psTrace(".psLib.dataManip.psMinimizeLM", 4, - "---- psMinimizeLM() begin ----\n"); - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "min->tol is %f\n", min->tol); - - for (i=0;in;i++) { - origParams->data.F32[i] = params->data.F32[i]; - } - - min->lastDelta = HUGE; - min->lastDelta = 12345.0; - min->iter = 0; - - while ((min->lastDelta > min->tol) && - (min->iter < min->maxIter)) { - psTrace(".psLib.dataManip.psMinimizeLM", 4, - "------------------------------------------------------\n"); - psTrace(".psLib.dataManip.psMinimizeLM", 4, - "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); - - min->value = func(deriv, params, coords); - - for (i=0;in;i++) { - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "params->data.F32[%d] is %f.\n", i, params->data.F32[i]); - } - for (i=0;in;i++) { - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "deriv->data.F32[%d] is %f.\n", i, deriv->data.F32[i]); - } - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "min->value is (%f)\n", min->value); - - for (i=0;in;i++) { - if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { - deriv->data.F32[i] = 0.0; - } - } - - // Calculate the BETA vector. - for (i=0;in;i++) { - beta->data.F64[i] = (float) deriv->data.F32[i]; - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "beta->data.F64[%d] is %f.\n", i, beta->data.F64[i]); - } - - // Calculate the ALPHA matrix. - for (j=0;jn;j++) { - for (k=0;kn;k++) { - alpha->data.F32[j][k] = deriv->data.F32[k] * - deriv->data.F32[j]; - } - } - - // Calculate the matrix A. - for (j=0;jn;j++) { - for (k=0;kn;k++) { - if (j == k) { - A->data.F64[j][k] = - (double) ((1.0 + lamda) * alpha->data.F32[j][k]); - } else { - A->data.F64[j][k] = (double) alpha->data.F32[j][k]; - } - } - } - for (j=0;jn;j++) { - psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n"); - for (k=0;kn;k++) { - psTrace(".psLib.dataManip.psMinimizeLM", 6, "%f ", A->data.F64[j][k]); - } - psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n"); - } - - // Solve A * alpha = Beta - aOut = psMatrixLUD(aOut, perm, A); - newParamsF64 = psMatrixLUSolve(newParamsF64, aOut, beta, perm); - - // Mask any masked parameters. - for (i=0;in;i++) { - if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { - newParamsF64->data.F64[i] = (double) origParams->data.F32[i]; - } - newParamsF32->data.F32[i] = (float) newParamsF64->data.F64[i]; - - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "newParamsF32->data.F32[%d] is %f.\n", i, newParamsF32->data.F32[i]); - NRparams->data.F32[i] = params->data.F32[i] - newParamsF32->data.F32[i]; - } - - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "Calling func() with new parameters:\n"); - for (i=0;in;i++) { - psTrace(".psLib.dataManip.psMinimizeLM", 6, - "NRparams->data.F32[%d] is %f.\n", i, NRparams->data.F32[i]); - } - - oldValue = min->value; - newValue = func(deriv, NRparams, coords); - psTrace(".psLib.dataManip.psMinimizeLM", 4, - "old/new values are (%f, %f)\n", oldValue, newValue); - for (i=0;in;i++) { - if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { - deriv->data.F32[i] = 0.0; - } - } - - if (oldValue > newValue) { - min->lastDelta = oldValue - newValue; - min->value = newValue; - - // No need to check the paramMask here since we already did so - // before the last function evaluation. - for (i=0;in;i++) { - // params->data.F32[i] = (float) newParamsF64->data.F64[i]; - params->data.F32[i] = (float) NRparams->data.F32[i]; - } - min->value = func(deriv, params, coords); - - lamda*= 0.1; - } else { - lamda*= 10.0; - } - psTrace(".psLib.dataManip.psMinimizeLM", 4, - "lamda is %f\n", lamda); - min->iter++; - } - psFree(beta); - psFree(perm); - psFree(newParamsF64); - psFree(newParamsF32); - psFree(origParams); - psFree(alpha); - psFree(A); - psFree(aOut); - psFree(deriv); - psFree(newDeriv); - - if ((min->iter < min->maxIter) || - (min->lastDelta <= min->tol)) { - return(true); - } - - psTrace(".psLib.dataManip.psMinimizeLM", 4, - "---- psMinimizeLM() end (false) ----\n"); - return(false); -} - -/****************************************************************************** -This routine will take an procedure which calculates an arbitrary function -and it's derivative and minimize the chi-squared match between that function -at the specified coords and the specified value at those coords. +psMinimizeLMChi2(): This routine will take an procedure which calculates an +arbitrary function and it's derivative and minimize the chi-squared match +between that function at the specified coords and the specified value at +those coords. XXX: Do this: @@ -1265,7 +389,7 @@ /****************************************************************************** - This routine must fit a polynomial of degree myPoly to the data points - (x, y) and return the coefficients of that polynomial, as well as the - error for each data point (yErr). +psVectorFitPolynomial1D(): This routine must fit a polynomial of degree +myPoly to the data points (x, y) and return the coefficients of that +polynomial, as well as the error for each data point (yErr). XXX: NOTE: yErr is currently ignored. @@ -1273,5 +397,6 @@ psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, const psVector* restrict x, - const psVector* restrict y, const psVector* restrict yErr) + const psVector* restrict y, + const psVector* restrict yErr) { int polyOrder = myPoly->n; @@ -1381,109 +506,4 @@ #define STEP_SIZE 0.10 -/****************************************************************************** - This routine takes as input an arbitrary function, and the parameter to - vary. This function produces as output a bracket [a, b, c] such that - f(b) is less than f(a) and f(b). - - Algorithm: XXX completely ad hoc: start with the user-supplied starting - parameter and call that b. Calculate a/c as a fractional amount - smaller/larger than b. Repeat this process until a local minimum is - found. - *****************************************************************************/ -psVector *p_psDetermineBracketOld(psVector *params, - int dim, - const psArray *coords, - psMinimizePowellFunc func) -{ - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; - int iter = 100; - float aDir = 0.0; - float cDir = 0.0; - float new_aDir = 0.0; - float new_cDir = 0.0; - psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - float stepSize = params->data.F32[dim] * STEP_SIZE; - // float initialParam = params->data.F32[dim]; - - if (stepSize == 0.0) { - stepSize = 1.0; - } - a = b = c = params->data.F32[0]; - a-= stepSize; - c+= stepSize; - - params->data.F32[dim] = a; - fa = func(params, coords); - params->data.F32[dim] = b; - fb = func(params, coords); - params->data.F32[dim] = c; - fc = func(params, coords); - if (fa < fb) { - aDir = -1; - } else { - aDir = 1; - } - - if (fc < fb) { - cDir = -1; - } else { - cDir = 1; - } - - while (iter > 0) { - if ((b < a) && (b < c)) { - bracket->data.F32[0] = a; - bracket->data.F32[1] = b; - bracket->data.F32[2] = c; - return(bracket); - } - stepSize*= (1.0 + stepSize); - a = a - stepSize; - c = c + stepSize; - - params->data.F32[dim] = a; - fa = func(params, coords); - params->data.F32[dim] = c; - fc = func(params, coords); - - //printf("HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc); - - if (fa < fb) { - new_aDir = -1; - } else { - new_aDir = 1; - } - - if (fc < fb) { - new_cDir = -1; - } else { - new_cDir = 1; - } - if ((new_aDir == 1) && (aDir == -1)) { - bracket->data.F32[0] = a; - bracket->data.F32[1] = b; - bracket->data.F32[2] = c; - return(bracket); - } - - if ((new_cDir == 1) && (cDir == -1)) { - bracket->data.F32[0] = a; - bracket->data.F32[1] = b; - bracket->data.F32[2] = c; - return(bracket); - } - aDir = new_aDir; - cDir = new_cDir; - iter--; - } - psFree(bracket); - return(NULL); -} - /****************************************************************************** This routine takes as input an arbitrary function, and the parameter to @@ -1666,87 +686,4 @@ } - -/****************************************************************************** - This routine must minimize a possibly multi-dimensional function - (several parameters) along a single dimension. - *****************************************************************************/ -bool psMinimize1DFunc(psMinimization *min, - psVector *params, - int dim, - const psArray *coords, - psMinimizePowellFunc func) -{ - psVector *bracket; - float a = 0.0; - float b = 0.0; - float c = 0.0; - float n = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; - float fn = 0.0; - // float initialParam = params->data.F32[dim]; - - bracket = p_psDetermineBracketOld(params, dim, coords, func); - if (bracket == NULL) { - psAbort(__func__, "Could not bracket minimum."); - } - - min->iter = 0; - while (min->iter < min->maxIter) { - min->iter++; - //printf("psMinimize1DFunc(): iteration %d\n", min->iter); - a = bracket->data.F32[0]; - b = bracket->data.F32[1]; - c = bracket->data.F32[2]; - - params->data.F32[dim] = a; - fa = func(params, coords); - params->data.F32[dim] = b; - fb = func(params, coords); - params->data.F32[dim] = c; - fc = func(params, coords); - //printf("Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc); - - // We determine which is the biggest segment in [a,b,c] then split - // that with the point n. - if ((b-a) > (c-b)) { - // This is the golden section formula - params->data.F32[dim] = n = a + (0.69 * (b-a)); - fn = func(params, coords); - if (fn > fb) { - // a = n, b = b, c = c - bracket->data.F32[0] = n; - } else { - // a = a, b = n, c = b - bracket->data.F32[1] = n; - bracket->data.F32[2] = b; - } - } else { - params->data.F32[dim] = n = b + (0.69 * (c-b)); - fn = func(params, coords); - if (fn > fb) { - // a = a, b = b, c = n - bracket->data.F32[2] = n; - } else { - // a = b, b = n, c = c - bracket->data.F32[0] = b; - bracket->data.F32[1] = n; - } - } - - if ((fabs(a-b) < min->tol) && - (fabs(b-c) < min->tol)) { - // psFree(bracket); - // XXX: is this line correct? - params->data.F32[dim] = bracket->data.F32[1]; - min->value = func(params, coords); - return(true); - } - } - - // psFree(bracket); - return(false); -} /****************************************************************************** @@ -2168,2 +1105,4 @@ return(psMinimizePowell(min, params, paramMask, coords, myPowellChi2Func)); } + +