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

Revision 7778, 56.9 KB checked in by aivar, 8 years ago (diff)

Fixed bug in NOISE_MLE case.

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 and global analysis code in
7   EcfGlobal.c.
8*/
9
10#include <math.h>
11#include <stdio.h>
12#include <stdlib.h>
13#include "EcfInternal.h"
14
15/* Predeclarations */
16int GCI_marquardt_step(float x[], float y[], int ndata,
17                                        noise_type noise, float sig[],
18                                        float param[], int paramfree[], int nparam,
19                                        restrain_type restrain,
20                                        void (*fitfunc)(float, float [], float *, float [], int),
21                                        float yfit[], float dy[],
22                                        float **covar, float **alpha, float *chisq,
23                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam);
24int GCI_marquardt_step_instr(float xincr, float y[],
25                                        int ndata, int fit_start, int fit_end,
26                                        float instr[], int ninstr,
27                                        noise_type noise, float sig[],
28                                        float param[], int paramfree[], int nparam,
29                                        restrain_type restrain,
30                                        void (*fitfunc)(float, float [], float *, float [], int),
31                                        float yfit[], float dy[],
32                                        float **covar, float **alpha, float *chisq,
33                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,
34                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
35                                        int *pfnvals_len, int *pdy_dparam_nparam_size);
36int GCI_marquardt_estimate_errors(float **alpha, int nparam, int mfit,
37                                                                  float d[], float **v, float interval);
38
39/* Globals */
40//static float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
41//static int fnvals_len=0, dy_dparam_nparam_size=0;
42// was Global, now thread safe
43
44/********************************************************************
45
46                                           SINGLE TRANSIENT FITTING
47
48                                                TRIPLE INTEGRAL METHOD
49
50 ********************************************************************/
51
52/* Start with an easy one: the three integral method.  This returns 0
53   on success, negative on error. */
54
55int GCI_triple_integral(float xincr, float y[],
56                                                int fit_start, int fit_end,
57                                                noise_type noise, float sig[],
58                                                float *Z, float *A, float *tau,
59                                                float *fitted, float *residuals,
60                                                float *chisq, int division)
61{
62        float d1, d2, d3, d12, d23;
63        float t0, dt, exp_dt_tau, exp_t0_tau;
64        int width;
65        int i;
66        float sigma2, res, chisq_local;
67
68        width = (fit_end - fit_start) / division;
69        if (width <= 0)
70                return -1;
71
72        t0 = fit_start * xincr;
73        dt = width * xincr;
74
75        d1 = d2 = d3 = 0;
76        for (i=fit_start; i<fit_start+width; i++)           d1 += y[i];
77        for (i=fit_start+width; i<fit_start+2*width; i++)   d2 += y[i];
78        for (i=fit_start+2*width; i<fit_start+3*width; i++) d3 += y[i];
79
80        // Those are raw sums, we now convert into areas */
81        d1 *= xincr;
82        d2 *= xincr;
83        d3 *= xincr;
84
85        d12 = d1 - d2;
86        d23 = d2 - d3;
87        if (d12 <= d23 || d23 <= 0)
88                return -2;
89
90        exp_dt_tau = d23 / d12;  /* exp(-dt/tau) */
91        *tau = -dt / log(exp_dt_tau);
92        exp_t0_tau = exp(-t0/(*tau));
93        *A = d12 / ((*tau) * exp_t0_tau * (1 - 2*exp_dt_tau + exp_dt_tau*exp_dt_tau));
94        *Z = (d1 - (*A) * (*tau) * exp_t0_tau * (1 - exp_dt_tau)) / dt;
95
96        /* Now calculate the fitted curve and chi-squared if wanted. */
97        if (fitted == NULL)
98                return 0;
99
100        for (i=0; i<fit_end; i++)
101                fitted[i] = (*Z) + (*A) * exp(-i*xincr/(*tau));
102
103        // OK, so now fitted contains our data for the timeslice of interest.
104        // We can calculate a chisq value and plot the graph, along with
105        // the residuals.
106
107        if (residuals == NULL && chisq == NULL)
108                return 0;
109
110        chisq_local = 0.0;
111        for (i=0; i<fit_start; i++) {
112                res = y[i]-fitted[i];
113                if (residuals != NULL)
114                        residuals[i] = res;
115        }
116
117        switch (noise) {
118        case NOISE_CONST:
119                /* Summation loop over all data */
120                for ( ; i<fit_end; i++) {
121                        res = y[i] - fitted[i];
122                        if (residuals != NULL)
123                                residuals[i] = res;
124                        chisq_local += res * res;
125                }
126
127                chisq_local /= (sig[0]*sig[0]);
128                break;
129
130        case NOISE_GIVEN:
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                        chisq_local += (res * res) / (sig[i] * sig[i]);
137                }
138                break;
139
140        case NOISE_POISSON_DATA:
141                /* Summation loop over all data */
142                for ( ; i<fit_end; i++) {
143                        res = y[i] - fitted[i];
144                        if (residuals != NULL)
145                                residuals[i] = res;
146                        /* don't let variance drop below 1 */
147                        sigma2 = (y[i] > 1 ? 1.0/y[i] : 1.0);
148                        chisq_local += res * res * sigma2;
149                }
150                break;
151
152        case NOISE_POISSON_FIT:
153        case NOISE_GAUSSIAN_FIT:   //  NOISE_GAUSSIAN_FIT and NOISE_MLE not implemented for triple integral
154        case NOISE_MLE:
155                /* Summation loop over all data */
156                for ( ; i<fit_end; i++) {
157                        res = y[i] - fitted[i];
158                        if (residuals != NULL)
159                                residuals[i] = res;
160                        /* don't let variance drop below 1 */
161                        sigma2 = (fitted[i] > 1 ? 1.0/fitted[i] : 1.0);
162                        chisq_local += res * res * sigma2;
163                }
164                break;
165
166        default:
167                return -3;
168                /* break; */   // (unreachable code)
169        }
170
171        if (chisq != NULL)
172                *chisq = chisq_local;
173
174        return 0;
175}
176
177
178int GCI_triple_integral_instr(float xincr, float y[],
179                                                          int fit_start, int fit_end,
180                                                          float instr[], int ninstr,
181                                                          noise_type noise, float sig[],
182                                                          float *Z, float *A, float *tau,
183                                                          float *fitted, float *residuals,
184                                                          float *chisq, int division)
185{
186        float d1, d2, d3, d12, d23;
187        float t0, dt, exp_dt_tau, exp_t0_tau;
188        int width;
189        int i, j;
190        float sigma2, res, chisq_local;
191        float sum, scaling;
192        int fitted_preconv_size=0;   // was static but now thread safe
193        float *fitted_preconv;   // was static but now thread safe
194
195        width = (fit_end - fit_start) / division;
196        if (width <= 0)
197                return -1;
198
199        t0 = fit_start * xincr;
200        dt = width * xincr;
201
202        d1 = d2 = d3 = 0;
203        for (i=fit_start; i<fit_start+width; i++)
204                d1 += y[i];
205        for (i=fit_start+width; i<fit_start+2*width; i++)
206                d2 += y[i];
207        for (i=fit_start+2*width; i<fit_start+3*width; i++)
208                d3 += y[i];
209
210        // Those are raw sums, we now convert into areas */
211        d1 *= xincr;
212        d2 *= xincr;
213        d3 *= xincr;
214
215        d12 = d1 - d2;
216        d23 = d2 - d3;
217        if (d12 <= d23 || d23 <= 0)
218                return -2;
219
220        exp_dt_tau = d23 / d12;  /* exp(-dt/tau) */
221        *tau = -dt / log(exp_dt_tau);
222        exp_t0_tau = exp(-t0/(*tau));
223        *A = d12 /
224                  ((*tau) * exp_t0_tau * (1 - 2*exp_dt_tau + exp_dt_tau*exp_dt_tau));
225        *Z = (d1 - (*A) * (*tau) * exp_t0_tau * (1 - exp_dt_tau)) / dt;
226
227        /* We now convolve with the instrument response to hopefully get a
228           slightly better fit.  We'll also scale by an appropriate
229           scale factor, which turns out to be:
230
231                  sum_{i=0}^{ninstr-1} instr[i]*exp(i*xincr/tau)
232
233           which should be only a little greater than the sum of the
234           instrument response values.
235        */
236
237        sum = scaling = 0;
238        for (i=0; i<ninstr; i++) {
239                sum += instr[i];
240                scaling += instr[i] * exp(i*xincr/(*tau));
241        }
242
243        scaling /= sum;  /* Make instrument response sum to 1.0 */
244        (*A) /= scaling;
245
246        /* Now calculate the fitted curve and chi-squared if wanted. */
247        if (fitted == NULL)
248                return 0;
249
250//      if (fitted_preconv_size < fit_end) {
251//              if (fitted_preconv_size > 0)
252//                      free(fitted_preconv);
253                if ((fitted_preconv = (float *) malloc(fit_end * sizeof(float)))
254                        == NULL)
255                        return -3;
256                else
257                        fitted_preconv_size = fit_end;
258//      }
259
260        for (i=0; i<fit_end; i++)
261                fitted_preconv[i] = (*A) * exp(-i*xincr/(*tau));
262
263        for (i=0; i<fit_end; i++) {
264                int convpts;
265
266                /* (Zero-basing everything in sight...)
267                   We wish to find fitted = fitted_preconv * instr, so explicitly:
268                          fitted[i] = sum_{j=0}^i fitted_preconv[i-j].instr[j]
269                   But instr[k]=0 for k >= ninstr, so we only need to sum:
270                          fitted[i] = sum_{j=0}^{min(ninstr-1,i)}
271                                              fitted_preconv[i-j].instr[j]
272                */
273
274                fitted[i] = 0;
275                convpts = (ninstr <= i) ? ninstr-1 : i;
276                for (j=0; j<=convpts; j++) {
277                        fitted[i] += fitted_preconv[i-j]*instr[j];
278                }
279
280                fitted[i] += *Z;
281        }
282
283        // OK, so now fitted contains our data for the timeslice of interest.
284        // We can calculate a chisq value and plot the graph, along with
285        // the residuals.
286
287        if (residuals == NULL && chisq == NULL)
288                return 0;
289
290        chisq_local = 0.0;
291        for (i=0; i<fit_start; i++) {
292                res = y[i]-fitted[i];
293                if (residuals != NULL)
294                        residuals[i] = res;
295        }
296
297        switch (noise) {
298        case NOISE_CONST:
299                /* Summation loop over all data */
300                for ( ; i<fit_end; i++) {
301                        res = y[i] - fitted[i];
302                        if (residuals != NULL)
303                                residuals[i] = res;
304                        chisq_local += res * res;
305                }
306
307                chisq_local /= (sig[0]*sig[0]);
308                break;
309
310        case NOISE_GIVEN:
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                        chisq_local += (res * res) / (sig[i] * sig[i]);
317                }
318                break;
319
320        case NOISE_POISSON_DATA:
321                /* Summation loop over all data */
322                for ( ; i<fit_end; i++) {
323                        res = y[i] - fitted[i];
324                        if (residuals != NULL)
325                                residuals[i] = res;
326                        /* don't let variance drop below 1 */
327                        sigma2 = (y[i] > 1 ? 1.0/y[i] : 1.0);
328                        chisq_local += res * res * sigma2;
329                }
330                break;
331
332        case NOISE_POISSON_FIT:
333        case NOISE_GAUSSIAN_FIT:   //  NOISE_GAUSSIAN_FIT and NOISE_MLE not implemented for triple integral
334        case NOISE_MLE:
335                /* Summation loop over all data */
336                for ( ; i<fit_end; i++) {
337                        res = y[i] - fitted[i];
338                        if (residuals != NULL)
339                                residuals[i] = res;
340                        /* don't let variance drop below 1 */
341                        sigma2 = (fitted[i] > 1 ? 1.0/fitted[i] : 1.0);
342                        chisq_local += res * res * sigma2;
343                }
344                break;
345
346        default:
347                return -3;
348                /* break; */   // (unreachable code)
349        }
350
351        if (chisq != NULL)
352                *chisq = chisq_local;
353
354        return 0;
355}
356
357int GCI_triple_integral_fitting_engine(float xincr, float y[], int fit_start, int fit_end,
358                                                          float instr[], int ninstr, noise_type noise, float sig[],
359                                                          float *Z, float *A, float *tau, float *fitted, float *residuals,
360                                                          float *chisq, float chisq_target)
361{
362        int tries=1, division=3;                 // the data
363        float local_chisq=3.0e38, oldChisq=3.0e38, oldZ, oldA, oldTau, *validFittedArray; // local_chisq a very high float but below oldChisq
364
365        if (fitted==NULL)   // we require chisq but have not supplied a "fitted" array so must malloc one
366        {
367                if ((validFittedArray = malloc(fit_end * sizeof(float)))== NULL) return (-1);
368        }
369        else validFittedArray = fitted;
370
371        if (instr==NULL)           // no instrument/prompt has been supplied
372        {
373                GCI_triple_integral(xincr, y, fit_start, fit_end, noise, sig,
374                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
375
376                while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS)
377                {
378                        oldChisq = local_chisq;
379                        oldZ = *Z;
380                        oldA = *A;
381                        oldTau = *tau;
382//                      division++;
383                        division+=division/3;
384                        tries++;
385                        GCI_triple_integral(xincr, y, fit_start, fit_end, noise, sig,
386                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
387                }
388        }
389        else
390        {
391                GCI_triple_integral_instr(xincr, y, fit_start, fit_end, instr, ninstr, noise, sig,
392                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
393
394                while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS)
395                {
396                        oldChisq = local_chisq;
397                        oldZ = *Z;
398                        oldA = *A;
399                        oldTau = *tau;
400//                      division++;
401                        division+=division/3;
402                        tries++;
403                        GCI_triple_integral_instr(xincr, y, fit_start, fit_end, instr, ninstr, noise, sig,
404                                                                Z, A, tau, validFittedArray, residuals, &local_chisq, division);
405
406                }
407        }
408
409        if (local_chisq>oldChisq)          // the previous fit was better
410        {
411                local_chisq = oldChisq;
412                *Z = oldZ;
413                *A = oldA;
414                *tau = oldTau;
415        }
416
417        if (chisq!=NULL) *chisq = local_chisq;
418
419        if (fitted==NULL)
420        {
421                free (validFittedArray);
422        }
423
424        return(tries);
425}
426
427/********************************************************************
428
429                                           SINGLE TRANSIENT FITTING
430
431                                          LEVENBERG-MARQUARDT METHOD
432
433 ********************************************************************/
434
435/* Now for the non-linear least squares fitting algorithms.
436
437   The process is:
438   - for Gaussian noise, use Levenberg-Marquardt directly
439   - for Poisson noise, use Levenberg-Marquardt to get an initial
440   estimate of the parameters assuming constant error variance, then
441   use amoeba to improve the estimate, assuming that the error
442   variance is proportional to the function value (with a minimum
443   variance of 15 to handle the case when the Poisson distribution
444   is not approximately Gaussian, so that the very noisy tails do
445   not inappropriately weight the solution).
446
447
448   This code contains two variants of the Levenberg-Marquardt method
449   for slightly different situations.  This attempts to reduce the
450   value chi^2 of a fit between a set of data points x[0..ndata-1],
451   y[0..ndata-1] and a nonlinear function dependent on nparam
452   coefficients param[0..nparam-1].  In the case that the x values are
453   equally spaced and start at zero, we can also handle convolution
454   with an instrument response instr[0..ninstr-1] and only look at the
455   data points from fit_start..fit_end-1.  The first variant does not
456   handle an instrument response and takes any values of
457   x[0..ndata-1].  The second variant takes an xincr and will handle
458   an instrument response if ninstr > 0.  The individual standard
459   deviations of the errors are determined by the value of noise: if
460   noise=NOISE_CONST, the standard deviations are constant, given by
461   sig[0]=*sig, if noise=NOISE_GIVEN, the standard deviations are
462   given by sig[0..ndata-1], if noise=NOISE_POISSON_DATA, the standard
463   deviations are taken to be given by Poisson noise, and the
464   variances are taken to be max(y[i],15), and if
465   noise=NOISE_POISSON_FIT, the variances are taken to be
466   max(yfit[i],15). If noise=NOISE_GAUSSIAN_FIT, the variances are taken to be
467   yfit[i] and if noise=NOISE_MLE then a optimisation is for the maximum
468   likelihood approximation (Laurence and Chromy in press 2010 and
469   G. Nishimura, and M. Tamura Phys Med Biol 2005).
470
471   The input array paramfree[0..nparam-1] indicates
472   by nonzero entries those components that should be held fixed at
473   their input values. The program returns current best-fit values for
474   the parameters param[0..nparam-1] and chi^2 = chisq.  The arrays
475   covar[0..nparam-1][0..nparam-1] and alpha[0..nparam-1][0..nparam-1]
476   are used as working space during most isterations.  Supply a
477   routine fitfunc(x,param,yfit,dy_dparam,nparam) that evaluates the
478   fitting function fitfunc and its derivatives dydy[1..nparam-1] with
479   respect to the fitting parameters param at x.  (See below for
480   information about zero offsets, though.)  The values of fitfunc,
481   modified by the instrument response, are returned in the array yfit
482   and the differences y - yfit in dy.  The first call _must_ provide
483   an initial guess for the parameters param and set alambda < 0 for
484   initialisation (which then sets alambda = 0.001).  If a step
485   succeeds, chisq becomes smaller and alambda decreases by a factor
486   of 10.  You must call this routine repeatedly until convergence is
487   achieved. Then make one final call with alambda=0 to perform
488   cleanups and so that covar[0..nparam-1][0..nparam-1] returns the
489   covariance matrix and alpha the curvature matrix.  (Parameters held
490   fixed will return zero covariances.)
491
492   One key extra piece which is particularly important in the
493   instrument response case.  The parameter param[0] is assumed to be
494   the zero offset of the signal, which applies before and after time
495   zero.  It thus simply contributes param[0]*sum(instr) to the signal
496   value rather than being convolved with the instrument response only
497   from time zero.  For this reason, the fitfunc should ignore
498   param[0], as the fitting routines will handle this offset.
499*/
500
501/* These two functions do the whole job */
502int GCI_marquardt(float x[], float y[], int ndata,
503                                  noise_type noise, float sig[],
504                                  float param[], int paramfree[], int nparam,
505                                  restrain_type restrain,
506                                  void (*fitfunc)(float, float [], float *, float [], int),
507                                  float *fitted, float *residuals,
508                                  float **covar, float **alpha, float *chisq,
509                                  float chisq_delta, float chisq_percent, float **erraxes)
510{
511        float alambda, ochisq;
512        int mfit;
513        float evals[MAXFIT];
514        int i, k, itst, itst_max;
515
516        float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];
517
518        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
519
520        mfit = 0;
521        for (i=0; i<nparam; i++) {
522                if (paramfree[i]) mfit++;
523        }
524
525        alambda = -1;
526        if (GCI_marquardt_step(x, y, ndata, noise, sig,
527                                                   param, paramfree, nparam, restrain,
528                                                   fitfunc, fitted, residuals,
529                                                   covar, alpha, chisq, &alambda,
530                                                   &mfit, &ochisq2, paramtry, beta, dparam) != 0) {
531                return -1;
532        }
533
534        k = 1;  /* Iteration counter */
535        itst = 0;
536        for (;;) {
537                k++;
538                if (k > MAXITERS) {
539                        return -2;
540                }
541
542                ochisq = *chisq;
543                if (GCI_marquardt_step(x, y, ndata, noise, sig,
544                                                           param, paramfree, nparam, restrain,
545                                                           fitfunc, fitted, residuals,
546                                                           covar, alpha, chisq, &alambda,
547                                                           &mfit, &ochisq2, paramtry, beta, dparam) != 0) {
548                        return -3;
549                }
550
551                if (*chisq > ochisq)
552                        itst = 0;
553                else if (ochisq - *chisq < chisq_delta)
554                        itst++;
555
556                if (itst < itst_max) continue;
557
558                /* Endgame */
559                alambda = 0.0;
560                if (GCI_marquardt_step(x, y, ndata, noise, sig,
561                                                           param, paramfree, nparam, restrain,
562                                                           fitfunc, fitted, residuals,
563                                                           covar, alpha, chisq, &alambda,
564                                                           &mfit, &ochisq2, paramtry, beta, dparam) != 0) {
565                        return -4;
566                }
567
568                if (erraxes == NULL)
569                        return k;
570
571                //TODO ARG
572                //if (GCI_marquardt_estimate_errors(alpha, nparam, mfit, evals,
573                //                                                                erraxes, chisq_percent) != 0) {
574                //      return -5;
575                //}
576
577                break;  /* We're done now */
578        }
579
580        return k;
581}
582
583#define do_frees \
584        if (fnvals) free(fnvals);\
585        if (dy_dparam_pure) GCI_ecf_free_matrix(dy_dparam_pure);\
586        if (dy_dparam_conv) GCI_ecf_free_matrix(dy_dparam_conv);
587
588int GCI_marquardt_instr(float xincr, float y[],
589                                        int ndata, int fit_start, int fit_end,
590                                        float instr[], int ninstr,
591                                        noise_type noise, float sig[],
592                                        float param[], int paramfree[], int nparam,
593                                        restrain_type restrain,
594                                        void (*fitfunc)(float, float [], float *, float [], int),
595                                        float *fitted, float *residuals,
596                                        float **covar, float **alpha, float *chisq,
597                                        float chisq_delta, float chisq_percent, float **erraxes)
598{
599        float alambda, ochisq;
600        int mfit, mfit2;
601        float evals[MAXFIT];
602        int i, k, itst, itst_max;
603
604        // The following are declared here to retain some optimisation by not repeatedly mallocing
605        // (only once per transient), but to remain thread safe.
606        // They are malloced by lower fns but at the end, freed by this fn.
607        // These vars were global or static before thread safety was introduced.
608        float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL;
609        int fnvals_len=0, dy_dparam_nparam_size=0;
610        float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];
611
612        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
613
614        mfit = 0;
615        for (i=0; i<nparam; i++) {
616                if (paramfree[i]) mfit++;
617        }
618
619        if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
620
621        alambda = -1;
622        if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
623                                                                 instr, ninstr, noise, sig,
624                                                                 param, paramfree, nparam, restrain,
625                                                                 fitfunc, fitted, residuals,
626                                                                 covar, alpha, chisq, &alambda,
627                                                                 &mfit2, &ochisq2, paramtry, beta, dparam,
628                                                                 &fnvals, &dy_dparam_pure, &dy_dparam_conv,
629                                                                 &fnvals_len, &dy_dparam_nparam_size) != 0) {
630                do_frees
631                return -1;
632        }
633
634        if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
635
636        k = 1;  /* Iteration counter */
637        itst = 0;
638        for (;;) {
639                k++;
640                if (k > MAXITERS) {
641                        do_frees
642                        return -2;
643                }
644
645                ochisq = *chisq;
646                if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
647                                                                         instr, ninstr, noise, sig,
648                                                                         param, paramfree, nparam, restrain,
649                                                                         fitfunc, fitted, residuals,
650                                                                         covar, alpha, chisq, &alambda,
651                                                                         &mfit2, &ochisq2, paramtry, beta, dparam,
652                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv,
653                                                                         &fnvals_len, &dy_dparam_nparam_size) != 0) {
654                        do_frees
655                        return -3;
656                }
657
658                if (ecf_exportParams) ecf_ExportParams (param, nparam, *chisq);
659
660                if (*chisq > ochisq)
661                        itst = 0;
662                else if (ochisq - *chisq < chisq_delta)
663                        itst++;
664
665                if (itst < itst_max) continue;
666
667                /* Endgame */
668                alambda = 0.0;
669                if (GCI_marquardt_step_instr(xincr, y, ndata, fit_start, fit_end,
670                                                                         instr, ninstr, noise, sig,
671                                                                         param, paramfree, nparam, restrain,
672                                                                         fitfunc, fitted, residuals,
673                                                                         covar, alpha, chisq, &alambda,
674                                                                         &mfit2, &ochisq2, paramtry, beta, dparam,
675                                                                         &fnvals, &dy_dparam_pure, &dy_dparam_conv,
676                                                                         &fnvals_len, &dy_dparam_nparam_size) != 0) {
677                        do_frees
678                        return -4;
679                }
680
681                if (erraxes == NULL){
682                        do_frees
683                        return k;
684                }
685
686//TODO ARG this estimate errors call was deleted in my latest version
687        //      if (GCI_marquardt_estimate_errors(alpha, nparam, mfit, evals,
688        //                                                                        erraxes, chisq_percent) != 0) {
689        //              do_frees
690        //              return -5;
691        //      }
692
693                break;  /* We're done now */
694        }
695
696        do_frees
697        return k;
698}
699#undef do_frees
700
701
702int GCI_marquardt_step(float x[], float y[], int ndata,
703                                        noise_type noise, float sig[],
704                                        float param[], int paramfree[], int nparam,
705                                        restrain_type restrain,
706                                        void (*fitfunc)(float, float [], float *, float [], int),
707                                        float yfit[], float dy[],
708                                        float **covar, float **alpha, float *chisq,
709                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam)
710{
711        int j, k, l, ret;
712//      static int mfit;   // was static but now thread safe
713//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe
714        int mfit = *pmfit;
715        float ochisq = *pochisq;
716
717        if (nparam > MAXFIT)
718                return -1;
719
720        /* Initialisation */
721        /* We assume we're given sensible starting values for param[] */
722        if (*alambda < 0.0) {
723                for (mfit=0, j=0; j<nparam; j++)
724                        if (paramfree[j])
725                                mfit++;
726
727                if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig,
728                                                                         param, paramfree, nparam, fitfunc,
729                                                                         yfit, dy,
730                                                                         alpha, beta, chisq, 0.0, *alambda) != 0)
731                        return -2;
732
733                *alambda = 0.001;
734                ochisq = (*chisq);
735                for (j=0; j<nparam; j++)
736                        paramtry[j] = param[j];
737        }
738
739        /* Alter linearised fitting matrix by augmenting diagonal elements */
740        for (j=0; j<mfit; j++) {
741                for (k=0; k<mfit; k++)
742                        covar[j][k] = alpha[j][k];
743
744                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
745                dparam[j] = beta[j];
746        }
747
748        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */
749        if (GCI_solve(covar, mfit, dparam) != 0)
750                return -1;
751
752        /* Once converged, evaluate covariance matrix */
753        if (*alambda == 0) {
754                if (GCI_marquardt_compute_fn_final(x, y, ndata, noise, sig,
755                        param, paramfree, nparam, fitfunc,
756                        yfit, dy, chisq) != 0)
757                    return -3;
758
759
760                for (j=0; j<mfit; j++)
761                    for (k=0; k<mfit; k++)
762                        covar[j][k] = alpha[j][k];
763                GCI_invert(covar, mfit);
764
765                if (mfit < nparam) {  /* no need to do this otherwise */
766                        GCI_covar_sort(covar, nparam, paramfree, mfit);
767                        GCI_covar_sort(alpha, nparam, paramfree, mfit);
768                }
769                return 0;
770        }
771
772        /* Did the trial succeed? */
773        for (j=0, l=0; l<nparam; l++)
774                if (paramfree[l])
775                        paramtry[l] = param[l] + dparam[j++];
776
777        if (restrain == ECF_RESTRAIN_DEFAULT)
778                ret = check_ecf_params (paramtry, nparam, fitfunc);
779        else
780                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
781
782        if (ret != 0) {
783                /* Bad parameters, increase alambda and return */
784                *alambda *= 10.0;
785                return 0;
786        }
787
788        if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig,
789                                                                 paramtry, paramfree, nparam, fitfunc,
790                                                                 yfit, dy, covar, dparam,
791                                                                 chisq, ochisq, *alambda) != 0)
792                return -2;
793
794        /* Success, accept the new solution */
795        if (*chisq < ochisq) {
796                *alambda *= 0.1;
797                ochisq = *chisq;
798                for (j=0; j<mfit; j++) {
799                        for (k=0; k<mfit; k++)
800                                alpha[j][k] = covar[j][k];
801                        beta[j] = dparam[j];
802                }
803                for (l=0; l<nparam; l++)
804                        param[l] = paramtry[l];
805        } else { /* Failure, increase alambda and return */
806                *alambda *= 10.0;
807                *chisq = ochisq;
808        }
809
810        return 0;
811}
812
813
814int GCI_marquardt_step_instr(float xincr, float y[],
815                                        int ndata, int fit_start, int fit_end,
816                                        float instr[], int ninstr,
817                                        noise_type noise, float sig[],
818                                        float param[], int paramfree[], int nparam,
819                                        restrain_type restrain,
820                                        void (*fitfunc)(float, float [], float *, float [], int),
821                                        float yfit[], float dy[],
822                                        float **covar, float **alpha, float *chisq,
823                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,
824                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
825                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
826{
827        int j, k, l, ret;
828//      static int mfit;   // was static but now thread safe
829//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe
830        //TODO ARG GCI_marquardt_step defines and uses int mfit = *pmfit;
831
832        if (nparam > MAXFIT)
833                return -10;
834        if (xincr <= 0)
835                return -11;
836        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
837                return -12;
838
839        /* Initialisation */
840        /* We assume we're given sensible starting values for param[] */
841        if (*alambda < 0.0) {
842
843                for ((*pmfit)=0, j=0; j<nparam; j++)
844                        if (paramfree[j])
845                                (*pmfit)++;
846
847                if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
848                                                                                   instr, ninstr, noise, sig,
849                                                                                   param, paramfree, nparam, fitfunc,
850                                                                                   yfit, dy, alpha, beta, chisq, 0.0,
851                                                                                   *alambda,
852                                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
853                                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0)
854                        return -2;
855
856                *alambda = 0.001;
857                *pochisq = *chisq;
858                for (j=0; j<nparam; j++)
859                        paramtry[j] = param[j];
860
861        }
862
863        /* Alter linearised fitting matrix by augmenting diagonal elements */
864        for (j=0; j<(*pmfit); j++) {
865                for (k=0; k<(*pmfit); k++)
866                        covar[j][k] = alpha[j][k];
867
868                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
869                dparam[j] = beta[j];
870        }
871
872        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */
873        if (GCI_solve(covar, *pmfit, dparam) != 0)
874                return -1;
875
876        /* Once converged, evaluate covariance matrix */
877        if (*alambda == 0) {
878                if (GCI_marquardt_compute_fn_final_instr(
879                        xincr, y, ndata, fit_start, fit_end,
880                        instr, ninstr, noise, sig,
881                        param, paramfree, nparam, fitfunc,
882                        yfit, dy, chisq,
883                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
884                        pfnvals_len, pdy_dparam_nparam_size) != 0)
885                    return -3;
886
887                for (j=0; j<(*pmfit); j++)
888                    for (k=0; k<(*pmfit); k++)
889                        covar[j][k] = alpha[j][k];
890                GCI_invert(covar, (*pmfit));
891
892
893                if (*pmfit < nparam) {  /* no need to do this otherwise */
894                        GCI_covar_sort(covar, nparam, paramfree, *pmfit);
895                        GCI_covar_sort(alpha, nparam, paramfree, *pmfit);
896                }
897                return 0;
898        }
899
900        /* Did the trial succeed? */
901        for (j=0, l=0; l<nparam; l++)
902                if (paramfree[l])
903                        paramtry[l] = param[l] + dparam[j++];
904
905        if (restrain == ECF_RESTRAIN_DEFAULT)
906                ret = check_ecf_params (paramtry, nparam, fitfunc);
907        else
908                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
909
910        if (ret != 0) {
911                /* Bad parameters, increase alambda and return */
912                *alambda *= 10.0;
913                return 0;
914        }
915
916        if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
917                                                                           instr, ninstr, noise, sig,
918                                                                           paramtry, paramfree, nparam, fitfunc,
919                                                                           yfit, dy, covar, dparam,
920                                                                           chisq, *pochisq, *alambda,
921                                                                           pfnvals, pdy_dparam_pure, pdy_dparam_conv,
922                                                                           pfnvals_len, pdy_dparam_nparam_size) != 0)
923                return -2;
924
925        /* Success, accept the new solution */
926        if (*chisq < *pochisq) {
927                *alambda *= 0.1;
928                *pochisq = *chisq;
929                for (j=0; j<(*pmfit); j++) {
930                        for (k=0; k<(*pmfit); k++)
931                                alpha[j][k] = covar[j][k];
932                        beta[j] = dparam[j];
933                }
934                for (l=0; l<nparam; l++)
935                        param[l] = paramtry[l];
936        } else { /* Failure, increase alambda and return */
937                *alambda *= 10.0;
938                *chisq = *pochisq;
939        }
940
941        return 0;
942}
943
944
945/* Used by GCI_marquardt to evaluate the linearised fitting matrix alpha
946   and vector beta and to calculate chi^2.      The equations involved are
947   given in section 15.5 of Numerical Recipes; basically:
948
949     \chi^2(param) = \sum_{i=1}^N ( (y_i-y(x_i;param)) / sigma_i )^2
950
951     beta_k = -1/2 (d/dparam_k)(chi^2)
952
953     alpha_kl = \sum_{i=1}^N (1/sigma_i^2) .
954                  (dy(x_i;param)/dparam_k) . (dy(x_i;param)/dparam_l)
955
956   (where all of the derivatives are partial).
957
958   If an instrument response is provided, we also take account of it
959   now.  We are given that:
960
961     observed(t) = response(t) * instr(t)
962
963   where response(t) is being fitted with fitfunc; it is also trivial to
964   show that (assuming that instr(t) is known and fixed, with no
965   dependencies on the param_k, the parameters being fitted):
966
967     (d/dparam_k) observed(t) = ((d/dparam_k) response(t)) * instr(t)
968
969   so we do not need to alter the response function in any way to
970   determined the fitted convolved response.
971
972   Again there are two variants of this function, corresponding to the
973   two variants of the Marquardt function.
974*/
975int GCI_marquardt_compute_fn(float x[], float y[], int ndata,
976                                         noise_type noise, float sig[],
977                                         float param[], int paramfree[], int nparam,
978                                         void (*fitfunc)(float, float [], float *, float [], int),
979                                         float yfit[], float dy[],
980                                         float **alpha, float beta[], float *chisq, float old_chisq,
981                                         float alambda)
982{
983        int i, j, k, l, m, mfit;
984        float wt, sig2i, y_ymod, dy_dparam[MAXBINS][MAXFIT];
985        float alpha_weight[MAXBINS];
986        float beta_weight[MAXBINS];
987        int q;
988        float weight;
989        int i_free;
990        int j_free;
991        float dot_product;
992        float beta_sum;
993        float dy_dparam_k_i;
994
995        for (j=0, mfit=0; j<nparam; j++)
996                if (paramfree[j])
997                        mfit++;
998
999        *chisq = 0.0f;
1000
1001        switch (noise) {
1002                case NOISE_CONST:
1003                {
1004                        float i_sig_0_squared = 1.0 / (sig[0] * sig[0]);
1005                        for (q = 0; q < ndata; ++q) {
1006                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1007                                yfit[q] += param[0];
1008                                dy_dparam[q][0] = 1.0;
1009                                dy[q] = y[q] - yfit[q];
1010                                weight = i_sig_0_squared;
1011                                alpha_weight[q] = weight; // 1 / (sig[0] * sig[0])
1012                                weight *= dy[q];
1013                                beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0])
1014                                weight *= dy[q];
1015                                *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0])
1016                        }
1017                        break;
1018                }
1019                case NOISE_GIVEN:
1020                {
1021                        for (q = 0; q < ndata; ++q) {
1022                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1023                                yfit[q] += param[0];
1024                                dy_dparam[q][0] = 1.0;
1025                                dy[q] = y[q] - yfit[q];
1026                                weight = 1.0f / (sig[q] * sig[q]);
1027                                alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
1028                                weight *= dy[q];
1029                                beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
1030                                weight *= dy[q];
1031                                *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
1032                        }
1033                        break;
1034                }
1035                case NOISE_POISSON_DATA:
1036                {
1037                        for (q = 0; q < ndata; ++q) {
1038                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1039                                yfit[q] += param[0];
1040                                dy_dparam[q][0] = 1.0;
1041                                dy[q] = y[q] - yfit[q];
1042                                weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
1043                                alpha_weight[q] = weight; // 1 / sig(q)
1044                                weight *= dy[q];
1045                                beta_weight[q] = weight; // dy[q] / sig(q)
1046                                weight *= dy[q];
1047                                *chisq += weight; // (dy[q] * dy[q]) / sig(q)
1048                        }
1049                        break;
1050                }
1051                case NOISE_POISSON_FIT:
1052                {
1053                        for (q = 0; q < ndata; ++q) {
1054                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1055                                yfit[q] += param[0];
1056                                dy_dparam[q][0] = 1.0;
1057                                dy[q] = y[q] - yfit[q];
1058                                weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1059                                alpha_weight[q] = weight; // 1 / sig(q)
1060                                weight *= dy[q];
1061                                beta_weight[q] = weight; // dy(q) / sig(q)
1062                                weight *= dy[q];
1063                                *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1064                        }
1065                        break;
1066                }
1067                case NOISE_GAUSSIAN_FIT:
1068                {
1069                        for (q = 0; q < ndata; ++q) {
1070                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1071                                yfit[q] += param[0];
1072                                dy_dparam[q][0] = 1.0;
1073                                dy[q] = y[q] - yfit[q];
1074                                weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1075                                alpha_weight[q] = weight; // 1 / sig(q)
1076                                weight *= dy[q];
1077                                beta_weight[q] = weight; // dy[q] / sig(q)
1078                                weight *= dy[q];
1079                                *chisq += weight; // dy[q] / sig(q)
1080                        }
1081                        break;
1082                }
1083                case NOISE_MLE:
1084                {
1085                        for (q = 0; q < ndata; ++q) {
1086                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1087                                yfit[q] += param[0];
1088                                dy_dparam[q][0] = 1.0;
1089                                dy[q] = y[q] - yfit[q];
1090                                weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f);
1091                                alpha_weight[q] = weight * y[q] / yfit[q];
1092                                beta_weight[q] = dy[q] * weight;
1093                                if (yfit[q] > 0.0) {
1094                                        *chisq += (0.0f == y[q])
1095                                                        ? 2.0 * yfit[q]
1096                                                        : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1097                                }
1098                        }
1099                        if (*chisq <= 0.0f) {
1100                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1101                        }
1102                        break;
1103                }
1104                default:
1105                {
1106                        return -3;
1107                }
1108        }
1109
1110        // Check if chi square worsened:
1111        if (0.0f != old_chisq && *chisq >= old_chisq) {
1112                // don't bother to set up the matrices for solution
1113                return 0;
1114        }
1115
1116        i_free = 0;
1117        // for all columns
1118        for (i = 0; i < nparam; ++i) {
1119                if (paramfree[i]) {
1120                        j_free = 0;
1121                        beta_sum = 0.0f;
1122                        // row loop, only need to consider lower triangle
1123                        for (j = 0; j <= i; ++j) {
1124                                if (paramfree[j]) {
1125                                        dot_product = 0.0f;
1126                                        if (0 == j_free) { // true only once for each outer loop i
1127                                                // for all data
1128                                                for (k = 0; k < ndata; ++k) {
1129                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k];
1130                                                        beta_sum += dy_dparam[k][i] * beta_weight[k]; // compute beta only once for each i
1131                                                }
1132                                        }
1133                                        else {
1134                                                // for all data
1135                                                for (k = 0; k < ndata; ++k) {
1136                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k];
1137                                                }
1138                                        } // k loop
1139                                               
1140                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; //TODO ARG w/n/b [i][j] more common usage?  row/column
1141                                        // if (i_free != j_free) { //TODO ARG this approach seemed slower at one time; c/b worth retesting:
1142                                        //     / is symmetric
1143                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1144                                        // }
1145                                        ++j_free;
1146                                }
1147                        } // j loop
1148                        beta[i_free] = beta_sum;
1149                        ++i_free;
1150                }
1151        } // i loop
1152       
1153        return 0;
1154}
1155
1156
1157/* And this is the variant which handles an instrument response. */
1158/* We assume that the function values are sensible. */
1159
1160int GCI_marquardt_compute_fn_instr(float xincr, float y[], int ndata,
1161                                   int fit_start, int fit_end,
1162                                   float instr[], int ninstr,
1163                                   noise_type noise, float sig[],
1164                                   float param[], int paramfree[], int nparam,
1165                                   void (*fitfunc)(float, float [], float *, float [], int),
1166                                   float yfit[], float dy[],
1167                                   float **alpha, float beta[], float *chisq, float old_chisq,
1168                                   float alambda,
1169                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1170                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1171{
1172        int i, j, k, l, m, mfit, ret;
1173        float wt, sig2i, y_ymod;
1174        float alpha_weight[MAXBINS];
1175        float beta_weight[MAXBINS];
1176        int q;
1177        float weight;
1178        int i_free;
1179        int j_free;
1180        float dot_product;
1181        float beta_sum;
1182        float dy_dparam_k_i;
1183       
1184        /* Are we initialising? */
1185        // Malloc the arrays that will get used again in this fit via the pointers passed in
1186        // They will be freed by the higher fn that declared them.
1187        if (alambda < 0) {
1188                /* do any necessary initialisation */
1189                if ((*pfnvals_len) < ndata) {  /* we will need this length for
1190                                                                                the final full computation */
1191                        if ((*pfnvals_len)) {
1192                                free((*pfnvals));
1193                                GCI_ecf_free_matrix((*pdy_dparam_pure));
1194                                GCI_ecf_free_matrix((*pdy_dparam_conv));
1195                                (*pfnvals_len) = 0;
1196                                (*pdy_dparam_nparam_size) = 0;
1197                        }
1198                } else if ((*pdy_dparam_nparam_size) < nparam) {
1199                        GCI_ecf_free_matrix((*pdy_dparam_pure));
1200                        GCI_ecf_free_matrix((*pdy_dparam_conv));
1201                        (*pdy_dparam_nparam_size) = 0;
1202                }
1203                if (! (*pfnvals_len)) {
1204                        if (((*pfnvals) = (float *) malloc(ndata * sizeof(float)))
1205                                == NULL)
1206                                return -1;
1207                        (*pfnvals_len) = ndata;
1208                }
1209                if (! (*pdy_dparam_nparam_size)) {
1210                        /* use (*pfnvals_len), not ndata here! */
1211                        if (((*pdy_dparam_pure) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
1212                                /* Can't keep (*pfnvals) around in this case */
1213                                free((*pfnvals));
1214                                (*pfnvals_len) = 0;
1215                                return -1;
1216                        }
1217                        if (((*pdy_dparam_conv) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
1218                                /* Can't keep (*pfnvals) around in this case */
1219                                free((*pfnvals));
1220                                free((*pdy_dparam_pure));
1221                                (*pfnvals_len) = 0;
1222                                return -1;
1223                        }
1224                        (*pdy_dparam_nparam_size) = nparam;
1225                }
1226        }
1227
1228        for (j=0, mfit=0; j<nparam; j++)
1229                if (paramfree[j]) mfit++;
1230
1231        /* Calculation of the fitting data will depend upon the type of
1232           noise and the type of instrument response */
1233
1234        /* Need to calculate unconvolved values all the way down to 0 for
1235           the instrument response case */
1236        if (ninstr > 0) {
1237                if (fitfunc == GCI_multiexp_lambda)
1238                        ret = multiexp_lambda_array(xincr, param, (*pfnvals),
1239                                                                                (*pdy_dparam_pure), fit_end, nparam);
1240                else if (fitfunc == GCI_multiexp_tau)
1241                        ret = multiexp_tau_array(xincr, param, (*pfnvals),
1242                                                                         (*pdy_dparam_pure), fit_end, nparam);
1243                else if (fitfunc == GCI_stretchedexp)
1244                        ret = stretchedexp_array(xincr, param, (*pfnvals),
1245                                                                         (*pdy_dparam_pure), fit_end, nparam);
1246                else
1247                        ret = -1;
1248
1249                if (ret < 0)
1250                        for (i=0; i<fit_end; i++)
1251                                (*fitfunc)(xincr*i, param, &(*pfnvals)[i],
1252                                                   (*pdy_dparam_pure)[i], nparam);
1253
1254                /* OK, we've got to convolve the model fit with the given
1255                   instrument response.  What we'll do here, then, is to
1256                   generate the whole model fit first, then do the convolution
1257                   with the instrument response.  Only after doing that will
1258                   we attempt to calculate the goodness of fit matrices.  This
1259                   means that we will be looping through all of the points
1260                   twice, which is not worth it if there is no convolution
1261                   necessary. */
1262
1263                for (i=fit_start; i<fit_end; i++) {
1264                        int convpts;
1265
1266                        /* We wish to find yfit = (*pfnvals) * instr, so explicitly:
1267                             yfit[i] = sum_{j=0}^i (*pfnvals)[i-j].instr[j]
1268                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1269                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1270                           (*pfnvals)[i-j].instr[j]
1271                        */
1272
1273                        /* Zero our adders */
1274                        yfit[i] = 0;
1275                        for (k=1; k<nparam; k++)
1276                                (*pdy_dparam_conv)[i][k] = 0;
1277
1278                        convpts = (ninstr <= i) ? ninstr-1 : i;
1279                        for (j=0; j<=convpts; j++) {
1280                                yfit[i] += (*pfnvals)[i-j] * instr[j];
1281                                for (k=1; k<nparam; k++)
1282                                        (*pdy_dparam_conv)[i][k] += (*pdy_dparam_pure)[i-j][k] * instr[j];
1283                        }
1284                }
1285        } else {
1286                /* Can go straight into the final arrays in this case */
1287                if (fitfunc == GCI_multiexp_lambda)
1288                        ret = multiexp_lambda_array(xincr, param, yfit,
1289                                                                                (*pdy_dparam_conv), fit_end, nparam);
1290                else if (fitfunc == GCI_multiexp_tau)
1291                        ret = multiexp_tau_array(xincr, param, yfit,
1292                                                                         (*pdy_dparam_conv), fit_end, nparam);
1293                else if (fitfunc == GCI_stretchedexp)
1294                        ret = stretchedexp_array(xincr, param, yfit,
1295                                                                         (*pdy_dparam_conv), fit_end, nparam);
1296                else
1297                        ret = -1;
1298
1299                if (ret < 0)
1300                        for (i=0; i<fit_end; i++)
1301                                (*fitfunc)(xincr*i, param, &yfit[i],
1302                                                   (*pdy_dparam_conv)[i], nparam);
1303        }
1304
1305        /* OK, now we've got our (possibly convolved) data, we can do the
1306           rest almost exactly as above. */
1307       
1308        *chisq = 0.0f;
1309       
1310        switch (noise) {
1311                case NOISE_CONST:
1312                {
1313                        for (q = fit_start; q < fit_end; ++q) {
1314                                (*pdy_dparam_conv)[q][0] = 1.0f;
1315                                yfit[q] += param[0];
1316                                dy[q] = y[q] - yfit[q];
1317                                weight = 1.0f / sig[0];
1318                                alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]);
1319                                weight *= dy[q];
1320                                beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]);
1321                                weight *= dy[q];
1322                                *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]);
1323                        }
1324                        break;
1325                }
1326                case NOISE_GIVEN:
1327                {
1328                        for (q = fit_start; q < fit_end; ++q) {
1329                                (*pdy_dparam_conv)[q][0] = 1.0f;
1330                                yfit[q] += param[0];
1331                                dy[q] = y[q] - yfit[q];
1332                                weight = 1.0f / (sig[q] * sig[q]);
1333                                alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
1334                                weight *= dy[q];
1335                                beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
1336                                weight *= dy[q];
1337                                *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
1338                        }
1339                        break;
1340                }
1341                case NOISE_POISSON_DATA:
1342                {
1343                        for (q = fit_start; q < fit_end; ++q) {
1344                                (*pdy_dparam_conv)[q][0] = 1.0f;
1345                                yfit[q] += param[0];
1346                                dy[q] = y[q] - yfit[q];
1347                                weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
1348                                alpha_weight[q] = weight; // 1 / sig(q)
1349                                weight *= dy[q];
1350                                beta_weight[q] = weight; // dy[q] / sig(q)
1351                                weight *= dy[q];
1352                                *chisq += weight; // (dy[q] * dy[q]) / sig(q)
1353                        }
1354                        break;
1355                }
1356                case NOISE_POISSON_FIT:
1357                {
1358                        for (q = fit_start; q < fit_end; ++q) {
1359                                (*pdy_dparam_conv)[q][0] = 1.0f;
1360                                yfit[q] += param[0];
1361                                dy[q] = y[q] - yfit[q];
1362                                weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1363                                alpha_weight[q] = weight; // 1 / sig(q)
1364                                weight *= dy[q];
1365                                beta_weight[q] = weight; // dy(q) / sig(q)
1366                                weight *= dy[q];
1367                                *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1368                        }
1369                        break;
1370                }
1371                case NOISE_GAUSSIAN_FIT:
1372                {
1373                        for (q = fit_start; q < fit_end; ++q) {
1374                                (*pdy_dparam_conv)[q][0] = 1.0f;
1375                                yfit[q] += param[0];
1376                                dy[q] = y[q] - yfit[q];
1377                                weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1378                                alpha_weight[q] = weight; // 1 / sig(q)
1379                                weight *= dy[q];
1380                                beta_weight[q] = weight; // dy[q] / sig(q)
1381                                weight *= dy[q];
1382                                *chisq += weight; // dy[q] / sig(q)
1383                        }
1384                        break;
1385                }
1386                case NOISE_MLE:
1387                {
1388                        for (q = fit_start; q < fit_end; ++q) {
1389                                (*pdy_dparam_conv)[q][0] = 1.0f;
1390                                yfit[q] += param[0];
1391                                dy[q] = y[q] - yfit[q];
1392                                weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f);
1393                                alpha_weight[q] = weight * y[q] / yfit[q];
1394                                beta_weight[q] = dy[q] * weight;
1395                                if (yfit[q] > 0.0) {
1396                                        *chisq += (0.0f == y[q])
1397                                                        ? 2.0 * yfit[q]
1398                                                        : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1399                                }
1400                        }
1401                        if (*chisq <= 0.0f) {
1402                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1403                        }
1404                        break;
1405                }
1406                default:
1407                {
1408                        return -3;
1409                }
1410        }
1411
1412        // Check if chi square worsened:
1413        if (0.0f != old_chisq && *chisq >= old_chisq) {
1414                // don't bother to set up the matrices for solution
1415                return 0;
1416        }
1417
1418        i_free = 0;
1419        // for all columns
1420        for (i = 0; i < nparam; ++i) {
1421                if (paramfree[i]) {
1422                        j_free = 0;
1423                        beta_sum = 0.0f;
1424                        // row loop, only need to consider lower triangle
1425                        for (j = 0; j <= i; ++j) {
1426                                if (paramfree[j]) {
1427                                        dot_product = 0.0f;
1428                                        if (0 == j_free) { // true only once for each outer loop i
1429                                                // for all data
1430                                                for (k = fit_start; k < fit_end; ++k) {
1431                                                        dy_dparam_k_i = (*pdy_dparam_conv)[k][i];
1432                                                        dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO ARG make it [i][k] and just *dy_dparam++ it.
1433                                                        beta_sum += dy_dparam_k_i * beta_weight[k];
1434                                                }
1435                                        }
1436                                        else {
1437                                                // for all data
1438                                                for (k = fit_start; k < fit_end; ++k) {
1439                                                        dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k];
1440                                                }
1441                                        } // k loop
1442                                       
1443                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1444                                        // if (i_free != j_free) {
1445                                        //     // matrix is symmetric
1446                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1447                                        // }
1448                                        ++j_free;
1449                                }
1450                        } // j loop
1451                        beta[i_free] = beta_sum;
1452                        ++i_free;
1453                }
1454        } // i loop
1455
1456        return 0;
1457}
1458
1459/* These two variants, used just before the Marquardt fitting
1460   functions terminate, compute the function values at all points,
1461   whether or not they are being fitted.  (All points are fitted in
1462   the non-instrument response variant.)  They also compute the
1463   residuals y - yfit at all of those points and compute a chi-squared
1464   value which is not modified at small data values in the POISSON
1465   noise models.  They do not calculate alpha or beta. */
1466
1467int GCI_marquardt_compute_fn_final(float x[], float y[], int ndata,
1468                                         noise_type noise, float sig[],
1469                                         float param[], int paramfree[], int nparam,
1470                                         void (*fitfunc)(float, float [], float *, float [], int),
1471                                         float yfit[], float dy[], float *chisq)
1472{
1473        int i, j, mfit;
1474        float sig2i, dy_dparam[MAXFIT];  /* dy_dparam needed for fitfunc */
1475
1476        for (j=0, mfit=0; j<nparam; j++)
1477                if (paramfree[j])
1478                        mfit++;
1479
1480        /* Calculation of the fitting data will depend upon the type of
1481           noise and the type of instrument response */
1482
1483        /* There's no convolution involved in this function.  This is then
1484           fairly straightforward, depending only upon the type of noise
1485           present.  Since we calculate the standard deviation at every
1486           point in a different way depending upon the type of noise, we
1487           will place our switch(noise) around the entire loop. */
1488
1489        switch (noise) {
1490        case NOISE_CONST:
1491                *chisq = 0.0;
1492                /* Summation loop over all data */
1493                for (i=0; i<ndata; i++) {
1494                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1495                        yfit[i] += param[0];
1496                        dy[i] = y[i] - yfit[i];
1497                        /* And find chi^2 */
1498                        *chisq += dy[i] * dy[i];
1499                }
1500
1501                /* Now divide everything by sigma^2 */
1502                sig2i = 1.0 / (sig[0] * sig[0]);
1503                *chisq *= sig2i;
1504                break;
1505
1506        case NOISE_GIVEN:  /* This is essentially the NR version */
1507                *chisq = 0.0;
1508                /* Summation loop over all data */
1509                for (i=0; i<ndata; i++) {
1510                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1511                        yfit[i] += param[0];
1512                        sig2i = 1.0 / (sig[i] * sig[i]);
1513                        dy[i] = y[i] - yfit[i];
1514                        /* And find chi^2 */
1515                        *chisq += dy[i] * dy[i] * sig2i;
1516                }
1517                break;
1518
1519        case NOISE_POISSON_DATA:
1520                *chisq = 0.0;
1521                /* Summation loop over all data */
1522                for (i=0; i<ndata; i++) {
1523                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1524                        yfit[i] += param[0];
1525                        /* we still don't let the variance drop below 1 */
1526                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1527                        dy[i] = y[i] - yfit[i];
1528                        /* And find chi^2 */
1529                        *chisq += dy[i] * dy[i] * sig2i;
1530                }
1531                break;
1532
1533        case NOISE_POISSON_FIT:
1534                *chisq = 0.0;
1535                // Summation loop over all data
1536                for (i=0; i<ndata; i++) {
1537                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1538                        yfit[i] += param[0];
1539                        // we still don't let the variance drop below 1
1540                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1541                        dy[i] = y[i] - yfit[i];
1542                        // And find chi^2
1543                        *chisq += dy[i] * dy[i] * sig2i;
1544                }
1545                break;
1546
1547        case NOISE_MLE:
1548                        *chisq = 0.0;
1549                        /* Summation loop over all data */
1550                        for (i=0; i<ndata; i++) {
1551                                (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1552                                yfit[i] += param[0];
1553//                              dy[i] = y[i] - yfit[i];
1554
1555                                /* And find chi^2 */
1556//                              sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1557//                              *chisq += dy[i] * dy[i] * sig2i;
1558                                if (yfit[i]<=0.0)
1559                                        ; // do nothing
1560                                else if (y[i]==0.0)
1561                                        *chisq += 2.0*yfit[i];   // to avoid NaN from log
1562                                else
1563                                        *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1564                        }
1565                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1566                break;
1567
1568        case NOISE_GAUSSIAN_FIT:
1569                *chisq = 0.0;
1570                // Summation loop over all data
1571                for (i=0; i<ndata; i++) {
1572                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1573                        yfit[i] += param[0];
1574                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1575                        dy[i] = y[i] - yfit[i];
1576                        // And find chi^2
1577                        *chisq += dy[i] * dy[i] * sig2i;
1578                }
1579                break;
1580
1581        default:
1582                return -3;
1583                /* break; */   // (unreachable code)
1584        }
1585
1586        return 0;
1587}
1588
1589
1590/* And this is the variant which handles an instrument response. */
1591/* We assume that the function values are sensible.  Note also that we
1592   have to be careful about which points which include in the
1593   chi-squared calculation.  Also, we are guaranteed that the
1594   initialisation of the convolution arrays has been performed. */
1595
1596int GCI_marquardt_compute_fn_final_instr(float xincr, float y[], int ndata,
1597                                   int fit_start, int fit_end,
1598                                   float instr[], int ninstr,
1599                                   noise_type noise, float sig[],
1600                                   float param[], int paramfree[], int nparam,
1601                                   void (*fitfunc)(float, float [], float *, float [], int),
1602                                   float yfit[], float dy[], float *chisq,
1603                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1604                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1605{
1606        int i, j, mfit, ret;
1607        float sig2i;
1608        float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
1609        int fnvals_len = *pfnvals_len;
1610        int dy_dparam_nparam_size = *pdy_dparam_nparam_size;
1611
1612        /* check the necessary initialisation for safety, bail out if
1613           broken */
1614        if ((fnvals_len < ndata) || (dy_dparam_nparam_size < nparam))
1615                return -1;
1616        fnvals = *pfnvals;
1617        dy_dparam_pure = *pdy_dparam_pure;
1618        dy_dparam_conv = *pdy_dparam_conv;
1619
1620        for (j=0, mfit=0; j<nparam; j++)
1621                if (paramfree[j]) mfit++;
1622
1623        /* Calculation of the fitting data will depend upon the type of
1624           noise and the type of instrument response */
1625
1626        /* Need to calculate unconvolved values all the way down to 0 for
1627           the instrument response case */
1628        if (ninstr > 0) {
1629                if (fitfunc == GCI_multiexp_lambda)
1630                        ret = multiexp_lambda_array(xincr, param, fnvals,
1631                                                                                dy_dparam_pure, ndata, nparam);
1632                else if (fitfunc == GCI_multiexp_tau)
1633                        ret = multiexp_tau_array(xincr, param, fnvals,
1634                                                                         dy_dparam_pure, ndata, nparam);
1635                else if (fitfunc == GCI_stretchedexp)
1636                        ret = stretchedexp_array(xincr, param, fnvals,
1637                                                                         dy_dparam_pure, ndata, nparam);
1638                else
1639                        ret = -1;
1640
1641                if (ret < 0)
1642                        for (i=0; i<ndata; i++)
1643                                (*fitfunc)(xincr*i, param, &fnvals[i],
1644                                                   dy_dparam_pure[i], nparam);
1645
1646                /* OK, we've got to convolve the model fit with the given
1647                   instrument response.  What we'll do here, then, is to
1648                   generate the whole model fit first, then do the convolution
1649                   with the instrument response.  Only after doing that will
1650                   we attempt to calculate the goodness of fit matrices.  This
1651                   means that we will be looping through all of the points
1652                   twice, which is not worth it if there is no convolution
1653                   necessary. */
1654
1655                for (i=0; i<ndata; i++) {
1656                        int convpts;
1657
1658                        /* We wish to find yfit = fnvals * instr, so explicitly:
1659                             yfit[i] = sum_{j=0}^i fnvals[i-j].instr[j]
1660                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1661                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1662                           fnvals[i-j].instr[j]
1663                        */
1664
1665                        /* Zero our adder; don't need to bother with dy_dparam
1666                           stuff here */
1667                        yfit[i] = 0;
1668
1669                        convpts = (ninstr <= i) ? ninstr-1 : i;
1670                        for (j=0; j<=convpts; j++)
1671                                yfit[i] += fnvals[i-j] * instr[j];
1672                }
1673        } else {
1674                /* Can go straight into the final arrays in this case */
1675                if (fitfunc == GCI_multiexp_lambda)
1676                        ret = multiexp_lambda_array(xincr, param, yfit,
1677                                                                                dy_dparam_conv, ndata, nparam);
1678                else if (fitfunc == GCI_multiexp_tau)
1679                        ret = multiexp_tau_array(xincr, param, yfit,
1680                                                                         dy_dparam_conv, ndata, nparam);
1681                else if (fitfunc == GCI_stretchedexp)
1682                        ret = stretchedexp_array(xincr, param, yfit,
1683                                                                         dy_dparam_conv, ndata, nparam);
1684                else
1685                        ret = -1;
1686
1687                if (ret < 0)
1688                        for (i=0; i<ndata; i++)
1689                                (*fitfunc)(xincr*i, param, &yfit[i],
1690                                                   dy_dparam_conv[i], nparam);
1691        }
1692
1693        /* OK, now we've got our (possibly convolved) data, we can do the
1694           rest almost exactly as above. */
1695
1696        switch (noise) {
1697        case NOISE_CONST:
1698                *chisq = 0.0;
1699                /* Summation loop over all data */
1700                for (i=0; i<fit_start; i++) {
1701                        yfit[i] += param[0];
1702                        dy[i] = y[i] - yfit[i];
1703                }
1704                for ( ; i<fit_end; i++) {
1705                        yfit[i] += param[0];
1706                        dy[i] = y[i] - yfit[i];
1707                        /* And find chi^2 */
1708                        *chisq += dy[i] * dy[i];
1709                }
1710                for ( ; i<ndata; i++) {
1711                        yfit[i] += param[0];
1712                        dy[i] = y[i] - yfit[i];
1713                }
1714
1715                /* Now divide chi-squared by sigma^2 */
1716                sig2i = 1.0 / (sig[0] * sig[0]);
1717                *chisq *= sig2i;
1718                break;
1719
1720        case NOISE_GIVEN:  /* This is essentially the NR version */
1721                *chisq = 0.0;
1722                /* Summation loop over all data */
1723                for (i=0; i<fit_start; i++) {
1724                        yfit[i] += param[0];
1725                        dy[i] = y[i] - yfit[i];
1726                }
1727                for ( ; i<fit_end; i++) {
1728                        yfit[i] += param[0];
1729                        dy[i] = y[i] - yfit[i];
1730                        /* And find chi^2 */
1731                        sig2i = 1.0 / (sig[i] * sig[i]);
1732                        *chisq += dy[i] * dy[i] * sig2i;
1733                }
1734                for ( ; i<ndata; i++) {
1735                        yfit[i] += param[0];
1736                        dy[i] = y[i] - yfit[i];
1737                }
1738                break;
1739
1740        case NOISE_POISSON_DATA:
1741                *chisq = 0.0;
1742                /* Summation loop over all data */
1743                for (i=0; i<fit_start; i++) {
1744                        yfit[i] += param[0];
1745                        dy[i] = y[i] - yfit[i];
1746                }
1747                for ( ; i<fit_end; i++) {
1748                        yfit[i] += param[0];
1749                        dy[i] = y[i] - yfit[i];
1750                        /* And find chi^2 */
1751                        /* we still don't let the variance drop below 1 */
1752                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1753                        *chisq += dy[i] * dy[i] * sig2i;
1754                }
1755                for (; i<ndata; i++) {
1756                        yfit[i] += param[0];
1757                        dy[i] = y[i] - yfit[i];
1758                }
1759                break;
1760
1761        case NOISE_POISSON_FIT:
1762                *chisq = 0.0;
1763                // Summation loop over all data
1764                for (i=0; i<fit_start; i++) {
1765                        yfit[i] += param[0];
1766                        dy[i] = y[i] - yfit[i];
1767                }
1768                for ( ; i<fit_end; i++) {
1769                        yfit[i] += param[0];
1770                        dy[i] = y[i] - yfit[i];
1771                        // And find chi^2
1772                        // we still don't let the variance drop below 1
1773                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1774                        *chisq += dy[i] * dy[i] * sig2i;
1775                }
1776                for ( ; i<ndata; i++) {
1777                        yfit[i] += param[0];
1778                        dy[i] = y[i] - yfit[i];
1779                }
1780                break;
1781
1782        case NOISE_MLE:                              // for the final chisq report a normal chisq measure for MLE
1783                *chisq = 0.0;
1784                // Summation loop over all data
1785                for (i=0; i<fit_start; i++) {
1786                        yfit[i] += param[0];
1787                        dy[i] = y[i] - yfit[i];
1788                }
1789                for ( ; i<fit_end; i++) {
1790                        yfit[i] += param[0];
1791                        dy[i] = y[i] - yfit[i];
1792                        // And find chi^2
1793//                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1794                        if (yfit[i]<=0.0)
1795                                ; // do nothing
1796                        else if (y[i]==0.0)
1797                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1798                        else
1799                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1800                }
1801                for ( ; i<ndata; i++) {
1802                        yfit[i] += param[0];
1803                        dy[i] = y[i] - yfit[i];
1804                }
1805                if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1806                break;
1807
1808        case NOISE_GAUSSIAN_FIT:
1809                *chisq = 0.0;
1810                // Summation loop over all data
1811                for (i=0; i<fit_start; i++) {
1812                        yfit[i] += param[0];
1813                        dy[i] = y[i] - yfit[i];
1814                }
1815                for ( ; i<fit_end; i++) {
1816                        yfit[i] += param[0];
1817                        dy[i] = y[i] - yfit[i];
1818                        // And find chi^2
1819                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1820                        *chisq += dy[i] * dy[i] * sig2i;
1821                }
1822                for ( ; i<ndata; i++) {
1823                        yfit[i] += param[0];
1824                        dy[i] = y[i] - yfit[i];
1825                }
1826                break;
1827
1828        default:
1829                return -3;
1830                /* break; */   // (unreachable code)
1831        }
1832
1833        return 0;
1834}
1835
1836
1837//********************************* GCI_marquardt_fitting_engine **********************************************************************
1838
1839// This returns the number of iterations or negative if an error occurred.
1840// passes all the data to the ecf routine, checks the returned chisq and re-fits if it is of benefit
1841// was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03
1842
1843int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end,
1844                                                float prompt[], int nprompt,
1845                                                noise_type noise, float sig[],
1846                                                float param[], int paramfree[],
1847                                           int nparam, restrain_type restrain,
1848                                           void (*fitfunc)(float, float [], float *, float [], int),
1849                                           float *fitted, float *residuals, float *chisq,
1850                                           float **covar, float **alpha, float **erraxes,
1851                                           float chisq_target, float chisq_delta, int chisq_percent)
1852{
1853        float oldChisq, local_chisq;
1854        int ret, tries=0;
1855
1856        if (ecf_exportParams) ecf_ExportParams_OpenFile ();
1857
1858        // All of the work is done by the ECF module
1859        ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1860                                                          prompt, nprompt, noise, sig,
1861                                                          param, paramfree, nparam, restrain, fitfunc,
1862                                                          fitted, residuals, covar, alpha, &local_chisq,
1863                                                          chisq_delta, chisq_percent, erraxes);
1864
1865        // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same
1866        // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr()
1867        oldChisq = 3.0e38;
1868        while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<MAXREFITS)
1869        {
1870                oldChisq = local_chisq;
1871                tries++;
1872                ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1873                                                          prompt, nprompt, noise, sig,
1874                                                          param, paramfree, nparam, restrain, fitfunc,
1875                                                          fitted, residuals, covar, alpha, &local_chisq,
1876                                                          chisq_delta, chisq_percent, erraxes);
1877        }
1878
1879        if (chisq!=NULL) *chisq = local_chisq;
1880
1881        if (ecf_exportParams) ecf_ExportParams_CloseFile ();
1882
1883        return ret;             // summed number of iterations
1884}
1885
1886/* Cleanup function */
1887void GCI_marquardt_cleanup(void)
1888{
1889/*      if (fnvals_len) {
1890                free(fnvals);
1891                GCI_ecf_free_matrix(dy_dparam_pure);
1892                GCI_ecf_free_matrix(dy_dparam_conv);
1893                fnvals_len = 0;
1894                dy_dparam_nparam_size = 0;
1895        }
1896*/
1897}
1898
1899
1900// Emacs settings:
1901// Local variables:
1902// mode: c
1903// c-basic-offset: 4
1904// tab-width: 4
1905// End:
Note: See TracBrowser for help on using the repository browser.