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

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

Added GPL license to all source code.

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