Context Navigation

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

Revision 7611, 99.7 KB checked in by aivar, 9 years ago (diff)
Added float chisq_delta parameter to GCI_marquart_global_generic_global_step.

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

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

Download in other formats: