Changeset 7720 for trunk/projects/slimcurve/src
 Timestamp:
 06/01/11 20:17:05 (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/projects/slimcurve/src/main/c/EcfGlobal.c
r7706 r7720 79 79 float *exp_conv[], 80 80 float yfit[], float dy[], 81 float **alpha, float *beta, float *chisq );81 float **alpha, float *beta, float *chisq, float old_chisq); 82 82 int GCI_marquardt_global_compute_exps_fn_final( 83 83 float xincr, float y[], … … 1164 1164 xincr, y, ndata, fit_start, fit_end, noise, sig, 1165 1165 ftype, param, paramfree, nparam, exp_conv, 1166 yfit, dy, alpha, beta, chisq ) != 0)1166 yfit, dy, alpha, beta, chisq, 0.0) != 0) 1167 1167 return 2; 1168 1168 … … 1219 1219 xincr, y, ndata, fit_start, fit_end, noise, sig, 1220 1220 ftype, paramtry, paramfree, nparam, exp_conv, 1221 yfit, dy, covar, dparam, chisq ) != 0)1221 yfit, dy, covar, dparam, chisq, ochisq) != 0) 1222 1222 return 2; 1223 1223 … … 1250 1250 float *exp_conv[], 1251 1251 float yfit[], float dy[], 1252 float **alpha, float *beta, float *chisq) 1252 float **alpha, float *beta, 1253 float *chisq, float old_chisq) 1253 1254 { 1254 1255 int i, j, k, l, m, mfit; 1255 float wt, sig2i, y_ymod, dy_dparam[MAXFIT]; 1256 float wt, sig2i, y_ymod; 1257 float dy_dparam[MAXBINS][MAXFIT]; 1258 float alpha_weight[MAXBINS]; 1259 float beta_weight[MAXBINS]; 1260 float weight; 1261 int i_free; 1262 int j_free; 1263 float dot_product; 1264 float beta_sum; 1265 float dy_dparam_k_i; 1256 1266 1257 1267 for (j=0, mfit=0; j<nparam; j++) … … 1259 1269 mfit++; 1260 1270 1261 /* Initialise (symmetric) alpha, beta */ 1262 for (j=0; j<mfit; j++) { 1263 for (k=0; k<=j; k++) 1264 alpha[j][k] = 0.0; 1265 beta[j] = 0.0; 1266 } 1267 1268 /* Calculation of the fitting data will depend upon the type of 1269 noise. Since there's no convolution involved here, this is 1270 very easy. */ 1271 *chisq = 0.0; 1271 1272 1272 1273 switch (ftype) { 1273 case FIT_GLOBAL_MULTIEXP: 1274 switch (noise) { 1275 case NOISE_CONST: 1276 *chisq = 0.0; 1277 /* Summation loop over all data */ 1278 for (i=fit_start; i<fit_end; i++) { 1279 yfit[i] = param[0]; /* Z */ 1280 dy_dparam[0] = 1.0; 1281 1282 for (j=1; j<nparam; j+=2) { 1283 yfit[i] += param[j] * exp_conv[j+1][i]; 1284 /* A_j . exp(x/tau_j) */ 1285 dy_dparam[j] = exp_conv[j+1][i]; 1286 /* exp(x/tau_j) */ 1287 dy_dparam[j+1] = param[j] * exp_conv[j][i]; 1288 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1274 case FIT_GLOBAL_MULTIEXP: 1275 switch (noise) { 1276 case NOISE_CONST: 1277 // loop over all data 1278 for (i = fit_start; i < fit_end; ++i) { 1279 // multiexponential fit 1280 yfit[i] = param[0]; /* Z */ 1281 dy_dparam[i][0] = 1.0; 1282 1283 for (j=1; j<nparam; j+=2) { 1284 yfit[i] += param[j] * exp_conv[j+1][i]; 1285 /* A_j . exp(x/tau_j) */ 1286 dy_dparam[i][j] = exp_conv[j+1][i]; 1287 /* exp(x/tau_j) */ 1288 dy_dparam[i][j+1] = param[j] * exp_conv[j][i]; 1289 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1290 } 1291 dy[i] = y[i]  yfit[i]; 1292 1293 // constant noise 1294 weight = 1.0f / sig[0]; 1295 alpha_weight[i] = weight; // 1 / (sig[0] * sig[0]); 1296 weight *= dy[i]; 1297 beta_weight[i] = weight; // dy[i] / (sig[0] * sig[0]); 1298 weight *= dy[i]; 1299 *chisq += weight; // (dy[i] * dy[i]) / (sig[0] * sig[0]); 1300 } 1301 break; 1302 case NOISE_GIVEN: 1303 // loop over all data 1304 for (i = fit_start; i < fit_end; ++i) { 1305 // multiexponential fit 1306 yfit[i] = param[0]; /* Z */ 1307 dy_dparam[i][0] = 1.0; 1308 1309 for (j=1; j<nparam; j+=2) { 1310 yfit[i] += param[j] * exp_conv[j+1][i]; 1311 /* A_j . exp(x/tau_j) */ 1312 dy_dparam[i][j] = exp_conv[j+1][i]; 1313 /* exp(x/tau_j) */ 1314 dy_dparam[i][j+1] = param[j] * exp_conv[j][i]; 1315 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1316 } 1317 dy[i] = y[i]  yfit[i]; 1318 1319 // given noise 1320 weight = 1.0f / (sig[i] * sig[i]); 1321 alpha_weight[i] = weight; // 1 / (sig[i] * sig[i]) 1322 weight *= dy[i]; 1323 beta_weight[i] = weight; // dy[i] / (sig[i] * sig[i]) 1324 weight *= dy[i]; 1325 *chisq += weight; // (dy[i] * dy[i]) / (sig[i] * sig[i]) 1326 } 1327 break; 1328 case NOISE_POISSON_DATA: 1329 // loop over all data 1330 for (i = fit_start; i < fit_end; ++i) { 1331 // multiexponential fit 1332 yfit[i] = param[0]; /* Z */ 1333 dy_dparam[i][0] = 1.0; 1334 1335 for (j=1; j<nparam; j+=2) { 1336 yfit[i] += param[j] * exp_conv[j+1][i]; 1337 /* A_j . exp(x/tau_j) */ 1338 dy_dparam[i][j] = exp_conv[j+1][i]; 1339 /* exp(x/tau_j) */ 1340 dy_dparam[i][j+1] = param[j] * exp_conv[j][i]; 1341 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1342 } 1343 dy[i] = y[i]  yfit[i]; 1344 1345 // poisson noise based on data 1346 weight = (y[i] > 15 ? 1.0f / y[i] : 1.0f / 15); 1347 alpha_weight[i] = weight; // 1 / sig(i) 1348 weight *= dy[i]; 1349 beta_weight[i] = weight; // dy[i] / sig(i) 1350 weight *= dy[i]; 1351 *chisq += weight; // (dy[i] * dy[i]) / sig(i) 1352 } 1353 break; 1354 case NOISE_POISSON_FIT: 1355 // loop over all data 1356 for (i = fit_start; i < fit_end; ++i) { 1357 // multiexponential fit 1358 yfit[i] = param[0]; /* Z */ 1359 dy_dparam[i][0] = 1.0; 1360 1361 for (j=1; j<nparam; j+=2) { 1362 yfit[i] += param[j] * exp_conv[j+1][i]; 1363 /* A_j . exp(x/tau_j) */ 1364 dy_dparam[i][j] = exp_conv[j+1][i]; 1365 /* exp(x/tau_j) */ 1366 dy_dparam[i][j+1] = param[j] * exp_conv[j][i]; 1367 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1368 } 1369 dy[i] = y[i]  yfit[i]; 1370 1371 // poisson noise based on fit 1372 weight = (yfit[i] > 15 ? 1.0f / yfit[i] : 1.0f / 15); 1373 alpha_weight[i] = weight; // 1 / sig(i) 1374 weight *= dy[i]; 1375 beta_weight[i] = weight; // dy(i) / sig(i) 1376 weight *= dy[i]; 1377 *chisq += weight; // (dy(i) * dy(i)) / sig(i) 1378 } 1379 break; 1380 case NOISE_GAUSSIAN_FIT: 1381 // loop over all data 1382 for (i = fit_start; i < fit_end; ++i) { 1383 // multiexponential fit 1384 yfit[i] = param[0]; /* Z */ 1385 dy_dparam[i][0] = 1.0; 1386 1387 for (j=1; j<nparam; j+=2) { 1388 yfit[i] += param[j] * exp_conv[j+1][i]; 1389 /* A_j . exp(x/tau_j) */ 1390 dy_dparam[i][j] = exp_conv[j+1][i]; 1391 /* exp(x/tau_j) */ 1392 dy_dparam[i][j+1] = param[j] * exp_conv[j][i]; 1393 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1394 } 1395 dy[i] = y[i]  yfit[i]; 1396 1397 // gaussian noise based on fit 1398 weight = (yfit[i] > 1.0f ? 1.0f / yfit[i] : 1.0f); 1399 alpha_weight[i] = weight; // 1 / sig(i) 1400 weight *= dy[i]; 1401 beta_weight[i] = weight; // dy[i] / sig(i) 1402 weight *= dy[i]; 1403 *chisq += weight; // dy[i] / sig(i) 1404 } 1405 break; 1406 case NOISE_MLE: 1407 // loop over all data 1408 for (i = fit_start; i < fit_end; ++i) { 1409 // multiexponential fit 1410 yfit[i] = param[0]; /* Z */ 1411 dy_dparam[i][0] = 1.0; 1412 1413 for (j=1; j<nparam; j+=2) { 1414 yfit[i] += param[j] * exp_conv[j+1][i]; 1415 /* A_j . exp(x/tau_j) */ 1416 dy_dparam[i][j] = exp_conv[j+1][i]; 1417 /* exp(x/tau_j) */ 1418 dy_dparam[i][j+1] = param[j] * exp_conv[j][i]; 1419 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1420 } 1421 dy[i] = y[i]  yfit[i]; 1422 1423 1424 // maximum likelihood estimation noise 1425 weight = (yfit[i] > 1 ? 1.0f / yfit[i] : 1.0f); 1426 alpha_weight[i] = weight * y[i] / yfit[i]; 1427 beta_weight[i] = dy[i] * weight; 1428 if (yfit[i] > 0.0) { 1429 *chisq += (0.0f == y[i]) 1430 ? 2.0 * yfit[i] 1431 : 2.0 * (yfit[i]  y[i])  2.0 * y[i] * log(yfit[i] / y[i]); 1432 } 1433 } 1434 if (*chisq <= 0.0f) { 1435 *chisq = 1.0e38f; // don't let chisq=0 through yfit being all ve 1436 } 1437 break; 1438 default: 1439 return 3; 1440 } 1441 break; 1442 case FIT_GLOBAL_STRETCHEDEXP: 1443 switch (noise) { 1444 case NOISE_CONST: 1445 // loop over all data 1446 for (i = fit_start; i < fit_end; ++i) { 1447 // stretched exponential fit 1448 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1449 dy[i] = y[i]  yfit[i]; 1450 1451 dy_dparam[i][0] = 1.0; 1452 dy_dparam[i][1] = exp_conv[1][i]; 1453 dy_dparam[i][2] = param[1] * exp_conv[2][i]; 1454 dy_dparam[i][3] = param[1] * exp_conv[3][i]; 1455 1456 // constant noise 1457 weight = 1.0f / sig[0]; 1458 alpha_weight[i] = weight; // 1 / (sig[0] * sig[0]); 1459 weight *= dy[i]; 1460 beta_weight[i] = weight; // dy[i] / (sig[0] * sig[0]); 1461 weight *= dy[i]; 1462 *chisq += weight; // (dy[i] * dy[i]) / (sig[0] * sig[0]); 1463 } 1464 break; 1465 case NOISE_GIVEN: 1466 // loop over all data 1467 for (i = fit_start; i < fit_end; ++i) { 1468 // stretched exponential fit 1469 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1470 dy[i] = y[i]  yfit[i]; 1471 1472 dy_dparam[i][0] = 1.0; 1473 dy_dparam[i][1] = exp_conv[1][i]; 1474 dy_dparam[i][2] = param[1] * exp_conv[2][i]; 1475 dy_dparam[i][3] = param[1] * exp_conv[3][i]; 1476 1477 // given noise 1478 weight = 1.0f / (sig[i] * sig[i]); 1479 alpha_weight[i] = weight; // 1 / (sig[i] * sig[i]) 1480 weight *= dy[i]; 1481 beta_weight[i] = weight; // dy[i] / (sig[i] * sig[i]) 1482 weight *= dy[i]; 1483 *chisq += weight; // (dy[i] * dy[i]) / (sig[i] * sig[i]) 1484 } 1485 break; 1486 case NOISE_POISSON_DATA: 1487 // loop over all data 1488 for (i = fit_start; i < fit_end; ++i) { 1489 // stretched exponential fit 1490 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1491 dy[i] = y[i]  yfit[i]; 1492 1493 dy_dparam[i][0] = 1.0; 1494 dy_dparam[i][1] = exp_conv[1][i]; 1495 dy_dparam[i][2] = param[1] * exp_conv[2][i]; 1496 dy_dparam[i][3] = param[1] * exp_conv[3][i]; 1497 1498 // poisson noise based on data 1499 weight = (y[i] > 15 ? 1.0f / y[i] : 1.0f / 15); 1500 alpha_weight[i] = weight; // 1 / sig(i) 1501 weight *= dy[i]; 1502 beta_weight[i] = weight; // dy[i] / sig(i) 1503 weight *= dy[i]; 1504 *chisq += weight; // (dy[i] * dy[i]) / sig(i) 1505 } 1506 break; 1507 case NOISE_POISSON_FIT: 1508 // loop over all data 1509 for (i = fit_start; i < fit_end; ++i) { 1510 // stretched exponential fit 1511 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1512 dy[i] = y[i]  yfit[i]; 1513 1514 dy_dparam[i][0] = 1.0; 1515 dy_dparam[i][1] = exp_conv[1][i]; 1516 dy_dparam[i][2] = param[1] * exp_conv[2][i]; 1517 dy_dparam[i][3] = param[1] * exp_conv[3][i]; 1518 1519 // poisson noise based on fit 1520 weight = (yfit[i] > 15 ? 1.0f / yfit[i] : 1.0f / 15); 1521 alpha_weight[i] = weight; // 1 / sig(i) 1522 weight *= dy[i]; 1523 beta_weight[i] = weight; // dy(i) / sig(i) 1524 weight *= dy[i]; 1525 *chisq += weight; // (dy(i) * dy(i)) / sig(i) 1526 } 1527 break; 1528 case NOISE_GAUSSIAN_FIT: 1529 // loop over all data 1530 for (i = fit_start; i < fit_end; ++i) { 1531 // stretched exponential fit 1532 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1533 dy[i] = y[i]  yfit[i]; 1534 1535 dy_dparam[i][0] = 1.0; 1536 dy_dparam[i][1] = exp_conv[1][i]; 1537 dy_dparam[i][2] = param[1] * exp_conv[2][i]; 1538 dy_dparam[i][3] = param[1] * exp_conv[3][i]; 1539 1540 // gaussian noise based on fit 1541 weight = (yfit[i] > 1.0f ? 1.0f / yfit[i] : 1.0f); 1542 alpha_weight[i] = weight; // 1 / sig(i) 1543 weight *= dy[i]; 1544 beta_weight[i] = weight; // dy[i] / sig(i) 1545 weight *= dy[i]; 1546 *chisq += weight; // dy[i] / sig(i) 1547 } 1548 break; 1549 case NOISE_MLE: 1550 // loop over all data 1551 for (i = fit_start; i < fit_end; ++i) { 1552 // stretched exponential fit 1553 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1554 dy[i] = y[i]  yfit[i]; 1555 1556 dy_dparam[i][0] = 1.0; 1557 dy_dparam[i][1] = exp_conv[1][i]; 1558 dy_dparam[i][2] = param[1] * exp_conv[2][i]; 1559 dy_dparam[i][3] = param[1] * exp_conv[3][i]; 1560 1561 // maximum likelihood estimation noise 1562 weight = (yfit[i] > 1 ? 1.0f / yfit[i] : 1.0f); 1563 alpha_weight[i] = weight * y[i] / yfit[i]; 1564 beta_weight[i] = dy[i] * weight; 1565 if (yfit[i] > 0.0) { 1566 *chisq += (0.0f == y[i]) 1567 ? 2.0 * yfit[i] 1568 : 2.0 * (yfit[i]  y[i])  2.0 * y[i] * log(yfit[i] / y[i]); 1569 } 1570 } 1571 if (*chisq <= 0.0f) { 1572 *chisq = 1.0e38f; // don't let chisq=0 through yfit being all ve 1573 } 1574 break; 1575 default: 1576 return 3; 1577 } 1578 break; 1579 default: 1580 dbgprintf(1, "compute_exps_fn: please update me!\n"); 1581 return 1; 1582 } 1583 1584 // Check if chi square worsened: 1585 if (0.0f != old_chisq && *chisq >= old_chisq) { 1586 // don't bother to set up the matrices for solution 1587 return 0; 1588 } 1589 1590 i_free = 0; 1591 // for all columns 1592 for (i = 0; i < nparam; ++i) { 1593 if (paramfree[i]) { 1594 j_free = 0; 1595 beta_sum = 0.0f; 1596 // row loop, only need to consider lower triangle 1597 for (j = 0; j <= i; ++j) { 1598 if (paramfree[j]) { 1599 dot_product = 0.0f; 1600 if (0 == j_free) { // true only once for each outer loop i 1601 // for all data 1602 for (k = fit_start; k < fit_end; ++k) { 1603 dy_dparam_k_i = dy_dparam[k][i]; 1604 dot_product += dy_dparam_k_i * dy_dparam[k][j] * alpha_weight[k]; //TODO ARG make it [i][k] and just *dy_dparam++ it. 1605 beta_sum += dy_dparam_k_i * beta_weight[k]; 1606 } 1607 } 1608 else { 1609 // for all data 1610 for (k = fit_start; k < fit_end; ++k) { 1611 dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k]; 1612 } 1613 } // k loop 1614 1615 alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 1616 // if (i_free != j_free) { 1617 // // matrix is symmetric 1618 // alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 1619 // } 1620 ++j_free; 1289 1621 } 1290 dy[i] = y[i]  yfit[i]; 1291 1292 for (j=0, l=0; l<nparam; l++) { 1293 if (paramfree[l]) { 1294 wt = dy_dparam[l]; 1295 for (k=0, m=0; m<=l; m++) 1296 if (paramfree[m]) 1297 alpha[j][k++] += wt * dy_dparam[m]; 1298 beta[j] += wt * dy[i]; 1299 j++; 1300 } 1301 } 1302 /* And find chi^2 */ 1303 *chisq += dy[i] * dy[i]; 1304 } 1305 1306 /* Now divide everything by sigma^2 */ 1307 sig2i = 1.0 / (sig[0] * sig[0]); 1308 *chisq *= sig2i; 1309 for (j=0; j<mfit; j++) { 1310 for (k=0; k<=j; k++) 1311 alpha[j][k] *= sig2i; 1312 beta[j] *= sig2i; 1313 } 1314 break; 1315 1316 case NOISE_GIVEN: /* This is essentially the NR version */ 1317 *chisq = 0.0; 1318 /* Summation loop over all data */ 1319 for (i=fit_start; i<fit_end; i++) { 1320 yfit[i] = param[0]; /* Z */ 1321 dy_dparam[0] = 1.0; 1322 1323 for (j=1; j<nparam; j+=2) { 1324 yfit[i] += param[j] * exp_conv[j+1][i]; 1325 /* A_j . exp(x/tau_j) */ 1326 dy_dparam[j] = exp_conv[j+1][i]; 1327 /* exp(x/tau_j) */ 1328 dy_dparam[j+1] = param[j] * exp_conv[j][i]; 1329 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1330 } 1331 sig2i = 1.0 / (sig[i] * sig[i]); 1332 dy[i] = y[i]  yfit[i]; 1333 1334 for (j=0, l=0; l<nparam; l++) { 1335 if (paramfree[l]) { 1336 wt = dy_dparam[l] * sig2i; 1337 for (k=0, m=0; m<=l; m++) 1338 if (paramfree[m]) 1339 alpha[j][k++] += wt * dy_dparam[m]; 1340 beta[j] += wt * dy[i]; 1341 j++; 1342 } 1343 } 1344 /* And find chi^2 */ 1345 *chisq += dy[i] * dy[i] * sig2i; 1346 } 1347 break; 1348 1349 case NOISE_POISSON_DATA: 1350 *chisq = 0.0; 1351 /* Summation loop over all data */ 1352 for (i=fit_start; i<fit_end; i++) { 1353 yfit[i] = param[0]; /* Z */ 1354 dy_dparam[0] = 1.0; 1355 1356 for (j=1; j<nparam; j+=2) { 1357 yfit[i] += param[j] * exp_conv[j+1][i]; 1358 /* A_j . exp(x/tau_j) */ 1359 dy_dparam[j] = exp_conv[j+1][i]; 1360 /* exp(x/tau_j) */ 1361 dy_dparam[j+1] = param[j] * exp_conv[j][i]; 1362 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1363 } 1364 sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15); 1365 dy[i] = y[i]  yfit[i]; 1366 1367 for (j=0, l=0; l<nparam; l++) { 1368 if (paramfree[l]) { 1369 wt = dy_dparam[l] * sig2i; 1370 for (k=0, m=0; m<=l; m++) 1371 if (paramfree[m]) 1372 alpha[j][k++] += wt * dy_dparam[m]; 1373 beta[j] += wt * dy[i]; 1374 j++; 1375 } 1376 } 1377 /* And find chi^2 */ 1378 *chisq += dy[i] * dy[i] * sig2i; 1379 } 1380 break; 1381 1382 case NOISE_POISSON_FIT: 1383 *chisq = 0.0; 1384 /* Summation loop over all data */ 1385 for (i=fit_start; i<fit_end; i++) { 1386 yfit[i] = param[0]; /* Z */ 1387 dy_dparam[0] = 1.0; 1388 1389 for (j=1; j<nparam; j+=2) { 1390 yfit[i] += param[j] * exp_conv[j+1][i]; 1391 /* A_j . exp(x/tau_j) */ 1392 dy_dparam[j] = exp_conv[j+1][i]; 1393 /* exp(x/tau_j) */ 1394 dy_dparam[j+1] = param[j] * exp_conv[j][i]; 1395 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1396 } 1397 sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15); 1398 dy[i] = y[i]  yfit[i]; 1399 1400 for (j=0, l=0; l<nparam; l++) { 1401 if (paramfree[l]) { 1402 wt = dy_dparam[l] * sig2i; 1403 for (k=0, m=0; m<=l; m++) 1404 if (paramfree[m]) 1405 alpha[j][k++] += wt * dy_dparam[m]; 1406 beta[j] += wt * dy[i]; 1407 j++; 1408 } 1409 } 1410 /* And find chi^2 */ 1411 *chisq += dy[i] * dy[i] * sig2i; 1412 } 1413 break; 1414 1415 case NOISE_GAUSSIAN_FIT: 1416 *chisq = 0.0; 1417 /* Summation loop over all data */ 1418 for (i=fit_start; i<fit_end; i++) { 1419 yfit[i] = param[0]; /* Z */ 1420 dy_dparam[0] = 1.0; 1421 1422 for (j=1; j<nparam; j+=2) { 1423 yfit[i] += param[j] * exp_conv[j+1][i]; 1424 /* A_j . exp(x/tau_j) */ 1425 dy_dparam[j] = exp_conv[j+1][i]; 1426 /* exp(x/tau_j) */ 1427 dy_dparam[j+1] = param[j] * exp_conv[j][i]; 1428 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1429 } 1430 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1431 dy[i] = y[i]  yfit[i]; 1432 1433 for (j=0, l=0; l<nparam; l++) { 1434 if (paramfree[l]) { 1435 wt = dy_dparam[l] * sig2i; 1436 for (k=0, m=0; m<=l; m++) 1437 if (paramfree[m]) 1438 alpha[j][k++] += wt * dy_dparam[m]; 1439 beta[j] += wt * dy[i]; 1440 j++; 1441 } 1442 } 1443 /* And find chi^2 */ 1444 *chisq += dy[i] * dy[i] * sig2i; 1445 } 1446 break; 1447 1448 case NOISE_MLE: 1449 *chisq = 0.0; 1450 /* Summation loop over all data */ 1451 for (i=fit_start; i<fit_end; i++) { 1452 yfit[i] = param[0]; /* Z */ 1453 dy_dparam[0] = 1.0; 1454 1455 for (j=1; j<nparam; j+=2) { 1456 yfit[i] += param[j] * exp_conv[j+1][i]; 1457 /* A_j . exp(x/tau_j) */ 1458 dy_dparam[j] = exp_conv[j+1][i]; 1459 /* exp(x/tau_j) */ 1460 dy_dparam[j+1] = param[j] * exp_conv[j][i]; 1461 /* A_j * x * exp(x/tau_j) / tau_j^2 */ 1462 } 1463 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1464 dy[i] = y[i]  yfit[i]; 1465 y_ymod=y[i]/yfit[i]; 1466 1467 for (j=0, l=0; l<nparam; l++) { 1468 if (paramfree[l]) { 1469 wt = dy_dparam[l] * sig2i; 1470 for (k=0, m=0; m<=l; m++) 1471 if (paramfree[m]) 1472 alpha[j][k++] += wt*dy_dparam[m]*y_ymod; // wt * dy_dparam[m]; 1473 beta[j] += wt * dy[i]; 1474 j++; 1475 } 1476 } 1477 // And find chi^2 1478 if (yfit[i]<=0.0) 1479 ; // do nothing 1480 else if (y[i]==0.0) 1481 *chisq += 2.0*yfit[i]; // to avoid NaN from log 1482 else 1483 *chisq += 2.0*(yfit[i]y[i])  2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; 1484 } 1485 if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all ve 1486 break; 1487 1488 default: 1489 return 3; 1490 /* break; */ // (unreachable code) 1491 } 1492 break; 1493 1494 case FIT_GLOBAL_STRETCHEDEXP: 1495 switch (noise) { 1496 case NOISE_CONST: 1497 *chisq = 0.0; 1498 /* Summation loop over all data */ 1499 for (i=fit_start; i<fit_end; i++) { 1500 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1501 dy[i] = y[i]  yfit[i]; 1502 1503 dy_dparam[0] = 1.0; 1504 dy_dparam[1] = exp_conv[1][i]; 1505 dy_dparam[2] = param[1] * exp_conv[2][i]; 1506 dy_dparam[3] = param[1] * exp_conv[3][i]; 1507 1508 for (j=0, l=0; l<nparam; l++) { 1509 if (paramfree[l]) { 1510 wt = dy_dparam[l]; 1511 for (k=0, m=0; m<=l; m++) 1512 if (paramfree[m]) 1513 alpha[j][k++] += wt * dy_dparam[m]; 1514 beta[j] += wt * dy[i]; 1515 j++; 1516 } 1517 } 1518 /* And find chi^2 */ 1519 *chisq += dy[i] * dy[i]; 1520 } 1521 1522 /* Now divide everything by sigma^2 */ 1523 sig2i = 1.0 / (sig[0] * sig[0]); 1524 *chisq *= sig2i; 1525 for (j=0; j<mfit; j++) { 1526 for (k=0; k<=j; k++) 1527 alpha[j][k] *= sig2i; 1528 beta[j] *= sig2i; 1529 } 1530 break; 1531 1532 case NOISE_GIVEN: /* This is essentially the NR version */ 1533 *chisq = 0.0; 1534 /* Summation loop over all data */ 1535 for (i=fit_start; i<fit_end; i++) { 1536 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1537 dy[i] = y[i]  yfit[i]; 1538 1539 dy_dparam[0] = 1.0; 1540 dy_dparam[1] = exp_conv[1][i]; 1541 dy_dparam[2] = param[1] * exp_conv[2][i]; 1542 dy_dparam[3] = param[1] * exp_conv[3][i]; 1543 1544 sig2i = 1.0 / (sig[i] * sig[i]); 1545 1546 for (j=0, l=0; l<nparam; l++) { 1547 if (paramfree[l]) { 1548 wt = dy_dparam[l] * sig2i; 1549 for (k=0, m=0; m<=l; m++) 1550 if (paramfree[m]) 1551 alpha[j][k++] += wt * dy_dparam[m]; 1552 beta[j] += wt * dy[i]; 1553 j++; 1554 } 1555 } 1556 /* And find chi^2 */ 1557 *chisq += dy[i] * dy[i] * sig2i; 1558 } 1559 break; 1560 1561 case NOISE_POISSON_DATA: 1562 *chisq = 0.0; 1563 /* Summation loop over all data */ 1564 for (i=fit_start; i<fit_end; i++) { 1565 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1566 dy[i] = y[i]  yfit[i]; 1567 1568 dy_dparam[0] = 1.0; 1569 dy_dparam[1] = exp_conv[1][i]; 1570 dy_dparam[2] = param[1] * exp_conv[2][i]; 1571 dy_dparam[3] = param[1] * exp_conv[3][i]; 1572 1573 sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15); 1574 1575 for (j=0, l=0; l<nparam; l++) { 1576 if (paramfree[l]) { 1577 wt = dy_dparam[l] * sig2i; 1578 for (k=0, m=0; m<=l; m++) 1579 if (paramfree[m]) 1580 alpha[j][k++] += wt * dy_dparam[m]; 1581 beta[j] += wt * dy[i]; 1582 j++; 1583 } 1584 } 1585 /* And find chi^2 */ 1586 *chisq += dy[i] * dy[i] * sig2i; 1587 } 1588 break; 1589 1590 case NOISE_POISSON_FIT: 1591 *chisq = 0.0; 1592 /* Summation loop over all data */ 1593 for (i=fit_start; i<fit_end; i++) { 1594 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1595 dy[i] = y[i]  yfit[i]; 1596 1597 dy_dparam[0] = 1.0; 1598 dy_dparam[1] = exp_conv[1][i]; 1599 dy_dparam[2] = param[1] * exp_conv[2][i]; 1600 dy_dparam[3] = param[1] * exp_conv[3][i]; 1601 1602 sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15); 1603 1604 for (j=0, l=0; l<nparam; l++) { 1605 if (paramfree[l]) { 1606 wt = dy_dparam[l] * sig2i; 1607 for (k=0, m=0; m<=l; m++) 1608 if (paramfree[m]) 1609 alpha[j][k++] += wt * dy_dparam[m]; 1610 beta[j] += wt * dy[i]; 1611 j++; 1612 } 1613 } 1614 /* And find chi^2 */ 1615 *chisq += dy[i] * dy[i] * sig2i; 1616 } 1617 break; 1618 1619 case NOISE_GAUSSIAN_FIT: 1620 *chisq = 0.0; 1621 /* Summation loop over all data */ 1622 for (i=fit_start; i<fit_end; i++) { 1623 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1624 dy[i] = y[i]  yfit[i]; 1625 1626 dy_dparam[0] = 1.0; 1627 dy_dparam[1] = exp_conv[1][i]; 1628 dy_dparam[2] = param[1] * exp_conv[2][i]; 1629 dy_dparam[3] = param[1] * exp_conv[3][i]; 1630 1631 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1632 1633 for (j=0, l=0; l<nparam; l++) { 1634 if (paramfree[l]) { 1635 wt = dy_dparam[l] * sig2i; 1636 for (k=0, m=0; m<=l; m++) 1637 if (paramfree[m]) 1638 alpha[j][k++] += wt * dy_dparam[m]; 1639 beta[j] += wt * dy[i]; 1640 j++; 1641 } 1642 } 1643 /* And find chi^2 */ 1644 *chisq += dy[i] * dy[i] * sig2i; 1645 } 1646 break; 1647 1648 case NOISE_MLE: 1649 *chisq = 0.0; 1650 /* Summation loop over all data */ 1651 for (i=fit_start; i<fit_end; i++) { 1652 yfit[i] = param[0] + param[1] * exp_conv[1][i]; 1653 dy[i] = y[i]  yfit[i]; 1654 1655 dy_dparam[0] = 1.0; 1656 dy_dparam[1] = exp_conv[1][i]; 1657 dy_dparam[2] = param[1] * exp_conv[2][i]; 1658 dy_dparam[3] = param[1] * exp_conv[3][i]; 1659 1660 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1661 y_ymod=y[i]/yfit[i]; 1662 1663 for (j=0, l=0; l<nparam; l++) { 1664 if (paramfree[l]) { 1665 wt = dy_dparam[l] * sig2i; 1666 for (k=0, m=0; m<=l; m++) 1667 if (paramfree[m]) 1668 alpha[j][k++] += wt*dy_dparam[m]*y_ymod; //wt * dy_dparam[m]; 1669 beta[j] += wt * dy[i]; 1670 j++; 1671 } 1672 } 1673 // And find chi^2 1674 if (yfit[i]<=0.0) 1675 ; // do nothing 1676 else if (y[i]==0.0) 1677 *chisq += 2.0*yfit[i]; // to avoid NaN from log 1678 else 1679 *chisq += 2.0*(yfit[i]y[i])  2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; 1680 } 1681 if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all ve 1682 break; 1683 1684 default: 1685 return 3; 1686 /* break; */ // (unreachable code) 1687 } 1688 break; 1689 1690 default: 1691 dbgprintf(1, "compute_exps_fn: please update me!\n"); 1692 return 1; 1693 } 1694 1695 /* Fill in the symmetric side */ 1696 for (j=1; j<mfit; j++) 1697 for (k=0; k<j; k++) 1698 alpha[k][j] = alpha[j][k]; 1622 } // j loop 1623 beta[i_free] = beta_sum; 1624 ++i_free; 1625 } 1626 } // i loop 1699 1627 1700 1628 return 0; … … 2463 2391 //// exp_conv, yfit[i], dy[i], 2464 2392 exp_conv, yfit[0], dy[0], 2465 alpha_scratch, beta_scratch, &chisq_trans[i] );2393 alpha_scratch, beta_scratch, &chisq_trans[i], 0.0); 2466 2394 2467 2395 if (ret != 0) {
Note: See TracChangeset
for help on using the changeset viewer.