Changeset 7719 for trunk/projects
 Timestamp:
 06/01/11 17:45:09 (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/projects/slimcurve/src/main/c/EcfSingle.c
r7715 r7719 984 984 int i, j, k, l, m, mfit; 985 985 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; 986 995 987 996 for (j=0, mfit=0; j<nparam; j++) … … 989 998 mfit++; 990 999 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 1161 1154 return 0; 1162 1155 } … … 1180 1173 int i, j, k, l, m, mfit, ret; 1181 1174 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 1183 1185 /* Are we initialising? */ 1184 1186 // Malloc the arrays that will get used again in this fit via the pointers passed in … … 1304 1306 /* OK, now we've got our (possibly convolved) data, we can do the 1305 1307 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 1321 1311 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 1467 1456 1468 1457 return 0;
Note: See TracChangeset
for help on using the changeset viewer.