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

Revision 7708, 66.4 KB checked in by aivar, 9 years ago (diff)

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

Tidied up a bit.

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

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
702//TODO ARG deleted in my latest version
703int GCI_marquardt_step(float x[], float y[], int ndata,
704                                        noise_type noise, float sig[],
705                                        float param[], int paramfree[], int nparam,
706                                        restrain_type restrain,
707                                        void (*fitfunc)(float, float [], float *, float [], int),
708                                        float yfit[], float dy[],
709                                        float **covar, float **alpha, float *chisq,
710                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam)
711{
712        int j, k, l, ret;
713//      static int mfit;   // was static but now thread safe
714//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe
715        int mfit = *pmfit;
716        float ochisq = *pochisq;
717
718        if (nparam > MAXFIT)
719                return -1;
720
721        /* Initialisation */
722        /* We assume we're given sensible starting values for param[] */
723        if (*alambda < 0.0) {
724                for (mfit=0, j=0; j<nparam; j++)
725                        if (paramfree[j])
726                                mfit++;
727
728                if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig,
729                                                                         param, paramfree, nparam, fitfunc,
730                                                                         yfit, dy,
731                                                                         alpha, beta, chisq, 0.0, *alambda) != 0)
732                        return -2;
733
734                *alambda = 0.001;
735                ochisq = (*chisq);
736                for (j=0; j<nparam; j++)
737                        paramtry[j] = param[j];
738        }
739
740        /* Alter linearised fitting matrix by augmenting diagonal elements */
741        for (j=0; j<mfit; j++) {
742                for (k=0; k<mfit; k++)
743                        covar[j][k] = alpha[j][k];
744
745                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
746                dparam[j] = beta[j];
747        }
748
749        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */
750        if (GCI_solve(covar, mfit, dparam) != 0)
751                return -1;
752
753        //TODO ARG GCI_gauss_jordan would invert the covar matrix as a side effect
754        /* Once converged, evaluate covariance matrix */
755        if (*alambda == 0) {
756                if (GCI_marquardt_compute_fn_final(x, y, ndata, noise, sig,
757                                                                                   param, paramfree, nparam, fitfunc,
758                                                                                   yfit, dy, chisq) != 0)
759                        return -3;
760                if (mfit < nparam) {  /* no need to do this otherwise */
761                        GCI_covar_sort(covar, nparam, paramfree, mfit);
762                        GCI_covar_sort(alpha, nparam, paramfree, mfit);
763                }
764                return 0;
765        }
766
767        /* Did the trial succeed? */
768        for (j=0, l=0; l<nparam; l++)
769                if (paramfree[l])
770                        paramtry[l] = param[l] + dparam[j++];
771
772        if (restrain == ECF_RESTRAIN_DEFAULT)
773                ret = check_ecf_params (paramtry, nparam, fitfunc);
774        else
775                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
776
777        if (ret != 0) {
778                /* Bad parameters, increase alambda and return */
779                *alambda *= 10.0;
780                return 0;
781        }
782
783        if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig,
784                                                                 paramtry, paramfree, nparam, fitfunc,
785                                                                 yfit, dy, covar, dparam,
786                                                                 chisq, ochisq, *alambda) != 0)
787                return -2;
788
789        /* Success, accept the new solution */
790        if (*chisq < ochisq) {
791                *alambda *= 0.1;
792                ochisq = *chisq;
793                for (j=0; j<mfit; j++) {
794                        for (k=0; k<mfit; k++)
795                                alpha[j][k] = covar[j][k];
796                        beta[j] = dparam[j];
797                }
798                for (l=0; l<nparam; l++)
799                        param[l] = paramtry[l];
800        } else { /* Failure, increase alambda and return */
801                *alambda *= 10.0;
802                *chisq = ochisq;
803        }
804
805        return 0;
806}
807
808
809int GCI_marquardt_step_instr(float xincr, float y[],
810                                        int ndata, int fit_start, int fit_end,
811                                        float instr[], int ninstr,
812                                        noise_type noise, float sig[],
813                                        float param[], int paramfree[], int nparam,
814                                        restrain_type restrain,
815                                        void (*fitfunc)(float, float [], float *, float [], int),
816                                        float yfit[], float dy[],
817                                        float **covar, float **alpha, float *chisq,
818                                        float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam,
819                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
820                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
821{
822        int j, k, l, ret;
823//      static int mfit;   // was static but now thread safe
824//      static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];   // was static but now thread safe
825
826        if (nparam > MAXFIT)
827                return -10;
828        if (xincr <= 0)
829                return -11;
830        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
831                return -12;
832
833        /* Initialisation */
834        /* We assume we're given sensible starting values for param[] */
835        if (*alambda < 0.0) {
836
837                for ((*pmfit)=0, j=0; j<nparam; j++)
838                        if (paramfree[j])
839                                (*pmfit)++;
840
841                if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
842                                                                                   instr, ninstr, noise, sig,
843                                                                                   param, paramfree, nparam, fitfunc,
844                                                                                   yfit, dy, alpha, beta, chisq, 0.0,
845                                                                                   *alambda,
846                                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
847                                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0)
848                        return -2;
849
850                *alambda = 0.001;
851                *pochisq = *chisq;
852                for (j=0; j<nparam; j++)
853                        paramtry[j] = param[j];
854
855        }
856
857        /* Alter linearised fitting matrix by augmenting diagonal elements */
858        for (j=0; j<(*pmfit); j++) {
859                for (k=0; k<(*pmfit); k++)
860                        covar[j][k] = alpha[j][k];
861
862                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
863                dparam[j] = beta[j];
864        }
865
866        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */
867        if (GCI_solve(covar, *pmfit, dparam) != 0)
868                return -1;
869
870//TODO ARG covar needs to get inverted; previously inverted as a side effect
871        /* Once converged, evaluate covariance matrix */
872        if (*alambda == 0) {
873                if (GCI_marquardt_compute_fn_final_instr(
874                                                                        xincr, y, ndata, fit_start, fit_end,
875                                                                        instr, ninstr, noise, sig,
876                                                                        param, paramfree, nparam, fitfunc,
877                                                                        yfit, dy, chisq,
878                                                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
879                                                                        pfnvals_len, pdy_dparam_nparam_size) != 0)
880                        return -3;
881                if (*pmfit < nparam) {  /* no need to do this otherwise */
882                        GCI_covar_sort(covar, nparam, paramfree, *pmfit);
883                        GCI_covar_sort(alpha, nparam, paramfree, *pmfit);
884                }
885                return 0;
886        }
887
888//TODO ARG c/b special case if dparam all zeroes; don't have to recalc alpha & beta
889        /* Did the trial succeed? */
890        for (j=0, l=0; l<nparam; l++)
891                if (paramfree[l])
892                        paramtry[l] = param[l] + dparam[j++];
893
894        if (restrain == ECF_RESTRAIN_DEFAULT)
895                ret = check_ecf_params (paramtry, nparam, fitfunc);
896        else
897                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
898
899        if (ret != 0) {
900                /* Bad parameters, increase alambda and return */
901                *alambda *= 10.0;
902                return 0;
903        }
904
905        if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
906                                                                           instr, ninstr, noise, sig,
907                                                                           paramtry, paramfree, nparam, fitfunc,
908                                                                           yfit, dy, covar, dparam,
909                                                                           chisq, *pochisq, *alambda,
910                                                                           pfnvals, pdy_dparam_pure, pdy_dparam_conv,
911                                                                           pfnvals_len, pdy_dparam_nparam_size) != 0)
912                return -2;
913
914        /* Success, accept the new solution */
915        if (*chisq < *pochisq) {
916                *alambda *= 0.1;
917                *pochisq = *chisq;
918                for (j=0; j<(*pmfit); j++) {
919                        for (k=0; k<(*pmfit); k++)
920                                alpha[j][k] = covar[j][k];
921                        beta[j] = dparam[j];
922                }
923                for (l=0; l<nparam; l++)
924                        param[l] = paramtry[l];
925        } else { /* Failure, increase alambda and return */
926                *alambda *= 10.0;
927                *chisq = *pochisq;
928        }
929
930        return 0;
931}
932
933
934/* Used by GCI_marquardt to evaluate the linearised fitting matrix alpha
935   and vector beta and to calculate chi^2.      The equations involved are
936   given in section 15.5 of Numerical Recipes; basically:
937
938     \chi^2(param) = \sum_{i=1}^N ( (y_i-y(x_i;param)) / sigma_i )^2
939
940     beta_k = -1/2 (d/dparam_k)(chi^2)
941
942     alpha_kl = \sum_{i=1}^N (1/sigma_i^2) .
943                  (dy(x_i;param)/dparam_k) . (dy(x_i;param)/dparam_l)
944
945   (where all of the derivatives are partial).
946
947   If an instrument response is provided, we also take account of it
948   now.  We are given that:
949
950     observed(t) = response(t) * instr(t)
951
952   where response(t) is being fitted with fitfunc; it is also trivial to
953   show that (assuming that instr(t) is known and fixed, with no
954   dependencies on the param_k, the parameters being fitted):
955
956     (d/dparam_k) observed(t) = ((d/dparam_k) response(t)) * instr(t)
957
958   so we do not need to alter the response function in any way to
959   determined the fitted convolved response.
960
961   Again there are two variants of this function, corresponding to the
962   two variants of the Marquardt function.
963*/
964
965//TODO ARG deleted in my version; needs to get de-NRed!!!!
966int GCI_marquardt_compute_fn(float x[], float y[], int ndata,
967                                         noise_type noise, float sig[],
968                                         float param[], int paramfree[], int nparam,
969                                         void (*fitfunc)(float, float [], float *, float [], int),
970                                         float yfit[], float dy[],
971                                         float **alpha, float beta[], float *chisq, float old_chisq,
972                                         float alambda)
973{
974        int i, j, k, l, m, mfit;
975        float wt, sig2i, y_ymod, dy_dparam[MAXFIT];
976
977        for (j=0, mfit=0; j<nparam; j++)
978                if (paramfree[j])
979                        mfit++;
980
981        //TODO ARG having this section { } of code here is just temporary;
982        // o'wise have to define these vars at top of function
983        {
984        float alpha_weight[MAXBINS];
985        float beta_weight[MAXBINS];
986        int q;
987        float weight;
988
989        int i_free;
990        int j_free;
991        float dot_product;
992        float beta_sum;
993        float dy_dparam_k_i;
994
995        *chisq = 0.0f;
996
997        switch (noise) {
998            case NOISE_CONST:
999            {
1000                for (q = 0; q < ndata; ++q) {
1001                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam);
1002                    yfit[q] += param[0];
1003                    dy[q] = y[q] - yfit[q];
1004                    weight = 1.0f; //TODO ARG this should be 1.0f / sig[0] !
1005                    alpha_weight[q] = weight; // 1
1006                    weight *= dy[q];
1007                    beta_weight[q] = weight; // dy[q]
1008                    weight *= dy[q];
1009                    *chisq += weight; // dy[q] * dy[q]
1010                }
1011                break;
1012            }
1013            case NOISE_GIVEN:
1014            {
1015                for (q = 0; q < ndata; ++q) {
1016                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam);
1017                    yfit[q] += param[0];
1018                    dy[q] = y[q] - yfit[q];
1019                    weight = 1.0f / (sig[q] * sig[q]);
1020                    alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
1021                    weight *= dy[q];
1022                    beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
1023                    weight *= dy[q];
1024                    *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
1025                }
1026                break;
1027            }
1028            case NOISE_POISSON_DATA:
1029            {
1030                for (q = 0; q < ndata; ++q) {
1031                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam);
1032                    yfit[q] += param[0];
1033                    dy[q] = y[q] - yfit[q];
1034                    weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
1035                    alpha_weight[q] = weight; // 1 / sig(q)
1036                    weight *= dy[q];
1037                    beta_weight[q] = weight; // dy[q] / sig(q)
1038                    weight *= dy[q];
1039                    *chisq += weight; // (dy[q] * dy[q]) / sig(q)
1040                }
1041                break;
1042            }
1043            case NOISE_POISSON_FIT:
1044            {
1045                for (q = 0; q < ndata; ++q) {
1046                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam);
1047                    yfit[q] += param[0];
1048                    dy[q] = y[q] - yfit[q];
1049                    weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1050                    alpha_weight[q] = weight; // 1 / sig(q)
1051                    weight *= dy[q];
1052                    beta_weight[q] = weight; // dy(q) / sig(q)
1053                    weight *= dy[q];
1054                    *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1055                }
1056                break;
1057            }
1058            case NOISE_GAUSSIAN_FIT:
1059            {
1060                 for (q = 0; q < ndata; ++q) {
1061                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam);
1062                    yfit[q] += param[0];
1063                    dy[q] = y[q] - yfit[q];
1064                    weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1065                    alpha_weight[q] = weight; // 1 / sig(q)
1066                    weight *= dy[q];
1067                    beta_weight[q] = weight; // dy[q] / sig(q)
1068                    weight *= dy[q];
1069                    *chisq += weight;
1070                 }
1071                 break;
1072           }
1073            case NOISE_MLE:
1074            {
1075                for (q = 0; q < ndata; ++q) {
1076                    (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam);
1077                    yfit[q] += param[0];
1078                    dy[q] = y[q] - yfit[q];
1079                    weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f);
1080                    alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus?
1081                    beta_weight[q] = dy[q] * weight;
1082                    *chisq += (0.0f == y[q])
1083                            ? 2.0 * yfit[q]
1084                            : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1085                }
1086                if (*chisq <= 0.0f) {
1087                    *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1088                }
1089                break;
1090            }
1091            default:
1092            {
1093                return -3;
1094            }
1095        }
1096
1097        // Check if chi square worsened:
1098        if (0.0f != old_chisq && *chisq >= old_chisq) {
1099            // don't bother to set up the matrices for solution
1100            return 0;
1101        }
1102
1103        i_free = 0;
1104        // for all columns
1105        for (i = 0; i < nparam; ++i) {
1106            if (paramfree[i]) {
1107                j_free = 0;
1108                beta_sum = 0.0f;
1109                // row loop, only need to consider lower triangle
1110                for (j = 0; j <= i; ++j) {
1111                    if (paramfree[j]) {
1112                        dot_product = 0.0f;
1113                        if (0 == j_free) {
1114                            // for all data
1115                            for (k = 0; k < ndata; ++k) {
1116                                //TODO ARG just to get this to compile, for now:
1117                  //TODO ARG from the _instr version:              dy_dparam_k_i = (*pdy_dparam_conv)[k][i];
1118                  //TODO ARG from the instr version              dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it.
1119                  //TODO ARG from the instr version              beta_sum += dy_dparam_k_i * beta_weight[k];
1120                            }
1121                        }
1122                        else {
1123                            // for all data
1124                            for (k = 0; k < ndata; ++k) {
1125                              //TODO ARG from the instr version:  dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k];
1126                            }
1127                        } // k loop
1128
1129                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1130                       // if (i_free != j_free) {
1131                       //     // matrix is symmetric
1132                       //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1133                       // }
1134                        ++j_free;
1135                    }
1136                } // j loop
1137                beta[i_free] = beta_sum;
1138                ++i_free;
1139            }
1140            //else printf("param not free %d\n", i);
1141        } // i loop
1142        }
1143
1144        return 0;
1145
1146
1147        /* Initialise (symmetric) alpha, beta */
1148        //TODO ARG FRI
1149        //for (j=0; j<mfit; j++) {
1150        //      for (k=0; k<=j; k++)
1151        //              alpha[j][k] = 0.0;
1152        //      beta[j] = 0.0;
1153        //}
1154
1155        /* Calculation of the fitting data will depend upon the type of
1156           noise and the type of instrument response */
1157
1158        /* There's no convolution involved in this function.  This is then
1159           fairly straightforward, depending only upon the type of noise
1160           present.  Since we calculate the standard deviation at every
1161           point in a different way depending upon the type of noise, we
1162           will place our switch(noise) around the entire loop. */
1163
1164        switch (noise) {
1165        case NOISE_CONST:
1166                *chisq = 0.0;
1167                /* Summation loop over all data */
1168                for (i=0; i<ndata; i++) {
1169                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1170                        yfit[i] += param[0];
1171                        dy_dparam[0] = 1;
1172                        dy[i] = y[i] - yfit[i];
1173                        for (j=0, l=0; l<nparam; l++) {
1174                                if (paramfree[l]) {
1175                                        wt = dy_dparam[l];       /* taken away the *sig2i from here */
1176                                        for (k=0, m=0; m<=l; m++)
1177                                                if (paramfree[m])
1178                                                        alpha[j][k++] += wt * dy_dparam[m];
1179                                        beta[j] += dy[i] * wt;
1180                                        j++;
1181                                }
1182                        }
1183                        /* And find chi^2 */
1184                        *chisq += dy[i] * dy[i];
1185                }
1186
1187                /* Now divide everything by sigma^2 */
1188                sig2i = 1.0 / (sig[0] * sig[0]);
1189                *chisq *= sig2i;
1190                for (j=0; j<mfit; j++) {
1191                        for (k=0; k<=j; k++)
1192                                alpha[j][k] *= sig2i;
1193                        beta[j] *= sig2i;
1194                }
1195                break;
1196
1197        case NOISE_GIVEN:  /* This is essentially the NR version */
1198                *chisq = 0.0;
1199                /* Summation loop over all data */
1200                for (i=0; i<ndata; i++) {
1201                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1202                        yfit[i] += param[0];
1203                        dy_dparam[0] = 1;
1204                        sig2i = 1.0 / (sig[i] * sig[i]);
1205                        dy[i] = y[i] - yfit[i];
1206                        for (j=0, l=0; l<nparam; l++) {
1207                                if (paramfree[l]) {
1208                                        wt = dy_dparam[l] * sig2i;
1209                                        for (k=0, m=0; m<=l; m++)
1210                                                if (paramfree[m])
1211                                                        alpha[j][k++] += wt * dy_dparam[m];
1212                                        beta[j] += wt * dy[i];
1213                                        j++;
1214                                }
1215                        }
1216                        /* And find chi^2 */
1217                        *chisq += dy[i] * dy[i] * sig2i;
1218                }
1219                break;
1220
1221        case NOISE_POISSON_DATA:
1222                *chisq = 0.0;
1223                /* Summation loop over all data */
1224                for (i=0; i<ndata; i++) {
1225                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1226                        yfit[i] += param[0];
1227                        dy_dparam[0] = 1;
1228                        sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15);
1229                        dy[i] = y[i] - yfit[i];
1230                        for (j=0, l=0; l<nparam; l++) {
1231                                if (paramfree[l]) {
1232                                        wt = dy_dparam[l] * sig2i;
1233                                        for (k=0, m=0; m<=l; m++)
1234                                                if (paramfree[m])
1235                                                        alpha[j][k++] += wt * dy_dparam[m];
1236                                        beta[j] += wt * dy[i];
1237                                        j++;
1238                                }
1239                        }
1240                        /* And find chi^2 */
1241                        *chisq += dy[i] * dy[i] * sig2i;
1242                }
1243                break;
1244
1245        case NOISE_POISSON_FIT:
1246                *chisq = 0.0;
1247                // Summation loop over all data
1248                for (i=0; i<ndata; i++) {
1249                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1250                        yfit[i] += param[0];
1251                        dy_dparam[0] = 1;
1252                        sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15);
1253                        dy[i] = y[i] - yfit[i];
1254                        for (j=0, l=0; l<nparam; l++) {
1255                                if (paramfree[l]) {
1256                                        wt = dy_dparam[l] * sig2i;
1257                                        for (k=0, m=0; m<=l; m++)
1258                                                if (paramfree[m])
1259                                                        alpha[j][k++] += wt * dy_dparam[m];
1260                                        beta[j] += wt * dy[i];
1261                                        j++;
1262                                }
1263                        }
1264                        // And find chi^2
1265                        *chisq += dy[i] * dy[i] * sig2i;
1266                }
1267                break;
1268
1269        case NOISE_GAUSSIAN_FIT:
1270                *chisq = 0.0;
1271                // Summation loop over all data
1272                for (i=0; i<ndata; i++) {
1273                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1274                        yfit[i] += param[0];
1275                        dy_dparam[0] = 1;
1276                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1277                        dy[i] = y[i] - yfit[i];
1278                        for (j=0, l=0; l<nparam; l++) {
1279                                if (paramfree[l]) {
1280                                        wt = dy_dparam[l] * sig2i;
1281                                        for (k=0, m=0; m<=l; m++)
1282                                                if (paramfree[m])
1283                                                        alpha[j][k++] += wt * dy_dparam[m];
1284                                        beta[j] += wt * dy[i];
1285                                        j++;
1286                                }
1287                        }
1288                        // And find chi^2
1289                        *chisq += dy[i] * dy[i] * sig2i;
1290                }
1291                break;
1292
1293        case NOISE_MLE:
1294                *chisq = 0.0;
1295                /* Summation loop over all data */
1296                for (i=0; i<ndata; i++) {
1297                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1298                        yfit[i] += param[0];
1299                        dy_dparam[0] = 1;
1300                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1301                        dy[i] = y[i] - yfit[i];
1302                        y_ymod=y[i]/yfit[i];
1303                        for (j=0, l=0; l<nparam; l++) {
1304                                if (paramfree[l]) {
1305                                        wt = dy_dparam[l] * sig2i;
1306                                        for (k=0, m=0; m<=l; m++)
1307                                                if (paramfree[m])
1308                                                        alpha[j][k++] += wt*dy_dparam[m]*y_ymod; //wt * dy_dparam[m];
1309                                        beta[j] += wt * dy[i];
1310                                        j++;
1311                                }
1312                        }
1313                        // And find chi^2
1314                        if (yfit[i]<=0.0)
1315                                ; // do nothing
1316                        else if (y[i]==0.0)
1317                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1318                        else
1319                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1320                }
1321                if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1322                break;
1323
1324        default:
1325                return -3;
1326                /* break; */   // (unreachable code)
1327        }
1328
1329        /* Fill in the symmetric side */
1330        for (j=1; j<mfit; j++)
1331                for (k=0; k<j; k++)
1332                        alpha[k][j] = alpha[j][k];
1333
1334        return 0;
1335}
1336
1337
1338/* And this is the variant which handles an instrument response. */
1339/* We assume that the function values are sensible. */
1340
1341int GCI_marquardt_compute_fn_instr(float xincr, float y[], int ndata,
1342                                   int fit_start, int fit_end,
1343                                   float instr[], int ninstr,
1344                                   noise_type noise, float sig[],
1345                                   float param[], int paramfree[], int nparam,
1346                                   void (*fitfunc)(float, float [], float *, float [], int),
1347                                   float yfit[], float dy[],
1348                                   float **alpha, float beta[], float *chisq, float old_chisq,
1349                                   float alambda,
1350                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1351                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1352{
1353        int i, j, k, l, m, mfit, ret;
1354        float wt, sig2i, y_ymod;
1355
1356        /* Are we initialising? */
1357        // Malloc the arrays that will get used again in this fit via the pointers passed in
1358        // They will be freed by the higher fn that declared them.
1359        if (alambda < 0) {
1360                /* do any necessary initialisation */
1361                if ((*pfnvals_len) < ndata) {  /* we will need this length for
1362                                                                          the final full computation */
1363                        if ((*pfnvals_len)) {
1364                                free((*pfnvals));
1365                                GCI_ecf_free_matrix((*pdy_dparam_pure));
1366                                GCI_ecf_free_matrix((*pdy_dparam_conv));
1367                                (*pfnvals_len) = 0;
1368                                (*pdy_dparam_nparam_size) = 0;
1369                        }
1370                } else if ((*pdy_dparam_nparam_size) < nparam) {
1371                        GCI_ecf_free_matrix((*pdy_dparam_pure));
1372                        GCI_ecf_free_matrix((*pdy_dparam_conv));
1373                        (*pdy_dparam_nparam_size) = 0;
1374                }
1375                if (! (*pfnvals_len)) {
1376                        if (((*pfnvals) = (float *) malloc(ndata * sizeof(float)))
1377                                == NULL)
1378                                return -1;
1379                        (*pfnvals_len) = ndata;
1380                }
1381                if (! (*pdy_dparam_nparam_size)) {
1382                        /* use (*pfnvals_len), not ndata here! */
1383                        if (((*pdy_dparam_pure) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
1384                                /* Can't keep (*pfnvals) around in this case */
1385                                free((*pfnvals));
1386                                (*pfnvals_len) = 0;
1387                                return -1;
1388                        }
1389                        if (((*pdy_dparam_conv) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
1390                                /* Can't keep (*pfnvals) around in this case */
1391                                free((*pfnvals));
1392                                free((*pdy_dparam_pure));
1393                                (*pfnvals_len) = 0;
1394                                return -1;
1395                        }
1396                        (*pdy_dparam_nparam_size) = nparam;
1397                }
1398        }
1399
1400        for (j=0, mfit=0; j<nparam; j++)
1401                if (paramfree[j]) mfit++;
1402
1403//TODO ARG first change:
1404//      /* Initialise (symmetric) alpha, beta */
1405//      for (j=0; j<mfit; j++) {
1406//              for (k=0; k<=j; k++)
1407//                      alpha[j][k] = 0.0;
1408//              beta[j] = 0.0;
1409//      }
1410
1411        /* Calculation of the fitting data will depend upon the type of
1412           noise and the type of instrument response */
1413
1414        /* Need to calculate unconvolved values all the way down to 0 for
1415           the instrument response case */
1416        if (ninstr > 0) {
1417                if (fitfunc == GCI_multiexp_lambda)
1418                        ret = multiexp_lambda_array(xincr, param, (*pfnvals),
1419                                                                                (*pdy_dparam_pure), fit_end, nparam);
1420                else if (fitfunc == GCI_multiexp_tau)
1421                        ret = multiexp_tau_array(xincr, param, (*pfnvals),
1422                                                                         (*pdy_dparam_pure), fit_end, nparam);
1423                else if (fitfunc == GCI_stretchedexp)
1424                        ret = stretchedexp_array(xincr, param, (*pfnvals),
1425                                                                         (*pdy_dparam_pure), fit_end, nparam);
1426                else
1427                        ret = -1;
1428
1429                if (ret < 0)
1430                        for (i=0; i<fit_end; i++)
1431                                (*fitfunc)(xincr*i, param, &(*pfnvals)[i],
1432                                                   (*pdy_dparam_pure)[i], nparam);
1433
1434                /* OK, we've got to convolve the model fit with the given
1435                   instrument response.  What we'll do here, then, is to
1436                   generate the whole model fit first, then do the convolution
1437                   with the instrument response.  Only after doing that will
1438                   we attempt to calculate the goodness of fit matrices.  This
1439                   means that we will be looping through all of the points
1440                   twice, which is not worth it if there is no convolution
1441                   necessary. */
1442
1443                for (i=fit_start; i<fit_end; i++) {
1444                        int convpts;
1445
1446                        /* We wish to find yfit = (*pfnvals) * instr, so explicitly:
1447                             yfit[i] = sum_{j=0}^i (*pfnvals)[i-j].instr[j]
1448                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1449                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1450                           (*pfnvals)[i-j].instr[j]
1451                        */
1452
1453                        /* Zero our adders */
1454                        yfit[i] = 0;
1455                        for (k=1; k<nparam; k++)
1456                                (*pdy_dparam_conv)[i][k] = 0;
1457
1458                        convpts = (ninstr <= i) ? ninstr-1 : i;
1459                        for (j=0; j<=convpts; j++) {
1460                                yfit[i] += (*pfnvals)[i-j] * instr[j];
1461                                for (k=1; k<nparam; k++)
1462                                        (*pdy_dparam_conv)[i][k] += (*pdy_dparam_pure)[i-j][k] * instr[j];
1463                        }
1464                }
1465        } else {
1466                /* Can go straight into the final arrays in this case */
1467                if (fitfunc == GCI_multiexp_lambda)
1468                        ret = multiexp_lambda_array(xincr, param, yfit,
1469                                                                                (*pdy_dparam_conv), fit_end, nparam);
1470                else if (fitfunc == GCI_multiexp_tau)
1471                        ret = multiexp_tau_array(xincr, param, yfit,
1472                                                                         (*pdy_dparam_conv), fit_end, nparam);
1473                else if (fitfunc == GCI_stretchedexp)
1474                        ret = stretchedexp_array(xincr, param, yfit,
1475                                                                         (*pdy_dparam_conv), fit_end, nparam);
1476                else
1477                        ret = -1;
1478
1479                if (ret < 0)
1480                        for (i=0; i<fit_end; i++)
1481                                (*fitfunc)(xincr*i, param, &yfit[i],
1482                                                   (*pdy_dparam_conv)[i], nparam);
1483        }
1484
1485        /* OK, now we've got our (possibly convolved) data, we can do the
1486           rest almost exactly as above. */
1487        //TODO ARG this new section of code is just temporary; o'wise have to define all these variables at the top of the function
1488        {
1489        float alpha_weight[MAXBINS];
1490        float beta_weight[MAXBINS];
1491        int q;
1492        float weight;
1493
1494        int i_free;
1495        int j_free;
1496        float dot_product;
1497        float beta_sum;
1498        float dy_dparam_k_i;
1499
1500        *chisq = 0.0f;
1501
1502                switch (noise) {
1503            case NOISE_CONST:
1504            {
1505                for (q = fit_start; q < fit_end; ++q) {
1506                    (*pdy_dparam_conv)[q][0] = 1.0f;
1507                    yfit[q] += param[0];
1508                    dy[q] = y[q] - yfit[q];
1509                    weight = 1.0f; //TODO ARG this should be 1.0f / sig[0] !
1510                    alpha_weight[q] = weight; // 1
1511                    weight *= dy[q];
1512                    beta_weight[q] = weight; // dy[q]
1513                    weight *= dy[q];
1514                    *chisq += weight; // dy[q] * dy[q]
1515                }
1516                break;
1517            }
1518            case NOISE_GIVEN:
1519            {
1520                for (q = fit_start; q < fit_end; ++q) {
1521                    (*pdy_dparam_conv)[q][0] = 1.0f;
1522                    yfit[q] += param[0];
1523                    dy[q] = y[q] - yfit[q];
1524                    weight = 1.0f / (sig[q] * sig[q]);
1525                    alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
1526                    weight *= dy[q];
1527                    beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
1528                    weight *= dy[q];
1529                    *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
1530                }
1531                break;
1532            }
1533            case NOISE_POISSON_DATA:
1534            {
1535                for (q = fit_start; q < fit_end; ++q) {
1536                    (*pdy_dparam_conv)[q][0] = 1.0f;
1537                    yfit[q] += param[0];
1538                    dy[q] = y[q] - yfit[q];
1539                    weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
1540                    alpha_weight[q] = weight; // 1 / sig(q)
1541                    weight *= dy[q];
1542                    beta_weight[q] = weight; // dy[q] / sig(q)
1543                    weight *= dy[q];
1544                    *chisq += weight; // (dy[q] * dy[q]) / sig(q)
1545                }
1546                break;
1547            }
1548            case NOISE_POISSON_FIT:
1549            {
1550                for (q = fit_start; q < fit_end; ++q) {
1551                    (*pdy_dparam_conv)[q][0] = 1.0f;
1552                    yfit[q] += param[0];
1553                    dy[q] = y[q] - yfit[q];
1554                    weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1555                    alpha_weight[q] = weight; // 1 / sig(q)
1556                    weight *= dy[q];
1557                    beta_weight[q] = weight; // dy(q) / sig(q)
1558                    weight *= dy[q];
1559                    *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1560                }
1561                break;
1562            }
1563            case NOISE_GAUSSIAN_FIT:
1564            {
1565                 for (q = fit_start; q < fit_end; ++q) {
1566                    (*pdy_dparam_conv)[q][0] = 1.0f;
1567                    yfit[q] += param[0];
1568                    dy[q] = y[q] - yfit[q];
1569                    weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1570                    alpha_weight[q] = weight; // 1 / sig(q)
1571                    weight *= dy[q];
1572                    beta_weight[q] = weight; // dy[q] / sig(q)
1573                    weight *= dy[q];
1574                    *chisq += weight;
1575                 }
1576                 break;
1577           }
1578            case NOISE_MLE:
1579            {
1580                for (q = fit_start; q < fit_end; ++q) {
1581                    (*pdy_dparam_conv)[q][0] = 1.0f;
1582                    yfit[q] += param[0];
1583                    dy[q] = y[q] - yfit[q];
1584                    weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f);
1585                    alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus?
1586                    beta_weight[q] = dy[q] * weight;
1587                    *chisq += (0.0f == y[q])
1588                            ? 2.0 * yfit[q]
1589                            : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1590                }
1591                if (*chisq <= 0.0f) {
1592                    *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1593                }
1594                break;
1595            }
1596            default:
1597            {
1598                return -3;
1599            }
1600        }
1601
1602        // Check if chi square worsened:
1603        if (0.0f != old_chisq && *chisq >= old_chisq) {
1604            // don't bother to set up the matrices for solution
1605            return 0;
1606        }
1607
1608        i_free = 0;
1609        // for all columns
1610        for (i = 0; i < nparam; ++i) {
1611            if (paramfree[i]) {
1612                j_free = 0;
1613                beta_sum = 0.0f;
1614                // row loop, only need to consider lower triangle
1615                for (j = 0; j <= i; ++j) {
1616                    if (paramfree[j]) {
1617                        dot_product = 0.0f;
1618                        if (0 == j_free) {
1619                            // for all data
1620                            for (k = fit_start; k < fit_end; ++k) {
1621                                dy_dparam_k_i = (*pdy_dparam_conv)[k][i];
1622                                dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it.
1623                                beta_sum += dy_dparam_k_i * beta_weight[k];
1624                            }
1625                        }
1626                        else {
1627                            // for all data
1628                            for (k = fit_start; k < fit_end; ++k) {
1629                                dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k];
1630                            }
1631                        } // k loop
1632
1633                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1634                       // if (i_free != j_free) {
1635                       //     // matrix is symmetric
1636                       //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1637                       // }
1638                        ++j_free;
1639                    }
1640                } // j loop
1641                beta[i_free] = beta_sum;
1642                ++i_free;
1643            }
1644            //else printf("param not free %d\n", i);
1645        } // i loop
1646        }
1647
1648        return 0;
1649}
1650
1651/* These two variants, used just before the Marquardt fitting
1652   functions terminate, compute the function values at all points,
1653   whether or not they are being fitted.  (All points are fitted in
1654   the non-instrument response variant.)  They also compute the
1655   residuals y - yfit at all of those points and compute a chi-squared
1656   value which is not modified at small data values in the POISSON
1657   noise models.  They do not calculate alpha or beta. */
1658
1659int GCI_marquardt_compute_fn_final(float x[], float y[], int ndata,
1660                                         noise_type noise, float sig[],
1661                                         float param[], int paramfree[], int nparam,
1662                                         void (*fitfunc)(float, float [], float *, float [], int),
1663                                         float yfit[], float dy[], float *chisq)
1664{
1665        int i, j, mfit;
1666        float sig2i, dy_dparam[MAXFIT];  /* dy_dparam needed for fitfunc */
1667
1668        for (j=0, mfit=0; j<nparam; j++)
1669                if (paramfree[j])
1670                        mfit++;
1671
1672        /* Calculation of the fitting data will depend upon the type of
1673           noise and the type of instrument response */
1674
1675        /* There's no convolution involved in this function.  This is then
1676           fairly straightforward, depending only upon the type of noise
1677           present.  Since we calculate the standard deviation at every
1678           point in a different way depending upon the type of noise, we
1679           will place our switch(noise) around the entire loop. */
1680
1681        switch (noise) {
1682        case NOISE_CONST:
1683                *chisq = 0.0;
1684                /* Summation loop over all data */
1685                for (i=0; i<ndata; i++) {
1686                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1687                        yfit[i] += param[0];
1688                        dy[i] = y[i] - yfit[i];
1689                        /* And find chi^2 */
1690                        *chisq += dy[i] * dy[i];
1691                }
1692
1693                /* Now divide everything by sigma^2 */
1694                sig2i = 1.0 / (sig[0] * sig[0]);
1695                *chisq *= sig2i;
1696                break;
1697
1698        case NOISE_GIVEN:  /* This is essentially the NR version */
1699                *chisq = 0.0;
1700                /* Summation loop over all data */
1701                for (i=0; i<ndata; i++) {
1702                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1703                        yfit[i] += param[0];
1704                        sig2i = 1.0 / (sig[i] * sig[i]);
1705                        dy[i] = y[i] - yfit[i];
1706                        /* And find chi^2 */
1707                        *chisq += dy[i] * dy[i] * sig2i;
1708                }
1709                break;
1710
1711        case NOISE_POISSON_DATA:
1712                *chisq = 0.0;
1713                /* Summation loop over all data */
1714                for (i=0; i<ndata; i++) {
1715                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1716                        yfit[i] += param[0];
1717                        /* we still don't let the variance drop below 1 */
1718                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1719                        dy[i] = y[i] - yfit[i];
1720                        /* And find chi^2 */
1721                        *chisq += dy[i] * dy[i] * sig2i;
1722                }
1723                break;
1724
1725        case NOISE_POISSON_FIT:
1726                *chisq = 0.0;
1727                // Summation loop over all data
1728                for (i=0; i<ndata; i++) {
1729                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1730                        yfit[i] += param[0];
1731                        // we still don't let the variance drop below 1
1732                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1733                        dy[i] = y[i] - yfit[i];
1734                        // And find chi^2
1735                        *chisq += dy[i] * dy[i] * sig2i;
1736                }
1737                break;
1738
1739        case NOISE_MLE:
1740                        *chisq = 0.0;
1741                        /* Summation loop over all data */
1742                        for (i=0; i<ndata; i++) {
1743                                (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1744                                yfit[i] += param[0];
1745//                              dy[i] = y[i] - yfit[i];
1746
1747                                /* And find chi^2 */
1748//                              sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1749//                              *chisq += dy[i] * dy[i] * sig2i;
1750                                if (yfit[i]<=0.0)
1751                                        ; // do nothing
1752                                else if (y[i]==0.0)
1753                                        *chisq += 2.0*yfit[i];   // to avoid NaN from log
1754                                else
1755                                        *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1756                        }
1757                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1758                break;
1759
1760        case NOISE_GAUSSIAN_FIT:
1761                *chisq = 0.0;
1762                // Summation loop over all data
1763                for (i=0; i<ndata; i++) {
1764                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1765                        yfit[i] += param[0];
1766                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1767                        dy[i] = y[i] - yfit[i];
1768                        // And find chi^2
1769                        *chisq += dy[i] * dy[i] * sig2i;
1770                }
1771                break;
1772
1773        default:
1774                return -3;
1775                /* break; */   // (unreachable code)
1776        }
1777
1778        return 0;
1779}
1780
1781
1782/* And this is the variant which handles an instrument response. */
1783/* We assume that the function values are sensible.  Note also that we
1784   have to be careful about which points which include in the
1785   chi-squared calculation.  Also, we are guaranteed that the
1786   initialisation of the convolution arrays has been performed. */
1787
1788int GCI_marquardt_compute_fn_final_instr(float xincr, float y[], int ndata,
1789                                   int fit_start, int fit_end,
1790                                   float instr[], int ninstr,
1791                                   noise_type noise, float sig[],
1792                                   float param[], int paramfree[], int nparam,
1793                                   void (*fitfunc)(float, float [], float *, float [], int),
1794                                   float yfit[], float dy[], float *chisq,
1795                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1796                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1797{
1798        int i, j, mfit, ret;
1799        float sig2i;
1800        float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
1801        int fnvals_len = *pfnvals_len;
1802        int dy_dparam_nparam_size = *pdy_dparam_nparam_size;
1803
1804        /* check the necessary initialisation for safety, bail out if
1805           broken */
1806        if ((fnvals_len < ndata) || (dy_dparam_nparam_size < nparam))
1807                return -1;
1808        fnvals = *pfnvals;
1809        dy_dparam_pure = *pdy_dparam_pure;
1810        dy_dparam_conv = *pdy_dparam_conv;
1811
1812        for (j=0, mfit=0; j<nparam; j++)
1813                if (paramfree[j]) mfit++;
1814
1815        /* Calculation of the fitting data will depend upon the type of
1816           noise and the type of instrument response */
1817
1818        /* Need to calculate unconvolved values all the way down to 0 for
1819           the instrument response case */
1820        if (ninstr > 0) {
1821                if (fitfunc == GCI_multiexp_lambda)
1822                        ret = multiexp_lambda_array(xincr, param, fnvals,
1823                                                                                dy_dparam_pure, ndata, nparam);
1824                else if (fitfunc == GCI_multiexp_tau)
1825                        ret = multiexp_tau_array(xincr, param, fnvals,
1826                                                                         dy_dparam_pure, ndata, nparam);
1827                else if (fitfunc == GCI_stretchedexp)
1828                        ret = stretchedexp_array(xincr, param, fnvals,
1829                                                                         dy_dparam_pure, ndata, nparam);
1830                else
1831                        ret = -1;
1832
1833                if (ret < 0)
1834                        for (i=0; i<ndata; i++)
1835                                (*fitfunc)(xincr*i, param, &fnvals[i],
1836                                                   dy_dparam_pure[i], nparam);
1837
1838                /* OK, we've got to convolve the model fit with the given
1839                   instrument response.  What we'll do here, then, is to
1840                   generate the whole model fit first, then do the convolution
1841                   with the instrument response.  Only after doing that will
1842                   we attempt to calculate the goodness of fit matrices.  This
1843                   means that we will be looping through all of the points
1844                   twice, which is not worth it if there is no convolution
1845                   necessary. */
1846
1847                for (i=0; i<ndata; i++) {
1848                        int convpts;
1849
1850                        /* We wish to find yfit = fnvals * instr, so explicitly:
1851                             yfit[i] = sum_{j=0}^i fnvals[i-j].instr[j]
1852                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1853                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1854                           fnvals[i-j].instr[j]
1855                        */
1856
1857                        /* Zero our adder; don't need to bother with dy_dparam
1858                           stuff here */
1859                        yfit[i] = 0;
1860
1861                        convpts = (ninstr <= i) ? ninstr-1 : i;
1862                        for (j=0; j<=convpts; j++)
1863                                yfit[i] += fnvals[i-j] * instr[j];
1864                }
1865        } else {
1866                /* Can go straight into the final arrays in this case */
1867                if (fitfunc == GCI_multiexp_lambda)
1868                        ret = multiexp_lambda_array(xincr, param, yfit,
1869                                                                                dy_dparam_conv, ndata, nparam);
1870                else if (fitfunc == GCI_multiexp_tau)
1871                        ret = multiexp_tau_array(xincr, param, yfit,
1872                                                                         dy_dparam_conv, ndata, nparam);
1873                else if (fitfunc == GCI_stretchedexp)
1874                        ret = stretchedexp_array(xincr, param, yfit,
1875                                                                         dy_dparam_conv, ndata, nparam);
1876                else
1877                        ret = -1;
1878
1879                if (ret < 0)
1880                        for (i=0; i<ndata; i++)
1881                                (*fitfunc)(xincr*i, param, &yfit[i],
1882                                                   dy_dparam_conv[i], nparam);
1883        }
1884
1885        /* OK, now we've got our (possibly convolved) data, we can do the
1886           rest almost exactly as above. */
1887
1888        switch (noise) {
1889        case NOISE_CONST:
1890                *chisq = 0.0;
1891                /* Summation loop over all data */
1892                for (i=0; i<fit_start; i++) {
1893                        yfit[i] += param[0];
1894                        dy[i] = y[i] - yfit[i];
1895                }
1896                for ( ; i<fit_end; i++) {
1897                        yfit[i] += param[0];
1898                        dy[i] = y[i] - yfit[i];
1899                        /* And find chi^2 */
1900                        *chisq += dy[i] * dy[i];
1901                }
1902                for ( ; i<ndata; i++) {
1903                        yfit[i] += param[0];
1904                        dy[i] = y[i] - yfit[i];
1905                }
1906
1907                /* Now divide chi-squared by sigma^2 */
1908                sig2i = 1.0 / (sig[0] * sig[0]);
1909                *chisq *= sig2i;
1910                break;
1911
1912        case NOISE_GIVEN:  /* This is essentially the NR version */
1913                *chisq = 0.0;
1914                /* Summation loop over all data */
1915                for (i=0; i<fit_start; i++) {
1916                        yfit[i] += param[0];
1917                        dy[i] = y[i] - yfit[i];
1918                }
1919                for ( ; i<fit_end; i++) {
1920                        yfit[i] += param[0];
1921                        dy[i] = y[i] - yfit[i];
1922                        /* And find chi^2 */
1923                        sig2i = 1.0 / (sig[i] * sig[i]);
1924                        *chisq += dy[i] * dy[i] * sig2i;
1925                }
1926                for ( ; i<ndata; i++) {
1927                        yfit[i] += param[0];
1928                        dy[i] = y[i] - yfit[i];
1929                }
1930                break;
1931
1932        case NOISE_POISSON_DATA:
1933                *chisq = 0.0;
1934                /* Summation loop over all data */
1935                for (i=0; i<fit_start; i++) {
1936                        yfit[i] += param[0];
1937                        dy[i] = y[i] - yfit[i];
1938                }
1939                for ( ; i<fit_end; i++) {
1940                        yfit[i] += param[0];
1941                        dy[i] = y[i] - yfit[i];
1942                        /* And find chi^2 */
1943                        /* we still don't let the variance drop below 1 */
1944                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1945                        *chisq += dy[i] * dy[i] * sig2i;
1946                }
1947                for (; i<ndata; i++) {
1948                        yfit[i] += param[0];
1949                        dy[i] = y[i] - yfit[i];
1950                }
1951                break;
1952
1953        case NOISE_POISSON_FIT:
1954                *chisq = 0.0;
1955                // Summation loop over all data
1956                for (i=0; i<fit_start; i++) {
1957                        yfit[i] += param[0];
1958                        dy[i] = y[i] - yfit[i];
1959                }
1960                for ( ; i<fit_end; i++) {
1961                        yfit[i] += param[0];
1962                        dy[i] = y[i] - yfit[i];
1963                        // And find chi^2
1964                        // we still don't let the variance drop below 1
1965                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1966                        *chisq += dy[i] * dy[i] * sig2i;
1967                }
1968                for ( ; i<ndata; i++) {
1969                        yfit[i] += param[0];
1970                        dy[i] = y[i] - yfit[i];
1971                }
1972                break;
1973
1974        case NOISE_MLE:                              // for the final chisq report a normal chisq measure for MLE
1975                *chisq = 0.0;
1976                // Summation loop over all data
1977                for (i=0; i<fit_start; i++) {
1978                        yfit[i] += param[0];
1979                        dy[i] = y[i] - yfit[i];
1980                }
1981                for ( ; i<fit_end; i++) {
1982                        yfit[i] += param[0];
1983                        dy[i] = y[i] - yfit[i];
1984                        // And find chi^2
1985//                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1986                        if (yfit[i]<=0.0)
1987                                ; // do nothing
1988                        else if (y[i]==0.0)
1989                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1990                        else
1991                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1992                }
1993                for ( ; i<ndata; i++) {
1994                        yfit[i] += param[0];
1995                        dy[i] = y[i] - yfit[i];
1996                }
1997                if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1998                break;
1999
2000        case NOISE_GAUSSIAN_FIT:
2001                *chisq = 0.0;
2002                // Summation loop over all data
2003                for (i=0; i<fit_start; i++) {
2004                        yfit[i] += param[0];
2005                        dy[i] = y[i] - yfit[i];
2006                }
2007                for ( ; i<fit_end; i++) {
2008                        yfit[i] += param[0];
2009                        dy[i] = y[i] - yfit[i];
2010                        // And find chi^2
2011                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
2012                        *chisq += dy[i] * dy[i] * sig2i;
2013                }
2014                for ( ; i<ndata; i++) {
2015                        yfit[i] += param[0];
2016                        dy[i] = y[i] - yfit[i];
2017                }
2018                break;
2019
2020        default:
2021                return -3;
2022                /* break; */   // (unreachable code)
2023        }
2024
2025        return 0;
2026}
2027
2028
2029//********************************* GCI_marquardt_fitting_engine **********************************************************************
2030
2031// This returns the number of iterations or negative if an error occurred.
2032// passes all the data to the ecf routine, checks the returned chisq and re-fits if it is of benefit
2033// was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03
2034
2035int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end,
2036                                                float prompt[], int nprompt,
2037                                                noise_type noise, float sig[],
2038                                                float param[], int paramfree[],
2039                                           int nparam, restrain_type restrain,
2040                                           void (*fitfunc)(float, float [], float *, float [], int),
2041                                           float *fitted, float *residuals, float *chisq,
2042                                           float **covar, float **alpha, float **erraxes,
2043                                           float chisq_target, float chisq_delta, int chisq_percent)
2044{
2045        float oldChisq, local_chisq;
2046        int ret, tries=0;
2047
2048        if (ecf_exportParams) ecf_ExportParams_OpenFile ();
2049
2050        // All of the work is done by the ECF module
2051        ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
2052                                                          prompt, nprompt, noise, sig,
2053                                                          param, paramfree, nparam, restrain, fitfunc,
2054                                                          fitted, residuals, covar, alpha, &local_chisq,
2055                                                          chisq_delta, chisq_percent, erraxes);
2056
2057        // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same
2058        // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr()
2059        oldChisq = 3.0e38;
2060        while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<MAXREFITS)
2061        {
2062                oldChisq = local_chisq;
2063                tries++;
2064                ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
2065                                                          prompt, nprompt, noise, sig,
2066                                                          param, paramfree, nparam, restrain, fitfunc,
2067                                                          fitted, residuals, covar, alpha, &local_chisq,
2068                                                          chisq_delta, chisq_percent, erraxes);
2069        }
2070
2071        if (chisq!=NULL) *chisq = local_chisq;
2072
2073        if (ecf_exportParams) ecf_ExportParams_CloseFile ();
2074
2075        return ret;             // summed number of iterations
2076}
2077
2078/* Cleanup function */
2079void GCI_marquardt_cleanup(void)
2080{
2081/*      if (fnvals_len) {
2082                free(fnvals);
2083                GCI_ecf_free_matrix(dy_dparam_pure);
2084                GCI_ecf_free_matrix(dy_dparam_conv);
2085                fnvals_len = 0;
2086                dy_dparam_nparam_size = 0;
2087        }
2088*/
2089}
2090
2091
2092// Emacs settings:
2093// Local variables:
2094// mode: c
2095// c-basic-offset: 4
2096// tab-width: 4
2097// End:
Note: See TracBrowser for help on using the repository browser.