Changeset 2202
- Timestamp:
- Oct 26, 2004, 1:14:04 PM (22 years ago)
- Location:
- trunk/psLib
- Files:
-
- 8 edited
-
src/dataManip/psMinimize.c (modified) (3 diffs)
-
src/dataManip/psMinimize.h (modified) (2 diffs)
-
src/math/psMinimize.c (modified) (3 diffs)
-
src/math/psMinimize.h (modified) (2 diffs)
-
src/sys/psTrace.h (modified) (1 diff)
-
src/sysUtils/psTrace.h (modified) (1 diff)
-
test/dataManip/tst_psMinimize06.c (modified) (4 diffs)
-
test/sysUtils/tst_psTrace.c (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataManip/psMinimize.c
r2197 r2202 9 9 * @author GLF, MHPCC 10 10 * 11 * @version $Revision: 1. 59$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-10-26 2 1:24:42$11 * @version $Revision: 1.60 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-10-26 23:14:04 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 393 393 394 394 /****************************************************************************** 395 psMinimizeLMChi2(): This routine will take an procedure which calculates an396 a rbitrary function and it's derivative and minimize the chi-squared match395 psMinimizeLMChi2(): This routine will take an procedure which calculates 396 an arbitrary function and it's derivative and minimize the chi-squared match 397 397 between that function at the specified coords and the specified value at 398 398 those coords. … … 417 417 *****************************************************************************/ 418 418 bool psMinimizeLMChi2(psMinimization *min, 419 psImage *covar, 419 420 psVector *params, 420 421 const psVector *paramMask, 421 psImage *covar,422 const ps Array *coords,423 const psVector * value,422 const psArray *x, 423 const psVector *y, 424 const psVector *yErr, 424 425 psMinimizeLMChi2Func func) 425 {426 PS_CHECK_NULL_PTR_RETURN_NULL(min);427 PS_CHECK_NULL_VECTOR_RETURN_NULL(params);428 PS_CHECK_EMPTY_VECTOR_RETURN_NULL(params);429 PS_CHECK_NULL_PTR_RETURN_NULL(coords);430 PS_CHECK_NULL_VECTOR_RETURN_NULL(value);431 PS_CHECK_EMPTY_VECTOR_RETURN_NULL(value);432 if (paramMask != NULL) {433 PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(params, paramMask);434 }435 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,436 "---- psMinimizeLMChi2() begin ----\n");437 int numData = value->n;438 int numParams = params->n;439 int i;440 int j;441 int k;442 int l;443 int n;444 int p;445 psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);446 psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);447 448 psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);449 psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);450 psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);451 452 psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);453 psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);454 psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);455 456 psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));457 for (i=0;i<numData;i++) {458 deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);459 }460 461 psVector *currValueVec = psVectorAlloc(value->n, PS_TYPE_F32);462 psVector *newValueVec = psVectorAlloc(value->n, PS_TYPE_F32);463 464 float currChi2 = 0.0;465 float newChi2 = 0.0;466 float lamda = 0.00005;467 468 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,469 "min->maxIter is %d\n", min->maxIter);470 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,471 "min->tol is %f\n", min->tol);472 473 for (p=0;p<numParams;p++) {474 origParams->data.F32[p] = params->data.F32[p];475 }476 477 min->lastDelta = HUGE;478 min->iter = 0;479 480 while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {481 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,482 "------------------------------------------------------\n");483 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,484 "Iteration %d. Delta is %f\n", min->iter, min->lastDelta);485 486 //487 // Calculate the current values and chi-squared of the function.488 //489 currChi2 = 0.0;490 for (n=0;n<numData;n++) {491 currValueVec->data.F32[n] = func(deriv[n],492 params,493 (psVector *) coords->data[n]);494 currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);495 }496 497 for (p=0;p<numParams;p++) {498 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,499 "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);500 }501 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,502 "Current chi-squared is (%f)\n", currChi2);503 504 //505 // Mask elements of the derivative for each data point.506 //507 for (p=0;p<numParams;p++) {508 if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {509 for (n=0;n<numData;n++) {510 (deriv[n])->data.F32[p] = 0.0;511 }512 }513 }514 515 //516 // Calculate the BETA vector.517 //518 for (p=0;p<numParams;p++) {519 beta->data.F64[p] = 0.0;520 for (n=0;n<numData;n++) {521 (beta->data.F64[p])+=522 (value->data.F32[n] - currValueVec->data.F32[n]) *523 (deriv[n])->data.F32[p];524 }525 // XXX: multiple by -1 here?526 (beta->data.F64[p])*= -1.0;527 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,528 "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);529 }530 531 //532 // Calculate the ALPHA matrix.533 //534 for (k=0;k<numParams;k++) {535 for (l=0;l<numParams;l++) {536 alpha->data.F32[k][l] = 0.0;537 for (n=0;n<numData;n++) {538 alpha->data.F32[k][l]+= (deriv[n])->data.F32[k] *539 (deriv[n])->data.F32[l];540 }541 }542 }543 544 //545 // Calculate the matrix A.546 //547 for (j=0;j<numParams;j++) {548 for (k=0;k<numParams;k++) {549 if (j == k) {550 A->data.F64[j][k] =551 (double) ((1.0 + lamda) * alpha->data.F32[j][k]);552 } else {553 A->data.F64[j][k] = (double) alpha->data.F32[j][k];554 }555 }556 }557 for (j=0;j<numParams;j++) {558 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");559 for (k=0;k<numParams;k++) {560 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);561 }562 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");563 }564 565 //566 // Solve A * alpha = Beta567 //568 aOut = psMatrixLUD(aOut, perm, A);569 paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);570 571 //572 // Mask any masked parameters.573 //574 for (i=0;i<numParams;i++) {575 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,576 "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);577 if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {578 newParams->data.F32[i] = origParams->data.F32[i];579 } else {580 newParams->data.F32[i] = params->data.F32[i] -581 (float) paramDeltasF64->data.F64[i];582 }583 }584 585 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,586 "Calling func() with new parameters:\n");587 for (i=0;i<numParams;i++) {588 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,589 "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);590 }591 592 593 //594 // Calculate new function values.595 //596 newChi2 = 0.0;597 for (n=0;n<numData;n++) {598 newValueVec->data.F32[n] = func(deriv[n], newParams,599 (psVector *) coords->data[n]);600 newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);601 }602 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,603 "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);604 605 //606 // If the new chi-squared is lower, then keep it.607 //608 if (currChi2 > newChi2) {609 min->lastDelta = currChi2 - newChi2;610 min->value = newChi2;611 612 // We already masked params.613 for (i=0;i<numParams;i++) {614 params->data.F32[i] = (float) newParams->data.F32[i];615 }616 lamda*= 0.1;617 } else {618 lamda*= 10.0;619 }620 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,621 "lamda is %f\n", lamda);622 min->iter++;623 }624 psFree(beta);625 psFree(perm);626 psFree(paramDeltasF64);627 psFree(origParams);628 psFree(newParams);629 psFree(alpha);630 psFree(A);631 psFree(aOut);632 for (i=0;i<numData;i++) {633 psFree(deriv[i]);634 }635 psFree(deriv);636 psFree(currValueVec);637 psFree(newValueVec);638 639 if ((min->iter < min->maxIter) ||640 (min->lastDelta <= min->tol)) {641 return(true);642 }643 644 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,645 "---- psMinimizeLMChi2() end (false) ----\n");646 return(false);647 }648 649 /******************************************************************************650 psMyMinimizeLMChi2(): This routine will take an procedure which calculates651 an arbitrary function and it's derivative and minimize the chi-squared match652 between that function at the specified coords and the specified value at653 those coords.654 655 XXX: Do this:656 After checking that all entries in the paramMask are 1 or 0, when657 forming the A matrix from alpha, try this:658 659 A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];660 661 XXX: This is very different from what is specified in the SDR. Must662 coordinate with IfA on new SDR.663 664 XXX: Do vector/image recycles.665 666 XXX: probably yErr will be part of the SDR.667 668 XXX: This must work for both F32 and F64. F32 is currently implemented.669 Note: since the LUD routines are only implemented in F64, then we670 will have to convert all F32 input vectors to F64 regardless. So,671 the F64 port might be.672 *****************************************************************************/673 bool psMyMinimizeLMChi2(psMinimization *min,674 psImage *covar,675 psVector *params,676 const psVector *paramMask,677 const psArray *x,678 const psVector *y,679 const psVector *yErr,680 psMyMinimizeLMChi2Func func)681 426 { 682 427 PS_CHECK_NULL_PTR_RETURN_NULL(min); -
trunk/psLib/src/dataManip/psMinimize.h
r2197 r2202 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.3 0$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-10-26 2 1:24:43$10 * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2004-10-26 23:14:04 $ 12 12 * 13 13 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 83 83 ); 84 84 85 86 /*87 Bug 203:88 89 typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params,90 const psVector *x);91 92 This function takes the current guess for the parameters for which we are93 trying to get the best values (params), and a single vector (x) of94 conditions from the array of x values fed into the minimiser.95 96 bool psMinimizeLMChi2(psMinimization *min,97 psImage *covar,98 psVector *params,99 const psVector *paramMask,100 const psArray *x,101 const psVector *y,102 const psVector *yErr,103 psMinimizeLMChi2Func func);104 105 This takes the minimization specs (min), returns the covariance matrix106 (covar), takes the best guess of initial parameters (params), the parameter107 mask (paramMask), and takes multiple vectors of conditions in an array (x),108 the corresponding measured values (y) and errors (yErr) and the function to109 fit (func).110 111 For example, for GRB afterglows, I have flux as a function of time and112 frequency. So I stuff into the "psArray *x" all my time and frequency values113 (so a whole heap of vectors of size 2), I have the measured values in "y" and114 errors in "yErr". Then each of the time-frequency pairs are passed to my model115 function with the current parameters, and the model function returns the flux116 for that time-frequency pair, and the derivative with respect to each of the117 parameters.118 119 This seems reasonable. The only thing we could change would be to have the120 function be defined:121 122 typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,123 const psVector *params,124 const psArray *x);125 126 So it would return the model value for each of the measurements at once, and127 return for each the derivatives (so that it returns a matrix).128 129 What do you think?130 131 132 133 I'm not sure I understand how LM chi-squared minimization will work with the134 following function that you define:135 136 typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,137 const psVector *params,138 const psArray *x);139 140 The chi-squared minimization algorithm, as defined in NR, requires that141 function to be minimized be evaluated at each data point, and that all142 derivatives, with respect to each parameter, be calculated at each data143 point. In the above, can I assume that144 145 x is an array of psVectors, with each vector corresponding to a single146 data point.147 148 The returned value has the same length as x. It contains the value of149 the function at each data point in x.150 151 deriv: an n-by-p matrix where "n" is the number of data points, and "p"152 is the number of parameters. The [i][j] element of this matrix153 holds the derivative of the function at the i-th data point with154 respect to the j-th parameter.155 156 157 158 */159 160 161 // XXX: What if any of these arguments are NULL?162 163 164 85 typedef 165 float (*psMinimizeLMChi2Func)(psVector*deriv,166 const psVector *params,167 const psVector *coords);86 psVector* (*psMinimizeLMChi2Func)(psImage *deriv, 87 const psVector *params, 88 const psArray *x); 168 89 169 90 bool psMinimizeLMChi2(psMinimization *min, 91 psImage *covar, 170 92 psVector *params, 171 93 const psVector *paramMask, 172 psImage *covar,173 const ps Array *coords,174 const psVector * value,94 const psArray *x, 95 const psVector *y, 96 const psVector *yErr, 175 97 psMinimizeLMChi2Func func); 176 177 178 typedef179 psVector* (*psMyMinimizeLMChi2Func)(psImage *deriv,180 const psVector *params,181 const psArray *x);182 183 bool psMyMinimizeLMChi2(psMinimization *min,184 psImage *covar,185 psVector *params,186 const psVector *paramMask,187 const psArray *x,188 const psVector *y,189 const psVector *yErr,190 psMyMinimizeLMChi2Func func);191 98 192 99 typedef -
trunk/psLib/src/math/psMinimize.c
r2197 r2202 9 9 * @author GLF, MHPCC 10 10 * 11 * @version $Revision: 1. 59$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-10-26 2 1:24:42$11 * @version $Revision: 1.60 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-10-26 23:14:04 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 393 393 394 394 /****************************************************************************** 395 psMinimizeLMChi2(): This routine will take an procedure which calculates an396 a rbitrary function and it's derivative and minimize the chi-squared match395 psMinimizeLMChi2(): This routine will take an procedure which calculates 396 an arbitrary function and it's derivative and minimize the chi-squared match 397 397 between that function at the specified coords and the specified value at 398 398 those coords. … … 417 417 *****************************************************************************/ 418 418 bool psMinimizeLMChi2(psMinimization *min, 419 psImage *covar, 419 420 psVector *params, 420 421 const psVector *paramMask, 421 psImage *covar,422 const ps Array *coords,423 const psVector * value,422 const psArray *x, 423 const psVector *y, 424 const psVector *yErr, 424 425 psMinimizeLMChi2Func func) 425 {426 PS_CHECK_NULL_PTR_RETURN_NULL(min);427 PS_CHECK_NULL_VECTOR_RETURN_NULL(params);428 PS_CHECK_EMPTY_VECTOR_RETURN_NULL(params);429 PS_CHECK_NULL_PTR_RETURN_NULL(coords);430 PS_CHECK_NULL_VECTOR_RETURN_NULL(value);431 PS_CHECK_EMPTY_VECTOR_RETURN_NULL(value);432 if (paramMask != NULL) {433 PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(params, paramMask);434 }435 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,436 "---- psMinimizeLMChi2() begin ----\n");437 int numData = value->n;438 int numParams = params->n;439 int i;440 int j;441 int k;442 int l;443 int n;444 int p;445 psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);446 psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);447 448 psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);449 psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);450 psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);451 452 psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);453 psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);454 psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);455 456 psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));457 for (i=0;i<numData;i++) {458 deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);459 }460 461 psVector *currValueVec = psVectorAlloc(value->n, PS_TYPE_F32);462 psVector *newValueVec = psVectorAlloc(value->n, PS_TYPE_F32);463 464 float currChi2 = 0.0;465 float newChi2 = 0.0;466 float lamda = 0.00005;467 468 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,469 "min->maxIter is %d\n", min->maxIter);470 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,471 "min->tol is %f\n", min->tol);472 473 for (p=0;p<numParams;p++) {474 origParams->data.F32[p] = params->data.F32[p];475 }476 477 min->lastDelta = HUGE;478 min->iter = 0;479 480 while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {481 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,482 "------------------------------------------------------\n");483 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,484 "Iteration %d. Delta is %f\n", min->iter, min->lastDelta);485 486 //487 // Calculate the current values and chi-squared of the function.488 //489 currChi2 = 0.0;490 for (n=0;n<numData;n++) {491 currValueVec->data.F32[n] = func(deriv[n],492 params,493 (psVector *) coords->data[n]);494 currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);495 }496 497 for (p=0;p<numParams;p++) {498 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,499 "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);500 }501 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,502 "Current chi-squared is (%f)\n", currChi2);503 504 //505 // Mask elements of the derivative for each data point.506 //507 for (p=0;p<numParams;p++) {508 if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {509 for (n=0;n<numData;n++) {510 (deriv[n])->data.F32[p] = 0.0;511 }512 }513 }514 515 //516 // Calculate the BETA vector.517 //518 for (p=0;p<numParams;p++) {519 beta->data.F64[p] = 0.0;520 for (n=0;n<numData;n++) {521 (beta->data.F64[p])+=522 (value->data.F32[n] - currValueVec->data.F32[n]) *523 (deriv[n])->data.F32[p];524 }525 // XXX: multiple by -1 here?526 (beta->data.F64[p])*= -1.0;527 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,528 "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);529 }530 531 //532 // Calculate the ALPHA matrix.533 //534 for (k=0;k<numParams;k++) {535 for (l=0;l<numParams;l++) {536 alpha->data.F32[k][l] = 0.0;537 for (n=0;n<numData;n++) {538 alpha->data.F32[k][l]+= (deriv[n])->data.F32[k] *539 (deriv[n])->data.F32[l];540 }541 }542 }543 544 //545 // Calculate the matrix A.546 //547 for (j=0;j<numParams;j++) {548 for (k=0;k<numParams;k++) {549 if (j == k) {550 A->data.F64[j][k] =551 (double) ((1.0 + lamda) * alpha->data.F32[j][k]);552 } else {553 A->data.F64[j][k] = (double) alpha->data.F32[j][k];554 }555 }556 }557 for (j=0;j<numParams;j++) {558 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");559 for (k=0;k<numParams;k++) {560 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);561 }562 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");563 }564 565 //566 // Solve A * alpha = Beta567 //568 aOut = psMatrixLUD(aOut, perm, A);569 paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);570 571 //572 // Mask any masked parameters.573 //574 for (i=0;i<numParams;i++) {575 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,576 "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);577 if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {578 newParams->data.F32[i] = origParams->data.F32[i];579 } else {580 newParams->data.F32[i] = params->data.F32[i] -581 (float) paramDeltasF64->data.F64[i];582 }583 }584 585 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,586 "Calling func() with new parameters:\n");587 for (i=0;i<numParams;i++) {588 psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,589 "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);590 }591 592 593 //594 // Calculate new function values.595 //596 newChi2 = 0.0;597 for (n=0;n<numData;n++) {598 newValueVec->data.F32[n] = func(deriv[n], newParams,599 (psVector *) coords->data[n]);600 newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);601 }602 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,603 "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);604 605 //606 // If the new chi-squared is lower, then keep it.607 //608 if (currChi2 > newChi2) {609 min->lastDelta = currChi2 - newChi2;610 min->value = newChi2;611 612 // We already masked params.613 for (i=0;i<numParams;i++) {614 params->data.F32[i] = (float) newParams->data.F32[i];615 }616 lamda*= 0.1;617 } else {618 lamda*= 10.0;619 }620 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,621 "lamda is %f\n", lamda);622 min->iter++;623 }624 psFree(beta);625 psFree(perm);626 psFree(paramDeltasF64);627 psFree(origParams);628 psFree(newParams);629 psFree(alpha);630 psFree(A);631 psFree(aOut);632 for (i=0;i<numData;i++) {633 psFree(deriv[i]);634 }635 psFree(deriv);636 psFree(currValueVec);637 psFree(newValueVec);638 639 if ((min->iter < min->maxIter) ||640 (min->lastDelta <= min->tol)) {641 return(true);642 }643 644 psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,645 "---- psMinimizeLMChi2() end (false) ----\n");646 return(false);647 }648 649 /******************************************************************************650 psMyMinimizeLMChi2(): This routine will take an procedure which calculates651 an arbitrary function and it's derivative and minimize the chi-squared match652 between that function at the specified coords and the specified value at653 those coords.654 655 XXX: Do this:656 After checking that all entries in the paramMask are 1 or 0, when657 forming the A matrix from alpha, try this:658 659 A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];660 661 XXX: This is very different from what is specified in the SDR. Must662 coordinate with IfA on new SDR.663 664 XXX: Do vector/image recycles.665 666 XXX: probably yErr will be part of the SDR.667 668 XXX: This must work for both F32 and F64. F32 is currently implemented.669 Note: since the LUD routines are only implemented in F64, then we670 will have to convert all F32 input vectors to F64 regardless. So,671 the F64 port might be.672 *****************************************************************************/673 bool psMyMinimizeLMChi2(psMinimization *min,674 psImage *covar,675 psVector *params,676 const psVector *paramMask,677 const psArray *x,678 const psVector *y,679 const psVector *yErr,680 psMyMinimizeLMChi2Func func)681 426 { 682 427 PS_CHECK_NULL_PTR_RETURN_NULL(min); -
trunk/psLib/src/math/psMinimize.h
r2197 r2202 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.3 0$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-10-26 2 1:24:43$10 * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2004-10-26 23:14:04 $ 12 12 * 13 13 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 83 83 ); 84 84 85 86 /*87 Bug 203:88 89 typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params,90 const psVector *x);91 92 This function takes the current guess for the parameters for which we are93 trying to get the best values (params), and a single vector (x) of94 conditions from the array of x values fed into the minimiser.95 96 bool psMinimizeLMChi2(psMinimization *min,97 psImage *covar,98 psVector *params,99 const psVector *paramMask,100 const psArray *x,101 const psVector *y,102 const psVector *yErr,103 psMinimizeLMChi2Func func);104 105 This takes the minimization specs (min), returns the covariance matrix106 (covar), takes the best guess of initial parameters (params), the parameter107 mask (paramMask), and takes multiple vectors of conditions in an array (x),108 the corresponding measured values (y) and errors (yErr) and the function to109 fit (func).110 111 For example, for GRB afterglows, I have flux as a function of time and112 frequency. So I stuff into the "psArray *x" all my time and frequency values113 (so a whole heap of vectors of size 2), I have the measured values in "y" and114 errors in "yErr". Then each of the time-frequency pairs are passed to my model115 function with the current parameters, and the model function returns the flux116 for that time-frequency pair, and the derivative with respect to each of the117 parameters.118 119 This seems reasonable. The only thing we could change would be to have the120 function be defined:121 122 typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,123 const psVector *params,124 const psArray *x);125 126 So it would return the model value for each of the measurements at once, and127 return for each the derivatives (so that it returns a matrix).128 129 What do you think?130 131 132 133 I'm not sure I understand how LM chi-squared minimization will work with the134 following function that you define:135 136 typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,137 const psVector *params,138 const psArray *x);139 140 The chi-squared minimization algorithm, as defined in NR, requires that141 function to be minimized be evaluated at each data point, and that all142 derivatives, with respect to each parameter, be calculated at each data143 point. In the above, can I assume that144 145 x is an array of psVectors, with each vector corresponding to a single146 data point.147 148 The returned value has the same length as x. It contains the value of149 the function at each data point in x.150 151 deriv: an n-by-p matrix where "n" is the number of data points, and "p"152 is the number of parameters. The [i][j] element of this matrix153 holds the derivative of the function at the i-th data point with154 respect to the j-th parameter.155 156 157 158 */159 160 161 // XXX: What if any of these arguments are NULL?162 163 164 85 typedef 165 float (*psMinimizeLMChi2Func)(psVector*deriv,166 const psVector *params,167 const psVector *coords);86 psVector* (*psMinimizeLMChi2Func)(psImage *deriv, 87 const psVector *params, 88 const psArray *x); 168 89 169 90 bool psMinimizeLMChi2(psMinimization *min, 91 psImage *covar, 170 92 psVector *params, 171 93 const psVector *paramMask, 172 psImage *covar,173 const ps Array *coords,174 const psVector * value,94 const psArray *x, 95 const psVector *y, 96 const psVector *yErr, 175 97 psMinimizeLMChi2Func func); 176 177 178 typedef179 psVector* (*psMyMinimizeLMChi2Func)(psImage *deriv,180 const psVector *params,181 const psArray *x);182 183 bool psMyMinimizeLMChi2(psMinimization *min,184 psImage *covar,185 psVector *params,186 const psVector *paramMask,187 const psArray *x,188 const psVector *y,189 const psVector *yErr,190 psMyMinimizeLMChi2Func func);191 98 192 99 typedef -
trunk/psLib/src/sys/psTrace.h
r1834 r2202 7 7 * 8 8 * @author Robert Lupton, Princeton University 9 * @author G eorge Gusciora, MHPCC9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $12 * @date $Date: 2004- 09-20 22:03:35$11 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-10-26 23:14:04 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii -
trunk/psLib/src/sysUtils/psTrace.h
r1834 r2202 7 7 * 8 8 * @author Robert Lupton, Princeton University 9 * @author G eorge Gusciora, MHPCC9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $12 * @date $Date: 2004- 09-20 22:03:35$11 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-10-26 23:14:04 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii -
trunk/psLib/test/dataManip/tst_psMinimize06.c
r2197 r2202 8 8 #define NUM_ITERATIONS 10000 9 9 #define ERR_TOL 0.0 10 #define N 510 #define N 20 11 11 #define MIN_VALUE 5.0 12 12 #define NUM_PARAMS 3 … … 24 24 25 25 *****************************************************************************/ 26 float myFunc(psVector *myDeriv, 27 psVector *myParams, 28 psVector *myCoords) 29 { 30 float sum = 0.0; 31 // float coordData = 0.0; 32 // float expData = 0.0; 33 int i; 34 35 if (myDeriv == NULL) { 36 myDeriv = psVectorAlloc(myParams->n, PS_TYPE_F32); 37 psError(__func__, "myDeriv is NULL.\n"); 38 } 39 40 // Simply test that coords were passed in correctly. 41 /* 42 for (i=0;i<N;i++) { 43 coordData = myCoords->data.F32[0]; 44 expData = (float) (i+10); 45 if (fabs(coordData - expData) > FLT_EPSILON) { 46 printf("ERROR(1): coordinate data was incorrectly passed to myFunc()\n"); 47 printf("ERROR(1): was (%f) should be (%f)\n", coordData, expData); 48 testStatus = false; 49 } 50 coordData = myCoords->data.F32[1]; 51 expData = (float) (i+3); 52 if (fabs(coordData - expData) > FLT_EPSILON) { 53 printf("ERROR(2): coordinate data was incorrectly passed to myFunc()\n"); 54 printf("ERROR(2): was (%f) should be (%f)\n", coordData, expData); 55 testStatus = false; 56 } 57 } 58 */ 59 60 sum = 0.0; 61 for (i=0;i<NUM_PARAMS;i++) { 62 sum+= (myParams->data.F32[i] - expectedParm[i]) * 63 (myParams->data.F32[i] - expectedParm[i]); 64 myDeriv->data.F32[i] = (2.0 * myParams->data.F32[i]) - 65 (2.0 * expectedParm[i]); 66 } 67 // for (i=0;i<NUM_PARAMS;i++) 68 // printf("HMMM: myParams->data.F32[%d] is %f\n", i, myParams->data.F32[i]); 69 // for (i=0;i<NUM_PARAMS;i++) 70 // printf("HMMM: myDeriv->data.F32[%d] is %f\n", i, myDeriv->data.F32[i]); 71 72 sum+= MIN_VALUE; 73 return(sum); 74 } 75 76 psVector *myFunc2(psImage *myDeriv, 77 psVector *myParams, 78 psArray *myCoords) 26 psVector *myFunc(psImage *myDeriv, 27 psVector *myParams, 28 psArray *myCoords) 79 29 { 80 30 psVector *sum = psVectorAlloc(myCoords->n, PS_TYPE_F32); … … 136 86 } 137 87 138 psM yMinimizeLMChi2(min,139 140 141 142 143 144 145 (psMyMinimizeLMChi2Func) myFunc2);88 psMinimizeLMChi2(min, 89 myCovar, 90 myParams, 91 NULL, 92 myCoords, 93 y, 94 NULL, 95 (psMinimizeLMChi2Func) myFunc); 146 96 147 97 printf("\nThe chi-squared is %f\n", min->value); … … 165 115 } 166 116 167 int t02()168 {169 int currentId = psMemGetId();170 int memLeaks = 0;171 int i = 0;172 psArray *myCoords;173 psVector *myParams;174 psVector *myParamMask;175 psImage *myCovar;176 psMinimization *min;177 psVector *y = psVectorAlloc(N, PS_TYPE_F32);178 179 psTraceSetLevel(".psLib", 0);180 t02();181 /**************************************************************************182 *************************************************************************/183 myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);184 myParamMask = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);185 min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL);186 myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);187 188 myCoords = psArrayAlloc(N);189 for (i=0;i<N;i++) {190 myCoords->data[i] = (psPTR *) psVectorAlloc(2, PS_TYPE_F32);191 ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) (i+10);192 ((psVector *) (myCoords->data[i]))->data.F32[1] = (float) (i+3);193 y->data.F32[i] = (float) i;194 }195 for (i=0;i<NUM_PARAMS;i++) {196 expectedParm[i] = 2.42 + (float) (2 * i);197 myParams->data.F32[i] = (float) i;198 myParams->data.F32[i] = expectedParm[i] * 1.3;199 myParams->data.F32[i] = (float) (5 + i);200 myParams->data.F32[i] = 0.0;201 myParamMask->data.U8[i] = 0;202 }203 204 psMinimizeLMChi2(min,205 myParams,206 NULL,207 myCovar,208 myCoords,209 y,210 (psMinimizeLMChi2Func) myFunc);211 212 printf("\nThe chi-squared is %f\n", min->value);213 for (i=0;i<NUM_PARAMS;i++) {214 printf("Parameter %d at the minimum is %f (expected: %f)\n", i,215 myParams->data.F32[i], expectedParm[i]);216 }217 218 psFree(myCoords);219 psFree(myParams);220 psFree(myParamMask);221 psFree(min);222 psFree(y);223 psFree(myCovar);224 225 psMemCheckCorruption(1);226 memLeaks = psMemCheckLeaks(currentId,NULL,NULL);227 if (0 != memLeaks) {228 psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);229 }230 return (!testStatus);231 }232 233 117 int main() 234 118 { 235 119 t01(); 236 // t02();237 120 } -
trunk/psLib/test/sysUtils/tst_psTrace.c
r1835 r2202 23 23 {testTrace00, 0, "psTraceSetLevel() and psTraceGetLevel()", 0, false}, 24 24 {testTrace01, 1, "psTraceSetLevel(): set multiple components in one call", 0, false}, 25 {testTrace02, 2, "psTraceSetLevel(): test static inheritance", 0, false},25 {testTrace02, 2, "psTraceSetLevel(): test static/dynamic inheritance", 0, false}, 26 26 {testTrace03, 3, "psTraceReset()", 0, false}, 27 27 {testTrace04, 4, "psTrace()", 0, false},
Note:
See TracChangeset
for help on using the changeset viewer.
