Changeset 7719 for trunk/projects


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

Fixed tabs. Moved newly added variables to start of method.

File:
1 edited

Legend:

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

    r7715 r7719  
    984984        int i, j, k, l, m, mfit; 
    985985        float wt, sig2i, y_ymod, dy_dparam[MAXBINS][MAXFIT]; 
     986        float alpha_weight[MAXBINS]; 
     987        float beta_weight[MAXBINS]; 
     988        int q; 
     989        float weight; 
     990        int i_free; 
     991        int j_free; 
     992        float dot_product; 
     993        float beta_sum; 
     994        float dy_dparam_k_i; 
    986995 
    987996        for (j=0, mfit=0; j<nparam; j++) 
     
    989998                        mfit++; 
    990999 
    991         //TODO ARG having this section { } of code here is just temporary; 
    992         // o'wise have to define these vars at top of function 
    993         { 
    994         float alpha_weight[MAXBINS]; 
    995         float beta_weight[MAXBINS]; 
    996         int q; 
    997         float weight; 
    998  
    999         int i_free; 
    1000         int j_free; 
    1001         float dot_product; 
    1002         float beta_sum; 
    1003         float dy_dparam_k_i; 
    1004  
    1005         *chisq = 0.0f; 
    1006  
    1007         switch (noise) { 
    1008             case NOISE_CONST: 
    1009             { 
    1010                 float i_sig_0_squared = 1.0 / (sig[0] * sig[0]); 
    1011                 for (q = 0; q < ndata; ++q) { 
    1012                     (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
    1013                     yfit[q] += param[0]; 
    1014                     dy_dparam[q][0] = 1.0; 
    1015                     dy[q] = y[q] - yfit[q]; 
    1016                     weight = i_sig_0_squared; 
    1017                     alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]) 
    1018                     weight *= dy[q]; 
    1019                     beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]) 
    1020                     weight *= dy[q]; 
    1021                     *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]) 
    1022                 } 
    1023                 break; 
    1024             } 
    1025             case NOISE_GIVEN: 
    1026             { 
    1027                 for (q = 0; q < ndata; ++q) { 
    1028                     (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
    1029                     yfit[q] += param[0]; 
    1030                     dy_dparam[q][0] = 1.0; 
    1031                     dy[q] = y[q] - yfit[q]; 
    1032                     weight = 1.0f / (sig[q] * sig[q]); 
    1033                     alpha_weight[q] = weight; // 1 / (sig[q] * sig[q]) 
    1034                     weight *= dy[q]; 
    1035                     beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q]) 
    1036                     weight *= dy[q]; 
    1037                     *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q]) 
    1038                 } 
    1039                 break; 
    1040             } 
    1041             case NOISE_POISSON_DATA: 
    1042             { 
    1043                 for (q = 0; q < ndata; ++q) { 
    1044                     (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
    1045                     yfit[q] += param[0]; 
    1046                     dy_dparam[q][0] = 1.0; 
    1047                     dy[q] = y[q] - yfit[q]; 
    1048                     weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 
    1049                     alpha_weight[q] = weight; // 1 / sig(q) 
    1050                     weight *= dy[q]; 
    1051                     beta_weight[q] = weight; // dy[q] / sig(q) 
    1052                     weight *= dy[q]; 
    1053                     *chisq += weight; // (dy[q] * dy[q]) / sig(q) 
    1054                 } 
    1055                 break; 
    1056             } 
    1057             case NOISE_POISSON_FIT: 
    1058             { 
    1059                 for (q = 0; q < ndata; ++q) { 
    1060                     (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
    1061                     yfit[q] += param[0]; 
    1062                     dy_dparam[q][0] = 1.0; 
    1063                     dy[q] = y[q] - yfit[q]; 
    1064                     weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 
    1065                     alpha_weight[q] = weight; // 1 / sig(q) 
    1066                     weight *= dy[q]; 
    1067                     beta_weight[q] = weight; // dy(q) / sig(q) 
    1068                     weight *= dy[q]; 
    1069                     *chisq += weight; // (dy(q) * dy(q)) / sig(q) 
    1070                 } 
    1071                 break; 
    1072             } 
    1073             case NOISE_GAUSSIAN_FIT: 
    1074             { 
    1075                  for (q = 0; q < ndata; ++q) { 
    1076                     (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
    1077                     yfit[q] += param[0]; 
    1078                     dy_dparam[q][0] = 1.0; 
    1079                     dy[q] = y[q] - yfit[q]; 
    1080                     weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 
    1081                     alpha_weight[q] = weight; // 1 / sig(q) 
    1082                     weight *= dy[q]; 
    1083                     beta_weight[q] = weight; // dy[q] / sig(q) 
    1084                     weight *= dy[q]; 
    1085                     *chisq += weight; // dy[q] / sig(q) 
    1086                  } 
    1087                  break; 
    1088            } 
    1089             case NOISE_MLE: 
    1090             { 
    1091                 for (q = 0; q < ndata; ++q) { 
    1092                     (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
    1093                     yfit[q] += param[0]; 
    1094                     dy_dparam[q][0] = 1.0; 
    1095                     dy[q] = y[q] - yfit[q]; 
    1096                     weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f); 
    1097                     alpha_weight[q] = weight * y[q] / yfit[q]; 
    1098                     beta_weight[q] = dy[q] * weight; 
    1099                     if (yfit[q] > 0.0) { 
    1100                         *chisq += (0.0f == y[q]) 
    1101                                 ? 2.0 * yfit[q] 
    1102                                 : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]); 
    1103                     } 
    1104                 } 
    1105                 if (*chisq <= 0.0f) { 
    1106                     *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
    1107                 } 
    1108                 break; 
    1109             } 
    1110             default: 
    1111             { 
    1112                 return -3; 
    1113             } 
    1114         } 
    1115  
    1116         // Check if chi square worsened: 
    1117         if (0.0f != old_chisq && *chisq >= old_chisq) { 
    1118             // don't bother to set up the matrices for solution 
    1119             return 0; 
    1120         } 
    1121  
    1122         i_free = 0; 
    1123         // for all columns 
    1124         for (i = 0; i < nparam; ++i) { 
    1125             if (paramfree[i]) { 
    1126                 j_free = 0; 
    1127                 beta_sum = 0.0f; 
    1128                 // row loop, only need to consider lower triangle 
    1129                 for (j = 0; j <= i; ++j) { 
    1130                     if (paramfree[j]) { 
    1131                         dot_product = 0.0f; 
    1132                         if (0 == j_free) { // true only once for each outer loop i 
    1133                             // for all data 
    1134                             for (k = 0; k < ndata; ++k) { 
    1135                                                                 dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k]; 
    1136                                                                 beta_sum += dy_dparam[k][i] * beta_weight[k]; // compute beta only once for each i 
    1137                             } 
    1138                         } 
    1139                         else { 
    1140                             // for all data 
    1141                             for (k = 0; k < ndata; ++k) { 
    1142                                                                 dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k]; 
    1143                             } 
    1144                         } // k loop 
    1145  
    1146                         alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; //TODO ARG w/n/b [i][j] more common usage?  row/column 
    1147                        // if (i_free != j_free) { //TODO ARG this approach seemed slower at one time; c/b worth retesting: 
    1148                        //     / is symmetric 
    1149                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 
    1150                        // } 
    1151                         ++j_free; 
    1152                     } 
    1153                 } // j loop 
    1154                 beta[i_free] = beta_sum; 
    1155                 ++i_free; 
    1156             } 
    1157         } // i loop 
    1158  
    1159         } 
    1160  
     1000        *chisq = 0.0f; 
     1001 
     1002        switch (noise) { 
     1003                case NOISE_CONST: 
     1004                { 
     1005                        float i_sig_0_squared = 1.0 / (sig[0] * sig[0]); 
     1006                        for (q = 0; q < ndata; ++q) { 
     1007                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
     1008                                yfit[q] += param[0]; 
     1009                                dy_dparam[q][0] = 1.0; 
     1010                                dy[q] = y[q] - yfit[q]; 
     1011                                weight = i_sig_0_squared; 
     1012                                alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]) 
     1013                                weight *= dy[q]; 
     1014                                beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]) 
     1015                                weight *= dy[q]; 
     1016                                *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]) 
     1017                        } 
     1018                        break; 
     1019                } 
     1020                case NOISE_GIVEN: 
     1021                { 
     1022                        for (q = 0; q < ndata; ++q) { 
     1023                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
     1024                                yfit[q] += param[0]; 
     1025                                dy_dparam[q][0] = 1.0; 
     1026                                dy[q] = y[q] - yfit[q]; 
     1027                                weight = 1.0f / (sig[q] * sig[q]); 
     1028                                alpha_weight[q] = weight; // 1 / (sig[q] * sig[q]) 
     1029                                weight *= dy[q]; 
     1030                                beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q]) 
     1031                                weight *= dy[q]; 
     1032                                *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q]) 
     1033                        } 
     1034                        break; 
     1035                } 
     1036                case NOISE_POISSON_DATA: 
     1037                { 
     1038                        for (q = 0; q < ndata; ++q) { 
     1039                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
     1040                                yfit[q] += param[0]; 
     1041                                dy_dparam[q][0] = 1.0; 
     1042                                dy[q] = y[q] - yfit[q]; 
     1043                                weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 
     1044                                alpha_weight[q] = weight; // 1 / sig(q) 
     1045                                weight *= dy[q]; 
     1046                                beta_weight[q] = weight; // dy[q] / sig(q) 
     1047                                weight *= dy[q]; 
     1048                                *chisq += weight; // (dy[q] * dy[q]) / sig(q) 
     1049                        } 
     1050                        break; 
     1051                } 
     1052                case NOISE_POISSON_FIT: 
     1053                { 
     1054                        for (q = 0; q < ndata; ++q) { 
     1055                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
     1056                                yfit[q] += param[0]; 
     1057                                dy_dparam[q][0] = 1.0; 
     1058                                dy[q] = y[q] - yfit[q]; 
     1059                                weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 
     1060                                alpha_weight[q] = weight; // 1 / sig(q) 
     1061                                weight *= dy[q]; 
     1062                                beta_weight[q] = weight; // dy(q) / sig(q) 
     1063                                weight *= dy[q]; 
     1064                                *chisq += weight; // (dy(q) * dy(q)) / sig(q) 
     1065                        } 
     1066                        break; 
     1067                } 
     1068                case NOISE_GAUSSIAN_FIT: 
     1069                { 
     1070                        for (q = 0; q < ndata; ++q) { 
     1071                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
     1072                                yfit[q] += param[0]; 
     1073                                dy_dparam[q][0] = 1.0; 
     1074                                dy[q] = y[q] - yfit[q]; 
     1075                                weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 
     1076                                alpha_weight[q] = weight; // 1 / sig(q) 
     1077                                weight *= dy[q]; 
     1078                                beta_weight[q] = weight; // dy[q] / sig(q) 
     1079                                weight *= dy[q]; 
     1080                                *chisq += weight; // dy[q] / sig(q) 
     1081                        } 
     1082                        break; 
     1083                } 
     1084                case NOISE_MLE: 
     1085                { 
     1086                        for (q = 0; q < ndata; ++q) { 
     1087                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 
     1088                                yfit[q] += param[0]; 
     1089                                dy_dparam[q][0] = 1.0; 
     1090                                dy[q] = y[q] - yfit[q]; 
     1091                                weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f); 
     1092                                alpha_weight[q] = weight * y[q] / yfit[q]; 
     1093                                beta_weight[q] = dy[q] * weight; 
     1094                                if (yfit[q] > 0.0) { 
     1095                                        *chisq += (0.0f == y[q]) 
     1096                                                        ? 2.0 * yfit[q] 
     1097                                                        : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]); 
     1098                                } 
     1099                        } 
     1100                        if (*chisq <= 0.0f) { 
     1101                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
     1102                        } 
     1103                        break; 
     1104                } 
     1105                default: 
     1106                { 
     1107                        return -3; 
     1108                } 
     1109        } 
     1110 
     1111        // Check if chi square worsened: 
     1112        if (0.0f != old_chisq && *chisq >= old_chisq) { 
     1113                // don't bother to set up the matrices for solution 
     1114                return 0; 
     1115        } 
     1116 
     1117        i_free = 0; 
     1118        // for all columns 
     1119        for (i = 0; i < nparam; ++i) { 
     1120                if (paramfree[i]) { 
     1121                        j_free = 0; 
     1122                        beta_sum = 0.0f; 
     1123                        // row loop, only need to consider lower triangle 
     1124                        for (j = 0; j <= i; ++j) { 
     1125                                if (paramfree[j]) { 
     1126                                        dot_product = 0.0f; 
     1127                                        if (0 == j_free) { // true only once for each outer loop i 
     1128                                                // for all data 
     1129                                                for (k = 0; k < ndata; ++k) { 
     1130                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k]; 
     1131                                                        beta_sum += dy_dparam[k][i] * beta_weight[k]; // compute beta only once for each i 
     1132                                                } 
     1133                                        } 
     1134                                        else { 
     1135                                                // for all data 
     1136                                                for (k = 0; k < ndata; ++k) { 
     1137                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k]; 
     1138                                                } 
     1139                                        } // k loop 
     1140                                                 
     1141                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; //TODO ARG w/n/b [i][j] more common usage?  row/column 
     1142                                        // if (i_free != j_free) { //TODO ARG this approach seemed slower at one time; c/b worth retesting: 
     1143                                        //     / is symmetric 
     1144                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 
     1145                                        // } 
     1146                                        ++j_free; 
     1147                                } 
     1148                        } // j loop 
     1149                        beta[i_free] = beta_sum; 
     1150                        ++i_free; 
     1151                } 
     1152        } // i loop 
     1153         
    11611154        return 0; 
    11621155} 
     
    11801173        int i, j, k, l, m, mfit, ret; 
    11811174        float wt, sig2i, y_ymod; 
    1182  
     1175        float alpha_weight[MAXBINS]; 
     1176        float beta_weight[MAXBINS]; 
     1177        int q; 
     1178        float weight; 
     1179        int i_free; 
     1180        int j_free; 
     1181        float dot_product; 
     1182        float beta_sum; 
     1183        float dy_dparam_k_i; 
     1184         
    11831185        /* Are we initialising? */ 
    11841186        // Malloc the arrays that will get used again in this fit via the pointers passed in 
     
    13041306        /* OK, now we've got our (possibly convolved) data, we can do the 
    13051307           rest almost exactly as above. */ 
    1306         //TODO ARG this new section of code is just temporary; o'wise have to define all these variables at the top of the function 
    1307         { 
    1308         float alpha_weight[MAXBINS]; 
    1309         float beta_weight[MAXBINS]; 
    1310         int q; 
    1311         float weight; 
    1312  
    1313         int i_free; 
    1314         int j_free; 
    1315         float dot_product; 
    1316         float beta_sum; 
    1317         float dy_dparam_k_i; 
    1318  
    1319         *chisq = 0.0f; 
    1320  
     1308         
     1309        *chisq = 0.0f; 
     1310         
    13211311        switch (noise) { 
    1322             case NOISE_CONST: 
    1323             { 
    1324                 for (q = fit_start; q < fit_end; ++q) { 
    1325                     (*pdy_dparam_conv)[q][0] = 1.0f; 
    1326                     yfit[q] += param[0]; 
    1327                     dy[q] = y[q] - yfit[q]; 
    1328                     weight = 1.0f / sig[0]; 
    1329                     alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]); 
    1330                     weight *= dy[q]; 
    1331                     beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]); 
    1332                     weight *= dy[q]; 
    1333                     *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]); 
    1334                 } 
    1335                 break; 
    1336             } 
    1337             case NOISE_GIVEN: 
    1338             { 
    1339                 for (q = fit_start; q < fit_end; ++q) { 
    1340                     (*pdy_dparam_conv)[q][0] = 1.0f; 
    1341                     yfit[q] += param[0]; 
    1342                     dy[q] = y[q] - yfit[q]; 
    1343                     weight = 1.0f / (sig[q] * sig[q]); 
    1344                     alpha_weight[q] = weight; // 1 / (sig[q] * sig[q]) 
    1345                     weight *= dy[q]; 
    1346                     beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q]) 
    1347                     weight *= dy[q]; 
    1348                     *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q]) 
    1349                 } 
    1350                 break; 
    1351             } 
    1352             case NOISE_POISSON_DATA: 
    1353             { 
    1354                 for (q = fit_start; q < fit_end; ++q) { 
    1355                     (*pdy_dparam_conv)[q][0] = 1.0f; 
    1356                     yfit[q] += param[0]; 
    1357                     dy[q] = y[q] - yfit[q]; 
    1358                     weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 
    1359                     alpha_weight[q] = weight; // 1 / sig(q) 
    1360                     weight *= dy[q]; 
    1361                     beta_weight[q] = weight; // dy[q] / sig(q) 
    1362                     weight *= dy[q]; 
    1363                     *chisq += weight; // (dy[q] * dy[q]) / sig(q) 
    1364                 } 
    1365                 break; 
    1366             } 
    1367             case NOISE_POISSON_FIT: 
    1368             { 
    1369                 for (q = fit_start; q < fit_end; ++q) { 
    1370                     (*pdy_dparam_conv)[q][0] = 1.0f; 
    1371                     yfit[q] += param[0]; 
    1372                     dy[q] = y[q] - yfit[q]; 
    1373                     weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 
    1374                     alpha_weight[q] = weight; // 1 / sig(q) 
    1375                     weight *= dy[q]; 
    1376                     beta_weight[q] = weight; // dy(q) / sig(q) 
    1377                     weight *= dy[q]; 
    1378                     *chisq += weight; // (dy(q) * dy(q)) / sig(q) 
    1379                 } 
    1380                 break; 
    1381             } 
    1382             case NOISE_GAUSSIAN_FIT: 
    1383             { 
    1384                  for (q = fit_start; q < fit_end; ++q) { 
    1385                     (*pdy_dparam_conv)[q][0] = 1.0f; 
    1386                     yfit[q] += param[0]; 
    1387                     dy[q] = y[q] - yfit[q]; 
    1388                     weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 
    1389                     alpha_weight[q] = weight; // 1 / sig(q) 
    1390                     weight *= dy[q]; 
    1391                     beta_weight[q] = weight; // dy[q] / sig(q) 
    1392                     weight *= dy[q]; 
    1393                     *chisq += weight; // dy[q] / sig(q) 
    1394                  } 
    1395                  break; 
    1396            } 
    1397             case NOISE_MLE: 
    1398             { 
    1399                 for (q = fit_start; q < fit_end; ++q) { 
    1400                     (*pdy_dparam_conv)[q][0] = 1.0f; 
    1401                     yfit[q] += param[0]; 
    1402                     dy[q] = y[q] - yfit[q]; 
    1403                     weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f); 
    1404                     alpha_weight[q] = weight * y[q] / yfit[q]; 
    1405                     beta_weight[q] = dy[q] * weight; 
    1406                     if (yfit[q] > 0.0) { 
    1407                         *chisq += (0.0f == y[q]) 
    1408                                 ? 2.0 * yfit[q] 
    1409                                 : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]); 
    1410                     } 
    1411                 } 
    1412                 if (*chisq <= 0.0f) { 
    1413                     *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
    1414                 } 
    1415                 break; 
    1416             } 
    1417             default: 
    1418             { 
    1419                 return -3; 
    1420             } 
    1421         } 
    1422  
    1423         // Check if chi square worsened: 
    1424         if (0.0f != old_chisq && *chisq >= old_chisq) { 
    1425             // don't bother to set up the matrices for solution 
    1426             return 0; 
    1427         } 
    1428  
    1429         i_free = 0; 
    1430         // for all columns 
    1431         for (i = 0; i < nparam; ++i) { 
    1432             if (paramfree[i]) { 
    1433                 j_free = 0; 
    1434                 beta_sum = 0.0f; 
    1435                 // row loop, only need to consider lower triangle 
    1436                 for (j = 0; j <= i; ++j) { 
    1437                     if (paramfree[j]) { 
    1438                         dot_product = 0.0f; 
    1439                         if (0 == j_free) { // true only once for each outer loop i 
    1440                             // for all data 
    1441                             for (k = fit_start; k < fit_end; ++k) { 
    1442                                 dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 
    1443                                 dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO ARG make it [i][k] and just *dy_dparam++ it. 
    1444                                 beta_sum += dy_dparam_k_i * beta_weight[k]; 
    1445                             } 
    1446                         } 
    1447                         else { 
    1448                             // for all data 
    1449                             for (k = fit_start; k < fit_end; ++k) { 
    1450                                 dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; 
    1451                             } 
    1452                         } // k loop 
    1453  
    1454                         alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 
    1455                        // if (i_free != j_free) { 
    1456                        //     // matrix is symmetric 
    1457                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 
    1458                        // } 
    1459                         ++j_free; 
    1460                     } 
    1461                 } // j loop 
    1462                 beta[i_free] = beta_sum; 
    1463                 ++i_free; 
    1464             } 
    1465         } // i loop 
    1466         } 
     1312                case NOISE_CONST: 
     1313                { 
     1314                        for (q = fit_start; q < fit_end; ++q) { 
     1315                                (*pdy_dparam_conv)[q][0] = 1.0f; 
     1316                                yfit[q] += param[0]; 
     1317                                dy[q] = y[q] - yfit[q]; 
     1318                                weight = 1.0f / sig[0]; 
     1319                                alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]); 
     1320                                weight *= dy[q]; 
     1321                                beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]); 
     1322                                weight *= dy[q]; 
     1323                                *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]); 
     1324                        } 
     1325                        break; 
     1326                } 
     1327                case NOISE_GIVEN: 
     1328                { 
     1329                        for (q = fit_start; q < fit_end; ++q) { 
     1330                                (*pdy_dparam_conv)[q][0] = 1.0f; 
     1331                                yfit[q] += param[0]; 
     1332                                dy[q] = y[q] - yfit[q]; 
     1333                                weight = 1.0f / (sig[q] * sig[q]); 
     1334                                alpha_weight[q] = weight; // 1 / (sig[q] * sig[q]) 
     1335                                weight *= dy[q]; 
     1336                                beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q]) 
     1337                                weight *= dy[q]; 
     1338                                *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q]) 
     1339                        } 
     1340                        break; 
     1341                } 
     1342                case NOISE_POISSON_DATA: 
     1343                { 
     1344                        for (q = fit_start; q < fit_end; ++q) { 
     1345                                (*pdy_dparam_conv)[q][0] = 1.0f; 
     1346                                yfit[q] += param[0]; 
     1347                                dy[q] = y[q] - yfit[q]; 
     1348                                weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 
     1349                                alpha_weight[q] = weight; // 1 / sig(q) 
     1350                                weight *= dy[q]; 
     1351                                beta_weight[q] = weight; // dy[q] / sig(q) 
     1352                                weight *= dy[q]; 
     1353                                *chisq += weight; // (dy[q] * dy[q]) / sig(q) 
     1354                        } 
     1355                        break; 
     1356                } 
     1357                case NOISE_POISSON_FIT: 
     1358                { 
     1359                        for (q = fit_start; q < fit_end; ++q) { 
     1360                                (*pdy_dparam_conv)[q][0] = 1.0f; 
     1361                                yfit[q] += param[0]; 
     1362                                dy[q] = y[q] - yfit[q]; 
     1363                                weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 
     1364                                alpha_weight[q] = weight; // 1 / sig(q) 
     1365                                weight *= dy[q]; 
     1366                                beta_weight[q] = weight; // dy(q) / sig(q) 
     1367                                weight *= dy[q]; 
     1368                                *chisq += weight; // (dy(q) * dy(q)) / sig(q) 
     1369                        } 
     1370                        break; 
     1371                } 
     1372                case NOISE_GAUSSIAN_FIT: 
     1373                { 
     1374                        for (q = fit_start; q < fit_end; ++q) { 
     1375                                (*pdy_dparam_conv)[q][0] = 1.0f; 
     1376                                yfit[q] += param[0]; 
     1377                                dy[q] = y[q] - yfit[q]; 
     1378                                weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 
     1379                                alpha_weight[q] = weight; // 1 / sig(q) 
     1380                                weight *= dy[q]; 
     1381                                beta_weight[q] = weight; // dy[q] / sig(q) 
     1382                                weight *= dy[q]; 
     1383                                *chisq += weight; // dy[q] / sig(q) 
     1384                        } 
     1385                        break; 
     1386                } 
     1387                case NOISE_MLE: 
     1388                { 
     1389                        for (q = fit_start; q < fit_end; ++q) { 
     1390                                (*pdy_dparam_conv)[q][0] = 1.0f; 
     1391                                yfit[q] += param[0]; 
     1392                                dy[q] = y[q] - yfit[q]; 
     1393                                weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f); 
     1394                                alpha_weight[q] = weight * y[q] / yfit[q]; 
     1395                                beta_weight[q] = dy[q] * weight; 
     1396                                if (yfit[q] > 0.0) { 
     1397                                        *chisq += (0.0f == y[q]) 
     1398                                                        ? 2.0 * yfit[q] 
     1399                                                        : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]); 
     1400                                } 
     1401                        } 
     1402                        if (*chisq <= 0.0f) { 
     1403                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
     1404                        } 
     1405                        break; 
     1406                } 
     1407                default: 
     1408                { 
     1409                        return -3; 
     1410                } 
     1411        } 
     1412 
     1413        // Check if chi square worsened: 
     1414        if (0.0f != old_chisq && *chisq >= old_chisq) { 
     1415                // don't bother to set up the matrices for solution 
     1416                return 0; 
     1417        } 
     1418 
     1419        i_free = 0; 
     1420        // for all columns 
     1421        for (i = 0; i < nparam; ++i) { 
     1422                if (paramfree[i]) { 
     1423                        j_free = 0; 
     1424                        beta_sum = 0.0f; 
     1425                        // row loop, only need to consider lower triangle 
     1426                        for (j = 0; j <= i; ++j) { 
     1427                                if (paramfree[j]) { 
     1428                                        dot_product = 0.0f; 
     1429                                        if (0 == j_free) { // true only once for each outer loop i 
     1430                                                // for all data 
     1431                                                for (k = fit_start; k < fit_end; ++k) { 
     1432                                                        dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 
     1433                                                        dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO ARG make it [i][k] and just *dy_dparam++ it. 
     1434                                                        beta_sum += dy_dparam_k_i * beta_weight[k]; 
     1435                                                } 
     1436                                        } 
     1437                                        else { 
     1438                                                // for all data 
     1439                                                for (k = fit_start; k < fit_end; ++k) { 
     1440                                                        dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; 
     1441                                                } 
     1442                                        } // k loop 
     1443                                         
     1444                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 
     1445                                        // if (i_free != j_free) { 
     1446                                        //     // matrix is symmetric 
     1447                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 
     1448                                        // } 
     1449                                        ++j_free; 
     1450                                } 
     1451                        } // j loop 
     1452                        beta[i_free] = beta_sum; 
     1453                        ++i_free; 
     1454                } 
     1455        } // i loop 
    14671456 
    14681457        return 0; 
Note: See TracChangeset for help on using the changeset viewer.