Changeset 7708 for trunk/projects


Ignore:
Timestamp:
05/20/11 22:24:58 (8 years ago)
Author:
aivar
Message:

Work in progress. Working on GCI_marquardt_compute_fn, code is commented out so it compiles.

Tidied up a bit.

Unexpected behavior in EcfUtil.c: Optimization looks good in generated assembly code but runs much slower.

Location:
trunk/projects/slim-curve/src/main/c
Files:
4 edited

Legend:

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

    r7706 r7708  
    1616#define MAXITERS 80 
    1717#define MAXREFITS 10 
     18#define MAXBINS 1024 /* Maximum number of lifetime bins; saves dynamic allocation of small arrays */ 
    1819 
    1920/* Functions from EcfSingle.c */ 
     
    2425                                         void (*fitfunc)(float, float [], float *, float [], int), 
    2526                                         float yfit[], float dy[], 
    26                                          float **alpha, float beta[], float *chisq, 
     27                                         float **alpha, float beta[], float *chisq, float old_chisq, 
    2728                                         float alambda); 
    2829int GCI_marquardt_compute_fn_instr(float xincr, float y[], int ndata, 
  • trunk/projects/slim-curve/src/main/c/EcfSingle.c

    r7706 r7708  
    44 
    55   This file contains functions for single transient analysis. 
    6    Utility code is found in EcfUtil.c. 
     6   Utility code is found in EcfUtil.c and global analysis code in 
     7   EcfGlobal.c. 
    78*/ 
    89 
     
    1314 
    1415/* Predeclarations */ 
     16int GCI_marquardt_step(float x[], float y[], int ndata, 
     17                                        noise_type noise, float sig[], 
     18                                        float param[], int paramfree[], int nparam, 
     19                                        restrain_type restrain, 
     20                                        void (*fitfunc)(float, float [], float *, float [], int), 
     21                                        float yfit[], float dy[], 
     22                                        float **covar, float **alpha, float *chisq, 
     23                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam); 
    1524int GCI_marquardt_step_instr(float xincr, float y[], 
    1625                                        int ndata, int fit_start, int fit_end, 
     
    2231                                        float yfit[], float dy[], 
    2332                                        float **covar, float **alpha, float *chisq, 
    24                                         float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,         
     33                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam, 
    2534                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, 
    2635                                        int *pfnvals_len, int *pdy_dparam_nparam_size); 
     36int GCI_marquardt_estimate_errors(float **alpha, int nparam, int mfit, 
     37                                                                  float d[], float **v, float interval); 
    2738 
    2839/* Globals */ 
     
    349360                                                          float *chisq, float chisq_target) 
    350361{ 
    351         int tries=1, division=3;                 // the data  
     362        int tries=1, division=3;                 // the data 
    352363        float local_chisq=3.0e38, oldChisq=3.0e38, oldZ, oldA, oldTau, *validFittedArray; // local_chisq a very high float but below oldChisq 
    353                  
     364 
    354365        if (fitted==NULL)   // we require chisq but have not supplied a "fitted" array so must malloc one 
    355366        { 
     
    357368        } 
    358369        else validFittedArray = fitted; 
    359          
     370 
    360371        if (instr==NULL)           // no instrument/prompt has been supplied 
    361372        { 
     
    363374                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division); 
    364375 
    365                 while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS)    
     376                while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS) 
    366377                { 
    367378                        oldChisq = local_chisq; 
     
    369380                        oldA = *A; 
    370381                        oldTau = *tau; 
    371 //                      division++;                       
    372                         division+=division/3;                     
     382//                      division++; 
     383                        division+=division/3; 
    373384                        tries++; 
    374385                        GCI_triple_integral(xincr, y, fit_start, fit_end, noise, sig, 
    375386                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division); 
    376                 }                                                  
     387                } 
    377388        } 
    378389        else 
     
    381392                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division); 
    382393 
    383                 while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS)    
     394                while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS) 
    384395                { 
    385                         oldChisq = local_chisq;  
     396                        oldChisq = local_chisq; 
    386397                        oldZ = *Z; 
    387398                        oldA = *A; 
    388399                        oldTau = *tau; 
    389 //                      division++;                       
     400//                      division++; 
    390401                        division+=division/3; 
    391402                        tries++; 
     
    393404                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division); 
    394405 
    395                 }                                                  
     406                } 
    396407        } 
    397408 
     
    403414                *tau = oldTau; 
    404415        } 
    405          
     416 
    406417        if (chisq!=NULL) *chisq = local_chisq; 
    407418 
    408         if (fitted==NULL)   
     419        if (fitted==NULL) 
    409420        { 
    410421                free (validFittedArray); 
     
    454465   noise=NOISE_POISSON_FIT, the variances are taken to be 
    455466   max(yfit[i],15). If noise=NOISE_GAUSSIAN_FIT, the variances are taken to be 
    456    yfit[i] and if noise=NOISE_MLE then a optimisation is for the maximum  
    457    likelihood approximation (Laurence and Chromy in press 2010 and  
     467   yfit[i] and if noise=NOISE_MLE then a optimisation is for the maximum 
     468   likelihood approximation (Laurence and Chromy in press 2010 and 
    458469   G. Nishimura, and M. Tamura Phys Med Biol 2005). 
    459     
     470 
    460471   The input array paramfree[0..nparam-1] indicates 
    461472   by nonzero entries those components that should be held fixed at 
     
    488499*/ 
    489500 
    490 /* This functions does the whole job */ 
    491  
     501/* These two functions do the whole job */ 
    492502int GCI_marquardt(float x[], float y[], int ndata, 
    493503                                  noise_type noise, float sig[], 
     
    499509                                  float chisq_delta, float chisq_percent, float **erraxes) 
    500510{ 
    501         // Empty fn to allow compile in TRI2 
     511        float alambda, ochisq; 
     512        int mfit; 
     513        float evals[MAXFIT]; 
     514        int i, k, itst, itst_max; 
     515 
     516        float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 
     517 
     518        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; 
     519 
     520        mfit = 0; 
     521        for (i=0; i<nparam; i++) { 
     522                if (paramfree[i]) mfit++; 
     523        } 
     524 
     525        alambda = -1; 
     526        if (GCI_marquardt_step(x, y, ndata, noise, sig, 
     527                                                   param, paramfree, nparam, restrain, 
     528                                                   fitfunc, fitted, residuals, 
     529                                                   covar, alpha, chisq, &alambda, 
     530                                                   &mfit, &ochisq2, paramtry, beta, dparam) != 0) { 
     531                return -1; 
     532        } 
     533 
     534        k = 1;  /* Iteration counter */ 
     535        itst = 0; 
     536        for (;;) { 
     537                k++; 
     538                if (k > MAXITERS) { 
     539                        return -2; 
     540                } 
     541 
     542                ochisq = *chisq; 
     543                if (GCI_marquardt_step(x, y, ndata, noise, sig, 
     544                                                           param, paramfree, nparam, restrain, 
     545                                                           fitfunc, fitted, residuals, 
     546                                                           covar, alpha, chisq, &alambda, 
     547                                                           &mfit, &ochisq2, paramtry, beta, dparam) != 0) { 
     548                        return -3; 
     549                } 
     550 
     551                if (*chisq > ochisq) 
     552                        itst = 0; 
     553                else if (ochisq - *chisq < chisq_delta) 
     554                        itst++; 
     555 
     556                if (itst < itst_max) continue; 
     557 
     558                /* Endgame */ 
     559                alambda = 0.0; 
     560                if (GCI_marquardt_step(x, y, ndata, noise, sig, 
     561                                                           param, paramfree, nparam, restrain, 
     562                                                           fitfunc, fitted, residuals, 
     563                                                           covar, alpha, chisq, &alambda, 
     564                                                           &mfit, &ochisq2, paramtry, beta, dparam) != 0) { 
     565                        return -4; 
     566                } 
     567 
     568                if (erraxes == NULL) 
     569                        return k; 
     570 
     571                //TODO ARG 
     572                //if (GCI_marquardt_estimate_errors(alpha, nparam, mfit, evals, 
     573                //                                                                erraxes, chisq_percent) != 0) { 
     574                //      return -5; 
     575                //} 
     576 
     577                break;  /* We're done now */ 
     578        } 
     579 
     580        return k; 
    502581} 
    503582 
     
    529608        float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL; 
    530609        int fnvals_len=0, dy_dparam_nparam_size=0; 
    531         float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];  
     610        float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 
    532611 
    533612        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; 
     
    545624                                                                 param, paramfree, nparam, restrain, 
    546625                                                                 fitfunc, fitted, residuals, 
    547                                                                  covar, alpha, chisq, &alambda,  
     626                                                                 covar, alpha, chisq, &alambda, 
    548627                                                                 &mfit2, &ochisq2, paramtry, beta, dparam, 
    549628                                                                 &fnvals, &dy_dparam_pure, &dy_dparam_conv, 
     
    569648                                                                         param, paramfree, nparam, restrain, 
    570649                                                                         fitfunc, fitted, residuals, 
    571                                                                          covar, alpha, chisq, &alambda,  
     650                                                                         covar, alpha, chisq, &alambda, 
    572651                                                                         &mfit2, &ochisq2, paramtry, beta, dparam, 
    573652                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv, 
     
    592671                                                                         param, paramfree, nparam, restrain, 
    593672                                                                         fitfunc, fitted, residuals, 
    594                                                                          covar, alpha, chisq, &alambda,  
     673                                                                         covar, alpha, chisq, &alambda, 
    595674                                                                         &mfit2, &ochisq2, paramtry, beta, dparam, 
    596675                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv, 
     
    605684                } 
    606685 
     686//TODO ARG this estimate errors call was deleted in my latest version 
     687        //      if (GCI_marquardt_estimate_errors(alpha, nparam, mfit, evals, 
     688        //                                                                        erraxes, chisq_percent) != 0) { 
     689        //              do_frees 
     690        //              return -5; 
     691        //      } 
     692 
    607693                break;  /* We're done now */ 
    608694        } 
     
    612698} 
    613699#undef do_frees 
     700 
     701 
     702//TODO ARG deleted in my latest version 
     703int GCI_marquardt_step(float x[], float y[], int ndata, 
     704                                        noise_type noise, float sig[], 
     705                                        float param[], int paramfree[], int nparam, 
     706                                        restrain_type restrain, 
     707                                        void (*fitfunc)(float, float [], float *, float [], int), 
     708                                        float yfit[], float dy[], 
     709                                        float **covar, float **alpha, float *chisq, 
     710                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam) 
     711{ 
     712        int j, k, l, ret; 
     713//      static int mfit;   // was static but now thread safe 
     714//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe 
     715        int mfit = *pmfit; 
     716        float ochisq = *pochisq; 
     717 
     718        if (nparam > MAXFIT) 
     719                return -1; 
     720 
     721        /* Initialisation */ 
     722        /* We assume we're given sensible starting values for param[] */ 
     723        if (*alambda < 0.0) { 
     724                for (mfit=0, j=0; j<nparam; j++) 
     725                        if (paramfree[j]) 
     726                                mfit++; 
     727 
     728                if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig, 
     729                                                                         param, paramfree, nparam, fitfunc, 
     730                                                                         yfit, dy, 
     731                                                                         alpha, beta, chisq, 0.0, *alambda) != 0) 
     732                        return -2; 
     733 
     734                *alambda = 0.001; 
     735                ochisq = (*chisq); 
     736                for (j=0; j<nparam; j++) 
     737                        paramtry[j] = param[j]; 
     738        } 
     739 
     740        /* Alter linearised fitting matrix by augmenting diagonal elements */ 
     741        for (j=0; j<mfit; j++) { 
     742                for (k=0; k<mfit; k++) 
     743                        covar[j][k] = alpha[j][k]; 
     744 
     745                covar[j][j] = alpha[j][j] * (1.0 + (*alambda)); 
     746                dparam[j] = beta[j]; 
     747        } 
     748 
     749        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */ 
     750        if (GCI_solve(covar, mfit, dparam) != 0) 
     751                return -1; 
     752 
     753        //TODO ARG GCI_gauss_jordan would invert the covar matrix as a side effect 
     754        /* Once converged, evaluate covariance matrix */ 
     755        if (*alambda == 0) { 
     756                if (GCI_marquardt_compute_fn_final(x, y, ndata, noise, sig, 
     757                                                                                   param, paramfree, nparam, fitfunc, 
     758                                                                                   yfit, dy, chisq) != 0) 
     759                        return -3; 
     760                if (mfit < nparam) {  /* no need to do this otherwise */ 
     761                        GCI_covar_sort(covar, nparam, paramfree, mfit); 
     762                        GCI_covar_sort(alpha, nparam, paramfree, mfit); 
     763                } 
     764                return 0; 
     765        } 
     766 
     767        /* Did the trial succeed? */ 
     768        for (j=0, l=0; l<nparam; l++) 
     769                if (paramfree[l]) 
     770                        paramtry[l] = param[l] + dparam[j++]; 
     771 
     772        if (restrain == ECF_RESTRAIN_DEFAULT) 
     773                ret = check_ecf_params (paramtry, nparam, fitfunc); 
     774        else 
     775                ret = check_ecf_user_params (paramtry, nparam, fitfunc); 
     776 
     777        if (ret != 0) { 
     778                /* Bad parameters, increase alambda and return */ 
     779                *alambda *= 10.0; 
     780                return 0; 
     781        } 
     782 
     783        if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig, 
     784                                                                 paramtry, paramfree, nparam, fitfunc, 
     785                                                                 yfit, dy, covar, dparam, 
     786                                                                 chisq, ochisq, *alambda) != 0) 
     787                return -2; 
     788 
     789        /* Success, accept the new solution */ 
     790        if (*chisq < ochisq) { 
     791                *alambda *= 0.1; 
     792                ochisq = *chisq; 
     793                for (j=0; j<mfit; j++) { 
     794                        for (k=0; k<mfit; k++) 
     795                                alpha[j][k] = covar[j][k]; 
     796                        beta[j] = dparam[j]; 
     797                } 
     798                for (l=0; l<nparam; l++) 
     799                        param[l] = paramtry[l]; 
     800        } else { /* Failure, increase alambda and return */ 
     801                *alambda *= 10.0; 
     802                *chisq = ochisq; 
     803        } 
     804 
     805        return 0; 
     806} 
    614807 
    615808 
     
    623816                                        float yfit[], float dy[], 
    624817                                        float **covar, float **alpha, float *chisq, 
    625                                         float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,         
     818                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam, 
    626819                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, 
    627820                                        int *pfnvals_len, int *pdy_dparam_nparam_size) 
     
    631824//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe 
    632825 
    633         /*if (nparam > MAXFIT) { 
    634             printf("GCI_marq_step_instr returns -10\n"); 
     826        if (nparam > MAXFIT) 
    635827                return -10; 
    636         } 
    637         if (xincr <= 0) { 
    638             printf("GCI_marq_step_instr returns -11\n"); 
     828        if (xincr <= 0) 
    639829                return -11; 
    640         } 
    641         if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) { 
    642             printf("GCI_marq_step_instr returns -12\n"); 
     830        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) 
    643831                return -12; 
    644         }*/ 
    645832 
    646833        /* Initialisation */ 
     
    656843                                                                                   param, paramfree, nparam, fitfunc, 
    657844                                                                                   yfit, dy, alpha, beta, chisq, 0.0, 
    658                                                                                    *alambda,     
     845                                                                                   *alambda, 
    659846                                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv, 
    660                                                                                         pfnvals_len, pdy_dparam_nparam_size) != 0) { 
     847                                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0) 
    661848                        return -2; 
    662                 } 
    663849 
    664850                *alambda = 0.001; 
     
    678864        } 
    679865 
    680         /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */ 
    681         if (GCI_solve(covar, *pmfit, dparam) != 0) { 
     866        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */ 
     867        if (GCI_solve(covar, *pmfit, dparam) != 0) 
    682868                return -1; 
    683         } 
    684  
    685         //TODO need to make sure covar gets inverted.  Previously the Gauss 
    686         // Jordan solution would invert covar as a side effect. 
    687  
     869 
     870//TODO ARG covar needs to get inverted; previously inverted as a side effect 
    688871        /* Once converged, evaluate covariance matrix */ 
    689872        if (*alambda == 0) { 
     
    692875                                                                        instr, ninstr, noise, sig, 
    693876                                                                        param, paramfree, nparam, fitfunc, 
    694                                                                         yfit, dy, chisq,         
     877                                                                        yfit, dy, chisq, 
    695878                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv, 
    696                                                                         pfnvals_len, pdy_dparam_nparam_size) != 0) { 
     879                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0) 
    697880                        return -3; 
    698                 } 
    699881                if (*pmfit < nparam) {  /* no need to do this otherwise */ 
    700882                        GCI_covar_sort(covar, nparam, paramfree, *pmfit); 
     
    704886        } 
    705887 
    706         /* Did the trial succeed? * 
    707         int allzeroes = 1; 
     888//TODO ARG c/b special case if dparam all zeroes; don't have to recalc alpha & beta 
     889        /* Did the trial succeed? */ 
    708890        for (j=0, l=0; l<nparam; l++) 
    709                 if (paramfree[l]) { 
    710 //TODO I don't think it ever recovers if all dparams are zero; alambda grows and grows, goes from int to inf!! 
    711 // could be an early way to detect we're done. 
    712 // at least, if delta params is all zeroes, there is no point in calculating a new alpha and beta. 
    713 //                    printf("paramtry[%d] was %f becomes %f, param[%d] or %f plus dparam[%d] or %f\n", l, paramtry[l], param[l] + dparam[j], l, param[l], j, dparam[j]); 
    714                     if (dparam[j] != 0.0) allzeroes = 0; 
     891                if (paramfree[l]) 
    715892                        paramtry[l] = param[l] + dparam[j++]; 
    716 } 
    717      //TODO experimental: 
    718      //  gave about a 10% speedup 
    719         if (allzeroes) return -12345; */ 
    720893 
    721894        if (restrain == ECF_RESTRAIN_DEFAULT) 
     
    729902                return 0; 
    730903        } 
     904 
    731905        if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end, 
    732906                                                                           instr, ninstr, noise, sig, 
     
    735909                                                                           chisq, *pochisq, *alambda, 
    736910                                                                           pfnvals, pdy_dparam_pure, pdy_dparam_conv, 
    737                                                                            pfnvals_len, pdy_dparam_nparam_size) != 0) { 
     911                                                                           pfnvals_len, pdy_dparam_nparam_size) != 0) 
    738912                return -2; 
    739         } 
    740913 
    741914        /* Success, accept the new solution */ 
    742915        if (*chisq < *pochisq) { 
    743 //            printf("success, alambda goes from %f to %f\n", *alambda, (*alambda)*0.1); 
    744916                *alambda *= 0.1; 
    745917                *pochisq = *chisq; 
     
    752924                        param[l] = paramtry[l]; 
    753925        } else { /* Failure, increase alambda and return */ 
    754 //TODO if alambda goes to inf, can we recover? 
    755 //            printf("failure, alambda goes from %f to %f\n", *alambda, (*alambda)*10.0); 
    756926                *alambda *= 10.0; 
    757927                *chisq = *pochisq; 
     
    789959   determined the fitted convolved response. 
    790960 
    791    This is the variant which handles an instrument response. */ 
    792  
     961   Again there are two variants of this function, corresponding to the 
     962   two variants of the Marquardt function. 
     963*/ 
     964 
     965//TODO ARG deleted in my version; needs to get de-NRed!!!! 
     966int GCI_marquardt_compute_fn(float x[], float y[], int ndata, 
     967                                         noise_type noise, float sig[], 
     968                                         float param[], int paramfree[], int nparam, 
     969                                         void (*fitfunc)(float, float [], float *, float [], int), 
     970                                         float yfit[], float dy[], 
     971                                         float **alpha, float beta[], float *chisq, float old_chisq, 
     972                                         float alambda) 
     973{ 
     974        int i, j, k, l, m, mfit; 
     975        float wt, sig2i, y_ymod, dy_dparam[MAXFIT]; 
     976 
     977        for (j=0, mfit=0; j<nparam; j++) 
     978                if (paramfree[j]) 
     979                        mfit++; 
     980 
     981        //TODO ARG having this section { } of code here is just temporary; 
     982        // o'wise have to define these vars at top of function 
     983        { 
     984        float alpha_weight[MAXBINS]; 
     985        float beta_weight[MAXBINS]; 
     986        int q; 
     987        float weight; 
     988 
     989        int i_free; 
     990        int j_free; 
     991        float dot_product; 
     992        float beta_sum; 
     993        float dy_dparam_k_i; 
     994 
     995        *chisq = 0.0f; 
     996 
     997        switch (noise) { 
     998            case NOISE_CONST: 
     999            { 
     1000                for (q = 0; q < ndata; ++q) { 
     1001                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 
     1002                    yfit[q] += param[0]; 
     1003                    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 
     1006                    weight *= dy[q]; 
     1007                    beta_weight[q] = weight; // dy[q] 
     1008                    weight *= dy[q]; 
     1009                    *chisq += weight; // dy[q] * dy[q] 
     1010                } 
     1011                break; 
     1012            } 
     1013            case NOISE_GIVEN: 
     1014            { 
     1015                for (q = 0; q < ndata; ++q) { 
     1016                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 
     1017                    yfit[q] += param[0]; 
     1018                    dy[q] = y[q] - yfit[q]; 
     1019                    weight = 1.0f / (sig[q] * sig[q]); 
     1020                    alpha_weight[q] = weight; // 1 / (sig[q] * sig[q]) 
     1021                    weight *= dy[q]; 
     1022                    beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q]) 
     1023                    weight *= dy[q]; 
     1024                    *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q]) 
     1025                } 
     1026                break; 
     1027            } 
     1028            case NOISE_POISSON_DATA: 
     1029            { 
     1030                for (q = 0; q < ndata; ++q) { 
     1031                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 
     1032                    yfit[q] += param[0]; 
     1033                    dy[q] = y[q] - yfit[q]; 
     1034                    weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 
     1035                    alpha_weight[q] = weight; // 1 / sig(q) 
     1036                    weight *= dy[q]; 
     1037                    beta_weight[q] = weight; // dy[q] / sig(q) 
     1038                    weight *= dy[q]; 
     1039                    *chisq += weight; // (dy[q] * dy[q]) / sig(q) 
     1040                } 
     1041                break; 
     1042            } 
     1043            case NOISE_POISSON_FIT: 
     1044            { 
     1045                for (q = 0; q < ndata; ++q) { 
     1046                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 
     1047                    yfit[q] += param[0]; 
     1048                    dy[q] = y[q] - yfit[q]; 
     1049                    weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 
     1050                    alpha_weight[q] = weight; // 1 / sig(q) 
     1051                    weight *= dy[q]; 
     1052                    beta_weight[q] = weight; // dy(q) / sig(q) 
     1053                    weight *= dy[q]; 
     1054                    *chisq += weight; // (dy(q) * dy(q)) / sig(q) 
     1055                } 
     1056                break; 
     1057            } 
     1058            case NOISE_GAUSSIAN_FIT: 
     1059            { 
     1060                 for (q = 0; q < ndata; ++q) { 
     1061                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 
     1062                    yfit[q] += param[0]; 
     1063                    dy[q] = y[q] - yfit[q]; 
     1064                    weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 
     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; 
     1070                 } 
     1071                 break; 
     1072           } 
     1073            case NOISE_MLE: 
     1074            { 
     1075                for (q = 0; q < ndata; ++q) { 
     1076                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 
     1077                    yfit[q] += param[0]; 
     1078                    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? 
     1081                    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]); 
     1085                } 
     1086                if (*chisq <= 0.0f) { 
     1087                    *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve 
     1088                } 
     1089                break; 
     1090            } 
     1091            default: 
     1092            { 
     1093                return -3; 
     1094            } 
     1095        } 
     1096 
     1097        // Check if chi square worsened: 
     1098        if (0.0f != old_chisq && *chisq >= old_chisq) { 
     1099            // don't bother to set up the matrices for solution 
     1100            return 0; 
     1101        } 
     1102 
     1103        i_free = 0; 
     1104        // for all columns 
     1105        for (i = 0; i < nparam; ++i) { 
     1106            if (paramfree[i]) { 
     1107                j_free = 0; 
     1108                beta_sum = 0.0f; 
     1109                // row loop, only need to consider lower triangle 
     1110                for (j = 0; j <= i; ++j) { 
     1111                    if (paramfree[j]) { 
     1112                        dot_product = 0.0f; 
     1113                        if (0 == j_free) { 
     1114                            // for all data 
     1115                            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]; 
     1120                            } 
     1121                        } 
     1122                        else { 
     1123                            // for all data 
     1124                            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]; 
     1126                            } 
     1127                        } // k loop 
     1128 
     1129                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 
     1130                       // if (i_free != j_free) { 
     1131                       //     // matrix is symmetric 
     1132                       //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 
     1133                       // } 
     1134                        ++j_free; 
     1135                    } 
     1136                } // j loop 
     1137                beta[i_free] = beta_sum; 
     1138                ++i_free; 
     1139            } 
     1140            //else printf("param not free %d\n", i); 
     1141        } // 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]; 
     1333 
     1334        return 0; 
     1335} 
     1336 
     1337 
     1338/* And this is the variant which handles an instrument response. */ 
    7931339/* We assume that the function values are sensible. */ 
    7941340 
     
    8551401                if (paramfree[j]) mfit++; 
    8561402 
     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 
    8571411        /* Calculation of the fitting data will depend upon the type of 
    8581412           noise and the type of instrument response */ 
     
    9311485        /* OK, now we've got our (possibly convolved) data, we can do the 
    9321486           rest almost exactly as above. */ 
     1487        //TODO ARG this new section of code is just temporary; o'wise have to define all these variables at the top of the function 
    9331488        { 
    934         float alpha_weight[256]; //TODO establish maximum # bins and use elsewhere (#define) 
    935         float beta_weight[256]; 
     1489        float alpha_weight[MAXBINS]; 
     1490        float beta_weight[MAXBINS]; 
    9361491        int q; 
    9371492        float weight; 
     
    9421497        float beta_sum; 
    9431498        float dy_dparam_k_i; 
    944         float *alpha_weight_ptr; 
    945         float *beta_weight_ptr; 
    9461499 
    9471500        *chisq = 0.0f; 
     
    9501503            case NOISE_CONST: 
    9511504            { 
    952                 float *alpha_ptr = &alpha_weight[fit_start]; 
    953                 float *beta_ptr = &beta_weight[fit_start]; 
    9541505                for (q = fit_start; q < fit_end; ++q) { 
    9551506                    (*pdy_dparam_conv)[q][0] = 1.0f; 
    9561507                    yfit[q] += param[0]; 
    9571508                    dy[q] = y[q] - yfit[q]; 
    958                     weight = 1.0f; 
    959                     *alpha_ptr++ = weight; // 
    960                     //alpha_weight[q] = weight; // 1 
     1509                    weight = 1.0f; //TODO ARG this should be 1.0f / sig[0] ! 
     1510                    alpha_weight[q] = weight; // 1 
    9611511                    weight *= dy[q]; 
    962                     //printf("beta weight %d %f is y %f - yfit %f\n", q, weight, y[q], yfit[q]); 
    963                     *beta_ptr++ = weight; // 
    964                     //beta_weight[q] = weight; // dy[q] 
     1512                    beta_weight[q] = weight; // dy[q] 
    9651513                    weight *= dy[q]; 
    9661514                    *chisq += weight; // dy[q] * dy[q] 
     
    10681616                    if (paramfree[j]) { 
    10691617                        dot_product = 0.0f; 
    1070                               alpha_weight_ptr = &alpha_weight[fit_start]; 
    10711618                        if (0 == j_free) { 
    1072                             beta_weight_ptr = &beta_weight[fit_start]; 
    10731619                            // for all data 
    10741620                            for (k = fit_start; k < fit_end; ++k) { 
    10751621                                dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 
    1076                                 dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * (*alpha_weight_ptr++); //alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it. 
    1077                                 beta_sum += dy_dparam_k_i * (*beta_weight_ptr++); ///beta_weight[k]; 
     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. 
     1623                                beta_sum += dy_dparam_k_i * beta_weight[k]; 
    10781624                            } 
    10791625                        } 
     
    10811627                            // for all data 
    10821628                            for (k = fit_start; k < fit_end; ++k) { 
    1083                                 dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * (*alpha_weight_ptr++); // alpha_weight[k]; 
     1629                                dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; 
    10841630                            } 
    10851631                        } // k loop 
     
    10981644            //else printf("param not free %d\n", i); 
    10991645        } // i loop 
    1100  
    1101  /*       printf("i_free is %d j_free %d\n", i_free, j_free); 
    1102         //TODO try padding with zeroes; leftovers would affect AX=B solution 
    1103         for (i = i_free; i < nparam; ++i) { 
    1104             for (j = 0; j < nparam; ++j) { 
    1105                 alpha[i][j] = alpha[j][i] = 0.0f; 
    1106             } 
    1107             beta[i] = 0.0f; 
    1108         }*/ 
    1109  
    1110         } 
     1646        } 
     1647 
    11111648        return 0; 
    11121649} 
     
    11201657   noise models.  They do not calculate alpha or beta. */ 
    11211658 
    1122 /* This is the variant which handles an instrument response. */ 
     1659int GCI_marquardt_compute_fn_final(float x[], float y[], int ndata, 
     1660                                         noise_type noise, float sig[], 
     1661                                         float param[], int paramfree[], int nparam, 
     1662                                         void (*fitfunc)(float, float [], float *, float [], int), 
     1663                                         float yfit[], float dy[], float *chisq) 
     1664{ 
     1665        int i, j, mfit; 
     1666        float sig2i, dy_dparam[MAXFIT];  /* dy_dparam needed for fitfunc */ 
     1667 
     1668        for (j=0, mfit=0; j<nparam; j++) 
     1669                if (paramfree[j]) 
     1670                        mfit++; 
     1671 
     1672        /* Calculation of the fitting data will depend upon the type of 
     1673           noise and the type of instrument response */ 
     1674 
     1675        /* There's no convolution involved in this function.  This is then 
     1676           fairly straightforward, depending only upon the type of noise 
     1677           present.  Since we calculate the standard deviation at every 
     1678           point in a different way depending upon the type of noise, we 
     1679           will place our switch(noise) around the entire loop. */ 
     1680 
     1681        switch (noise) { 
     1682        case NOISE_CONST: 
     1683                *chisq = 0.0; 
     1684                /* Summation loop over all data */ 
     1685                for (i=0; i<ndata; i++) { 
     1686                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 
     1687                        yfit[i] += param[0]; 
     1688                        dy[i] = y[i] - yfit[i]; 
     1689                        /* And find chi^2 */ 
     1690                        *chisq += dy[i] * dy[i]; 
     1691                } 
     1692 
     1693                /* Now divide everything by sigma^2 */ 
     1694                sig2i = 1.0 / (sig[0] * sig[0]); 
     1695                *chisq *= sig2i; 
     1696                break; 
     1697 
     1698        case NOISE_GIVEN:  /* This is essentially the NR version */ 
     1699                *chisq = 0.0; 
     1700                /* Summation loop over all data */ 
     1701                for (i=0; i<ndata; i++) { 
     1702                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 
     1703                        yfit[i] += param[0]; 
     1704                        sig2i = 1.0 / (sig[i] * sig[i]); 
     1705                        dy[i] = y[i] - yfit[i]; 
     1706                        /* And find chi^2 */ 
     1707                        *chisq += dy[i] * dy[i] * sig2i; 
     1708                } 
     1709                break; 
     1710 
     1711        case NOISE_POISSON_DATA: 
     1712                *chisq = 0.0; 
     1713                /* Summation loop over all data */ 
     1714                for (i=0; i<ndata; i++) { 
     1715                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 
     1716                        yfit[i] += param[0]; 
     1717                        /* we still don't let the variance drop below 1 */ 
     1718                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0); 
     1719                        dy[i] = y[i] - yfit[i]; 
     1720                        /* And find chi^2 */ 
     1721                        *chisq += dy[i] * dy[i] * sig2i; 
     1722                } 
     1723                break; 
     1724 
     1725        case NOISE_POISSON_FIT: 
     1726                *chisq = 0.0; 
     1727                // Summation loop over all data 
     1728                for (i=0; i<ndata; i++) { 
     1729                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 
     1730                        yfit[i] += param[0]; 
     1731                        // we still don't let the variance drop below 1 
     1732                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 
     1733                        dy[i] = y[i] - yfit[i]; 
     1734                        // And find chi^2 
     1735                        *chisq += dy[i] * dy[i] * sig2i; 
     1736                } 
     1737                break; 
     1738 
     1739        case NOISE_MLE: 
     1740                        *chisq = 0.0; 
     1741                        /* Summation loop over all data */ 
     1742                        for (i=0; i<ndata; i++) { 
     1743                                (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 
     1744                                yfit[i] += param[0]; 
     1745//                              dy[i] = y[i] - yfit[i]; 
     1746 
     1747                                /* And find chi^2 */ 
     1748//                              sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 
     1749//                              *chisq += dy[i] * dy[i] * sig2i; 
     1750                                if (yfit[i]<=0.0) 
     1751                                        ; // do nothing 
     1752                                else if (y[i]==0.0) 
     1753                                        *chisq += 2.0*yfit[i];   // to avoid NaN from log 
     1754                                else 
     1755                                        *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; 
     1756                        } 
     1757                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve 
     1758                break; 
     1759 
     1760        case NOISE_GAUSSIAN_FIT: 
     1761                *chisq = 0.0; 
     1762                // Summation loop over all data 
     1763                for (i=0; i<ndata; i++) { 
     1764                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 
     1765                        yfit[i] += param[0]; 
     1766                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 
     1767                        dy[i] = y[i] - yfit[i]; 
     1768                        // And find chi^2 
     1769                        *chisq += dy[i] * dy[i] * sig2i; 
     1770                } 
     1771                break; 
     1772 
     1773        default: 
     1774                return -3; 
     1775                /* break; */   // (unreachable code) 
     1776        } 
     1777 
     1778        return 0; 
     1779} 
     1780 
     1781 
     1782/* And this is the variant which handles an instrument response. */ 
    11231783/* We assume that the function values are sensible.  Note also that we 
    11241784   have to be careful about which points which include in the 
     
    11321792                                   float param[], int paramfree[], int nparam, 
    11331793                                   void (*fitfunc)(float, float [], float *, float [], int), 
    1134                                    float yfit[], float dy[], float *chisq,       
     1794                                   float yfit[], float dy[], float *chisq, 
    11351795                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, 
    11361796                                        int *pfnvals_len, int *pdy_dparam_nparam_size) 
     
    11871847                for (i=0; i<ndata; i++) { 
    11881848                        int convpts; 
    1189                  
     1849 
    11901850                        /* We wish to find yfit = fnvals * instr, so explicitly: 
    11911851                             yfit[i] = sum_{j=0}^i fnvals[i-j].instr[j] 
     
    12221882                                                   dy_dparam_conv[i], nparam); 
    12231883        } 
    1224           
     1884 
    12251885        /* OK, now we've got our (possibly convolved) data, we can do the 
    12261886           rest almost exactly as above. */ 
     
    12931953        case NOISE_POISSON_FIT: 
    12941954                *chisq = 0.0; 
    1295                 // Summation loop over all data  
     1955                // Summation loop over all data 
    12961956                for (i=0; i<fit_start; i++) { 
    12971957                        yfit[i] += param[0]; 
     
    13011961                        yfit[i] += param[0]; 
    13021962                        dy[i] = y[i] - yfit[i]; 
    1303                         // And find chi^2  
    1304                         // we still don't let the variance drop below 1  
     1963                        // And find chi^2 
     1964                        // we still don't let the variance drop below 1 
    13051965                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 
    13061966                        *chisq += dy[i] * dy[i] * sig2i; 
     
    13111971                } 
    13121972                break; 
    1313                  
     1973 
    13141974        case NOISE_MLE:                              // for the final chisq report a normal chisq measure for MLE 
    13151975                *chisq = 0.0; 
    1316                 // Summation loop over all data  
     1976                // Summation loop over all data 
    13171977                for (i=0; i<fit_start; i++) { 
    13181978                        yfit[i] += param[0]; 
     
    13221982                        yfit[i] += param[0]; 
    13231983                        dy[i] = y[i] - yfit[i]; 
    1324                         // And find chi^2  
     1984                        // And find chi^2 
    13251985//                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 
    13261986                        if (yfit[i]<=0.0) 
     
    13402000        case NOISE_GAUSSIAN_FIT: 
    13412001                *chisq = 0.0; 
    1342                 // Summation loop over all data  
     2002                // Summation loop over all data 
    13432003                for (i=0; i<fit_start; i++) { 
    13442004                        yfit[i] += param[0]; 
     
    13482008                        yfit[i] += param[0]; 
    13492009                        dy[i] = y[i] - yfit[i]; 
    1350                         // And find chi^2  
     2010                        // And find chi^2 
    13512011                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 
    13522012                        *chisq += dy[i] * dy[i] * sig2i; 
     
    13572017                } 
    13582018                break; 
    1359          
     2019 
    13602020        default: 
    13612021                return -3; 
     
    13732033// was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03 
    13742034 
    1375 int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end,  
     2035int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end, 
    13762036                                                float prompt[], int nprompt, 
    13772037                                                noise_type noise, float sig[], 
     
    13882048        if (ecf_exportParams) ecf_ExportParams_OpenFile (); 
    13892049 
    1390         // All of the work is done by the ECF module  
     2050        // All of the work is done by the ECF module 
    13912051        ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end, 
    13922052                                                          prompt, nprompt, noise, sig, 
     
    13942054                                                          fitted, residuals, covar, alpha, &local_chisq, 
    13952055                                                          chisq_delta, chisq_percent, erraxes); 
    1396                                                            
     2056 
    13972057        // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same 
    13982058        // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr() 
    13992059        oldChisq = 3.0e38; 
    1400         while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<(MAXREFITS*3)) 
     2060        while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<MAXREFITS) 
    14012061        { 
    1402                 oldChisq = local_chisq;  
     2062                oldChisq = local_chisq; 
    14032063                tries++; 
    14042064                ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end, 
     
    14072067                                                          fitted, residuals, covar, alpha, &local_chisq, 
    14082068                                                          chisq_delta, chisq_percent, erraxes); 
    1409         }                                                  
    1410    /*     if (local_chisq<=chisq_target) { 
    1411             printf("local_chisq %f target %f\n", local_chisq, chisq_target); 
    1412         } 
    1413         if (local_chisq>=oldChisq) { 
    1414             printf("local_chisq %f oldChisq %f\n", local_chisq, oldChisq); 
    1415         } 
    1416         if (tries>=(MAXREFITS*3)) { 
    1417             printf("tries is %d\n", tries); 
    1418             printf("local_chisq %f oldChisq %f\n", local_chisq, oldChisq); 
    1419         } 
    1420  
    1421         show_timing();*/ 
     2069        } 
    14222070 
    14232071        if (chisq!=NULL) *chisq = local_chisq; 
     
    14412089} 
    14422090 
     2091 
    14432092// Emacs settings: 
    14442093// Local variables: 
  • trunk/projects/slim-curve/src/main/c/EcfUtil.c

    r7706 r7708  
    2222int ECF_debug = 0; 
    2323 
     24//TODO ARG 
     25//#define SPEEDUP1 1 
     26//#define SPEEDUP2 1 
     27//#define SPEEDUP3 1 
    2428 
    2529/******************************************************************** 
     
    176180    } 
    177181 
     182    #ifdef SPEEDUP1 
     183    //TODO ARG this is actually slower! 
     184    // swap rows 
     185    float *ptr1; 
     186    float *ptr2; 
     187    if (pivotRow != col) { 
     188        // swap elements in a matrix 
     189        ptr1 = &a[col][0]; 
     190        ptr2 = &a[pivotRow][0]; 
     191        for (i = 0; i < n; ++i) { 
     192            float temp; 
     193            SWAP(*ptr1, *ptr2); 
     194            ++ptr1; 
     195            ++ptr2; 
     196        } 
     197 
     198        // swap elements in order vector 
     199        { 
     200            int temp; 
     201            SWAP(order[col], order[pivotRow]); 
     202        } 
     203    } 
     204    #else 
    178205    // swap rows 
    179206    if (pivotRow != col) { 
     
    190217        } 
    191218    } 
     219    #endif 
    192220} 
    193221 
     
    210238    int kCol; 
    211239    float sum; 
    212      
     240 
    213241    // initialize ordering vector 
     242    #ifdef SPEEDUP2 
     243    //TODO ARG this is *slightly* slower 
     244    int *order_ptr = order; 
     245    for (i = 0; i < n; ++i) 
     246    { 
     247        *order_ptr++ = i; 
     248    } 
     249    #else 
    214250    for (i = 0; i < n; ++i) 
    215251    { 
    216252        order[i] = i; 
    217253    } 
     254    #endif 
    218255 
    219256    // pivot first column 
     
    227264 
    228265    // compute first row of upper 
     266    #ifdef SPEEDUP3 
     267    //TODO ARG this is *much* slower!!! 
     268    //  Note compiler probably realizes a[0] is a constant anyway 
     269    inverse = 1.0 / a[0][0]; 
     270    float *a_0_ptr = a[0]; 
     271    for (i = 1; i < n; ++i) { 
     272        *a_0_ptr++ *= inverse; 
     273    } 
     274    #else 
    229275    inverse = 1.0 / a[0][0]; 
    230276    for (i = 1; i < n; ++i) { 
    231277        a[0][i] *= inverse; 
    232278    } 
     279    #endif 
    233280 
    234281    // continue computing columns of lowers then rows of uppers 
     
    377424int GCI_solve(float **a, int n, float *b) 
    378425{ 
    379     return GCI_solve_Gaussian(a, n, b); 
    380     //return GCI_solve_lu_decomp(a, n, b); 
     426    //return GCI_solve_Gaussian(a, n, b); 
     427    return GCI_solve_lu_decomp(a, n, b); 
    381428} 
    382429 
     
    387434int GCI_invert(float **a, int n) 
    388435{ 
    389     return GCI_invert_Gaussian(a, n); 
    390     //return GCI_invert_lu_decomp(a, n); 
     436    //return GCI_invert_Gaussian(a, n); 
     437    return GCI_invert_lu_decomp(a, n); 
    391438} 
    392439 
  • trunk/projects/slim-curve/src/main/c/ecf.h

    r7387 r7708  
    44#ifndef _GCI_ECF 
    55#define _GCI_ECF 
    6  
    7 #define ELIMINATE_NR 1 //TODO ARG 9/14/10 
    86 
    97/* #defines which are publically needed */ 
Note: See TracChangeset for help on using the changeset viewer.