source: trunk/projects/slim-curve/src/main/c/EcfSingle.c @ 7706

Revision 7706, 47.4 KB checked in by aivar, 9 years ago (diff)

Changed GCI_gauss_jordan to be GCI_solve and GCI_invert. The current default CGI_solve uses Gaussian Elimination. Added inversion method. Added Lower/Upper Decomposition code for both solve and invert, but it is commented out.

Line 
1/* Based on the 2003 version of the ECF library.  This has been
2   modified to remove modified Numeric Recipes code. Also, this
3   takes account of the fact that we may be handling Poisson noise.
4
5   This file contains functions for single transient analysis.
6   Utility code is found in EcfUtil.c.
7*/
8
9#include <math.h>
10#include <stdio.h>
11#include <stdlib.h>
12#include "EcfInternal.h"
13
14/* Predeclarations */
15int GCI_marquardt_step_instr(float xincr, float y[],
16                                        int ndata, int fit_start, int fit_end,
17                                        float instr[], int ninstr,
18                                        noise_type noise, float sig[],
19                                        float param[], int paramfree[], int nparam,
20                                        restrain_type restrain,
21                                        void (*fitfunc)(float, float [], float *, float [], int),
22                                        float yfit[], float dy[],
23                                        float **covar, float **alpha, float *chisq,
24                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,       
25                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
26                                        int *pfnvals_len, int *pdy_dparam_nparam_size);
27
28/* Globals */
29//static float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
30//static int fnvals_len=0, dy_dparam_nparam_size=0;
31// was Global, now thread safe
32
33/********************************************************************
34
35                                           SINGLE TRANSIENT FITTING
36
37                                                TRIPLE INTEGRAL METHOD
38
39 ********************************************************************/
40
41/* Start with an easy one: the three integral method.  This returns 0
42   on success, negative on error. */
43
44int GCI_triple_integral(float xincr, float y[],
45                                                int fit_start, int fit_end,
46                                                noise_type noise, float sig[],
47                                                float *Z, float *A, float *tau,
48                                                float *fitted, float *residuals,
49                                                float *chisq, int division)
50{
51        float d1, d2, d3, d12, d23;
52        float t0, dt, exp_dt_tau, exp_t0_tau;
53        int width;
54        int i;
55        float sigma2, res, chisq_local;
56
57        width = (fit_end - fit_start) / division;
58        if (width <= 0)
59                return -1;
60
61        t0 = fit_start * xincr;
62        dt = width * xincr;
63
64        d1 = d2 = d3 = 0;
65        for (i=fit_start; i<fit_start+width; i++)           d1 += y[i];
66        for (i=fit_start+width; i<fit_start+2*width; i++)   d2 += y[i];
67        for (i=fit_start+2*width; i<fit_start+3*width; i++) d3 += y[i];
68
69        // Those are raw sums, we now convert into areas */
70        d1 *= xincr;
71        d2 *= xincr;
72        d3 *= xincr;
73
74        d12 = d1 - d2;
75        d23 = d2 - d3;
76        if (d12 <= d23 || d23 <= 0)
77                return -2;
78
79        exp_dt_tau = d23 / d12;  /* exp(-dt/tau) */
80        *tau = -dt / log(exp_dt_tau);
81        exp_t0_tau = exp(-t0/(*tau));
82        *A = d12 / ((*tau) * exp_t0_tau * (1 - 2*exp_dt_tau + exp_dt_tau*exp_dt_tau));
83        *Z = (d1 - (*A) * (*tau) * exp_t0_tau * (1 - exp_dt_tau)) / dt;
84
85        /* Now calculate the fitted curve and chi-squared if wanted. */
86        if (fitted == NULL)
87                return 0;
88
89        for (i=0; i<fit_end; i++)
90                fitted[i] = (*Z) + (*A) * exp(-i*xincr/(*tau));
91
92        // OK, so now fitted contains our data for the timeslice of interest.
93        // We can calculate a chisq value and plot the graph, along with
94        // the residuals.
95
96        if (residuals == NULL && chisq == NULL)
97                return 0;
98
99        chisq_local = 0.0;
100        for (i=0; i<fit_start; i++) {
101                res = y[i]-fitted[i];
102                if (residuals != NULL)
103                        residuals[i] = res;
104        }
105
106        switch (noise) {
107        case NOISE_CONST:
108                /* Summation loop over all data */
109                for ( ; i<fit_end; i++) {
110                        res = y[i] - fitted[i];
111                        if (residuals != NULL)
112                                residuals[i] = res;
113                        chisq_local += res * res;
114                }
115
116                chisq_local /= (sig[0]*sig[0]);
117                break;
118
119        case NOISE_GIVEN:
120                /* Summation loop over all data */
121                for ( ; i<fit_end; i++) {
122                        res = y[i] - fitted[i];
123                        if (residuals != NULL)
124                                residuals[i] = res;
125                        chisq_local += (res * res) / (sig[i] * sig[i]);
126                }
127                break;
128
129        case NOISE_POISSON_DATA:
130                /* Summation loop over all data */
131                for ( ; i<fit_end; i++) {
132                        res = y[i] - fitted[i];
133                        if (residuals != NULL)
134                                residuals[i] = res;
135                        /* don't let variance drop below 1 */
136                        sigma2 = (y[i] > 1 ? 1.0/y[i] : 1.0);
137                        chisq_local += res * res * sigma2;
138                }
139                break;
140
141        case NOISE_POISSON_FIT:
142        case NOISE_GAUSSIAN_FIT:   //  NOISE_GAUSSIAN_FIT and NOISE_MLE not implemented for triple integral
143        case NOISE_MLE:
144                /* Summation loop over all data */
145                for ( ; i<fit_end; i++) {
146                        res = y[i] - fitted[i];
147                        if (residuals != NULL)
148                                residuals[i] = res;
149                        /* don't let variance drop below 1 */
150                        sigma2 = (fitted[i] > 1 ? 1.0/fitted[i] : 1.0);
151                        chisq_local += res * res * sigma2;
152                }
153                break;
154
155        default:
156                return -3;
157                /* break; */   // (unreachable code)
158        }
159
160        if (chisq != NULL)
161                *chisq = chisq_local;
162
163        return 0;
164}
165
166
167int GCI_triple_integral_instr(float xincr, float y[],
168                                                          int fit_start, int fit_end,
169                                                          float instr[], int ninstr,
170                                                          noise_type noise, float sig[],
171                                                          float *Z, float *A, float *tau,
172                                                          float *fitted, float *residuals,
173                                                          float *chisq, int division)
174{
175        float d1, d2, d3, d12, d23;
176        float t0, dt, exp_dt_tau, exp_t0_tau;
177        int width;
178        int i, j;
179        float sigma2, res, chisq_local;
180        float sum, scaling;
181        int fitted_preconv_size=0;   // was static but now thread safe
182        float *fitted_preconv;   // was static but now thread safe
183
184        width = (fit_end - fit_start) / division;
185        if (width <= 0)
186                return -1;
187
188        t0 = fit_start * xincr;
189        dt = width * xincr;
190
191        d1 = d2 = d3 = 0;
192        for (i=fit_start; i<fit_start+width; i++)
193                d1 += y[i];
194        for (i=fit_start+width; i<fit_start+2*width; i++)
195                d2 += y[i];
196        for (i=fit_start+2*width; i<fit_start+3*width; i++)
197                d3 += y[i];
198
199        // Those are raw sums, we now convert into areas */
200        d1 *= xincr;
201        d2 *= xincr;
202        d3 *= xincr;
203
204        d12 = d1 - d2;
205        d23 = d2 - d3;
206        if (d12 <= d23 || d23 <= 0)
207                return -2;
208
209        exp_dt_tau = d23 / d12;  /* exp(-dt/tau) */
210        *tau = -dt / log(exp_dt_tau);
211        exp_t0_tau = exp(-t0/(*tau));
212        *A = d12 /
213                  ((*tau) * exp_t0_tau * (1 - 2*exp_dt_tau + exp_dt_tau*exp_dt_tau));
214        *Z = (d1 - (*A) * (*tau) * exp_t0_tau * (1 - exp_dt_tau)) / dt;
215
216        /* We now convolve with the instrument response to hopefully get a
217           slightly better fit.  We'll also scale by an appropriate
218           scale factor, which turns out to be:
219
220                  sum_{i=0}^{ninstr-1} instr[i]*exp(i*xincr/tau)
221
222           which should be only a little greater than the sum of the
223           instrument response values.
224        */
225
226        sum = scaling = 0;
227        for (i=0; i<ninstr; i++) {
228                sum += instr[i];
229                scaling += instr[i] * exp(i*xincr/(*tau));
230        }
231
232        scaling /= sum;  /* Make instrument response sum to 1.0 */
233        (*A) /= scaling;
234
235        /* Now calculate the fitted curve and chi-squared if wanted. */
236        if (fitted == NULL)
237                return 0;
238
239//      if (fitted_preconv_size < fit_end) {
240//              if (fitted_preconv_size > 0)
241//                      free(fitted_preconv);
242                if ((fitted_preconv = (float *) malloc(fit_end * sizeof(float)))
243                        == NULL)
244                        return -3;
245                else
246                        fitted_preconv_size = fit_end;
247//      }
248
249        for (i=0; i<fit_end; i++)
250                fitted_preconv[i] = (*A) * exp(-i*xincr/(*tau));
251
252        for (i=0; i<fit_end; i++) {
253                int convpts;
254
255                /* (Zero-basing everything in sight...)
256                   We wish to find fitted = fitted_preconv * instr, so explicitly:
257                          fitted[i] = sum_{j=0}^i fitted_preconv[i-j].instr[j]
258                   But instr[k]=0 for k >= ninstr, so we only need to sum:
259                          fitted[i] = sum_{j=0}^{min(ninstr-1,i)}
260                                              fitted_preconv[i-j].instr[j]
261                */
262
263                fitted[i] = 0;
264                convpts = (ninstr <= i) ? ninstr-1 : i;
265                for (j=0; j<=convpts; j++) {
266                        fitted[i] += fitted_preconv[i-j]*instr[j];
267                }
268
269                fitted[i] += *Z;
270        }
271
272        // OK, so now fitted contains our data for the timeslice of interest.
273        // We can calculate a chisq value and plot the graph, along with
274        // the residuals.
275
276        if (residuals == NULL && chisq == NULL)
277                return 0;
278
279        chisq_local = 0.0;
280        for (i=0; i<fit_start; i++) {
281                res = y[i]-fitted[i];
282                if (residuals != NULL)
283                        residuals[i] = res;
284        }
285
286        switch (noise) {
287        case NOISE_CONST:
288                /* Summation loop over all data */
289                for ( ; i<fit_end; i++) {
290                        res = y[i] - fitted[i];
291                        if (residuals != NULL)
292                                residuals[i] = res;
293                        chisq_local += res * res;
294                }
295
296                chisq_local /= (sig[0]*sig[0]);
297                break;
298
299        case NOISE_GIVEN:
300                /* Summation loop over all data */
301                for ( ; i<fit_end; i++) {
302                        res = y[i] - fitted[i];
303                        if (residuals != NULL)
304                                residuals[i] = res;
305                        chisq_local += (res * res) / (sig[i] * sig[i]);
306                }
307                break;
308
309        case NOISE_POISSON_DATA:
310                /* Summation loop over all data */
311                for ( ; i<fit_end; i++) {
312                        res = y[i] - fitted[i];
313                        if (residuals != NULL)
314                                residuals[i] = res;
315                        /* don't let variance drop below 1 */
316                        sigma2 = (y[i] > 1 ? 1.0/y[i] : 1.0);
317                        chisq_local += res * res * sigma2;
318                }
319                break;
320
321        case NOISE_POISSON_FIT:
322        case NOISE_GAUSSIAN_FIT:   //  NOISE_GAUSSIAN_FIT and NOISE_MLE not implemented for triple integral
323        case NOISE_MLE:
324                /* Summation loop over all data */
325                for ( ; i<fit_end; i++) {
326                        res = y[i] - fitted[i];
327                        if (residuals != NULL)
328                                residuals[i] = res;
329                        /* don't let variance drop below 1 */
330                        sigma2 = (fitted[i] > 1 ? 1.0/fitted[i] : 1.0);
331                        chisq_local += res * res * sigma2;
332                }
333                break;
334
335        default:
336                return -3;
337                /* break; */   // (unreachable code)
338        }
339
340        if (chisq != NULL)
341                *chisq = chisq_local;
342
343        return 0;
344}
345
346int GCI_triple_integral_fitting_engine(float xincr, float y[], int fit_start, int fit_end,
347                                                          float instr[], int ninstr, noise_type noise, float sig[],
348                                                          float *Z, float *A, float *tau, float *fitted, float *residuals,
349                                                          float *chisq, float chisq_target)
350{
351        int tries=1, division=3;                 // the data
352        float local_chisq=3.0e38, oldChisq=3.0e38, oldZ, oldA, oldTau, *validFittedArray; // local_chisq a very high float but below oldChisq
353               
354        if (fitted==NULL)   // we require chisq but have not supplied a "fitted" array so must malloc one
355        {
356                if ((validFittedArray = malloc(fit_end * sizeof(float)))== NULL) return (-1);
357        }
358        else validFittedArray = fitted;
359       
360        if (instr==NULL)           // no instrument/prompt has been supplied
361        {
362                GCI_triple_integral(xincr, y, fit_start, fit_end, noise, sig,
363                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
364
365                while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS)   
366                {
367                        oldChisq = local_chisq;
368                        oldZ = *Z;
369                        oldA = *A;
370                        oldTau = *tau;
371//                      division++;                     
372                        division+=division/3;                   
373                        tries++;
374                        GCI_triple_integral(xincr, y, fit_start, fit_end, noise, sig,
375                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
376                }                                                 
377        }
378        else
379        {
380                GCI_triple_integral_instr(xincr, y, fit_start, fit_end, instr, ninstr, noise, sig,
381                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
382
383                while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS)   
384                {
385                        oldChisq = local_chisq; 
386                        oldZ = *Z;
387                        oldA = *A;
388                        oldTau = *tau;
389//                      division++;                     
390                        division+=division/3;
391                        tries++;
392                        GCI_triple_integral_instr(xincr, y, fit_start, fit_end, instr, ninstr, noise, sig,
393                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
394
395                }                                                 
396        }
397
398        if (local_chisq>oldChisq)          // the previous fit was better
399        {
400                local_chisq = oldChisq;
401                *Z = oldZ;
402                *A = oldA;
403                *tau = oldTau;
404        }
405       
406        if (chisq!=NULL) *chisq = local_chisq;
407
408        if (fitted==NULL) 
409        {
410                free (validFittedArray);
411        }
412
413        return(tries);
414}
415
416/********************************************************************
417
418                                           SINGLE TRANSIENT FITTING
419
420                                          LEVENBERG-MARQUARDT METHOD
421
422 ********************************************************************/
423
424/* Now for the non-linear least squares fitting algorithms.
425
426   The process is:
427   - for Gaussian noise, use Levenberg-Marquardt directly
428   - for Poisson noise, use Levenberg-Marquardt to get an initial
429   estimate of the parameters assuming constant error variance, then
430   use amoeba to improve the estimate, assuming that the error
431   variance is proportional to the function value (with a minimum
432   variance of 15 to handle the case when the Poisson distribution
433   is not approximately Gaussian, so that the very noisy tails do
434   not inappropriately weight the solution).
435
436
437   This code contains two variants of the Levenberg-Marquardt method
438   for slightly different situations.  This attempts to reduce the
439   value chi^2 of a fit between a set of data points x[0..ndata-1],
440   y[0..ndata-1] and a nonlinear function dependent on nparam
441   coefficients param[0..nparam-1].  In the case that the x values are
442   equally spaced and start at zero, we can also handle convolution
443   with an instrument response instr[0..ninstr-1] and only look at the
444   data points from fit_start..fit_end-1.  The first variant does not
445   handle an instrument response and takes any values of
446   x[0..ndata-1].  The second variant takes an xincr and will handle
447   an instrument response if ninstr > 0.  The individual standard
448   deviations of the errors are determined by the value of noise: if
449   noise=NOISE_CONST, the standard deviations are constant, given by
450   sig[0]=*sig, if noise=NOISE_GIVEN, the standard deviations are
451   given by sig[0..ndata-1], if noise=NOISE_POISSON_DATA, the standard
452   deviations are taken to be given by Poisson noise, and the
453   variances are taken to be max(y[i],15), and if
454   noise=NOISE_POISSON_FIT, the variances are taken to be
455   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
458   G. Nishimura, and M. Tamura Phys Med Biol 2005).
459   
460   The input array paramfree[0..nparam-1] indicates
461   by nonzero entries those components that should be held fixed at
462   their input values. The program returns current best-fit values for
463   the parameters param[0..nparam-1] and chi^2 = chisq.  The arrays
464   covar[0..nparam-1][0..nparam-1] and alpha[0..nparam-1][0..nparam-1]
465   are used as working space during most isterations.  Supply a
466   routine fitfunc(x,param,yfit,dy_dparam,nparam) that evaluates the
467   fitting function fitfunc and its derivatives dydy[1..nparam-1] with
468   respect to the fitting parameters param at x.  (See below for
469   information about zero offsets, though.)  The values of fitfunc,
470   modified by the instrument response, are returned in the array yfit
471   and the differences y - yfit in dy.  The first call _must_ provide
472   an initial guess for the parameters param and set alambda < 0 for
473   initialisation (which then sets alambda = 0.001).  If a step
474   succeeds, chisq becomes smaller and alambda decreases by a factor
475   of 10.  You must call this routine repeatedly until convergence is
476   achieved. Then make one final call with alambda=0 to perform
477   cleanups and so that covar[0..nparam-1][0..nparam-1] returns the
478   covariance matrix and alpha the curvature matrix.  (Parameters held
479   fixed will return zero covariances.)
480
481   One key extra piece which is particularly important in the
482   instrument response case.  The parameter param[0] is assumed to be
483   the zero offset of the signal, which applies before and after time
484   zero.  It thus simply contributes param[0]*sum(instr) to the signal
485   value rather than being convolved with the instrument response only
486   from time zero.  For this reason, the fitfunc should ignore
487   param[0], as the fitting routines will handle this offset.
488*/
489
490/* This functions does the whole job */
491
492int GCI_marquardt(float x[], float y[], int ndata,
493                                  noise_type noise, float sig[],
494                                  float param[], int paramfree[], int nparam,
495                                  restrain_type restrain,
496                                  void (*fitfunc)(float, float [], float *, float [], int),
497                                  float *fitted, float *residuals,
498                                  float **covar, float **alpha, float *chisq,
499                                  float chisq_delta, float chisq_percent, float **erraxes)
500{
501        // Empty fn to allow compile in TRI2
502}
503
504#define do_frees \
505        if (fnvals) free(fnvals);\
506        if (dy_dparam_pure) GCI_ecf_free_matrix(dy_dparam_pure);\
507        if (dy_dparam_conv) GCI_ecf_free_matrix(dy_dparam_conv);
508
509int GCI_marquardt_instr(float xincr, float y[],
510                                        int ndata, int fit_start, int fit_end,
511                                        float instr[], int ninstr,
512                                        noise_type noise, float sig[],
513                                        float param[], int paramfree[], int nparam,
514                                        restrain_type restrain,
515                                        void (*fitfunc)(float, float [], float *, float [], int),
516                                        float *fitted, float *residuals,
517                                        float **covar, float **alpha, float *chisq,
518                                        float chisq_delta, float chisq_percent, float **erraxes)
519{
520        float alambda, ochisq;
521        int mfit, mfit2;
522        float evals[MAXFIT];
523        int i, k, itst, itst_max;
524
525        // The following are declared here to retain some optimisation by not repeatedly mallocing
526        // (only once per transient), but to remain thread safe.
527        // They are malloced by lower fns but at the end, freed by this fn.
528        // These vars were global or static before thread safety was introduced.
529        float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL;
530        int fnvals_len=0, dy_dparam_nparam_size=0;
531        float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 
532
533        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
534
535        mfit = 0;
536        for (i=0; i<nparam; i++) {
537                if (paramfree[i]) mfit++;
538        }
539
540        if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
541
542        alambda = -1;
543        if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
544                                                                 instr, ninstr, noise, sig,
545                                                                 param, paramfree, nparam, restrain,
546                                                                 fitfunc, fitted, residuals,
547                                                                 covar, alpha, chisq, &alambda, 
548                                                                 &mfit2, &ochisq2, paramtry, beta, dparam,
549                                                                 &fnvals, &dy_dparam_pure, &dy_dparam_conv,
550                                                                 &fnvals_len, &dy_dparam_nparam_size) != 0) {
551                do_frees
552                return -1;
553        }
554
555        if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
556
557        k = 1;  /* Iteration counter */
558        itst = 0;
559        for (;;) {
560                k++;
561                if (k > MAXITERS) {
562                        do_frees
563                        return -2;
564                }
565
566                ochisq = *chisq;
567                if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
568                                                                         instr, ninstr, noise, sig,
569                                                                         param, paramfree, nparam, restrain,
570                                                                         fitfunc, fitted, residuals,
571                                                                         covar, alpha, chisq, &alambda, 
572                                                                         &mfit2, &ochisq2, paramtry, beta, dparam,
573                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv,
574                                                                         &fnvals_len, &dy_dparam_nparam_size) != 0) {
575                        do_frees
576                        return -3;
577                }
578
579                if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
580
581                if (*chisq > ochisq)
582                        itst = 0;
583                else if (ochisq - *chisq < chisq_delta)
584                        itst++;
585
586                if (itst < itst_max) continue;
587
588                /* Endgame */
589                alambda = 0.0;
590                if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
591                                                                         instr, ninstr, noise, sig,
592                                                                         param, paramfree, nparam, restrain,
593                                                                         fitfunc, fitted, residuals,
594                                                                         covar, alpha, chisq, &alambda, 
595                                                                         &mfit2, &ochisq2, paramtry, beta, dparam,
596                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv,
597                                                                         &fnvals_len, &dy_dparam_nparam_size) != 0) {
598                        do_frees
599                        return -4;
600                }
601
602                if (erraxes == NULL){
603                        do_frees
604                        return k;
605                }
606
607                break;  /* We're done now */
608        }
609
610        do_frees
611        return k;
612}
613#undef do_frees
614
615
616int GCI_marquardt_step_instr(float xincr, float y[],
617                                        int ndata, int fit_start, int fit_end,
618                                        float instr[], int ninstr,
619                                        noise_type noise, float sig[],
620                                        float param[], int paramfree[], int nparam,
621                                        restrain_type restrain,
622                                        void (*fitfunc)(float, float [], float *, float [], int),
623                                        float yfit[], float dy[],
624                                        float **covar, float **alpha, float *chisq,
625                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,       
626                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
627                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
628{
629        int j, k, l, ret;
630//      static int mfit;   // was static but now thread safe
631//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe
632
633        /*if (nparam > MAXFIT) {
634            printf("GCI_marq_step_instr returns -10\n");
635                return -10;
636        }
637        if (xincr <= 0) {
638            printf("GCI_marq_step_instr returns -11\n");
639                return -11;
640        }
641        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) {
642            printf("GCI_marq_step_instr returns -12\n");
643                return -12;
644        }*/
645
646        /* Initialisation */
647        /* We assume we're given sensible starting values for param[] */
648        if (*alambda < 0.0) {
649
650                for ((*pmfit)=0, j=0; j<nparam; j++)
651                        if (paramfree[j])
652                                (*pmfit)++;
653
654                if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
655                                                                                   instr, ninstr, noise, sig,
656                                                                                   param, paramfree, nparam, fitfunc,
657                                                                                   yfit, dy, alpha, beta, chisq, 0.0,
658                                                                                   *alambda,   
659                                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
660                                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0) {
661                        return -2;
662                }
663
664                *alambda = 0.001;
665                *pochisq = *chisq;
666                for (j=0; j<nparam; j++)
667                        paramtry[j] = param[j];
668
669        }
670
671        /* Alter linearised fitting matrix by augmenting diagonal elements */
672        for (j=0; j<(*pmfit); j++) {
673                for (k=0; k<(*pmfit); k++)
674                        covar[j][k] = alpha[j][k];
675
676                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
677                dparam[j] = beta[j];
678        }
679
680        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */
681        if (GCI_solve(covar, *pmfit, dparam) != 0) {
682                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
688        /* Once converged, evaluate covariance matrix */
689        if (*alambda == 0) {
690                if (GCI_marquardt_compute_fn_final_instr(
691                                                                        xincr, y, ndata, fit_start, fit_end,
692                                                                        instr, ninstr, noise, sig,
693                                                                        param, paramfree, nparam, fitfunc,
694                                                                        yfit, dy, chisq,       
695                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
696                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0) {
697                        return -3;
698                }
699                if (*pmfit < nparam) {  /* no need to do this otherwise */
700                        GCI_covar_sort(covar, nparam, paramfree, *pmfit);
701                        GCI_covar_sort(alpha, nparam, paramfree, *pmfit);
702                }
703                return 0;
704        }
705
706        /* Did the trial succeed? *
707        int allzeroes = 1;
708        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;
715                        paramtry[l] = param[l] + dparam[j++];
716}
717     //TODO experimental:
718     //  gave about a 10% speedup
719        if (allzeroes) return -12345; */
720
721        if (restrain == ECF_RESTRAIN_DEFAULT)
722                ret = check_ecf_params (paramtry, nparam, fitfunc);
723        else
724                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
725
726        if (ret != 0) {
727                /* Bad parameters, increase alambda and return */
728                *alambda *= 10.0;
729                return 0;
730        }
731        if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
732                                                                           instr, ninstr, noise, sig,
733                                                                           paramtry, paramfree, nparam, fitfunc,
734                                                                           yfit, dy, covar, dparam,
735                                                                           chisq, *pochisq, *alambda,
736                                                                           pfnvals, pdy_dparam_pure, pdy_dparam_conv,
737                                                                           pfnvals_len, pdy_dparam_nparam_size) != 0) {
738                return -2;
739        }
740
741        /* Success, accept the new solution */
742        if (*chisq < *pochisq) {
743//            printf("success, alambda goes from %f to %f\n", *alambda, (*alambda)*0.1);
744                *alambda *= 0.1;
745                *pochisq = *chisq;
746                for (j=0; j<(*pmfit); j++) {
747                        for (k=0; k<(*pmfit); k++)
748                                alpha[j][k] = covar[j][k];
749                        beta[j] = dparam[j];
750                }
751                for (l=0; l<nparam; l++)
752                        param[l] = paramtry[l];
753        } 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);
756                *alambda *= 10.0;
757                *chisq = *pochisq;
758        }
759
760        return 0;
761}
762
763
764/* Used by GCI_marquardt to evaluate the linearised fitting matrix alpha
765   and vector beta and to calculate chi^2.      The equations involved are
766   given in section 15.5 of Numerical Recipes; basically:
767
768     \chi^2(param) = \sum_{i=1}^N ( (y_i-y(x_i;param)) / sigma_i )^2
769
770     beta_k = -1/2 (d/dparam_k)(chi^2)
771
772     alpha_kl = \sum_{i=1}^N (1/sigma_i^2) .
773                  (dy(x_i;param)/dparam_k) . (dy(x_i;param)/dparam_l)
774
775   (where all of the derivatives are partial).
776
777   If an instrument response is provided, we also take account of it
778   now.  We are given that:
779
780     observed(t) = response(t) * instr(t)
781
782   where response(t) is being fitted with fitfunc; it is also trivial to
783   show that (assuming that instr(t) is known and fixed, with no
784   dependencies on the param_k, the parameters being fitted):
785
786     (d/dparam_k) observed(t) = ((d/dparam_k) response(t)) * instr(t)
787
788   so we do not need to alter the response function in any way to
789   determined the fitted convolved response.
790
791   This is the variant which handles an instrument response. */
792
793/* We assume that the function values are sensible. */
794
795int GCI_marquardt_compute_fn_instr(float xincr, float y[], int ndata,
796                                   int fit_start, int fit_end,
797                                   float instr[], int ninstr,
798                                   noise_type noise, float sig[],
799                                   float param[], int paramfree[], int nparam,
800                                   void (*fitfunc)(float, float [], float *, float [], int),
801                                   float yfit[], float dy[],
802                                   float **alpha, float beta[], float *chisq, float old_chisq,
803                                   float alambda,
804                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
805                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
806{
807        int i, j, k, l, m, mfit, ret;
808        float wt, sig2i, y_ymod;
809
810        /* Are we initialising? */
811        // Malloc the arrays that will get used again in this fit via the pointers passed in
812        // They will be freed by the higher fn that declared them.
813        if (alambda < 0) {
814                /* do any necessary initialisation */
815                if ((*pfnvals_len) < ndata) {  /* we will need this length for
816                                                                          the final full computation */
817                        if ((*pfnvals_len)) {
818                                free((*pfnvals));
819                                GCI_ecf_free_matrix((*pdy_dparam_pure));
820                                GCI_ecf_free_matrix((*pdy_dparam_conv));
821                                (*pfnvals_len) = 0;
822                                (*pdy_dparam_nparam_size) = 0;
823                        }
824                } else if ((*pdy_dparam_nparam_size) < nparam) {
825                        GCI_ecf_free_matrix((*pdy_dparam_pure));
826                        GCI_ecf_free_matrix((*pdy_dparam_conv));
827                        (*pdy_dparam_nparam_size) = 0;
828                }
829                if (! (*pfnvals_len)) {
830                        if (((*pfnvals) = (float *) malloc(ndata * sizeof(float)))
831                                == NULL)
832                                return -1;
833                        (*pfnvals_len) = ndata;
834                }
835                if (! (*pdy_dparam_nparam_size)) {
836                        /* use (*pfnvals_len), not ndata here! */
837                        if (((*pdy_dparam_pure) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
838                                /* Can't keep (*pfnvals) around in this case */
839                                free((*pfnvals));
840                                (*pfnvals_len) = 0;
841                                return -1;
842                        }
843                        if (((*pdy_dparam_conv) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
844                                /* Can't keep (*pfnvals) around in this case */
845                                free((*pfnvals));
846                                free((*pdy_dparam_pure));
847                                (*pfnvals_len) = 0;
848                                return -1;
849                        }
850                        (*pdy_dparam_nparam_size) = nparam;
851                }
852        }
853
854        for (j=0, mfit=0; j<nparam; j++)
855                if (paramfree[j]) mfit++;
856
857        /* Calculation of the fitting data will depend upon the type of
858           noise and the type of instrument response */
859
860        /* Need to calculate unconvolved values all the way down to 0 for
861           the instrument response case */
862        if (ninstr > 0) {
863                if (fitfunc == GCI_multiexp_lambda)
864                        ret = multiexp_lambda_array(xincr, param, (*pfnvals),
865                                                                                (*pdy_dparam_pure), fit_end, nparam);
866                else if (fitfunc == GCI_multiexp_tau)
867                        ret = multiexp_tau_array(xincr, param, (*pfnvals),
868                                                                         (*pdy_dparam_pure), fit_end, nparam);
869                else if (fitfunc == GCI_stretchedexp)
870                        ret = stretchedexp_array(xincr, param, (*pfnvals),
871                                                                         (*pdy_dparam_pure), fit_end, nparam);
872                else
873                        ret = -1;
874
875                if (ret < 0)
876                        for (i=0; i<fit_end; i++)
877                                (*fitfunc)(xincr*i, param, &(*pfnvals)[i],
878                                                   (*pdy_dparam_pure)[i], nparam);
879
880                /* OK, we've got to convolve the model fit with the given
881                   instrument response.  What we'll do here, then, is to
882                   generate the whole model fit first, then do the convolution
883                   with the instrument response.  Only after doing that will
884                   we attempt to calculate the goodness of fit matrices.  This
885                   means that we will be looping through all of the points
886                   twice, which is not worth it if there is no convolution
887                   necessary. */
888
889                for (i=fit_start; i<fit_end; i++) {
890                        int convpts;
891
892                        /* We wish to find yfit = (*pfnvals) * instr, so explicitly:
893                             yfit[i] = sum_{j=0}^i (*pfnvals)[i-j].instr[j]
894                           But instr[k]=0 for k >= ninstr, so we only need to sum:
895                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
896                           (*pfnvals)[i-j].instr[j]
897                        */
898
899                        /* Zero our adders */
900                        yfit[i] = 0;
901                        for (k=1; k<nparam; k++)
902                                (*pdy_dparam_conv)[i][k] = 0;
903
904                        convpts = (ninstr <= i) ? ninstr-1 : i;
905                        for (j=0; j<=convpts; j++) {
906                                yfit[i] += (*pfnvals)[i-j] * instr[j];
907                                for (k=1; k<nparam; k++)
908                                        (*pdy_dparam_conv)[i][k] += (*pdy_dparam_pure)[i-j][k] * instr[j];
909                        }
910                }
911        } else {
912                /* Can go straight into the final arrays in this case */
913                if (fitfunc == GCI_multiexp_lambda)
914                        ret = multiexp_lambda_array(xincr, param, yfit,
915                                                                                (*pdy_dparam_conv), fit_end, nparam);
916                else if (fitfunc == GCI_multiexp_tau)
917                        ret = multiexp_tau_array(xincr, param, yfit,
918                                                                         (*pdy_dparam_conv), fit_end, nparam);
919                else if (fitfunc == GCI_stretchedexp)
920                        ret = stretchedexp_array(xincr, param, yfit,
921                                                                         (*pdy_dparam_conv), fit_end, nparam);
922                else
923                        ret = -1;
924
925                if (ret < 0)
926                        for (i=0; i<fit_end; i++)
927                                (*fitfunc)(xincr*i, param, &yfit[i],
928                                                   (*pdy_dparam_conv)[i], nparam);
929        }
930
931        /* OK, now we've got our (possibly convolved) data, we can do the
932           rest almost exactly as above. */
933        {
934        float alpha_weight[256]; //TODO establish maximum # bins and use elsewhere (#define)
935        float beta_weight[256];
936        int q;
937        float weight;
938
939        int i_free;
940        int j_free;
941        float dot_product;
942        float beta_sum;
943        float dy_dparam_k_i;
944        float *alpha_weight_ptr;
945        float *beta_weight_ptr;
946
947        *chisq = 0.0f;
948
949                switch (noise) {
950            case NOISE_CONST:
951            {
952                float *alpha_ptr = &alpha_weight[fit_start];
953                float *beta_ptr = &beta_weight[fit_start];
954                for (q = fit_start; q < fit_end; ++q) {
955                    (*pdy_dparam_conv)[q][0] = 1.0f;
956                    yfit[q] += param[0];
957                    dy[q] = y[q] - yfit[q];
958                    weight = 1.0f;
959                    *alpha_ptr++ = weight; //
960                    //alpha_weight[q] = weight; // 1
961                    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]
965                    weight *= dy[q];
966                    *chisq += weight; // dy[q] * dy[q]
967                }
968                break;
969            }
970            case NOISE_GIVEN:
971            {
972                for (q = fit_start; q < fit_end; ++q) {
973                    (*pdy_dparam_conv)[q][0] = 1.0f;
974                    yfit[q] += param[0];
975                    dy[q] = y[q] - yfit[q];
976                    weight = 1.0f / (sig[q] * sig[q]);
977                    alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
978                    weight *= dy[q];
979                    beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
980                    weight *= dy[q];
981                    *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
982                }
983                break;
984            }
985            case NOISE_POISSON_DATA:
986            {
987                for (q = fit_start; q < fit_end; ++q) {
988                    (*pdy_dparam_conv)[q][0] = 1.0f;
989                    yfit[q] += param[0];
990                    dy[q] = y[q] - yfit[q];
991                    weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
992                    alpha_weight[q] = weight; // 1 / sig(q)
993                    weight *= dy[q];
994                    beta_weight[q] = weight; // dy[q] / sig(q)
995                    weight *= dy[q];
996                    *chisq += weight; // (dy[q] * dy[q]) / sig(q)
997                }
998                break;
999            }
1000            case NOISE_POISSON_FIT:
1001            {
1002                for (q = fit_start; q < fit_end; ++q) {
1003                    (*pdy_dparam_conv)[q][0] = 1.0f;
1004                    yfit[q] += param[0];
1005                    dy[q] = y[q] - yfit[q];
1006                    weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1007                    alpha_weight[q] = weight; // 1 / sig(q)
1008                    weight *= dy[q];
1009                    beta_weight[q] = weight; // dy(q) / sig(q)
1010                    weight *= dy[q];
1011                    *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1012                }
1013                break;
1014            }
1015            case NOISE_GAUSSIAN_FIT:
1016            {
1017                 for (q = fit_start; q < fit_end; ++q) {
1018                    (*pdy_dparam_conv)[q][0] = 1.0f;
1019                    yfit[q] += param[0];
1020                    dy[q] = y[q] - yfit[q];
1021                    weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1022                    alpha_weight[q] = weight; // 1 / sig(q)
1023                    weight *= dy[q];
1024                    beta_weight[q] = weight; // dy[q] / sig(q)
1025                    weight *= dy[q];
1026                    *chisq += weight;
1027                 }
1028                 break;
1029           }
1030            case NOISE_MLE:
1031            {
1032                for (q = fit_start; q < fit_end; ++q) {
1033                    (*pdy_dparam_conv)[q][0] = 1.0f;
1034                    yfit[q] += param[0];
1035                    dy[q] = y[q] - yfit[q];
1036                    weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f);
1037                    alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus?
1038                    beta_weight[q] = dy[q] * weight;
1039                    *chisq += (0.0f == y[q])
1040                            ? 2.0 * yfit[q]
1041                            : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1042                }
1043                if (*chisq <= 0.0f) {
1044                    *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1045                }
1046                break;
1047            }
1048            default:
1049            {
1050                return -3;
1051            }
1052        }
1053
1054        // Check if chi square worsened:
1055        if (0.0f != old_chisq && *chisq >= old_chisq) {
1056            // don't bother to set up the matrices for solution
1057            return 0;
1058        }
1059
1060        i_free = 0;
1061        // for all columns
1062        for (i = 0; i < nparam; ++i) {
1063            if (paramfree[i]) {
1064                j_free = 0;
1065                beta_sum = 0.0f;
1066                // row loop, only need to consider lower triangle
1067                for (j = 0; j <= i; ++j) {
1068                    if (paramfree[j]) {
1069                        dot_product = 0.0f;
1070                              alpha_weight_ptr = &alpha_weight[fit_start];
1071                        if (0 == j_free) {
1072                            beta_weight_ptr = &beta_weight[fit_start];
1073                            // for all data
1074                            for (k = fit_start; k < fit_end; ++k) {
1075                                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];
1078                            }
1079                        }
1080                        else {
1081                            // for all data
1082                            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];
1084                            }
1085                        } // k loop
1086
1087                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1088                       // if (i_free != j_free) {
1089                       //     // matrix is symmetric
1090                       //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1091                       // }
1092                        ++j_free;
1093                    }
1094                } // j loop
1095                beta[i_free] = beta_sum;
1096                ++i_free;
1097            }
1098            //else printf("param not free %d\n", i);
1099        } // 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        }
1111        return 0;
1112}
1113
1114/* These two variants, used just before the Marquardt fitting
1115   functions terminate, compute the function values at all points,
1116   whether or not they are being fitted.  (All points are fitted in
1117   the non-instrument response variant.)  They also compute the
1118   residuals y - yfit at all of those points and compute a chi-squared
1119   value which is not modified at small data values in the POISSON
1120   noise models.  They do not calculate alpha or beta. */
1121
1122/* This is the variant which handles an instrument response. */
1123/* We assume that the function values are sensible.  Note also that we
1124   have to be careful about which points which include in the
1125   chi-squared calculation.  Also, we are guaranteed that the
1126   initialisation of the convolution arrays has been performed. */
1127
1128int GCI_marquardt_compute_fn_final_instr(float xincr, float y[], int ndata,
1129                                   int fit_start, int fit_end,
1130                                   float instr[], int ninstr,
1131                                   noise_type noise, float sig[],
1132                                   float param[], int paramfree[], int nparam,
1133                                   void (*fitfunc)(float, float [], float *, float [], int),
1134                                   float yfit[], float dy[], float *chisq,     
1135                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1136                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1137{
1138        int i, j, mfit, ret;
1139        float sig2i;
1140        float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
1141        int fnvals_len = *pfnvals_len;
1142        int dy_dparam_nparam_size = *pdy_dparam_nparam_size;
1143
1144        /* check the necessary initialisation for safety, bail out if
1145           broken */
1146        if ((fnvals_len < ndata) || (dy_dparam_nparam_size < nparam))
1147                return -1;
1148        fnvals = *pfnvals;
1149        dy_dparam_pure = *pdy_dparam_pure;
1150        dy_dparam_conv = *pdy_dparam_conv;
1151
1152        for (j=0, mfit=0; j<nparam; j++)
1153                if (paramfree[j]) mfit++;
1154
1155        /* Calculation of the fitting data will depend upon the type of
1156           noise and the type of instrument response */
1157
1158        /* Need to calculate unconvolved values all the way down to 0 for
1159           the instrument response case */
1160        if (ninstr > 0) {
1161                if (fitfunc == GCI_multiexp_lambda)
1162                        ret = multiexp_lambda_array(xincr, param, fnvals,
1163                                                                                dy_dparam_pure, ndata, nparam);
1164                else if (fitfunc == GCI_multiexp_tau)
1165                        ret = multiexp_tau_array(xincr, param, fnvals,
1166                                                                         dy_dparam_pure, ndata, nparam);
1167                else if (fitfunc == GCI_stretchedexp)
1168                        ret = stretchedexp_array(xincr, param, fnvals,
1169                                                                         dy_dparam_pure, ndata, nparam);
1170                else
1171                        ret = -1;
1172
1173                if (ret < 0)
1174                        for (i=0; i<ndata; i++)
1175                                (*fitfunc)(xincr*i, param, &fnvals[i],
1176                                                   dy_dparam_pure[i], nparam);
1177
1178                /* OK, we've got to convolve the model fit with the given
1179                   instrument response.  What we'll do here, then, is to
1180                   generate the whole model fit first, then do the convolution
1181                   with the instrument response.  Only after doing that will
1182                   we attempt to calculate the goodness of fit matrices.  This
1183                   means that we will be looping through all of the points
1184                   twice, which is not worth it if there is no convolution
1185                   necessary. */
1186
1187                for (i=0; i<ndata; i++) {
1188                        int convpts;
1189               
1190                        /* We wish to find yfit = fnvals * instr, so explicitly:
1191                             yfit[i] = sum_{j=0}^i fnvals[i-j].instr[j]
1192                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1193                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1194                           fnvals[i-j].instr[j]
1195                        */
1196
1197                        /* Zero our adder; don't need to bother with dy_dparam
1198                           stuff here */
1199                        yfit[i] = 0;
1200
1201                        convpts = (ninstr <= i) ? ninstr-1 : i;
1202                        for (j=0; j<=convpts; j++)
1203                                yfit[i] += fnvals[i-j] * instr[j];
1204                }
1205        } else {
1206                /* Can go straight into the final arrays in this case */
1207                if (fitfunc == GCI_multiexp_lambda)
1208                        ret = multiexp_lambda_array(xincr, param, yfit,
1209                                                                                dy_dparam_conv, ndata, nparam);
1210                else if (fitfunc == GCI_multiexp_tau)
1211                        ret = multiexp_tau_array(xincr, param, yfit,
1212                                                                         dy_dparam_conv, ndata, nparam);
1213                else if (fitfunc == GCI_stretchedexp)
1214                        ret = stretchedexp_array(xincr, param, yfit,
1215                                                                         dy_dparam_conv, ndata, nparam);
1216                else
1217                        ret = -1;
1218
1219                if (ret < 0)
1220                        for (i=0; i<ndata; i++)
1221                                (*fitfunc)(xincr*i, param, &yfit[i],
1222                                                   dy_dparam_conv[i], nparam);
1223        }
1224         
1225        /* OK, now we've got our (possibly convolved) data, we can do the
1226           rest almost exactly as above. */
1227
1228        switch (noise) {
1229        case NOISE_CONST:
1230                *chisq = 0.0;
1231                /* Summation loop over all data */
1232                for (i=0; i<fit_start; i++) {
1233                        yfit[i] += param[0];
1234                        dy[i] = y[i] - yfit[i];
1235                }
1236                for ( ; i<fit_end; i++) {
1237                        yfit[i] += param[0];
1238                        dy[i] = y[i] - yfit[i];
1239                        /* And find chi^2 */
1240                        *chisq += dy[i] * dy[i];
1241                }
1242                for ( ; i<ndata; i++) {
1243                        yfit[i] += param[0];
1244                        dy[i] = y[i] - yfit[i];
1245                }
1246
1247                /* Now divide chi-squared by sigma^2 */
1248                sig2i = 1.0 / (sig[0] * sig[0]);
1249                *chisq *= sig2i;
1250                break;
1251
1252        case NOISE_GIVEN:  /* This is essentially the NR version */
1253                *chisq = 0.0;
1254                /* Summation loop over all data */
1255                for (i=0; i<fit_start; i++) {
1256                        yfit[i] += param[0];
1257                        dy[i] = y[i] - yfit[i];
1258                }
1259                for ( ; i<fit_end; i++) {
1260                        yfit[i] += param[0];
1261                        dy[i] = y[i] - yfit[i];
1262                        /* And find chi^2 */
1263                        sig2i = 1.0 / (sig[i] * sig[i]);
1264                        *chisq += dy[i] * dy[i] * sig2i;
1265                }
1266                for ( ; i<ndata; i++) {
1267                        yfit[i] += param[0];
1268                        dy[i] = y[i] - yfit[i];
1269                }
1270                break;
1271
1272        case NOISE_POISSON_DATA:
1273                *chisq = 0.0;
1274                /* Summation loop over all data */
1275                for (i=0; i<fit_start; i++) {
1276                        yfit[i] += param[0];
1277                        dy[i] = y[i] - yfit[i];
1278                }
1279                for ( ; i<fit_end; i++) {
1280                        yfit[i] += param[0];
1281                        dy[i] = y[i] - yfit[i];
1282                        /* And find chi^2 */
1283                        /* we still don't let the variance drop below 1 */
1284                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1285                        *chisq += dy[i] * dy[i] * sig2i;
1286                }
1287                for (; i<ndata; i++) {
1288                        yfit[i] += param[0];
1289                        dy[i] = y[i] - yfit[i];
1290                }
1291                break;
1292
1293        case NOISE_POISSON_FIT:
1294                *chisq = 0.0;
1295                // Summation loop over all data
1296                for (i=0; i<fit_start; i++) {
1297                        yfit[i] += param[0];
1298                        dy[i] = y[i] - yfit[i];
1299                }
1300                for ( ; i<fit_end; i++) {
1301                        yfit[i] += param[0];
1302                        dy[i] = y[i] - yfit[i];
1303                        // And find chi^2
1304                        // we still don't let the variance drop below 1
1305                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1306                        *chisq += dy[i] * dy[i] * sig2i;
1307                }
1308                for ( ; i<ndata; i++) {
1309                        yfit[i] += param[0];
1310                        dy[i] = y[i] - yfit[i];
1311                }
1312                break;
1313               
1314        case NOISE_MLE:                              // for the final chisq report a normal chisq measure for MLE
1315                *chisq = 0.0;
1316                // Summation loop over all data
1317                for (i=0; i<fit_start; i++) {
1318                        yfit[i] += param[0];
1319                        dy[i] = y[i] - yfit[i];
1320                }
1321                for ( ; i<fit_end; i++) {
1322                        yfit[i] += param[0];
1323                        dy[i] = y[i] - yfit[i];
1324                        // And find chi^2
1325//                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1326                        if (yfit[i]<=0.0)
1327                                ; // do nothing
1328                        else if (y[i]==0.0)
1329                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1330                        else
1331                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1332                }
1333                for ( ; i<ndata; i++) {
1334                        yfit[i] += param[0];
1335                        dy[i] = y[i] - yfit[i];
1336                }
1337                if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1338                break;
1339
1340        case NOISE_GAUSSIAN_FIT:
1341                *chisq = 0.0;
1342                // Summation loop over all data
1343                for (i=0; i<fit_start; i++) {
1344                        yfit[i] += param[0];
1345                        dy[i] = y[i] - yfit[i];
1346                }
1347                for ( ; i<fit_end; i++) {
1348                        yfit[i] += param[0];
1349                        dy[i] = y[i] - yfit[i];
1350                        // And find chi^2
1351                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1352                        *chisq += dy[i] * dy[i] * sig2i;
1353                }
1354                for ( ; i<ndata; i++) {
1355                        yfit[i] += param[0];
1356                        dy[i] = y[i] - yfit[i];
1357                }
1358                break;
1359       
1360        default:
1361                return -3;
1362                /* break; */   // (unreachable code)
1363        }
1364
1365        return 0;
1366}
1367
1368
1369//********************************* GCI_marquardt_fitting_engine **********************************************************************
1370
1371// This returns the number of iterations or negative if an error occurred.
1372// passes all the data to the ecf routine, checks the returned chisq and re-fits if it is of benefit
1373// was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03
1374
1375int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end, 
1376                                                float prompt[], int nprompt,
1377                                                noise_type noise, float sig[],
1378                                                float param[], int paramfree[],
1379                                           int nparam, restrain_type restrain,
1380                                           void (*fitfunc)(float, float [], float *, float [], int),
1381                                           float *fitted, float *residuals, float *chisq,
1382                                           float **covar, float **alpha, float **erraxes,
1383                                           float chisq_target, float chisq_delta, int chisq_percent)
1384{
1385        float oldChisq, local_chisq;
1386        int ret, tries=0;
1387
1388        if (ecf_exportParams) ecf_ExportParams_OpenFile ();
1389
1390        // All of the work is done by the ECF module
1391        ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1392                                                          prompt, nprompt, noise, sig,
1393                                                          param, paramfree, nparam, restrain, fitfunc,
1394                                                          fitted, residuals, covar, alpha, &local_chisq,
1395                                                          chisq_delta, chisq_percent, erraxes);
1396                                                         
1397        // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same
1398        // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr()
1399        oldChisq = 3.0e38;
1400        while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<(MAXREFITS*3))
1401        {
1402                oldChisq = local_chisq; 
1403                tries++;
1404                ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1405                                                          prompt, nprompt, noise, sig,
1406                                                          param, paramfree, nparam, restrain, fitfunc,
1407                                                          fitted, residuals, covar, alpha, &local_chisq,
1408                                                          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();*/
1422
1423        if (chisq!=NULL) *chisq = local_chisq;
1424
1425        if (ecf_exportParams) ecf_ExportParams_CloseFile ();
1426
1427        return ret;             // summed number of iterations
1428}
1429
1430/* Cleanup function */
1431void GCI_marquardt_cleanup(void)
1432{
1433/*      if (fnvals_len) {
1434                free(fnvals);
1435                GCI_ecf_free_matrix(dy_dparam_pure);
1436                GCI_ecf_free_matrix(dy_dparam_conv);
1437                fnvals_len = 0;
1438                dy_dparam_nparam_size = 0;
1439        }
1440*/
1441}
1442
1443// Emacs settings:
1444// Local variables:
1445// mode: c
1446// c-basic-offset: 4
1447// tab-width: 4
1448// End:
Note: See TracBrowser for help on using the repository browser.