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

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

Added GPL license to all source code.

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