Changeset 7715 for trunk/projects


Ignore:
Timestamp:
05/27/11 20:40:45 (9 years ago)
Author:
aivar
Message:

Got GCI_marquardt_compute_fn working, in testing it agrees with the old code. Computing covar as the inverse of the alpha matrix.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r7708 r7715  
    700700 
    701701 
    702 //TODO ARG deleted in my latest version 
    703702int GCI_marquardt_step(float x[], float y[], int ndata, 
    704703                                        noise_type noise, float sig[], 
    … …  
    751750                return -1; 
    752751 
    753         //TODO ARG GCI_gauss_jordan would invert the covar matrix as a side effect 
    754752        /* Once converged, evaluate covariance matrix */ 
    755753        if (*alambda == 0) { 
    756754                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 
    760765                if (mfit < nparam) {  /* no need to do this otherwise */ 
    761766                        GCI_covar_sort(covar, nparam, paramfree, mfit); 
    … …  
    823828//      static int mfit;   // was static but now thread safe 
    824829//      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; 
    825831 
    826832        if (nparam > MAXFIT) 
    … …  
    868874                return -1; 
    869875 
    870 //TODO ARG covar needs to get inverted; previously inverted as a side effect 
    871876        /* Once converged, evaluate covariance matrix */ 
    872877        if (*alambda == 0) { 
    873878                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 
    881893                if (*pmfit < nparam) {  /* no need to do this otherwise */ 
    882894                        GCI_covar_sort(covar, nparam, paramfree, *pmfit); 
    … …  
    962974   two variants of the Marquardt function. 
    963975*/ 
    964  
    965 //TODO ARG deleted in my version; needs to get de-NRed!!!! 
    966976int GCI_marquardt_compute_fn(float x[], float y[], int ndata, 
    967977                                         noise_type noise, float sig[], 
    … …  
    973983{ 
    974984        int i, j, k, l, m, mfit; 
    975         float wt, sig2i, y_ymod, dy_dparam[MAXFIT]; 
     985        float wt, sig2i, y_ymod, dy_dparam[MAXBINS][MAXFIT]; 
    976986 
    977987        for (j=0, mfit=0; j<nparam; j++) 
    … …  
    9981008            case NOISE_CONST: 
    9991009            { 
     1010                float i_sig_0_squared = 1.0 / (sig[0] * sig[0]); 
    10001011                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); 
    10021013                    yfit[q] += param[0]; 
     1014                    dy_dparam[q][0] = 1.0; 
    10031015                    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]) 
    10061018                    weight *= dy[q]; 
    1007                     beta_weight[q] = weight; // dy[q] 
     1019                    beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]) 
    10081020                    weight *= dy[q]; 
    1009                     *chisq += weight; // dy[q] * dy[q] 
     1021                    *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]) 
    10101022                } 
    10111023                break; 
    … …  
    10141026            { 
    10151027                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); 
    10171029                    yfit[q] += param[0]; 
     1030                    dy_dparam[q][0] = 1.0; 
    10181031                    dy[q] = y[q] - yfit[q]; 
    10191032                    weight = 1.0f / (sig[q] * sig[q]); 
    … …  
    10291042            { 
    10301043                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); 
    10321045                    yfit[q] += param[0]; 
     1046                    dy_dparam[q][0] = 1.0; 
    10331047                    dy[q] = y[q] - yfit[q]; 
    10341048                    weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 
    … …  
    10441058            { 
    10451059                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); 
    10471061                    yfit[q] += param[0]; 
     1062                    dy_dparam[q][0] = 1.0; 
    10481063                    dy[q] = y[q] - yfit[q]; 
    10491064                    weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 
    … …  
    10591074            { 
    10601075                 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); 
    10621077                    yfit[q] += param[0]; 
     1078                    dy_dparam[q][0] = 1.0; 
    10631079                    dy[q] = y[q] - yfit[q]; 
    10641080                    weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 
    … …  
    10671083                    beta_weight[q] = weight; // dy[q] / sig(q) 
    10681084                    weight *= dy[q]; 
    1069                     *chisq += weight; 
     1085                    *chisq += weight; // dy[q] / sig(q) 
    10701086                 } 
    10711087                 break; 
    … …  
    10741090            { 
    10751091                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); 
    10771093                    yfit[q] += param[0]; 
     1094                    dy_dparam[q][0] = 1.0; 
    10781095                    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]; 
    10811098                    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                    } 
    10851104                } 
    1086                 if (*chisq <= 0.0f) { 
     1105                if (*chisq <= 0.0f) { 
    10871106                    *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
    10881107                } 
    … …  
    11111130                    if (paramfree[j]) { 
    11121131                        dot_product = 0.0f; 
    1113                         if (0 == j_free) { 
     1132                        if (0 == j_free) { // true only once for each outer loop i 
    11141133                            // for all data 
    11151134                            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 
    11201137                            } 
    11211138                        } 
    … …  
    11231140                            // for all data 
    11241141                            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]; 
    11261143                            } 
    11271144                        } // k loop 
    11281145 
    1129                         alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 
    1130                        // if (i_free != j_free) { 
    1131                        //     // matrix is symmetric 
     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 
    11321149                       //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 
    11331150                       // } 
    … …  
    11381155                ++i_free; 
    11391156            } 
    1140             //else printf("param not free %d\n", i); 
    11411157        } // 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        } 
    13331160 
    13341161        return 0; 
    … …  
    13601187                /* do any necessary initialisation */ 
    13611188                if ((*pfnvals_len) < ndata) {  /* we will need this length for 
    1362                                                                           the final full computation */ 
     1189                                                                                the final full computation */ 
    13631190                        if ((*pfnvals_len)) { 
    13641191                                free((*pfnvals)); 
    … …  
    14011228                if (paramfree[j]) mfit++; 
    14021229 
    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  
    14111230        /* Calculation of the fitting data will depend upon the type of 
    14121231           noise and the type of instrument response */ 
    … …  
    15001319        *chisq = 0.0f; 
    15011320 
    1502                 switch (noise) { 
     1321        switch (noise) { 
    15031322            case NOISE_CONST: 
    15041323            { 
    … …  
    15071326                    yfit[q] += param[0]; 
    15081327                    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]); 
    15111330                    weight *= dy[q]; 
    1512                     beta_weight[q] = weight; // dy[q] 
     1331                    beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]); 
    15131332                    weight *= dy[q]; 
    1514                     *chisq += weight; // dy[q] * dy[q] 
     1333                    *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]); 
    15151334                } 
    15161335                break; 
    … …  
    15721391                    beta_weight[q] = weight; // dy[q] / sig(q) 
    15731392                    weight *= dy[q]; 
    1574                     *chisq += weight; 
     1393                    *chisq += weight; // dy[q] / sig(q) 
    15751394                 } 
    15761395                 break; 
    … …  
    15831402                    dy[q] = y[q] - yfit[q]; 
    15841403                    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]; 
    15861405                    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) { 
    15921413                    *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
    15931414                } 
    … …  
    16161437                    if (paramfree[j]) { 
    16171438                        dot_product = 0.0f; 
    1618                         if (0 == j_free) { 
     1439                        if (0 == j_free) { // true only once for each outer loop i 
    16191440                            // for all data 
    16201441                            for (k = fit_start; k < fit_end; ++k) { 
    16211442                                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. 
    16231444                                beta_sum += dy_dparam_k_i * beta_weight[k]; 
    16241445                            } 
    … …  
    16421463                ++i_free; 
    16431464            } 
    1644             //else printf("param not free %d\n", i); 
    16451465        } // i loop 
    16461466        } 
Note: See TracChangeset for help on using the changeset viewer.