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

Revision 7803, 57.7 KB checked in by paulbarber, 8 years ago (diff)

Fixed some memory leaks and some compile warnings.

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        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */
894        if (GCI_solve(covar, *pmfit, dparam) != 0)
895                return -1;
896
897        /* Once converged, evaluate covariance matrix */
898        if (*alambda == 0) {
899                if (GCI_marquardt_compute_fn_final_instr(
900                        xincr, y, ndata, fit_start, fit_end,
901                        instr, ninstr, noise, sig,
902                        param, paramfree, nparam, fitfunc,
903                        yfit, dy, chisq,
904                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
905                        pfnvals_len, pdy_dparam_nparam_size) != 0)
906                    return -3;
907
908                for (j=0; j<(*pmfit); j++)
909                    for (k=0; k<(*pmfit); k++)
910                        covar[j][k] = alpha[j][k];
911                GCI_invert(covar, (*pmfit));
912
913
914                if (*pmfit < nparam) {  /* no need to do this otherwise */
915                        GCI_covar_sort(covar, nparam, paramfree, *pmfit);
916                        GCI_covar_sort(alpha, nparam, paramfree, *pmfit);
917                }
918                return 0;
919        }
920
921        /* Did the trial succeed? */
922        for (j=0, l=0; l<nparam; l++)
923                if (paramfree[l])
924                        paramtry[l] = param[l] + dparam[j++];
925
926        if (restrain == ECF_RESTRAIN_DEFAULT)
927                ret = check_ecf_params (paramtry, nparam, fitfunc);
928        else
929                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
930
931        if (ret != 0) {
932                /* Bad parameters, increase alambda and return */
933                *alambda *= 10.0;
934                return 0;
935        }
936
937        if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end,
938                                                                           instr, ninstr, noise, sig,
939                                                                           paramtry, paramfree, nparam, fitfunc,
940                                                                           yfit, dy, covar, dparam,
941                                                                           chisq, *pochisq, *alambda,
942                                                                           pfnvals, pdy_dparam_pure, pdy_dparam_conv,
943                                                                           pfnvals_len, pdy_dparam_nparam_size) != 0)
944                return -2;
945
946        /* Success, accept the new solution */
947        if (*chisq < *pochisq) {
948                *alambda *= 0.1;
949                *pochisq = *chisq;
950                for (j=0; j<(*pmfit); j++) {
951                        for (k=0; k<(*pmfit); k++)
952                                alpha[j][k] = covar[j][k];
953                        beta[j] = dparam[j];
954                }
955                for (l=0; l<nparam; l++)
956                        param[l] = paramtry[l];
957        } else { /* Failure, increase alambda and return */
958                *alambda *= 10.0;
959                *chisq = *pochisq;
960        }
961
962        return 0;
963}
964
965
966/* Used by GCI_marquardt to evaluate the linearised fitting matrix alpha
967   and vector beta and to calculate chi^2.      The equations involved are
968   given in section 15.5 of Numerical Recipes; basically:
969
970     \chi^2(param) = \sum_{i=1}^N ( (y_i-y(x_i;param)) / sigma_i )^2
971
972     beta_k = -1/2 (d/dparam_k)(chi^2)
973
974     alpha_kl = \sum_{i=1}^N (1/sigma_i^2) .
975                  (dy(x_i;param)/dparam_k) . (dy(x_i;param)/dparam_l)
976
977   (where all of the derivatives are partial).
978
979   If an instrument response is provided, we also take account of it
980   now.  We are given that:
981
982     observed(t) = response(t) * instr(t)
983
984   where response(t) is being fitted with fitfunc; it is also trivial to
985   show that (assuming that instr(t) is known and fixed, with no
986   dependencies on the param_k, the parameters being fitted):
987
988     (d/dparam_k) observed(t) = ((d/dparam_k) response(t)) * instr(t)
989
990   so we do not need to alter the response function in any way to
991   determined the fitted convolved response.
992
993   Again there are two variants of this function, corresponding to the
994   two variants of the Marquardt function.
995*/
996int GCI_marquardt_compute_fn(float x[], float y[], int ndata,
997                                         noise_type noise, float sig[],
998                                         float param[], int paramfree[], int nparam,
999                                         void (*fitfunc)(float, float [], float *, float [], int),
1000                                         float yfit[], float dy[],
1001                                         float **alpha, float beta[], float *chisq, float old_chisq,
1002                                         float alambda)
1003{
1004        int i, j, k, l, m, mfit;
1005        float wt, sig2i, y_ymod, dy_dparam[MAXBINS][MAXFIT];
1006        float alpha_weight[MAXBINS];
1007        float beta_weight[MAXBINS];
1008        int q;
1009        float weight;
1010        int i_free;
1011        int j_free;
1012        float dot_product;
1013        float beta_sum;
1014        float dy_dparam_k_i;
1015
1016        for (j=0, mfit=0; j<nparam; j++)
1017                if (paramfree[j])
1018                        mfit++;
1019
1020        *chisq = 0.0f;
1021
1022        switch (noise) {
1023                case NOISE_CONST:
1024                {
1025                        float i_sig_0_squared = 1.0 / (sig[0] * sig[0]);
1026                        for (q = 0; q < ndata; ++q) {
1027                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1028                                yfit[q] += param[0];
1029                                dy_dparam[q][0] = 1.0;
1030                                dy[q] = y[q] - yfit[q];
1031                                weight = i_sig_0_squared;
1032                                alpha_weight[q] = weight; // 1 / (sig[0] * sig[0])
1033                                weight *= dy[q];
1034                                beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0])
1035                                weight *= dy[q];
1036                                *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0])
1037                        }
1038                        break;
1039                }
1040                case NOISE_GIVEN:
1041                {
1042                        for (q = 0; q < ndata; ++q) {
1043                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1044                                yfit[q] += param[0];
1045                                dy_dparam[q][0] = 1.0;
1046                                dy[q] = y[q] - yfit[q];
1047                                weight = 1.0f / (sig[q] * sig[q]);
1048                                alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
1049                                weight *= dy[q];
1050                                beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
1051                                weight *= dy[q];
1052                                *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
1053                        }
1054                        break;
1055                }
1056                case NOISE_POISSON_DATA:
1057                {
1058                        for (q = 0; q < ndata; ++q) {
1059                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1060                                yfit[q] += param[0];
1061                                dy_dparam[q][0] = 1.0;
1062                                dy[q] = y[q] - yfit[q];
1063                                weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
1064                                alpha_weight[q] = weight; // 1 / sig(q)
1065                                weight *= dy[q];
1066                                beta_weight[q] = weight; // dy[q] / sig(q)
1067                                weight *= dy[q];
1068                                *chisq += weight; // (dy[q] * dy[q]) / sig(q)
1069                        }
1070                        break;
1071                }
1072                case NOISE_POISSON_FIT:
1073                {
1074                        for (q = 0; q < ndata; ++q) {
1075                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1076                                yfit[q] += param[0];
1077                                dy_dparam[q][0] = 1.0;
1078                                dy[q] = y[q] - yfit[q];
1079                                weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1080                                alpha_weight[q] = weight; // 1 / sig(q)
1081                                weight *= dy[q];
1082                                beta_weight[q] = weight; // dy(q) / sig(q)
1083                                weight *= dy[q];
1084                                *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1085                        }
1086                        break;
1087                }
1088                case NOISE_GAUSSIAN_FIT:
1089                {
1090                        for (q = 0; q < ndata; ++q) {
1091                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1092                                yfit[q] += param[0];
1093                                dy_dparam[q][0] = 1.0;
1094                                dy[q] = y[q] - yfit[q];
1095                                weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1096                                alpha_weight[q] = weight; // 1 / sig(q)
1097                                weight *= dy[q];
1098                                beta_weight[q] = weight; // dy[q] / sig(q)
1099                                weight *= dy[q];
1100                                *chisq += weight; // dy[q] / sig(q)
1101                        }
1102                        break;
1103                }
1104                case NOISE_MLE:
1105                {
1106                        for (q = 0; q < ndata; ++q) {
1107                                (*fitfunc)(x[q], param, &yfit[q], dy_dparam[q], nparam);
1108                                yfit[q] += param[0];
1109                                dy_dparam[q][0] = 1.0;
1110                                dy[q] = y[q] - yfit[q];
1111                                weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f);
1112                                alpha_weight[q] = weight * y[q] / yfit[q];
1113                                beta_weight[q] = dy[q] * weight;
1114                                if (yfit[q] > 0.0) {
1115                                        *chisq += (0.0f == y[q])
1116                                                        ? 2.0 * yfit[q]
1117                                                        : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1118                                }
1119                        }
1120                        if (*chisq <= 0.0f) {
1121                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1122                        }
1123                        break;
1124                }
1125                default:
1126                {
1127                        return -3;
1128                }
1129        }
1130
1131        // Check if chi square worsened:
1132        if (0.0f != old_chisq && *chisq >= old_chisq) {
1133                // don't bother to set up the matrices for solution
1134                return 0;
1135        }
1136
1137        i_free = 0;
1138        // for all columns
1139        for (i = 0; i < nparam; ++i) {
1140                if (paramfree[i]) {
1141                        j_free = 0;
1142                        beta_sum = 0.0f;
1143                        // row loop, only need to consider lower triangle
1144                        for (j = 0; j <= i; ++j) {
1145                                if (paramfree[j]) {
1146                                        dot_product = 0.0f;
1147                                        if (0 == j_free) { // true only once for each outer loop i
1148                                                // for all data
1149                                                for (k = 0; k < ndata; ++k) {
1150                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k];
1151                                                        beta_sum += dy_dparam[k][i] * beta_weight[k]; // compute beta only once for each i
1152                                                }
1153                                        }
1154                                        else {
1155                                                // for all data
1156                                                for (k = 0; k < ndata; ++k) {
1157                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k];
1158                                                }
1159                                        } // k loop
1160                                               
1161                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; //TODO ARG w/n/b [i][j] more common usage?  row/column
1162                                        // if (i_free != j_free) { //TODO ARG this approach seemed slower at one time; c/b worth retesting:
1163                                        //     / is symmetric
1164                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1165                                        // }
1166                                        ++j_free;
1167                                }
1168                        } // j loop
1169                        beta[i_free] = beta_sum;
1170                        ++i_free;
1171                }
1172        } // i loop
1173       
1174        return 0;
1175}
1176
1177
1178/* And this is the variant which handles an instrument response. */
1179/* We assume that the function values are sensible. */
1180
1181int GCI_marquardt_compute_fn_instr(float xincr, float y[], int ndata,
1182                                   int fit_start, int fit_end,
1183                                   float instr[], int ninstr,
1184                                   noise_type noise, float sig[],
1185                                   float param[], int paramfree[], int nparam,
1186                                   void (*fitfunc)(float, float [], float *, float [], int),
1187                                   float yfit[], float dy[],
1188                                   float **alpha, float beta[], float *chisq, float old_chisq,
1189                                   float alambda,
1190                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1191                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1192{
1193        int i, j, k, l, m, mfit, ret;
1194        float wt, sig2i, y_ymod;
1195        float alpha_weight[MAXBINS];
1196        float beta_weight[MAXBINS];
1197        int q;
1198        float weight;
1199        int i_free;
1200        int j_free;
1201        float dot_product;
1202        float beta_sum;
1203        float dy_dparam_k_i;
1204       
1205        /* Are we initialising? */
1206        // Malloc the arrays that will get used again in this fit via the pointers passed in
1207        // They will be freed by the higher fn that declared them.
1208        if (alambda < 0) {
1209                /* do any necessary initialisation */
1210                if ((*pfnvals_len) < ndata) {  /* we will need this length for
1211                                                                                the final full computation */
1212                        if ((*pfnvals_len)) {
1213                                free((*pfnvals));
1214                                GCI_ecf_free_matrix((*pdy_dparam_pure));
1215                                GCI_ecf_free_matrix((*pdy_dparam_conv));
1216                                (*pfnvals_len) = 0;
1217                                (*pdy_dparam_nparam_size) = 0;
1218                        }
1219                } else if ((*pdy_dparam_nparam_size) < nparam) {
1220                        GCI_ecf_free_matrix((*pdy_dparam_pure));
1221                        GCI_ecf_free_matrix((*pdy_dparam_conv));
1222                        (*pdy_dparam_nparam_size) = 0;
1223                }
1224                if (! (*pfnvals_len)) {
1225                        if (((*pfnvals) = (float *) malloc(ndata * sizeof(float)))
1226                                == NULL)
1227                                return -1;
1228                        (*pfnvals_len) = ndata;
1229                }
1230                if (! (*pdy_dparam_nparam_size)) {
1231                        /* use (*pfnvals_len), not ndata here! */
1232                        if (((*pdy_dparam_pure) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
1233                                /* Can't keep (*pfnvals) around in this case */
1234                                free((*pfnvals));
1235                                (*pfnvals_len) = 0;
1236                                return -1;
1237                        }
1238                        if (((*pdy_dparam_conv) = GCI_ecf_matrix((*pfnvals_len), nparam)) == NULL) {
1239                                /* Can't keep (*pfnvals) around in this case */
1240                                free((*pfnvals));
1241                                free((*pdy_dparam_pure));
1242                                (*pfnvals_len) = 0;
1243                                return -1;
1244                        }
1245                        (*pdy_dparam_nparam_size) = nparam;
1246                }
1247        }
1248
1249        for (j=0, mfit=0; j<nparam; j++)
1250                if (paramfree[j]) mfit++;
1251
1252        /* Calculation of the fitting data will depend upon the type of
1253           noise and the type of instrument response */
1254
1255        /* Need to calculate unconvolved values all the way down to 0 for
1256           the instrument response case */
1257        if (ninstr > 0) {
1258                if (fitfunc == GCI_multiexp_lambda)
1259                        ret = multiexp_lambda_array(xincr, param, (*pfnvals),
1260                                                                                (*pdy_dparam_pure), fit_end, nparam);
1261                else if (fitfunc == GCI_multiexp_tau)
1262                        ret = multiexp_tau_array(xincr, param, (*pfnvals),
1263                                                                         (*pdy_dparam_pure), fit_end, nparam);
1264                else if (fitfunc == GCI_stretchedexp)
1265                        ret = stretchedexp_array(xincr, param, (*pfnvals),
1266                                                                         (*pdy_dparam_pure), fit_end, nparam);
1267                else
1268                        ret = -1;
1269
1270                if (ret < 0)
1271                        for (i=0; i<fit_end; i++)
1272                                (*fitfunc)(xincr*i, param, &(*pfnvals)[i],
1273                                                   (*pdy_dparam_pure)[i], nparam);
1274
1275                /* OK, we've got to convolve the model fit with the given
1276                   instrument response.  What we'll do here, then, is to
1277                   generate the whole model fit first, then do the convolution
1278                   with the instrument response.  Only after doing that will
1279                   we attempt to calculate the goodness of fit matrices.  This
1280                   means that we will be looping through all of the points
1281                   twice, which is not worth it if there is no convolution
1282                   necessary. */
1283
1284                for (i=fit_start; i<fit_end; i++) {
1285                        int convpts;
1286
1287                        /* We wish to find yfit = (*pfnvals) * instr, so explicitly:
1288                             yfit[i] = sum_{j=0}^i (*pfnvals)[i-j].instr[j]
1289                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1290                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1291                           (*pfnvals)[i-j].instr[j]
1292                        */
1293
1294                        /* Zero our adders */
1295                        yfit[i] = 0;
1296                        for (k=1; k<nparam; k++)
1297                                (*pdy_dparam_conv)[i][k] = 0;
1298
1299                        convpts = (ninstr <= i) ? ninstr-1 : i;
1300                        for (j=0; j<=convpts; j++) {
1301                                yfit[i] += (*pfnvals)[i-j] * instr[j];
1302                                for (k=1; k<nparam; k++)
1303                                        (*pdy_dparam_conv)[i][k] += (*pdy_dparam_pure)[i-j][k] * instr[j];
1304                        }
1305                }
1306        } else {
1307                /* Can go straight into the final arrays in this case */
1308                if (fitfunc == GCI_multiexp_lambda)
1309                        ret = multiexp_lambda_array(xincr, param, yfit,
1310                                                                                (*pdy_dparam_conv), fit_end, nparam);
1311                else if (fitfunc == GCI_multiexp_tau)
1312                        ret = multiexp_tau_array(xincr, param, yfit,
1313                                                                         (*pdy_dparam_conv), fit_end, nparam);
1314                else if (fitfunc == GCI_stretchedexp)
1315                        ret = stretchedexp_array(xincr, param, yfit,
1316                                                                         (*pdy_dparam_conv), fit_end, nparam);
1317                else
1318                        ret = -1;
1319
1320                if (ret < 0)
1321                        for (i=0; i<fit_end; i++)
1322                                (*fitfunc)(xincr*i, param, &yfit[i],
1323                                                   (*pdy_dparam_conv)[i], nparam);
1324        }
1325
1326        /* OK, now we've got our (possibly convolved) data, we can do the
1327           rest almost exactly as above. */
1328       
1329        *chisq = 0.0f;
1330       
1331        switch (noise) {
1332                case NOISE_CONST:
1333                {
1334                        for (q = fit_start; q < fit_end; ++q) {
1335                                (*pdy_dparam_conv)[q][0] = 1.0f;
1336                                yfit[q] += param[0];
1337                                dy[q] = y[q] - yfit[q];
1338                                weight = 1.0f / sig[0];
1339                                alpha_weight[q] = weight; // 1 / (sig[0] * sig[0]);
1340                                weight *= dy[q];
1341                                beta_weight[q] = weight; // dy[q] / (sig[0] * sig[0]);
1342                                weight *= dy[q];
1343                                *chisq += weight; // (dy[q] * dy[q]) / (sig[0] * sig[0]);
1344                        }
1345                        break;
1346                }
1347                case NOISE_GIVEN:
1348                {
1349                        for (q = fit_start; q < fit_end; ++q) {
1350                                (*pdy_dparam_conv)[q][0] = 1.0f;
1351                                yfit[q] += param[0];
1352                                dy[q] = y[q] - yfit[q];
1353                                weight = 1.0f / (sig[q] * sig[q]);
1354                                alpha_weight[q] = weight; // 1 / (sig[q] * sig[q])
1355                                weight *= dy[q];
1356                                beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q])
1357                                weight *= dy[q];
1358                                *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q])
1359                        }
1360                        break;
1361                }
1362                case NOISE_POISSON_DATA:
1363                {
1364                        for (q = fit_start; q < fit_end; ++q) {
1365                                (*pdy_dparam_conv)[q][0] = 1.0f;
1366                                yfit[q] += param[0];
1367                                dy[q] = y[q] - yfit[q];
1368                                weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15);
1369                                alpha_weight[q] = weight; // 1 / sig(q)
1370                                weight *= dy[q];
1371                                beta_weight[q] = weight; // dy[q] / sig(q)
1372                                weight *= dy[q];
1373                                *chisq += weight; // (dy[q] * dy[q]) / sig(q)
1374                        }
1375                        break;
1376                }
1377                case NOISE_POISSON_FIT:
1378                {
1379                        for (q = fit_start; q < fit_end; ++q) {
1380                                (*pdy_dparam_conv)[q][0] = 1.0f;
1381                                yfit[q] += param[0];
1382                                dy[q] = y[q] - yfit[q];
1383                                weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15);
1384                                alpha_weight[q] = weight; // 1 / sig(q)
1385                                weight *= dy[q];
1386                                beta_weight[q] = weight; // dy(q) / sig(q)
1387                                weight *= dy[q];
1388                                *chisq += weight; // (dy(q) * dy(q)) / sig(q)
1389                        }
1390                        break;
1391                }
1392                case NOISE_GAUSSIAN_FIT:
1393                {
1394                        for (q = fit_start; q < fit_end; ++q) {
1395                                (*pdy_dparam_conv)[q][0] = 1.0f;
1396                                yfit[q] += param[0];
1397                                dy[q] = y[q] - yfit[q];
1398                                weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f);
1399                                alpha_weight[q] = weight; // 1 / sig(q)
1400                                weight *= dy[q];
1401                                beta_weight[q] = weight; // dy[q] / sig(q)
1402                                weight *= dy[q];
1403                                *chisq += weight; // dy[q] / sig(q)
1404                        }
1405                        break;
1406                }
1407                case NOISE_MLE:
1408                {
1409                        for (q = fit_start; q < fit_end; ++q) {
1410                                (*pdy_dparam_conv)[q][0] = 1.0f;
1411                                yfit[q] += param[0];
1412                                dy[q] = y[q] - yfit[q];
1413                                weight = (yfit[q] > 1 ? 1.0f / yfit[q] : 1.0f);
1414                                alpha_weight[q] = weight * y[q] / yfit[q];
1415                                beta_weight[q] = dy[q] * weight;
1416                                if (yfit[q] > 0.0) {
1417                                        *chisq += (0.0f == y[q])
1418                                                        ? 2.0 * yfit[q]
1419                                                        : 2.0 * (yfit[q] - y[q]) - 2.0 * y[q] * log(yfit[q] / y[q]);
1420                                }
1421                        }
1422                        if (*chisq <= 0.0f) {
1423                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1424                        }
1425                        break;
1426                }
1427                default:
1428                {
1429                        return -3;
1430                }
1431        }
1432
1433        // Check if chi square worsened:
1434        if (0.0f != old_chisq && *chisq >= old_chisq) {
1435                // don't bother to set up the matrices for solution
1436                return 0;
1437        }
1438
1439        i_free = 0;
1440        // for all columns
1441        for (i = 0; i < nparam; ++i) {
1442                if (paramfree[i]) {
1443                        j_free = 0;
1444                        beta_sum = 0.0f;
1445                        // row loop, only need to consider lower triangle
1446                        for (j = 0; j <= i; ++j) {
1447                                if (paramfree[j]) {
1448                                        dot_product = 0.0f;
1449                                        if (0 == j_free) { // true only once for each outer loop i
1450                                                // for all data
1451                                                for (k = fit_start; k < fit_end; ++k) {
1452                                                        dy_dparam_k_i = (*pdy_dparam_conv)[k][i];
1453                                                        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.
1454                                                        beta_sum += dy_dparam_k_i * beta_weight[k];
1455                                                }
1456                                        }
1457                                        else {
1458                                                // for all data
1459                                                for (k = fit_start; k < fit_end; ++k) {
1460                                                        dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k];
1461                                                }
1462                                        } // k loop
1463                                       
1464                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1465                                        // if (i_free != j_free) {
1466                                        //     // matrix is symmetric
1467                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1468                                        // }
1469                                        ++j_free;
1470                                }
1471                        } // j loop
1472                        beta[i_free] = beta_sum;
1473                        ++i_free;
1474                }
1475        } // i loop
1476
1477        return 0;
1478}
1479
1480/* These two variants, used just before the Marquardt fitting
1481   functions terminate, compute the function values at all points,
1482   whether or not they are being fitted.  (All points are fitted in
1483   the non-instrument response variant.)  They also compute the
1484   residuals y - yfit at all of those points and compute a chi-squared
1485   value which is not modified at small data values in the POISSON
1486   noise models.  They do not calculate alpha or beta. */
1487
1488int GCI_marquardt_compute_fn_final(float x[], float y[], int ndata,
1489                                         noise_type noise, float sig[],
1490                                         float param[], int paramfree[], int nparam,
1491                                         void (*fitfunc)(float, float [], float *, float [], int),
1492                                         float yfit[], float dy[], float *chisq)
1493{
1494        int i, j, mfit;
1495        float sig2i, dy_dparam[MAXFIT];  /* dy_dparam needed for fitfunc */
1496
1497        for (j=0, mfit=0; j<nparam; j++)
1498                if (paramfree[j])
1499                        mfit++;
1500
1501        /* Calculation of the fitting data will depend upon the type of
1502           noise and the type of instrument response */
1503
1504        /* There's no convolution involved in this function.  This is then
1505           fairly straightforward, depending only upon the type of noise
1506           present.  Since we calculate the standard deviation at every
1507           point in a different way depending upon the type of noise, we
1508           will place our switch(noise) around the entire loop. */
1509
1510        switch (noise) {
1511        case NOISE_CONST:
1512                *chisq = 0.0;
1513                /* Summation loop over all data */
1514                for (i=0; i<ndata; i++) {
1515                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1516                        yfit[i] += param[0];
1517                        dy[i] = y[i] - yfit[i];
1518                        /* And find chi^2 */
1519                        *chisq += dy[i] * dy[i];
1520                }
1521
1522                /* Now divide everything by sigma^2 */
1523                sig2i = 1.0 / (sig[0] * sig[0]);
1524                *chisq *= sig2i;
1525                break;
1526
1527        case NOISE_GIVEN:  /* This is essentially the NR version */
1528                *chisq = 0.0;
1529                /* Summation loop over all data */
1530                for (i=0; i<ndata; i++) {
1531                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1532                        yfit[i] += param[0];
1533                        sig2i = 1.0 / (sig[i] * sig[i]);
1534                        dy[i] = y[i] - yfit[i];
1535                        /* And find chi^2 */
1536                        *chisq += dy[i] * dy[i] * sig2i;
1537                }
1538                break;
1539
1540        case NOISE_POISSON_DATA:
1541                *chisq = 0.0;
1542                /* Summation loop over all data */
1543                for (i=0; i<ndata; i++) {
1544                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1545                        yfit[i] += param[0];
1546                        /* we still don't let the variance drop below 1 */
1547                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1548                        dy[i] = y[i] - yfit[i];
1549                        /* And find chi^2 */
1550                        *chisq += dy[i] * dy[i] * sig2i;
1551                }
1552                break;
1553
1554        case NOISE_POISSON_FIT:
1555                *chisq = 0.0;
1556                // Summation loop over all data
1557                for (i=0; i<ndata; i++) {
1558                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1559                        yfit[i] += param[0];
1560                        // we still don't let the variance drop below 1
1561                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1562                        dy[i] = y[i] - yfit[i];
1563                        // And find chi^2
1564                        *chisq += dy[i] * dy[i] * sig2i;
1565                }
1566                break;
1567
1568        case NOISE_MLE:
1569                        *chisq = 0.0;
1570                        /* Summation loop over all data */
1571                        for (i=0; i<ndata; i++) {
1572                                (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1573                                yfit[i] += param[0];
1574//                              dy[i] = y[i] - yfit[i];
1575
1576                                /* And find chi^2 */
1577//                              sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1578//                              *chisq += dy[i] * dy[i] * sig2i;
1579                                if (yfit[i]<=0.0)
1580                                        ; // do nothing
1581                                else if (y[i]==0.0)
1582                                        *chisq += 2.0*yfit[i];   // to avoid NaN from log
1583                                else
1584                                        *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1585                        }
1586                        if (*chisq <= 0.0) *chisq = 1.0e38; // don't let chisq=0 through yfit being all -ve
1587                break;
1588
1589        case NOISE_GAUSSIAN_FIT:
1590                *chisq = 0.0;
1591                // Summation loop over all data
1592                for (i=0; i<ndata; i++) {
1593                        (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam);
1594                        yfit[i] += param[0];
1595                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1596                        dy[i] = y[i] - yfit[i];
1597                        // And find chi^2
1598                        *chisq += dy[i] * dy[i] * sig2i;
1599                }
1600                break;
1601
1602        default:
1603                return -3;
1604                /* break; */   // (unreachable code)
1605        }
1606
1607        return 0;
1608}
1609
1610
1611/* And this is the variant which handles an instrument response. */
1612/* We assume that the function values are sensible.  Note also that we
1613   have to be careful about which points which include in the
1614   chi-squared calculation.  Also, we are guaranteed that the
1615   initialisation of the convolution arrays has been performed. */
1616
1617int GCI_marquardt_compute_fn_final_instr(float xincr, float y[], int ndata,
1618                                   int fit_start, int fit_end,
1619                                   float instr[], int ninstr,
1620                                   noise_type noise, float sig[],
1621                                   float param[], int paramfree[], int nparam,
1622                                   void (*fitfunc)(float, float [], float *, float [], int),
1623                                   float yfit[], float dy[], float *chisq,
1624                                        float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
1625                                        int *pfnvals_len, int *pdy_dparam_nparam_size)
1626{
1627        int i, j, mfit, ret;
1628        float sig2i;
1629        float *fnvals, **dy_dparam_pure, **dy_dparam_conv;
1630        int fnvals_len = *pfnvals_len;
1631        int dy_dparam_nparam_size = *pdy_dparam_nparam_size;
1632
1633        /* check the necessary initialisation for safety, bail out if
1634           broken */
1635        if ((fnvals_len < ndata) || (dy_dparam_nparam_size < nparam))
1636                return -1;
1637        fnvals = *pfnvals;
1638        dy_dparam_pure = *pdy_dparam_pure;
1639        dy_dparam_conv = *pdy_dparam_conv;
1640
1641        for (j=0, mfit=0; j<nparam; j++)
1642                if (paramfree[j]) mfit++;
1643
1644        /* Calculation of the fitting data will depend upon the type of
1645           noise and the type of instrument response */
1646
1647        /* Need to calculate unconvolved values all the way down to 0 for
1648           the instrument response case */
1649        if (ninstr > 0) {
1650                if (fitfunc == GCI_multiexp_lambda)
1651                        ret = multiexp_lambda_array(xincr, param, fnvals,
1652                                                                                dy_dparam_pure, ndata, nparam);
1653                else if (fitfunc == GCI_multiexp_tau)
1654                        ret = multiexp_tau_array(xincr, param, fnvals,
1655                                                                         dy_dparam_pure, ndata, nparam);
1656                else if (fitfunc == GCI_stretchedexp)
1657                        ret = stretchedexp_array(xincr, param, fnvals,
1658                                                                         dy_dparam_pure, ndata, nparam);
1659                else
1660                        ret = -1;
1661
1662                if (ret < 0)
1663                        for (i=0; i<ndata; i++)
1664                                (*fitfunc)(xincr*i, param, &fnvals[i],
1665                                                   dy_dparam_pure[i], nparam);
1666
1667                /* OK, we've got to convolve the model fit with the given
1668                   instrument response.  What we'll do here, then, is to
1669                   generate the whole model fit first, then do the convolution
1670                   with the instrument response.  Only after doing that will
1671                   we attempt to calculate the goodness of fit matrices.  This
1672                   means that we will be looping through all of the points
1673                   twice, which is not worth it if there is no convolution
1674                   necessary. */
1675
1676                for (i=0; i<ndata; i++) {
1677                        int convpts;
1678
1679                        /* We wish to find yfit = fnvals * instr, so explicitly:
1680                             yfit[i] = sum_{j=0}^i fnvals[i-j].instr[j]
1681                           But instr[k]=0 for k >= ninstr, so we only need to sum:
1682                             yfit[i] = sum_{j=0}^{min(ninstr-1,i)}
1683                           fnvals[i-j].instr[j]
1684                        */
1685
1686                        /* Zero our adder; don't need to bother with dy_dparam
1687                           stuff here */
1688                        yfit[i] = 0;
1689
1690                        convpts = (ninstr <= i) ? ninstr-1 : i;
1691                        for (j=0; j<=convpts; j++)
1692                                yfit[i] += fnvals[i-j] * instr[j];
1693                }
1694        } else {
1695                /* Can go straight into the final arrays in this case */
1696                if (fitfunc == GCI_multiexp_lambda)
1697                        ret = multiexp_lambda_array(xincr, param, yfit,
1698                                                                                dy_dparam_conv, ndata, nparam);
1699                else if (fitfunc == GCI_multiexp_tau)
1700                        ret = multiexp_tau_array(xincr, param, yfit,
1701                                                                         dy_dparam_conv, ndata, nparam);
1702                else if (fitfunc == GCI_stretchedexp)
1703                        ret = stretchedexp_array(xincr, param, yfit,
1704                                                                         dy_dparam_conv, ndata, nparam);
1705                else
1706                        ret = -1;
1707
1708                if (ret < 0)
1709                        for (i=0; i<ndata; i++)
1710                                (*fitfunc)(xincr*i, param, &yfit[i],
1711                                                   dy_dparam_conv[i], nparam);
1712        }
1713
1714        /* OK, now we've got our (possibly convolved) data, we can do the
1715           rest almost exactly as above. */
1716
1717        switch (noise) {
1718        case NOISE_CONST:
1719                *chisq = 0.0;
1720                /* Summation loop over all data */
1721                for (i=0; i<fit_start; i++) {
1722                        yfit[i] += param[0];
1723                        dy[i] = y[i] - yfit[i];
1724                }
1725                for ( ; i<fit_end; i++) {
1726                        yfit[i] += param[0];
1727                        dy[i] = y[i] - yfit[i];
1728                        /* And find chi^2 */
1729                        *chisq += dy[i] * dy[i];
1730                }
1731                for ( ; i<ndata; i++) {
1732                        yfit[i] += param[0];
1733                        dy[i] = y[i] - yfit[i];
1734                }
1735
1736                /* Now divide chi-squared by sigma^2 */
1737                sig2i = 1.0 / (sig[0] * sig[0]);
1738                *chisq *= sig2i;
1739                break;
1740
1741        case NOISE_GIVEN:  /* This is essentially the NR version */
1742                *chisq = 0.0;
1743                /* Summation loop over all data */
1744                for (i=0; i<fit_start; i++) {
1745                        yfit[i] += param[0];
1746                        dy[i] = y[i] - yfit[i];
1747                }
1748                for ( ; i<fit_end; i++) {
1749                        yfit[i] += param[0];
1750                        dy[i] = y[i] - yfit[i];
1751                        /* And find chi^2 */
1752                        sig2i = 1.0 / (sig[i] * sig[i]);
1753                        *chisq += dy[i] * dy[i] * sig2i;
1754                }
1755                for ( ; i<ndata; i++) {
1756                        yfit[i] += param[0];
1757                        dy[i] = y[i] - yfit[i];
1758                }
1759                break;
1760
1761        case NOISE_POISSON_DATA:
1762                *chisq = 0.0;
1763                /* Summation loop over all data */
1764                for (i=0; i<fit_start; i++) {
1765                        yfit[i] += param[0];
1766                        dy[i] = y[i] - yfit[i];
1767                }
1768                for ( ; i<fit_end; i++) {
1769                        yfit[i] += param[0];
1770                        dy[i] = y[i] - yfit[i];
1771                        /* And find chi^2 */
1772                        /* we still don't let the variance drop below 1 */
1773                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1774                        *chisq += dy[i] * dy[i] * sig2i;
1775                }
1776                for (; i<ndata; i++) {
1777                        yfit[i] += param[0];
1778                        dy[i] = y[i] - yfit[i];
1779                }
1780                break;
1781
1782        case NOISE_POISSON_FIT:
1783                *chisq = 0.0;
1784                // Summation loop over all data
1785                for (i=0; i<fit_start; i++) {
1786                        yfit[i] += param[0];
1787                        dy[i] = y[i] - yfit[i];
1788                }
1789                for ( ; i<fit_end; i++) {
1790                        yfit[i] += param[0];
1791                        dy[i] = y[i] - yfit[i];
1792                        // And find chi^2
1793                        // we still don't let the variance drop below 1
1794                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1795                        *chisq += dy[i] * dy[i] * sig2i;
1796                }
1797                for ( ; i<ndata; i++) {
1798                        yfit[i] += param[0];
1799                        dy[i] = y[i] - yfit[i];
1800                }
1801                break;
1802
1803        case NOISE_MLE:                              // for the final chisq report a normal chisq measure for MLE
1804                *chisq = 0.0;
1805                // Summation loop over all data
1806                for (i=0; i<fit_start; i++) {
1807                        yfit[i] += param[0];
1808                        dy[i] = y[i] - yfit[i];
1809                }
1810                for ( ; i<fit_end; i++) {
1811                        yfit[i] += param[0];
1812                        dy[i] = y[i] - yfit[i];
1813                        // And find chi^2
1814//                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1815                        if (yfit[i]<=0.0)
1816                                ; // do nothing
1817                        else if (y[i]==0.0)
1818                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1819                        else
1820                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1821                }
1822                for ( ; i<ndata; i++) {
1823                        yfit[i] += param[0];
1824                        dy[i] = y[i] - yfit[i];
1825                }
1826                if (*chisq <= 0.0) *chisq = 1.0e38; // don't let chisq=0 through yfit being all -ve
1827                break;
1828
1829        case NOISE_GAUSSIAN_FIT:
1830                *chisq = 0.0;
1831                // Summation loop over all data
1832                for (i=0; i<fit_start; i++) {
1833                        yfit[i] += param[0];
1834                        dy[i] = y[i] - yfit[i];
1835                }
1836                for ( ; i<fit_end; i++) {
1837                        yfit[i] += param[0];
1838                        dy[i] = y[i] - yfit[i];
1839                        // And find chi^2
1840                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1841                        *chisq += dy[i] * dy[i] * sig2i;
1842                }
1843                for ( ; i<ndata; i++) {
1844                        yfit[i] += param[0];
1845                        dy[i] = y[i] - yfit[i];
1846                }
1847                break;
1848
1849        default:
1850                return -3;
1851                /* break; */   // (unreachable code)
1852        }
1853
1854        return 0;
1855}
1856
1857
1858//********************************* GCI_marquardt_fitting_engine **********************************************************************
1859
1860// This returns the number of iterations or negative if an error occurred.
1861// passes all the data to the ecf routine, checks the returned chisq and re-fits if it is of benefit
1862// was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03
1863
1864int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end,
1865                                                float prompt[], int nprompt,
1866                                                noise_type noise, float sig[],
1867                                                float param[], int paramfree[],
1868                                           int nparam, restrain_type restrain,
1869                                           void (*fitfunc)(float, float [], float *, float [], int),
1870                                           float *fitted, float *residuals, float *chisq,
1871                                           float **covar, float **alpha, float **erraxes,
1872                                           float chisq_target, float chisq_delta, int chisq_percent)
1873{
1874        float oldChisq, local_chisq;
1875        int ret, tries=0;
1876
1877        if (ecf_exportParams) ecf_ExportParams_OpenFile ();
1878
1879        // All of the work is done by the ECF module
1880        ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1881                                                          prompt, nprompt, noise, sig,
1882                                                          param, paramfree, nparam, restrain, fitfunc,
1883                                                          fitted, residuals, covar, alpha, &local_chisq,
1884                                                          chisq_delta, chisq_percent, erraxes);
1885
1886        // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same
1887        // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr()
1888        oldChisq = 3.0e38;
1889        while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<MAXREFITS)
1890        {
1891                oldChisq = local_chisq;
1892                tries++;
1893                ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end,
1894                                                          prompt, nprompt, noise, sig,
1895                                                          param, paramfree, nparam, restrain, fitfunc,
1896                                                          fitted, residuals, covar, alpha, &local_chisq,
1897                                                          chisq_delta, chisq_percent, erraxes);
1898        }
1899
1900        if (chisq!=NULL) *chisq = local_chisq;
1901
1902        if (ecf_exportParams) ecf_ExportParams_CloseFile ();
1903
1904        return ret;             // summed number of iterations
1905}
1906
1907/* Cleanup function */
1908void GCI_marquardt_cleanup(void)
1909{
1910/*      if (fnvals_len) {
1911                free(fnvals);
1912                GCI_ecf_free_matrix(dy_dparam_pure);
1913                GCI_ecf_free_matrix(dy_dparam_conv);
1914                fnvals_len = 0;
1915                dy_dparam_nparam_size = 0;
1916        }
1917*/
1918}
1919
1920
1921// Emacs settings:
1922// Local variables:
1923// mode: c
1924// c-basic-offset: 4
1925// tab-width: 4
1926// End:
Note: See TracBrowser for help on using the repository browser.