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

Revision 7720, 99.5 KB checked in by aivar, 9 years ago (diff)

Updated GCI_marquardt_global_compute_exps_fn.

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