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

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

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

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