Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 1921) +++ trunk/psLib/src/math/psMinimize.c (revision 1945) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-28 23:27:37 $ + * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-03 23:35:47 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -682,5 +682,4 @@ psScalar *fScalar; - // XXX: These assignments appear too simple to warrant code and // variable declarations. I retain them here to maintain coherence @@ -714,5 +713,7 @@ fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32, 3, tmpScalar); + f->data.F64[i] = (double) fScalar->data.F32; + psFree(fScalar); psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6, @@ -724,5 +725,7 @@ // coefficients of the Chebyshev polynomial: NR 5.8.7. fac = 2.0/((float) n); - for (j=0;jn;j++) { sum = 0.0; for (k=0;kcoeff[j] = fac * sum; } // XXX: Must free memory. + psFree(f); psFree(x32); psFree(y32); psFree(tmpScalar); - psFree(fScalar); + return(myPoly); } @@ -745,6 +750,4 @@ polynomial of degree myPoly to the data points (x, y) and return the coefficients of that polynomial. - -XXX: yErr is currently ignored. XXX: must add type F32 (currently F64 only). @@ -795,9 +798,9 @@ // Initialize data structures. - for (i = 0; i < (polyOrder); i++) { + for (i = 0; i < polyOrder; i++) { B->data.F64[i] = 0.0; coeffs->data.F64[i] = 0.0; outPerm->data.F64[i] = 0.0; - for (j = 0; j < (polyOrder); j++) { + for (j = 0; j < polyOrder; j++) { A->data.F64[i][j] = 0.0; ALUD->data.F64[i][j] = 0.0; @@ -809,14 +812,32 @@ // Build the B and A data structs. - for (i = 0; i < X->n; i++) { - p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); - - for (k = 0; k < (polyOrder); k++) { - B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k]; - } - - for (k = 0; k < (polyOrder); k++) { - for (j = 0; j < (polyOrder); j++) { - A->data.F64[k][j] += xSums->data.F64[k + j]; + if (yErr == NULL) { + for (i = 0; i < X->n; i++) { + p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); + + for (k = 0; k < polyOrder; k++) { + B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k]; + } + + for (k = 0; k < polyOrder; k++) { + for (j = 0; j < polyOrder; j++) { + A->data.F64[k][j] += xSums->data.F64[k + j]; + } + } + } + } else { + for (i = 0; i < X->n; i++) { + p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); + + for (k = 0; k < polyOrder; k++) { + B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k] / + yErr->data.F64[i]; + } + + for (k = 0; k < polyOrder; k++) { + for (j = 0; j < polyOrder; j++) { + A->data.F64[k][j] += xSums->data.F64[k + j] / + yErr->data.F64[i]; + } } } @@ -826,5 +847,5 @@ coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - for (k = 0; k < (polyOrder); k++) { + for (k = 0; k < polyOrder; k++) { myPoly->coeff[k] = coeffs->data.F64[k]; // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]); @@ -853,6 +874,4 @@ myPoly to the data points (x, y) and return the coefficients of that polynomial. - -XXX: yErr is currently ignored. XXX: must add type F32 (currently F64 only). @@ -870,7 +889,7 @@ psVector *myYErr = NULL; - PS_CHECK_NULL_1DPOLY(myPoly); - PS_CHECK_NULL_VECTOR(y); - PS_CHECK_EMPTY_VECTOR(y); + PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0); + PS_CHECK_NULL_VECTOR_ACTION(y, 0); + PS_CHECK_EMPTY_VECTOR_ACTION(y, 0); // If yErr==NULL, set all errors equal. @@ -880,9 +899,9 @@ if (y->type.type == PS_TYPE_F32) { - for (i=0;in;i++) { + for (i=0;in;i++) { myYErr->data.F32[i] = 1.0; } } else if (y->type.type == PS_TYPE_F64) { - for (i=0;in;i++) { + for (i=0;in;i++) { myYErr->data.F64[i] = 1.0; } @@ -893,5 +912,5 @@ // If x==NULL, create an myX vector with x values set to (0:n), and if - // this is a CHebyshev polynomial, we must scale to (-1:1). + // this is a Chebyshev polynomial, we must scale to (-1:1). // XXX: Verify that this is the correct action. @@ -902,5 +921,5 @@ if (y->type.type == PS_TYPE_F32) { if (myPoly->type == PS_POLYNOMIAL_ORD) { - for (i=0;in;i++) { + for (i=0;in;i++) { myX->data.F32[i] = (float) i; } @@ -909,5 +928,5 @@ float max = (float) (y->n - 1); - for (i=0;in;i++) { + for (i=0;in;i++) { myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) / (0.5 * (max - min)); @@ -916,5 +935,5 @@ } else if (y->type.type == PS_TYPE_F64) { if (myPoly->type == PS_POLYNOMIAL_ORD) { - for (i=0;in;i++) { + for (i=0;in;i++) { myX->data.F64[i] = (float) i; } @@ -923,5 +942,5 @@ double max = (double) (y->n - 1); - for (i=0;in;i++) { + for (i=0;in;i++) { myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) / (0.5 * (max - min)); @@ -1161,15 +1180,19 @@ /****************************************************************************** - This routine must minimize a possibly multi-dimensional function - along a vector defined by line. - -XXX: Use a p_psName(). +This routine takes as input a possibly multi-dimensional function, along +with an initial guess at the parameters of that function and vector "line" +of the same size as the parameter vector. It will minimize the function +along that vector anr returns the offset along that vector at which the +minimum is determined. + +XXX: This routine is not very efficient in terms of total evaluations of the +function. *****************************************************************************/ -float psLineMin(psMinimization *min, - psVector *params, - psVector *line, - const psVector *paramMask, - const psArray *coords, - psMinimizePowellFunc func) +float p_psLineMin(psMinimization *min, + psVector *params, + psVector *line, + const psVector *paramMask, + const psArray *coords, + psMinimizePowellFunc func) { psVector *bracket; @@ -1190,9 +1213,9 @@ int null = 0; - psTrace(".psLib.dataManip.psLineMin", 4, - "---- psLineMin() begin ----\n"); - psTrace(".psLib.dataManip.psLineMin", 6, + psTrace(".psLib.dataManip.p_psLineMin", 4, + "---- p_psLineMin() begin ----\n"); + psTrace(".psLib.dataManip.p_psLineMin", 6, "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.psLineMin", 6, + psTrace(".psLib.dataManip.p_psLineMin", 6, "min->tol is %f\n", min->tol); @@ -1202,5 +1225,5 @@ if (line->data.F32[i] >= FLT_EPSILON) { null = 1; - psTrace(".psLib.dataManip.psLineMin", 4, + psTrace(".psLib.dataManip.p_psLineMin", 4, "line->data.F32[%d] is %f\n", i, line->data.F32[i]); } @@ -1211,5 +1234,5 @@ if (line->data.F32[i] >= FLT_EPSILON) { null = 1; - psTrace(".psLib.dataManip.psLineMin", 4, + psTrace(".psLib.dataManip.p_psLineMin", 4, "line->data.F32[%d] is %f\n", i, line->data.F32[i]); } @@ -1218,6 +1241,6 @@ if (null == 0) { - psTrace(".psLib.dataManip.psLineMin", 2, - "psLineMin() called with zero line vector.\n"); + psTrace(".psLib.dataManip.p_psLineMin", 2, + "p_psLineMin() called with zero line vector.\n"); return(0.0); } @@ -1229,5 +1252,5 @@ for (i=0;in;i++) { - psTrace(".psLib.dataManip.psLineMin", 6, + psTrace(".psLib.dataManip.p_psLineMin", 6, "params, paramMask, line [%d] is (%f %d %f)\n", i, params->data.F32[i], @@ -1248,6 +1271,6 @@ while (min->iter < min->maxIter) { min->iter++; - psTrace(".psLib.dataManip.psLineMin", 6, - "psLineMin(): iteration %d\n", min->iter); + psTrace(".psLib.dataManip.p_psLineMin", 6, + "p_psLineMin(): iteration %d\n", min->iter); a = bracket->data.F32[0]; @@ -1269,5 +1292,5 @@ fb = func(tmpb, coords); fc = func(tmpc, coords); - psTrace(".psLib.dataManip.psLineMin", 6, + psTrace(".psLib.dataManip.p_psLineMin", 6, "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc); @@ -1307,10 +1330,9 @@ } } - psTrace(".psLib.dataManip.psLineMin", 6, + psTrace(".psLib.dataManip.p_psLineMin", 6, "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]); mul = bracket->data.F32[1]; - if ((fabs(a-b) < min->tol) && - (fabs(b-c) < min->tol)) { + if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) { for (i=0;in;i++) { if (paramMask->data.U8[i] == 0) { @@ -1324,6 +1346,6 @@ psFree(tmpc); psFree(tmpn); - psTrace(".psLib.dataManip.psLineMin", 4, - "---- psLineMin() end.a (%f) ----\n", mul); + psTrace(".psLib.dataManip.p_psLineMin", 4, + "---- p_psLineMin() end.a (%f) ----\n", mul); return(mul); } @@ -1336,6 +1358,6 @@ psFree(tmpn); - psTrace(".psLib.dataManip.psLineMin", 4, - "---- psLineMin() end.b (0.0) ----\n"); + psTrace(".psLib.dataManip.p_psLineMin", 4, + "---- p_psLineMin() end.b (0.0) ----\n"); return(0.0); } @@ -1429,5 +1451,5 @@ dummyMin.maxIter = 100; dummyMin.tol = 0.01; - mul = psLineMin(&dummyMin, Q, v[i], paramMask, coords, func); + mul = p_psLineMin(&dummyMin, Q, v[i], paramMask, coords, func); if (fabs(dummyMin.value - currFuncVal) > biggestDiff) { biggestDiff = fabs(dummyMin.value - currFuncVal); @@ -1460,5 +1482,5 @@ } - mul = psLineMin(min, params, u, paramMask, coords, func); + mul = p_psLineMin(min, params, u, paramMask, coords, func); // 6: