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

Revision 7611, 99.7 KB checked in by aivar, 9 years ago (diff)

Added float chisq_delta parameter to GCI_marquart_global_generic_global_step.

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);
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)
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_gauss_jordan solves Ax=b rather than AX=B */
1184        if (GCI_gauss_jordan(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) != 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, float *chisq)
1253{
1254        int i, j, k, l, m, mfit;
1255        float wt, sig2i, y_ymod, dy_dparam[MAXFIT];
1256
1257        for (j=0, mfit=0; j<nparam; j++)
1258                if (paramfree[j])
1259                        mfit++;
1260
1261        /* Initialise (symmetric) alpha, beta */
1262        for (j=0; j<mfit; j++) {
1263                for (k=0; k<=j; k++)
1264                        alpha[j][k] = 0.0;
1265                beta[j] = 0.0;
1266        }
1267
1268        /* Calculation of the fitting data will depend upon the type of
1269           noise.  Since there's no convolution involved here, this is
1270           very easy. */
1271
1272        switch (ftype) {
1273        case FIT_GLOBAL_MULTIEXP:
1274                switch (noise) {
1275                case NOISE_CONST:
1276                        *chisq = 0.0;
1277                        /* Summation loop over all data */
1278                        for (i=fit_start; i<fit_end; i++) {
1279                                yfit[i] = param[0];  /* Z */
1280                                dy_dparam[0] = 1.0;
1281
1282                                for (j=1; j<nparam; j+=2) {
1283                                        yfit[i] += param[j] * exp_conv[j+1][i];
1284                                        /* A_j . exp(-x/tau_j) */
1285                                        dy_dparam[j] = exp_conv[j+1][i];
1286                                        /* exp(-x/tau_j) */
1287                                        dy_dparam[j+1] = param[j] * exp_conv[j][i];
1288                                        /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1289                                }
1290                                dy[i] = y[i] - yfit[i];
1291
1292                                for (j=0, l=0; l<nparam; l++) {
1293                                        if (paramfree[l]) {
1294                                                wt = dy_dparam[l];
1295                                                for (k=0, m=0; m<=l; m++)
1296                                                        if (paramfree[m])
1297                                                                alpha[j][k++] += wt * dy_dparam[m];
1298                                                beta[j] += wt * dy[i];
1299                                                j++;
1300                                        }
1301                                }
1302                                /* And find chi^2 */
1303                                *chisq += dy[i] * dy[i];
1304                        }
1305
1306                        /* Now divide everything by sigma^2 */
1307                        sig2i = 1.0 / (sig[0] * sig[0]);
1308                        *chisq *= sig2i;
1309                        for (j=0; j<mfit; j++) {
1310                                for (k=0; k<=j; k++)
1311                                        alpha[j][k] *= sig2i;
1312                                beta[j] *= sig2i;
1313                        }
1314                        break;
1315
1316                case NOISE_GIVEN:  /* This is essentially the NR version */
1317                        *chisq = 0.0;
1318                        /* Summation loop over all data */
1319                        for (i=fit_start; i<fit_end; i++) {
1320                                yfit[i] = param[0];  /* Z */
1321                                dy_dparam[0] = 1.0;
1322
1323                                for (j=1; j<nparam; j+=2) {
1324                                        yfit[i] += param[j] * exp_conv[j+1][i];
1325                                        /* A_j . exp(-x/tau_j) */
1326                                        dy_dparam[j] = exp_conv[j+1][i];
1327                                        /* exp(-x/tau_j) */
1328                                        dy_dparam[j+1] = param[j] * exp_conv[j][i];
1329                                            /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1330                                }
1331                                sig2i = 1.0 / (sig[i] * sig[i]);
1332                                dy[i] = y[i] - yfit[i];
1333
1334                                for (j=0, l=0; l<nparam; l++) {
1335                                        if (paramfree[l]) {
1336                                                wt = dy_dparam[l] * sig2i;
1337                                                for (k=0, m=0; m<=l; m++)
1338                                                        if (paramfree[m])
1339                                                                alpha[j][k++] += wt * dy_dparam[m];
1340                                                beta[j] += wt * dy[i];
1341                                                j++;
1342                                        }
1343                                }
1344                                /* And find chi^2 */
1345                                *chisq += dy[i] * dy[i] * sig2i;
1346                        }
1347                        break;
1348
1349                case NOISE_POISSON_DATA:
1350                        *chisq = 0.0;
1351                        /* Summation loop over all data */
1352                        for (i=fit_start; i<fit_end; i++) {
1353                                yfit[i] = param[0];  /* Z */
1354                                dy_dparam[0] = 1.0;
1355
1356                                for (j=1; j<nparam; j+=2) {
1357                                        yfit[i] += param[j] * exp_conv[j+1][i];
1358                                            /* A_j . exp(-x/tau_j) */
1359                                        dy_dparam[j] = exp_conv[j+1][i];
1360                                            /* exp(-x/tau_j) */
1361                                        dy_dparam[j+1] = param[j] * exp_conv[j][i];
1362                                            /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1363                                }
1364                                sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15);
1365                                dy[i] = y[i] - yfit[i];
1366
1367                                for (j=0, l=0; l<nparam; l++) {
1368                                        if (paramfree[l]) {
1369                                                wt = dy_dparam[l] * sig2i;
1370                                                for (k=0, m=0; m<=l; m++)
1371                                                        if (paramfree[m])
1372                                                                alpha[j][k++] += wt * dy_dparam[m];
1373                                                beta[j] += wt * dy[i];
1374                                                j++;
1375                                        }
1376                                }
1377                                /* And find chi^2 */
1378                                *chisq += dy[i] * dy[i] * sig2i;
1379                        }
1380                        break;
1381
1382                case NOISE_POISSON_FIT:
1383                        *chisq = 0.0;
1384                        /* Summation loop over all data */
1385                        for (i=fit_start; i<fit_end; i++) {
1386                                yfit[i] = param[0];  /* Z */
1387                                dy_dparam[0] = 1.0;
1388
1389                                for (j=1; j<nparam; j+=2) {
1390                                        yfit[i] += param[j] * exp_conv[j+1][i];
1391                                            /* A_j . exp(-x/tau_j) */
1392                                        dy_dparam[j] = exp_conv[j+1][i];
1393                                            /* exp(-x/tau_j) */
1394                                        dy_dparam[j+1] = param[j] * exp_conv[j][i];
1395                                            /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1396                                }
1397                                sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15);
1398                                dy[i] = y[i] - yfit[i];
1399
1400                                for (j=0, l=0; l<nparam; l++) {
1401                                        if (paramfree[l]) {
1402                                                wt = dy_dparam[l] * sig2i;
1403                                                for (k=0, m=0; m<=l; m++)
1404                                                        if (paramfree[m])
1405                                                                alpha[j][k++] += wt * dy_dparam[m];
1406                                                beta[j] += wt * dy[i];
1407                                                j++;
1408                                        }
1409                                }
1410                                /* And find chi^2 */
1411                                *chisq += dy[i] * dy[i] * sig2i;
1412                        }
1413                        break;
1414
1415                case NOISE_GAUSSIAN_FIT:
1416                        *chisq = 0.0;
1417                        /* Summation loop over all data */
1418                        for (i=fit_start; i<fit_end; i++) {
1419                                yfit[i] = param[0];  /* Z */
1420                                dy_dparam[0] = 1.0;
1421
1422                                for (j=1; j<nparam; j+=2) {
1423                                        yfit[i] += param[j] * exp_conv[j+1][i];
1424                                            /* A_j . exp(-x/tau_j) */
1425                                        dy_dparam[j] = exp_conv[j+1][i];
1426                                            /* exp(-x/tau_j) */
1427                                        dy_dparam[j+1] = param[j] * exp_conv[j][i];
1428                                            /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1429                                }
1430                                sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1431                                dy[i] = y[i] - yfit[i];
1432
1433                                for (j=0, l=0; l<nparam; l++) {
1434                                        if (paramfree[l]) {
1435                                                wt = dy_dparam[l] * sig2i;
1436                                                for (k=0, m=0; m<=l; m++)
1437                                                        if (paramfree[m])
1438                                                                alpha[j][k++] += wt * dy_dparam[m];
1439                                                beta[j] += wt * dy[i];
1440                                                j++;
1441                                        }
1442                                }
1443                                /* And find chi^2 */
1444                                *chisq += dy[i] * dy[i] * sig2i;
1445                        }
1446                        break;
1447
1448                case NOISE_MLE:
1449                        *chisq = 0.0;
1450                        /* Summation loop over all data */
1451                        for (i=fit_start; i<fit_end; i++) {
1452                                yfit[i] = param[0];  /* Z */
1453                                dy_dparam[0] = 1.0;
1454
1455                                for (j=1; j<nparam; j+=2) {
1456                                        yfit[i] += param[j] * exp_conv[j+1][i];
1457                                            /* A_j . exp(-x/tau_j) */
1458                                        dy_dparam[j] = exp_conv[j+1][i];
1459                                            /* exp(-x/tau_j) */
1460                                        dy_dparam[j+1] = param[j] * exp_conv[j][i];
1461                                            /* A_j * x * exp(-x/tau_j) / tau_j^2 */
1462                                }
1463                                sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1464                                dy[i] = y[i] - yfit[i];
1465                                y_ymod=y[i]/yfit[i];
1466
1467                                for (j=0, l=0; l<nparam; l++) {
1468                                        if (paramfree[l]) {
1469                                                wt = dy_dparam[l] * sig2i;
1470                                                for (k=0, m=0; m<=l; m++)
1471                                                        if (paramfree[m])
1472                                                                alpha[j][k++] += wt*dy_dparam[m]*y_ymod; // wt * dy_dparam[m];
1473                                                beta[j] += wt * dy[i];
1474                                                j++;
1475                                        }
1476                                }
1477                                // And find chi^2
1478                                if (yfit[i]<=0.0)
1479                                        ; // do nothing
1480                                else if (y[i]==0.0)
1481                                        *chisq += 2.0*yfit[i];   // to avoid NaN from log
1482                                else
1483                                        *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1484                        }
1485                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1486                        break;
1487
1488                default:
1489                        return -3;
1490                        /* break; */   // (unreachable code)
1491                }
1492                break;
1493
1494        case FIT_GLOBAL_STRETCHEDEXP:
1495                switch (noise) {
1496                case NOISE_CONST:
1497                        *chisq = 0.0;
1498                        /* Summation loop over all data */
1499                        for (i=fit_start; i<fit_end; i++) {
1500                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1501                                dy[i] = y[i] - yfit[i];
1502
1503                                dy_dparam[0] = 1.0;
1504                                dy_dparam[1] = exp_conv[1][i];
1505                                dy_dparam[2] = param[1] * exp_conv[2][i];
1506                                dy_dparam[3] = param[1] * exp_conv[3][i];
1507
1508                                for (j=0, l=0; l<nparam; l++) {
1509                                        if (paramfree[l]) {
1510                                                wt = dy_dparam[l];
1511                                                for (k=0, m=0; m<=l; m++)
1512                                                        if (paramfree[m])
1513                                                                alpha[j][k++] += wt * dy_dparam[m];
1514                                                beta[j] += wt * dy[i];
1515                                                j++;
1516                                        }
1517                                }
1518                                /* And find chi^2 */
1519                                *chisq += dy[i] * dy[i];
1520                        }
1521
1522                        /* Now divide everything by sigma^2 */
1523                        sig2i = 1.0 / (sig[0] * sig[0]);
1524                        *chisq *= sig2i;
1525                        for (j=0; j<mfit; j++) {
1526                                for (k=0; k<=j; k++)
1527                                        alpha[j][k] *= sig2i;
1528                                beta[j] *= sig2i;
1529                        }
1530                        break;
1531
1532                case NOISE_GIVEN:  /* This is essentially the NR version */
1533                        *chisq = 0.0;
1534                        /* Summation loop over all data */
1535                        for (i=fit_start; i<fit_end; i++) {
1536                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1537                                dy[i] = y[i] - yfit[i];
1538
1539                                dy_dparam[0] = 1.0;
1540                                dy_dparam[1] = exp_conv[1][i];
1541                                dy_dparam[2] = param[1] * exp_conv[2][i];
1542                                dy_dparam[3] = param[1] * exp_conv[3][i];
1543
1544                                sig2i = 1.0 / (sig[i] * sig[i]);
1545
1546                                for (j=0, l=0; l<nparam; l++) {
1547                                        if (paramfree[l]) {
1548                                                wt = dy_dparam[l] * sig2i;
1549                                                for (k=0, m=0; m<=l; m++)
1550                                                        if (paramfree[m])
1551                                                                alpha[j][k++] += wt * dy_dparam[m];
1552                                                beta[j] += wt * dy[i];
1553                                                j++;
1554                                        }
1555                                }
1556                                /* And find chi^2 */
1557                                *chisq += dy[i] * dy[i] * sig2i;
1558                        }
1559                        break;
1560
1561                case NOISE_POISSON_DATA:
1562                        *chisq = 0.0;
1563                        /* Summation loop over all data */
1564                        for (i=fit_start; i<fit_end; i++) {
1565                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1566                                dy[i] = y[i] - yfit[i];
1567
1568                                dy_dparam[0] = 1.0;
1569                                dy_dparam[1] = exp_conv[1][i];
1570                                dy_dparam[2] = param[1] * exp_conv[2][i];
1571                                dy_dparam[3] = param[1] * exp_conv[3][i];
1572
1573                                sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15);
1574
1575                                for (j=0, l=0; l<nparam; l++) {
1576                                        if (paramfree[l]) {
1577                                                wt = dy_dparam[l] * sig2i;
1578                                                for (k=0, m=0; m<=l; m++)
1579                                                        if (paramfree[m])
1580                                                                alpha[j][k++] += wt * dy_dparam[m];
1581                                                beta[j] += wt * dy[i];
1582                                                j++;
1583                                        }
1584                                }
1585                                /* And find chi^2 */
1586                                *chisq += dy[i] * dy[i] * sig2i;
1587                        }
1588                        break;
1589
1590                case NOISE_POISSON_FIT:
1591                        *chisq = 0.0;
1592                        /* Summation loop over all data */
1593                        for (i=fit_start; i<fit_end; i++) {
1594                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1595                                dy[i] = y[i] - yfit[i];
1596
1597                                dy_dparam[0] = 1.0;
1598                                dy_dparam[1] = exp_conv[1][i];
1599                                dy_dparam[2] = param[1] * exp_conv[2][i];
1600                                dy_dparam[3] = param[1] * exp_conv[3][i];
1601
1602                                sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15);
1603
1604                                for (j=0, l=0; l<nparam; l++) {
1605                                        if (paramfree[l]) {
1606                                                wt = dy_dparam[l] * sig2i;
1607                                                for (k=0, m=0; m<=l; m++)
1608                                                        if (paramfree[m])
1609                                                                alpha[j][k++] += wt * dy_dparam[m];
1610                                                beta[j] += wt * dy[i];
1611                                                j++;
1612                                        }
1613                                }
1614                                /* And find chi^2 */
1615                                *chisq += dy[i] * dy[i] * sig2i;
1616                        }
1617                        break;
1618
1619                case NOISE_GAUSSIAN_FIT:
1620                        *chisq = 0.0;
1621                        /* Summation loop over all data */
1622                        for (i=fit_start; i<fit_end; i++) {
1623                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1624                                dy[i] = y[i] - yfit[i];
1625
1626                                dy_dparam[0] = 1.0;
1627                                dy_dparam[1] = exp_conv[1][i];
1628                                dy_dparam[2] = param[1] * exp_conv[2][i];
1629                                dy_dparam[3] = param[1] * exp_conv[3][i];
1630
1631                                sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1632
1633                                for (j=0, l=0; l<nparam; l++) {
1634                                        if (paramfree[l]) {
1635                                                wt = dy_dparam[l] * sig2i;
1636                                                for (k=0, m=0; m<=l; m++)
1637                                                        if (paramfree[m])
1638                                                                alpha[j][k++] += wt * dy_dparam[m];
1639                                                beta[j] += wt * dy[i];
1640                                                j++;
1641                                        }
1642                                }
1643                                /* And find chi^2 */
1644                                *chisq += dy[i] * dy[i] * sig2i;
1645                        }
1646                        break;
1647
1648                case NOISE_MLE:
1649                        *chisq = 0.0;
1650                        /* Summation loop over all data */
1651                        for (i=fit_start; i<fit_end; i++) {
1652                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1653                                dy[i] = y[i] - yfit[i];
1654
1655                                dy_dparam[0] = 1.0;
1656                                dy_dparam[1] = exp_conv[1][i];
1657                                dy_dparam[2] = param[1] * exp_conv[2][i];
1658                                dy_dparam[3] = param[1] * exp_conv[3][i];
1659
1660                                sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1661                                y_ymod=y[i]/yfit[i];
1662
1663                                for (j=0, l=0; l<nparam; l++) {
1664                                        if (paramfree[l]) {
1665                                                wt = dy_dparam[l] * sig2i;
1666                                                for (k=0, m=0; m<=l; m++)
1667                                                        if (paramfree[m])
1668                                                                alpha[j][k++] += wt*dy_dparam[m]*y_ymod; //wt * dy_dparam[m];
1669                                                beta[j] += wt * dy[i];
1670                                                j++;
1671                                        }
1672                                }
1673                                // And find chi^2
1674                                if (yfit[i]<=0.0)
1675                                        ; // do nothing
1676                                else if (y[i]==0.0)
1677                                        *chisq += 2.0*yfit[i];   // to avoid NaN from log
1678                                else
1679                                        *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1680                        }
1681                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1682                        break;
1683
1684                default:
1685                        return -3;
1686                        /* break; */   // (unreachable code)
1687                }
1688                break;
1689
1690        default:
1691                dbgprintf(1, "compute_exps_fn: please update me!\n");
1692                return -1;
1693        }
1694
1695        /* Fill in the symmetric side */
1696        for (j=1; j<mfit; j++)
1697                for (k=0; k<j; k++)
1698                        alpha[k][j] = alpha[j][k];
1699
1700        return 0;
1701}
1702
1703
1704/* And this is a final variant which computes the true chi-squared
1705   values and the full fit, as in EcfSingle.c */
1706int GCI_marquardt_global_compute_exps_fn_final(
1707                        float xincr, float y[],
1708                        int ndata, int fit_start, int fit_end,
1709                        noise_type noise, float sig[], int ftype,
1710                        float param[], int paramfree[], int nparam,
1711                        float *exp_conv[],
1712                        float yfit[], float dy[], float *chisq)
1713{
1714        int i, j, mfit;
1715        float sig2i;
1716
1717        for (j=0, mfit=0; j<nparam; j++)
1718                if (paramfree[j])
1719                        mfit++;
1720
1721        /* Calculation of the fitting data will depend upon the type of
1722           noise.  Since there's no convolution involved here, this is
1723           very easy. */
1724
1725        switch (ftype) {
1726        case FIT_GLOBAL_MULTIEXP:
1727                switch (noise) {
1728                case NOISE_CONST:
1729                        *chisq = 0.0;
1730                        /* Summation loop over all data */
1731                        for (i=0; i<ndata; i++) {
1732                                yfit[i] = param[0];  /* Z */
1733                                for (j=1; j<nparam; j+=2) {
1734                                        yfit[i] += param[j] * exp_conv[j+1][i];
1735                                        /* A_j . exp(-x/tau_j) */
1736                                }
1737                                dy[i] = y[i] - yfit[i];
1738
1739                                /* And find chi^2 */
1740                                if (i >= fit_start && i < fit_end)
1741                                        *chisq += dy[i] * dy[i];
1742                        }
1743
1744                        /* Now divide by sigma^2 */
1745                        sig2i = 1.0 / (sig[0] * sig[0]);
1746                        *chisq *= sig2i;
1747                        break;
1748
1749                case NOISE_GIVEN:  /* This is essentially the NR version */
1750                        *chisq = 0.0;
1751                        /* Summation loop over all data */
1752                        for (i=0; i<ndata; i++) {
1753                                yfit[i] = param[0];  /* Z */
1754                                for (j=1; j<nparam; j+=2) {
1755                                        yfit[i] += param[j] * exp_conv[j+1][i];
1756                                        /* A_j . exp(-x/tau_j) */
1757                                }
1758                                dy[i] = y[i] - yfit[i];
1759
1760                                /* And find chi^2 */
1761                                if (i >= fit_start && i < fit_end) {
1762                                        sig2i = 1.0 / (sig[i] * sig[i]);
1763                                        *chisq += dy[i] * dy[i] * sig2i;
1764                                }
1765                        }
1766                        break;
1767
1768                case NOISE_POISSON_DATA:
1769                        *chisq = 0.0;
1770                        /* Summation loop over all data */
1771                        for (i=0; i<ndata; i++) {
1772                                yfit[i] = param[0];  /* Z */
1773                                for (j=1; j<nparam; j+=2) {
1774                                        yfit[i] += param[j] * exp_conv[j+1][i];
1775                                            /* A_j . exp(-x/tau_j) */
1776                                }
1777                                dy[i] = y[i] - yfit[i];
1778
1779                                /* And find chi^2 */
1780                                if (i >= fit_start && i < fit_end) {
1781                                        /* we still don't let the variance drop below 1 */
1782                                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1783                                        *chisq += dy[i] * dy[i] * sig2i;
1784                                }
1785                        }
1786                        break;
1787
1788                case NOISE_POISSON_FIT:
1789                        *chisq = 0.0;
1790                        /* Summation loop over all data */
1791                        for (i=0; i<ndata; i++) {
1792                                yfit[i] = param[0];  /* Z */
1793                                for (j=1; j<nparam; j+=2) {
1794                                        yfit[i] += param[j] * exp_conv[j+1][i];
1795                                            /* A_j . exp(-x/tau_j) */
1796                                }
1797                                dy[i] = y[i] - yfit[i];
1798
1799                                /* And find chi^2 */
1800                                if (i >= fit_start && i < fit_end) {
1801                                        /* we still don't let the variance drop below 1 */
1802                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1803                                        *chisq += dy[i] * dy[i] * sig2i;
1804                                }
1805                        }
1806                        break;
1807
1808                case NOISE_MLE:
1809                        *chisq = 0.0;
1810                        /* Summation loop over all data */
1811                        for (i=0; i<ndata; i++) {
1812                                yfit[i] = param[0];  /* Z */
1813                                for (j=1; j<nparam; j+=2) {
1814                                        yfit[i] += param[j] * exp_conv[j+1][i];
1815                                            /* A_j . exp(-x/tau_j) */
1816                                }
1817                                dy[i] = y[i] - yfit[i];
1818
1819                                /* And find chi^2 */
1820                                if (i >= fit_start && i < fit_end) {
1821//                                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1822//                                      *chisq += dy[i] * dy[i] * sig2i;
1823                                        if (yfit[i]<=0.0)
1824                                                ; // do nothing
1825                                        else if (y[i]==0.0)
1826                                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1827                                        else
1828                                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1829                                }
1830                        }
1831                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1832                        break;
1833
1834
1835                case NOISE_GAUSSIAN_FIT:
1836                        *chisq = 0.0;
1837                        /* Summation loop over all data */
1838                        for (i=0; i<ndata; i++) {
1839                                yfit[i] = param[0];  /* Z */
1840                                for (j=1; j<nparam; j+=2) {
1841                                        yfit[i] += param[j] * exp_conv[j+1][i];
1842                                            /* A_j . exp(-x/tau_j) */
1843                                }
1844                                dy[i] = y[i] - yfit[i];
1845
1846                                /* And find chi^2 */
1847                                if (i >= fit_start && i < fit_end) {
1848                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1849                                        *chisq += dy[i] * dy[i] * sig2i;
1850                                }
1851                        }
1852                        break;
1853
1854                default:
1855                        return -3;
1856                        /* break; */   // (unreachable code)
1857                }
1858                break;
1859
1860        case FIT_GLOBAL_STRETCHEDEXP:
1861                switch (noise) {
1862                case NOISE_CONST:
1863                        *chisq = 0.0;
1864                        /* Summation loop over all data */
1865                        for (i=0; i<ndata; i++) {
1866                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1867                                dy[i] = y[i] - yfit[i];
1868
1869                                /* And find chi^2 */
1870                                if (i >= fit_start && i < fit_end)
1871                                        *chisq += dy[i] * dy[i];
1872                        }
1873
1874                        /* Now divide by sigma^2 */
1875                        sig2i = 1.0 / (sig[0] * sig[0]);
1876                        *chisq *= sig2i;
1877                        break;
1878
1879                case NOISE_GIVEN:  /* This is essentially the NR version */
1880                        *chisq = 0.0;
1881                        /* Summation loop over all data */
1882                        for (i=0; i<ndata; i++) {
1883                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1884                                dy[i] = y[i] - yfit[i];
1885
1886                                /* And find chi^2 */
1887                                if (i >= fit_start && i < fit_end) {
1888                                        sig2i = 1.0 / (sig[i] * sig[i]);
1889                                        *chisq += dy[i] * dy[i] * sig2i;
1890                                }
1891                        }
1892                        break;
1893
1894                case NOISE_POISSON_DATA:
1895                        *chisq = 0.0;
1896                        /* Summation loop over all data */
1897                        for (i=0; i<ndata; i++) {
1898                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1899                                dy[i] = y[i] - yfit[i];
1900
1901                                /* And find chi^2 */
1902                                if (i >= fit_start && i < fit_end) {
1903                                        /* we still don't let the variance drop below 1 */
1904                                        sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0);
1905                                        *chisq += dy[i] * dy[i] * sig2i;
1906                                }
1907                        }
1908                        break;
1909
1910                case NOISE_POISSON_FIT:
1911                        *chisq = 0.0;
1912                        /* Summation loop over all data */
1913                        for (i=0; i<ndata; i++) {
1914                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1915                                dy[i] = y[i] - yfit[i];
1916
1917                                /* And find chi^2 */
1918                                if (i >= fit_start && i < fit_end) {
1919                                        /* we still don't let the variance drop below 1 */
1920                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1921                                        *chisq += dy[i] * dy[i] * sig2i;
1922                                }
1923                        }
1924                        break;
1925
1926                case NOISE_MLE:
1927                        *chisq = 0.0;
1928                        /* Summation loop over all data */
1929                        for (i=0; i<ndata; i++) {
1930                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1931                                dy[i] = y[i] - yfit[i];
1932
1933                                /* And find chi^2 */
1934                                if (i >= fit_start && i < fit_end) {
1935//                                      sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1936//                                      *chisq += dy[i] * dy[i] * sig2i;
1937                                        if (yfit[i]<=0.0)
1938                                                ; // do nothing
1939                                        else if (y[i]==0.0)
1940                                                *chisq += 2.0*yfit[i];   // to avoid NaN from log
1941                                        else
1942                                                *chisq += 2.0*(yfit[i]-y[i]) - 2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i;
1943                                }
1944                        }
1945                        if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all -ve
1946                        break;
1947
1948                case NOISE_GAUSSIAN_FIT:
1949                        *chisq = 0.0;
1950                        /* Summation loop over all data */
1951                        for (i=0; i<ndata; i++) {
1952                                yfit[i] = param[0] + param[1] * exp_conv[1][i];
1953                                dy[i] = y[i] - yfit[i];
1954
1955                                /* And find chi^2 */
1956                                if (i >= fit_start && i < fit_end) {
1957                                        /* we still don't let the variance drop below 1 */
1958                                        sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0);
1959                                        *chisq += dy[i] * dy[i] * sig2i;
1960                                }
1961                        }
1962                        break;
1963
1964                default:
1965                        return -3;
1966                        /* break; */   // (unreachable code)
1967                }
1968                break;
1969
1970        default:
1971                dbgprintf(1, "compute_exps_fn: please update me!\n");
1972                return -1;
1973        }
1974
1975        return 0;
1976}
1977
1978
1979/* This one does the actual global fitting for multiexponential taus.
1980   It is basically similar to the above do_fit_single function, except
1981   that now we handle the extra intricacies involved in global
1982   fitting, in particular, the much larger alpha matrix is handled in
1983   a special way. */
1984
1985int GCI_marquardt_global_exps_do_fit_instr(
1986                float xincr, float **trans, int ndata, int ntrans,
1987                int fit_start, int fit_end, float instr[], int ninstr,
1988                noise_type noise, float sig[], int ftype,
1989                float **param, int paramfree[], int nparam, restrain_type restrain,
1990                float chisq_delta, float exp_pure[], float *exp_conv[],
1991                float **fitted, float **residuals,
1992                float **covar_scratch, float **alpha_scratch,
1993                float *chisq_trans, float *chisq_global,
1994                int drop_bad_transients)
1995{
1996        float alambda, ochisq_global, *ochisq_trans;
1997        int i, k, itst, itst_max;
1998        int ret;
1999
2000        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
2001
2002        /* If there are no global parameters being fitted, we simply fit
2003           each local set. */
2004        switch (ftype) {
2005        case FIT_GLOBAL_MULTIEXP:
2006                for (i=2; i<nparam; i+=2) {
2007                        if (paramfree[i]) {
2008                                i = 0; /* sentinel value */
2009                                break;
2010                        }
2011                }
2012                break;
2013
2014        case FIT_GLOBAL_STRETCHEDEXP:
2015                for (i=2; i<nparam; i++) {
2016                        if (paramfree[i]) {
2017                                i = 0; /* sentinel value */
2018                                break;
2019                        }
2020                }
2021                break;
2022
2023        default:
2024                dbgprintf(1, "exps_do_fit_instr: please update me!\n");
2025                return -1;
2026        }
2027
2028        if (i > 0) { /* no globals to fit */
2029                if (GCI_marquardt_global_exps_calculate_exps_instr(
2030                                xincr, ndata, instr, ninstr, ftype, param[0], nparam,
2031                                exp_pure, exp_conv) != 0)
2032                        return -2;
2033
2034                *chisq_global = 0;
2035
2036                for (i=0; i<ntrans; i++) {
2037                        if (drop_bad_transients && chisq_trans[i] < 0)
2038                                continue;
2039
2040                        ret = GCI_marquardt_global_exps_do_fit_single(
2041                                        xincr, trans[i], ndata, fit_start, fit_end,
2042                                        noise, sig, ftype, param[i], paramfree, nparam, restrain,
2043////                                    exp_conv, fitted[i], residuals[i],
2044                                        chisq_delta, exp_conv, fitted[0], residuals[0],
2045                                        covar_scratch, alpha_scratch, &chisq_trans[i]);
2046                        if (ret < 0) {
2047                                if (drop_bad_transients) {
2048                                        dbgprintf(1, "In do_fit_instr, do_fit_single returned %d "
2049                                                          "for transient %d, dropping it\n", ret, i);
2050                                        chisq_trans[i] = -1;
2051                                } else {
2052                                        dbgprintf(1, "In do_fit_instr, do_fit_single returned %d "
2053                                                          "for transient %d\n", ret, i);
2054                                        return -10 + ret;
2055                                }
2056                        } else {
2057                                *chisq_global += chisq_trans[i];
2058                        }
2059                }
2060                return 0;
2061        }
2062
2063        /* If there are no free local variables to fit, we still do the
2064           global fitting, but we have to be a little careful in some of
2065           the later routines */
2066
2067        /* Now allocate all of the arrays we will need. */
2068
2069        if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL)
2070                return -1;
2071
2072        /* We now begin our standard Marquardt loop, with several
2073           modifications */
2074        alambda = -1;
2075        ret = GCI_marquardt_global_exps_global_step(
2076                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2077                                instr, ninstr, noise, sig, ftype,
2078                                param, paramfree, nparam, restrain, exp_pure, exp_conv,
2079                                fitted, residuals, chisq_trans, chisq_global,
2080                                alpha_scratch, &alambda, drop_bad_transients);
2081        if (ret != 0) {
2082                dbgprintf(1, "In do_fit_instr, first global_step returned %d\n", ret);
2083                if (ret != -1) {
2084                        /* Wasn't a memory error, so unallocate arrays */
2085                        alambda = 0.0;
2086                        GCI_marquardt_global_exps_global_step(
2087                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2088                                instr, ninstr, noise, sig, ftype,
2089                                param, paramfree, nparam, restrain, exp_pure, exp_conv,
2090                                fitted, residuals, chisq_trans, chisq_global,
2091                                alpha_scratch, &alambda, drop_bad_transients);
2092                }
2093                free(ochisq_trans);
2094                return ret;
2095        }
2096
2097        k = 1;  /* Iteration counter */
2098        itst = 0;
2099        for (;;) {
2100                dbgprintf(3, "In do_fit_instr, beginning iteration %d:\n", k);
2101                dbgprintf(3, " itst = %d, chisq_global = %.4f\n", itst, *chisq_global);
2102
2103                k++;
2104                if (k > MAXITERS) {
2105                        return -2;
2106                }
2107
2108                ochisq_global = *chisq_global;
2109                for (i=0; i<ntrans; i++)
2110                        ochisq_trans[i] = chisq_trans[i];
2111
2112                ret = GCI_marquardt_global_exps_global_step(
2113                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2114                                        instr, ninstr, noise, sig, ftype,
2115                                        param, paramfree, nparam, restrain, exp_pure, exp_conv,
2116                                        fitted, residuals, chisq_trans, chisq_global,
2117                                        alpha_scratch, &alambda, drop_bad_transients);
2118                if (ret != 0) {
2119                        dbgprintf(1, "In do_fit_instr, second global_step returned %d\n",
2120                                          ret);
2121                        /* Unallocate arrays */
2122                        alambda = 0.0;
2123                        GCI_marquardt_global_exps_global_step(
2124                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2125                                instr, ninstr, noise, sig, ftype,
2126                                param, paramfree, nparam, restrain, exp_pure, exp_conv,
2127                                fitted, residuals, chisq_trans, chisq_global,
2128                                alpha_scratch, &alambda, drop_bad_transients);
2129                        free(ochisq_trans);
2130                        return ret;
2131                }
2132
2133                if (*chisq_global > ochisq_global)
2134                        itst = 0;
2135                else {
2136                        /* Let's try this approach; I really don't know what will
2137                           be best */
2138                        float maxdiff;
2139
2140                        maxdiff = 0.0;
2141                        for (i=0; i<ntrans; i++)
2142                                if (ochisq_trans[i] - chisq_trans[i] > maxdiff)
2143                                        maxdiff = ochisq_trans[i] - chisq_trans[i];
2144
2145                        if (maxdiff < chisq_delta)
2146                                itst++;
2147                        dbgprintf(3, "In do_fit_instr, maxdiff = %.3f:\n", maxdiff);
2148                }
2149
2150                if (itst < itst_max) continue;
2151
2152                /* Endgame */
2153                alambda = 0.0;
2154                ret = GCI_marquardt_global_exps_global_step(
2155                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2156                                        instr, ninstr, noise, sig, ftype,
2157                                        param, paramfree, nparam, restrain, exp_pure, exp_conv,
2158                                        fitted, residuals, chisq_trans, chisq_global,
2159                                        alpha_scratch, &alambda, drop_bad_transients);
2160                if (ret != 0) {
2161                        dbgprintf(1, "In do_fit_instr, final global_step returned %d\n",
2162                                          ret);
2163                        free(ochisq_trans);
2164                        return ret;
2165                }
2166
2167                free(ochisq_trans);
2168                return k;  /* We're done now */
2169        }
2170}
2171
2172
2173/* And this one is basically a specialised GCI_marquardt_instr_step
2174   for the global fitting setup. */
2175int GCI_marquardt_global_exps_global_step(
2176                                float xincr, float **trans,
2177                                int ndata, int ntrans, int fit_start, int fit_end,
2178                                float instr[], int ninstr,
2179                                noise_type noise, float sig[], int ftype,
2180                                float **param, int paramfree[], int nparam,
2181                                restrain_type restrain,
2182                                float exp_pure[], float *exp_conv[],
2183                                float **yfit, float **dy,
2184                                float *chisq_trans, float *chisq_global,
2185                                float **alpha_scratch, float *alambda,
2186                                int drop_bad_transients)
2187{
2188        int i, j, ret;
2189        static global_matrix alpha, covar;
2190        static global_vector beta, dparam;
2191        static float **paramtry;
2192        static int mfit_local, mfit_global;
2193        static int gindex[MAXFIT], lindex[MAXFIT];
2194        static float ochisq_global, *ochisq_trans;
2195        static void (*fitfunc)(float, float [], float *, float [], int);
2196        static int initialised=0;
2197
2198        if (nparam > MAXFIT)
2199                return -10;
2200        if (xincr <= 0)
2201                return -11;
2202        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
2203                return -12;
2204
2205        /* Initialisation */
2206        /* We assume we're given sensible starting values for param[] */
2207        if (*alambda < 0.0) {
2208                /* Start by allocating lots of variables we will need */
2209                mfit_local = mfit_global = 0;
2210
2211                switch (ftype) {
2212                case FIT_GLOBAL_MULTIEXP:
2213                        fitfunc = GCI_multiexp_tau;
2214
2215                        /* We know that all of param[2i], the taus, are the global
2216                           variables, and that the param[2i+1], the A's, are the
2217                           local variables, along with param[0] = Z.  We stored
2218                           the indices of the local and global free variables in
2219                           lindex and gindex respectively. */
2220                        if (paramfree[0]) {
2221                                lindex[mfit_local++] = 0;
2222                        }
2223                        for (i=1; i<nparam; i+=2) {
2224                                if (paramfree[i])
2225                                        lindex[mfit_local++] = i;
2226                        }
2227                        for (i=2; i<nparam; i+=2) {
2228                                if (paramfree[i])
2229                                        gindex[mfit_global++] = i;
2230                        }
2231                        break;
2232
2233                case FIT_GLOBAL_STRETCHEDEXP:
2234                        fitfunc = GCI_stretchedexp;
2235
2236                        if (paramfree[0])
2237                                lindex[mfit_local++] = 0;  /* Z */
2238                        if (paramfree[1])
2239                                lindex[mfit_local++] = 1;  /* A */
2240                        if (paramfree[2])
2241                                gindex[mfit_global++] = 2;  /* tau */
2242                        if (paramfree[3])
2243                                gindex[mfit_global++] = 3;  /* h */
2244                        break;
2245
2246                default:
2247                        dbgprintf(1, "exps_global_step: please update me!\n");
2248                        return -1;
2249                }
2250
2251                if (initialised) {
2252                        GCI_free_global_matrix(&alpha);   GCI_free_global_matrix(&covar);
2253                        GCI_ecf_free_matrix(paramtry);        GCI_free_global_vector(&beta);
2254                        GCI_free_global_vector(&dparam);  free(ochisq_trans);
2255                        initialised = 0;
2256                }
2257
2258                if (GCI_alloc_global_matrix(&alpha, mfit_global, mfit_local, ntrans)
2259                        != 0)
2260                        return -1;
2261
2262                if (GCI_alloc_global_matrix(&covar, mfit_global, mfit_local, ntrans)
2263                        != 0) {
2264                        GCI_free_global_matrix(&alpha);
2265                        return -1;
2266                }
2267
2268                if ((paramtry = GCI_ecf_matrix(ntrans, nparam)) == NULL) {
2269                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2270                        return -1;
2271                }
2272
2273                if (GCI_alloc_global_vector(&beta, mfit_global, mfit_local, ntrans)
2274                        != 0) {
2275                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2276                        GCI_ecf_free_matrix(paramtry);
2277                        return -1;
2278                }
2279
2280                if (GCI_alloc_global_vector(&dparam, mfit_global, mfit_local, ntrans)
2281                        != 0) {
2282                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2283                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
2284                        return -1;
2285                }
2286
2287                if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float)))
2288                        == NULL) {
2289                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2290                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
2291                        GCI_free_global_vector(&dparam);
2292                        return -1;
2293                }
2294
2295                initialised = 1;
2296
2297                if (GCI_marquardt_global_compute_global_exps_fn(
2298                                        xincr, trans, ndata, ntrans,
2299                                        fit_start, fit_end, instr, ninstr, noise, sig, ftype,
2300                                        param, paramfree, nparam, mfit_global, mfit_local,
2301                                        gindex, lindex, exp_pure, exp_conv,
2302                                        yfit, dy, alpha, beta, alpha_scratch,
2303                                        chisq_trans, chisq_global, drop_bad_transients) != 0)
2304                        return -2;
2305
2306                *alambda = 0.001;
2307                ochisq_global = *chisq_global;
2308                for (i=0; i<ntrans; i++)
2309                        ochisq_trans[i] = chisq_trans[i];
2310
2311                /* Initialise paramtry to param */
2312                for (i=0; i<ntrans; i++) {
2313                        for (j=0; j<nparam; j++)
2314                                paramtry[i][j] = param[i][j];
2315                }
2316        }
2317
2318        /* Once converged, evaluate covariance matrix */
2319        if (*alambda == 0) {
2320                if (GCI_marquardt_global_compute_global_exps_fn_final(
2321                                        xincr, trans, ndata, ntrans,
2322                                        fit_start, fit_end, instr, ninstr, noise, sig, ftype,
2323                                        param, paramfree, nparam, mfit_global, mfit_local,
2324                                        gindex, lindex, exp_pure, exp_conv, yfit, dy,
2325                                        chisq_trans, chisq_global, drop_bad_transients) != 0)
2326                        return -3;
2327                /* Don't need to do this here; if we wished to, we'd have to
2328                   move this code (the "if (*alambda == 0)" block) to after
2329                   the Gauss-Jordan call.  We'd also need to rewrite it for
2330                   our situation.... */ 
2331                //      if (mfit < nparam) {  /* no need to do this otherwise */
2332                //              GCI_covar_sort(covar, nparam, paramfree, mfit);
2333                //              GCI_covar_sort(alpha, nparam, paramfree, mfit);
2334                //      }
2335                GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
2336                GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
2337                GCI_free_global_vector(&dparam); free(ochisq_trans);
2338                initialised = 0;
2339                return 0;
2340        }
2341
2342        /* Alter linearised fitting matrix by augmenting diagonal
2343           elements. */
2344        GCI_copy_global_matrix(covar, alpha, mfit_global, mfit_local, ntrans);
2345        GCI_copy_global_vector(dparam, beta, mfit_global, mfit_local, ntrans);
2346        for (j=0; j<mfit_global; j++)
2347                covar.P[j][j] *= 1.0 + (*alambda);
2348        for (i=0; i<ntrans; i++)
2349                for (j=0; j<mfit_local; j++)
2350                        covar.S[i][j][j] *= 1.0 + (*alambda);
2351
2352        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */
2353        if (GCI_marquardt_global_solve_eqn(covar, dparam,
2354                                                                           mfit_global, mfit_local, ntrans) != 0)
2355                return -3;
2356
2357        /* Did the trial succeed?  Modify param by dparam... */
2358        for (i=0; i<ntrans; i++) {
2359                for (j=0; j<mfit_global; j++)
2360                        paramtry[i][gindex[j]] = param[i][gindex[j]] + dparam.global[j];
2361                for (j=0; j<mfit_local; j++)
2362                        paramtry[i][lindex[j]] =
2363                                param[i][lindex[j]] + dparam.local[i*mfit_local + j];
2364        }
2365
2366        for (i=0; i<ntrans; i++) {
2367                if (drop_bad_transients && chisq_trans[i] < 0)
2368                        continue;
2369
2370                if (restrain == ECF_RESTRAIN_DEFAULT)
2371                        ret = check_ecf_params (paramtry[i], nparam, fitfunc);
2372                else
2373                        ret = check_ecf_user_params (paramtry[i], nparam, fitfunc);
2374
2375                if (ret != 0) {
2376                        /* Bad parameters, increase alambda and return */
2377                        *alambda *= 10.0;
2378                        return 0;
2379                }
2380        }
2381
2382        if (GCI_marquardt_global_compute_global_exps_fn(
2383                                xincr, trans, ndata, ntrans,
2384                                fit_start, fit_end, instr, ninstr, noise, sig, ftype,
2385                                paramtry, paramfree, nparam, mfit_global, mfit_local,
2386                                gindex, lindex, exp_pure, exp_conv,
2387                                yfit, dy, covar, dparam, alpha_scratch,
2388                                chisq_trans, chisq_global, drop_bad_transients) != 0)
2389                return -2;
2390
2391        /* Success, accept the new solution */
2392        if (*chisq_global < ochisq_global) {
2393                *alambda *= 0.1;
2394                ochisq_global = *chisq_global;
2395                for (i=0; i<ntrans; i++)
2396                        ochisq_trans[i] = chisq_trans[i];
2397                GCI_copy_global_matrix(alpha, covar, mfit_global, mfit_local, ntrans);
2398                GCI_copy_global_vector(beta, dparam, mfit_global, mfit_local, ntrans);
2399                for (i=0; i<ntrans; i++) {
2400                        for (j=0; j<nparam; j++)
2401                                param[i][j] = paramtry[i][j];
2402                }
2403        } else { /* Failure, increase alambda and return */
2404                *alambda *= 10.0;
2405                *chisq_global = ochisq_global;
2406                for (i=0; i<ntrans; i++)
2407                        chisq_trans[i] = ochisq_trans[i];
2408        }
2409
2410        return 0;
2411}
2412
2413
2414/* Here we use alpha only for scratch space */
2415int GCI_marquardt_global_compute_global_exps_fn(
2416                float xincr, float **trans, int ndata, int ntrans,
2417                int fit_start, int fit_end, float instr[], int ninstr,
2418                noise_type noise, float sig[], int ftype,
2419                float **param, int paramfree[], int nparam,
2420                int mfit_global, int mfit_local, int gindex[], int lindex[],
2421                float exp_pure[], float *exp_conv[],
2422                float **yfit, float **dy, global_matrix alpha, global_vector beta,
2423                float **alpha_scratch, float *chisq_trans, float *chisq_global,
2424                int drop_bad_transients)
2425{
2426        int i, j, k, ret;
2427        float beta_scratch[MAXFIT];  /* scratch space */
2428
2429        /* Calculate the exponential array once only */
2430        if (GCI_marquardt_global_exps_calculate_exps_instr(
2431                        xincr, ndata, instr, ninstr, ftype, param[0], nparam,
2432                        exp_pure, exp_conv) != 0)
2433                return -1;
2434
2435        /* We initialise P and beta_global to zero; the others don't
2436           matter, as they will be totally overwritten */
2437        for (i=0; i<mfit_global; i++) {
2438                for (j=0; j<mfit_global; j++)
2439                        alpha.P[i][j] = 0;
2440                beta.global[i] = 0;
2441        }
2442        *chisq_global = 0.0;
2443
2444        for (i=0; i<ntrans; i++) {
2445                if (drop_bad_transients && chisq_trans[i] < 0) {
2446                        for (j=0; j<mfit_global; j++)
2447                                for (k=0; k<mfit_local; k++)
2448                                        alpha.Q[j][i*mfit_local + k] = 0.0;
2449                        for (j=0; j<mfit_local; j++) {
2450                                /* Make this component of S an identity matrix and of
2451                                   beta zero */
2452                                for (k=0; k<mfit_local; k++)
2453                                        alpha.S[i][j][k] = (j == k) ? 1.0 : 0.0;
2454                                beta.local[i*mfit_local + j] = 0;
2455                        }
2456                        continue;
2457                }
2458                       
2459                /* This transient is fine! */
2460                ret = GCI_marquardt_global_compute_exps_fn(
2461                                        xincr, trans[i], ndata, fit_start, fit_end, noise, sig,
2462                                        ftype, param[i], paramfree, nparam,
2463////                                    exp_conv, yfit[i], dy[i],
2464                                        exp_conv, yfit[0], dy[0],
2465                                        alpha_scratch, beta_scratch, &chisq_trans[i]);
2466
2467                if (ret != 0) {
2468                        if (drop_bad_transients) {
2469                                dbgprintf(3, "In compute_global_exps_fn, "
2470                                                  "compute_exps_fn returned %d for transient %d\n",
2471                                                  ret, i);
2472                                chisq_trans[i] = -1;
2473                                continue;
2474                        } else {
2475                                dbgprintf(1, "In compute_global_exps_fn, "
2476                                                  "compute_exps_fn returned %d for transient %d\n",
2477                                                  ret, i);
2478                                return -2;
2479                        }
2480                }
2481
2482                /* So now have to populate alpha and beta with the contents of
2483                   alpha_scratch and beta_scratch. */
2484
2485                for (j=0; j<mfit_global; j++) {
2486                        for (k=0; k<mfit_global; k++)
2487                                alpha.P[j][k] += alpha_scratch[gindex[j]][gindex[k]];
2488                        for (k=0; k<mfit_local; k++)
2489                                alpha.Q[j][i*mfit_local + k] =
2490                                        alpha_scratch[gindex[j]][lindex[k]];
2491                        beta.global[j] += beta_scratch[gindex[j]];
2492                }
2493                for (j=0; j<mfit_local; j++) {
2494                        for (k=0; k<mfit_local; k++)
2495                                alpha.S[i][j][k] = alpha_scratch[lindex[j]][lindex[k]];
2496                        beta.local[i*mfit_local + j] = beta_scratch[lindex[j]];
2497                }
2498
2499                *chisq_global += chisq_trans[i];
2500        }
2501
2502        return 0;
2503}
2504
2505
2506/* The final variant */
2507int GCI_marquardt_global_compute_global_exps_fn_final(
2508                float xincr, float **trans, int ndata, int ntrans,
2509                int fit_start, int fit_end, float instr[], int ninstr,
2510                noise_type noise, float sig[], int ftype,
2511                float **param, int paramfree[], int nparam,
2512                int mfit_global, int mfit_local, int gindex[], int lindex[],
2513                float exp_pure[], float *exp_conv[],
2514                float **yfit, float **dy,
2515                float *chisq_trans, float *chisq_global, int drop_bad_transients)
2516{
2517        int i, ret;
2518
2519        /* Calculate the exponential array once only */
2520        if (GCI_marquardt_global_exps_calculate_exps_instr(
2521                        xincr, ndata, instr, ninstr, ftype, param[0], nparam,
2522                        exp_pure, exp_conv) != 0)
2523                return -1;
2524
2525        *chisq_global = 0.0;
2526
2527        for (i=0; i<ntrans; i++) {
2528                if (drop_bad_transients && chisq_trans[i] < 0)
2529                        continue;
2530                       
2531                /* This transient is fine! */
2532                ret = GCI_marquardt_global_compute_exps_fn_final(
2533                                        xincr, trans[i], ndata, fit_start, fit_end, noise, sig,
2534                                        ftype, param[i], paramfree, nparam,
2535////                                    exp_conv, yfit[i], dy[i], &chisq_trans[i]);
2536                                        exp_conv, yfit[0], dy[0], &chisq_trans[i]);
2537
2538                if (ret != 0) {
2539                        if (drop_bad_transients) {
2540                                dbgprintf(3, "In compute_global_exps_fn_final, "
2541                                                  "compute_exps_fn_final returned %d "
2542                                                  "for transient %d\n",
2543                                                  ret, i);
2544                                chisq_trans[i] = -1;
2545                                continue;
2546                        } else {
2547                                dbgprintf(1, "In compute_global_exps_fn_final, "
2548                                                  "compute_exps_fn_final returned %d "
2549                                                  "for transient %d\n",
2550                                                  ret, i);
2551                                return -2;
2552                        }
2553                }
2554
2555                *chisq_global += chisq_trans[i];
2556        }
2557
2558        return 0;
2559}
2560
2561
2562/* This function solves the equation Ax=b, where A is the alpha
2563   matrix, which has the form:
2564
2565     A = (P Q)
2566         (R S)
2567
2568   Here P is an mfit_global x mfit_global square matrix, S is a block
2569   diagonal matrix with ntrans blocks, each of size mfit_local x
2570   mfit_local, and Q and R are the right sizes to make the whole
2571   matrix square.  We solve it by inverting the matrix A using the
2572   formulae given in Numerical Recipes section 2.7, then multiplying
2573   the inverse by b to get x.  We are not too concerned about
2574   accuracy, as this does not affect the solution found, only the
2575   route taken to find it.
2576
2577   Numerical Recipes, section 2.7, notes that A^{-1} is given by:
2578
2579      (P' Q')
2580      (R' S')
2581
2582   with:
2583
2584      P' = (P - Q.S^{-1}.R)^{-1}
2585      Q' = - (P - Q.S^{-1}.R)^{-1} . (Q.S^{-1})
2586      R' = - (S^{-1}.R) . (P - Q.S^{-1}.R)^{-1}
2587      S' = S^{-1} + (S^{-1}.R) . (P - Q.S^{-1}.R)^{-1} . (Q.S^{-1})
2588
2589   We also make use of the fact that A is symmetric, so in particular,
2590   (S^{-1}.R) = (Q.S^{-1})^T and R' = Q'^T.
2591
2592   We are given A as a global_matrix and b as a global_vector.  This
2593   function destroys the original matrix and returns the solution in
2594   place of b.
2595*/
2596
2597int GCI_marquardt_global_solve_eqn(global_matrix A, global_vector b,
2598                        int mfit_global, int mfit_local, int ntrans)
2599{
2600        int row, col, block, i, j;
2601        float x_temp[MAXFIT], x_temp2[MAXFIT];
2602        static float **QS;
2603        static global_vector x;
2604        static int saved_global=0, saved_local=0, saved_ntrans=0;
2605
2606        /* If no local parameters, just do a straight Gauss-Jordan method */
2607        if (mfit_local == 0) {
2608                if (GCI_gauss_jordan(A.P, mfit_global, b.global) != 0)
2609                        return -2;
2610                return 0;
2611        }
2612
2613        /* Allocate arrays if necessary */
2614        if ((saved_global != mfit_global) || (saved_local != mfit_local) ||
2615                (saved_ntrans != ntrans)) {
2616                if (saved_global > 0) {
2617                        GCI_ecf_free_matrix(QS);
2618                        GCI_free_global_vector(&x);
2619                        saved_global = 0;
2620                }
2621                if ((QS = GCI_ecf_matrix(mfit_global, mfit_local*ntrans)) == NULL)
2622                        return -1;
2623                if (GCI_alloc_global_vector(&x, mfit_global, mfit_local, ntrans)
2624                        != 0) {
2625                        GCI_ecf_free_matrix(QS);
2626                        return -1;
2627                }
2628                saved_global = mfit_global;
2629                saved_local = mfit_local;
2630                saved_ntrans = ntrans;
2631        }
2632
2633        /* Start by inverting S */
2634        for (block=0; block<ntrans; block++)
2635                if (GCI_gauss_jordan(A.S[block], mfit_local, NULL) != 0)
2636                        return -2;
2637
2638        /* Calculate Q.S^{-1} */
2639        for (row=0; row<mfit_global; row++)
2640                for (block=0; block<ntrans; block++)
2641                        for (i=0; i<mfit_local; i++) {
2642                                QS[row][block*mfit_local + i] = 0;
2643                                for (j=0; j<mfit_local; j++)
2644                                        QS[row][block*mfit_local + i] +=
2645                                                A.Q[row][block*mfit_local + j] * A.S[block][j][i];
2646                        }
2647
2648        /* Now find P - Q.S^{-1}.R */
2649        for (row=0; row<mfit_global; row++)
2650                for (col=0; col<mfit_global; col++)
2651                        for (i=0; i<ntrans*mfit_local; i++)
2652                                A.P[row][col] -= QS[row][i] * A.Q[col][i];  /* Q = R^T */
2653
2654        /* And invert it to get P' */
2655        if (GCI_gauss_jordan(A.P, mfit_global, NULL) != 0)
2656                return -3;
2657
2658        /* Now overwrite Q with Q' */
2659        for (row=0; row<mfit_global; row++)
2660                for (col=0; col<ntrans*mfit_local; col++) {
2661                        A.Q[row][col] = 0;
2662                        for (i=0; i<mfit_global; i++)
2663                                A.Q[row][col] -= A.P[row][i] * QS[i][col];  /* P contains P' */
2664                }
2665
2666        /* Finally, we can solve to find x */
2667        /* We have x.global = P'.(b.global) + Q'.(b.local)
2668           and     x.local  = R'.(b.global) + S'.(b.local)
2669           We do x.global first. */
2670        for (row=0; row<mfit_global; row++) {
2671                x.global[row] = 0;
2672                for (i=0; i<mfit_global; i++)
2673                        x.global[row] += A.P[row][i] * b.global[i];
2674                /* Recall that Q now contains Q' */
2675                for (i=0; i < ntrans * mfit_local; i++)
2676                        x.global[row] += A.Q[row][i] * b.local[i];
2677        }
2678
2679        /* Now x_local; the R'.b_global component first, recalling that R'
2680           is Q' transposed and that Q' is stored in Q. */
2681        for (row=0; row < ntrans * mfit_local; row++) {
2682                x.local[row] = 0;
2683                for (i=0; i<mfit_global; i++)
2684                        x.local[row] += A.Q[i][row] * b.global[i];
2685        }
2686
2687        /* Now S' = S^{-1} + (S^{-1}.R).(P-Q.S^{-1}.R)^{-1}.(Q.S^{-1}).
2688           We first handle the S^{-1} term, then the remaining term.  */
2689        for (block=0; block<ntrans; block++)
2690                for (row=0; row<mfit_local; row++) {
2691                        for (j=0; j<mfit_local; j++)
2692                                x.local[block*mfit_local + row] +=
2693                                        A.S[block][row][j] * b.local[block*mfit_local + j];
2694                }
2695
2696        /* For the remaining term, we have an x.local[row] contribution of
2697
2698              sum_j sum_k sum_l
2699                  (S^{-1}.R)_{row,j} . (P-Q.S^{-1}.R)^{-1}_{j,k} .
2700                  (Q.S^{-1})_{k,l} . b.local_{l}
2701
2702           In order to save computations, we calculate the matrices once
2703           only.  We start with (Q.S^{-1}) . b_local, which is an
2704           mfit_global x 1 array, and store this in x_temp; premultiplying
2705           this by the square matrix P' again gives an mfit_global x 1
2706           array, which goes in x_temp2, then premultiplying this by
2707           (S^{-1}.R) gives an (ntrans * mfit_local) x 1 array which is
2708           added directly onto x_local.  Recall also that S^{-1}.R is the
2709           transpose of Q.S^{-1}, which is currently stored in QS, and
2710           that the middle term is currently stored in P. */
2711
2712        for (row=0; row<mfit_global; row++) {
2713                x_temp[row] = 0;
2714                for (i=0; i < ntrans*mfit_local; i++)
2715                        x_temp[row] += QS[row][i] * b.local[i];
2716        }
2717        for (row=0; row<mfit_global; row++) {
2718                x_temp2[row] = 0;
2719                for (i=0; i<mfit_global; i++)
2720                        x_temp2[row] += A.P[row][i] * x_temp[i];
2721        }
2722        /* Again, S^{-1}.R is the transpose of Q.S^{-1} */
2723        for (row=0; row < ntrans * mfit_local; row++) {
2724                for (i=0; i<mfit_global; i++)
2725                        x.local[row] += QS[i][row] * x_temp2[i];
2726        }
2727
2728        /* And we're done, once we've copied x into b */
2729        GCI_copy_global_vector(b, x, mfit_global, mfit_local, ntrans);
2730
2731        return 0;
2732}
2733
2734
2735/*         *****        GENERIC FUNCTION CODE        *****        */
2736
2737/* These functions are essentially the same as the above functions
2738   GCI_marquardt_global_exps_instr and
2739   GCI_marquardt_global_exps_do_fit_instr and the latter's dependents,
2740   except that this version takes an arbitrary function and a list of
2741   global parameters.  This function is designed to be called from
2742   external code.  It is nowhere near as efficient as the streamlined
2743   code above for globally fitting taus for multi-exponential models.
2744   Also, this function must be provided with meaningful starting
2745   estimates for all parameters.
2746   */
2747   
2748int GCI_marquardt_global_generic_instr(float xincr, float **trans,
2749                                        int ndata, int ntrans, int fit_start, int fit_end,
2750                                        float instr[], int ninstr,
2751                                        noise_type noise, float sig[],
2752                                        float **param, int paramfree[], int nparam, int gparam[],
2753                                        restrain_type restrain, float chisq_delta,
2754                                        void (*fitfunc)(float, float [], float *, float [], int),
2755                                        float **fitted, float **residuals,
2756                                        float chisq_trans[], float *chisq_global, int *df)
2757{
2758        float **covar, **alpha, *scaled_instr, instrsum;
2759        int i, ret;
2760        int mlocal, mglobal;
2761
2762        /* Some basic parameter checks */
2763        if (xincr <= 0) return -1;
2764        if (ntrans < 1) return -1;
2765        if (ndata < 1)  return -1;
2766        if (fit_start < 0 || fit_end > ndata) return -1;
2767        if (ninstr < 1) return -1;
2768        if (nparam < 1) return -1;
2769
2770        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL)
2771                return -2;
2772
2773        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
2774                GCI_ecf_free_matrix(covar);
2775                return -3;
2776        }
2777
2778        if ((scaled_instr = (float *) malloc(ninstr * sizeof(float))) == NULL) {
2779                GCI_ecf_free_matrix(covar);
2780                GCI_ecf_free_matrix(alpha);
2781                return -4;
2782        }
2783
2784        /* Scale the instrument response */
2785        for (i=0, instrsum=0; i<ninstr; i++)
2786                instrsum += instr[i];
2787        if (instrsum == 0) return -6;
2788        for (i=0; i<ninstr; i++)
2789                scaled_instr[i] = instr[i] / instrsum;
2790
2791        /* Now call the global fitting function. */
2792        ret = GCI_marquardt_global_generic_do_fit_instr(
2793                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2794                        scaled_instr, ninstr, noise, sig,
2795                        param, paramfree, nparam, gparam, restrain, chisq_delta,
2796                        fitfunc, fitted, residuals, covar, alpha,
2797                        chisq_trans, chisq_global);
2798
2799        GCI_ecf_free_matrix(covar);
2800        GCI_ecf_free_matrix(alpha);
2801        free(scaled_instr);
2802        GCI_marquardt_cleanup();
2803
2804        if (ret < 0) {
2805                dbgprintf(1, "Fit failed, ret = %d\n", ret);
2806                return -10 + ret;
2807        }
2808
2809        dbgprintf(2, "Fit succeeded, ret = %d\n", ret);
2810
2811        /* Before we return, calculate the number of degrees of freedom */
2812        /* The number of degrees of freedom is given by:
2813              d.f. = ntrans * ((fit_end - fit_start) - # free local parameters)
2814                         - # free global parameters
2815        */
2816
2817        mglobal = mlocal = 0;
2818        for (i=0; i<nparam; i++)
2819                if (paramfree[i]) {
2820                        if (gparam[i]) mglobal++;
2821                        else mlocal++;
2822                }
2823
2824        *df = ntrans * ((fit_end - fit_start) - mlocal) - mglobal;
2825
2826        return ret;
2827}
2828
2829
2830int GCI_marquardt_global_generic_do_fit_instr(
2831                float xincr, float **trans, int ndata, int ntrans,
2832                int fit_start, int fit_end, float instr[], int ninstr,
2833                noise_type noise, float sig[],
2834                float **param, int paramfree[], int nparam, int gparam[],
2835                restrain_type restrain, float chisq_delta,
2836                void (*fitfunc)(float, float [], float *, float [], int),
2837                float **fitted, float **residuals,
2838                float **covar_scratch, float **alpha_scratch,
2839                float *chisq_trans, float *chisq_global)
2840{
2841        // PRB 03/07 Although **fitted and **residuals are provided only one "transient" is required and used, fitted[0] and residuals[0]
2842       
2843        float alambda, ochisq_global, *ochisq_trans;
2844        int i, k, itst, itst_max;
2845        int ret;
2846
2847        itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6;
2848
2849        /* If there are no global parameters being fitted, we simply fit
2850           each local set. */
2851        for (i=0; i<nparam; i++) {
2852                if (gparam[i] && paramfree[i]) {
2853                        i = -1; /* sentinel value */
2854                        break;
2855                }
2856        }
2857
2858        if (i >= 0) { /* no globals to fit */
2859                *chisq_global = 0;
2860
2861                for (i=0; i<ntrans; i++) {
2862                        ret = GCI_marquardt_instr(xincr, trans[i],
2863                                        ndata, fit_start, fit_end, instr, ninstr, noise, sig,
2864                                        param[i], paramfree, nparam, restrain, fitfunc,
2865                                        fitted[0], residuals[0], covar_scratch, alpha_scratch,
2866                                        &chisq_trans[i], chisq_delta, 0, NULL);
2867                        if (ret < 0) {
2868                                dbgprintf(1, "In do_fit_instr, marquardt_instr returned %d "
2869                                                  "for transient %d\n", ret, i);
2870                                return -10 + ret;
2871                        } else {
2872                                *chisq_global += chisq_trans[i];
2873                        }
2874                }
2875                return 0;
2876        }
2877
2878        /* If there are no free local variables to fit, we still do the
2879           global fitting, but we have to be a little careful in some of
2880           the later routines */
2881
2882        /* Now allocate all of the arrays we will need. */
2883
2884        if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL)
2885                return -1;
2886
2887        /* We now begin our standard Marquardt loop, with several
2888           modifications */
2889        alambda = -1;
2890        ret = GCI_marquardt_global_generic_global_step(
2891                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2892                                instr, ninstr, noise, sig,
2893                                param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2894                                fitted, residuals, chisq_trans, chisq_global,
2895                                alpha_scratch, &alambda);
2896        if (ret != 0) {
2897                dbgprintf(1, "In do_fit_instr, first global_step returned %d\n", ret);
2898                if (ret != -1) {
2899                        /* Wasn't a memory error, so unallocate arrays */
2900                        alambda = 0.0;
2901                        GCI_marquardt_global_generic_global_step(
2902                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2903                                instr, ninstr, noise, sig,
2904                                param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2905                                fitted, residuals, chisq_trans, chisq_global,
2906                                alpha_scratch, &alambda);
2907                }
2908                free(ochisq_trans);
2909                return ret;
2910        }
2911
2912        k = 1;  /* Iteration counter */
2913        itst = 0;
2914        for (;;) {
2915                dbgprintf(3, "In do_fit_instr, beginning iteration %d:\n", k);
2916                dbgprintf(3, " itst = %d, chisq_global = %.4f\n", itst, *chisq_global);
2917
2918                k++;
2919                if (k > MAXITERS) {
2920                        return -2;
2921                }
2922
2923                ochisq_global = *chisq_global;
2924                for (i=0; i<ntrans; i++)
2925                        ochisq_trans[i] = chisq_trans[i];
2926
2927                ret = GCI_marquardt_global_generic_global_step(
2928                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2929                                        instr, ninstr, noise, sig,
2930                                        param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2931                                        fitted, residuals, chisq_trans, chisq_global,
2932                                        alpha_scratch, &alambda);
2933                if (ret != 0) {
2934                        dbgprintf(1, "In do_fit_instr, second global_step returned %d\n",
2935                                          ret);
2936                        /* Unallocate arrays */
2937                        alambda = 0.0;
2938                        GCI_marquardt_global_generic_global_step(
2939                                xincr, trans, ndata, ntrans, fit_start, fit_end,
2940                                instr, ninstr, noise, sig,
2941                                param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2942                                fitted, residuals, chisq_trans, chisq_global,
2943                                alpha_scratch, &alambda);
2944                        free(ochisq_trans);
2945                        return ret;
2946                }
2947
2948                if (*chisq_global > ochisq_global)
2949                        itst = 0;
2950                else {
2951                        /* Let's try this approach; I really don't know what will
2952                           be best */
2953                        float maxdiff;
2954
2955                        maxdiff = 0.0;
2956                        for (i=0; i<ntrans; i++)
2957                                if (ochisq_trans[i] - chisq_trans[i] > maxdiff)
2958                                        maxdiff = ochisq_trans[i] - chisq_trans[i];
2959
2960                        if (maxdiff < chisq_delta)
2961                                itst++;
2962                        dbgprintf(3, "In do_fit_instr, maxdiff = %.3f:\n", maxdiff);
2963                }
2964
2965                if (itst < itst_max) continue;
2966
2967                /* Endgame */
2968                alambda = 0.0;
2969                ret = GCI_marquardt_global_generic_global_step(
2970                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
2971                                        instr, ninstr, noise, sig,
2972                                        param, paramfree, nparam, gparam, restrain, chisq_delta, fitfunc,
2973                                        fitted, residuals, chisq_trans, chisq_global,
2974                                        alpha_scratch, &alambda);
2975                if (ret != 0) {
2976                        dbgprintf(1, "In do_fit_instr, final global_step returned %d\n",
2977                                          ret);
2978                        free(ochisq_trans);
2979                        return ret;
2980                }
2981
2982                free(ochisq_trans);
2983                return k;  /* We're done now */
2984        }
2985}
2986
2987
2988/* And this one is basically a specialised GCI_marquardt_instr_step
2989   for the global fitting setup. */
2990
2991#define do_frees \
2992        if (fnvals) free(fnvals);\
2993        if (dy_dparam_pure) GCI_ecf_free_matrix(dy_dparam_pure);\
2994        if (dy_dparam_conv) GCI_ecf_free_matrix(dy_dparam_conv);
2995
2996int GCI_marquardt_global_generic_global_step(
2997                                float xincr, float **trans,
2998                                int ndata, int ntrans, int fit_start, int fit_end,
2999                                float instr[], int ninstr,
3000                                noise_type noise, float sig[],
3001                                float **param, int paramfree[], int nparam, int gparam[],
3002                                restrain_type restrain, float chisq_delta,
3003                                void (*fitfunc)(float, float [], float *, float [], int),
3004                                float **yfit, float **dy,
3005                                float *chisq_trans, float *chisq_global,
3006                                float **alpha_scratch, float *alambda)
3007{
3008        int i, j, ret;
3009        static global_matrix alpha, covar;
3010        static global_vector beta, dparam;
3011        static float **paramtry;
3012        static int mfit_local, mfit_global;
3013        static int gindex[MAXFIT], lindex[MAXFIT];
3014        static float ochisq_global, *ochisq_trans;
3015        static int initialised=0;
3016
3017        // The following are declared here to retain some optimisation by not repeatedly mallocing
3018        // (only once per transient), but to remain thread safe.
3019        // They are malloced by lower fns but at the end, freed by this fn.
3020        // These vars were global or static before thread safety was introduced.
3021        float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL;
3022        int fnvals_len=0, dy_dparam_nparam_size=0;
3023
3024        if (nparam > MAXFIT)
3025                return -10;
3026        if (xincr <= 0)
3027                return -11;
3028        if (fit_start < 0 || fit_start > fit_end || fit_end > ndata)
3029                return -12;
3030
3031        /* Initialisation */
3032        /* We assume we're given sensible starting values for param[] */
3033        if (*alambda < 0.0) {
3034                /* Start by allocating lots of variables we will need */
3035                mfit_local = mfit_global = 0;
3036
3037                for (i=0; i<nparam; i++) {
3038                        if (paramfree[i]) {
3039                                if (gparam[i])
3040                                        gindex[mfit_global++] = i;
3041                                else
3042                                        lindex[mfit_local++] = i;
3043                        }
3044                }
3045
3046                if (initialised) {
3047                        GCI_free_global_matrix(&alpha);   GCI_free_global_matrix(&covar);
3048                        GCI_ecf_free_matrix(paramtry);        GCI_free_global_vector(&beta);
3049                        GCI_free_global_vector(&dparam);  free(ochisq_trans);
3050                        initialised = 0;
3051                }
3052
3053                if (GCI_alloc_global_matrix(&alpha, mfit_global, mfit_local, ntrans)
3054                        != 0)
3055                        return -1;
3056
3057                if (GCI_alloc_global_matrix(&covar, mfit_global, mfit_local, ntrans)
3058                        != 0) {
3059                        GCI_free_global_matrix(&alpha);
3060                        return -1;
3061                }
3062
3063                if ((paramtry = GCI_ecf_matrix(ntrans, nparam)) == NULL) {
3064                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3065                        return -1;
3066                }
3067
3068                if (GCI_alloc_global_vector(&beta, mfit_global, mfit_local, ntrans)
3069                        != 0) {
3070                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3071                        GCI_ecf_free_matrix(paramtry);
3072                        return -1;
3073                }
3074
3075                if (GCI_alloc_global_vector(&dparam, mfit_global, mfit_local, ntrans)
3076                        != 0) {
3077                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3078                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
3079                        return -1;
3080                }
3081
3082                if ((ochisq_trans = (float *) malloc(ntrans * sizeof(float)))
3083                        == NULL) {
3084                        GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3085                        GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
3086                        GCI_free_global_vector(&dparam);
3087                        return -1;
3088                }
3089
3090                initialised = 1;
3091
3092                if (GCI_marquardt_global_compute_global_generic_fn(
3093                                        xincr, trans, ndata, ntrans,
3094                                        fit_start, fit_end, instr, ninstr, noise, sig,
3095                                        param, paramfree, nparam, gparam,
3096                                        mfit_global, mfit_local, gindex, lindex, fitfunc,
3097                                        yfit, dy, alpha, beta, alpha_scratch,
3098                                        chisq_trans, chisq_global, *alambda,
3099                                        &fnvals, &dy_dparam_pure, &dy_dparam_conv,
3100                                        &fnvals_len, &dy_dparam_nparam_size) != 0)
3101                        return -2;
3102
3103                *alambda = 0.001;
3104                ochisq_global = *chisq_global;
3105                for (i=0; i<ntrans; i++)
3106                        ochisq_trans[i] = chisq_trans[i];
3107
3108                /* Initialise paramtry to param */
3109                for (i=0; i<ntrans; i++) {
3110                        for (j=0; j<nparam; j++)
3111                                paramtry[i][j] = param[i][j];
3112                }
3113        }
3114
3115        /* Once converged, evaluate covariance matrix */
3116        if (*alambda == 0) {
3117                if (GCI_marquardt_global_compute_global_generic_fn_final(
3118                                        xincr, trans, ndata, ntrans,
3119                                        fit_start, fit_end, instr, ninstr, noise, sig,
3120                                        param, paramfree, nparam, gparam,
3121                                        mfit_global, mfit_local, gindex, lindex, fitfunc,
3122                                        yfit, dy, chisq_trans, chisq_global,
3123                                        &fnvals, &dy_dparam_pure, &dy_dparam_conv,
3124                                        &fnvals_len, &dy_dparam_nparam_size) != 0)
3125                        return -3;
3126                /* Don't need to do this here; if we wished to, we'd have to
3127                   move this code (the "if (*alambda == 0)" block) to after
3128                   the Gauss-Jordan call.  We'd also need to rewrite it for
3129                   our situation.... */ 
3130                //      if (mfit < nparam) {  /* no need to do this otherwise */
3131                //              GCI_covar_sort(covar, nparam, paramfree, mfit);
3132                //              GCI_covar_sort(alpha, nparam, paramfree, mfit);
3133                //      }
3134                GCI_free_global_matrix(&alpha);  GCI_free_global_matrix(&covar);
3135                GCI_ecf_free_matrix(paramtry);       GCI_free_global_vector(&beta);
3136                GCI_free_global_vector(&dparam); free(ochisq_trans);
3137                initialised = 0;
3138                return 0;
3139        }
3140
3141        /* Alter linearised fitting matrix by augmenting diagonal
3142           elements. */
3143        GCI_copy_global_matrix(covar, alpha, mfit_global, mfit_local, ntrans);
3144        GCI_copy_global_vector(dparam, beta, mfit_global, mfit_local, ntrans);
3145        for (j=0; j<mfit_global; j++)
3146                covar.P[j][j] *= 1.0 + (*alambda);
3147        for (i=0; i<ntrans; i++)
3148                for (j=0; j<mfit_local; j++)
3149                        covar.S[i][j][j] *= 1.0 + (*alambda);
3150
3151        /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */
3152        if (GCI_marquardt_global_solve_eqn(covar, dparam,
3153                                                                           mfit_global, mfit_local, ntrans) != 0)
3154                return -3;
3155
3156        /* Did the trial succeed?  Modify param by dparam... */
3157        for (i=0; i<ntrans; i++) {
3158                for (j=0; j<mfit_global; j++)
3159                        paramtry[i][gindex[j]] = param[i][gindex[j]] + dparam.global[j];
3160                for (j=0; j<mfit_local; j++)
3161                        paramtry[i][lindex[j]] =
3162                                param[i][lindex[j]] + dparam.local[i*mfit_local + j];
3163        }
3164
3165        for (i=0; i<ntrans; i++) {
3166                if (restrain == ECF_RESTRAIN_DEFAULT)
3167                        ret = check_ecf_params (paramtry[i], nparam, fitfunc);
3168                else
3169                        ret = check_ecf_user_params (paramtry[i], nparam, fitfunc);
3170
3171                if (ret != 0) {
3172                        /* Bad parameters, increase alambda and return */
3173                        *alambda *= 10.0;
3174                        return 0;
3175                }
3176        }
3177
3178        if (GCI_marquardt_global_compute_global_generic_fn(
3179                                xincr, trans, ndata, ntrans,
3180                                fit_start, fit_end, instr, ninstr, noise, sig,
3181                                paramtry, paramfree, nparam, gparam,
3182                                mfit_global, mfit_local, gindex, lindex, fitfunc,
3183                                yfit, dy, covar, dparam, alpha_scratch,
3184                                chisq_trans, chisq_global, *alambda,
3185                                &fnvals, &dy_dparam_pure, &dy_dparam_conv,
3186                                &fnvals_len, &dy_dparam_nparam_size) != 0)
3187                return -2;
3188
3189        /* Success, accept the new solution */
3190        if (*chisq_global < ochisq_global) {
3191                *alambda *= 0.1;
3192                ochisq_global = *chisq_global;
3193                for (i=0; i<ntrans; i++)
3194                        ochisq_trans[i] = chisq_trans[i];
3195                GCI_copy_global_matrix(alpha, covar, mfit_global, mfit_local, ntrans);
3196                GCI_copy_global_vector(beta, dparam, mfit_global, mfit_local, ntrans);
3197                for (i=0; i<ntrans; i++) {
3198                        for (j=0; j<nparam; j++)
3199                                param[i][j] = paramtry[i][j];
3200                }
3201        } else { /* Failure, increase alambda and return */
3202                *alambda *= 10.0;
3203                *chisq_global = ochisq_global;
3204                for (i=0; i<ntrans; i++)
3205                        chisq_trans[i] = ochisq_trans[i];
3206        }
3207
3208        return 0;
3209}
3210
3211
3212/* Here we use alpha only for scratch space */
3213int GCI_marquardt_global_compute_global_generic_fn(
3214                float xincr, float **trans, int ndata, int ntrans,
3215                int fit_start, int fit_end, float instr[], int ninstr,
3216                noise_type noise, float sig[],
3217                float **param, int paramfree[], int nparam, int gparam[],
3218                int mfit_global, int mfit_local, int gindex[], int lindex[],
3219                void (*fitfunc)(float, float [], float *, float [], int),
3220                float **yfit, float **dy, global_matrix alpha, global_vector beta,
3221                float **alpha_scratch, float *chisq_trans, float *chisq_global,
3222                float alambda, 
3223                float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
3224                int *pfnvals_len, int *pdy_dparam_nparam_size)
3225{
3226        int i, j, k, ret;
3227        float beta_scratch[MAXFIT];  /* scratch space */
3228
3229        /* We initialise P and beta_global to zero; the others don't
3230           matter, as they will be totally overwritten */
3231        for (i=0; i<mfit_global; i++) {
3232                for (j=0; j<mfit_global; j++)
3233                        alpha.P[i][j] = 0;
3234                beta.global[i] = 0;
3235        }
3236        *chisq_global = 0.0;
3237
3238        for (i=0; i<ntrans; i++) {
3239                /* Only pass the true alambda, used for initialisation, for
3240                   the first transient */
3241                ret = GCI_marquardt_compute_fn_instr(
3242                                        xincr, trans[i], ndata, fit_start, fit_end,
3243                                        instr, ninstr, noise, sig,
3244                                        param[i], paramfree, nparam, fitfunc,
3245////                                    yfit[i], dy[i], alpha_scratch, beta_scratch,
3246                                        yfit[0], dy[0], alpha_scratch, beta_scratch,
3247                                        &chisq_trans[i], (i == 0) ? alambda : 0.0,
3248                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
3249                                        pfnvals_len, pdy_dparam_nparam_size);
3250
3251                if (ret != 0) {
3252                        dbgprintf(1, "In compute_global_generic_fn, "
3253                                          "compute_fn_instr returned %d for transient %d\n",
3254                                          ret, i);
3255                        return -2;
3256                }
3257
3258                /* So now have to populate alpha and beta with the contents of
3259                   alpha_scratch and beta_scratch. */
3260
3261                for (j=0; j<mfit_global; j++) {
3262                        for (k=0; k<mfit_global; k++)
3263                                alpha.P[j][k] += alpha_scratch[gindex[j]][gindex[k]];
3264                        for (k=0; k<mfit_local; k++)
3265                                alpha.Q[j][i*mfit_local + k] =
3266                                        alpha_scratch[gindex[j]][lindex[k]];
3267                        beta.global[j] += beta_scratch[gindex[j]];
3268                }
3269                for (j=0; j<mfit_local; j++) {
3270                        for (k=0; k<mfit_local; k++)
3271                                alpha.S[i][j][k] = alpha_scratch[lindex[j]][lindex[k]];
3272                        beta.local[i*mfit_local + j] = beta_scratch[lindex[j]];
3273                }
3274
3275                *chisq_global += chisq_trans[i];
3276        }
3277
3278        return 0;
3279}
3280
3281
3282/* And the final variant */
3283int GCI_marquardt_global_compute_global_generic_fn_final(
3284                float xincr, float **trans, int ndata, int ntrans,
3285                int fit_start, int fit_end, float instr[], int ninstr,
3286                noise_type noise, float sig[],
3287                float **param, int paramfree[], int nparam, int gparam[],
3288                int mfit_global, int mfit_local, int gindex[], int lindex[],
3289                void (*fitfunc)(float, float [], float *, float [], int),
3290                float **yfit, float **dy,
3291                float *chisq_trans, float *chisq_global,       
3292                float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv,
3293                int *pfnvals_len, int *pdy_dparam_nparam_size)
3294{
3295        int i, ret;
3296
3297        *chisq_global = 0.0;
3298
3299        for (i=0; i<ntrans; i++) {
3300                /* Only pass the true alambda, used for initialisation, for
3301                   the first transient */
3302                ret = GCI_marquardt_compute_fn_final_instr(
3303                                        xincr, trans[i], ndata, fit_start, fit_end,
3304                                        instr, ninstr, noise, sig,
3305                                        param[i], paramfree, nparam, fitfunc,
3306//                                      yfit[i], dy[i], &chisq_trans[i]);
3307                                        yfit[0], dy[0], &chisq_trans[i],
3308                                        pfnvals, pdy_dparam_pure, pdy_dparam_conv,
3309                                        pfnvals_len, pdy_dparam_nparam_size);
3310
3311                if (ret != 0) {
3312                        dbgprintf(1, "In compute_global_generic_fn_final, "
3313                                          "compute_fn_final_instr returned %d for transient %d\n",
3314                                          ret, i);
3315                        return -2;
3316                }
3317
3318                *chisq_global += chisq_trans[i];
3319        }
3320
3321        return 0;
3322}
3323
3324
3325// Emacs settings:
3326// Local variables:
3327// mode: c
3328// c-basic-offset: 4
3329// tab-width: 4
3330// End:
Note: See TracBrowser for help on using the repository browser.