Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 1780) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 1781) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-10 02:52:02 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-10 22:25:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -441,4 +441,14 @@ float xSum = 1.0; + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "---- Calling p_psOrdPolynomial1DEval(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "Polynomial order is %d\n", myPoly->n); + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]); + } + + if (NULL == myPoly) { psAbort(__func__, "psPolynomial1DEval(): myPoly is NULL\n"); @@ -455,5 +465,5 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] != 0) { + if (myPoly->mask[loop_x] == 0) { polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -535,5 +545,5 @@ ySum = xSum; for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += ySum * myPoly->coeff[loop_x][loop_y]; ySum *= y; @@ -564,5 +574,5 @@ for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += myPoly->coeff[loop_x][loop_y] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -639,5 +649,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; zSum *= z; @@ -675,5 +685,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_x][loop_y][loop_z] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -768,5 +778,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; zSum *= z; @@ -811,5 +821,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * psPolynomial1DEval(w, chebPolys[loop_w]) * @@ -1129,5 +1139,5 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] != 0) { + if (myPoly->mask[loop_x] == 0) { polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -1210,5 +1220,5 @@ ySum = xSum; for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += ySum * myPoly->coeff[loop_x][loop_y]; ySum *= y; @@ -1239,5 +1249,5 @@ for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += myPoly->coeff[loop_x][loop_y] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -1313,5 +1323,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; zSum *= z; @@ -1349,5 +1359,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_x][loop_y][loop_z] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -1444,5 +1454,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; zSum *= z; @@ -1487,5 +1497,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * psPolynomial1DEval(w, chebPolys[loop_w]) * @@ -1608,9 +1618,12 @@ tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); width = (max - min) / ((float) numSplines); - tmpDomain = min; - for (i=0;idomains)[0] = min; + tmpDomain = min+width; + for (i=1;idomains)[i] = tmpDomain; tmpDomain+= width; } + (tmp->domains)[numSplines] = max; return(tmp); @@ -1661,4 +1674,7 @@ int mid; + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Calling VectorBinDisect(%f)\n", x); + if ((x < bins[0]) || (x > bins[numBins-1])) { @@ -1670,18 +1686,26 @@ min = 0; - max = numBins-1; - mid = (max-min)/2; - - while (min+1 < max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max); - mid = (max-min)/2; - - if (x < bins[mid]) { - max = mid; + max = numBins-2; + mid = ((max+1)-min)/2; + + while (min != max) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", + min, mid, max, x, bins[mid]); + + if (x == bins[mid]) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Exiting VectorBinDisect(): bin %d\n", min); + return(min); + } else if (x < bins[mid]) { + max = mid-1; } else { min = mid; } - } - + mid = ((max+1)+min)/2; + } + + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Exiting VectorBinDisect(): bin %d\n", min); return(min); } @@ -1733,4 +1757,7 @@ float slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Calling psSpline1DGen()\n"); + if (data->n != (1 + spline->n)) { psLogMsg(__func__, PS_LOG_WARN, @@ -1744,8 +1771,16 @@ slope = (data->data.F32[i+1] - data->data.F32[i]) / (spline->domains[i+1] - spline->domains[i]); - (spline->spline[0])->coeff[0] = data->data.F32[i] - + (spline->spline[i])->coeff[0] = data->data.F32[i] - (slope * spline->domains[i]); - (spline->spline[0])->coeff[1] = slope; - } + //printf("HMMM (%f, %f, %f)\n", data->data.F32[i], slope, spline->domains[i]); + + (spline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- spline %d coeffs are (%f, %f)\n", i, + (spline->spline[i])->coeff[0], + (spline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Exiting psSpline1DGen()\n"); return((psSpline1D *) spline); } else if (4 == (spline->spline[0])->n) { Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 1780) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 1781) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-10 02:52:02 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-10 22:25:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -845,5 +845,16 @@ } /****************************************************************************** +This routine will take an procedure which calculates an arbitrary function +and it's derivative and minimize it. + GUS + +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, @@ -857,4 +868,5 @@ psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32); psVector *newParams = 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_F32); @@ -869,7 +881,17 @@ float lamda = 1.0; + for (i=0;in;i++) { + origParams->data.F32[i] = params->data.F32[i]; + } + min->lastDelta = -HUGE; min->iter = 0; min->value = func(deriv, params, coords); + for (i=0;in;i++) { + if ((paramMask != NULL) && + (paramMask->data.U8[i] != 0)) { + deriv->data.F32[i] = 0.0; + } + } while ((min->lastDelta > min->tol) && @@ -897,10 +919,26 @@ newParams = psMatrixLUSolve(newParams, aOut, beta, perm); + for (i=0;in;i++) { + if ((paramMask != NULL) && + (paramMask->data.U8[i] != 0)) { + newParams->data.F32[i] = origParams->data.F32[i]; + } + } + oldValue = min->value; newValue = func(newDeriv, newParams, coords); + 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] = newParams->data.F32[i]; @@ -918,4 +956,5 @@ psFree(perm); psFree(newParams); + psFree(origParams); psFree(alpha); psFree(A); @@ -1769,5 +1808,5 @@ The basic idea is calculate chi-squared for a set of params/coords/errors. -This functions uses globale variables to receive the function pointer, the +This functions uses global variables to receive the function pointer, the data values, and the data errors. *****************************************************************************/ Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 1780) +++ /trunk/psLib/src/math/psMinimize.c (revision 1781) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-10 02:52:02 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-10 22:25:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -845,5 +845,16 @@ } /****************************************************************************** +This routine will take an procedure which calculates an arbitrary function +and it's derivative and minimize it. + GUS + +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, @@ -857,4 +868,5 @@ psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32); psVector *newParams = 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_F32); @@ -869,7 +881,17 @@ float lamda = 1.0; + for (i=0;in;i++) { + origParams->data.F32[i] = params->data.F32[i]; + } + min->lastDelta = -HUGE; min->iter = 0; min->value = func(deriv, params, coords); + for (i=0;in;i++) { + if ((paramMask != NULL) && + (paramMask->data.U8[i] != 0)) { + deriv->data.F32[i] = 0.0; + } + } while ((min->lastDelta > min->tol) && @@ -897,10 +919,26 @@ newParams = psMatrixLUSolve(newParams, aOut, beta, perm); + for (i=0;in;i++) { + if ((paramMask != NULL) && + (paramMask->data.U8[i] != 0)) { + newParams->data.F32[i] = origParams->data.F32[i]; + } + } + oldValue = min->value; newValue = func(newDeriv, newParams, coords); + 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] = newParams->data.F32[i]; @@ -918,4 +956,5 @@ psFree(perm); psFree(newParams); + psFree(origParams); psFree(alpha); psFree(A); @@ -1769,5 +1808,5 @@ The basic idea is calculate chi-squared for a set of params/coords/errors. -This functions uses globale variables to receive the function pointer, the +This functions uses global variables to receive the function pointer, the data values, and the data errors. *****************************************************************************/ Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 1780) +++ /trunk/psLib/src/math/psPolynomial.c (revision 1781) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-10 02:52:02 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-10 22:25:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -441,4 +441,14 @@ float xSum = 1.0; + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "---- Calling p_psOrdPolynomial1DEval(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "Polynomial order is %d\n", myPoly->n); + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]); + } + + if (NULL == myPoly) { psAbort(__func__, "psPolynomial1DEval(): myPoly is NULL\n"); @@ -455,5 +465,5 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] != 0) { + if (myPoly->mask[loop_x] == 0) { polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -535,5 +545,5 @@ ySum = xSum; for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += ySum * myPoly->coeff[loop_x][loop_y]; ySum *= y; @@ -564,5 +574,5 @@ for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += myPoly->coeff[loop_x][loop_y] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -639,5 +649,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; zSum *= z; @@ -675,5 +685,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_x][loop_y][loop_z] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -768,5 +778,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; zSum *= z; @@ -811,5 +821,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * psPolynomial1DEval(w, chebPolys[loop_w]) * @@ -1129,5 +1139,5 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] != 0) { + if (myPoly->mask[loop_x] == 0) { polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -1210,5 +1220,5 @@ ySum = xSum; for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += ySum * myPoly->coeff[loop_x][loop_y]; ySum *= y; @@ -1239,5 +1249,5 @@ for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += myPoly->coeff[loop_x][loop_y] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -1313,5 +1323,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; zSum *= z; @@ -1349,5 +1359,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_x][loop_y][loop_z] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -1444,5 +1454,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; zSum *= z; @@ -1487,5 +1497,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * psPolynomial1DEval(w, chebPolys[loop_w]) * @@ -1608,9 +1618,12 @@ tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); width = (max - min) / ((float) numSplines); - tmpDomain = min; - for (i=0;idomains)[0] = min; + tmpDomain = min+width; + for (i=1;idomains)[i] = tmpDomain; tmpDomain+= width; } + (tmp->domains)[numSplines] = max; return(tmp); @@ -1661,4 +1674,7 @@ int mid; + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Calling VectorBinDisect(%f)\n", x); + if ((x < bins[0]) || (x > bins[numBins-1])) { @@ -1670,18 +1686,26 @@ min = 0; - max = numBins-1; - mid = (max-min)/2; - - while (min+1 < max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max); - mid = (max-min)/2; - - if (x < bins[mid]) { - max = mid; + max = numBins-2; + mid = ((max+1)-min)/2; + + while (min != max) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", + min, mid, max, x, bins[mid]); + + if (x == bins[mid]) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Exiting VectorBinDisect(): bin %d\n", min); + return(min); + } else if (x < bins[mid]) { + max = mid-1; } else { min = mid; } - } - + mid = ((max+1)+min)/2; + } + + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Exiting VectorBinDisect(): bin %d\n", min); return(min); } @@ -1733,4 +1757,7 @@ float slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Calling psSpline1DGen()\n"); + if (data->n != (1 + spline->n)) { psLogMsg(__func__, PS_LOG_WARN, @@ -1744,8 +1771,16 @@ slope = (data->data.F32[i+1] - data->data.F32[i]) / (spline->domains[i+1] - spline->domains[i]); - (spline->spline[0])->coeff[0] = data->data.F32[i] - + (spline->spline[i])->coeff[0] = data->data.F32[i] - (slope * spline->domains[i]); - (spline->spline[0])->coeff[1] = slope; - } + //printf("HMMM (%f, %f, %f)\n", data->data.F32[i], slope, spline->domains[i]); + + (spline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- spline %d coeffs are (%f, %f)\n", i, + (spline->spline[i])->coeff[0], + (spline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Exiting psSpline1DGen()\n"); return((psSpline1D *) spline); } else if (4 == (spline->spline[0])->n) { Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 1780) +++ /trunk/psLib/src/math/psSpline.c (revision 1781) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-10 02:52:02 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-10 22:25:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -441,4 +441,14 @@ float xSum = 1.0; + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "---- Calling p_psOrdPolynomial1DEval(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "Polynomial order is %d\n", myPoly->n); + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4, + "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]); + } + + if (NULL == myPoly) { psAbort(__func__, "psPolynomial1DEval(): myPoly is NULL\n"); @@ -455,5 +465,5 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] != 0) { + if (myPoly->mask[loop_x] == 0) { polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -535,5 +545,5 @@ ySum = xSum; for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += ySum * myPoly->coeff[loop_x][loop_y]; ySum *= y; @@ -564,5 +574,5 @@ for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += myPoly->coeff[loop_x][loop_y] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -639,5 +649,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; zSum *= z; @@ -675,5 +685,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_x][loop_y][loop_z] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -768,5 +778,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; zSum *= z; @@ -811,5 +821,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * psPolynomial1DEval(w, chebPolys[loop_w]) * @@ -1129,5 +1139,5 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { - if (myPoly->mask[loop_x] != 0) { + if (myPoly->mask[loop_x] == 0) { polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -1210,5 +1220,5 @@ ySum = xSum; for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += ySum * myPoly->coeff[loop_x][loop_y]; ySum *= y; @@ -1239,5 +1249,5 @@ for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { - if (myPoly->mask[loop_x][loop_y] != 0) { + if (myPoly->mask[loop_x][loop_y] == 0) { polySum += myPoly->coeff[loop_x][loop_y] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -1313,5 +1323,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; zSum *= z; @@ -1349,5 +1359,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_x][loop_y][loop_z] * psPolynomial1DEval(x, chebPolys[loop_x]) * @@ -1444,5 +1454,5 @@ zSum = ySum; for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; zSum *= z; @@ -1487,5 +1497,5 @@ for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { - if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) { + if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) { polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] * psPolynomial1DEval(w, chebPolys[loop_w]) * @@ -1608,9 +1618,12 @@ tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); width = (max - min) / ((float) numSplines); - tmpDomain = min; - for (i=0;idomains)[0] = min; + tmpDomain = min+width; + for (i=1;idomains)[i] = tmpDomain; tmpDomain+= width; } + (tmp->domains)[numSplines] = max; return(tmp); @@ -1661,4 +1674,7 @@ int mid; + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Calling VectorBinDisect(%f)\n", x); + if ((x < bins[0]) || (x > bins[numBins-1])) { @@ -1670,18 +1686,26 @@ min = 0; - max = numBins-1; - mid = (max-min)/2; - - while (min+1 < max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max); - mid = (max-min)/2; - - if (x < bins[mid]) { - max = mid; + max = numBins-2; + mid = ((max+1)-min)/2; + + while (min != max) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", + min, mid, max, x, bins[mid]); + + if (x == bins[mid]) { + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Exiting VectorBinDisect(): bin %d\n", min); + return(min); + } else if (x < bins[mid]) { + max = mid-1; } else { min = mid; } - } - + mid = ((max+1)+min)/2; + } + + psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, + "---- Exiting VectorBinDisect(): bin %d\n", min); return(min); } @@ -1733,4 +1757,7 @@ float slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Calling psSpline1DGen()\n"); + if (data->n != (1 + spline->n)) { psLogMsg(__func__, PS_LOG_WARN, @@ -1744,8 +1771,16 @@ slope = (data->data.F32[i+1] - data->data.F32[i]) / (spline->domains[i+1] - spline->domains[i]); - (spline->spline[0])->coeff[0] = data->data.F32[i] - + (spline->spline[i])->coeff[0] = data->data.F32[i] - (slope * spline->domains[i]); - (spline->spline[0])->coeff[1] = slope; - } + //printf("HMMM (%f, %f, %f)\n", data->data.F32[i], slope, spline->domains[i]); + + (spline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- spline %d coeffs are (%f, %f)\n", i, + (spline->spline[i])->coeff[0], + (spline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Exiting psSpline1DGen()\n"); return((psSpline1D *) spline); } else if (4 == (spline->spline[0])->n) { Index: /trunk/psLib/test/dataManip/tst_psFunc04.c =================================================================== --- /trunk/psLib/test/dataManip/tst_psFunc04.c (revision 1780) +++ /trunk/psLib/test/dataManip/tst_psFunc04.c (revision 1781) @@ -11,5 +11,5 @@ #include "psFunctions.h" -#define N 10 +#define N 8 int main() @@ -23,5 +23,5 @@ psVector *data; - psTraceSetLevel(".", 9); + psTraceSetLevel(".", 0); /****************************************************************************/ /****************************************************************************/ @@ -33,6 +33,11 @@ tmpSpline = psSpline1DAlloc(N, 1, 1.0, 10.0); + for (i=0;idata.F32[i] = (float) tmpSpline->domains[i]; + printf("tmpSpline->domains[%d] is %f\n", i, tmpSpline->domains[i]); + } + + for (i=0;idata.F32[i] = tmpSpline->domains[i]; }