source: branches/maven/projects/slim-curve/src/main/c/EcfUtil.c @ 7387

Revision 7387, 20.8 KB checked in by aivar, 9 years ago (diff)

A version of Paul Barber's ECF library with functions used in SLIM Plugin. These have been rewritten to remove NR dependencies. Functions not used and other code that still uses NR all removed.

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//
115// Developer: Henry Guennadi Levkin
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
208
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    if (NULL == raw)
220    {
221        return NULL;
222    }
223    float **row = (float **) raw;
224    float *data = (float *) (row + nrows);
225    int i;
226    for (i = 0; i < nrows; ++i)
227    {
228        row[i] = data;
229        data += ncols;
230    }
231    return row;
232}
233
234/* Frees allocated float matrix.
235 */
236void GCI_ecf_free_matrix(float **m)
237{
238    if (NULL != m)
239    {
240        free(m);
241    }
242}
243
244
245/********************************************************************
246
247                                FITTING FUNCTION CALCULATING FUNCTIONS
248
249*********************************************************************/
250
251
252/* Some match functions for use in the above routines
253
254   Prototypes: func(float x, float param[], float *y,
255                    float dy_dparam[], int nparam)
256   Inputs:  x            data point
257            param[]      function parameters
258            nparam       number of parameters (may or may not be
259                         used by the function)
260   Outputs: y            function value at this point
261            dy_dparam[]  vector of dy/dparam_k values at this point
262
263   There are also variants which calculate the values for a whole
264   array of x values in one go, which is significantly faster, as
265   there are many fewer calls to exp().  Note that this is only
266   reliable if tau (or lambda) is positive.  These functions are
267   called from above in certain situations to speed things up.  Their
268   prototype is:
269
270     void fitfunc_array(float xincr, float param[],
271                        float *y, float **dy_dparam, int nx, int nparam)
272
273   where the fitfunc will be evaluated for x=i*xincr, with i=0, 1,
274   ..., nx-1.  The results will be placed in the y[] array and the
275   dy_dparam[nx][nparam] array (which is assumed to have been allocated by
276   GCI_matrix).
277*/
278
279/* This one produces multiexponentials using lambdas:
280
281      y(x) = param[0] + param[1]*exp(-param[2]*x) +
282               param[3]*exp(-param[4]*x) + ...
283
284   This gives:
285
286      dy/dparam_0 = 1
287      dy/dparam_1 = exp(-param[2]*x)
288      dy/dparam_2 = -param[1]*x*exp(-param[2]*x)
289
290   and similarly for dy/dparam_3 and dy/dparam_4, etc.
291
292   As discussed above, though, we ignore the param[0] term.
293*/
294
295void GCI_multiexp_lambda(float x, float param[],
296                                                 float *y, float dy_dparam[], int nparam)
297{
298        int i;
299        float ex;
300
301        *y = 0;
302
303        for (i=1; i<nparam-1; i+=2) {
304                dy_dparam[i] = ex = exp(-param[i+1] * x);
305                ex *= param[i];
306                *y += ex;
307                dy_dparam[i+1] = -ex * x;
308        }
309}
310
311
312int multiexp_lambda_array(float xincr, float param[],
313                                                  float *y, float **dy_dparam, int nx, int nparam)
314{
315        int i, j;
316        float ex;
317        double exincr[MAXFIT];  /* exp(-lambda*xincr) */
318        double excur[MAXFIT];   /* exp(-lambda*x)     */
319
320        if (xincr <= 0) return -1;
321
322        for (j=1; j<nparam-1; j+=2) {
323                if (param[j+1] < 0) return -1;
324                excur[j] = 1.0;
325                exincr[j] = exp(-(double) param[j+1] * xincr);
326        }
327
328        for (i=0; i<nx; i++) {
329                y[i] = 0;
330                for (j=1; j<nparam-1; j+=2) {
331                        dy_dparam[i][j] = ex = excur[j];
332                        ex *= param[j];
333                        y[i] += ex;
334                        dy_dparam[i][j+1] = -ex * xincr * i;
335                        /* And ready for next loop... */
336                        excur[j] *= exincr[j];
337                }
338        }
339
340        return 0;
341}
342
343
344/* This one produces multiexponentials using taus:
345
346      y(x) = param[0] + param[1]*exp(-x/param[2]) +
347               param[3]*exp(-x/param[4]) + ...
348
349   This gives:
350
351      dy/dparam_0 = 1
352      dy/dparam_1 = exp(-x/param[2])
353      dy/dparam_2 = param[1]*x*exp(-x/param[2]) / param[2]^2
354
355   and similarly for dy/dparam_3 and dy/dparam_4, etc.
356
357   Again, we ignore the param[0] term.
358*/
359
360void GCI_multiexp_tau(float x, float param[],
361                                          float *y, float dy_dparam[], int nparam)
362{
363        int i;
364        float ex, xa;
365
366        *y = 0;
367
368        for (i=1; i<nparam-1; i+=2) {
369                xa = x / param[i+1];
370                dy_dparam[i] = ex = exp(-xa);
371                ex *= param[i];
372                *y += ex;
373                dy_dparam[i+1] = ex * xa / param[i+1];
374        }
375}
376
377
378int multiexp_tau_array(float xincr, float param[],
379                                           float *y, float **dy_dparam, int nx, int nparam)
380{
381        int i, j;
382        float ex;
383        float a2[MAXFIT];       /* 1/(param[j]*param[j]) for taus */
384        double exincr[MAXFIT];  /* exp(-xincr/tau) */
385        double excur[MAXFIT];   /* exp(-x/tau)     */
386
387        if (xincr <= 0) return -1;
388
389        for (j=1; j<nparam-1; j+=2) {
390                if (param[j+1] < 0) return -1;
391                excur[j] = 1.0;
392                exincr[j] = exp(-xincr / (double) param[j+1]);
393                a2[j] = 1 / (param[j+1] * param[j+1]);
394        }
395
396        for (i=0; i<nx; i++) {
397                y[i] = 0;
398                for (j=1; j<nparam-1; j+=2) {
399                        dy_dparam[i][j] = ex = excur[j];
400                        ex *= param[j];
401                        y[i] += ex;
402                        dy_dparam[i][j+1] = ex * xincr * i * a2[j];
403                        /* And ready for next loop... */
404                        excur[j] *= exincr[j];
405                }
406        }
407
408        return 0;
409}
410
411
412/* And this one produces stretched exponentials:
413
414      y(x) = Z + A exp(-(x/tau)^(1/h))
415
416   or translated into C-speak:
417
418      y(x) = param[0] + param[1]*exp(-(x/param[2])^(1/param[3]))
419
420   This gives:
421
422      dy/dparam_0 = 1
423      dy/dparam_1 = exp(-(x/param[2])^(1/param[3]))
424      dy/dparam_2 = param[1]*exp(-(x/param[2])^(1/param[3])) *
425                      (x/param[2])^(1/param[3]) * (1/param[2]*param[3])
426      dy/dparam_3 = param[1]*exp(-(x/param[2])^(1/param[3])) *
427                      (x/param[2])^(1/param[3]) *
428                      (1/param[3]^2)*log(x/param[2])
429*/
430
431void GCI_stretchedexp(float x, float param[],
432                                          float *y, float dy_dparam[], int nparam)
433{
434        float ex, xa, lxa, xah;
435
436        if (x > 0) {
437                xa = x / param[2];         /* xa = x/param[2] */
438                lxa = log(xa);             /* lxa = log(x/param[2]) */
439                xah = exp(lxa / param[3]); /* xah = exp(log(x/param[2])/param[3])
440                                                  = (x/param[2])^(1/param[3]) */
441                dy_dparam[1] = ex = exp(-xah);
442                                           /* ex = exp(-(x/param[2])^(1/param[3])) */
443                ex *= param[1];            /* ex = param[1] *
444                                                     exp(-(x/param[2])^(1/param[3])) */
445                *y = ex;                   /* y is now correct */
446                ex *= xah / param[3];      /* ex = param[1] * exp(...) *
447                                                      (x/param[2])^(1/param[3]) *
448                                                      1/param[3]   */
449                dy_dparam[2] = ex / param[2];
450                dy_dparam[3] = ex * lxa / param[3];
451        } else if (x > -1e-10) {  // Almost zero
452                *y = param[1];
453                dy_dparam[1] = 1;
454                dy_dparam[2] = dy_dparam[3] = 0;
455        } else {
456                /* Ouch: x<0 */
457                fprintf(stderr, "Can't have x < 0 in stretched exponential!!\n");
458        }
459}
460
461/* This is actually essentially the same as the other version; because
462   of the power (x/tau)^(1/h), we cannot use the same tricks as for
463   the multiexponential case, unfortunately. */
464
465int stretchedexp_array(float xincr, float param[],
466                                           float *y, float **dy_dparam, int nx, int nparam)
467{
468        int i;
469        float ex, lxa, xah, a2inv, a3inv;
470        double xa, xaincr;
471
472        if (xincr == 0) return -1;
473        xa=0;
474        xaincr = xincr / param[2];
475        a2inv = 1/param[2];
476        a3inv = 1/param[3];
477       
478        /* When x=0 */
479        y[0] = param[1];
480        dy_dparam[0][1] = 1;
481        dy_dparam[0][2] = dy_dparam[0][3] = 0;
482
483        for (i=1; i<nx; i++) {
484                xa += xaincr;       /* xa = (xincr*i)/param[2] */
485                lxa = log(xa);      /* lxa = log(x/param[2]) */
486                xah = exp(lxa * a3inv);  /* xah = exp(log(x/param[2])/param[3])
487                                                = (x/param[2])^(1/param[3]) */
488                dy_dparam[i][1] = ex = exp(-xah);
489                                    /* ex = exp(-(x/param[2])^(1/param[3])) */
490                ex *= param[1];     /* ex = param[1]*exp(-(x/param[2])^(1/param[3])) */
491                y[i] = ex;          /* y is now correct */
492                ex *= xah * a3inv;  /* ex = param[1] * exp(...) *
493                                              (x/param[2])^(1/param[3]) * 1/param[3] */
494                dy_dparam[i][2] = ex * a2inv;
495                dy_dparam[i][3] = ex * lxa * a3inv;
496        }
497
498        return 0;
499}
500
501
502/********************************************************************
503
504                           CHECKING AND RESTRICTED FITTING ROUTINES
505
506*********************************************************************/
507
508
509/* This is the default routine */
510#define MIN_Z -1e5
511#define MIN_Z_FACTOR 0.4
512#define MAX_Z 1e10  /* Silly */
513#define MIN_A 0
514#define MAX_A 1e10  /* Silly */
515#define MIN_TAU 0.001  /* Works for lambda too */
516#define MAX_TAU 1000
517#define MIN_H 1
518#define MAX_H 10
519
520int check_ecf_params (float param[], int nparam,
521                                        void (*fitfunc)(float, float [], float *, float [], int))
522{
523        if (fitfunc == GCI_multiexp_lambda || fitfunc == GCI_multiexp_tau) {
524                switch (nparam) {
525                case 3:
526                        if (param[0] < MIN_Z ||
527                                param[0] < -MIN_Z_FACTOR * fabs(param[1]) ||
528                                param[0] > MAX_Z)
529                                return -21;
530                        if (param[1] < MIN_A || param[1] > MAX_A)
531                                return -22;
532                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
533                                return -23;
534                        break;
535
536                case 5:
537                        if (param[0] < MIN_Z ||
538                                param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3])) ||
539                                param[0] > MAX_Z)
540                                return -21;
541                        if (param[1] < MIN_A || param[1] > MAX_A)
542                                return -22;
543                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
544                                return -23;
545                        if (param[3] < MIN_A || param[3] > MAX_A)
546                                return -24;
547                        if (param[4] < MIN_TAU || param[4] > MAX_TAU)
548                                return -25;
549                        break;
550
551                case 7:
552                        if (param[0] < MIN_Z ||
553                                param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3]) +
554                                                                                        fabs(param[5])) ||
555                                param[0] > MAX_Z)
556                                return -21;
557                        if (param[1] < MIN_A || param[1] > MAX_A)
558                                return -22;
559                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
560                                return -23;
561                        if (param[3] < MIN_A || param[3] > MAX_A)
562                                return -24;
563                        if (param[4] < MIN_TAU || param[4] > MAX_TAU)
564                                return -25;
565                        if (param[5] < MIN_A || param[5] > MAX_A)
566                                return -26;
567                        if (param[6] < MIN_TAU || param[6] > MAX_TAU)
568                                return -27;
569                        break;
570                }
571        } else if (fitfunc == GCI_stretchedexp) {
572                if (param[0] < MIN_Z || param[0] < -MIN_Z_FACTOR * fabs(param[1]) ||
573                        param[0] > MAX_Z)
574                        return -21;
575                if (param[1] < MIN_A || param[1] > MAX_A)
576                        return -22;
577                if (param[2] < MIN_TAU || param[2] > MAX_TAU)
578                        return -23;
579                if (param[3] < MIN_H || param[3] > MAX_H)
580                        return -24;
581        }
582        return 0;
583}
584
585
586/* For the user-specified version, we have some global variables to
587   store the settings between invocations. */
588
589static int restrain_nparam=0;  /* How many parameters have we set up? */
590static int restrain_restraining[MAXFIT];  /* Do we check parameter i? */
591static float restrain_minval[MAXFIT]; /* Minimum acceptable parameter value */
592static float restrain_maxval[MAXFIT]; /* Maximum acceptable parameter value */
593
594int GCI_set_restrain_limits(int nparam, int restrain[],
595                                                        float minval[], float maxval[])
596{
597        int i;
598
599        if (nparam < 0 || nparam > MAXFIT)
600                return -1;
601
602        /* We're going to be doing something, so clear the memory */
603        restrain_nparam = 0;
604
605        for (i=0; i<nparam; i++) {
606                if (restrain[i]) {
607                        restrain_restraining[i] = 1;
608
609                        if (minval[i] > maxval[i])
610                                return -2;
611                        restrain_minval[i] = minval[i];
612                        restrain_maxval[i] = maxval[i];
613                } else
614                        restrain_restraining[i] = 0;
615        }
616
617        restrain_nparam = nparam;
618        return 0;
619}
620
621/* original version
622int check_ecf_user_params (float param[], int nparam,
623                                        void (*fitfunc)(float, float [], float *, float [], int))
624{
625        int i;
626
627        if (restrain_nparam != nparam) {
628                dbgprintf(0, "Using user-defined parameter restraining with "
629                                  "wrong number of parameters:\n"
630                                  "actual nparam = %d, user restraining nparam = %d\n"
631                                  "Defaulting to standard tests\n", nparam, restrain_nparam);
632                return check_ecf_params(param, nparam, fitfunc);
633        }
634
635        for (i=0; i<nparam; i++) {
636                if (restrain_restraining[i]) {
637                        if ((param[i] < restrain_minval[i]) ||
638                                (param[i] > restrain_maxval[i]))
639                                return -21;
640                }
641        }
642
643        return 0;
644}
645*/
646// new version from J Gilbey 31.03.03
647int check_ecf_user_params (float param[], int nparam,
648                                        void (*fitfunc)(float, float [], float *, float [], int))
649{
650        int i;
651
652
653        if (restrain_nparam != nparam) {
654                dbgprintf(0, "Using user-defined parameter restraining with "
655                                  "wrong number of parameters:\n"
656                                  "actual nparam = %d, user restraining nparam = %d\n"
657                                  "Defaulting to standard tests\n", nparam, restrain_nparam);
658                return check_ecf_params(param, nparam, fitfunc);
659        }
660
661
662        for (i=0; i<nparam; i++) {
663                if (restrain_restraining[i]) {
664                        if (param[i] < restrain_minval[i])
665                                param[i] = restrain_minval[i];
666                        else if (param[i] > restrain_maxval[i])
667                                param[i] = restrain_maxval[i];
668                }
669        }
670
671
672        return 0;
673}
674
675
676/********************************************************************
677
678                                                          MISCELLANY
679
680*********************************************************************/
681
682
683/* From filters */
684int ECF_Find_Float_Max (float data[], int np, float *max_val)
685{
686        int i, maxi;
687        float *ptr=data;
688
689        if (np < 1)        /* daft input */
690                return -1;
691
692        maxi = 0;
693        *max_val = *ptr;
694
695        for (i=1; i<np; i++)
696                if (*++ptr>*max_val) {
697                        *max_val = *ptr;
698                        maxi = i;
699                }
700
701        return maxi;
702}
703
704
705/* Debugging version of printf */
706#ifdef _CVI_
707int dbgprintf(int dbg_level, const char *format, ...)
708{
709        int ret = 0;
710        char msg[1000];
711        va_list va;
712
713        if (ECF_debug >= dbg_level) {
714                va_start(va, format);
715                ret = vsprintf(msg, format, va);  // Here's hoping, no vsnprintf ...
716                MessagePopup("Debug info", msg);
717                va_end(va);
718        }
719
720        return ret;
721}
722#else
723int dbgprintf(int dbg_level, const char *format, ...)
724{
725        int ret = 0;
726        va_list va;
727
728        if (ECF_debug >= dbg_level) {
729                va_start(va, format);
730                ret = vfprintf(stderr, format, va);
731                va_end(va);
732        }
733
734        return ret;
735}
736#endif
737
738#ifdef TESTCHISQ
739int main(int ac, char **av)
740{
741        int i, nu, ret;
742        int debug = 0;
743        float root, chisq;
744        float percents[4] = { 0.50, 0.68, 0.90, 0.95 };
745
746        for (i=0; i<4; i++) {
747                chisq = percents[i];
748                if (debug)
749                        fprintf(stderr, "Finding x for chisq=%.2f:\n", chisq);
750                printf("static float chisq%d[MAXFIT+1] = { 0", (int)(100*chisq + 0.1));
751                for (nu=1; nu<=20; nu++) {
752                        ret = GCI_chisq(nu, chisq, &root);
753                        if (ret != 0)
754                                fprintf(stderr, "error %d with nu=%d, chisq=%.2f\n",
755                                                ret, nu, chisq);
756                        if (debug)
757                                fprintf(stderr, "nu=%2d  x=%6.2f  chisq(%d,%.3f)=%6.4f\n",
758                                                nu, root, nu, root, GCI_gammp(0.5*nu, 0.5*root));
759                        printf(",%s%.2f", nu % 5 == 0 ? "\n                                   " : " ", root);
760                }
761                if (debug) fprintf(stderr, "\n");
762                printf(" };\n");
763        }
764}
765#endif
766
767//***************************************** ExportParams ***********************************************/
768
769void ECF_ExportParams_start (char path[])
770{
771        ecf_exportParams = 1;
772        if (path) strcpy(ecf_exportParams_path, path);
773}
774
775void ECF_ExportParams_stop (void)
776{
777        ecf_exportParams = 0;
778}
779
780static FILE *ecf_exportFileStream;
781
782void ecf_ExportParams_OpenFile (void)
783{
784        ecf_exportFileStream = fopen(ecf_exportParams_path, "a"); 
785}
786
787void ecf_ExportParams_CloseFile (void)
788{
789        fprintf(ecf_exportFileStream, "\n");
790        fclose(ecf_exportFileStream);
791}
792
793void ecf_ExportParams (float param[], int nparam, float chisq)
794{
795        int i;
796       
797        for (i=0; i<nparam; i++) fprintf(ecf_exportFileStream, "%g, ", param[i]);
798        fprintf(ecf_exportFileStream, "%g\n", chisq);
799}
800
801// Emacs settings:
802// Local variables:
803// mode: c
804// c-basic-offset: 4
805// tab-width: 4
806// End:
Note: See TracBrowser for help on using the repository browser.