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

Revision 7809, 57.9 KB checked in by paulbarber, 8 years ago (diff)

Better comments on why the convolution is actually correct.
Fixed one possible crash point if there are no params to fit, now just returns the fitted fn with params provided.

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