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

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

Fixed some memory leaks and some compile warnings.

Line 
1/*
2This file is part of the SLIM-curve package for exponential curve fitting of spectral lifetime data.
3
4Copyright (c) 2010, 2011, Gray Institute University of Oxford & UW-Madison LOCI.
5
6    This program is free software: you can redistribute it and/or modify
7    it under the terms of the GNU General Public License as published by
8    the Free Software Foundation, either version 3 of the License, or
9    (at your option) any later version.
10
11    This program is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14    GNU General Public License for more details.
15
16    You should have received a copy of the GNU General Public License
17    along with this program.  If not, see <http://www.gnu.org/licenses/>.
18 */
19
20/* The 2003 version of the ECF library.  This takes account of the fact that we may be
21   handling Poisson noise.
22
23   This file contains functions for global analysis, referring to a
24   little code from the single transient analysis functions.  Utility
25   code is found in EcfUtil.c and single transient analysis code in
26   EcfSingle.c.
27*/
28
29#include <stdio.h>
30#include <stdlib.h>
31#include <math.h>
32#include "EcfInternal.h"
33
34typedef struct {
35        float **P, **Q, ***S;
36} global_matrix;
37
38typedef struct {
39        float *global, *local;
40} global_vector;
41
42/* Predeclarations */
43
44int GCI_alloc_global_matrix(global_matrix *m,
45                                                        int global, int local, int ntrans);
46void GCI_free_global_matrix(global_matrix *m);
47void GCI_copy_global_matrix(global_matrix dest, global_matrix src,
48                                                        int global, int local, int ntrans);
49int GCI_alloc_global_vector(global_vector *v,
50                                                        int global, int local, int ntrans);
51void GCI_free_global_vector(global_vector *v);
52void GCI_copy_global_vector(global_vector dest, global_vector src,
53                                                        int global, int local, int ntrans);
54int GCI_marquardt_global_exps_est_globals_instr(
55                        float xincr, float **trans, int ndata, int ntrans,
56                        int fit_start, int fit_end, float instr[], int ninstr,
57                        noise_type noise, float sig[], int ftype,
58                        float **param, int paramfree[], int nparam, float gparam[],
59                        restrain_type restrain, float chisq_delta,
60                        float fitted[], float residuals[],
61                        float **covar, float **alpha, float *chisq_global);
62int GCI_marquardt_global_exps_est_params_instr(
63                        float xincr, float **trans, int ndata, int ntrans,
64                        int fit_start, int fit_end, float instr[], int ninstr,
65                        noise_type noise, float sig[], int ftype,
66                        float **param, int paramfree[], int nparam, restrain_type restrain, float chisq_delta,
67                        float exp_pure[], float *exp_conv[],
68                        float **fitted, float **residuals,
69                        float **covar, float **alpha, float chisq_trans[],
70                        int drop_bad_transients);
71int GCI_marquardt_global_exps_calculate_exps_instr(
72                        float xincr, int ndata, float instr[], int ninstr,
73                        int ftype, float param[], int nparam,
74                        float exp_pure[], float *exp_conv[]);
75int GCI_marquardt_global_exps_do_fit_single(
76                float xincr, float y[], int ndata, int fit_start, int fit_end,
77                noise_type noise, float sig[], int ftype,
78                float param[], int paramfree[], int nparam, restrain_type restrain, float chisq_delta,
79                float *exp_conv[], float *fitted, float *residuals,
80                float **covar, float **alpha, float *chisq_trans);
81int GCI_marquardt_global_exps_single_step(
82                                float xincr, float y[],
83                                int ndata, int fit_start, int fit_end,
84                                noise_type noise, float sig[], int ftype,
85                                float param[], int paramfree[], int nparam,
86                                restrain_type restrain,
87                                float *exp_conv[],
88                                float yfit[], float dy[],
89                                float **covar, float **alpha, float *chisq,
90                                float *alambda);
91int GCI_marquardt_global_compute_exps_fn(
92                        float xincr, float y[],
93                        int ndata, int fit_start, int fit_end,
94                        noise_type noise, float sig[], int ftype,
95                        float param[], int paramfree[], int nparam,
96                        float *exp_conv[],
97                        float yfit[], float dy[],
98                        float **alpha, float *beta, float *chisq, float old_chisq);
99int GCI_marquardt_global_compute_exps_fn_final(
100                        float xincr, float y[],
101                        int ndata, int fit_start, int fit_end,
102                        noise_type noise, float sig[], int ftype,
103                        float param[], int paramfree[], int nparam,
104                        float *exp_conv[],
105                        float yfit[], float dy[], float *chisq);
106int GCI_marquardt_global_exps_do_fit_instr(
107                float xincr, float **trans, int ndata, int ntrans,
108                int fit_start, int fit_end, float instr[], int ninstr,
109                noise_type noise, float sig[], int ftype,
110                float **param, int paramfree[], int nparam, restrain_type restrain, float chisq_delta,
111                float exp_pure[], float *exp_conv[],
112                float **fitted, float **residuals,
113                float **covar_scratch, float **alpha_scratch,
114                float *chisq_trans, float *chisq_global,
115                int drop_bad_transients);
116int GCI_marquardt_global_exps_global_step(
117                                float xincr, float **trans,
118                                int ndata, int ntrans, int fit_start, int fit_end,
119                                float instr[], int ninstr,
120                                noise_type noise, float sig[], int ftype,
121                                float **param, int paramfree[], int nparam,
122                                restrain_type restrain,
123                                float exp_pure[], float *exp_conv[],
124                                float **yfit, float **dy,
125                                float *chisq_trans, float *chisq_global,
126                                float **alpha_scratch, float *alambda,
127                                int drop_bad_transients);
128int GCI_marquardt_global_compute_global_exps_fn(
129                float xincr, float **trans, int ndata, int ntrans,
130                int fit_start, int fit_end, float instr[], int ninstr,
131                noise_type noise, float sig[], int ftype,
132                float **param, int paramfree[], int nparam,
133                int mfit_global, int mfit_local, int gindex[], int lindex[],
134                float exp_pure[], float *exp_conv[],
135                float **yfit, float **dy, global_matrix alpha, global_vector beta,
136                float **alpha_scratch, float *chisq_trans, float *chisq_global,
137                int drop_bad_transients);
138int GCI_marquardt_global_compute_global_exps_fn_final(
139                float xincr, float **trans, int ndata, int ntrans,
140                int fit_start, int fit_end, float instr[], int ninstr,
141                noise_type noise, float sig[], int ftype,
142                float **param, int paramfree[], int nparam,
143                int mfit_global, int mfit_local, int gindex[], int lindex[],
144                float exp_pure[], float *exp_conv[],
145                float **yfit, float **dy,
146                float *chisq_trans, float *chisq_global, int drop_bad_transients);
147
148int GCI_marquardt_global_solve_eqn(global_matrix A, global_vector b,
149                        int mfit_global, int mfit_local, int ntrans);
150
151int GCI_marquardt_global_generic_do_fit_instr(
152                float xincr, float **trans, int ndata, int ntrans,
153                int fit_start, int fit_end, float instr[], int ninstr,
154                noise_type noise, float sig[],
155                float **param, int paramfree[], int nparam, int gparam[],
156                restrain_type restrain, float chisq_delta,
157                void (*fitfunc)(float, float [], float *, float [], int),
158                float **fitted, float **residuals,
159                float **covar_scratch, float **alpha_scratch,
160                float *chisq_trans, float *chisq_global);
161int GCI_marquardt_global_generic_global_step(
162                                float xincr, float **trans,
163                                int ndata, int ntrans, int fit_start, int fit_end,
164                                float instr[], int ninstr,
165                                noise_type noise, float sig[],
166                                float **param, int paramfree[], int nparam, int gparam[],
167                                restrain_type restrain, float chisq_delta,
168                                void (*fitfunc)(float, float [], float *, float [], int),
169                                float **yfit, float **dy,
170                                float *chisq_trans, float *chisq_global,
171                                float **alpha_scratch, float *alambda);
172int GCI_marquardt_global_compute_global_generic_fn(
173                float xincr, float **trans, int ndata, int ntrans,
174                int fit_start, int fit_end, float instr[], int ninstr,
175                noise_type noise, float sig[],
176                float **param, int paramfree[], int nparam, int gparam[],
177                int mfit_global, int mfit_local, int gindex[], int lindex[],
178                void (*fitfunc)(float, float [], float *, float [], int),
179                float **yfit, float **dy, global_matrix alpha, global_vector beta,
180                float **alpha_scratch, float *chisq_trans, float *chisq_global,
181                float alambda,
182                float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
183                int *pfnvals_len, int *pdy_dparam_nparam_size);
184int GCI_marquardt_global_compute_global_generic_fn_final(
185                float xincr, float **trans, int ndata, int ntrans,
186                int fit_start, int fit_end, float instr[], int ninstr,
187                noise_type noise, float sig[],
188                float **param, int paramfree[], int nparam, int gparam[],
189                int mfit_global, int mfit_local, int gindex[], int lindex[],
190                void (*fitfunc)(float, float [], float *, float [], int),
191                float **yfit, float **dy,
192                float *chisq_trans, float *chisq_global,
193                float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
194                int *pfnvals_len, int *pdy_dparam_nparam_size);
195
196
197/********************************************************************
198
199                                                   GLOBAL ANALYSIS
200
201 ********************************************************************/
202
203/* We only work with the case of multiple transients, all with the
204   same instrument/prompt response, the same xincr, the same number of
205   points, and so on.  The recommended fitting algorithm is a
206   three-step process:
207
208   (1) Sum the transients and use this to get initial estimates for
209       the global parameters we are estimating.
210
211   (2) Fixing these parameters, perform a Marquardt fit on each of the
212       transients.  In our cases, this will be fairly efficient, as we
213       will not need to repeatedly calculate the exponential decay.
214
215   (3) Finally, perform the global fit.  There's lots of interesting
216       maths here which we'll discuss when we get to it.  Note that
217       again we will not need to repeatedly calculate the exponential
218       decays for each transient, which will hopefully make the
219       process significantly faster.
220
221   We provide special code to perform these steps in the case of
222   multiexponential and stretched exponential fits where we are aiming
223   to globally fit all of the taus and the h parameter (in the
224   stretched exponential case); these can be performed far more
225   efficiently than the general case.  We also provide code to perform
226   step (3) in the general case; steps (1) and (2) will have to be
227   handled on a case-by-case basis by the calling code.  We also
228   provide a version of step (3) to handle the case of arbitrary
229   x data with no instrument response.
230*/
231
232/*         *****             UTILITY CODE             *****        */
233
234/* First, we define functions to handle the data structures we will be
235   using later on to store the alpha and covar matrices and the beta
236   etc. vectors for global analysis */
237
238int GCI_alloc_global_matrix(global_matrix *m,
239                                                        int global, int local, int ntrans)
240{
241        if (global <= 0 || local < 0 || ntrans <= 0)
242                return -2;
243
244        if ((m->P = GCI_ecf_matrix(global, global)) == NULL)
245                return -1;
246        if (local > 0) {
247                if ((m->Q = GCI_ecf_matrix(global, ntrans*local)) == NULL) {
248                        GCI_ecf_free_matrix(m->P);
249                        return -1;
250                }
251                if ((m->S = GCI_ecf_matrix_array(ntrans, local, local)) == NULL) {
252                        GCI_ecf_free_matrix(m->P);
253                        GCI_ecf_free_matrix(m->Q);
254                        return -1;
255                }
256        } else {
257                m->Q = NULL;
258                m->S = NULL;
259        }
260
261        return 0;
262}
263
264void GCI_free_global_matrix(global_matrix *m)
265{
266        GCI_ecf_free_matrix(m->P);
267        if (m->Q != NULL) {
268                GCI_ecf_free_matrix(m->Q);
269                GCI_ecf_free_matrix_array(m->S);
270        }
271}
272
273void GCI_copy_global_matrix(global_matrix dest, global_matrix src,
274                                                        int global, int local, int ntrans)
275{
276        int i, j, k;
277
278        for (i=0; i<global; i++)
279                for (j=0; j<global; j++)
280                        dest.P[i][j] = src.P[i][j];
281
282        if (local > 0) {
283                for (i=0; i<global; i++)
284                        for (j=0; j<ntrans*local; j++)
285                                dest.Q[i][j] = src.Q[i][j];
286
287                for (i=0; i<ntrans; i++)
288                        for (j=0; j<local; j++)
289                                for (k=0; k<local; k++)
290                                        dest.S[i][j][k] = src.S[i][j][k];
291        }
292}
293
294
295int GCI_alloc_global_vector(global_vector *v,
296                                                        int global, int local, int ntrans)
297{
298        if (global <= 0 || local < 0 || ntrans <= 0)
299                return -2;
300
301        if ((v->global = (float *) malloc(global * sizeof(float))) == NULL)
302                return -1;
303        if (local > 0) {
304                if ((v->local =
305                         (float *) malloc(ntrans * local * sizeof(float))) == NULL) {
306                        free(v->global);
307                        return -1;
308                }
309        } else {
310                v->local = NULL;
311        }
312
313        return 0;
314}
315
316void GCI_free_global_vector(global_vector *v)
317{
318        free(v->global);
319        if (v->local != NULL)
320                free(v->local);
321}
322
323void GCI_copy_global_vector(global_vector dest, global_vector src,
324                                                        int global, int local, int ntrans)
325{
326        int i;
327
328        for (i=0; i<global; i++)
329                dest.global[i] = src.global[i];
330
331        if (local > 0) {
332                for (i=0; i<ntrans*local; i++)
333                        dest.local[i] = src.local[i];
334        }
335}
336
337
338/*         *****          EXPONENTIALS CODE           *****        */
339
340/* Now the code for performing a global fit for multiexponential taus
341   and stretched exponentials.  This is the function which is called
342   from external programs. */
343
344/* I don't even want to _contemplate_ error axes for this!  It would
345   be computationally very messy, as there would be very large numbers
346   of large vectors involved. */
347
348int GCI_marquardt_global_exps_instr(float xincr, float **trans,
349                                        int ndata, int ntrans, int fit_start, int fit_end,
350                                        float instr[], int ninstr,
351                                        noise_type noise, float sig[], int ftype,
352                                        float **param, int paramfree[], int nparam,
353                                        restrain_type restrain, float chisq_delta,
354                                        float **fitted, float **residuals,
355                                        float chisq_trans[], float *chisq_global, int *df,
356                                        int drop_bad_transients)
357{
358        float **covar, **alpha, *scaled_instr, instrsum;
359        int i, j, ret;
360        int mlocal, mglobal;
361        float gparam[MAXFIT];
362        float *exp_pure, *exp_conv[MAXFIT];
363//      double time;
364
365//      time=Timer();
366
367        /* Some basic parameter checks */
368        if (xincr <= 0) return -1;
369        if (ntrans < 1) return -1;
370        if (ndata < 1)  return -1;
371        if (fit_start < 0 || fit_end > ndata) return -1;
372//      if (ninstr < 1) return -1;
373        if (nparam < 1) return -1;
374
375        switch (ftype) {
376        case FIT_GLOBAL_MULTIEXP:
377                if (nparam % 2 != 1) {
378                        dbgprintf(1, "global fitting: "
379                                          "multiexp needs odd number of parameters\n");
380                        return -1;
381                }
382                break;
383
384        case FIT_GLOBAL_STRETCHEDEXP:
385                if (nparam != 4) {
386                        dbgprintf(1, "global fitting: "
387                                          "stretched exp needs precisely 4 parameters\n");
388                        return -1;
389                }
390                break;
391
392        default:
393                dbgprintf(1, "global fitting: unknown fitting type\n");
394                return -1;
395        }
396
397        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL)
398                return -2;
399
400        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
401                GCI_ecf_free_matrix(covar);
402                return -3;
403        }
404
405        if ((scaled_instr = (float *) malloc(ninstr * sizeof(float))) == NULL && ninstr>1) {
406                GCI_ecf_free_matrix(covar);
407                GCI_ecf_free_matrix(alpha);
408                return -4;
409        }
410
411        /* Also allocate space for the exp_pure and exp_conv arrays */
412        if ((exp_conv[0] = (float *) malloc(nparam * ndata * sizeof(float)))
413                == NULL) {
414                GCI_ecf_free_matrix(covar);
415                GCI_ecf_free_matrix(alpha);
416                free(scaled_instr);
417                return -5;
418        }
419
420        exp_pure = exp_conv[0];
421        for (i=1; i<nparam; i++)
422                exp_conv[i] = exp_conv[0] + i * ndata;
423
424        /* Scale the instrument response */
425        for (i=0, instrsum=0; i<ninstr; i++)
426                instrsum += instr[i];
427        if (instrsum == 0) instrsum=1.0;  //return -6;
428        for (i=0; i<ninstr; i++)
429                scaled_instr[i] = instr[i] / instrsum;
430
431//      printf("that took: %f secs.\n", Timer()-time); time=Timer();
432        dbgprintf(2, "About to enter step (1)\n");
433
434        /* Step (1): estimate the global taus */
435        ret = GCI_marquardt_global_exps_est_globals_instr(
436                        xincr, trans, ndata, ntrans, fit_start, fit_end,
437                        scaled_instr, ninstr, noise, sig,
438                        ftype, param, paramfree, nparam, gparam, restrain, chisq_delta,
439                        fitted[0], residuals[0], covar, alpha, chisq_global);
440
441        if (ret != 0) {
442                dbgprintf(1, "Step (1) failed, ret = %d\n", ret);
443                GCI_ecf_free_matrix(covar);
444                GCI_ecf_free_matrix(alpha);
445                free(scaled_instr);
446                free(exp_conv[0]);
447                return -10 + ret;
448        }
449
450//      printf("that took: %f secs.\n", Timer()-time); time=Timer();
451        /* Copy the estimated global taus to the parameters array */
452
453        switch (ftype)
454        {
455                case FIT_GLOBAL_MULTIEXP:
456                        for (i=2; i<nparam; i+=2)
457                                for (j=0; j<ntrans; j++)
458                                        param[j][i] = gparam[i];
459                break;
460
461                case FIT_GLOBAL_STRETCHEDEXP:
462                        for (i=2; i<nparam; i++)  /* param 2 = tau, param 3 = h */
463                                for (j=0; j<ntrans; j++)
464                                        param[j][i] = gparam[i];
465                break;
466
467                default:
468                        dbgprintf(1, "global_exps_instr: please update me!\n");
469                        GCI_ecf_free_matrix(covar);
470                        GCI_ecf_free_matrix(alpha);
471                        free(scaled_instr);
472                        free(exp_conv[0]);
473                        return -1;
474        }
475
476//      printf("that took: %f secs.\n", Timer()-time); time=Timer();
477        dbgprintf(2, "About to enter step (2)\n");
478
479        /* Step (2): use these taus to estimate initial values for all of
480           the individual transient parameters */
481        ret = GCI_marquardt_global_exps_est_params_instr(
482                        xincr, trans, ndata, ntrans, fit_start, fit_end,
483                        scaled_instr, ninstr, noise, sig, ftype,
484                        param, paramfree, nparam, restrain, chisq_delta,
485                        exp_pure, exp_conv, fitted, residuals, covar, alpha,
486                        chisq_trans, drop_bad_transients);
487
488        if (ret != 0) {
489                dbgprintf(1, "Step (2) failed, ret = %d\n", ret);
490                GCI_ecf_free_matrix(covar);
491                GCI_ecf_free_matrix(alpha);
492                free(scaled_instr);
493                free(exp_conv[0]);
494                return -20 + ret;
495        }
496
497//      printf("that took: %f secs.\n", Timer()-time); time=Timer();
498        dbgprintf(2, "About to enter step (3)\n");
499
500        /* Step (3): now that we have estimates for initial values for all
501           parameters, we can do the global Marquardt fitting.  Note that
502           covar and alpha are only provided as scratch space. */
503        ret = GCI_marquardt_global_exps_do_fit_instr(
504                        xincr, trans, ndata, ntrans, fit_start, fit_end,
505                        scaled_instr, ninstr, noise, sig, ftype,
506                        param, paramfree, nparam, restrain, chisq_delta,
507                        exp_pure, exp_conv, fitted, residuals, covar, alpha,
508                        chisq_trans, chisq_global, drop_bad_transients);
509
510        GCI_ecf_free_matrix(covar);
511        GCI_ecf_free_matrix(alpha);
512        free(scaled_instr);
513        free(exp_conv[0]);
514
515        if (ret < 0) {
516                dbgprintf(1, "Step (3) failed, ret = %d\n", ret);
517                return -30 + ret;
518        }
519
520//      printf("that took: %f secs.\n", Timer()-time); time=Timer();
521        dbgprintf(2, "Step (3) succeeded, ret = %d\n", ret);
522
523        /* Before we return, calculate the number of degrees of freedom */
524        /* The number of degrees of freedom is given by:
525              d.f. = ntrans * ((fit_end - fit_start) - # free local parameters)
526                         - # free global parameters
527        */
528
529        if (drop_bad_transients) {
530                *df = 0;
531                for (i=0; i<ntrans; i++) {
532                        if (chisq_trans[i] > 0)
533                                (*df)++;
534                }
535        } else
536                *df = ntrans;
537
538        mglobal = mlocal = 0;
539
540        switch (ftype)
541        {
542        case FIT_GLOBAL_MULTIEXP:
543                for (i=2; i<nparam; i+=2)
544                if (paramfree[i]) mglobal++;
545                for (i=1; i<nparam; i+=2)
546                        if (paramfree[i]) mlocal++;
547                if (paramfree[0]) mlocal++;
548                break;
549
550        case FIT_GLOBAL_STRETCHEDEXP:
551                if (paramfree[0]) mlocal++;
552                if (paramfree[1]) mlocal++;
553                if (paramfree[2]) mglobal++;
554                if (paramfree[3]) mglobal++;
555                break;
556
557        default:
558                dbgprintf(1, "global_exps_instr: please update me!\n");
559                return -1;
560        }
561
562        *df *= ((fit_end - fit_start) - mlocal);
563        *df -= mglobal;
564
565//      printf("mlocal %d\nmglobal %d\ndf %d\n", mlocal, mglobal, *df);
566
567//      printf("that took: %f secs.\n", Timer()-time); time=Timer();
568        return ret;
569}
570
571
572int GCI_marquardt_global_exps_est_globals_instr(
573                        float xincr, float **trans, int ndata, int ntrans,
574                        int fit_start, int fit_end, float instr[], int ninstr,
575                        noise_type noise, float sig[], int ftype,
576                        float **param, int paramfree[], int nparam, float gparam[],
577                        restrain_type restrain, float chisq_delta,
578                        float fitted[], float residuals[],
579                        float **covar, float **alpha, float *chisq_global)
580{
581        int i, j, ret, nparamfree;
582        float *summed, *tptr;
583        int data_start;
584        float Z, A, tau;
585        void (*fitfunc)(float, float [], float *, float [], int);
586
587        if ((summed = (float *) calloc(ndata, sizeof(float))) == NULL)
588                return -1;
589
590        for (i=0; i<ntrans; i++) {
591                tptr = trans[i];
592                for (j=0; j<ndata; j++)
593                        summed[j] += tptr[j];
594        }
595
596        /* This code is now lifted from fitting.c, appropriately modified */
597        data_start = fit_start + ECF_Find_Float_Max(&summed[fit_start],
598                                                                                                fit_end - fit_start, &A);
599//      ret = GCI_triple_integral_instr(xincr, summed, data_start, fit_end,
600//                                                                      instr, ninstr, noise, sig,
601//                                                                      &Z, &A, &tau, NULL, NULL, NULL);
602
603        ret = GCI_triple_integral_fitting_engine(xincr, summed, data_start, fit_end,
604                                                                                        instr, ninstr, noise, sig,
605                                                                                        &Z, &A, &tau, NULL, NULL, NULL, (float)1.5*(fit_end-fit_start-3));
606
607        dbgprintf(3, "In est_globals_instr, triple integral ret = %d\n", ret);
608
609        if (ret < 0) {
610                Z = 0;
611                ECF_Find_Float_Max(&summed[fit_start], fit_end - fit_start, &A);
612                tau = 2.0;
613        }
614
615
616        /* We set gparam[] to be an array which holds initial estimates of
617           the parameters for the _sum_ of all the transients, for those
618           parameters which are fixed and therefore "known"; the rest will
619           be estimated later.  It doesn't matter if we set a few other
620           values, as these will be overwritten later.  We could also
621           merge this switch() with the next one, but then the code would
622           possibly be a little less easy to follow, so we won't. */
623
624        switch (ftype) {
625        case FIT_GLOBAL_MULTIEXP:
626                /* Z */
627                if (! paramfree[0])
628                {
629                        gparam[0] = 0;
630                        for (j=0; j<ntrans; j++) gparam[0] += param[j][0];
631                }
632
633                /* A's */
634                for (i=1; i<nparam; i++)
635                        if (! paramfree[i])
636                        {
637                                gparam[i] = 0;
638                                for (j=0; j<ntrans; j++) gparam[i] += param[j][i];
639                        }
640
641                /* taus last (was first) */
642                for (i=2; i<nparam; i+=2)
643                        gparam[i] = param[0][i];
644
645                break;
646
647        case FIT_GLOBAL_STRETCHEDEXP:
648                /* Z */
649                if (! paramfree[0]) {
650                        gparam[0] = 0;
651                        for (j=0; j<ntrans; j++) gparam[0] += param[j][0];
652                }
653
654                /* A */
655                if (! paramfree[1]) {
656                        gparam[1] = 0;
657                        for (j=0; j<ntrans; j++) gparam[1] += param[j][1];
658                }
659
660                /* tau and h last (were first) */
661                for (i=2; i<nparam; i++)
662                        gparam[i] = param[0][i];
663
664                break;
665
666        default:
667                dbgprintf(1, "global_exps_est_globals_instr: please update me!\n");
668                free(summed);
669                return -1;
670        }
671
672
673        /* Now we can set any non-fixed parameters to meaningful initial
674           estimates */
675        switch (ftype) {
676        case FIT_GLOBAL_MULTIEXP:
677                fitfunc = GCI_multiexp_tau;
678
679                switch (nparam) {
680                case 3:
681                        if (paramfree[0]) gparam[0] = Z;
682                        if (paramfree[1]) gparam[1] = A;
683                        if (paramfree[2]) gparam[2] = tau;
684                        break;
685
686                case 5:
687                        if (paramfree[0]) gparam[0] = Z;
688                        if (paramfree[1]) gparam[1] = A*3/4;
689                        if (paramfree[2]) gparam[2] = tau;
690                        if (paramfree[3]) gparam[3] = A*1/4;
691                        if (paramfree[4]) gparam[4] = tau*2/3;
692                        break;
693
694                default:
695                        if (nparam<7) {
696                                free(summed);
697                                return -2;
698                        }
699                        if (paramfree[0]) gparam[0] = Z;
700                        if (paramfree[1]) gparam[1] = A*3/4;
701                        if (paramfree[2]) gparam[2] = tau;
702                        if (paramfree[3]) gparam[3] = A*1/6;
703                        if (paramfree[4]) gparam[4] = tau*2/3;
704                        if (paramfree[5]) gparam[5] = A*1/6;
705                        if (paramfree[6]) gparam[6] = tau/3;
706                        for (i=7; i<nparam; i+=2) {
707                                if (paramfree[i]) gparam[i] = A/i;
708                                if (paramfree[i+1]) gparam[i+1] = tau/i;
709                        }
710                        break;
711                }
712                break;
713
714        case FIT_GLOBAL_STRETCHEDEXP:
715                fitfunc = GCI_stretchedexp;
716
717                if (paramfree[0]) gparam[0] = Z;
718                if (paramfree[1]) gparam[1] = A;
719                if (paramfree[2]) gparam[2] = tau;
720                if (paramfree[3]) gparam[3] = 1.5;  /* h */
721                break;
722
723        default:
724                dbgprintf(1, "est_globals_instr: please update me!\n");
725                free(summed);
726                return -1;
727        }
728
729        for (i=0, nparamfree=0; i<nparam; i++) if (paramfree[i]) nparamfree++;
730
731        /* Note that the only values in the gparam array which are of
732           interest are the taus and h for stretched exp, so we don't need
733           to rescale Z and the A's back again */
734//      ret = GCI_marquardt_instr(xincr, summed, ndata, fit_start, fit_end,
735//                                                        instr, ninstr, noise, sig,
736//                                                        gparam, paramfree, nparam, restrain, fitfunc,
737//                                                        fitted, residuals, covar, alpha, chisq_global,
738//                                                        0, NULL);
739
740        ret = GCI_marquardt_fitting_engine(xincr, summed, ndata, fit_start, fit_end,
741                                                          instr, ninstr, noise, sig,
742                                                          gparam, paramfree, nparam, restrain, fitfunc,
743                                                          fitted, residuals, chisq_global, covar, alpha,
744                                                          NULL, (float)1.5*(fit_end-fit_start-nparamfree), chisq_delta, 0);
745
746        dbgprintf(3, "In est_globals_instr, marquardt ret = %d\n", ret);
747
748        free(summed);
749
750        if (ret < 0)
751                return -3;
752        else
753                return 0;
754}
755
756
757int GCI_marquardt_global_exps_est_params_instr(
758                        float xincr, float **trans, int ndata, int ntrans,
759                        int fit_start, int fit_end, float instr[], int ninstr,
760                        noise_type noise, float sig[], int ftype,
761                        float **param, int paramfree[], int nparam, restrain_type restrain, float chisq_delta,
762                        float exp_pure[], float *exp_conv[],
763                        float **fitted, float **residuals,
764                        float **covar, float **alpha, float chisq_trans[],
765                        int drop_bad_transients)
766{
767        int i, j, sortkey[MAXFIT], paramfree_local[MAXFIT], tempi, ret;
768        int data_start;
769        float Z, A, tau;
770
771        /* We begin by estimating the non-tau parameters for each
772           transient.  We do this by performing a triple-integral fit and
773           using the Z and A resulting as a basis for our initial
774           parameters.  In the case of multiexponential fitting, we also
775           sort the taus into decreasing order, and assume that the
776           largest tau is the most significant component. */
777
778        // PRB 03/07 Although **fitted and **residuals are provided only one "transient" is required and used, fitted[0] and residuals[0]
779
780        if (ftype ==  FIT_GLOBAL_MULTIEXP) {
781                /* Initialise */
782                for (i=2; i<nparam; i+=2)
783                        sortkey[i] = i;
784
785                /* Bubblesort :-) */
786                for (i=2; i<nparam; i+=2)
787                        for (j=2*(nparam/2); j>i; j-=2)  /* nparam is odd */
788                                if (param[0][sortkey[j]] > param[0][sortkey[j-2]]) {
789                                        tempi = sortkey[j];
790                                        sortkey[j] = sortkey[j-2];
791                                        sortkey[j-2] = tempi;
792                                }
793        }
794
795        dbgprintf(3, "In est_params_instr, parameters initialised to:\n");
796
797        for (i=0; i<ntrans; i++) {
798                /* This code is now lifted from fitting.c, appropriately modified */
799                data_start = fit_start + ECF_Find_Float_Max(&trans[i][fit_start],
800                                                                                                        fit_end - fit_start, &A);
801//              ret = GCI_triple_integral_instr(xincr, trans[i], data_start, fit_end,
802//                                                                              instr, ninstr, noise, sig,
803//                                                                              &Z, &A, &tau, NULL, NULL, NULL);
804
805                ret = GCI_triple_integral_fitting_engine(xincr, trans[i], data_start, fit_end,
806                                                                                        instr, ninstr, noise, sig,
807                                                                                        &Z, &A, &tau, NULL, NULL, NULL, (float)1.5*(fit_end-fit_start-3));
808                if (ret < 0) {
809                        Z = 0;
810                        ECF_Find_Float_Max(&trans[i][fit_start], fit_end - fit_start, &A);
811                }
812
813                switch (ftype) {
814                case FIT_GLOBAL_MULTIEXP:
815                        switch (nparam) {
816                        case 3:
817                                if (paramfree[0]) param[i][0] = Z;
818                                if (paramfree[1]) param[i][1] = A;
819                                break;
820
821                        case 5:
822                                if (paramfree[0]) param[i][0] = Z;
823                                if (paramfree[sortkey[2]-1]) param[i][sortkey[2]-1] = A*3/4;
824                                if (paramfree[sortkey[4]-1]) param[i][sortkey[4]-1] = A*1/4;
825                                break;
826
827                        default:
828                                if (nparam<7) { /* only actually need to do this once */
829                                        return -1;
830                                }
831                                if (paramfree[0]) param[i][0] = Z;
832                                if (paramfree[sortkey[2]-1]) param[i][sortkey[2]-1] = A*3/4;
833                                if (paramfree[sortkey[4]-1]) param[i][sortkey[4]-1] = A*1/6;
834                                if (paramfree[sortkey[6]-1]) param[i][sortkey[6]-1] = A*1/6;
835                                /* this is all pretty meaningless from here on, anyway */
836                                for (j=8; j<nparam; j+=2) {
837                                        if (paramfree[sortkey[j]-1]) param[i][sortkey[j]-1] = A/j;
838                                }
839                                break;
840                        }
841                        break;
842
843                case FIT_GLOBAL_STRETCHEDEXP:
844                        if (paramfree[0]) param[i][0] = Z;
845                        if (paramfree[1]) param[i][1] = A;
846                        break;
847
848                default:
849                        dbgprintf(1, "est_params_instr: please update me!\n");
850                        return -1;
851                }
852
853                if (ECF_debug >= 3) {
854                        for (j=0; j<nparam; j++)
855                                dbgprintf(3, "param[%d][%d] = %.4g\n", i, j, param[i][j]);
856                }
857        }
858
859
860        /* OK, the initial parameters are set up, we now do a Marquardt
861           fit on all of the transients simultaneously to get more decent
862           initial A and Z values.  But we do this manually without
863           recalculating the exponentials repeatedly.  Furthermore, since
864           the instrument response is convolved linearly with the
865           exponentials, we can get by doing the convolution once only as
866           well, making further major time savings. */
867
868        if (GCI_marquardt_global_exps_calculate_exps_instr(
869                        xincr, ndata, instr, ninstr, ftype, param[0], nparam,
870                        exp_pure, exp_conv) != 0)
871                return -2;
872
873        /* Now we can do a Marquardt fit on each of the transients */
874
875        /* Create a paramfree[] array which fixes all of the taus */
876        switch (ftype) {
877        case FIT_GLOBAL_MULTIEXP:
878                paramfree_local[0] = paramfree[0];  /* Z */
879                for (i=1; i<nparam; i+=2)
880                        paramfree_local[i] = paramfree[i];  /* the A's */
881                for (i=2; i<nparam; i+=2)
882                        paramfree_local[i] = 0;   /* the taus */
883                break;
884
885        case FIT_GLOBAL_STRETCHEDEXP:
886                paramfree_local[0] = paramfree[0];  /* Z */
887                paramfree_local[1] = paramfree[1];  /* A */
888                paramfree_local[2] = 0;  /* tau */
889                paramfree_local[3] = 0;  /* h */
890                break;
891
892        default:
893                dbgprintf(1, "est_params_instr: please update me!\n");
894                return -1;
895        }
896
897        dbgprintf(3, "In est_params_instr, after do_fit_single, "
898                          "parameters initialised to:\n");
899        for (i=0; i<ntrans; i++) {
900                ret = GCI_marquardt_global_exps_do_fit_single(
901                                        xincr, trans[i], ndata, fit_start, fit_end,
902                                        noise, sig, ftype,
903                                        param[i], paramfree_local, nparam, restrain, chisq_delta, exp_conv,
904                                        fitted[0], residuals[0], covar, alpha, &chisq_trans[i]);
905
906                if (ret < 0) {
907                        if (drop_bad_transients) {
908                                dbgprintf(2, "In est_params_instr, transient %d gave "
909                                                  "do_fit_single return val %d; dropping it\n",
910                                                  i, ret);
911                                chisq_trans[i] = -1;
912                                continue;
913                        } else {
914                                dbgprintf(1, "In est_params_instr, transient %d gave "
915                                                  "do_fit_single return val %d\n", i, ret);
916                                return -10 + ret;
917                        }
918                }
919
920                /* We try a second time with these parameters if we got
921                   nonsense */
922                if (chisq_trans[i] > 20 * (fit_end - fit_start)) {
923                        ret = GCI_marquardt_global_exps_do_fit_single(
924                                        xincr, trans[i], ndata, fit_start, fit_end,
925                                        noise, sig, ftype,
926                                        param[i], paramfree_local, nparam, restrain, chisq_delta, exp_conv,
927                                        fitted[0], residuals[0], covar, alpha, &chisq_trans[i]);
928
929                        /* Improved? */
930                        if (ret < 0 || chisq_trans[i] > 20 * (fit_end - fit_start)) {
931                                if (drop_bad_transients) {
932                                        dbgprintf(2, "In est_params_instr, transient %d gave "
933                                                          "do_fit_single return val %d, chisq value %.3f; "
934                                                          "dropping it\n", i, ret, chisq_trans[i]);
935                                        chisq_trans[i] = -1;
936                                        continue;
937                                } else {
938                                        dbgprintf(1, "In est_params_instr, transient %d gave "
939                                                          "do_fit_single return val %d, "
940                                                          "chisq value %.3f\n", i, ret, chisq_trans[i]);
941                                        return -10 + ret;
942                                }
943                        }
944
945                        if (ECF_debug >= 3) {
946                                for (j=0; j<nparam; j++)
947                                        dbgprintf(3, "param[%d][%d] = %.4g\n", i, j, param[i][j]);
948                        }
949                }
950        }
951
952        return 0;
953}
954
955
956/* This finds values of exp(-x/tau) and x*exp(-x/tau)/tau^2, which are
957   needed later for the multiexponential case, and finds the
958   equivalents in the stretched exponential case. */
959// should also now handle no instrument response, i.e. instr=NULL
960
961int GCI_marquardt_global_exps_calculate_exps_instr(
962                        float xincr, int ndata, float instr[], int ninstr,
963                        int ftype, float param[], int nparam,
964                        float exp_pure[], float *exp_conv[])
965{
966        int i, j, k;
967        int convpts;
968        double excur;  /* exp(-x/tau) */
969        double exincr; /* exp(-xincr/tau) */
970        float *expi;
971        float ex, lxa, xah, a2inv, a3inv;  /* for stetched exp */
972        double xa, xaincr;
973
974        switch (ftype) {
975        case FIT_GLOBAL_MULTIEXP:
976                /* Not quite the most efficient way to do this, but not bad */
977                /* First we calculate the exp(-x/tau) array */
978                for (i=2; i<nparam; i+=2) {
979                        expi = exp_conv[i];
980
981                        excur = 1.0;
982                        exincr = exp(-xincr/(double)param[i]);
983                        for (j=0; j<ndata; j++) {
984                                exp_pure[j] = excur;
985                                excur *= exincr;
986
987                                /* And convolve the exponentials with the instrument response if possible */
988                                expi[j] = 0;
989                                convpts = (ninstr <= j) ? ninstr-1 : j;
990                                if (convpts<=0 || instr==NULL) expi[j] = exp_pure[j];
991                                else for (k=0; k<=convpts; k++) expi[j] += exp_pure[j-k]*instr[k];
992                        }
993                }
994
995                /* Now we repeat the exercise for x*exp(-x/tau) / tau^2 */
996                for (i=2; i<nparam; i+=2) {
997                        expi = exp_conv[i-1];
998
999                        excur = 1.0 / (param[i]*param[i]);  /* 1/tau^2 */
1000                        exincr = exp(-xincr/(double)param[i]);
1001                        for (j=0; j<ndata; j++) {
1002                                exp_pure[j] = (xincr*i) * excur; /* x*exp(-x/tau) / tau^2 */
1003                                excur *= exincr;
1004
1005                                /* And convolve the exponentials with the instrument response if possible */
1006                                expi[j] = 0;
1007                                convpts = (ninstr <= j) ? ninstr-1 : j;
1008                                if (convpts<=0 || instr==NULL) expi[j] = exp_pure[j];
1009                                else for (k=0; k<=convpts; k++) expi[j] += exp_pure[j-k]*instr[k];
1010                        }
1011                }
1012                break;
1013
1014        case FIT_GLOBAL_STRETCHEDEXP:
1015                /* We start by essentially repeating the stretchedexp_array
1016                   function */
1017                xa=0;
1018                xaincr = xincr / param[2];
1019                a2inv = 1/param[2];
1020                a3inv = 1/param[3];
1021
1022                /* When x=0 */
1023                exp_conv[1][0] = 1.0;
1024                exp_conv[2][0] = exp_conv[3][0] = 0.0;
1025
1026                for (i=1; i<ndata; i++) {
1027                        xa += xaincr;       /* xa = (xincr*i)/param[2] */
1028                        lxa = log(xa);      /* lxa = log(x/param[2]) */
1029                        xah = exp(lxa * a3inv);  /* xah = exp(log(x/param[2])/param[3])
1030                                                    = (x/param[2])^(1/param[3]) */
1031                        exp_conv[1][i] = ex = exp(-xah);
1032                                            /* ex = exp(-(x/param[2])^(1/param[3])) */
1033                        ex *= xah * a3inv;  /* ex = exp(...) * (x/param[2])^(1/param[3]) *
1034                                                      1/param[3] */
1035                        exp_conv[2][i] = ex * a2inv;
1036                        exp_conv[3][i] = ex * lxa * a3inv;
1037                }
1038
1039                if (ninstr>0 && instr!=NULL)   // else exp_conv already contains the pure data
1040                {
1041                        /* Now convolve these with the instrument response */
1042                        for (i=1; i<4; i++)
1043                        {
1044                                expi = exp_conv[i];
1045
1046                                for (j=0; j<ndata; j++)
1047                                        exp_pure[j] = expi[j];  /* save the array in temp storage */
1048
1049                                for (j=0; j<ndata; j++)
1050                                {
1051                                        expi[j] = 0;
1052                                        convpts = (ninstr <= j) ? ninstr-1 : j;
1053                                        for (k=0; k<=convpts; k++)
1054                                                expi[j] += exp_pure[j-k]*instr[k];
1055                                }
1056                        }
1057                }
1058                break;
1059
1060        default:
1061                dbgprintf(1, "calculate_exps_instr: please update me!\n");
1062                return -1;
1063        }
1064
1065        return 0;
1066}
1067
1068
1069/* This is just like the normal GCI_marquardt_instr, except that it
1070   is designed for the multiexp case where we provide the exponential
1071   decays in advance, and where we don't care about error axes */
1072int GCI_marquardt_global_exps_do_fit_single(
1073                float xincr, float y[], int ndata, int fit_start, int fit_end,
1074                noise_type noise, float sig[], int ftype,
1075                float param[], int paramfree[], int nparam, restrain_type restrain,
1076                float chisq_delta, float *exp_conv[], float *fitted, float *residuals,
1077                float **covar, float **alpha, float *chisq)
1078{
1079        float alambda, ochisq;
1080        int mfit;
1081        int i, k, itst, itst_max;
1082
1083        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
1084
1085        mfit = 0;
1086        for (i=0; i<nparam; i++) {
1087                if (paramfree[i])
1088                        mfit++;
1089        }
1090
1091        alambda = -1;
1092        if (GCI_marquardt_global_exps_single_step(
1093                                xincr, y, ndata, fit_start, fit_end,
1094                                noise, sig, ftype, param, paramfree, nparam, restrain,
1095                                exp_conv, fitted, residuals, covar, alpha,
1096                                chisq, &alambda) != 0) {
1097                return -1;
1098        }
1099
1100        k = 1;  /* Iteration counter */
1101        itst = 0;
1102        for (;;) {
1103                k++;
1104                if (k > MAXITERS) {
1105                        return -2;
1106                }
1107
1108                ochisq = *chisq;
1109                if (GCI_marquardt_global_exps_single_step(
1110                                        xincr, y, ndata, fit_start, fit_end,
1111                                        noise, sig, ftype, param, paramfree, nparam, restrain,
1112                                        exp_conv, fitted, residuals, covar, alpha,
1113                                        chisq, &alambda) != 0) {
1114                        return -3;
1115                }
1116
1117                if (*chisq > ochisq)
1118                        itst = 0;
1119                else if (ochisq - *chisq < chisq_delta)
1120                        itst++;
1121
1122                if (itst < itst_max) continue;
1123
1124                /* Endgame; this also handles correcting the chi-squared values */
1125                alambda=0.0;
1126                if (GCI_marquardt_global_exps_single_step(
1127                                        xincr, y, ndata, fit_start, fit_end,
1128                                        noise, sig, ftype, param, paramfree, nparam, restrain,
1129                                        exp_conv, fitted, residuals, covar, alpha,
1130                                        chisq, &alambda) != 0) {
1131                        return -4;
1132                }
1133
1134                return k;  /* We're done now */
1135        }
1136}
1137
1138
1139/* And this one is basically a specialised GCI_marquardt_instr_step */
1140int GCI_marquardt_global_exps_single_step(
1141                                float xincr, float y[],
1142                                int ndata, int fit_start, int fit_end,
1143                                noise_type noise, float sig[], int ftype,
1144                                float param[], int paramfree[], int nparam,
1145                                restrain_type restrain,
1146                                float *exp_conv[],
1147                                float yfit[], float dy[],
1148                                float **covar, float **alpha, float *chisq,
1149                                float *alambda)
1150{
1151        int j, k, l, ret;
1152        static int mfit;
1153        static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT];
1154        static void (*fitfunc)(float, float [], float *, float [], int);
1155
1156        if (nparam > MAXFIT)
1157                return -10;
1158        if (xincr <= 0)
1159                return -11;
1160        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
1161                return -12;
1162
1163        /* Initialisation */
1164        /* We assume we're given sensible starting values for param[] */
1165        if (*alambda < 0.0) {
1166                switch (ftype) {
1167                case FIT_GLOBAL_MULTIEXP:
1168                        fitfunc = GCI_multiexp_tau;
1169                        break;
1170
1171                case FIT_GLOBAL_STRETCHEDEXP:
1172                        fitfunc = GCI_stretchedexp;
1173                        break;
1174
1175                default:
1176                        dbgprintf(1, "exps_single_step: please update me!\n");
1177                        return -1;
1178                }
1179
1180                for (mfit=0, j=0; j<nparam; j++)
1181                        if (paramfree[j])
1182                                mfit++;
1183
1184                if (GCI_marquardt_global_compute_exps_fn(
1185                                        xincr, y, ndata, fit_start, fit_end, noise, sig,
1186                                        ftype, param, paramfree, nparam, exp_conv,
1187                                        yfit, dy, alpha, beta, chisq, 0.0) != 0)
1188                        return -2;
1189
1190                *alambda = 0.001;
1191                ochisq = *chisq;
1192                for (j=0; j<nparam; j++)
1193                        paramtry[j] = param[j];
1194        }
1195
1196        /* Alter linearised fitting matrix by augmenting diagonal elements */
1197        for (j=0; j<mfit; j++) {
1198                for (k=0; k<mfit; k++)
1199                        covar[j][k] = alpha[j][k];
1200                covar[j][j] = alpha[j][j] * (1.0 + (*alambda));
1201                dparam[j] = beta[j];
1202        }
1203
1204        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */
1205        if (GCI_solve(covar, mfit, dparam) != 0)
1206                return -1;
1207
1208        /* Once converged, calculate corrected chi-squared values */
1209        if (*alambda == 0) {
1210                if (GCI_marquardt_global_compute_exps_fn_final(
1211                                        xincr, y, ndata, fit_start, fit_end, noise, sig,
1212                                        ftype, param, paramfree, nparam, exp_conv,
1213                                        yfit, dy, chisq) != 0)
1214                        return -4;
1215                /* Irrelevant */
1216                // if (mfit < nparam) {  /* no need to do this otherwise */
1217                //      GCI_covar_sort(covar, nparam, paramfree, mfit);
1218                //      GCI_covar_sort(alpha, nparam, paramfree, mfit);
1219                // }
1220                return 0;
1221        }
1222
1223        /* Did the trial succeed? */
1224        for (j=0, l=0; l<nparam; l++)
1225                if (paramfree[l])
1226                        paramtry[l] = param[l] + dparam[j++];
1227
1228        if (restrain == ECF_RESTRAIN_DEFAULT)
1229                ret = check_ecf_params (paramtry, nparam, fitfunc);
1230        else
1231                ret = check_ecf_user_params (paramtry, nparam, fitfunc);
1232
1233        if (ret != 0) {
1234                /* Bad parameters, increase alambda and return */
1235                *alambda *= 10.0;
1236                return 0;
1237        }
1238
1239        if (GCI_marquardt_global_compute_exps_fn(
1240                                xincr, y, ndata, fit_start, fit_end, noise, sig,
1241                                ftype, paramtry, paramfree, nparam, exp_conv,
1242                                yfit, dy, covar, dparam, chisq, ochisq) != 0)
1243                return -2;
1244
1245        /* Success, accept the new solution */
1246        if (*chisq < ochisq) {
1247                *alambda *= 0.1;
1248                ochisq = *chisq;
1249                for (j=0; j<mfit; j++) {
1250                        for (k=0; k<mfit; k++)
1251                                alpha[j][k] = covar[j][k];
1252                        beta[j] = dparam[j];
1253                }
1254                for (l=0; l<nparam; l++)
1255                        param[l] = paramtry[l];
1256        } else { /* Failure, increase alambda and return */
1257                *alambda *= 10.0;
1258                *chisq = ochisq;
1259        }
1260
1261        return 0;
1262}
1263
1264
1265/* This is a streamlined GCI_marquardt_compute_fn_instr */
1266int GCI_marquardt_global_compute_exps_fn(
1267                        float xincr, float y[],
1268                        int ndata, int fit_start, int fit_end,
1269                        noise_type noise, float sig[], int ftype,
1270                        float param[], int paramfree[], int nparam,
1271                        float *exp_conv[],
1272                        float yfit[], float dy[],
1273                        float **alpha, float *beta,
1274                        float *chisq, float old_chisq)
1275{
1276        int i, j, k, l, m, mfit;
1277        float wt, sig2i, y_ymod;
1278        float dy_dparam[MAXBINS][MAXFIT];
1279        float alpha_weight[MAXBINS];
1280        float beta_weight[MAXBINS];
1281        float weight;
1282        int i_free;
1283        int j_free;
1284        float dot_product;
1285        float beta_sum;
1286        float dy_dparam_k_i;
1287
1288        for (j=0, mfit=0; j<nparam; j++)
1289                if (paramfree[j])
1290                        mfit++;
1291
1292        *chisq = 0.0;
1293
1294        switch (ftype) {
1295                case FIT_GLOBAL_MULTIEXP:
1296                        switch (noise) {
1297                                case NOISE_CONST:
1298                                        // loop over all data
1299                                        for (i = fit_start; i < fit_end; ++i) {
1300                                                // multi-exponential fit
1301                                                yfit[i] = param[0];  /* Z */
1302                                                dy_dparam[i][0] = 1.0;
1303
1304                                                for (j=1; j<nparam; j+=2) {
1305                                                        yfit[i] += param[j] * exp_conv[j+1][i];
1306                                                        /* A_j . exp(-x/tau_j) */
1307                                                        dy_dparam[i][j] = exp_conv[j+1][i];
1308                                                        /* exp(-x/tau_j) */
1309                                                        dy_dparam[i][j+1] = param[j] * exp_conv[j][i];
1310                                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1311                                                }
1312                                                dy[i] = y[i] - yfit[i];
1313
1314                                                // constant noise
1315                                                weight = 1.0f / sig[0];
1316                                                alpha_weight[i] = weight; // 1 / (sig[0] * sig[0]);
1317                                                weight *= dy[i];
1318                                                beta_weight[i] = weight; // dy[i] / (sig[0] * sig[0]);
1319                                                weight *= dy[i];
1320                                                *chisq += weight; // (dy[i] * dy[i]) / (sig[0] * sig[0]);
1321                                        }
1322                                        break;
1323                                case NOISE_GIVEN:
1324                                        // loop over all data
1325                                        for (i = fit_start; i < fit_end; ++i) {
1326                                                // multi-exponential fit
1327                                                yfit[i] = param[0];  /* Z */
1328                                                dy_dparam[i][0] = 1.0;
1329
1330                                                for (j=1; j<nparam; j+=2) {
1331                                                        yfit[i] += param[j] * exp_conv[j+1][i];
1332                                                        /* A_j . exp(-x/tau_j) */
1333                                                        dy_dparam[i][j] = exp_conv[j+1][i];
1334                                                        /* exp(-x/tau_j) */
1335                                                        dy_dparam[i][j+1] = param[j] * exp_conv[j][i];
1336                                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1337                                                }
1338                                                dy[i] = y[i] - yfit[i];
1339
1340                                                // given noise
1341                                                weight = 1.0f / (sig[i] * sig[i]);
1342                                                alpha_weight[i] = weight; // 1 / (sig[i] * sig[i])
1343                                                weight *= dy[i];
1344                                                beta_weight[i] = weight; // dy[i] / (sig[i] * sig[i])
1345                                                weight *= dy[i];
1346                                                *chisq += weight; // (dy[i] * dy[i]) / (sig[i] * sig[i])
1347                                        }
1348                                        break;
1349                                case NOISE_POISSON_DATA:
1350                                        // loop over all data
1351                                        for (i = fit_start; i < fit_end; ++i) {
1352                                                // multi-exponential fit
1353                                                yfit[i] = param[0];  /* Z */
1354                                                dy_dparam[i][0] = 1.0;
1355
1356                                                for (j=1; j<nparam; j+=2) {
1357                                                        yfit[i] += param[j] * exp_conv[j+1][i];
1358                                                        /* A_j . exp(-x/tau_j) */
1359                                                        dy_dparam[i][j] = exp_conv[j+1][i];
1360                                                        /* exp(-x/tau_j) */
1361                                                        dy_dparam[i][j+1] = param[j] * exp_conv[j][i];
1362                                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1363                                                }
1364                                                dy[i] = y[i] - yfit[i];
1365
1366                                                // poisson noise based on data
1367                                                weight = (y[i] > 15 ? 1.0f / y[i] : 1.0f / 15);
1368                                                alpha_weight[i] = weight; // 1 / sig(i)
1369                                                weight *= dy[i];
1370                                                beta_weight[i] = weight; // dy[i] / sig(i)
1371                                                weight *= dy[i];
1372                                                *chisq += weight; // (dy[i] * dy[i]) / sig(i)
1373                                        }
1374                                        break;
1375                                case NOISE_POISSON_FIT:
1376                                        // loop over all data
1377                                        for (i = fit_start; i < fit_end; ++i) {
1378                                                // multi-exponential fit
1379                                                yfit[i] = param[0];  /* Z */
1380                                                dy_dparam[i][0] = 1.0;
1381
1382                                                for (j=1; j<nparam; j+=2) {
1383                                                        yfit[i] += param[j] * exp_conv[j+1][i];
1384                                                        /* A_j . exp(-x/tau_j) */
1385                                                        dy_dparam[i][j] = exp_conv[j+1][i];
1386                                                        /* exp(-x/tau_j) */
1387                                                        dy_dparam[i][j+1] = param[j] * exp_conv[j][i];
1388                                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1389                                                }
1390                                                dy[i] = y[i] - yfit[i];
1391
1392                                                // poisson noise based on fit
1393                                                weight = (yfit[i] > 15 ? 1.0f / yfit[i] : 1.0f / 15);
1394                                                alpha_weight[i] = weight; // 1 / sig(i)
1395                                                weight *= dy[i];
1396                                                beta_weight[i] = weight; // dy(i) / sig(i)
1397                                                weight *= dy[i];
1398                                                *chisq += weight; // (dy(i) * dy(i)) / sig(i)
1399                                        }
1400                                        break;
1401                                case NOISE_GAUSSIAN_FIT:
1402                                        // loop over all data
1403                                        for (i = fit_start; i < fit_end; ++i) {
1404                                                // multi-exponential fit
1405                                                yfit[i] = param[0];  /* Z */
1406                                                dy_dparam[i][0] = 1.0;
1407
1408                                                for (j=1; j<nparam; j+=2) {
1409                                                        yfit[i] += param[j] * exp_conv[j+1][i];
1410                                                        /* A_j . exp(-x/tau_j) */
1411                                                        dy_dparam[i][j] = exp_conv[j+1][i];
1412                                                        /* exp(-x/tau_j) */
1413                                                        dy_dparam[i][j+1] = param[j] * exp_conv[j][i];
1414                                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1415                                                }
1416                                                dy[i] = y[i] - yfit[i];
1417
1418                                                // gaussian noise based on fit
1419                                                weight = (yfit[i] > 1.0f ? 1.0f / yfit[i] : 1.0f);
1420                                                alpha_weight[i] = weight; // 1 / sig(i)
1421                                                weight *= dy[i];
1422                                                beta_weight[i] = weight; // dy[i] / sig(i)
1423                                                weight *= dy[i];
1424                                                *chisq += weight; // dy[i] / sig(i)
1425                                        }
1426                                        break;
1427                                case NOISE_MLE:
1428                                        // loop over all data
1429                                        for (i = fit_start; i < fit_end; ++i) {
1430                                                // multi-exponential fit
1431                                                yfit[i] = param[0];  /* Z */
1432                                                dy_dparam[i][0] = 1.0;
1433
1434                                                for (j=1; j<nparam; j+=2) {
1435                                                        yfit[i] += param[j] * exp_conv[j+1][i];
1436                                                        /* A_j . exp(-x/tau_j) */
1437                                                        dy_dparam[i][j] = exp_conv[j+1][i];
1438                                                        /* exp(-x/tau_j) */
1439                                                        dy_dparam[i][j+1] = param[j] * exp_conv[j][i];
1440                                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1441                                                }
1442                                                dy[i] = y[i] - yfit[i];
1443
1444
1445                                                // maximum likelihood estimation noise
1446                                                weight = (yfit[i] > 1 ? 1.0f / yfit[i] : 1.0f);
1447                                                alpha_weight[i] = weight * y[i] / yfit[i];
1448                                                beta_weight[i] = dy[i] * weight;
1449                                                if (yfit[i] > 0.0) {
1450                                                        *chisq += (0.0f == y[i])
1451                                                                        ? 2.0 * yfit[i]
1452                                                                        : 2.0 * (yfit[i] - y[i]) - 2.0 * y[i] * log(yfit[i] / y[i]);
1453                                                }
1454                                        }
1455                                        if (*chisq <= 0.0f) {
1456                                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1457                                        }
1458                                        break;
1459                                default:
1460                                        return -3;
1461                        }
1462                        break;
1463                case FIT_GLOBAL_STRETCHEDEXP:
1464                        switch (noise) {
1465                                case NOISE_CONST:
1466                                        // loop over all data
1467                                        for (i = fit_start; i < fit_end; ++i) {
1468                                                // stretched exponential fit
1469                                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1470                                                dy[i] = y[i] - yfit[i];
1471
1472                                                dy_dparam[i][0] = 1.0;
1473                                                dy_dparam[i][1] = exp_conv[1][i];
1474                                                dy_dparam[i][2] = param[1] * exp_conv[2][i];
1475                                                dy_dparam[i][3] = param[1] * exp_conv[3][i];
1476
1477                                                // constant noise
1478                                                weight = 1.0f / sig[0];
1479                                                alpha_weight[i] = weight; // 1 / (sig[0] * sig[0]);
1480                                                weight *= dy[i];
1481                                                beta_weight[i] = weight; // dy[i] / (sig[0] * sig[0]);
1482                                                weight *= dy[i];
1483                                                *chisq += weight; // (dy[i] * dy[i]) / (sig[0] * sig[0]);
1484                                        }
1485                                        break;
1486                                case NOISE_GIVEN:
1487                                        // loop over all data
1488                                        for (i = fit_start; i < fit_end; ++i) {
1489                                                // stretched exponential fit
1490                                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1491                                                dy[i] = y[i] - yfit[i];
1492
1493                                                dy_dparam[i][0] = 1.0;
1494                                                dy_dparam[i][1] = exp_conv[1][i];
1495                                                dy_dparam[i][2] = param[1] * exp_conv[2][i];
1496                                                dy_dparam[i][3] = param[1] * exp_conv[3][i];
1497
1498                                                // given noise
1499                                                weight = 1.0f / (sig[i] * sig[i]);
1500                                                alpha_weight[i] = weight; // 1 / (sig[i] * sig[i])
1501                                                weight *= dy[i];
1502                                                beta_weight[i] = weight; // dy[i] / (sig[i] * sig[i])
1503                                                weight *= dy[i];
1504                                                *chisq += weight; // (dy[i] * dy[i]) / (sig[i] * sig[i])
1505                                        }
1506                                        break;
1507                                case NOISE_POISSON_DATA:
1508                                        // loop over all data
1509                                        for (i = fit_start; i < fit_end; ++i) {
1510                                                // stretched exponential fit
1511                                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1512                                                dy[i] = y[i] - yfit[i];
1513
1514                                                dy_dparam[i][0] = 1.0;
1515                                                dy_dparam[i][1] = exp_conv[1][i];
1516                                                dy_dparam[i][2] = param[1] * exp_conv[2][i];
1517                                                dy_dparam[i][3] = param[1] * exp_conv[3][i];
1518
1519                                                // poisson noise based on data
1520                                                weight = (y[i] > 15 ? 1.0f / y[i] : 1.0f / 15);
1521                                                alpha_weight[i] = weight; // 1 / sig(i)
1522                                                weight *= dy[i];
1523                                                beta_weight[i] = weight; // dy[i] / sig(i)
1524                                                weight *= dy[i];
1525                                                *chisq += weight; // (dy[i] * dy[i]) / sig(i)
1526                                        }
1527                                        break;
1528                                case NOISE_POISSON_FIT:
1529                                        // loop over all data
1530                                        for (i = fit_start; i < fit_end; ++i) {
1531                                                // stretched exponential fit
1532                                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1533                                                dy[i] = y[i] - yfit[i];
1534
1535                                                dy_dparam[i][0] = 1.0;
1536                                                dy_dparam[i][1] = exp_conv[1][i];
1537                                                dy_dparam[i][2] = param[1] * exp_conv[2][i];
1538                                                dy_dparam[i][3] = param[1] * exp_conv[3][i];
1539
1540                                                // poisson noise based on fit
1541                                                weight = (yfit[i] > 15 ? 1.0f / yfit[i] : 1.0f / 15);
1542                                                alpha_weight[i] = weight; // 1 / sig(i)
1543                                                weight *= dy[i];
1544                                                beta_weight[i] = weight; // dy(i) / sig(i)
1545                                                weight *= dy[i];
1546                                                *chisq += weight; // (dy(i) * dy(i)) / sig(i)
1547                                        }
1548                                        break;
1549                                case NOISE_GAUSSIAN_FIT:
1550                                        // loop over all data
1551                                        for (i = fit_start; i < fit_end; ++i) {
1552                                                // stretched exponential fit
1553                                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1554                                                dy[i] = y[i] - yfit[i];
1555
1556                                                dy_dparam[i][0] = 1.0;
1557                                                dy_dparam[i][1] = exp_conv[1][i];
1558                                                dy_dparam[i][2] = param[1] * exp_conv[2][i];
1559                                                dy_dparam[i][3] = param[1] * exp_conv[3][i];
1560
1561                                                // gaussian noise based on fit
1562                                                weight = (yfit[i] > 1.0f ? 1.0f / yfit[i] : 1.0f);
1563                                                alpha_weight[i] = weight; // 1 / sig(i)
1564                                                weight *= dy[i];
1565                                                beta_weight[i] = weight; // dy[i] / sig(i)
1566                                                weight *= dy[i];
1567                                                *chisq += weight; // dy[i] / sig(i)
1568                                        }
1569                                        break;
1570                                case NOISE_MLE:
1571                                        // loop over all data
1572                                        for (i = fit_start; i < fit_end; ++i) {
1573                                                // stretched exponential fit
1574                                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1575                                                dy[i] = y[i] - yfit[i];
1576
1577                                                dy_dparam[i][0] = 1.0;
1578                                                dy_dparam[i][1] = exp_conv[1][i];
1579                                                dy_dparam[i][2] = param[1] * exp_conv[2][i];
1580                                                dy_dparam[i][3] = param[1] * exp_conv[3][i];
1581
1582                                                // maximum likelihood estimation noise
1583                                                weight = (yfit[i] > 1 ? 1.0f / yfit[i] : 1.0f);
1584                                                alpha_weight[i] = weight * y[i] / yfit[i];
1585                                                beta_weight[i] = dy[i] * weight;
1586                                                if (yfit[i] > 0.0) {
1587                                                        *chisq += (0.0f == y[i])
1588                                                                        ? 2.0 * yfit[i]
1589                                                                        : 2.0 * (yfit[i] - y[i]) - 2.0 * y[i] * log(yfit[i] / y[i]);
1590                                                }
1591                                        }
1592                                        if (*chisq <= 0.0f) {
1593                                                *chisq = 1.0e38f; // don't let chisq=0 through yfit being all -ve
1594                                        }
1595                                        break;
1596                                default:
1597                                        return -3;
1598                        }
1599                        break;
1600                default:
1601                        dbgprintf(1, "compute_exps_fn: please update me!\n");
1602                        return -1;
1603        }
1604
1605        // Check if chi square worsened:
1606        if (0.0f != old_chisq && *chisq >= old_chisq) {
1607                // don't bother to set up the matrices for solution
1608                return 0;
1609        }
1610
1611        i_free = 0;
1612        // for all columns
1613        for (i = 0; i < nparam; ++i) {
1614                if (paramfree[i]) {
1615                        j_free = 0;
1616                        beta_sum = 0.0f;
1617                        // row loop, only need to consider lower triangle
1618                        for (j = 0; j <= i; ++j) {
1619                                if (paramfree[j]) {
1620                                        dot_product = 0.0f;
1621                                        if (0 == j_free) { // true only once for each outer loop i
1622                                                // for all data
1623                                                for (k = fit_start; k < fit_end; ++k) {
1624                                                        dy_dparam_k_i = dy_dparam[k][i];
1625                                                        dot_product += dy_dparam_k_i * dy_dparam[k][j] * alpha_weight[k]; //TODO ARG make it [i][k] and just *dy_dparam++ it.
1626                                                        beta_sum += dy_dparam_k_i * beta_weight[k];
1627                                                }
1628                                        }
1629                                        else {
1630                                                // for all data
1631                                                for (k = fit_start; k < fit_end; ++k) {
1632                                                        dot_product += dy_dparam[k][i] * dy_dparam[k][j] * alpha_weight[k];
1633                                                }
1634                                        } // k loop
1635
1636                                        alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product;
1637                                        // if (i_free != j_free) {
1638                                        //     // matrix is symmetric
1639                                        //     alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!!
1640                                        // }
1641                                        ++j_free;
1642                                }
1643                        } // j loop
1644                        beta[i_free] = beta_sum;
1645                        ++i_free;
1646                }
1647        } // i loop
1648
1649        return 0;
1650}
1651
1652
1653/* And this is a final variant which computes the true chi-squared
1654   values and the full fit, as in EcfSingle.c */
1655int GCI_marquardt_global_compute_exps_fn_final(
1656                        float xincr, float y[],
1657                        int ndata, int fit_start, int fit_end,
1658                        noise_type noise, float sig[], int ftype,
1659                        float param[], int paramfree[], int nparam,
1660                        float *exp_conv[],
1661                        float yfit[], float dy[], float *chisq)
1662{
1663        int i, j, mfit;
1664        float sig2i;
1665
1666        for (j=0, mfit=0; j<nparam; j++)
1667                if (paramfree[j])
1668                        mfit++;
1669
1670        /* Calculation of the fitting data will depend upon the type of
1671           noise.  Since there's no convolution involved here, this is
1672           very easy. */
1673
1674        switch (ftype) {
1675        case FIT_GLOBAL_MULTIEXP:
1676                switch (noise) {
1677                case NOISE_CONST:
1678                        *chisq = 0.0;
1679                        /* Summation loop over all data */
1680                        for (i=0; i<ndata; i++) {
1681                                yfit[i] = param[0];  /* Z */
1682                                for (j=1; j<nparam; j+=2) {
1683                                        yfit[i] += param[j] * exp_conv[j+1][i];
1684                                        /* A_j . exp(-x/tau_j) */
1685                                }
1686                                dy[i] = y[i] - yfit[i];
1687
1688                                /* And find chi^2 */
1689                                if (i >= fit_start && i < fit_end)
1690                                        *chisq += dy[i] * dy[i];
1691                        }
1692
1693                        /* Now divide by sigma^2 */
1694                        sig2i = 1.0 / (sig[0] * sig[0]);
1695                        *chisq *= sig2i;
1696                        break;
1697
1698                case NOISE_GIVEN:  /* This is essentially the NR version */
1699                        *chisq = 0.0;
1700                        /* Summation loop over all data */
1701                        for (i=0; i<ndata; i++) {
1702                                yfit[i] = param[0];  /* Z */
1703                                for (j=1; j<nparam; j+=2) {
1704                                        yfit[i] += param[j] * exp_conv[j+1][i];
1705                                        /* A_j . exp(-x/tau_j) */
1706                                }
1707                                dy[i] = y[i] - yfit[i];
1708
1709                                /* And find chi^2 */
1710                                if (i >= fit_start && i < fit_end) {
1711                                        sig2i = 1.0 / (sig[i] * sig[i]);
1712                                        *chisq += dy[i] * dy[i] * sig2i;
1713                                }
1714                        }
1715                        break;
1716
1717                case NOISE_POISSON_DATA:
1718                        *chisq = 0.0;
1719                        /* Summation loop over all data */
1720                        for (i=0; i<ndata; i++) {
1721                                yfit[i] = param[0];  /* Z */
1722                                for (j=1; j<nparam; j+=2) {
1723                                        yfit[i] += param[j] * exp_conv[j+1][i];
1724                                            /* A_j . exp(-x/tau_j) */
1725                                }
1726                                dy[i] = y[i] - yfit[i];
1727
1728                                /* And find chi^2 */
1729                                if (i >= fit_start && i < fit_end) {
1730                                        /* we still don't let the variance drop below 1 */
1731                                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1732                                        *chisq += dy[i] * dy[i] * sig2i;
1733                                }
1734                        }
1735                        break;
1736
1737                case NOISE_POISSON_FIT:
1738                        *chisq = 0.0;
1739                        /* Summation loop over all data */
1740                        for (i=0; i<ndata; i++) {
1741                                yfit[i] = param[0];  /* Z */
1742                                for (j=1; j<nparam; j+=2) {
1743                                        yfit[i] += param[j] * exp_conv[j+1][i];
1744                                            /* A_j . exp(-x/tau_j) */
1745                                }
1746                                dy[i] = y[i] - yfit[i];
1747
1748                                /* And find chi^2 */
1749                                if (i >= fit_start && i < fit_end) {
1750                                        /* we still don't let the variance drop below 1 */
1751                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1752                                        *chisq += dy[i] * dy[i] * sig2i;
1753                                }
1754                        }
1755                        break;
1756
1757                case NOISE_MLE:
1758                        *chisq = 0.0;
1759                        /* Summation loop over all data */
1760                        for (i=0; i<ndata; i++) {
1761                                yfit[i] = param[0];  /* Z */
1762                                for (j=1; j<nparam; j+=2) {
1763                                        yfit[i] += param[j] * exp_conv[j+1][i];
1764                                            /* A_j . exp(-x/tau_j) */
1765                                }
1766                                dy[i] = y[i] - yfit[i];
1767
1768                                /* And find chi^2 */
1769                                if (i >= fit_start && i < fit_end) {
1770//                                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1771//                                      *chisq += dy[i] * dy[i] * sig2i;
1772                                        if (yfit[i]<=0.0)
1773                                                ; // do nothing
1774                                        else if (y[i]==0.0)
1775                                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1776                                        else
1777                                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1778                                }
1779                        }
1780                        if (*chisq <= 0.0) *chisq = 1.0e38; // don't let chisq=0 through yfit being all -ve
1781                        break;
1782
1783
1784                case NOISE_GAUSSIAN_FIT:
1785                        *chisq = 0.0;
1786                        /* Summation loop over all data */
1787                        for (i=0; i<ndata; i++) {
1788                                yfit[i] = param[0];  /* Z */
1789                                for (j=1; j<nparam; j+=2) {
1790                                        yfit[i] += param[j] * exp_conv[j+1][i];
1791                                            /* A_j . exp(-x/tau_j) */
1792                                }
1793                                dy[i] = y[i] - yfit[i];
1794
1795                                /* And find chi^2 */
1796                                if (i >= fit_start && i < fit_end) {
1797                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1798                                        *chisq += dy[i] * dy[i] * sig2i;
1799                                }
1800                        }
1801                        break;
1802
1803                default:
1804                        return -3;
1805                        /* break; */   // (unreachable code)
1806                }
1807                break;
1808
1809        case FIT_GLOBAL_STRETCHEDEXP:
1810                switch (noise) {
1811                case NOISE_CONST:
1812                        *chisq = 0.0;
1813                        /* Summation loop over all data */
1814                        for (i=0; i<ndata; i++) {
1815                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1816                                dy[i] = y[i] - yfit[i];
1817
1818                                /* And find chi^2 */
1819                                if (i >= fit_start && i < fit_end)
1820                                        *chisq += dy[i] * dy[i];
1821                        }
1822
1823                        /* Now divide by sigma^2 */
1824                        sig2i = 1.0 / (sig[0] * sig[0]);
1825                        *chisq *= sig2i;
1826                        break;
1827
1828                case NOISE_GIVEN:  /* This is essentially the NR version */
1829                        *chisq = 0.0;
1830                        /* Summation loop over all data */
1831                        for (i=0; i<ndata; i++) {
1832                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1833                                dy[i] = y[i] - yfit[i];
1834
1835                                /* And find chi^2 */
1836                                if (i >= fit_start && i < fit_end) {
1837                                        sig2i = 1.0 / (sig[i] * sig[i]);
1838                                        *chisq += dy[i] * dy[i] * sig2i;
1839                                }
1840                        }
1841                        break;
1842
1843                case NOISE_POISSON_DATA:
1844                        *chisq = 0.0;
1845                        /* Summation loop over all data */
1846                        for (i=0; i<ndata; i++) {
1847                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1848                                dy[i] = y[i] - yfit[i];
1849
1850                                /* And find chi^2 */
1851                                if (i >= fit_start && i < fit_end) {
1852                                        /* we still don't let the variance drop below 1 */
1853                                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1854                                        *chisq += dy[i] * dy[i] * sig2i;
1855                                }
1856                        }
1857                        break;
1858
1859                case NOISE_POISSON_FIT:
1860                        *chisq = 0.0;
1861                        /* Summation loop over all data */
1862                        for (i=0; i<ndata; i++) {
1863                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1864                                dy[i] = y[i] - yfit[i];
1865
1866                                /* And find chi^2 */
1867                                if (i >= fit_start && i < fit_end) {
1868                                        /* we still don't let the variance drop below 1 */
1869                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1870                                        *chisq += dy[i] * dy[i] * sig2i;
1871                                }
1872                        }
1873                        break;
1874
1875                case NOISE_MLE:
1876                        *chisq = 0.0;
1877                        /* Summation loop over all data */
1878                        for (i=0; i<ndata; i++) {
1879                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1880                                dy[i] = y[i] - yfit[i];
1881
1882                                /* And find chi^2 */
1883                                if (i >= fit_start && i < fit_end) {
1884//                                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1885//                                      *chisq += dy[i] * dy[i] * sig2i;
1886                                        if (yfit[i]<=0.0)
1887                                                ; // do nothing
1888                                        else if (y[i]==0.0)
1889                                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1890                                        else
1891                                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1892                                }
1893                        }
1894                        if (*chisq <= 0.0) *chisq = 1.0e38; // don't let chisq=0 through yfit being all -ve
1895                        break;
1896
1897                case NOISE_GAUSSIAN_FIT:
1898                        *chisq = 0.0;
1899                        /* Summation loop over all data */
1900                        for (i=0; i<ndata; i++) {
1901                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1902                                dy[i] = y[i] - yfit[i];
1903
1904                                /* And find chi^2 */
1905                                if (i >= fit_start && i < fit_end) {
1906                                        /* we still don't let the variance drop below 1 */
1907                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1908                                        *chisq += dy[i] * dy[i] * sig2i;
1909                                }
1910                        }
1911                        break;
1912
1913                default:
1914                        return -3;
1915                        /* break; */   // (unreachable code)
1916                }
1917                break;
1918
1919        default:
1920                dbgprintf(1, "compute_exps_fn: please update me!\n");
1921                return -1;
1922        }
1923
1924        return 0;
1925}
1926
1927
1928/* This one does the actual global fitting for multiexponential taus.
1929   It is basically similar to the above do_fit_single function, except
1930   that now we handle the extra intricacies involved in global
1931   fitting, in particular, the much larger alpha matrix is handled in
1932   a special way. */
1933
1934int GCI_marquardt_global_exps_do_fit_instr(
1935                float xincr, float **trans, int ndata, int ntrans,
1936                int fit_start, int fit_end, float instr[], int ninstr,
1937                noise_type noise, float sig[], int ftype,
1938                float **param, int paramfree[], int nparam, restrain_type restrain,
1939                float chisq_delta, float exp_pure[], float *exp_conv[],
1940                float **fitted, float **residuals,
1941                float **covar_scratch, float **alpha_scratch,
1942                float *chisq_trans, float *chisq_global,
1943                int drop_bad_transients)
1944{
1945        float alambda, ochisq_global, *ochisq_trans;
1946        int i, k, itst, itst_max;
1947        int ret;
1948
1949        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
1950
1951        /* If there are no global parameters being fitted, we simply fit
1952           each local set. */
1953        switch (ftype) {
1954        case FIT_GLOBAL_MULTIEXP:
1955                for (i=2; i<nparam; i+=2) {
1956                        if (paramfree[i]) {
1957                                i = 0; /* sentinel value */
1958                                break;
1959                        }
1960                }
1961                break;
1962
1963        case FIT_GLOBAL_STRETCHEDEXP:
1964                for (i=2; i<nparam; i++) {
1965                        if (paramfree[i]) {
1966                                i = 0; /* sentinel value */
1967                                break;
1968                        }
1969                }
1970                break;
1971
1972        default:
1973                dbgprintf(1, "exps_do_fit_instr: please update me!\n");
1974                return -1;
1975        }
1976
1977        if (i > 0) { /* no globals to fit */
1978                if (GCI_marquardt_global_exps_calculate_exps_instr(
1979                                xincr, ndata, instr, ninstr, ftype, param[0], nparam,
1980                                exp_pure, exp_conv) != 0)
1981                        return -2;
1982
1983                *chisq_global = 0;
1984
1985                for (i=0; i<ntrans; i++) {
1986                        if (drop_bad_transients && chisq_trans[i] < 0)
1987                                continue;
1988
1989                        ret = GCI_marquardt_global_exps_do_fit_single(
1990                                        xincr, trans[i], ndata, fit_start, fit_end,
1991                                        noise, sig, ftype, param[i], paramfree, nparam, restrain,
1992////                                    exp_conv, fitted[i], residuals[i],
1993                                        chisq_delta, exp_conv, fitted[0], residuals[0],
1994                                        covar_scratch, alpha_scratch, &chisq_trans[i]);
1995                        if (ret < 0) {
1996                                if (drop_bad_transients) {
1997                                        dbgprintf(1, "In do_fit_instr, do_fit_single returned %d "
1998                                                          "for transient %d, dropping it\n", ret, i);
1999                                        chisq_trans[i] = -1;
2000                                } else {
2001                                        dbgprintf(1, "In do_fit_instr, do_fit_single returned %d "
2002                                                          "for transient %d\n", ret, i);
2003                                        return -10 + ret;
2004                                }
2005                        } else {
2006                                *chisq_global += chisq_trans[i];
2007                        }
2008                }
2009                return 0;
2010        }
2011
2012        /* If there are no free local variables to fit, we still do the
2013           global fitting, but we have to be a little careful in some of
2014           the later routines */
2015
2016        /* Now allocate all of the arrays we will need. */
2017
2018        if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL)
2019                return -1;
2020
2021        /* We now begin our standard Marquardt loop, with several
2022           modifications */
2023        alambda = -1;
2024        ret = GCI_marquardt_global_exps_global_step(
2025                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2026                                instr, ninstr, noise, sig, ftype,
2027                                param, paramfree, nparam, restrain, exp_pure, exp_conv,
2028                                fitted, residuals, chisq_trans, chisq_global,
2029                                alpha_scratch, &alambda, drop_bad_transients);
2030        if (ret != 0) {
2031                dbgprintf(1, "In do_fit_instr, first global_step returned %d\n", ret);
2032                if (ret != -1) {
2033                        /* Wasn't a memory error, so unallocate arrays */
2034                        alambda = 0.0;
2035                        GCI_marquardt_global_exps_global_step(
2036                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2037                                instr, ninstr, noise, sig, ftype,
2038                                param, paramfree, nparam, restrain, exp_pure, exp_conv,
2039                                fitted, residuals, chisq_trans, chisq_global,
2040                                alpha_scratch, &alambda, drop_bad_transients);
2041                }
2042                free(ochisq_trans);
2043                return ret;
2044        }
2045
2046        k = 1;  /* Iteration counter */
2047        itst = 0;
2048        for (;;) {
2049                dbgprintf(3, "In do_fit_instr, beginning iteration %d:\n", k);
2050                dbgprintf(3, " itst = %d, chisq_global = %.4f\n", itst, *chisq_global);
2051
2052                k++;
2053                if (k > MAXITERS) {
2054                        free(ochisq_trans);
2055                        return -2;
2056                }
2057
2058                ochisq_global = *chisq_global;
2059                for (i=0; i<ntrans; i++)
2060                        ochisq_trans[i] = chisq_trans[i];
2061
2062                ret = GCI_marquardt_global_exps_global_step(
2063                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2064                                        instr, ninstr, noise, sig, ftype,
2065                                        param, paramfree, nparam, restrain, exp_pure, exp_conv,
2066                                        fitted, residuals, chisq_trans, chisq_global,
2067                                        alpha_scratch, &alambda, drop_bad_transients);
2068                if (ret != 0) {
2069                        dbgprintf(1, "In do_fit_instr, second global_step returned %d\n",
2070                                          ret);
2071                        /* Unallocate arrays */
2072                        alambda = 0.0;
2073                        GCI_marquardt_global_exps_global_step(
2074                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2075                                instr, ninstr, noise, sig, ftype,
2076                                param, paramfree, nparam, restrain, exp_pure, exp_conv,
2077                                fitted, residuals, chisq_trans, chisq_global,
2078                                alpha_scratch, &alambda, drop_bad_transients);
2079                        free(ochisq_trans);
2080                        return ret;
2081                }
2082
2083                if (*chisq_global > ochisq_global)
2084                        itst = 0;
2085                else {
2086                        /* Let's try this approach; I really don't know what will
2087                           be best */
2088                        float maxdiff;
2089
2090                        maxdiff = 0.0;
2091                        for (i=0; i<ntrans; i++)
2092                                if (ochisq_trans[i] - chisq_trans[i] > maxdiff)
2093                                        maxdiff = ochisq_trans[i] - chisq_trans[i];
2094
2095                        if (maxdiff < chisq_delta)
2096                                itst++;
2097                        dbgprintf(3, "In do_fit_instr, maxdiff = %.3f:\n", maxdiff);
2098                }
2099
2100                if (itst < itst_max) continue;
2101
2102                /* Endgame */
2103                alambda = 0.0;
2104                ret = GCI_marquardt_global_exps_global_step(
2105                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2106                                        instr, ninstr, noise, sig, ftype,
2107                                        param, paramfree, nparam, restrain, exp_pure, exp_conv,
2108                                        fitted, residuals, chisq_trans, chisq_global,
2109                                        alpha_scratch, &alambda, drop_bad_transients);
2110                if (ret != 0) {
2111                        dbgprintf(1, "In do_fit_instr, final global_step returned %d\n",
2112                                          ret);
2113                        free(ochisq_trans);
2114                        return ret;
2115                }
2116
2117                free(ochisq_trans);
2118                return k;  /* We're done now */
2119        }
2120}
2121
2122
2123/* And this one is basically a specialised GCI_marquardt_instr_step
2124   for the global fitting setup. */
2125int GCI_marquardt_global_exps_global_step(
2126                                float xincr, float **trans,
2127                                int ndata, int ntrans, int fit_start, int fit_end,
2128                                float instr[], int ninstr,
2129                                noise_type noise, float sig[], int ftype,
2130                                float **param, int paramfree[], int nparam,
2131                                restrain_type restrain,
2132                                float exp_pure[], float *exp_conv[],
2133                                float **yfit, float **dy,
2134                                float *chisq_trans, float *chisq_global,
2135                                float **alpha_scratch, float *alambda,
2136                                int drop_bad_transients)
2137{
2138        int i, j, ret;
2139        static global_matrix alpha, covar;
2140        static global_vector beta, dparam;
2141        static float **paramtry;
2142        static int mfit_local, mfit_global;
2143        static int gindex[MAXFIT], lindex[MAXFIT];
2144        static float ochisq_global, *ochisq_trans;
2145        static void (*fitfunc)(float, float [], float *, float [], int);
2146        static int initialised=0;
2147
2148        if (nparam > MAXFIT)
2149                return -10;
2150        if (xincr <= 0)
2151                return -11;
2152        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
2153                return -12;
2154
2155        /* Initialisation */
2156        /* We assume we're given sensible starting values for param[] */
2157        if (*alambda < 0.0) {
2158                /* Start by allocating lots of variables we will need */
2159                mfit_local = mfit_global = 0;
2160
2161                switch (ftype) {
2162                case FIT_GLOBAL_MULTIEXP:
2163                        fitfunc = GCI_multiexp_tau;
2164
2165                        /* We know that all of param[2i], the taus, are the global
2166                           variables, and that the param[2i+1], the A's, are the
2167                           local variables, along with param[0] = Z.  We stored
2168                           the indices of the local and global free variables in
2169                           lindex and gindex respectively. */
2170                        if (paramfree[0]) {
2171                                lindex[mfit_local++] = 0;
2172                        }
2173                        for (i=1; i<nparam; i+=2) {
2174                                if (paramfree[i])
2175                                        lindex[mfit_local++] = i;
2176                        }
2177                        for (i=2; i<nparam; i+=2) {
2178                                if (paramfree[i])
2179                                        gindex[mfit_global++] = i;
2180                        }
2181                        break;
2182
2183                case FIT_GLOBAL_STRETCHEDEXP:
2184                        fitfunc = GCI_stretchedexp;
2185
2186                        if (paramfree[0])
2187                                lindex[mfit_local++] = 0;  /* Z */
2188                        if (paramfree[1])
2189                                lindex[mfit_local++] = 1;  /* A */
2190                        if (paramfree[2])
2191                                gindex[mfit_global++] = 2;  /* tau */
2192                        if (paramfree[3])
2193                                gindex[mfit_global++] = 3;  /* h */
2194                        break;
2195
2196                default:
2197                        dbgprintf(1, "exps_global_step: please update me!\n");
2198                        return -1;
2199                }
2200
2201                if (initialised) {
2202                        GCI_free_global_matrix(&alpha);   GCI_free_global_matrix(&covar);
2203                        GCI_ecf_free_matrix(paramtry);        GCI_free_global_vector(&beta);
2204                        GCI_free_global_vector(&dparam);  free(ochisq_trans);
2205                        initialised = 0;
2206                }
2207
2208                if (GCI_alloc_global_matrix(&alpha, mfit_global, mfit_local, ntrans)
2209                        != 0)
2210                        return -1;
2211
2212                if (GCI_alloc_global_matrix(&covar, mfit_global, mfit_local, ntrans)
2213                        != 0) {
2214                        GCI_free_global_matrix(&alpha);
2215                        return -1;
2216                }
2217
2218                if ((paramtry = GCI_ecf_matrix(ntrans, nparam)) == NULL) {
2219                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2220                        return -1;
2221                }
2222
2223                if (GCI_alloc_global_vector(&beta, mfit_global, mfit_local, ntrans)
2224                        != 0) {
2225                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2226                        GCI_ecf_free_matrix(paramtry);
2227                        return -1;
2228                }
2229
2230                if (GCI_alloc_global_vector(&dparam, mfit_global, mfit_local, ntrans)
2231                        != 0) {
2232                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2233                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
2234                        return -1;
2235                }
2236
2237                if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float)))
2238                        == NULL) {
2239                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2240                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
2241                        GCI_free_global_vector(&dparam);
2242                        return -1;
2243                }
2244
2245                initialised = 1;
2246
2247                if (GCI_marquardt_global_compute_global_exps_fn(
2248                                        xincr, trans, ndata, ntrans,
2249                                        fit_start, fit_end, instr, ninstr, noise, sig, ftype,
2250                                        param, paramfree, nparam, mfit_global, mfit_local,
2251                                        gindex, lindex, exp_pure, exp_conv,
2252                                        yfit, dy, alpha, beta, alpha_scratch,
2253                                        chisq_trans, chisq_global, drop_bad_transients) != 0)
2254                        return -2;
2255
2256                *alambda = 0.001;
2257                ochisq_global = *chisq_global;
2258                for (i=0; i<ntrans; i++)
2259                        ochisq_trans[i] = chisq_trans[i];
2260
2261                /* Initialise paramtry to param */
2262                for (i=0; i<ntrans; i++) {
2263                        for (j=0; j<nparam; j++)
2264                                paramtry[i][j] = param[i][j];
2265                }
2266        }
2267
2268        /* Once converged, evaluate covariance matrix */
2269        if (*alambda == 0) {
2270                if (GCI_marquardt_global_compute_global_exps_fn_final(
2271                                        xincr, trans, ndata, ntrans,
2272                                        fit_start, fit_end, instr, ninstr, noise, sig, ftype,
2273                                        param, paramfree, nparam, mfit_global, mfit_local,
2274                                        gindex, lindex, exp_pure, exp_conv, yfit, dy,
2275                                        chisq_trans, chisq_global, drop_bad_transients) != 0)
2276                        return -3;
2277                /* Don't need to do this here; if we wished to, we'd have to
2278                   move this code (the "if (*alambda == 0)" block) to after
2279                   the Gauss-Jordan call.  We'd also need to rewrite it for
2280                   our situation.... */
2281                //      if (mfit < nparam) {  /* no need to do this otherwise */
2282                //              GCI_covar_sort(covar, nparam, paramfree, mfit);
2283                //              GCI_covar_sort(alpha, nparam, paramfree, mfit);
2284                //      }
2285                GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2286                GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
2287                GCI_free_global_vector(&dparam); free(ochisq_trans);
2288                initialised = 0;
2289                return 0;
2290        }
2291
2292        /* Alter linearised fitting matrix by augmenting diagonal
2293           elements. */
2294        GCI_copy_global_matrix(covar, alpha, mfit_global, mfit_local, ntrans);
2295        GCI_copy_global_vector(dparam, beta, mfit_global, mfit_local, ntrans);
2296        for (j=0; j<mfit_global; j++)
2297                covar.P[j][j] *= 1.0 + (*alambda);
2298        for (i=0; i<ntrans; i++)
2299                for (j=0; j<mfit_local; j++)
2300                        covar.S[i][j][j] *= 1.0 + (*alambda);
2301
2302        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */
2303        if (GCI_marquardt_global_solve_eqn(covar, dparam,
2304                                                                           mfit_global, mfit_local, ntrans) != 0)
2305                return -3;
2306
2307        /* Did the trial succeed?  Modify param by dparam... */
2308        for (i=0; i<ntrans; i++) {
2309                for (j=0; j<mfit_global; j++)
2310                        paramtry[i][gindex[j]] = param[i][gindex[j]] + dparam.global[j];
2311                for (j=0; j<mfit_local; j++)
2312                        paramtry[i][lindex[j]] =
2313                                param[i][lindex[j]] + dparam.local[i*mfit_local + j];
2314        }
2315
2316        for (i=0; i<ntrans; i++) {
2317                if (drop_bad_transients && chisq_trans[i] < 0)
2318                        continue;
2319
2320                if (restrain == ECF_RESTRAIN_DEFAULT)
2321                        ret = check_ecf_params (paramtry[i], nparam, fitfunc);
2322                else
2323                        ret = check_ecf_user_params (paramtry[i], nparam, fitfunc);
2324
2325                if (ret != 0) {
2326                        /* Bad parameters, increase alambda and return */
2327                        *alambda *= 10.0;
2328                        return 0;
2329                }
2330        }
2331
2332        if (GCI_marquardt_global_compute_global_exps_fn(
2333                                xincr, trans, ndata, ntrans,
2334                                fit_start, fit_end, instr, ninstr, noise, sig, ftype,
2335                                paramtry, paramfree, nparam, mfit_global, mfit_local,
2336                                gindex, lindex, exp_pure, exp_conv,
2337                                yfit, dy, covar, dparam, alpha_scratch,
2338                                chisq_trans, chisq_global, drop_bad_transients) != 0)
2339                return -2;
2340
2341        /* Success, accept the new solution */
2342        if (*chisq_global < ochisq_global) {
2343                *alambda *= 0.1;
2344                ochisq_global = *chisq_global;
2345                for (i=0; i<ntrans; i++)
2346                        ochisq_trans[i] = chisq_trans[i];
2347                GCI_copy_global_matrix(alpha, covar, mfit_global, mfit_local, ntrans);
2348                GCI_copy_global_vector(beta, dparam, mfit_global, mfit_local, ntrans);
2349                for (i=0; i<ntrans; i++) {
2350                        for (j=0; j<nparam; j++)
2351                                param[i][j] = paramtry[i][j];
2352                }
2353        } else { /* Failure, increase alambda and return */
2354                *alambda *= 10.0;
2355                *chisq_global = ochisq_global;
2356                for (i=0; i<ntrans; i++)
2357                        chisq_trans[i] = ochisq_trans[i];
2358        }
2359
2360        return 0;
2361}
2362
2363
2364/* Here we use alpha only for scratch space */
2365int GCI_marquardt_global_compute_global_exps_fn(
2366                float xincr, float **trans, int ndata, int ntrans,
2367                int fit_start, int fit_end, float instr[], int ninstr,
2368                noise_type noise, float sig[], int ftype,
2369                float **param, int paramfree[], int nparam,
2370                int mfit_global, int mfit_local, int gindex[], int lindex[],
2371                float exp_pure[], float *exp_conv[],
2372                float **yfit, float **dy, global_matrix alpha, global_vector beta,
2373                float **alpha_scratch, float *chisq_trans, float *chisq_global,
2374                int drop_bad_transients)
2375{
2376        int i, j, k, ret;
2377        float beta_scratch[MAXFIT];  /* scratch space */
2378
2379        /* Calculate the exponential array once only */
2380        if (GCI_marquardt_global_exps_calculate_exps_instr(
2381                        xincr, ndata, instr, ninstr, ftype, param[0], nparam,
2382                        exp_pure, exp_conv) != 0)
2383                return -1;
2384
2385        /* We initialise P and beta_global to zero; the others don't
2386           matter, as they will be totally overwritten */
2387        for (i=0; i<mfit_global; i++) {
2388                for (j=0; j<mfit_global; j++)
2389                        alpha.P[i][j] = 0;
2390                beta.global[i] = 0;
2391        }
2392        *chisq_global = 0.0;
2393
2394        for (i=0; i<ntrans; i++) {
2395                if (drop_bad_transients && chisq_trans[i] < 0) {
2396                        for (j=0; j<mfit_global; j++)
2397                                for (k=0; k<mfit_local; k++)
2398                                        alpha.Q[j][i*mfit_local + k] = 0.0;
2399                        for (j=0; j<mfit_local; j++) {
2400                                /* Make this component of S an identity matrix and of
2401                                   beta zero */
2402                                for (k=0; k<mfit_local; k++)
2403                                        alpha.S[i][j][k] = (j == k) ? 1.0 : 0.0;
2404                                beta.local[i*mfit_local + j] = 0;
2405                        }
2406                        continue;
2407                }
2408
2409                /* This transient is fine! */
2410                ret = GCI_marquardt_global_compute_exps_fn(
2411                                        xincr, trans[i], ndata, fit_start, fit_end, noise, sig,
2412                                        ftype, param[i], paramfree, nparam,
2413////                                    exp_conv, yfit[i], dy[i],
2414                                        exp_conv, yfit[0], dy[0],
2415                                        alpha_scratch, beta_scratch, &chisq_trans[i], 0.0);
2416
2417                if (ret != 0) {
2418                        if (drop_bad_transients) {
2419                                dbgprintf(3, "In compute_global_exps_fn, "
2420                                                  "compute_exps_fn returned %d for transient %d\n",
2421                                                  ret, i);
2422                                chisq_trans[i] = -1;
2423                                continue;
2424                        } else {
2425                                dbgprintf(1, "In compute_global_exps_fn, "
2426                                                  "compute_exps_fn returned %d for transient %d\n",
2427                                                  ret, i);
2428                                return -2;
2429                        }
2430                }
2431
2432                /* So now have to populate alpha and beta with the contents of
2433                   alpha_scratch and beta_scratch. */
2434
2435                for (j=0; j<mfit_global; j++) {
2436                        for (k=0; k<mfit_global; k++)
2437                                alpha.P[j][k] += alpha_scratch[gindex[j]][gindex[k]];
2438                        for (k=0; k<mfit_local; k++)
2439                                alpha.Q[j][i*mfit_local + k] =
2440                                        alpha_scratch[gindex[j]][lindex[k]];
2441                        beta.global[j] += beta_scratch[gindex[j]];
2442                }
2443                for (j=0; j<mfit_local; j++) {
2444                        for (k=0; k<mfit_local; k++)
2445                                alpha.S[i][j][k] = alpha_scratch[lindex[j]][lindex[k]];
2446                        beta.local[i*mfit_local + j] = beta_scratch[lindex[j]];
2447                }
2448
2449                *chisq_global += chisq_trans[i];
2450        }
2451
2452        return 0;
2453}
2454
2455
2456/* The final variant */
2457int GCI_marquardt_global_compute_global_exps_fn_final(
2458                float xincr, float **trans, int ndata, int ntrans,
2459                int fit_start, int fit_end, float instr[], int ninstr,
2460                noise_type noise, float sig[], int ftype,
2461                float **param, int paramfree[], int nparam,
2462                int mfit_global, int mfit_local, int gindex[], int lindex[],
2463                float exp_pure[], float *exp_conv[],
2464                float **yfit, float **dy,
2465                float *chisq_trans, float *chisq_global, int drop_bad_transients)
2466{
2467        int i, ret;
2468
2469        /* Calculate the exponential array once only */
2470        if (GCI_marquardt_global_exps_calculate_exps_instr(
2471                        xincr, ndata, instr, ninstr, ftype, param[0], nparam,
2472                        exp_pure, exp_conv) != 0)
2473                return -1;
2474
2475        *chisq_global = 0.0;
2476
2477        for (i=0; i<ntrans; i++) {
2478                if (drop_bad_transients && chisq_trans[i] < 0)
2479                        continue;
2480
2481                /* This transient is fine! */
2482                ret = GCI_marquardt_global_compute_exps_fn_final(
2483                                        xincr, trans[i], ndata, fit_start, fit_end, noise, sig,
2484                                        ftype, param[i], paramfree, nparam,
2485////                                    exp_conv, yfit[i], dy[i], &chisq_trans[i]);
2486                                        exp_conv, yfit[0], dy[0], &chisq_trans[i]);
2487
2488                if (ret != 0) {
2489                        if (drop_bad_transients) {
2490                                dbgprintf(3, "In compute_global_exps_fn_final, "
2491                                                  "compute_exps_fn_final returned %d "
2492                                                  "for transient %d\n",
2493                                                  ret, i);
2494                                chisq_trans[i] = -1;
2495                                continue;
2496                        } else {
2497                                dbgprintf(1, "In compute_global_exps_fn_final, "
2498                                                  "compute_exps_fn_final returned %d "
2499                                                  "for transient %d\n",
2500                                                  ret, i);
2501                                return -2;
2502                        }
2503                }
2504
2505                *chisq_global += chisq_trans[i];
2506        }
2507
2508        return 0;
2509}
2510
2511
2512/* This function solves the equation Ax=b, where A is the alpha
2513   matrix, which has the form:
2514
2515     A = (P Q)
2516         (R S)
2517
2518   Here P is an mfit_global x mfit_global square matrix, S is a block
2519   diagonal matrix with ntrans blocks, each of size mfit_local x
2520   mfit_local, and Q and R are the right sizes to make the whole
2521   matrix square.  We solve it by inverting the matrix A using the
2522   formulae given in Numerical Recipes section 2.7, then multiplying
2523   the inverse by b to get x.  We are not too concerned about
2524   accuracy, as this does not affect the solution found, only the
2525   route taken to find it.
2526
2527   Numerical Recipes, section 2.7, notes that A^{-1} is given by:
2528
2529      (P' Q')
2530      (R' S')
2531
2532   with:
2533
2534      P' = (P - Q.S^{-1}.R)^{-1}
2535      Q' = - (P - Q.S^{-1}.R)^{-1} . (Q.S^{-1})
2536      R' = - (S^{-1}.R) . (P - Q.S^{-1}.R)^{-1}
2537      S' = S^{-1} + (S^{-1}.R) . (P - Q.S^{-1}.R)^{-1} . (Q.S^{-1})
2538
2539   We also make use of the fact that A is symmetric, so in particular,
2540   (S^{-1}.R) = (Q.S^{-1})^T and R' = Q'^T.
2541
2542   We are given A as a global_matrix and b as a global_vector.  This
2543   function destroys the original matrix and returns the solution in
2544   place of b.
2545*/
2546
2547int GCI_marquardt_global_solve_eqn(global_matrix A, global_vector b,
2548                        int mfit_global, int mfit_local, int ntrans)
2549{
2550        int row, col, block, i, j;
2551        float x_temp[MAXFIT], x_temp2[MAXFIT];
2552        static float **QS;
2553        static global_vector x;
2554        static int saved_global=0, saved_local=0, saved_ntrans=0;
2555
2556        /* If no local parameters, just do a straight matrix solution */
2557        if (mfit_local == 0) {
2558                if (GCI_solve(A.P, mfit_global, b.global) != 0)
2559                        return -2;
2560                return 0;
2561        }
2562
2563        /* Allocate arrays if necessary */
2564        if ((saved_global != mfit_global) || (saved_local != mfit_local) ||
2565                (saved_ntrans != ntrans)) {
2566                if (saved_global > 0) {
2567                        GCI_ecf_free_matrix(QS);
2568                        GCI_free_global_vector(&x);
2569                        saved_global = 0;
2570                }
2571                if ((QS = GCI_ecf_matrix(mfit_global, mfit_local*ntrans)) == NULL)
2572                        return -1;
2573                if (GCI_alloc_global_vector(&x, mfit_global, mfit_local, ntrans)
2574                        != 0) {
2575                        GCI_ecf_free_matrix(QS);
2576                        return -1;
2577                }
2578                saved_global = mfit_global;
2579                saved_local = mfit_local;
2580                saved_ntrans = ntrans;
2581        }
2582
2583        /* Start by inverting S */
2584        for (block=0; block<ntrans; block++)
2585                if (GCI_invert(A.S[block], mfit_local) != 0)
2586                        return -2;
2587
2588        /* Calculate Q.S^{-1} */
2589        for (row=0; row<mfit_global; row++)
2590                for (block=0; block<ntrans; block++)
2591                        for (i=0; i<mfit_local; i++) {
2592                                QS[row][block*mfit_local + i] = 0;
2593                                for (j=0; j<mfit_local; j++)
2594                                        QS[row][block*mfit_local + i] +=
2595                                                A.Q[row][block*mfit_local + j] * A.S[block][j][i];
2596                        }
2597
2598        /* Now find P - Q.S^{-1}.R */
2599        for (row=0; row<mfit_global; row++)
2600                for (col=0; col<mfit_global; col++)
2601                        for (i=0; i<ntrans*mfit_local; i++)
2602                                A.P[row][col] -= QS[row][i] * A.Q[col][i];  /* Q = R^T */
2603
2604        /* And invert it to get P' */
2605        if (GCI_invert(A.P, mfit_global) != 0)
2606                return -3;
2607
2608        /* Now overwrite Q with Q' */
2609        for (row=0; row<mfit_global; row++)
2610                for (col=0; col<ntrans*mfit_local; col++) {
2611                        A.Q[row][col] = 0;
2612                        for (i=0; i<mfit_global; i++)
2613                                A.Q[row][col] -= A.P[row][i] * QS[i][col];  /* P contains P' */
2614                }
2615
2616        /* Finally, we can solve to find x */
2617        /* We have x.global = P'.(b.global) + Q'.(b.local)
2618           and     x.local  = R'.(b.global) + S'.(b.local)
2619           We do x.global first. */
2620        for (row=0; row<mfit_global; row++) {
2621                x.global[row] = 0;
2622                for (i=0; i<mfit_global; i++)
2623                        x.global[row] += A.P[row][i] * b.global[i];
2624                /* Recall that Q now contains Q' */
2625                for (i=0; i < ntrans * mfit_local; i++)
2626                        x.global[row] += A.Q[row][i] * b.local[i];
2627        }
2628
2629        /* Now x_local; the R'.b_global component first, recalling that R'
2630           is Q' transposed and that Q' is stored in Q. */
2631        for (row=0; row < ntrans * mfit_local; row++) {
2632                x.local[row] = 0;
2633                for (i=0; i<mfit_global; i++)
2634                        x.local[row] += A.Q[i][row] * b.global[i];
2635        }
2636
2637        /* Now S' = S^{-1} + (S^{-1}.R).(P-Q.S^{-1}.R)^{-1}.(Q.S^{-1}).
2638           We first handle the S^{-1} term, then the remaining term.  */
2639        for (block=0; block<ntrans; block++)
2640                for (row=0; row<mfit_local; row++) {
2641                        for (j=0; j<mfit_local; j++)
2642                                x.local[block*mfit_local + row] +=
2643                                        A.S[block][row][j] * b.local[block*mfit_local + j];
2644                }
2645
2646        /* For the remaining term, we have an x.local[row] contribution of
2647
2648              sum_j sum_k sum_l
2649                  (S^{-1}.R)_{row,j} . (P-Q.S^{-1}.R)^{-1}_{j,k} .
2650                  (Q.S^{-1})_{k,l} . b.local_{l}
2651
2652           In order to save computations, we calculate the matrices once
2653           only.  We start with (Q.S^{-1}) . b_local, which is an
2654           mfit_global x 1 array, and store this in x_temp; premultiplying
2655           this by the square matrix P' again gives an mfit_global x 1
2656           array, which goes in x_temp2, then premultiplying this by
2657           (S^{-1}.R) gives an (ntrans * mfit_local) x 1 array which is
2658           added directly onto x_local.  Recall also that S^{-1}.R is the
2659           transpose of Q.S^{-1}, which is currently stored in QS, and
2660           that the middle term is currently stored in P. */
2661
2662        for (row=0; row<mfit_global; row++) {
2663                x_temp[row] = 0;
2664                for (i=0; i < ntrans*mfit_local; i++)
2665                        x_temp[row] += QS[row][i] * b.local[i];
2666        }
2667        for (row=0; row<mfit_global; row++) {
2668                x_temp2[row] = 0;
2669                for (i=0; i<mfit_global; i++)
2670                        x_temp2[row] += A.P[row][i] * x_temp[i];
2671        }
2672        /* Again, S^{-1}.R is the transpose of Q.S^{-1} */
2673        for (row=0; row < ntrans * mfit_local; row++) {
2674                for (i=0; i<mfit_global; i++)
2675                        x.local[row] += QS[i][row] * x_temp2[i];
2676        }
2677
2678        /* And we're done, once we've copied x into b */
2679        GCI_copy_global_vector(b, x, mfit_global, mfit_local, ntrans);
2680
2681        return 0;
2682}
2683
2684
2685/*         *****        GENERIC FUNCTION CODE        *****        */
2686
2687/* These functions are essentially the same as the above functions
2688   GCI_marquardt_global_exps_instr and
2689   GCI_marquardt_global_exps_do_fit_instr and the latter's dependents,
2690   except that this version takes an arbitrary function and a list of
2691   global parameters.  This function is designed to be called from
2692   external code.  It is nowhere near as efficient as the streamlined
2693   code above for globally fitting taus for multi-exponential models.
2694   Also, this function must be provided with meaningful starting
2695   estimates for all parameters.
2696   */
2697
2698int GCI_marquardt_global_generic_instr(float xincr, float **trans,
2699                                        int ndata, int ntrans, int fit_start, int fit_end,
2700                                        float instr[], int ninstr,
2701                                        noise_type noise, float sig[],
2702                                        float **param, int paramfree[], int nparam, int gparam[],
2703                                        restrain_type restrain, float chisq_delta,
2704                                        void (*fitfunc)(float, float [], float *, float [], int),
2705                                        float **fitted, float **residuals,
2706                                        float chisq_trans[], float *chisq_global, int *df)
2707{
2708        float **covar, **alpha, *scaled_instr, instrsum;
2709        int i, ret;
2710        int mlocal, mglobal;
2711
2712        /* Some basic parameter checks */
2713        if (xincr <= 0) return -1;
2714        if (ntrans < 1) return -1;
2715        if (ndata < 1)  return -1;
2716        if (fit_start < 0 || fit_end > ndata) return -1;
2717        if (ninstr < 1) return -1;
2718        if (nparam < 1) return -1;
2719
2720        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL)
2721                return -2;
2722
2723        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
2724                GCI_ecf_free_matrix(covar);
2725                return -3;
2726        }
2727
2728        if ((scaled_instr = (float *) malloc(ninstr * sizeof(float))) == NULL) {
2729                GCI_ecf_free_matrix(covar);
2730                GCI_ecf_free_matrix(alpha);
2731                return -4;
2732        }
2733
2734        /* Scale the instrument response */
2735        for (i=0, instrsum=0; i<ninstr; i++)
2736                instrsum += instr[i];
2737        if (instrsum == 0) {
2738                GCI_ecf_free_matrix(covar);
2739                GCI_ecf_free_matrix(alpha);
2740                free(scaled_instr);
2741                return -6;
2742        }
2743       
2744        for (i=0; i<ninstr; i++)
2745                scaled_instr[i] = instr[i] / instrsum;
2746
2747        /* Now call the global fitting function. */
2748        ret = GCI_marquardt_global_generic_do_fit_instr(
2749                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2750                        scaled_instr, ninstr, noise, sig,
2751                        param, paramfree, nparam, gparam, restrain, chisq_delta,
2752                        fitfunc, fitted, residuals, covar, alpha,
2753                        chisq_trans, chisq_global);
2754
2755        GCI_ecf_free_matrix(covar);
2756        GCI_ecf_free_matrix(alpha);
2757        free(scaled_instr);
2758        GCI_marquardt_cleanup();
2759
2760        if (ret < 0) {
2761                dbgprintf(1, "Fit failed, ret = %d\n", ret);
2762                GCI_ecf_free_matrix(covar);
2763                GCI_ecf_free_matrix(alpha);
2764                free(scaled_instr);
2765                return -10 + ret;
2766        }
2767
2768        dbgprintf(2, "Fit succeeded, ret = %d\n", ret);
2769
2770        /* Before we return, calculate the number of degrees of freedom */
2771        /* The number of degrees of freedom is given by:
2772              d.f. = ntrans * ((fit_end - fit_start) - # free local parameters)
2773                         - # free global parameters
2774        */
2775
2776        mglobal = mlocal = 0;
2777        for (i=0; i<nparam; i++)
2778                if (paramfree[i]) {
2779                        if (gparam[i]) mglobal++;
2780                        else mlocal++;
2781                }
2782
2783        *df = ntrans * ((fit_end - fit_start) - mlocal) - mglobal;
2784
2785        GCI_ecf_free_matrix(covar);
2786        GCI_ecf_free_matrix(alpha);
2787        free(scaled_instr);
2788
2789        return ret;
2790}
2791
2792
2793int GCI_marquardt_global_generic_do_fit_instr(
2794                float xincr, float **trans, int ndata, int ntrans,
2795                int fit_start, int fit_end, float instr[], int ninstr,
2796                noise_type noise, float sig[],
2797                float **param, int paramfree[], int nparam, int gparam[],
2798                restrain_type restrain, float chisq_delta,
2799                void (*fitfunc)(float, float [], float *, float [], int),
2800                float **fitted, float **residuals,
2801                float **covar_scratch, float **alpha_scratch,
2802                float *chisq_trans, float *chisq_global)
2803{
2804        // PRB 03/07 Although **fitted and **residuals are provided only one "transient" is required and used, fitted[0] and residuals[0]
2805
2806        float alambda, ochisq_global, *ochisq_trans;
2807        int i, k, itst, itst_max;
2808        int ret;
2809
2810        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
2811
2812        /* If there are no global parameters being fitted, we simply fit
2813           each local set. */
2814        for (i=0; i<nparam; i++) {
2815                if (gparam[i] && paramfree[i]) {
2816                        i = -1; /* sentinel value */
2817                        break;
2818                }
2819        }
2820
2821        if (i >= 0) { /* no globals to fit */
2822                *chisq_global = 0;
2823
2824                for (i=0; i<ntrans; i++) {
2825                        ret = GCI_marquardt_instr(xincr, trans[i],
2826                                        ndata, fit_start, fit_end, instr, ninstr, noise, sig,
2827                                        param[i], paramfree, nparam, restrain, fitfunc,
2828                                        fitted[0], residuals[0], covar_scratch, alpha_scratch,
2829                                        &chisq_trans[i], chisq_delta, 0, NULL);
2830                        if (ret < 0) {
2831                                dbgprintf(1, "In do_fit_instr, marquardt_instr returned %d "
2832                                                  "for transient %d\n", ret, i);
2833                                return -10 + ret;
2834                        } else {
2835                                *chisq_global += chisq_trans[i];
2836                        }
2837                }
2838                return 0;
2839        }
2840
2841        /* If there are no free local variables to fit, we still do the
2842           global fitting, but we have to be a little careful in some of
2843           the later routines */
2844
2845        /* Now allocate all of the arrays we will need. */
2846
2847        if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL)
2848                return -1;
2849
2850        /* We now begin our standard Marquardt loop, with several
2851           modifications */
2852        alambda = -1;
2853        ret = GCI_marquardt_global_generic_global_step(
2854                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2855                                instr, ninstr, noise, sig,
2856                                param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2857                                fitted, residuals, chisq_trans, chisq_global,
2858                                alpha_scratch, &alambda);
2859        if (ret != 0) {
2860                dbgprintf(1, "In do_fit_instr, first global_step returned %d\n", ret);
2861                if (ret != -1) {
2862                        /* Wasn't a memory error, so unallocate arrays */
2863                        alambda = 0.0;
2864                        GCI_marquardt_global_generic_global_step(
2865                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2866                                instr, ninstr, noise, sig,
2867                                param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2868                                fitted, residuals, chisq_trans, chisq_global,
2869                                alpha_scratch, &alambda);
2870                }
2871                free(ochisq_trans);
2872                return ret;
2873        }
2874
2875        k = 1;  /* Iteration counter */
2876        itst = 0;
2877        for (;;) {
2878                dbgprintf(3, "In do_fit_instr, beginning iteration %d:\n", k);
2879                dbgprintf(3, " itst = %d, chisq_global = %.4f\n", itst, *chisq_global);
2880
2881                k++;
2882                if (k > MAXITERS) {
2883                        free(ochisq_trans);
2884                        return -2;
2885                }
2886
2887                ochisq_global = *chisq_global;
2888                for (i=0; i<ntrans; i++)
2889                        ochisq_trans[i] = chisq_trans[i];
2890
2891                ret = GCI_marquardt_global_generic_global_step(
2892                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2893                                        instr, ninstr, noise, sig,
2894                                        param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2895                                        fitted, residuals, chisq_trans, chisq_global,
2896                                        alpha_scratch, &alambda);
2897                if (ret != 0) {
2898                        dbgprintf(1, "In do_fit_instr, second global_step returned %d\n",
2899                                          ret);
2900                        /* Unallocate arrays */
2901                        alambda = 0.0;
2902                        GCI_marquardt_global_generic_global_step(
2903                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2904                                instr, ninstr, noise, sig,
2905                                param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2906                                fitted, residuals, chisq_trans, chisq_global,
2907                                alpha_scratch, &alambda);
2908                        free(ochisq_trans);
2909                        return ret;
2910                }
2911
2912                if (*chisq_global > ochisq_global)
2913                        itst = 0;
2914                else {
2915                        /* Let's try this approach; I really don't know what will
2916                           be best */
2917                        float maxdiff;
2918
2919                        maxdiff = 0.0;
2920                        for (i=0; i<ntrans; i++)
2921                                if (ochisq_trans[i] - chisq_trans[i] > maxdiff)
2922                                        maxdiff = ochisq_trans[i] - chisq_trans[i];
2923
2924                        if (maxdiff < chisq_delta)
2925                                itst++;
2926                        dbgprintf(3, "In do_fit_instr, maxdiff = %.3f:\n", maxdiff);
2927                }
2928
2929                if (itst < itst_max) continue;
2930
2931                /* Endgame */
2932                alambda = 0.0;
2933                ret = GCI_marquardt_global_generic_global_step(
2934                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2935                                        instr, ninstr, noise, sig,
2936                                        param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2937                                        fitted, residuals, chisq_trans, chisq_global,
2938                                        alpha_scratch, &alambda);
2939                if (ret != 0) {
2940                        dbgprintf(1, "In do_fit_instr, final global_step returned %d\n",
2941                                          ret);
2942                        free(ochisq_trans);
2943                        return ret;
2944                }
2945
2946                free(ochisq_trans);
2947                return k;  /* We're done now */
2948        }
2949}
2950
2951
2952/* And this one is basically a specialised GCI_marquardt_instr_step
2953   for the global fitting setup. */
2954
2955#define do_frees \
2956        if (fnvals) free(fnvals);\
2957        if (dy_dparam_pure) GCI_ecf_free_matrix(dy_dparam_pure);\
2958        if (dy_dparam_conv) GCI_ecf_free_matrix(dy_dparam_conv);
2959
2960int GCI_marquardt_global_generic_global_step(
2961                                float xincr, float **trans,
2962                                int ndata, int ntrans, int fit_start, int fit_end,
2963                                float instr[], int ninstr,
2964                                noise_type noise, float sig[],
2965                                float **param, int paramfree[], int nparam, int gparam[],
2966                                restrain_type restrain, float chisq_delta,
2967                                void (*fitfunc)(float, float [], float *, float [], int),
2968                                float **yfit, float **dy,
2969                                float *chisq_trans, float *chisq_global,
2970                                float **alpha_scratch, float *alambda)
2971{
2972        int i, j, ret;
2973        static global_matrix alpha, covar;
2974        static global_vector beta, dparam;
2975        static float **paramtry;
2976        static int mfit_local, mfit_global;
2977        static int gindex[MAXFIT], lindex[MAXFIT];
2978        static float ochisq_global, *ochisq_trans;
2979        static int initialised=0;
2980
2981        // The following are declared here to retain some optimisation by not repeatedly mallocing
2982        // (only once per transient), but to remain thread safe.
2983        // They are malloced by lower fns but at the end, freed by this fn.
2984        // These vars were global or static before thread safety was introduced.
2985        float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL;
2986        int fnvals_len=0, dy_dparam_nparam_size=0;
2987
2988        if (nparam > MAXFIT)
2989                return -10;
2990        if (xincr <= 0)
2991                return -11;
2992        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
2993                return -12;
2994
2995        /* Initialisation */
2996        /* We assume we're given sensible starting values for param[] */
2997        if (*alambda < 0.0) {
2998                /* Start by allocating lots of variables we will need */
2999                mfit_local = mfit_global = 0;
3000
3001                for (i=0; i<nparam; i++) {
3002                        if (paramfree[i]) {
3003                                if (gparam[i])
3004                                        gindex[mfit_global++] = i;
3005                                else
3006                                        lindex[mfit_local++] = i;
3007                        }
3008                }
3009
3010                if (initialised) {
3011                        GCI_free_global_matrix(&alpha);   GCI_free_global_matrix(&covar);
3012                        GCI_ecf_free_matrix(paramtry);        GCI_free_global_vector(&beta);
3013                        GCI_free_global_vector(&dparam);  free(ochisq_trans);
3014                        initialised = 0;
3015                }
3016
3017                if (GCI_alloc_global_matrix(&alpha, mfit_global, mfit_local, ntrans)
3018                        != 0)
3019                        return -1;
3020
3021                if (GCI_alloc_global_matrix(&covar, mfit_global, mfit_local, ntrans)
3022                        != 0) {
3023                        GCI_free_global_matrix(&alpha);
3024                        return -1;
3025                }
3026
3027                if ((paramtry = GCI_ecf_matrix(ntrans, nparam)) == NULL) {
3028                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3029                        return -1;
3030                }
3031
3032                if (GCI_alloc_global_vector(&beta, mfit_global, mfit_local, ntrans)
3033                        != 0) {
3034                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3035                        GCI_ecf_free_matrix(paramtry);
3036                        return -1;
3037                }
3038
3039                if (GCI_alloc_global_vector(&dparam, mfit_global, mfit_local, ntrans)
3040                        != 0) {
3041                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3042                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
3043                        return -1;
3044                }
3045
3046                if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float)))
3047                        == NULL) {
3048                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3049                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
3050                        GCI_free_global_vector(&dparam);
3051                        return -1;
3052                }
3053
3054                initialised = 1;
3055
3056                if (GCI_marquardt_global_compute_global_generic_fn(
3057                                        xincr, trans, ndata, ntrans,
3058                                        fit_start, fit_end, instr, ninstr, noise, sig,
3059                                        param, paramfree, nparam, gparam,
3060                                        mfit_global, mfit_local, gindex, lindex, fitfunc,
3061                                        yfit, dy, alpha, beta, alpha_scratch,
3062                                        chisq_trans, chisq_global, *alambda,
3063                                        &fnvals, &dy_dparam_pure, &dy_dparam_conv,
3064                                        &fnvals_len, &dy_dparam_nparam_size) != 0)
3065                        return -2;
3066
3067                *alambda = 0.001;
3068                ochisq_global = *chisq_global;
3069                for (i=0; i<ntrans; i++)
3070                        ochisq_trans[i] = chisq_trans[i];
3071
3072                /* Initialise paramtry to param */
3073                for (i=0; i<ntrans; i++) {
3074                        for (j=0; j<nparam; j++)
3075                                paramtry[i][j] = param[i][j];
3076                }
3077        }
3078
3079        /* Once converged, evaluate covariance matrix */
3080        if (*alambda == 0) {
3081                if (GCI_marquardt_global_compute_global_generic_fn_final(
3082                                        xincr, trans, ndata, ntrans,
3083                                        fit_start, fit_end, instr, ninstr, noise, sig,
3084                                        param, paramfree, nparam, gparam,
3085                                        mfit_global, mfit_local, gindex, lindex, fitfunc,
3086                                        yfit, dy, chisq_trans, chisq_global,
3087                                        &fnvals, &dy_dparam_pure, &dy_dparam_conv,
3088                                        &fnvals_len, &dy_dparam_nparam_size) != 0)
3089                        return -3;
3090                /* Don't need to do this here; if we wished to, we'd have to
3091                   move this code (the "if (*alambda == 0)" block) to after
3092                   the Gauss-Jordan call.  We'd also need to rewrite it for
3093                   our situation.... */
3094                //      if (mfit < nparam) {  /* no need to do this otherwise */
3095                //              GCI_covar_sort(covar, nparam, paramfree, mfit);
3096                //              GCI_covar_sort(alpha, nparam, paramfree, mfit);
3097                //      }
3098                GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3099                GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
3100                GCI_free_global_vector(&dparam); free(ochisq_trans);
3101                initialised = 0;
3102                return 0;
3103        }
3104
3105        /* Alter linearised fitting matrix by augmenting diagonal
3106           elements. */
3107        GCI_copy_global_matrix(covar, alpha, mfit_global, mfit_local, ntrans);
3108        GCI_copy_global_vector(dparam, beta, mfit_global, mfit_local, ntrans);
3109        for (j=0; j<mfit_global; j++)
3110                covar.P[j][j] *= 1.0 + (*alambda);
3111        for (i=0; i<ntrans; i++)
3112                for (j=0; j<mfit_local; j++)
3113                        covar.S[i][j][j] *= 1.0 + (*alambda);
3114
3115        /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */
3116        if (GCI_marquardt_global_solve_eqn(covar, dparam,
3117                                                                           mfit_global, mfit_local, ntrans) != 0)
3118                return -3;
3119
3120        /* Did the trial succeed?  Modify param by dparam... */
3121        for (i=0; i<ntrans; i++) {
3122                for (j=0; j<mfit_global; j++)
3123                        paramtry[i][gindex[j]] = param[i][gindex[j]] + dparam.global[j];
3124                for (j=0; j<mfit_local; j++)
3125                        paramtry[i][lindex[j]] =
3126                                param[i][lindex[j]] + dparam.local[i*mfit_local + j];
3127        }
3128
3129        for (i=0; i<ntrans; i++) {
3130                if (restrain == ECF_RESTRAIN_DEFAULT)
3131                        ret = check_ecf_params (paramtry[i], nparam, fitfunc);
3132                else
3133                        ret = check_ecf_user_params (paramtry[i], nparam, fitfunc);
3134
3135                if (ret != 0) {
3136                        /* Bad parameters, increase alambda and return */
3137                        *alambda *= 10.0;
3138                        return 0;
3139                }
3140        }
3141
3142        if (GCI_marquardt_global_compute_global_generic_fn(
3143                                xincr, trans, ndata, ntrans,
3144                                fit_start, fit_end, instr, ninstr, noise, sig,
3145                                paramtry, paramfree, nparam, gparam,
3146                                mfit_global, mfit_local, gindex, lindex, fitfunc,
3147                                yfit, dy, covar, dparam, alpha_scratch,
3148                                chisq_trans, chisq_global, *alambda,
3149                                &fnvals, &dy_dparam_pure, &dy_dparam_conv,
3150                                &fnvals_len, &dy_dparam_nparam_size) != 0)
3151                return -2;
3152
3153        /* Success, accept the new solution */
3154        if (*chisq_global < ochisq_global) {
3155                *alambda *= 0.1;
3156                ochisq_global = *chisq_global;
3157                for (i=0; i<ntrans; i++)
3158                        ochisq_trans[i] = chisq_trans[i];
3159                GCI_copy_global_matrix(alpha, covar, mfit_global, mfit_local, ntrans);
3160                GCI_copy_global_vector(beta, dparam, mfit_global, mfit_local, ntrans);
3161                for (i=0; i<ntrans; i++) {
3162                        for (j=0; j<nparam; j++)
3163                                param[i][j] = paramtry[i][j];
3164                }
3165        } else { /* Failure, increase alambda and return */
3166                *alambda *= 10.0;
3167                *chisq_global = ochisq_global;
3168                for (i=0; i<ntrans; i++)
3169                        chisq_trans[i] = ochisq_trans[i];
3170        }
3171
3172        return 0;
3173}
3174
3175
3176/* Here we use alpha only for scratch space */
3177int GCI_marquardt_global_compute_global_generic_fn(
3178                float xincr, float **trans, int ndata, int ntrans,
3179                int fit_start, int fit_end, float instr[], int ninstr,
3180                noise_type noise, float sig[],
3181                float **param, int paramfree[], int nparam, int gparam[],
3182                int mfit_global, int mfit_local, int gindex[], int lindex[],
3183                void (*fitfunc)(float, float [], float *, float [], int),
3184                float **yfit, float **dy, global_matrix alpha, global_vector beta,
3185                float **alpha_scratch, float *chisq_trans, float *chisq_global,
3186                float alambda,
3187                float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
3188                int *pfnvals_len, int *pdy_dparam_nparam_size)
3189{
3190        int i, j, k, ret;
3191        float beta_scratch[MAXFIT];  /* scratch space */
3192
3193        /* We initialise P and beta_global to zero; the others don't
3194           matter, as they will be totally overwritten */
3195        for (i=0; i<mfit_global; i++) {
3196                for (j=0; j<mfit_global; j++)
3197                        alpha.P[i][j] = 0;
3198                beta.global[i] = 0;
3199        }
3200        *chisq_global = 0.0;
3201
3202        for (i=0; i<ntrans; i++) {
3203                /* Only pass the true alambda, used for initialisation, for
3204                   the first transient */
3205                ret = GCI_marquardt_compute_fn_instr(
3206                                        xincr, trans[i], ndata, fit_start, fit_end,
3207                                        instr, ninstr, noise, sig,
3208                                        param[i], paramfree, nparam, fitfunc,
3209////                                    yfit[i], dy[i], alpha_scratch, beta_scratch,
3210                                        yfit[0], dy[0], alpha_scratch, beta_scratch,
3211                                        &chisq_trans[i], 0.0f, (i == 0) ? alambda : 0.0, //TODO ARG added 0.0f here for new old_chisq parameter
3212                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
3213                                        pfnvals_len, pdy_dparam_nparam_size);
3214
3215                if (ret != 0) {
3216                        dbgprintf(1, "In compute_global_generic_fn, "
3217                                          "compute_fn_instr returned %d for transient %d\n",
3218                                          ret, i);
3219                        return -2;
3220                }
3221
3222                /* So now have to populate alpha and beta with the contents of
3223                   alpha_scratch and beta_scratch. */
3224
3225                for (j=0; j<mfit_global; j++) {
3226                        for (k=0; k<mfit_global; k++)
3227                                alpha.P[j][k] += alpha_scratch[gindex[j]][gindex[k]];
3228                        for (k=0; k<mfit_local; k++)
3229                                alpha.Q[j][i*mfit_local + k] =
3230                                        alpha_scratch[gindex[j]][lindex[k]];
3231                        beta.global[j] += beta_scratch[gindex[j]];
3232                }
3233                for (j=0; j<mfit_local; j++) {
3234                        for (k=0; k<mfit_local; k++)
3235                                alpha.S[i][j][k] = alpha_scratch[lindex[j]][lindex[k]];
3236                        beta.local[i*mfit_local + j] = beta_scratch[lindex[j]];
3237                }
3238
3239                *chisq_global += chisq_trans[i];
3240        }
3241
3242        return 0;
3243}
3244
3245
3246/* And the final variant */
3247int GCI_marquardt_global_compute_global_generic_fn_final(
3248                float xincr, float **trans, int ndata, int ntrans,
3249                int fit_start, int fit_end, float instr[], int ninstr,
3250                noise_type noise, float sig[],
3251                float **param, int paramfree[], int nparam, int gparam[],
3252                int mfit_global, int mfit_local, int gindex[], int lindex[],
3253                void (*fitfunc)(float, float [], float *, float [], int),
3254                float **yfit, float **dy,
3255                float *chisq_trans, float *chisq_global,
3256                float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
3257                int *pfnvals_len, int *pdy_dparam_nparam_size)
3258{
3259        int i, ret;
3260
3261        *chisq_global = 0.0;
3262
3263        for (i=0; i<ntrans; i++) {
3264                /* Only pass the true alambda, used for initialisation, for
3265                   the first transient */
3266                ret = GCI_marquardt_compute_fn_final_instr(
3267                                        xincr, trans[i], ndata, fit_start, fit_end,
3268                                        instr, ninstr, noise, sig,
3269                                        param[i], paramfree, nparam, fitfunc,
3270//                                      yfit[i], dy[i], &chisq_trans[i]);
3271                                        yfit[0], dy[0], &chisq_trans[i],
3272                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
3273                                        pfnvals_len, pdy_dparam_nparam_size);
3274
3275                if (ret != 0) {
3276                        dbgprintf(1, "In compute_global_generic_fn_final, "
3277                                          "compute_fn_final_instr returned %d for transient %d\n",
3278                                          ret, i);
3279                        return -2;
3280                }
3281
3282                *chisq_global += chisq_trans[i];
3283        }
3284
3285        return 0;
3286}
3287
3288
3289// Emacs settings:
3290// Local variables:
3291// mode: c
3292// c-basic-offset: 4
3293// tab-width: 4
3294// End:
Note: See TracBrowser for help on using the repository browser.