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

Revision 7606, 46.4 KB checked in by paulbarber, 9 years ago (diff)

To compile in TRI2, moved declarations to top of relevant block.
Added empty GCI_marquardt.
Added matrix array routinnes.
Added Global and SPA files.

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