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

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

To compile in TRI2, moved declarations to top of relevant block.

Line
1//#include <ansi_c.h>
2/* The 2010 version of the ECF library.  This has basically been
3   completely rewritten to avoid license issues.
4   Also, this takes account of the fact that we may be
5   handling Poisson noise.
6
7   This file contains utility code for various functions needed by the
8   ECF library, and which are not particularly specific to the single
9   or global fitting routines.
10*/
11
12#include <stdio.h>
13#include <stdlib.h>
14#include <stdarg.h>
15#include <math.h>
16#ifdef _CVI_
17#include <userint.h>
18#endif
19#include "EcfInternal.h"  /* For #defines */
20
21int ECF_debug = 0;
22
23
24/********************************************************************
25
26                           GAUSS-JORDAN AND COVAR-SORTING ROUTINES
27
28*********************************************************************/
29
30
31/* Linear equation solution by Gauss-Jordan elimination.
32   a[0..n-1][0..n-1] is the input matrix,
33   b[0..n] is input containing the single right-hand side vector.
34   On output, a is replaced by its matrix inverse and b is replaced by
35   the corresponding set of solution vectors. */
36
37#define SWAP(a,b) { temp=(a); (a)=(b); (b)=temp; }
38
39/*
40 * Based on LinearEquationSolving from linearEquations.c by Henry Guennadi Levkin.
41 */
42int GCI_gauss_jordan(float **a, int n, float *b)
43{
44    float max;
45    float tmp;
46
47    int i, j, k, m;
48
49    // base row of matrix
50    for (k = 0; k < n - 1; ++k)
51    {
52        // search for line with max element
53        max = fabs(a[k][k]);
54        m = k;
55        for (i = k + 1; i < n; ++i)
56        {
57            if (max < fabs(a[i][k])) // row i col k
58            {
59                max = a[i][k];
60                m = i;
61            }
62        }
63
64        // permutation of base line (index k) and max element line (index m)
65        if (m != k)
66        {
67            for (i = k; i < n; ++i)
68            {
69                tmp = a[k][i];
70                a[k][i] = a[m][i];
71                a[m][i] = tmp;
72            }
73            tmp = b[k];
74            b[k] = b[m];
75            b[m] = tmp;
76        }
77
78        if (0.0 == a[k][k])
79        {
80            return -2; // singular matrix
81        }
82
83        // triangulation of matrix with coefficients
84        for (j = k + 1; j < n; ++j) // current row of matrix
85        {
86            tmp = -a[j][k] / a[k][k];
87            for (i = k; i < n; ++i)
88            {
89                a[j][i] += tmp * a[k][i];
90            }
91            b[j] += tmp * b[k]; // free member recalculation
92        }
93    }
94
95    for (k = n - 1; k >= 0; --k)
96    {
97        for (i = k + 1; i < n; ++i)
98        {
99            b[k] -= a[k][i] * b[i];
100        }
101        b[k] /= a[k][k];
102    }
103    return 0;
104}
105
106//==============================================================================
107// return 1 if system not solving
108// nDim - system dimension
109// pfMatr - matrix with coefficients
110// pfVect - vector with free members
111// pfSolution - vector with system solution
112// pfMatr becames trianglular after function call
113// pfVect changes after function call
114//
116//
117//==============================================================================
118int LinearEquationsSolving(int nDim, double* pfMatr, double* pfVect, double* pfSolution)
119{
120  double fMaxElem;
121  double fAcc;
122
123  int i , j, k, m;
124
125
126  for(k=0; k<(nDim-1); k++) // base row of matrix
127  {
128    // search of line with max element
129    fMaxElem = fabs( pfMatr[k*nDim + k] );
130    m = k;
131    for(i=k+1; i<nDim; i++)
132    {
133      if(fMaxElem < fabs(pfMatr[i*nDim + k]) )
134      {
135        fMaxElem = pfMatr[i*nDim + k];
136        m = i;
137      }
138    }
139
140    // permutation of base line (index k) and max element line(index m)
141    if(m != k)
142    {
143      for(i=k; i<nDim; i++)
144      {
145        fAcc               = pfMatr[k*nDim + i];
146        pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
147        pfMatr[m*nDim + i] = fAcc;
148      }
149      fAcc = pfVect[k];
150      pfVect[k] = pfVect[m];
151      pfVect[m] = fAcc;
152    }
153
154    if( pfMatr[k*nDim + k] == 0.) return 1; // needs improvement !!!
155
156    // triangulation of matrix with coefficients
157    for(j=(k+1); j<nDim; j++) // current row of matrix
158    {
159      fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
160      for(i=k; i<nDim; i++)
161      {
162        pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
163      }
164      pfVect[j] = pfVect[j] + fAcc*pfVect[k]; // free member recalculation
165    }
166  }
167
168  for(k=(nDim-1); k>=0; k--)
169  {
170    pfSolution[k] = pfVect[k];
171    for(i=(k+1); i<nDim; i++)
172    {
173      pfSolution[k] -= (pfMatr[k*nDim + i]*pfSolution[i]);
174    }
175    pfSolution[k] = pfSolution[k] / pfMatr[k*nDim + k];
176  }
177
178  return 0;
179}
180
181
182void GCI_covar_sort(float **covar, int nparam, int paramfree[], int mfit)
183{
184        int i, j, k;
185        float temp;
186
187        for (i=mfit; i<nparam; i++)
188                for (j=0; j<=i; j++)
189                        covar[i][j] = covar[j][i] = 0.0;
190
191        k = mfit-1;
192        for (j=nparam-1; j>=0; j--)
193        {
194                if (paramfree[j])
195                {
196                        for (i=0; i<nparam; i++) SWAP(covar[i][k], covar[i][j]);
197                        for (i=0; i<nparam; i++) SWAP(covar[k][i], covar[j][i]);
198                        k--;
199                }
200        }
201}
202#undef SWAP
203
204
205/********************************************************************
206
207                                          ARRAY ALLOCATION ROUTINES
208http://www.nr.com/public-domain.html
209*********************************************************************/
210
211/* This function allocates a float matrix with subscript range
212 * m[0..nrows-1][0..ncols-1]
213 */
214float **GCI_ecf_matrix(long nrows, long ncols)
215{
216    int row_size = nrows * sizeof(float *);
217    int data_size = nrows * ncols * sizeof(float);
218    unsigned char *raw = malloc(row_size + data_size);
219    float **row = (float **) raw;
220    float *data = (float *) (row + nrows);
221    int i;
222
223        if (NULL == raw)
224    {
225        return NULL;
226    }
227
228        for (i = 0; i < nrows; ++i)
229    {
230        row[i] = data;
231        data += ncols;
232    }
233    return row;
234}
235
236/* Frees allocated float matrix.
237 */
238void GCI_ecf_free_matrix(float **m)
239{
240    if (NULL != m)
241    {
242        free(m);
243    }
244}
245
246float ***GCI_ecf_matrix_array(long nblocks, long nrows, long ncols)
247/* allocate a float matrix array with range
248   marr[0..nblocks][0..nrows][0..ncols] */
249{
250        long i;
251        float ***marr;
252
253        /* allocate pointers to blocks */
254        if ((marr = (float ***) malloc(nblocks * sizeof(float **))) == NULL)
255                return NULL;
256
257        /* allocate blocks (= pointers to rows) and set pointers to them */
258        if ((marr[0] = (float **) malloc(nblocks * nrows * sizeof(float *)))
259                == NULL) {
260                free(marr);
261                return NULL;
262        }
263
264        for (i=1; i<nblocks; i++)
265                marr[i] = marr[i-1] + nrows;
266
267        /* allocate rows (= pointers to column entries) and set pointers to them */
268        if ((marr[0][0] = (float *)malloc(nblocks * nrows * ncols * sizeof(float)))
269                == NULL) {
270                free(marr[0]);
271                free(marr);
272                return NULL;
273        }
274
275        /* This sneaky loop does the whole lot!  For note that
276           marr[block][row] = marr[0][block*nrows + row] */
277        for (i=1; i < nblocks * nrows; i++)
278                marr[0][i] = marr[0][i-1] + ncols;
279
280        return marr;
281}
282
283
284void GCI_ecf_free_matrix_array(float ***marr)
285{
286        if (marr != NULL) {
287                free(marr[0][0]);
288                free(marr[0]);
289                free(marr);
290        }
291}
292
293/********************************************************************
294
295                                FITTING FUNCTION CALCULATING FUNCTIONS
296
297*********************************************************************/
298
299
300/* Some match functions for use in the above routines
301
302   Prototypes: func(float x, float param[], float *y,
303                    float dy_dparam[], int nparam)
304   Inputs:  x            data point
305            param[]      function parameters
306            nparam       number of parameters (may or may not be
307                         used by the function)
308   Outputs: y            function value at this point
309            dy_dparam[]  vector of dy/dparam_k values at this point
310
311   There are also variants which calculate the values for a whole
312   array of x values in one go, which is significantly faster, as
313   there are many fewer calls to exp().  Note that this is only
314   reliable if tau (or lambda) is positive.  These functions are
315   called from above in certain situations to speed things up.  Their
316   prototype is:
317
318     void fitfunc_array(float xincr, float param[],
319                        float *y, float **dy_dparam, int nx, int nparam)
320
321   where the fitfunc will be evaluated for x=i*xincr, with i=0, 1,
322   ..., nx-1.  The results will be placed in the y[] array and the
323   dy_dparam[nx][nparam] array (which is assumed to have been allocated by
324   GCI_matrix).
325*/
326
327/* This one produces multiexponentials using lambdas:
328
329      y(x) = param[0] + param[1]*exp(-param[2]*x) +
330               param[3]*exp(-param[4]*x) + ...
331
332   This gives:
333
334      dy/dparam_0 = 1
335      dy/dparam_1 = exp(-param[2]*x)
336      dy/dparam_2 = -param[1]*x*exp(-param[2]*x)
337
338   and similarly for dy/dparam_3 and dy/dparam_4, etc.
339
340   As discussed above, though, we ignore the param[0] term.
341*/
342
343void GCI_multiexp_lambda(float x, float param[],
344                                                 float *y, float dy_dparam[], int nparam)
345{
346        int i;
347        float ex;
348
349        *y = 0;
350
351        for (i=1; i<nparam-1; i+=2) {
352                dy_dparam[i] = ex = exp(-param[i+1] * x);
353                ex *= param[i];
354                *y += ex;
355                dy_dparam[i+1] = -ex * x;
356        }
357}
358
359
360int multiexp_lambda_array(float xincr, float param[],
361                                                  float *y, float **dy_dparam, int nx, int nparam)
362{
363        int i, j;
364        float ex;
365        double exincr[MAXFIT];  /* exp(-lambda*xincr) */
366        double excur[MAXFIT];   /* exp(-lambda*x)     */
367
368        if (xincr <= 0) return -1;
369
370        for (j=1; j<nparam-1; j+=2) {
371                if (param[j+1] < 0) return -1;
372                excur[j] = 1.0;
373                exincr[j] = exp(-(double) param[j+1] * xincr);
374        }
375
376        for (i=0; i<nx; i++) {
377                y[i] = 0;
378                for (j=1; j<nparam-1; j+=2) {
379                        dy_dparam[i][j] = ex = excur[j];
380                        ex *= param[j];
381                        y[i] += ex;
382                        dy_dparam[i][j+1] = -ex * xincr * i;
383                        /* And ready for next loop... */
384                        excur[j] *= exincr[j];
385                }
386        }
387
388        return 0;
389}
390
391
392/* This one produces multiexponentials using taus:
393
394      y(x) = param[0] + param[1]*exp(-x/param[2]) +
395               param[3]*exp(-x/param[4]) + ...
396
397   This gives:
398
399      dy/dparam_0 = 1
400      dy/dparam_1 = exp(-x/param[2])
401      dy/dparam_2 = param[1]*x*exp(-x/param[2]) / param[2]^2
402
403   and similarly for dy/dparam_3 and dy/dparam_4, etc.
404
405   Again, we ignore the param[0] term.
406*/
407
408void GCI_multiexp_tau(float x, float param[],
409                                          float *y, float dy_dparam[], int nparam)
410{
411        int i;
412        float ex, xa;
413
414        *y = 0;
415
416        for (i=1; i<nparam-1; i+=2) {
417                xa = x / param[i+1];
418                dy_dparam[i] = ex = exp(-xa);
419                ex *= param[i];
420                *y += ex;
421                dy_dparam[i+1] = ex * xa / param[i+1];
422        }
423}
424
425
426int multiexp_tau_array(float xincr, float param[],
427                                           float *y, float **dy_dparam, int nx, int nparam)
428{
429        int i, j;
430        float ex;
431        float a2[MAXFIT];       /* 1/(param[j]*param[j]) for taus */
432        double exincr[MAXFIT];  /* exp(-xincr/tau) */
433        double excur[MAXFIT];   /* exp(-x/tau)     */
434
435        if (xincr <= 0) return -1;
436
437        for (j=1; j<nparam-1; j+=2) {
438                if (param[j+1] < 0) return -1;
439                excur[j] = 1.0;
440                exincr[j] = exp(-xincr / (double) param[j+1]);
441                a2[j] = 1 / (param[j+1] * param[j+1]);
442        }
443
444        for (i=0; i<nx; i++) {
445                y[i] = 0;
446                for (j=1; j<nparam-1; j+=2) {
447                        dy_dparam[i][j] = ex = excur[j];
448                        ex *= param[j];
449                        y[i] += ex;
450                        dy_dparam[i][j+1] = ex * xincr * i * a2[j];
451                        /* And ready for next loop... */
452                        excur[j] *= exincr[j];
453                }
454        }
455
456        return 0;
457}
458
459
460/* And this one produces stretched exponentials:
461
462      y(x) = Z + A exp(-(x/tau)^(1/h))
463
464   or translated into C-speak:
465
466      y(x) = param[0] + param[1]*exp(-(x/param[2])^(1/param[3]))
467
468   This gives:
469
470      dy/dparam_0 = 1
471      dy/dparam_1 = exp(-(x/param[2])^(1/param[3]))
472      dy/dparam_2 = param[1]*exp(-(x/param[2])^(1/param[3])) *
473                      (x/param[2])^(1/param[3]) * (1/param[2]*param[3])
474      dy/dparam_3 = param[1]*exp(-(x/param[2])^(1/param[3])) *
475                      (x/param[2])^(1/param[3]) *
476                      (1/param[3]^2)*log(x/param[2])
477*/
478
479void GCI_stretchedexp(float x, float param[],
480                                          float *y, float dy_dparam[], int nparam)
481{
482        float ex, xa, lxa, xah;
483
484        if (x > 0) {
485                xa = x / param[2];         /* xa = x/param[2] */
486                lxa = log(xa);             /* lxa = log(x/param[2]) */
487                xah = exp(lxa / param[3]); /* xah = exp(log(x/param[2])/param[3])
488                                                  = (x/param[2])^(1/param[3]) */
489                dy_dparam[1] = ex = exp(-xah);
490                                           /* ex = exp(-(x/param[2])^(1/param[3])) */
491                ex *= param[1];            /* ex = param[1] *
492                                                     exp(-(x/param[2])^(1/param[3])) */
493                *y = ex;                   /* y is now correct */
494                ex *= xah / param[3];      /* ex = param[1] * exp(...) *
495                                                      (x/param[2])^(1/param[3]) *
496                                                      1/param[3]   */
497                dy_dparam[2] = ex / param[2];
498                dy_dparam[3] = ex * lxa / param[3];
499        } else if (x > -1e-10) {  // Almost zero
500                *y = param[1];
501                dy_dparam[1] = 1;
502                dy_dparam[2] = dy_dparam[3] = 0;
503        } else {
504                /* Ouch: x<0 */
505                fprintf(stderr, "Can't have x < 0 in stretched exponential!!\n");
506        }
507}
508
509/* This is actually essentially the same as the other version; because
510   of the power (x/tau)^(1/h), we cannot use the same tricks as for
511   the multiexponential case, unfortunately. */
512
513int stretchedexp_array(float xincr, float param[],
514                                           float *y, float **dy_dparam, int nx, int nparam)
515{
516        int i;
517        float ex, lxa, xah, a2inv, a3inv;
518        double xa, xaincr;
519
520        if (xincr == 0) return -1;
521        xa=0;
522        xaincr = xincr / param[2];
523        a2inv = 1/param[2];
524        a3inv = 1/param[3];
525
526        /* When x=0 */
527        y[0] = param[1];
528        dy_dparam[0][1] = 1;
529        dy_dparam[0][2] = dy_dparam[0][3] = 0;
530
531        for (i=1; i<nx; i++) {
532                xa += xaincr;       /* xa = (xincr*i)/param[2] */
533                lxa = log(xa);      /* lxa = log(x/param[2]) */
534                xah = exp(lxa * a3inv);  /* xah = exp(log(x/param[2])/param[3])
535                                                = (x/param[2])^(1/param[3]) */
536                dy_dparam[i][1] = ex = exp(-xah);
537                                    /* ex = exp(-(x/param[2])^(1/param[3])) */
538                ex *= param[1];     /* ex = param[1]*exp(-(x/param[2])^(1/param[3])) */
539                y[i] = ex;          /* y is now correct */
540                ex *= xah * a3inv;  /* ex = param[1] * exp(...) *
541                                              (x/param[2])^(1/param[3]) * 1/param[3] */
542                dy_dparam[i][2] = ex * a2inv;
543                dy_dparam[i][3] = ex * lxa * a3inv;
544        }
545
546        return 0;
547}
548
549
550/********************************************************************
551
552                           CHECKING AND RESTRICTED FITTING ROUTINES
553
554*********************************************************************/
555
556
557/* This is the default routine */
558#define MIN_Z -1e5
559#define MIN_Z_FACTOR 0.4
560#define MAX_Z 1e10  /* Silly */
561#define MIN_A 0
562#define MAX_A 1e10  /* Silly */
563#define MIN_TAU 0.001  /* Works for lambda too */
564#define MAX_TAU 1000
565#define MIN_H 1
566#define MAX_H 10
567
568int check_ecf_params (float param[], int nparam,
569                                        void (*fitfunc)(float, float [], float *, float [], int))
570{
571        if (fitfunc == GCI_multiexp_lambda || fitfunc == GCI_multiexp_tau) {
572                switch (nparam) {
573                case 3:
574                        if (param[0] < MIN_Z ||
575                                param[0] < -MIN_Z_FACTOR * fabs(param[1]) ||
576                                param[0] > MAX_Z)
577                                return -21;
578                        if (param[1] < MIN_A || param[1] > MAX_A)
579                                return -22;
580                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
581                                return -23;
582                        break;
583
584                case 5:
585                        if (param[0] < MIN_Z ||
586                                param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3])) ||
587                                param[0] > MAX_Z)
588                                return -21;
589                        if (param[1] < MIN_A || param[1] > MAX_A)
590                                return -22;
591                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
592                                return -23;
593                        if (param[3] < MIN_A || param[3] > MAX_A)
594                                return -24;
595                        if (param[4] < MIN_TAU || param[4] > MAX_TAU)
596                                return -25;
597                        break;
598
599                case 7:
600                        if (param[0] < MIN_Z ||
601                                param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3]) +
602                                                                                        fabs(param[5])) ||
603                                param[0] > MAX_Z)
604                                return -21;
605                        if (param[1] < MIN_A || param[1] > MAX_A)
606                                return -22;
607                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
608                                return -23;
609                        if (param[3] < MIN_A || param[3] > MAX_A)
610                                return -24;
611                        if (param[4] < MIN_TAU || param[4] > MAX_TAU)
612                                return -25;
613                        if (param[5] < MIN_A || param[5] > MAX_A)
614                                return -26;
615                        if (param[6] < MIN_TAU || param[6] > MAX_TAU)
616                                return -27;
617                        break;
618                }
619        } else if (fitfunc == GCI_stretchedexp) {
620                if (param[0] < MIN_Z || param[0] < -MIN_Z_FACTOR * fabs(param[1]) ||
621                        param[0] > MAX_Z)
622                        return -21;
623                if (param[1] < MIN_A || param[1] > MAX_A)
624                        return -22;
625                if (param[2] < MIN_TAU || param[2] > MAX_TAU)
626                        return -23;
627                if (param[3] < MIN_H || param[3] > MAX_H)
628                        return -24;
629        }
630        return 0;
631}
632
633
634/* For the user-specified version, we have some global variables to
635   store the settings between invocations. */
636
637static int restrain_nparam=0;  /* How many parameters have we set up? */
638static int restrain_restraining[MAXFIT];  /* Do we check parameter i? */
639static float restrain_minval[MAXFIT]; /* Minimum acceptable parameter value */
640static float restrain_maxval[MAXFIT]; /* Maximum acceptable parameter value */
641
642int GCI_set_restrain_limits(int nparam, int restrain[],
643                                                        float minval[], float maxval[])
644{
645        int i;
646
647        if (nparam < 0 || nparam > MAXFIT)
648                return -1;
649
650        /* We're going to be doing something, so clear the memory */
651        restrain_nparam = 0;
652
653        for (i=0; i<nparam; i++) {
654                if (restrain[i]) {
655                        restrain_restraining[i] = 1;
656
657                        if (minval[i] > maxval[i])
658                                return -2;
659                        restrain_minval[i] = minval[i];
660                        restrain_maxval[i] = maxval[i];
661                } else
662                        restrain_restraining[i] = 0;
663        }
664
665        restrain_nparam = nparam;
666        return 0;
667}
668
669/* original version
670int check_ecf_user_params (float param[], int nparam,
671                                        void (*fitfunc)(float, float [], float *, float [], int))
672{
673        int i;
674
675        if (restrain_nparam != nparam) {
676                dbgprintf(0, "Using user-defined parameter restraining with "
677                                  "wrong number of parameters:\n"
678                                  "actual nparam = %d, user restraining nparam = %d\n"
679                                  "Defaulting to standard tests\n", nparam, restrain_nparam);
680                return check_ecf_params(param, nparam, fitfunc);
681        }
682
683        for (i=0; i<nparam; i++) {
684                if (restrain_restraining[i]) {
685                        if ((param[i] < restrain_minval[i]) ||
686                                (param[i] > restrain_maxval[i]))
687                                return -21;
688                }
689        }
690
691        return 0;
692}
693*/
694// new version from J Gilbey 31.03.03
695int check_ecf_user_params (float param[], int nparam,
696                                        void (*fitfunc)(float, float [], float *, float [], int))
697{
698        int i;
699
700
701        if (restrain_nparam != nparam) {
702                dbgprintf(0, "Using user-defined parameter restraining with "
703                                  "wrong number of parameters:\n"
704                                  "actual nparam = %d, user restraining nparam = %d\n"
705                                  "Defaulting to standard tests\n", nparam, restrain_nparam);
706                return check_ecf_params(param, nparam, fitfunc);
707        }
708
709
710        for (i=0; i<nparam; i++) {
711                if (restrain_restraining[i]) {
712                        if (param[i] < restrain_minval[i])
713                                param[i] = restrain_minval[i];
714                        else if (param[i] > restrain_maxval[i])
715                                param[i] = restrain_maxval[i];
716                }
717        }
718
719
720        return 0;
721}
722
723
724/********************************************************************
725
726                                                          MISCELLANY
727
728*********************************************************************/
729
730
731/* From filters */
732int ECF_Find_Float_Max (float data[], int np, float *max_val)
733{
734        int i, maxi;
735        float *ptr=data;
736
737        if (np < 1)        /* daft input */
738                return -1;
739
740        maxi = 0;
741        *max_val = *ptr;
742
743        for (i=1; i<np; i++)
744                if (*++ptr>*max_val) {
745                        *max_val = *ptr;
746                        maxi = i;
747                }
748
749        return maxi;
750}
751
752
753/* Debugging version of printf */
754#ifdef _CVI_
755int dbgprintf(int dbg_level, const char *format, ...)
756{
757        int ret = 0;
758        char msg[1000];
759        va_list va;
760
761        if (ECF_debug >= dbg_level) {
762                va_start(va, format);
763                ret = vsprintf(msg, format, va);  // Here's hoping, no vsnprintf ...
764                MessagePopup("Debug info", msg);
765                va_end(va);
766        }
767
768        return ret;
769}
770#else
771int dbgprintf(int dbg_level, const char *format, ...)
772{
773        int ret = 0;
774        va_list va;
775
776        if (ECF_debug >= dbg_level) {
777                va_start(va, format);
778                ret = vfprintf(stderr, format, va);
779                va_end(va);
780        }
781
782        return ret;
783}
784#endif
785
786#ifdef TESTCHISQ
787int main(int ac, char **av)
788{
789        int i, nu, ret;
790        int debug = 0;
791        float root, chisq;
792        float percents[4] = { 0.50, 0.68, 0.90, 0.95 };
793
794        for (i=0; i<4; i++) {
795                chisq = percents[i];
796                if (debug)
797                        fprintf(stderr, "Finding x for chisq=%.2f:\n", chisq);
798                printf("static float chisq%d[MAXFIT+1] = { 0", (int)(100*chisq + 0.1));
799                for (nu=1; nu<=20; nu++) {
800                        ret = GCI_chisq(nu, chisq, &root);
801                        if (ret != 0)
802                                fprintf(stderr, "error %d with nu=%d, chisq=%.2f\n",
803                                                ret, nu, chisq);
804                        if (debug)
805                                fprintf(stderr, "nu=%2d  x=%6.2f  chisq(%d,%.3f)=%6.4f\n",
806                                                nu, root, nu, root, GCI_gammp(0.5*nu, 0.5*root));
807                        printf(",%s%.2f", nu % 5 == 0 ? "\n                                   " : " ", root);
808                }
809                if (debug) fprintf(stderr, "\n");
810                printf(" };\n");
811        }
812}
813#endif
814
815//***************************************** ExportParams ***********************************************/
816
817void ECF_ExportParams_start (char path[])
818{
819        ecf_exportParams = 1;
820        if (path) strcpy(ecf_exportParams_path, path);
821}
822
823void ECF_ExportParams_stop (void)
824{
825        ecf_exportParams = 0;
826}
827
828static FILE *ecf_exportFileStream;
829
830void ecf_ExportParams_OpenFile (void)
831{
832        ecf_exportFileStream = fopen(ecf_exportParams_path, "a");
833}
834
835void ecf_ExportParams_CloseFile (void)
836{
837        fprintf(ecf_exportFileStream, "\n");
838        fclose(ecf_exportFileStream);
839}
840
841void ecf_ExportParams (float param[], int nparam, float chisq)
842{
843        int i;
844
845        for (i=0; i<nparam; i++) fprintf(ecf_exportFileStream, "%g, ", param[i]);
846        fprintf(ecf_exportFileStream, "%g\n", chisq);
847}
848
849// Emacs settings:
850// Local variables:
851// mode: c
852// c-basic-offset: 4
853// tab-width: 4
854// End:
Note: See TracBrowser for help on using the repository browser.