Context Navigation

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

Revision 7803, 100.7 KB checked in by paulbarber, 8 years ago (diff)
Fixed some memory leaks and some compile warnings.

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: