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

Revision 7606, 99.7 KB checked in by paulbarber, 9 years ago (diff)

To compile in TRI2, moved declarations to top of relevant block.
Added empty GCI_marquardt.
Added matrix array routinnes.
Added Global and SPA files.

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