Changeset 7715 for trunk/projects
- Timestamp:
- 05/27/11 20:40:45 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/projects/slim-curve/src/main/c/EcfSingle.c
r7708 r7715 700 700 701 701 702 //TODO ARG deleted in my latest version703 702 int GCI_marquardt_step(float x[], float y[], int ndata, 704 703 noise_type noise, float sig[], … … 751 750 return -1; 752 751 753 //TODO ARG GCI_gauss_jordan would invert the covar matrix as a side effect754 752 /* Once converged, evaluate covariance matrix */ 755 753 if (*alambda == 0) { 756 754 if (GCI_marquardt_compute_fn_final(x, y, ndata, noise, sig, 757 param, paramfree, nparam, fitfunc, 758 yfit, dy, chisq) != 0) 759 return -3; 755 param, paramfree, nparam, fitfunc, 756 yfit, dy, chisq) != 0) 757 return -3; 758 759 760 for (j=0; j<mfit; j++) 761 for (k=0; k<mfit; k++) 762 covar[j][k] = alpha[j][k]; 763 GCI_invert(covar, mfit); 764 760 765 if (mfit < nparam) { /* no need to do this otherwise */ 761 766 GCI_covar_sort(covar, nparam, paramfree, mfit); … … 823 828 // static int mfit; // was static but now thread safe 824 829 // static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; // was static but now thread safe 830 //TODO ARG GCI_marquardt_step defines and uses int mfit = *pmfit; 825 831 826 832 if (nparam > MAXFIT) … … 868 874 return -1; 869 875 870 //TODO ARG covar needs to get inverted; previously inverted as a side effect871 876 /* Once converged, evaluate covariance matrix */ 872 877 if (*alambda == 0) { 873 878 if (GCI_marquardt_compute_fn_final_instr( 874 xincr, y, ndata, fit_start, fit_end, 875 instr, ninstr, noise, sig, 876 param, paramfree, nparam, fitfunc, 877 yfit, dy, chisq, 878 pfnvals, pdy_dparam_pure, pdy_dparam_conv, 879 pfnvals_len, pdy_dparam_nparam_size) != 0) 880 return -3; 879 xincr, y, ndata, fit_start, fit_end, 880 instr, ninstr, noise, sig, 881 param, paramfree, nparam, fitfunc, 882 yfit, dy, chisq, 883 pfnvals, pdy_dparam_pure, pdy_dparam_conv, 884 pfnvals_len, pdy_dparam_nparam_size) != 0) 885 return -3; 886 887 for (j=0; j<(*pmfit); j++) 888 for (k=0; k<(*pmfit); k++) 889 covar[j][k] = alpha[j][k]; 890 GCI_invert(covar, (*pmfit)); 891 892 881 893 if (*pmfit < nparam) { /* no need to do this otherwise */ 882 894 GCI_covar_sort(covar, nparam, paramfree, *pmfit); … … 962 974 two variants of the Marquardt function. 963 975 */ 964 965 //TODO ARG deleted in my version; needs to get de-NRed!!!!966 976 int GCI_marquardt_compute_fn(float x[], float y[], int ndata, 967 977 noise_type noise, float sig[], … … 973 983 { 974 984 int i, j, k, l, m, mfit; 975 float wt, sig2i, y_ymod, dy_dparam[MAX FIT];985 float wt, sig2i, y_ymod, dy_dparam[MAXBINS][MAXFIT]; 976 986 977 987 for (j=0, mfit=0; j<nparam; j++) … … 998 1008 case NOISE_CONST: 999 1009 { 1010 float i_sig_0_squared = 1.0 / (sig[0] * sig[0]); 1000 1011 for (q = 0; q < ndata; ++q) { 1001 (*fitfunc)(x[q], param, &yfit[q], dy_dparam , nparam);1012 (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 1002 1013 yfit[q] += param[0]; 1014 dy_dparam[q][0] = 1.0; 1003 1015 dy[q] = y[q] - yfit[q]; 1004 weight = 1.0f; //TODO ARG this should be 1.0f / sig[0] !1005 alpha_weight[q] = weight; // 1 1016 weight = i_sig_0_squared; 1017 alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]) 1006 1018 weight *= dy[q]; 1007 beta_weight[q] = weight; // dy[q] 1019 beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]) 1008 1020 weight *= dy[q]; 1009 *chisq += weight; // dy[q] * dy[q]1021 *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]) 1010 1022 } 1011 1023 break; … … 1014 1026 { 1015 1027 for (q = 0; q < ndata; ++q) { 1016 (*fitfunc)(x[q], param, &yfit[q], dy_dparam , nparam);1028 (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 1017 1029 yfit[q] += param[0]; 1030 dy_dparam[q][0] = 1.0; 1018 1031 dy[q] = y[q] - yfit[q]; 1019 1032 weight = 1.0f / (sig[q] * sig[q]); … … 1029 1042 { 1030 1043 for (q = 0; q < ndata; ++q) { 1031 (*fitfunc)(x[q], param, &yfit[q], dy_dparam , nparam);1044 (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 1032 1045 yfit[q] += param[0]; 1046 dy_dparam[q][0] = 1.0; 1033 1047 dy[q] = y[q] - yfit[q]; 1034 1048 weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); … … 1044 1058 { 1045 1059 for (q = 0; q < ndata; ++q) { 1046 (*fitfunc)(x[q], param, &yfit[q], dy_dparam , nparam);1060 (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 1047 1061 yfit[q] += param[0]; 1062 dy_dparam[q][0] = 1.0; 1048 1063 dy[q] = y[q] - yfit[q]; 1049 1064 weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); … … 1059 1074 { 1060 1075 for (q = 0; q < ndata; ++q) { 1061 (*fitfunc)(x[q], param, &yfit[q], dy_dparam , nparam);1076 (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 1062 1077 yfit[q] += param[0]; 1078 dy_dparam[q][0] = 1.0; 1063 1079 dy[q] = y[q] - yfit[q]; 1064 1080 weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); … … 1067 1083 beta_weight[q] = weight; // dy[q] / sig(q) 1068 1084 weight *= dy[q]; 1069 *chisq += weight; 1085 *chisq += weight; // dy[q] / sig(q) 1070 1086 } 1071 1087 break; … … 1074 1090 { 1075 1091 for (q = 0; q < ndata; ++q) { 1076 (*fitfunc)(x[q], param, &yfit[q], dy_dparam , nparam);1092 (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam); 1077 1093 yfit[q] += param[0]; 1094 dy_dparam[q][0] = 1.0; 1078 1095 dy[q] = y[q] - yfit[q]; 1079 weight = (yfit[q] > 1 ? 1.0f / yfit[ i] : 1.0f);1080 alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus?1096 weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f); 1097 alpha_weight[q] = weight * y[q] / yfit[q]; 1081 1098 beta_weight[q] = dy[q] * weight; 1082 *chisq += (0.0f == y[q]) 1083 ? 2.0 * yfit[q] 1084 : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]); 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 } 1085 1104 } 1086 1105 if (*chisq <= 0.0f) { 1087 1106 *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 1088 1107 } … … 1111 1130 if (paramfree[j]) { 1112 1131 dot_product = 0.0f; 1113 if (0 == j_free) { 1132 if (0 == j_free) { // true only once for each outer loop i 1114 1133 // for all data 1115 1134 for (k = 0; k < ndata; ++k) { 1116 //TODO ARG just to get this to compile, for now: 1117 //TODO ARG from the _instr version: dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 1118 //TODO ARG from the instr version dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it. 1119 //TODO ARG from the instr version beta_sum += dy_dparam_k_i * beta_weight[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 1120 1137 } 1121 1138 } … … 1123 1140 // for all data 1124 1141 for (k = 0; k < ndata; ++k) { 1125 //TODO ARG from the instr version: dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k];1142 dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k]; 1126 1143 } 1127 1144 } // k loop 1128 1145 1129 alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 1130 // if (i_free != j_free) { 1131 // / / matrixis symmetric1146 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 1132 1149 // alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 1133 1150 // } … … 1138 1155 ++i_free; 1139 1156 } 1140 //else printf("param not free %d\n", i);1141 1157 } // i loop 1142 } 1143 1144 return 0; 1145 1146 1147 /* Initialise (symmetric) alpha, beta */ 1148 //TODO ARG FRI 1149 //for (j=0; j<mfit; j++) { 1150 // for (k=0; k<=j; k++) 1151 // alpha[j][k] = 0.0; 1152 // beta[j] = 0.0; 1153 //} 1154 1155 /* Calculation of the fitting data will depend upon the type of 1156 noise and the type of instrument response */ 1157 1158 /* There's no convolution involved in this function. This is then 1159 fairly straightforward, depending only upon the type of noise 1160 present. Since we calculate the standard deviation at every 1161 point in a different way depending upon the type of noise, we 1162 will place our switch(noise) around the entire loop. */ 1163 1164 switch (noise) { 1165 case NOISE_CONST: 1166 *chisq = 0.0; 1167 /* Summation loop over all data */ 1168 for (i=0; i<ndata; i++) { 1169 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1170 yfit[i] += param[0]; 1171 dy_dparam[0] = 1; 1172 dy[i] = y[i] - yfit[i]; 1173 for (j=0, l=0; l<nparam; l++) { 1174 if (paramfree[l]) { 1175 wt = dy_dparam[l]; /* taken away the *sig2i from here */ 1176 for (k=0, m=0; m<=l; m++) 1177 if (paramfree[m]) 1178 alpha[j][k++] += wt * dy_dparam[m]; 1179 beta[j] += dy[i] * wt; 1180 j++; 1181 } 1182 } 1183 /* And find chi^2 */ 1184 *chisq += dy[i] * dy[i]; 1185 } 1186 1187 /* Now divide everything by sigma^2 */ 1188 sig2i = 1.0 / (sig[0] * sig[0]); 1189 *chisq *= sig2i; 1190 for (j=0; j<mfit; j++) { 1191 for (k=0; k<=j; k++) 1192 alpha[j][k] *= sig2i; 1193 beta[j] *= sig2i; 1194 } 1195 break; 1196 1197 case NOISE_GIVEN: /* This is essentially the NR version */ 1198 *chisq = 0.0; 1199 /* Summation loop over all data */ 1200 for (i=0; i<ndata; i++) { 1201 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1202 yfit[i] += param[0]; 1203 dy_dparam[0] = 1; 1204 sig2i = 1.0 / (sig[i] * sig[i]); 1205 dy[i] = y[i] - yfit[i]; 1206 for (j=0, l=0; l<nparam; l++) { 1207 if (paramfree[l]) { 1208 wt = dy_dparam[l] * sig2i; 1209 for (k=0, m=0; m<=l; m++) 1210 if (paramfree[m]) 1211 alpha[j][k++] += wt * dy_dparam[m]; 1212 beta[j] += wt * dy[i]; 1213 j++; 1214 } 1215 } 1216 /* And find chi^2 */ 1217 *chisq += dy[i] * dy[i] * sig2i; 1218 } 1219 break; 1220 1221 case NOISE_POISSON_DATA: 1222 *chisq = 0.0; 1223 /* Summation loop over all data */ 1224 for (i=0; i<ndata; i++) { 1225 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1226 yfit[i] += param[0]; 1227 dy_dparam[0] = 1; 1228 sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15); 1229 dy[i] = y[i] - yfit[i]; 1230 for (j=0, l=0; l<nparam; l++) { 1231 if (paramfree[l]) { 1232 wt = dy_dparam[l] * sig2i; 1233 for (k=0, m=0; m<=l; m++) 1234 if (paramfree[m]) 1235 alpha[j][k++] += wt * dy_dparam[m]; 1236 beta[j] += wt * dy[i]; 1237 j++; 1238 } 1239 } 1240 /* And find chi^2 */ 1241 *chisq += dy[i] * dy[i] * sig2i; 1242 } 1243 break; 1244 1245 case NOISE_POISSON_FIT: 1246 *chisq = 0.0; 1247 // Summation loop over all data 1248 for (i=0; i<ndata; i++) { 1249 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1250 yfit[i] += param[0]; 1251 dy_dparam[0] = 1; 1252 sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15); 1253 dy[i] = y[i] - yfit[i]; 1254 for (j=0, l=0; l<nparam; l++) { 1255 if (paramfree[l]) { 1256 wt = dy_dparam[l] * sig2i; 1257 for (k=0, m=0; m<=l; m++) 1258 if (paramfree[m]) 1259 alpha[j][k++] += wt * dy_dparam[m]; 1260 beta[j] += wt * dy[i]; 1261 j++; 1262 } 1263 } 1264 // And find chi^2 1265 *chisq += dy[i] * dy[i] * sig2i; 1266 } 1267 break; 1268 1269 case NOISE_GAUSSIAN_FIT: 1270 *chisq = 0.0; 1271 // Summation loop over all data 1272 for (i=0; i<ndata; i++) { 1273 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1274 yfit[i] += param[0]; 1275 dy_dparam[0] = 1; 1276 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1277 dy[i] = y[i] - yfit[i]; 1278 for (j=0, l=0; l<nparam; l++) { 1279 if (paramfree[l]) { 1280 wt = dy_dparam[l] * sig2i; 1281 for (k=0, m=0; m<=l; m++) 1282 if (paramfree[m]) 1283 alpha[j][k++] += wt * dy_dparam[m]; 1284 beta[j] += wt * dy[i]; 1285 j++; 1286 } 1287 } 1288 // And find chi^2 1289 *chisq += dy[i] * dy[i] * sig2i; 1290 } 1291 break; 1292 1293 case NOISE_MLE: 1294 *chisq = 0.0; 1295 /* Summation loop over all data */ 1296 for (i=0; i<ndata; i++) { 1297 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1298 yfit[i] += param[0]; 1299 dy_dparam[0] = 1; 1300 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1301 dy[i] = y[i] - yfit[i]; 1302 y_ymod=y[i]/yfit[i]; 1303 for (j=0, l=0; l<nparam; l++) { 1304 if (paramfree[l]) { 1305 wt = dy_dparam[l] * sig2i; 1306 for (k=0, m=0; m<=l; m++) 1307 if (paramfree[m]) 1308 alpha[j][k++] += wt*dy_dparam[m]*y_ymod; //wt * dy_dparam[m]; 1309 beta[j] += wt * dy[i]; 1310 j++; 1311 } 1312 } 1313 // And find chi^2 1314 if (yfit[i]<=0.0) 1315 ; // do nothing 1316 else if (y[i]==0.0) 1317 *chisq += 2.0*yfit[i]; // to avoid NaN from log 1318 else 1319 *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; 1320 } 1321 if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve 1322 break; 1323 1324 default: 1325 return -3; 1326 /* break; */ // (unreachable code) 1327 } 1328 1329 /* Fill in the symmetric side */ 1330 for (j=1; j<mfit; j++) 1331 for (k=0; k<j; k++) 1332 alpha[k][j] = alpha[j][k]; 1158 1159 } 1333 1160 1334 1161 return 0; … … 1360 1187 /* do any necessary initialisation */ 1361 1188 if ((*pfnvals_len) < ndata) { /* we will need this length for 1362 1189 the final full computation */ 1363 1190 if ((*pfnvals_len)) { 1364 1191 free((*pfnvals)); … … 1401 1228 if (paramfree[j]) mfit++; 1402 1229 1403 //TODO ARG first change:1404 // /* Initialise (symmetric) alpha, beta */1405 // for (j=0; j<mfit; j++) {1406 // for (k=0; k<=j; k++)1407 // alpha[j][k] = 0.0;1408 // beta[j] = 0.0;1409 // }1410 1411 1230 /* Calculation of the fitting data will depend upon the type of 1412 1231 noise and the type of instrument response */ … … 1500 1319 *chisq = 0.0f; 1501 1320 1502 1321 switch (noise) { 1503 1322 case NOISE_CONST: 1504 1323 { … … 1507 1326 yfit[q] += param[0]; 1508 1327 dy[q] = y[q] - yfit[q]; 1509 weight = 1.0f ; //TODO ARG this should be 1.0f / sig[0] !1510 alpha_weight[q] = weight; // 1 1328 weight = 1.0f / sig[0]; 1329 alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]); 1511 1330 weight *= dy[q]; 1512 beta_weight[q] = weight; // dy[q] 1331 beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]); 1513 1332 weight *= dy[q]; 1514 *chisq += weight; // dy[q] * dy[q]1333 *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]); 1515 1334 } 1516 1335 break; … … 1572 1391 beta_weight[q] = weight; // dy[q] / sig(q) 1573 1392 weight *= dy[q]; 1574 *chisq += weight; 1393 *chisq += weight; // dy[q] / sig(q) 1575 1394 } 1576 1395 break; … … 1583 1402 dy[q] = y[q] - yfit[q]; 1584 1403 weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f); 1585 alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus?1404 alpha_weight[q] = weight * y[q] / yfit[q]; 1586 1405 beta_weight[q] = dy[q] * weight; 1587 *chisq += (0.0f == y[q]) 1588 ? 2.0 * yfit[q] 1589 : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]); 1590 } 1591 if (*chisq <= 0.0f) { 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) { 1592 1413 *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 1593 1414 } … … 1616 1437 if (paramfree[j]) { 1617 1438 dot_product = 0.0f; 1618 if (0 == j_free) { 1439 if (0 == j_free) { // true only once for each outer loop i 1619 1440 // for all data 1620 1441 for (k = fit_start; k < fit_end; ++k) { 1621 1442 dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 1622 dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it.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. 1623 1444 beta_sum += dy_dparam_k_i * beta_weight[k]; 1624 1445 } … … 1642 1463 ++i_free; 1643 1464 } 1644 //else printf("param not free %d\n", i);1645 1465 } // i loop 1646 1466 }
Note: See TracChangeset
for help on using the changeset viewer.