source: branches/maven/projects/slim-curve/src/main/c/EcfSingle.c @ 7387

Revision 7387, 46.0 KB checked in by aivar, 9 years ago (diff)

A version of Paul Barber's ECF library with functions used in SLIM Plugin. These have been rewritten to remove NR dependencies. Functions not used and other code that still uses NR all removed.

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#define do_frees \
493        if (fnvals) free(fnvals);\
494        if (dy_dparam_pure) GCI_ecf_free_matrix(dy_dparam_pure);\
495        if (dy_dparam_conv) GCI_ecf_free_matrix(dy_dparam_conv);
496
497int GCI_marquardt_instr(float xincr, float y[],
498                                        int ndata, int fit_start, int fit_end,
499                                        float instr[], int ninstr,
500                                        noise_type noise, float sig[],
501                                        float param[], int paramfree[], int nparam,
502                                        restrain_type restrain,
503                                        void (*fitfunc)(float, float [], float *, float [], int),
504                                        float *fitted, float *residuals,
505                                        float **covar, float **alpha, float *chisq,
506                                        float chisq_delta, float chisq_percent, float **erraxes)
507{
508        float alambda, ochisq;
509        int mfit, mfit2;
510        float evals[MAXFIT];
511        int i, k, itst, itst_max;
512
513        // The following are declared here to retain some optimisation by not repeatedly mallocing
514        // (only once per transient), but to remain thread safe.
515        // They are malloced by lower fns but at the end, freed by this fn.
516        // These vars were global or static before thread safety was introduced.
517        float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL;
518        int fnvals_len=0, dy_dparam_nparam_size=0;
519        float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 
520
521        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
522
523        mfit = 0;
524        for (i=0; i<nparam; i++) {
525                if (paramfree[i]) mfit++;
526        }
527
528        if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
529
530        alambda = -1;
531        if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
532                                                                 instr, ninstr, noise, sig,
533                                                                 param, paramfree, nparam, restrain,
534                                                                 fitfunc, fitted, residuals,
535                                                                 covar, alpha, chisq, &alambda, 
536                                                                 &mfit2, &ochisq2, paramtry, beta, dparam,
537                                                                 &fnvals, &dy_dparam_pure, &dy_dparam_conv,
538                                                                 &fnvals_len, &dy_dparam_nparam_size) != 0) {
539                do_frees
540                return -1;
541        }
542
543        if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
544
545        k = 1;  /* Iteration counter */
546        itst = 0;
547        for (;;) {
548                k++;
549                if (k > MAXITERS) {
550                        do_frees
551                        return -2;
552                }
553
554                ochisq = *chisq;
555                if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
556                                                                         instr, ninstr, noise, sig,
557                                                                         param, paramfree, nparam, restrain,
558                                                                         fitfunc, fitted, residuals,
559                                                                         covar, alpha, chisq, &alambda, 
560                                                                         &mfit2, &ochisq2, paramtry, beta, dparam,
561                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv,
562                                                                         &fnvals_len, &dy_dparam_nparam_size) != 0) {
563                        do_frees
564                        return -3;
565                }
566
567                if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
568
569                if (*chisq > ochisq)
570                        itst = 0;
571                else if (ochisq - *chisq < chisq_delta)
572                        itst++;
573
574                if (itst < itst_max) continue;
575
576                /* Endgame */
577                alambda = 0.0;
578                if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
579                                                                         instr, ninstr, noise, sig,
580                                                                         param, paramfree, nparam, restrain,
581                                                                         fitfunc, fitted, residuals,
582                                                                         covar, alpha, chisq, &alambda, 
583                                                                         &mfit2, &ochisq2, paramtry, beta, dparam,
584                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv,
585                                                                         &fnvals_len, &dy_dparam_nparam_size) != 0) {
586                        do_frees
587                        return -4;
588                }
589
590                if (erraxes == NULL){
591                        do_frees
592                        return k;
593                }
594
595                break;  /* We're done now */
596        }
597
598        do_frees
599        return k;
600}
601#undef do_frees
602
603
604int GCI_gauss_jordan(float **a, int n, float *b);
605
606int GCI_marquardt_step_instr(float xincr, float y[],
607                                        int ndata, int fit_start, int fit_end,
608                                        float instr[], int ninstr,
609                                        noise_type noise, float sig[],
610                                        float param[], int paramfree[], int nparam,
611                                        restrain_type restrain,
612                                        void (*fitfunc)(float, float [], float *, float [], int),
613                                        float yfit[], float dy[],
614                                        float **covar, float **alpha, float *chisq,
615                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,       
616                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
617                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
618{
619        int j, k, l, ret;
620//      static int mfit;   // was static but now thread safe
621//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe
622
623        /*if (nparam > MAXFIT) {
624            printf("GCI_marq_step_instr returns -10\n");
625                return -10;
626        }
627        if (xincr <= 0) {
628            printf("GCI_marq_step_instr returns -11\n");
629                return -11;
630        }
631        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata) {
632            printf("GCI_marq_step_instr returns -12\n");
633                return -12;
634        }*/
635
636        /* Initialisation */
637        /* We assume we're given sensible starting values for param[] */
638        if (*alambda < 0.0) {
639
640                for ((*pmfit)=0, j=0; j<nparam; j++)
641                        if (paramfree[j])
642                                (*pmfit)++;
643    global_chisq = 0.0f;
644                if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
645                                                                                   instr, ninstr, noise, sig,
646                                                                                   param, paramfree, nparam, fitfunc,
647                                                                                   yfit, dy, alpha, beta, chisq,
648                                                                                   *alambda,   
649                                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
650                                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0) {
651                        return -2;
652                }
653
654                *alambda = 0.001;
655                *pochisq = *chisq;
656                for (j=0; j<nparam; j++)
657                        paramtry[j] = param[j];
658
659        }
660
661        /* Alter linearised fitting matrix by augmenting diagonal elements */
662        for (j=0; j<(*pmfit); j++) {
663                for (k=0; k<(*pmfit); k++)
664                        covar[j][k] = alpha[j][k];
665
666                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
667                dparam[j] = beta[j];
668        }
669
670        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */
671        if (GCI_gauss_jordan(covar, *pmfit, dparam) != 0) {
672                return -1;
673        }
674       
675
676        /* Once converged, evaluate covariance matrix */
677        if (*alambda == 0) {
678                if (GCI_marquardt_compute_fn_final_instr(
679                                                                        xincr, y, ndata, fit_start, fit_end,
680                                                                        instr, ninstr, noise, sig,
681                                                                        param, paramfree, nparam, fitfunc,
682                                                                        yfit, dy, chisq,       
683                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
684                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0) {
685                        return -3;
686                }
687                if (*pmfit < nparam) {  /* no need to do this otherwise */
688                        GCI_covar_sort(covar, nparam, paramfree, *pmfit);
689                        GCI_covar_sort(alpha, nparam, paramfree, *pmfit);
690                }
691                return 0;
692        }
693
694        /* Did the trial succeed? */
695        for (j=0, l=0; l<nparam; l++)
696                if (paramfree[l])
697                        paramtry[l] = param[l] + dparam[j++];
698
699        if (restrain == ECF_RESTRAIN_DEFAULT)
700                ret = check_ecf_params (paramtry, nparam, fitfunc);
701        else
702                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
703
704        if (ret != 0) {
705                /* Bad parameters, increase alambda and return */
706                *alambda *= 10.0;
707                return 0;
708        }
709    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
710        if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
711                                                                           instr, ninstr, noise, sig,
712                                                                           paramtry, paramfree, nparam, fitfunc,
713                                                                           yfit, dy, covar, dparam,
714                                                                           chisq, *alambda,     
715                                                                           pfnvals, pdy_dparam_pure, pdy_dparam_conv,
716                                                                           pfnvals_len, pdy_dparam_nparam_size) != 0) {
717                return -2;
718        }
719
720        /* Success, accept the new solution */
721        if (*chisq < *pochisq) {
722                *alambda *= 0.1;
723                *pochisq = *chisq;
724                for (j=0; j<(*pmfit); j++) {
725                        for (k=0; k<(*pmfit); k++)
726                                alpha[j][k] = covar[j][k];
727                        beta[j] = dparam[j];
728                }
729                for (l=0; l<nparam; l++)
730                        param[l] = paramtry[l];
731        } else { /* Failure, increase alambda and return */
732                *alambda *= 10.0;
733                *chisq = *pochisq;
734        }
735
736        return 0;
737}
738
739
740/* Used by GCI_marquardt to evaluate the linearised fitting matrix alpha
741   and vector beta and to calculate chi^2.      The equations involved are
742   given in section 15.5 of Numerical Recipes; basically:
743
744     \chi^2(param) = \sum_{i=1}^N ( (y_i-y(x_i;param)) / sigma_i )^2
745
746     beta_k = -1/2 (d/dparam_k)(chi^2)
747
748     alpha_kl = \sum_{i=1}^N (1/sigma_i^2) .
749                  (dy(x_i;param)/dparam_k) . (dy(x_i;param)/dparam_l)
750
751   (where all of the derivatives are partial).
752
753   If an instrument response is provided, we also take account of it
754   now.  We are given that:
755
756     observed(t) = response(t) * instr(t)
757
758   where response(t) is being fitted with fitfunc; it is also trivial to
759   show that (assuming that instr(t) is known and fixed, with no
760   dependencies on the param_k, the parameters being fitted):
761
762     (d/dparam_k) observed(t) = ((d/dparam_k) response(t)) * instr(t)
763
764   so we do not need to alter the response function in any way to
765   determined the fitted convolved response.
766
767   This is the variant which handles an instrument response. */
768
769/* We assume that the function values are sensible. */
770
771int GCI_marquardt_compute_fn_instr(float xincr, float y[], int ndata,
772                                   int fit_start, int fit_end,
773                                   float instr[], int ninstr,
774                                   noise_type noise, float sig[],
775                                   float param[], int paramfree[], int nparam,
776                                   void (*fitfunc)(float, float [], float *, float [], int),
777                                   float yfit[], float dy[],
778                                   float **alpha, float beta[], float *chisq,
779                                   float alambda,
780                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
781                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
782{
783        int i, j, k, l, m, mfit, ret;
784        float wt, sig2i, y_ymod;
785
786        /* Are we initialising? */
787        // Malloc the arrays that will get used again in this fit via the pointers passed in
788        // They will be freed by the higher fn that declared them.
789        if (alambda < 0) {
790                /* do any necessary initialisation */
791                if ((*pfnvals_len) < ndata) {  /* we will need this length for
792                                                                          the final full computation */
793                        if ((*pfnvals_len)) {
794                                free((*pfnvals));
795                                GCI_ecf_free_matrix((*pdy_dparam_pure));
796                                GCI_ecf_free_matrix((*pdy_dparam_conv));
797                                (*pfnvals_len) = 0;
798                                (*pdy_dparam_nparam_size) = 0;
799                        }
800                } else if ((*pdy_dparam_nparam_size) < nparam) {
801                        GCI_ecf_free_matrix((*pdy_dparam_pure));
802                        GCI_ecf_free_matrix((*pdy_dparam_conv));
803                        (*pdy_dparam_nparam_size) = 0;
804                }
805                if (! (*pfnvals_len)) {
806                        if (((*pfnvals) = (float *) malloc(ndata * sizeof(float)))
807                                == NULL)
808                                return -1;
809                        (*pfnvals_len) = ndata;
810                }
811                if (! (*pdy_dparam_nparam_size)) {
812                        /* use (*pfnvals_len), not ndata here! */
813                        if (((*pdy_dparam_pure) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
814                                /* Can't keep (*pfnvals) around in this case */
815                                free((*pfnvals));
816                                (*pfnvals_len) = 0;
817                                return -1;
818                        }
819                        if (((*pdy_dparam_conv) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
820                                /* Can't keep (*pfnvals) around in this case */
821                                free((*pfnvals));
822                                free((*pdy_dparam_pure));
823                                (*pfnvals_len) = 0;
824                                return -1;
825                        }
826                        (*pdy_dparam_nparam_size) = nparam;
827                }
828        }
829
830        for (j=0, mfit=0; j<nparam; j++)
831                if (paramfree[j]) mfit++;
832
833        /* Calculation of the fitting data will depend upon the type of
834           noise and the type of instrument response */
835
836        /* Need to calculate unconvolved values all the way down to 0 for
837           the instrument response case */
838        if (ninstr > 0) {
839                if (fitfunc == GCI_multiexp_lambda)
840                        ret = multiexp_lambda_array(xincr, param, (*pfnvals),
841                                                                                (*pdy_dparam_pure), fit_end, nparam);
842                else if (fitfunc == GCI_multiexp_tau)
843                        ret = multiexp_tau_array(xincr, param, (*pfnvals),
844                                                                         (*pdy_dparam_pure), fit_end, nparam);
845                else if (fitfunc == GCI_stretchedexp)
846                        ret = stretchedexp_array(xincr, param, (*pfnvals),
847                                                                         (*pdy_dparam_pure), fit_end, nparam);
848                else
849                        ret = -1;
850
851                if (ret < 0)
852                        for (i=0; i<fit_end; i++)
853                                (*fitfunc)(xincr*i, param, &(*pfnvals)[i],
854                                                   (*pdy_dparam_pure)[i], nparam);
855
856                /* OK, we've got to convolve the model fit with the given
857                   instrument response.  What we'll do here, then, is to
858                   generate the whole model fit first, then do the convolution
859                   with the instrument response.  Only after doing that will
860                   we attempt to calculate the goodness of fit matrices.  This
861                   means that we will be looping through all of the points
862                   twice, which is not worth it if there is no convolution
863                   necessary. */
864
865                for (i=fit_start; i<fit_end; i++) {
866                        int convpts;
867
868                        /* We wish to find yfit = (*pfnvals) * instr, so explicitly:
869                             yfit[i] = sum_{j=0}^i (*pfnvals)[i-j].instr[j]
870                           But instr[k]=0 for k >= ninstr, so we only need to sum:
871                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
872                           (*pfnvals)[i-j].instr[j]
873                        */
874
875                        /* Zero our adders */
876                        yfit[i] = 0;
877                        for (k=1; k<nparam; k++)
878                                (*pdy_dparam_conv)[i][k] = 0;
879
880                        convpts = (ninstr <= i) ? ninstr-1 : i;
881                        for (j=0; j<=convpts; j++) {
882                                yfit[i] += (*pfnvals)[i-j] * instr[j];
883                                for (k=1; k<nparam; k++)
884                                        (*pdy_dparam_conv)[i][k] += (*pdy_dparam_pure)[i-j][k] * instr[j];
885                        }
886                }
887        } else {
888                /* Can go straight into the final arrays in this case */
889                if (fitfunc == GCI_multiexp_lambda)
890                        ret = multiexp_lambda_array(xincr, param, yfit,
891                                                                                (*pdy_dparam_conv), fit_end, nparam);
892                else if (fitfunc == GCI_multiexp_tau)
893                        ret = multiexp_tau_array(xincr, param, yfit,
894                                                                         (*pdy_dparam_conv), fit_end, nparam);
895                else if (fitfunc == GCI_stretchedexp)
896                        ret = stretchedexp_array(xincr, param, yfit,
897                                                                         (*pdy_dparam_conv), fit_end, nparam);
898                else
899                        ret = -1;
900
901                if (ret < 0)
902                        for (i=0; i<fit_end; i++)
903                                (*fitfunc)(xincr*i, param, &yfit[i],
904                                                   (*pdy_dparam_conv)[i], nparam);
905        }
906
907        /* OK, now we've got our (possibly convolved) data, we can do the
908           rest almost exactly as above. */
909
910        float alpha_weight[256]; //TODO establish maximum # bins and use elsewhere (#define)
911        float beta_weight[256];
912
913        *chisq = 0.0f;
914        int q;
915        float weight;
916        switch (noise) {
917            case NOISE_CONST:
918            {
919                float *alpha_ptr = &alpha_weight[fit_start];
920                float *beta_ptr = &beta_weight[fit_start];
921                for (q = fit_start; q < fit_end; ++q) {
922                    (*pdy_dparam_conv)[q][0] = 1.0f;
923                    yfit[q] += param[0];
924                    dy[q] = y[q] - yfit[q];
925                    weight = 1.0f;
926                    *alpha_ptr++ = weight; //
927                    //alpha_weight[q] = weight; // 1
928                    weight *= dy[q];
929                    //printf("beta weight %d %f is y %f - yfit %f\n", q, weight, y[q], yfit[q]);
930                    *beta_ptr++ = weight; //
931                    //beta_weight[q] = weight; // dy[q]
932                    weight *= dy[q];
933                    *chisq += weight; // dy[q] * dy[q]
934                }
935                break;
936            }
937            case NOISE_GIVEN:
938            {
939                for (q = fit_start; q < fit_end; ++q) {
940                    (*pdy_dparam_conv)[q][0] = 1.0f;
941                    yfit[q] += param[0];
942                    dy[q] = y[q] - yfit[q];
943                    weight = 1.0f / (sig[q] * sig[q]);
944                    alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
945                    weight *= dy[q];
946                    beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
947                    weight *= dy[q];
948                    *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
949                }
950                break;
951            }
952            case NOISE_POISSON_DATA:
953            {
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 = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
959                    alpha_weight[q] = weight; // 1 / sig(q)
960                    weight *= dy[q];
961                    beta_weight[q] = weight; // dy[q] / sig(q)
962                    weight *= dy[q];
963                    *chisq += weight; // (dy[q] * dy[q]) / sig(q)
964                }
965                break;
966            }
967            case NOISE_POISSON_FIT:
968            {
969                for (q = fit_start; q < fit_end; ++q) {
970                    (*pdy_dparam_conv)[q][0] = 1.0f;
971                    yfit[q] += param[0];
972                    dy[q] = y[q] - yfit[q];
973                    weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
974                    alpha_weight[q] = weight; // 1 / sig(q)
975                    weight *= dy[q];
976                    beta_weight[q] = weight; // dy(q) / sig(q)
977                    weight *= dy[q];
978                    *chisq += weight; // (dy(q) * dy(q)) / sig(q)
979                }
980                break;
981            }
982            case NOISE_GAUSSIAN_FIT:
983            {
984                 for (q = fit_start; q < fit_end; ++q) {
985                    (*pdy_dparam_conv)[q][0] = 1.0f;
986                    yfit[q] += param[0];
987                    dy[q] = y[q] - yfit[q];
988                    weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
989                    alpha_weight[q] = weight; // 1 / sig(q)
990                    weight *= dy[q];
991                    beta_weight[q] = weight; // dy[q] / sig(q)
992                    weight *= dy[q];
993                    *chisq += weight;
994                 }
995                 break;
996           }
997            case NOISE_MLE:
998            {
999                for (q = fit_start; q < fit_end; ++q) {
1000                    (*pdy_dparam_conv)[q][0] = 1.0f;
1001                    yfit[q] += param[0];
1002                    dy[q] = y[q] - yfit[q];
1003                    weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f);
1004                    alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus?
1005                    beta_weight[q] = dy[q] * weight;
1006                    *chisq += (0.0f == y[q])
1007                            ? 2.0 * yfit[q]
1008                            : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1009                }
1010                if (*chisq <= 0.0f) {
1011                    *chisq = 1.0e308f; // don't let chisq=0 through yfit being all -ve
1012                }
1013                break;
1014            }
1015            default:
1016            {
1017                return -3;
1018            }
1019        }
1020
1021        // Check if chi square worsened:
1022        if (0.0f != global_chisq && *chisq >= global_chisq) { //TODO pass in the old chi square as a parameter
1023            // don't bother to set up the matrices for solution
1024            return 0;
1025        }
1026
1027        int i_free;
1028        int j_free;
1029        float dot_product;
1030        float beta_sum;
1031        float dy_dparam_k_i;
1032        float *alpha_weight_ptr;
1033        float *beta_weight_ptr;
1034
1035        i_free = 0;
1036        // for all columns
1037        for (i = 0; i < nparam; ++i) {
1038            if (paramfree[i]) {
1039                j_free = 0;
1040                beta_sum = 0.0f;
1041                // row loop, only need to consider lower triangle
1042                for (j = 0; j <= i; ++j) {
1043                    if (paramfree[j]) {
1044                        dot_product = 0.0f;
1045                              alpha_weight_ptr = &alpha_weight[fit_start];
1046                        if (0 == j_free) {
1047                            beta_weight_ptr = &beta_weight[fit_start];
1048                            // for all data
1049                            for (k = fit_start; k < fit_end; ++k) {
1050                                dy_dparam_k_i = (*pdy_dparam_conv)[k][i];
1051                                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.
1052                                beta_sum += dy_dparam_k_i * (*beta_weight_ptr++); ///beta_weight[k];
1053                            }
1054                        }
1055                        else {
1056                            // for all data
1057                            for (k = fit_start; k < fit_end; ++k) {
1058                                dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * (*alpha_weight_ptr++); // alpha_weight[k];
1059                            }
1060                        } // k loop
1061
1062                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1063                       // if (i_free != j_free) {
1064                       //     // matrix is symmetric
1065                       //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1066                       // }
1067                        ++j_free;
1068                    }
1069                } // j loop
1070                beta[i_free] = beta_sum;
1071                ++i_free;
1072            }
1073        } // i loop
1074
1075        return 0;
1076}
1077
1078/* These two variants, used just before the Marquardt fitting
1079   functions terminate, compute the function values at all points,
1080   whether or not they are being fitted.  (All points are fitted in
1081   the non-instrument response variant.)  They also compute the
1082   residuals y - yfit at all of those points and compute a chi-squared
1083   value which is not modified at small data values in the POISSON
1084   noise models.  They do not calculate alpha or beta. */
1085
1086/* This is the variant which handles an instrument response. */
1087/* We assume that the function values are sensible.  Note also that we
1088   have to be careful about which points which include in the
1089   chi-squared calculation.  Also, we are guaranteed that the
1090   initialisation of the convolution arrays has been performed. */
1091
1092int GCI_marquardt_compute_fn_final_instr(float xincr, float y[], int ndata,
1093                                   int fit_start, int fit_end,
1094                                   float instr[], int ninstr,
1095                                   noise_type noise, float sig[],
1096                                   float param[], int paramfree[], int nparam,
1097                                   void (*fitfunc)(float, float [], float *, float [], int),
1098                                   float yfit[], float dy[], float *chisq,     
1099                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1100                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1101{
1102        int i, j, mfit, ret;
1103        float sig2i;
1104        float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
1105        int fnvals_len = *pfnvals_len;
1106        int dy_dparam_nparam_size = *pdy_dparam_nparam_size;
1107
1108        /* check the necessary initialisation for safety, bail out if
1109           broken */
1110        if ((fnvals_len < ndata) || (dy_dparam_nparam_size < nparam))
1111                return -1;
1112        fnvals = *pfnvals;
1113        dy_dparam_pure = *pdy_dparam_pure;
1114        dy_dparam_conv = *pdy_dparam_conv;
1115
1116        for (j=0, mfit=0; j<nparam; j++)
1117                if (paramfree[j]) mfit++;
1118
1119        /* Calculation of the fitting data will depend upon the type of
1120           noise and the type of instrument response */
1121
1122        /* Need to calculate unconvolved values all the way down to 0 for
1123           the instrument response case */
1124        if (ninstr > 0) {
1125                if (fitfunc == GCI_multiexp_lambda)
1126                        ret = multiexp_lambda_array(xincr, param, fnvals,
1127                                                                                dy_dparam_pure, ndata, nparam);
1128                else if (fitfunc == GCI_multiexp_tau)
1129                        ret = multiexp_tau_array(xincr, param, fnvals,
1130                                                                         dy_dparam_pure, ndata, nparam);
1131                else if (fitfunc == GCI_stretchedexp)
1132                        ret = stretchedexp_array(xincr, param, fnvals,
1133                                                                         dy_dparam_pure, ndata, nparam);
1134                else
1135                        ret = -1;
1136
1137                if (ret < 0)
1138                        for (i=0; i<ndata; i++)
1139                                (*fitfunc)(xincr*i, param, &fnvals[i],
1140                                                   dy_dparam_pure[i], nparam);
1141
1142                /* OK, we've got to convolve the model fit with the given
1143                   instrument response.  What we'll do here, then, is to
1144                   generate the whole model fit first, then do the convolution
1145                   with the instrument response.  Only after doing that will
1146                   we attempt to calculate the goodness of fit matrices.  This
1147                   means that we will be looping through all of the points
1148                   twice, which is not worth it if there is no convolution
1149                   necessary. */
1150
1151                for (i=0; i<ndata; i++) {
1152                        int convpts;
1153               
1154                        /* We wish to find yfit = fnvals * instr, so explicitly:
1155                             yfit[i] = sum_{j=0}^i fnvals[i-j].instr[j]
1156                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1157                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1158                           fnvals[i-j].instr[j]
1159                        */
1160
1161                        /* Zero our adder; don't need to bother with dy_dparam
1162                           stuff here */
1163                        yfit[i] = 0;
1164
1165                        convpts = (ninstr <= i) ? ninstr-1 : i;
1166                        for (j=0; j<=convpts; j++)
1167                                yfit[i] += fnvals[i-j] * instr[j];
1168                }
1169        } else {
1170                /* Can go straight into the final arrays in this case */
1171                if (fitfunc == GCI_multiexp_lambda)
1172                        ret = multiexp_lambda_array(xincr, param, yfit,
1173                                                                                dy_dparam_conv, ndata, nparam);
1174                else if (fitfunc == GCI_multiexp_tau)
1175                        ret = multiexp_tau_array(xincr, param, yfit,
1176                                                                         dy_dparam_conv, ndata, nparam);
1177                else if (fitfunc == GCI_stretchedexp)
1178                        ret = stretchedexp_array(xincr, param, yfit,
1179                                                                         dy_dparam_conv, ndata, nparam);
1180                else
1181                        ret = -1;
1182
1183                if (ret < 0)
1184                        for (i=0; i<ndata; i++)
1185                                (*fitfunc)(xincr*i, param, &yfit[i],
1186                                                   dy_dparam_conv[i], nparam);
1187        }
1188         
1189        /* OK, now we've got our (possibly convolved) data, we can do the
1190           rest almost exactly as above. */
1191
1192        switch (noise) {
1193        case NOISE_CONST:
1194                *chisq = 0.0;
1195                /* Summation loop over all data */
1196                for (i=0; i<fit_start; i++) {
1197                        yfit[i] += param[0];
1198                        dy[i] = y[i] - yfit[i];
1199                }
1200                for ( ; i<fit_end; i++) {
1201                        yfit[i] += param[0];
1202                        dy[i] = y[i] - yfit[i];
1203                        /* And find chi^2 */
1204                        *chisq += dy[i] * dy[i];
1205                }
1206                for ( ; i<ndata; i++) {
1207                        yfit[i] += param[0];
1208                        dy[i] = y[i] - yfit[i];
1209                }
1210
1211                /* Now divide chi-squared by sigma^2 */
1212                sig2i = 1.0 / (sig[0] * sig[0]);
1213                *chisq *= sig2i;
1214                break;
1215
1216        case NOISE_GIVEN:  /* This is essentially the NR version */
1217                *chisq = 0.0;
1218                /* Summation loop over all data */
1219                for (i=0; i<fit_start; i++) {
1220                        yfit[i] += param[0];
1221                        dy[i] = y[i] - yfit[i];
1222                }
1223                for ( ; i<fit_end; i++) {
1224                        yfit[i] += param[0];
1225                        dy[i] = y[i] - yfit[i];
1226                        /* And find chi^2 */
1227                        sig2i = 1.0 / (sig[i] * sig[i]);
1228                        *chisq += dy[i] * dy[i] * sig2i;
1229                }
1230                for ( ; i<ndata; i++) {
1231                        yfit[i] += param[0];
1232                        dy[i] = y[i] - yfit[i];
1233                }
1234                break;
1235
1236        case NOISE_POISSON_DATA:
1237                *chisq = 0.0;
1238                /* Summation loop over all data */
1239                for (i=0; i<fit_start; i++) {
1240                        yfit[i] += param[0];
1241                        dy[i] = y[i] - yfit[i];
1242                }
1243                for ( ; i<fit_end; i++) {
1244                        yfit[i] += param[0];
1245                        dy[i] = y[i] - yfit[i];
1246                        /* And find chi^2 */
1247                        /* we still don't let the variance drop below 1 */
1248                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1249                        *chisq += dy[i] * dy[i] * sig2i;
1250                }
1251                for (; i<ndata; i++) {
1252                        yfit[i] += param[0];
1253                        dy[i] = y[i] - yfit[i];
1254                }
1255                break;
1256
1257        case NOISE_POISSON_FIT:
1258                *chisq = 0.0;
1259                // Summation loop over all data
1260                for (i=0; i<fit_start; i++) {
1261                        yfit[i] += param[0];
1262                        dy[i] = y[i] - yfit[i];
1263                }
1264                for ( ; i<fit_end; i++) {
1265                        yfit[i] += param[0];
1266                        dy[i] = y[i] - yfit[i];
1267                        // And find chi^2
1268                        // we still don't let the variance drop below 1
1269                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1270                        *chisq += dy[i] * dy[i] * sig2i;
1271                }
1272                for ( ; i<ndata; i++) {
1273                        yfit[i] += param[0];
1274                        dy[i] = y[i] - yfit[i];
1275                }
1276                break;
1277               
1278        case NOISE_MLE:                              // for the final chisq report a normal chisq measure for MLE
1279                *chisq = 0.0;
1280                // Summation loop over all data
1281                for (i=0; i<fit_start; i++) {
1282                        yfit[i] += param[0];
1283                        dy[i] = y[i] - yfit[i];
1284                }
1285                for ( ; i<fit_end; i++) {
1286                        yfit[i] += param[0];
1287                        dy[i] = y[i] - yfit[i];
1288                        // And find chi^2
1289//                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1290                        if (yfit[i]<=0.0)
1291                                ; // do nothing
1292                        else if (y[i]==0.0)
1293                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1294                        else
1295                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1296                }
1297                for ( ; i<ndata; i++) {
1298                        yfit[i] += param[0];
1299                        dy[i] = y[i] - yfit[i];
1300                }
1301                if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1302                break;
1303
1304        case NOISE_GAUSSIAN_FIT:
1305                *chisq = 0.0;
1306                // Summation loop over all data
1307                for (i=0; i<fit_start; i++) {
1308                        yfit[i] += param[0];
1309                        dy[i] = y[i] - yfit[i];
1310                }
1311                for ( ; i<fit_end; i++) {
1312                        yfit[i] += param[0];
1313                        dy[i] = y[i] - yfit[i];
1314                        // And find chi^2
1315                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1316                        *chisq += dy[i] * dy[i] * sig2i;
1317                }
1318                for ( ; i<ndata; i++) {
1319                        yfit[i] += param[0];
1320                        dy[i] = y[i] - yfit[i];
1321                }
1322                break;
1323       
1324        default:
1325                return -3;
1326                /* break; */   // (unreachable code)
1327        }
1328
1329        return 0;
1330}
1331
1332
1333//********************************* GCI_marquardt_fitting_engine **********************************************************************
1334
1335// This returns the number of iterations or negative if an error occurred.
1336// passes all the data to the ecf routine, checks the returned chisq and re-fits if it is of benefit
1337// was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03
1338
1339int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end, 
1340                                                float prompt[], int nprompt,
1341                                                noise_type noise, float sig[],
1342                                                float param[], int paramfree[],
1343                                           int nparam, restrain_type restrain,
1344                                           void (*fitfunc)(float, float [], float *, float [], int),
1345                                           float *fitted, float *residuals, float *chisq,
1346                                           float **covar, float **alpha, float **erraxes,
1347                                           float chisq_target, float chisq_delta, int chisq_percent)
1348{
1349        float oldChisq, local_chisq;
1350        int ret, tries=0;
1351
1352        if (ecf_exportParams) ecf_ExportParams_OpenFile ();
1353
1354        // All of the work is done by the ECF module
1355        ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1356                                                          prompt, nprompt, noise, sig,
1357                                                          param, paramfree, nparam, restrain, fitfunc,
1358                                                          fitted, residuals, covar, alpha, &local_chisq,
1359                                                          chisq_delta, chisq_percent, erraxes);
1360                                                         
1361        // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same
1362        // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr()
1363        oldChisq = 3.0e38;
1364        while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<(MAXREFITS*3))
1365        {
1366                oldChisq = local_chisq; 
1367                tries++;
1368                ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1369                                                          prompt, nprompt, noise, sig,
1370                                                          param, paramfree, nparam, restrain, fitfunc,
1371                                                          fitted, residuals, covar, alpha, &local_chisq,
1372                                                          chisq_delta, chisq_percent, erraxes);
1373        }                                                 
1374   /*     if (local_chisq<=chisq_target) {
1375            printf("local_chisq %f target %f\n", local_chisq, chisq_target);
1376        }
1377        if (local_chisq>=oldChisq) {
1378            printf("local_chisq %f oldChisq %f\n", local_chisq, oldChisq);
1379        }
1380        if (tries>=(MAXREFITS*3)) {
1381            printf("tries is %d\n", tries);
1382            printf("local_chisq %f oldChisq %f\n", local_chisq, oldChisq);
1383        }
1384
1385        show_timing();*/
1386
1387        if (chisq!=NULL) *chisq = local_chisq;
1388
1389        if (ecf_exportParams) ecf_ExportParams_CloseFile ();
1390
1391        return ret;             // summed number of iterations
1392}
1393
1394/* Cleanup function */
1395void GCI_marquardt_cleanup(void)
1396{
1397/*      if (fnvals_len) {
1398                free(fnvals);
1399                GCI_ecf_free_matrix(dy_dparam_pure);
1400                GCI_ecf_free_matrix(dy_dparam_conv);
1401                fnvals_len = 0;
1402                dy_dparam_nparam_size = 0;
1403        }
1404*/
1405}
1406
1407// Emacs settings:
1408// Local variables:
1409// mode: c
1410// c-basic-offset: 4
1411// tab-width: 4
1412// End:
Note: See TracBrowser for help on using the repository browser.