source: trunk/projects/slim-curve/src/main/c/EcfSPA.c @ 7781

Revision 7781, 23.3 KB checked in by paulbarber, 8 years ago (diff)

Added GPL license to all source code.

Line 
1/*
2This file is part of the SLIM-curve package for exponential curve fitting of spectral lifetime data.
3
4Copyright (c) 2010, 2011, Gray Institute University of Oxford & UW-Madison LOCI.
5
6    This program is free software: you can redistribute it and/or modify
7    it under the terms of the GNU General Public License as published by
8    the Free Software Foundation, either version 3 of the License, or
9    (at your option) any later version.
10
11    This program is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14    GNU General Public License for more details.
15
16    You should have received a copy of the GNU General Public License
17    along with this program.  If not, see <http://www.gnu.org/licenses/>.
18 */
19
20/* The 2003 version of the ECF library.  This takes account of the fact that we may be
21   handling Poisson noise.
22
23   This file contains functions for support plane analysis.
24*/
25
26#include <stdio.h>
27#include <stdlib.h>
28#include "EcfInternal.h"
29
30
31/* The support plane analysis code calls the appropriate Marquardt
32   function to perform a fit for each parameter value set we are
33   interested in.
34
35   The variants of these functions are as follows:
36
37   GCI_SPA_1D_*  performs a 1-dimensional support plane analysis.  It
38                 takes the same fitting paramters as the corresponding
39                 standard Marquardt function, without a few of the
40                 transient-specific parameters, and with the following
41                 extra parameters:
42
43                 int spa_param      which parameter are we analysing?
44                                 int spa_nvalues    how many different parameter
45                                    values will we calculate
46                                    chi-squared values for?
47                                 float spa_low      the lowest parameter value
48                 float spa_high     the highest parameter value
49                                 float chisq[]      the resulting chisq values
50
51                                 In the case of the global analysis for exponentials
52                                 function, there is an extra parameter, int df[],
53                                 which gives the number of degrees of freedom for each
54                                 value; this is relevant if the drop_bad_transients
55                                 parameter is set.  If an individual fit fails, the
56                                 corresponding chisq[] value is set to -1.  Also, the
57                                 chisq[] array runs from chisq[0], corresponding to
58                                 spa_low to chisq[spa_nvalues-1], corresponding to
59                                 spa_high.
60
61   GCI_SPA_2D_*  performs a 2-dimensional support plane analysis.  It
62                 is the same as the 1-dimensional variant, except that
63                 now the extra parameters are:
64
65                 int spa_param1, int spa_nvalues1,
66                                 float spa_low1, float spa_high1,
67                 int spa_param1, int spa_nvalues1,
68                                 float spa_low1, float spa_high1,
69                                 float chisq[][], (int df[][])
70
71                                 Here chisq[i][j] will correspond to the i-th value of
72                                 param1 and the j-th value of param2; chisq[][] itself
73                                 will have been allocated by the GCI_ecf_matrix function.
74
75   In all cases where the original function must be provided with
76   reasonable initial estimates for the parameters, so must these
77   functions.
78*/
79
80
81int GCI_SPA_1D_marquardt(
82                                float x[], float y[], int ndata,
83                                noise_type noise, float sig[],
84                                float param[], int paramfree[], int nparam,
85                                restrain_type restrain, float chisq_delta,
86                                void (*fitfunc)(float, float [], float *, float [], int),
87                                int spa_param, int spa_nvalues,
88                                float spa_low, float spa_high,
89                                float chisq[], void (*progressfunc)(float))
90{
91        int i, j, ret;
92        float *fitted, *residuals, **covar, **alpha;
93        float param_copy[MAXFIT];
94        int paramfree_copy[MAXFIT];
95       
96        if (spa_param < 0 || spa_param >= nparam)  /* and so nparam > 0, too */
97                return -1;
98        if (spa_nvalues < 2)
99                return -2;
100       
101        if ((fitted = (float *) malloc(ndata * sizeof(float))) == NULL)
102                return -3;
103        if ((residuals = (float *) malloc(ndata * sizeof(float))) == NULL) {
104                free(fitted);
105                return -3;
106        }
107        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL) {
108                free(fitted);
109                free(residuals);
110                return -3;
111        }
112        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
113                free(fitted);
114                free(residuals);
115                GCI_ecf_free_matrix(covar);
116                return -3;
117        }
118
119        for (j=0; j<nparam; j++)
120                paramfree_copy[j] = paramfree[j];
121        paramfree_copy[spa_param] = 0;  /* we fix the parameter we are
122                                                                           analysing */
123
124        for (i=0; i<spa_nvalues; i++) {
125                /* Initialise parameter array each time */
126                for (j=0; j<nparam; j++)
127                        param_copy[j] = param[j];
128                /* Set the parameter we are analysing */
129                param_copy[spa_param] =
130                        spa_low + (spa_high - spa_low) * i / (spa_nvalues - 1);
131
132                ret = GCI_marquardt(x, y, ndata, noise, sig,
133                                                        param_copy, paramfree_copy, nparam, restrain,
134                                                        fitfunc, fitted, residuals, covar, alpha,
135                                                        &chisq[i], chisq_delta, 0, NULL);
136                       
137                progressfunc ((float)i/(float)(spa_nvalues-1));
138
139                if (ret < 0)
140                        chisq[i] = -1;
141        }
142
143        free(fitted);
144        free(residuals);
145        GCI_ecf_free_matrix(covar);
146        GCI_ecf_free_matrix(alpha);
147
148        return 0;
149}
150
151
152int GCI_SPA_2D_marquardt(
153                                float x[], float y[], int ndata,
154                                noise_type noise, float sig[],
155                                float param[], int paramfree[], int nparam,
156                                restrain_type restrain, float chisq_delta,
157                                void (*fitfunc)(float, float [], float *, float [], int),
158                                int spa_param1, int spa_nvalues1,
159                                float spa_low1, float spa_high1,
160                                int spa_param2, int spa_nvalues2,
161                                float spa_low2, float spa_high2,
162                                float **chisq, void (*progressfunc)(float))
163{
164        int i1, i2, j, ret, progress, total;
165        float *fitted, *residuals, **covar, **alpha;
166        float param_copy[MAXFIT];
167        int paramfree_copy[MAXFIT];
168       
169        if (spa_param1 < 0 || spa_param1 >= nparam)  /* and so nparam > 0, too */
170                return -1;
171        if (spa_param2 < 0 || spa_param2 >= nparam)
172                return -1;
173        if (spa_param1 == spa_param2)
174                return -1;
175        if (spa_nvalues1 < 2 || spa_nvalues2 < 2)
176                return -2;
177       
178        if ((fitted = (float *) malloc(ndata * sizeof(float))) == NULL)
179                return -3;
180        if ((residuals = (float *) malloc(ndata * sizeof(float))) == NULL) {
181                free(fitted);
182                return -3;
183        }
184        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL) {
185                free(fitted);
186                free(residuals);
187                return -3;
188        }
189        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
190                free(fitted);
191                free(residuals);
192                GCI_ecf_free_matrix(covar);
193                return -3;
194        }
195
196        for (j=0; j<nparam; j++)
197                paramfree_copy[j] = paramfree[j];
198        paramfree_copy[spa_param1] = 0;
199        paramfree_copy[spa_param2] = 0;  /* we fix the parameters we are
200                                                                                analysing */
201
202        progress = 0;
203        total    = spa_nvalues1*spa_nvalues2; 
204        for (i1=0; i1<spa_nvalues1; i1++) {
205                for (i2=0; i2<spa_nvalues2; i2++) {
206                        /* Initialise parameter array each time */
207                        for (j=0; j<nparam; j++)
208                                param_copy[j] = param[j];
209                        /* Set the parameters we are analysing */
210                        param_copy[spa_param1] =
211                                spa_low1 + (spa_high1 - spa_low1) * i1 / (spa_nvalues1 - 1);
212                        param_copy[spa_param2] =
213                                spa_low2 + (spa_high2 - spa_low2) * i2 / (spa_nvalues2 - 1);
214
215                        ret = GCI_marquardt(x, y, ndata, noise, sig,
216                                                                param_copy, paramfree_copy, nparam, restrain,
217                                                                fitfunc, fitted, residuals, covar, alpha,
218                                                                &chisq[i1][i2], chisq_delta, 0, NULL);
219                       
220                        progress++;
221                        progressfunc ((float)progress/((float)(total-1)));
222               
223                        if (ret < 0)
224                                chisq[i1][i2] = -1;
225                }
226        }
227
228        free(fitted);
229        free(residuals);
230        GCI_ecf_free_matrix(covar);
231        GCI_ecf_free_matrix(alpha);
232
233        return 0;
234}
235
236
237int GCI_SPA_1D_marquardt_instr(
238                                float xincr, float y[],
239                                int ndata, int fit_start, int fit_end,
240                                float instr[], int ninstr,
241                                noise_type noise, float sig[],
242                                float param[], int paramfree[], int nparam,
243                                restrain_type restrain, float chisq_delta,
244                                void (*fitfunc)(float, float [], float *, float [], int),
245                                int spa_param, int spa_nvalues,
246                                float spa_low, float spa_high,
247                                float chisq[], float chisq_target, void (*progressfunc)(float))
248{
249        int i, j, ret;
250        float *fitted, *residuals, **covar, **alpha;
251        float param_copy[MAXFIT];
252        int paramfree_copy[MAXFIT];
253       
254        if (spa_param < 0 || spa_param >= nparam)  /* and so nparam > 0, too */
255                return -1;
256        if (spa_nvalues < 2)
257                return -2;
258       
259        if ((fitted = (float *) malloc(ndata * sizeof(float))) == NULL)
260                return -3;
261        if ((residuals = (float *) malloc(ndata * sizeof(float))) == NULL) {
262                free(fitted);
263                return -3;
264        }
265        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL) {
266                free(fitted);
267                free(residuals);
268                return -3;
269        }
270        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
271                free(fitted);
272                free(residuals);
273                GCI_ecf_free_matrix(covar);
274                return -3;
275        }
276
277        for (j=0; j<nparam; j++)
278                paramfree_copy[j] = paramfree[j];
279        paramfree_copy[spa_param] = 0;  /* we fix the parameter we are
280                                                                           analysing */
281
282        for (i=0; i<spa_nvalues; i++) {
283                /* Initialise parameter array each time */
284                for (j=0; j<nparam; j++)
285                        param_copy[j] = param[j];
286                /* Set the parameter we are analysing */
287                param_copy[spa_param] =
288                        spa_low + (spa_high - spa_low) * i / (spa_nvalues - 1);
289
290                ret = GCI_marquardt_fitting_engine(xincr, y, ndata, fit_start, fit_end, instr, ninstr,
291                                                                                        noise, NULL, 
292                                                                                        param_copy, paramfree_copy, nparam, restrain,
293                                                                                        fitfunc, fitted, residuals, &chisq[i],
294                                                                                        covar, alpha, NULL, chisq_target, chisq_delta, 0);
295                       
296                progressfunc ((float)i/(float)(spa_nvalues-1));
297               
298                if (ret < 0)
299                        chisq[i] = -1;
300        }
301
302        free(fitted);
303        free(residuals);
304        GCI_ecf_free_matrix(covar);
305        GCI_ecf_free_matrix(alpha);
306        GCI_marquardt_cleanup();
307
308        return 0;
309}
310
311
312int GCI_SPA_2D_marquardt_instr(
313                                float xincr, float y[],
314                                int ndata, int fit_start, int fit_end,
315                                float instr[], int ninstr,
316                                noise_type noise, float sig[],
317                                float param[], int paramfree[], int nparam,
318                                restrain_type restrain, float chisq_delta,
319                                void (*fitfunc)(float, float [], float *, float [], int),
320                                int spa_param1, int spa_nvalues1,
321                                float spa_low1, float spa_high1,
322                                int spa_param2, int spa_nvalues2,
323                                float spa_low2, float spa_high2,
324                                float **chisq, float chisq_target, void (*progressfunc)(float))
325{
326        int i1, i2, j, ret, progress, total;
327        float *fitted, *residuals, **covar, **alpha;
328        float param_copy[MAXFIT];
329        int paramfree_copy[MAXFIT];
330       
331        if (spa_param1 < 0 || spa_param1 >= nparam)  /* and so nparam > 0, too */
332                return -1;
333        if (spa_param2 < 0 || spa_param2 >= nparam)
334                return -1;
335        if (spa_param1 == spa_param2)
336                return -1;
337        if (spa_nvalues1 < 2 || spa_nvalues2 < 2)
338                return -2;
339       
340        if ((fitted = (float *) malloc(ndata * sizeof(float))) == NULL)
341                return -3;
342        if ((residuals = (float *) malloc(ndata * sizeof(float))) == NULL) {
343                free(fitted);
344                return -3;
345        }
346        if ((covar = GCI_ecf_matrix(nparam, nparam)) == NULL) {
347                free(fitted);
348                free(residuals);
349                return -3;
350        }
351        if ((alpha = GCI_ecf_matrix(nparam, nparam)) == NULL) {
352                free(fitted);
353                free(residuals);
354                GCI_ecf_free_matrix(covar);
355                return -3;
356        }
357
358        for (j=0; j<nparam; j++)
359                paramfree_copy[j] = paramfree[j];
360        paramfree_copy[spa_param1] = 0;
361        paramfree_copy[spa_param2] = 0;  /* we fix the parameters we are
362                                                                                analysing */
363        progress = 0;
364        total    = spa_nvalues1*spa_nvalues2; 
365        for (i1=0; i1<spa_nvalues1; i1++) {
366                for (i2=0; i2<spa_nvalues2; i2++) {
367                        /* Initialise parameter array each time */
368                        for (j=0; j<nparam; j++)
369                                param_copy[j] = param[j];
370                        /* Set the parameters we are analysing */
371                        param_copy[spa_param1] =
372                                spa_low1 + (spa_high1 - spa_low1) * i1 / (spa_nvalues1 - 1);
373                        param_copy[spa_param2] =
374                                spa_low2 + (spa_high2 - spa_low2) * i2 / (spa_nvalues2 - 1);
375
376                        ret = GCI_marquardt_fitting_engine(xincr, y, ndata, fit_start, fit_end, instr, ninstr,
377                                                                        noise, NULL,
378                                                                        param_copy, paramfree_copy, nparam, restrain,
379                                                                        fitfunc, fitted, residuals, &chisq[i1][i2],
380                                                                        covar, alpha, NULL, chisq_target, chisq_delta, 0);
381                       
382                        progress++;
383                        progressfunc ((float)progress/((float)(total-1)));
384               
385                        if (ret < 0)
386                                chisq[i1][i2] = -1;
387       
388               
389                }
390        }
391
392        free(fitted);
393        free(residuals);
394        GCI_ecf_free_matrix(covar);
395        GCI_ecf_free_matrix(alpha);
396        GCI_marquardt_cleanup();
397
398        return 0;
399}
400
401
402int GCI_SPA_1D_marquardt_global_exps_instr(
403                                        float xincr, float **trans,
404                                        int ndata, int ntrans, int fit_start, int fit_end,
405                                        float instr[], int ninstr,
406                                        noise_type noise, float sig[], int ftype,
407                                        float **param, int paramfree[], int nparam,
408                                        restrain_type restrain, float chisq_delta, int drop_bad_transients,
409                                        int spa_param, int spa_nvalues,
410                                        float spa_low, float spa_high,
411                                        float chisq_global[], int df[], void (*progressfunc)(float))
412{
413        int i, j, k, ret;
414        float **param_copy;
415        float **fitted, **residuals, *chisq_trans;
416        int paramfree_copy[MAXFIT];
417       
418        if (spa_param < 0 || spa_param >= nparam)  /* and so nparam > 0, too */
419                return -1;
420        if (spa_nvalues < 2)
421                return -2;
422        if (ntrans < 1)
423                return -3;
424
425        if ((param_copy = GCI_ecf_matrix(ntrans, nparam)) == NULL)
426                return -4;
427        if ((fitted = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
428                GCI_ecf_free_matrix(param_copy);
429                return -4;
430        }
431        if ((residuals = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
432                GCI_ecf_free_matrix(param_copy);
433                GCI_ecf_free_matrix(fitted);
434                return -4;
435        }
436        if ((chisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL) {
437                GCI_ecf_free_matrix(param_copy);
438                GCI_ecf_free_matrix(fitted);
439                GCI_ecf_free_matrix(residuals);
440                return -4;
441        }
442
443        /* We set up the paramfree array, and also count the number of
444           free parameters for the degrees of freedom calculation */
445        for (j=0; j<nparam; j++)
446                paramfree_copy[j] = paramfree[j];
447        paramfree_copy[spa_param] = 0;  /* we fix the parameter we are
448                                                                           analysing */
449
450        /* We only need to initialise this parameter array once at the
451           beginning, in case there are any global variables. */
452        for (j=0; j<ntrans; j++) {
453                for (k=0; k<nparam; k++)
454                        param_copy[j][k] = param[j][k];
455        }
456
457        for (i=0; i<spa_nvalues; i++) {
458                /* Set the parameter we are analysing */
459                for (j=0; j<ntrans; j++) {
460                        param_copy[j][spa_param] =
461                                spa_low + (spa_high - spa_low) * i / (spa_nvalues - 1);
462                }
463
464                ret = GCI_marquardt_global_exps_instr(
465                                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
466                                                        instr, ninstr, noise, sig, ftype,
467                                                        param_copy, paramfree_copy, nparam, restrain, chisq_delta,
468                                                        fitted, residuals, chisq_trans,
469                                                        &chisq_global[i], &df[i], drop_bad_transients);
470                       
471                progressfunc ((float)i/(float)(spa_nvalues-1));
472
473                if (ret < 0)
474                        chisq_global[i] = -1;
475        }
476
477        GCI_ecf_free_matrix(param_copy);
478        GCI_ecf_free_matrix(fitted);
479        GCI_ecf_free_matrix(residuals);
480        free(chisq_trans);
481        GCI_marquardt_cleanup();
482
483        return 0;
484}
485
486
487int GCI_SPA_2D_marquardt_global_exps_instr(
488                                        float xincr, float **trans,
489                                        int ndata, int ntrans, int fit_start, int fit_end,
490                                        float instr[], int ninstr,
491                                        noise_type noise, float sig[], int ftype,
492                                        float **param, int paramfree[], int nparam,
493                                        restrain_type restrain, float chisq_delta, int drop_bad_transients,
494                                        int spa_param1, int spa_nvalues1,
495                                        float spa_low1, float spa_high1,
496                                        int spa_param2, int spa_nvalues2,
497                                        float spa_low2, float spa_high2,
498                                        float **chisq_global, int **df, void (*progressfunc)(float))
499{
500        int i1, i2, j, k, ret, progress, total;
501        float **param_copy;
502        float **fitted, **residuals, *chisq_trans;
503        int paramfree_copy[MAXFIT];
504       
505        if (spa_param1 < 0 || spa_param1 >= nparam)  /* and so nparam > 0, too */
506                return -1;
507        if (spa_param2 < 0 || spa_param2 >= nparam)
508                return -1;
509        if (spa_param1 == spa_param2)
510                return -1;
511        if (spa_nvalues1 < 2 || spa_nvalues2 < 2)
512                return -2;
513        if (ntrans < 1)
514                return -3;
515
516        if ((param_copy = GCI_ecf_matrix(ntrans, nparam)) == NULL)
517                return -4;
518        if ((fitted = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
519                GCI_ecf_free_matrix(param_copy);
520                return -4;
521        }
522        if ((residuals = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
523                GCI_ecf_free_matrix(param_copy);
524                GCI_ecf_free_matrix(fitted);
525                return -4;
526        }
527        if ((chisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL) {
528                GCI_ecf_free_matrix(param_copy);
529                GCI_ecf_free_matrix(fitted);
530                GCI_ecf_free_matrix(residuals);
531                return -4;
532        }
533
534        /* We set up the paramfree array, and also count the number of
535           free parameters for the degrees of freedom calculation */
536        for (j=0; j<nparam; j++)
537                paramfree_copy[j] = paramfree[j];
538        paramfree_copy[spa_param1] = 0;
539        paramfree_copy[spa_param2] = 0;  /* we fix the parameters we are
540                                                                                analysing */
541
542        /* We only need to initialise this parameter array once at the
543           beginning, in case there are any global variables; all of the
544           other estimates are done by the ECF code itself */
545        for (j=0; j<ntrans; j++) {
546                for (k=0; k<nparam; k++)
547                        param_copy[j][k] = param[j][k];
548        }
549
550        progress = 0;
551        total    = spa_nvalues1*spa_nvalues2; 
552        for (i1=0; i1<spa_nvalues1; i1++) {
553                for (i2=0; i2<spa_nvalues2; i2++) {
554                        /* Set the parameters we are analysing */
555                        for (j=0; j<ntrans; j++) {
556                                param_copy[j][spa_param1] = spa_low1 +
557                                        (spa_high1 - spa_low1) * i1 / (spa_nvalues1 - 1);
558                                param_copy[j][spa_param2] = spa_low2 +
559                                        (spa_high2 - spa_low2) * i2 / (spa_nvalues2 - 1);
560                        }
561
562                        ret = GCI_marquardt_global_exps_instr(
563                                                        xincr, trans, ndata, ntrans, fit_start, fit_end,
564                                                        instr, ninstr, noise, sig, ftype,
565                                                        param_copy, paramfree_copy, nparam, restrain, chisq_delta,
566                                                        fitted, residuals, chisq_trans,
567                                                        &chisq_global[i1][i2], &df[i1][i2],
568                                                        drop_bad_transients);
569                       
570                        progress++;
571                        progressfunc ((float)progress/((float)(total-1)));
572               
573                        if (ret < 0)
574                                chisq_global[i1][i2] = -1;
575                }
576        }
577
578        GCI_ecf_free_matrix(param_copy);
579        GCI_ecf_free_matrix(fitted);
580        GCI_ecf_free_matrix(residuals);
581        free(chisq_trans);
582        GCI_marquardt_cleanup();
583
584        return 0;
585}
586
587
588int GCI_SPA_1D_marquardt_global_generic_instr(
589                                        float xincr, float **trans,
590                                        int ndata, int ntrans, int fit_start, int fit_end,
591                                        float instr[], int ninstr,
592                                        noise_type noise, float sig[],
593                                        float **param, int paramfree[], int nparam, int gparam[],
594                                        restrain_type restrain, float chisq_delta,
595                                        void (*fitfunc)(float, float [], float *, float [], int),
596                                        int spa_param, int spa_nvalues,
597                                        float spa_low, float spa_high,
598                                        float chisq_global[], int df[], void (*progressfunc)(float))
599{
600        int i, j, k, ret;
601        float **param_copy;
602        float **fitted, **residuals, *chisq_trans;
603        int paramfree_copy[MAXFIT];
604       
605        if (spa_param < 0 || spa_param >= nparam)  /* and so nparam > 0, too */
606                return -1;
607        if (spa_nvalues < 2)
608                return -2;
609        if (ntrans < 1)
610                return -3;
611
612        if ((param_copy = GCI_ecf_matrix(ntrans, nparam)) == NULL)
613                return -4;
614        if ((fitted = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
615                GCI_ecf_free_matrix(param_copy);
616                return -4;
617        }
618        if ((residuals = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
619                GCI_ecf_free_matrix(param_copy);
620                GCI_ecf_free_matrix(fitted);
621                return -4;
622        }
623        if ((chisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL) {
624                GCI_ecf_free_matrix(param_copy);
625                GCI_ecf_free_matrix(fitted);
626                GCI_ecf_free_matrix(residuals);
627                return -4;
628        }
629
630        /* We set up the paramfree array, and also count the number of
631           free parameters for the degrees of freedom calculation */
632        for (j=0; j<nparam; j++)
633                paramfree_copy[j] = paramfree[j];
634        paramfree_copy[spa_param] = 0;  /* we fix the parameter we are
635                                                                           analysing */
636
637        for (i=0; i<spa_nvalues; i++) {
638                /* Initialise parameter array each time */
639                for (j=0; j<ntrans; j++) {
640                        for (k=0; k<nparam; k++)
641                                param_copy[j][k] = param[j][k];
642                        /* Set the parameter we are analysing */
643                        param_copy[j][spa_param] =
644                                spa_low + (spa_high - spa_low) * i / (spa_nvalues - 1);
645                }
646
647                ret = GCI_marquardt_global_generic_instr(
648                                                xincr, trans, ndata, ntrans, fit_start, fit_end,
649                                                instr, ninstr, noise, sig,
650                                                param_copy, paramfree_copy, nparam, gparam, restrain, chisq_delta,
651                                                fitfunc, fitted, residuals, chisq_trans,
652                                                &chisq_global[i], &df[i]);
653                       
654                progressfunc ((float)i/(float)(spa_nvalues-1));
655
656                if (ret < 0)
657                        chisq_global[i] = -1;
658        }
659
660        GCI_ecf_free_matrix(param_copy);
661        GCI_ecf_free_matrix(fitted);
662        GCI_ecf_free_matrix(residuals);
663        free(chisq_trans);
664        GCI_marquardt_cleanup();
665
666        return 0;
667}
668
669
670int GCI_SPA_2D_marquardt_global_generic_instr(
671                                        float xincr, float **trans,
672                                        int ndata, int ntrans, int fit_start, int fit_end,
673                                        float instr[], int ninstr,
674                                        noise_type noise, float sig[],
675                                        float **param, int paramfree[], int nparam, int gparam[],
676                                        restrain_type restrain, float chisq_delta,
677                                        void (*fitfunc)(float, float [], float *, float [], int),
678                                        int spa_param1, int spa_nvalues1,
679                                        float spa_low1, float spa_high1,
680                                        int spa_param2, int spa_nvalues2,
681                                        float spa_low2, float spa_high2,
682                                        float **chisq_global, int **df, void (*progressfunc)(float))
683{
684        int i1, i2, j, k, ret, progress, total;
685        float **param_copy;
686        float **fitted, **residuals, *chisq_trans;
687        int paramfree_copy[MAXFIT];
688       
689        if (spa_param1 < 0 || spa_param1 >= nparam)  /* and so nparam > 0, too */
690                return -1;
691        if (spa_param2 < 0 || spa_param2 >= nparam)
692                return -1;
693        if (spa_param1 == spa_param2)
694                return -1;
695        if (spa_nvalues1 < 2 || spa_nvalues2 < 2)
696                return -2;
697        if (ntrans < 1)
698                return -3;
699
700        if ((param_copy = GCI_ecf_matrix(ntrans, nparam)) == NULL)
701                return -4;
702        if ((fitted = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
703                GCI_ecf_free_matrix(param_copy);
704                return -4;
705        }
706        if ((residuals = GCI_ecf_matrix(ntrans, ndata)) == NULL) {
707                GCI_ecf_free_matrix(param_copy);
708                GCI_ecf_free_matrix(fitted);
709                return -4;
710        }
711        if ((chisq_trans = (float *) malloc(ntrans * sizeof(float))) == NULL) {
712                GCI_ecf_free_matrix(param_copy);
713                GCI_ecf_free_matrix(fitted);
714                GCI_ecf_free_matrix(residuals);
715                return -4;
716        }
717
718        /* We set up the paramfree array, and also count the number of
719           free parameters for the degrees of freedom calculation */
720        for (j=0; j<nparam; j++)
721                paramfree_copy[j] = paramfree[j];
722        paramfree_copy[spa_param1] = 0;
723        paramfree_copy[spa_param2] = 0;  /* we fix the parameters we are
724                                                                                analysing */
725
726        progress = 0;
727        total    = spa_nvalues1*spa_nvalues2; 
728        for (i1=0; i1<spa_nvalues1; i1++) {
729                for (i2=0; i2<spa_nvalues2; i2++) {
730                        /* Initialise parameter array each time */
731                        for (j=0; j<ntrans; j++) {
732                                for (k=0; k<nparam; k++)
733                                        param_copy[j][k] = param[j][k];
734                                /* Set the parameters we are analysing */
735                                param_copy[j][spa_param1] = spa_low1 +
736                                        (spa_high1 - spa_low1) * i1 / (spa_nvalues1 - 1);
737                                param_copy[j][spa_param2] = spa_low2 +
738                                        (spa_high2 - spa_low2) * i2 / (spa_nvalues2 - 1);
739                        }
740
741                        ret = GCI_marquardt_global_generic_instr(
742                                                xincr, trans, ndata, ntrans, fit_start, fit_end,
743                                                instr, ninstr, noise, sig,
744                                                param_copy, paramfree_copy, nparam, gparam, restrain, chisq_delta,
745                                                fitfunc, fitted, residuals, chisq_trans,
746                                                &chisq_global[i1][i2], &df[i1][i2]);
747                       
748                        progress++;
749                        progressfunc ((float)progress/((float)(total-1)));
750               
751                        if (ret < 0)
752                                chisq_global[i1][i2] = -1;
753                }
754        }
755
756        GCI_ecf_free_matrix(param_copy);
757        GCI_ecf_free_matrix(fitted);
758        GCI_ecf_free_matrix(residuals);
759        free(chisq_trans);
760        GCI_marquardt_cleanup();
761
762        return 0;
763}
764
765
766// Emacs settings:
767// Local variables:
768// mode: c
769// c-basic-offset: 4
770// tab-width: 4
771// End:
Note: See TracBrowser for help on using the repository browser.