Ignore:
Timestamp:
06/01/11 20:17:05 (9 years ago)
Author:
aivar
Message:

Updated GCI_marquardt_global_compute_exps_fn.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/projects/slim-curve/src/main/c/EcfGlobal.c

    r7706 r7720  
    7979                        float *exp_conv[], 
    8080                        float yfit[], float dy[], 
    81                         float **alpha, float *beta, float *chisq); 
     81                        float **alpha, float *beta, float *chisq, float old_chisq); 
    8282int GCI_marquardt_global_compute_exps_fn_final( 
    8383                        float xincr, float y[], 
     
    11641164                                        xincr, y, ndata, fit_start, fit_end, noise, sig, 
    11651165                                        ftype, param, paramfree, nparam, exp_conv, 
    1166                                         yfit, dy, alpha, beta, chisq) != 0) 
     1166                                        yfit, dy, alpha, beta, chisq, 0.0) != 0) 
    11671167                        return -2; 
    11681168 
     
    12191219                                xincr, y, ndata, fit_start, fit_end, noise, sig, 
    12201220                                ftype, paramtry, paramfree, nparam, exp_conv, 
    1221                                 yfit, dy, covar, dparam, chisq) != 0) 
     1221                                yfit, dy, covar, dparam, chisq, ochisq) != 0) 
    12221222                return -2; 
    12231223 
     
    12501250                        float *exp_conv[], 
    12511251                        float yfit[], float dy[], 
    1252                         float **alpha, float *beta, float *chisq) 
     1252                        float **alpha, float *beta, 
     1253                        float *chisq, float old_chisq) 
    12531254{ 
    12541255        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; 
    12561266 
    12571267        for (j=0, mfit=0; j<nparam; j++) 
     
    12591269                        mfit++; 
    12601270 
    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; 
    12711272 
    12721273        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                                                // multi-exponential 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                                                // multi-exponential 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                                                // multi-exponential 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                                                // multi-exponential 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                                                // multi-exponential 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                                                // multi-exponential 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; 
    12891621                                } 
    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 
    16991627 
    17001628        return 0; 
     
    24632391////                                    exp_conv, yfit[i], dy[i], 
    24642392                                        exp_conv, yfit[0], dy[0], 
    2465                                         alpha_scratch, beta_scratch, &chisq_trans[i]); 
     2393                                        alpha_scratch, beta_scratch, &chisq_trans[i], 0.0); 
    24662394 
    24672395                if (ret != 0) { 
Note: See TracChangeset for help on using the changeset viewer.