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

Revision 7780, 27.2 KB checked in by aivar, 8 years ago (diff)

Tidied up.

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 <string.h>
16#include <math.h>
17#ifdef _CVI_
18#include <userint.h>
19#endif
20#include "EcfInternal.h"  /* For #defines */
21
22int ECF_debug = 0;
23
24/********************************************************************
25
26                           EQUATION-SOLVING AND COVAR-SORTING ROUTINES
27
28*********************************************************************/
29
30#define SWAP(a,b) { temp=(a); (a)=(b); (b)=temp; }
31
32/* Linear equation solution of Ax = b by Gaussian elimination.
33   A is the n x n input matrix, b is the right-hand side vector, length n.
34   On output, b is replaced by the corresponding set of solution vectors
35   and A is trashed.
36 */
37int GCI_solve_Gaussian(float **a, int n, float *b)
38{
39    float max;
40    float temp;
41    float *pivotInverse = (float *)malloc(n * sizeof(float));
42    int i, j, k, m;
43
44    // base row of matrix
45    for (k = 0; k < n - 1; ++k)
46    {
47        // search for line with max element
48        max = fabs(a[k][k]);
49        m = k;
50        for (i = k + 1; i < n; ++i)
51        {
52            if (max < fabs(a[i][k])) // row i col k
53            {
54                max = a[i][k];
55                m = i;
56            }
57        }
58
59        // permutation of base line (index k) and max element line (index m)
60        if (m != k)
61        {
62            for (i = k; i < n; ++i)
63            {
64                SWAP(a[k][i], a[m][i]);
65            }
66            SWAP(b[k], b[m]);
67        }
68
69        if (0.0 == a[k][k])
70        {
71            free(pivotInverse);
72            return -2; // singular matrix
73        }
74
75        // triangulation of matrix with coefficients
76        pivotInverse[k] = 1.0 / a[k][k];
77        for (j = k + 1; j < n; ++j) // current row of matrix
78        {
79            // want "temp = -a[j][k] / a[k][k]"
80            temp = -a[j][k] * pivotInverse[k];
81            for (i = k; i < n; ++i)
82            {
83                a[j][i] += temp * a[k][i];
84            }
85            b[j] += temp * b[k]; // free member recalculation
86        }
87    }
88    // precalculate last pivot inverse
89    pivotInverse[n - 1] = 1.0 / a[n - 1][n - 1];
90
91    for (k = n - 1; k >= 0; --k)
92    {
93        for (i = k + 1; i < n; ++i)
94        {
95            b[k] -= a[k][i] * b[i];
96        }
97        b[k] *= pivotInverse[k];
98    }
99
100    free(pivotInverse);
101    return 0;
102}
103
104/* Matrix inversion by Gaussian elimination.
105   A is the n x n input matrix.
106   On output, A is replaced by its matrix inverse.
107   Returns 0 upon success, -2 if matrix is singular.
108 */
109int GCI_invert_Gaussian(float **a, int n)
110{
111    int returnValue = 0;
112    float **identity = GCI_ecf_matrix(n, n);
113    float **work = GCI_ecf_matrix(n, n);
114    int i, j, k;
115
116    for (j = 0; j < n; ++j) {
117        // find inverse by columns
118        for (i = 0; i < n; ++i) {
119            identity[j][i] = 0.0;
120            // need a fresh copy of matrix a
121            for (k = 0; k < n; ++k) {
122                work[k][i] = a[k][i];
123            }
124        }
125        identity[j][j] = 1.0;
126        returnValue = GCI_solve_Gaussian(work, n, identity[j]);
127        if (returnValue < 0) {
128            return returnValue;
129        }
130    }
131
132    // copy over results
133    for (j = 0; j < n; ++j) {
134        for (i = 0; i < n; ++i) {
135            a[j][i] = identity[j][i];
136        }
137    }
138
139    GCI_ecf_free_matrix(identity);
140    GCI_ecf_free_matrix(work);
141    return returnValue;
142}
143
144/* Pivots matrix on given column.  Also keeps track of order.
145 */
146void pivot(float **a, int n, int *order, int col)
147{
148    int pivotRow;
149    float maxValue;
150    float rowValue;
151    int i;
152
153    // find row with maximum element value in col, below diagonal
154    pivotRow = col;
155    maxValue = fabs(a[col][col]);
156    for (i = col + 1; i < n; ++i) {
157        rowValue = fabs(a[i][col]);
158        if (rowValue > maxValue) {
159            pivotRow = i;
160            maxValue = rowValue;
161        }
162    }
163
164    // swap rows
165    if (pivotRow != col) {
166        // swap elements in a matrix
167        for (i = 0; i < n; ++i) {
168            float temp;
169            SWAP(a[col][i], a[pivotRow][i]);
170        }
171
172        // swap elements in order vector
173        {
174            int temp;
175            SWAP(order[col], order[pivotRow]);
176        }
177    }
178}
179
180/*
181  Performs an in-place Crout lower/upper decomposition of n x n matrix A.
182  Values on or below diagonals are lowers, values about the
183  diagonal are uppers, with an implicit 1.0 value for the
184  diagonals.
185  Returns 0 upon success, -2 if matrix is singular.
186
187  Based on _Applied Numerical Analysis_, Fourth Edition, Gerald & Wheatley
188  Sections 2.5 & 2.14.
189 */
190int lu_decomp(float **a, int n, int *order)
191{
192    int i;
193    float inverse;
194    int jCol;
195    int iRow;
196    int kCol;
197    float sum;
198
199    // initialize ordering vector
200    for (i = 0; i < n; ++i)
201    {
202        order[i] = i;
203    }
204
205    // pivot first column
206    pivot(a, n, order, 0);
207
208    // check for singularity
209    if (0.0 == a[0][0])
210    {
211        return -2;
212    }
213
214    // compute first row of upper
215    inverse = 1.0 / a[0][0];
216    for (i = 1; i < n; ++i) {
217        a[0][i] *= inverse;
218    }
219
220    // continue computing columns of lowers then rows of uppers
221    for (jCol = 1; jCol < n - 1; ++jCol) {
222        // compute column of lowers
223        for (iRow = jCol; iRow < n; ++iRow) {
224            sum = 0.0;
225            for (kCol = 0; kCol < jCol; ++kCol) {
226                sum += a[iRow][kCol] * a[kCol][jCol];
227            }
228            a[iRow][jCol] -= sum;
229        }
230
231        // find pivot for row
232        pivot(a, n, order, jCol);
233        if (0.0 == a[jCol][jCol])
234        {
235            return -2;
236        }
237
238        // build row of uppers
239        inverse = 1.0 / a[jCol][jCol];
240        for (kCol = jCol + 1; kCol < n; ++kCol) {
241            sum = 0.0;
242            for (iRow = 0; iRow < jCol; ++iRow) {
243                sum += a[jCol][iRow] * a[iRow][kCol];
244            }
245            a[jCol][kCol] -= sum;
246            a[jCol][kCol] *= inverse;
247        }
248    }
249
250    // get remaining lower
251    sum = 0.0;
252    for (kCol = 0; kCol < n - 1; ++kCol) {
253        sum += a[n - 1][kCol] * a[kCol][n - 1];
254    }
255    a[n - 1][n - 1] -= sum;
256    return 0;
257}
258
259/*
260 Given a LU decomposition of an n x n matrix A, and an order vector
261 specifying any reordering done during pivoting, solves the equation
262 Ax = b.  Returns result in b.
263 */
264//TODO check for lu[i][i] != 0?
265int solve_lu(float **lu, int n, float *b, int *order)
266{
267    int startIndex;
268    int index;
269    float temp;
270    int iRow;
271    int jCol;
272    int nvbl;
273    float sum;
274
275    // rearrange the elements of the b vector in place.
276    startIndex = order[0];
277    index = startIndex;
278    temp = b[index];
279    while (1) {
280        int nextIndex = order[index];
281        if (nextIndex == startIndex) {
282            b[index] = temp;
283            break;
284        }
285        b[index] = b[nextIndex];
286        index = nextIndex;
287    }
288
289    // compute the b' vector
290    b[0] /= lu[0][0];
291    for (iRow = 1; iRow < n; ++iRow) {
292        sum = 0.0;
293        for (jCol = 0; jCol < iRow; ++jCol) {
294            sum += lu[iRow][jCol] * b[jCol];
295        }
296        b[iRow] -= sum;
297        b[iRow] /= lu[iRow][iRow];
298    }
299
300    // get the solution, we have b[n-1] already
301    for (iRow = 1; iRow < n; ++iRow) { // iRow goes from 1 to n-1
302        nvbl = n - iRow - 1;           // nvbl goes from n-2 to 0
303        sum = 0.0f;
304        for (jCol = nvbl + 1; jCol < n; ++jCol) {
305            sum += lu[nvbl][jCol] * b[jCol];
306        }
307        b[nvbl] -= sum;
308    }
309    return 0;
310}
311
312/* Linear equation solution of Ax = b by lower/upper decomposition.
313   A is the n x n input max, b is the right-hand side vector, length n.
314   On output, b is replaced by the corresponding set of solution vectors
315   and A is trashed.
316 */
317int GCI_solve_lu_decomp(float **a, int n, float *b)
318{
319    int *order = (int *) malloc(n * sizeof(int));
320    int return_value = lu_decomp(a, n, order);
321    if (return_value >= 0) {
322        return_value = solve_lu(a, n, b, order);
323    }
324    free(order);
325    return return_value;
326}
327
328/* Matrix inversion by lower/upper decomposition.
329   A is the n x n input matrix.
330   On output, a is replaced by its matrix inverse..
331 */
332int GCI_invert_lu_decomp(float **a, int n)
333{
334    int returnValue;
335    int *order = (int *) malloc(n * sizeof(int));
336    float **identity = GCI_ecf_matrix(n, n);
337    int i, j;
338
339    returnValue = lu_decomp(a, n, order);
340    if (returnValue >= 0) {
341        for (j = 0; j < n; ++j) {
342            // find inverse by columns
343            for (i = 0; i < n; ++i) {
344                identity[j][i] = 0.0;
345            }
346            identity[j][j] = 1.0;
347            solve_lu(a, n, identity[j], order);
348        }
349        for (j = 0; j < n; ++j) {
350            for (i = 0; i < n; ++i) {
351                a[j][i] = identity[j][i];
352            }
353        }
354    }
355
356    free(order);
357    GCI_ecf_free_matrix(identity);
358    return returnValue;
359}
360
361/* Linear equation solution of Ax = b..
362   A is the n x n input max, b is the right-hand side vector, length n.
363   On output, b is replaced by the corresponding set of solution vectors
364   and A is trashed.
365 */
366int GCI_solve(float **a, int n, float *b)
367{
368    return GCI_solve_Gaussian(a, n, b);
369    //return GCI_solve_lu_decomp(a, n, b);
370}
371
372/* Matrix inversion.
373   A is the n x n input matrix.
374   On output, a is replaced by its matrix inverse..
375 */
376int GCI_invert(float **a, int n)
377{
378    return GCI_invert_Gaussian(a, n);
379    //return GCI_invert_lu_decomp(a, n);
380}
381
382
383void GCI_covar_sort(float **covar, int nparam, int paramfree[], int mfit)
384{
385        int i, j, k;
386        float temp;
387
388        for (i=mfit; i<nparam; i++)
389                for (j=0; j<=i; j++)
390                        covar[i][j] = covar[j][i] = 0.0;
391
392        k = mfit-1;
393        for (j=nparam-1; j>=0; j--)
394        {
395                if (paramfree[j])
396                {
397                        for (i=0; i<nparam; i++) SWAP(covar[i][k], covar[i][j]);
398                        for (i=0; i<nparam; i++) SWAP(covar[k][i], covar[j][i]);
399                        k--;
400                }
401        }
402}
403#undef SWAP
404
405
406/********************************************************************
407
408                                          ARRAY ALLOCATION ROUTINES
409http://www.nr.com/public-domain.html
410*********************************************************************/
411
412/* This function allocates a float matrix with subscript range
413 * m[0..nrows-1][0..ncols-1]
414 */
415float **GCI_ecf_matrix(long nrows, long ncols)
416{
417    int row_size = nrows * sizeof(float *);
418    int data_size = nrows * ncols * sizeof(float);
419    unsigned char *raw = malloc(row_size + data_size);
420    float **row = (float **) raw;
421    float *data = (float *) (row + nrows);
422    int i;
423
424        if (NULL == raw)
425    {
426        return NULL;
427    }
428
429        for (i = 0; i < nrows; ++i)
430    {
431        row[i] = data;
432        data += ncols;
433    }
434    return row;
435}
436
437/* Frees allocated float matrix.
438 */
439void GCI_ecf_free_matrix(float **m)
440{
441    if (NULL != m)
442    {
443        free(m);
444    }
445}
446
447float ***GCI_ecf_matrix_array(long nblocks, long nrows, long ncols)
448/* allocate a float matrix array with range
449   marr[0..nblocks][0..nrows][0..ncols] */
450{
451        long i;
452        float ***marr;
453
454        /* allocate pointers to blocks */
455        if ((marr = (float ***) malloc(nblocks * sizeof(float **))) == NULL)
456                return NULL;
457
458        /* allocate blocks (= pointers to rows) and set pointers to them */
459        if ((marr[0] = (float **) malloc(nblocks * nrows * sizeof(float *)))
460                == NULL) {
461                free(marr);
462                return NULL;
463        }
464
465        for (i=1; i<nblocks; i++)
466                marr[i] = marr[i-1] + nrows;
467
468        /* allocate rows (= pointers to column entries) and set pointers to them */
469        if ((marr[0][0] = (float *)malloc(nblocks * nrows * ncols * sizeof(float)))
470                == NULL) {
471                free(marr[0]);
472                free(marr);
473                return NULL;
474        }
475
476        /* This sneaky loop does the whole lot!  For note that
477           marr[block][row] = marr[0][block*nrows + row] */
478        for (i=1; i < nblocks * nrows; i++)
479                marr[0][i] = marr[0][i-1] + ncols;
480
481        return marr;
482}
483
484
485void GCI_ecf_free_matrix_array(float ***marr)
486{
487        if (marr != NULL) {
488                free(marr[0][0]);
489                free(marr[0]);
490                free(marr);
491        }
492}
493
494/********************************************************************
495
496                                FITTING FUNCTION CALCULATING FUNCTIONS
497
498*********************************************************************/
499
500
501/* Some match functions for use in the above routines
502
503   Prototypes: func(float x, float param[], float *y,
504                    float dy_dparam[], int nparam)
505   Inputs:  x            data point
506            param[]      function parameters
507            nparam       number of parameters (may or may not be
508                         used by the function)
509   Outputs: y            function value at this point
510            dy_dparam[]  vector of dy/dparam_k values at this point
511
512   There are also variants which calculate the values for a whole
513   array of x values in one go, which is significantly faster, as
514   there are many fewer calls to exp().  Note that this is only
515   reliable if tau (or lambda) is positive.  These functions are
516   called from above in certain situations to speed things up.  Their
517   prototype is:
518
519     void fitfunc_array(float xincr, float param[],
520                        float *y, float **dy_dparam, int nx, int nparam)
521
522   where the fitfunc will be evaluated for x=i*xincr, with i=0, 1,
523   ..., nx-1.  The results will be placed in the y[] array and the
524   dy_dparam[nx][nparam] array (which is assumed to have been allocated by
525   GCI_matrix).
526*/
527
528/* This one produces multiexponentials using lambdas:
529
530      y(x) = param[0] + param[1]*exp(-param[2]*x) +
531               param[3]*exp(-param[4]*x) + ...
532
533   This gives:
534
535      dy/dparam_0 = 1
536      dy/dparam_1 = exp(-param[2]*x)
537      dy/dparam_2 = -param[1]*x*exp(-param[2]*x)
538
539   and similarly for dy/dparam_3 and dy/dparam_4, etc.
540
541   As discussed above, though, we ignore the param[0] term.
542*/
543
544void GCI_multiexp_lambda(float x, float param[],
545                                                 float *y, float dy_dparam[], int nparam)
546{
547        int i;
548        float ex;
549
550        *y = 0;
551
552        for (i=1; i<nparam-1; i+=2) {
553                dy_dparam[i] = ex = exp(-param[i+1] * x);
554                ex *= param[i];
555                *y += ex;
556                dy_dparam[i+1] = -ex * x;
557        }
558}
559
560
561int multiexp_lambda_array(float xincr, float param[],
562                                                  float *y, float **dy_dparam, int nx, int nparam)
563{
564        int i, j;
565        float ex;
566        double exincr[MAXFIT];  /* exp(-lambda*xincr) */
567        double excur[MAXFIT];   /* exp(-lambda*x)     */
568
569        if (xincr <= 0) return -1;
570
571        for (j=1; j<nparam-1; j+=2) {
572                if (param[j+1] < 0) return -1;
573                excur[j] = 1.0;
574                exincr[j] = exp(-(double) param[j+1] * xincr);
575        }
576
577        for (i=0; i<nx; i++) {
578                y[i] = 0;
579                for (j=1; j<nparam-1; j+=2) {
580                        dy_dparam[i][j] = ex = excur[j];
581                        ex *= param[j];
582                        y[i] += ex;
583                        dy_dparam[i][j+1] = -ex * xincr * i;
584                        /* And ready for next loop... */
585                        excur[j] *= exincr[j];
586                }
587        }
588
589        return 0;
590}
591
592
593/* This one produces multiexponentials using taus:
594
595      y(x) = param[0] + param[1]*exp(-x/param[2]) +
596               param[3]*exp(-x/param[4]) + ...
597
598   This gives:
599
600      dy/dparam_0 = 1
601      dy/dparam_1 = exp(-x/param[2])
602      dy/dparam_2 = param[1]*x*exp(-x/param[2]) / param[2]^2
603
604   and similarly for dy/dparam_3 and dy/dparam_4, etc.
605
606   Again, we ignore the param[0] term.
607*/
608
609void GCI_multiexp_tau(float x, float param[],
610                                          float *y, float dy_dparam[], int nparam)
611{
612        int i;
613        float ex, xa;
614
615        *y = 0;
616
617        for (i=1; i<nparam-1; i+=2) {
618                xa = x / param[i+1];
619                dy_dparam[i] = ex = exp(-xa);
620                ex *= param[i];
621                *y += ex;
622                dy_dparam[i+1] = ex * xa / param[i+1];
623        }
624}
625
626
627int multiexp_tau_array(float xincr, float param[],
628                                           float *y, float **dy_dparam, int nx, int nparam)
629{
630        int i, j;
631        float ex;
632        float a2[MAXFIT];       /* 1/(param[j]*param[j]) for taus */
633        double exincr[MAXFIT];  /* exp(-xincr/tau) */
634        double excur[MAXFIT];   /* exp(-x/tau)     */
635
636        if (xincr <= 0) return -1;
637
638        for (j=1; j<nparam-1; j+=2) {
639                if (param[j+1] < 0) return -1;
640                excur[j] = 1.0;
641                exincr[j] = exp(-xincr / (double) param[j+1]);
642                a2[j] = 1 / (param[j+1] * param[j+1]);
643        }
644
645        for (i=0; i<nx; i++) {
646                y[i] = 0;
647                for (j=1; j<nparam-1; j+=2) {
648                        dy_dparam[i][j] = ex = excur[j];
649                        ex *= param[j];
650                        y[i] += ex;
651                        dy_dparam[i][j+1] = ex * xincr * i * a2[j];
652                        /* And ready for next loop... */
653                        excur[j] *= exincr[j];
654                }
655        }
656
657        return 0;
658}
659
660
661/* And this one produces stretched exponentials:
662
663      y(x) = Z + A exp(-(x/tau)^(1/h))
664
665   or translated into C-speak:
666
667      y(x) = param[0] + param[1]*exp(-(x/param[2])^(1/param[3]))
668
669   This gives:
670
671      dy/dparam_0 = 1
672      dy/dparam_1 = exp(-(x/param[2])^(1/param[3]))
673      dy/dparam_2 = param[1]*exp(-(x/param[2])^(1/param[3])) *
674                      (x/param[2])^(1/param[3]) * (1/param[2]*param[3])
675      dy/dparam_3 = param[1]*exp(-(x/param[2])^(1/param[3])) *
676                      (x/param[2])^(1/param[3]) *
677                      (1/param[3]^2)*log(x/param[2])
678*/
679
680void GCI_stretchedexp(float x, float param[],
681                                          float *y, float dy_dparam[], int nparam)
682{
683        float ex, xa, lxa, xah;
684
685        if (x > 0) {
686                xa = x / param[2];         /* xa = x/param[2] */
687                lxa = log(xa);             /* lxa = log(x/param[2]) */
688                xah = exp(lxa / param[3]); /* xah = exp(log(x/param[2])/param[3])
689                                                  = (x/param[2])^(1/param[3]) */
690                dy_dparam[1] = ex = exp(-xah);
691                                           /* ex = exp(-(x/param[2])^(1/param[3])) */
692                ex *= param[1];            /* ex = param[1] *
693                                                     exp(-(x/param[2])^(1/param[3])) */
694                *y = ex;                   /* y is now correct */
695                ex *= xah / param[3];      /* ex = param[1] * exp(...) *
696                                                      (x/param[2])^(1/param[3]) *
697                                                      1/param[3]   */
698                dy_dparam[2] = ex / param[2];
699                dy_dparam[3] = ex * lxa / param[3];
700        } else if (x > -1e-10) {  // Almost zero
701                *y = param[1];
702                dy_dparam[1] = 1;
703                dy_dparam[2] = dy_dparam[3] = 0;
704        } else {
705                /* Ouch: x<0 */
706                fprintf(stderr, "Can't have x < 0 in stretched exponential!!\n");
707        }
708}
709
710/* This is actually essentially the same as the other version; because
711   of the power (x/tau)^(1/h), we cannot use the same tricks as for
712   the multiexponential case, unfortunately. */
713
714int stretchedexp_array(float xincr, float param[],
715                                           float *y, float **dy_dparam, int nx, int nparam)
716{
717        int i;
718        float ex, lxa, xah, a2inv, a3inv;
719        double xa, xaincr;
720
721        if (xincr == 0) return -1;
722        xa=0;
723        xaincr = xincr / param[2];
724        a2inv = 1/param[2];
725        a3inv = 1/param[3];
726
727        /* When x=0 */
728        y[0] = param[1];
729        dy_dparam[0][1] = 1;
730        dy_dparam[0][2] = dy_dparam[0][3] = 0;
731
732        for (i=1; i<nx; i++) {
733                xa += xaincr;       /* xa = (xincr*i)/param[2] */
734                lxa = log(xa);      /* lxa = log(x/param[2]) */
735                xah = exp(lxa * a3inv);  /* xah = exp(log(x/param[2])/param[3])
736                                                = (x/param[2])^(1/param[3]) */
737                dy_dparam[i][1] = ex = exp(-xah);
738                                    /* ex = exp(-(x/param[2])^(1/param[3])) */
739                ex *= param[1];     /* ex = param[1]*exp(-(x/param[2])^(1/param[3])) */
740                y[i] = ex;          /* y is now correct */
741                ex *= xah * a3inv;  /* ex = param[1] * exp(...) *
742                                              (x/param[2])^(1/param[3]) * 1/param[3] */
743                dy_dparam[i][2] = ex * a2inv;
744                dy_dparam[i][3] = ex * lxa * a3inv;
745        }
746
747        return 0;
748}
749
750
751/********************************************************************
752
753                           CHECKING AND RESTRICTED FITTING ROUTINES
754
755*********************************************************************/
756
757
758/* This is the default routine */
759#define MIN_Z -1e5
760#define MIN_Z_FACTOR 0.4
761#define MAX_Z 1e10  /* Silly */
762#define MIN_A 0
763#define MAX_A 1e10  /* Silly */
764#define MIN_TAU 0.001  /* Works for lambda too */
765#define MAX_TAU 1000
766#define MIN_H 1
767#define MAX_H 10
768
769int check_ecf_params (float param[], int nparam,
770                                        void (*fitfunc)(float, float [], float *, float [], int))
771{
772        if (fitfunc == GCI_multiexp_lambda || fitfunc == GCI_multiexp_tau) {
773                switch (nparam) {
774                case 3:
775                        if (param[0] < MIN_Z ||
776                                param[0] < -MIN_Z_FACTOR * fabs(param[1]) ||
777                                param[0] > MAX_Z)
778                                return -21;
779                        if (param[1] < MIN_A || param[1] > MAX_A)
780                                return -22;
781                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
782                                return -23;
783                        break;
784
785                case 5:
786                        if (param[0] < MIN_Z ||
787                                param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3])) ||
788                                param[0] > MAX_Z)
789                                return -21;
790                        if (param[1] < MIN_A || param[1] > MAX_A)
791                                return -22;
792                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
793                                return -23;
794                        if (param[3] < MIN_A || param[3] > MAX_A)
795                                return -24;
796                        if (param[4] < MIN_TAU || param[4] > MAX_TAU)
797                                return -25;
798                        break;
799
800                case 7:
801                        if (param[0] < MIN_Z ||
802                                param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3]) +
803                                                                                        fabs(param[5])) ||
804                                param[0] > MAX_Z)
805                                return -21;
806                        if (param[1] < MIN_A || param[1] > MAX_A)
807                                return -22;
808                        if (param[2] < MIN_TAU || param[2] > MAX_TAU)
809                                return -23;
810                        if (param[3] < MIN_A || param[3] > MAX_A)
811                                return -24;
812                        if (param[4] < MIN_TAU || param[4] > MAX_TAU)
813                                return -25;
814                        if (param[5] < MIN_A || param[5] > MAX_A)
815                                return -26;
816                        if (param[6] < MIN_TAU || param[6] > MAX_TAU)
817                                return -27;
818                        break;
819                }
820        } else if (fitfunc == GCI_stretchedexp) {
821                if (param[0] < MIN_Z || param[0] < -MIN_Z_FACTOR * fabs(param[1]) ||
822                        param[0] > MAX_Z)
823                        return -21;
824                if (param[1] < MIN_A || param[1] > MAX_A)
825                        return -22;
826                if (param[2] < MIN_TAU || param[2] > MAX_TAU)
827                        return -23;
828                if (param[3] < MIN_H || param[3] > MAX_H)
829                        return -24;
830        }
831        return 0;
832}
833
834
835/* For the user-specified version, we have some global variables to
836   store the settings between invocations. */
837
838static int restrain_nparam=0;  /* How many parameters have we set up? */
839static int restrain_restraining[MAXFIT];  /* Do we check parameter i? */
840static float restrain_minval[MAXFIT]; /* Minimum acceptable parameter value */
841static float restrain_maxval[MAXFIT]; /* Maximum acceptable parameter value */
842
843int GCI_set_restrain_limits(int nparam, int restrain[],
844                                                        float minval[], float maxval[])
845{
846        int i;
847
848        if (nparam < 0 || nparam > MAXFIT)
849                return -1;
850
851        /* We're going to be doing something, so clear the memory */
852        restrain_nparam = 0;
853
854        for (i=0; i<nparam; i++) {
855                if (restrain[i]) {
856                        restrain_restraining[i] = 1;
857
858                        if (minval[i] > maxval[i])
859                                return -2;
860                        restrain_minval[i] = minval[i];
861                        restrain_maxval[i] = maxval[i];
862                } else
863                        restrain_restraining[i] = 0;
864        }
865
866        restrain_nparam = nparam;
867        return 0;
868}
869
870/* original version
871int check_ecf_user_params (float param[], int nparam,
872                                        void (*fitfunc)(float, float [], float *, float [], int))
873{
874        int i;
875
876        if (restrain_nparam != nparam) {
877                dbgprintf(0, "Using user-defined parameter restraining with "
878                                  "wrong number of parameters:\n"
879                                  "actual nparam = %d, user restraining nparam = %d\n"
880                                  "Defaulting to standard tests\n", nparam, restrain_nparam);
881                return check_ecf_params(param, nparam, fitfunc);
882        }
883
884        for (i=0; i<nparam; i++) {
885                if (restrain_restraining[i]) {
886                        if ((param[i] < restrain_minval[i]) ||
887                                (param[i] > restrain_maxval[i]))
888                                return -21;
889                }
890        }
891
892        return 0;
893}
894*/
895// new version from J Gilbey 31.03.03
896int check_ecf_user_params (float param[], int nparam,
897                                        void (*fitfunc)(float, float [], float *, float [], int))
898{
899        int i;
900
901
902        if (restrain_nparam != nparam) {
903                dbgprintf(0, "Using user-defined parameter restraining with "
904                                  "wrong number of parameters:\n"
905                                  "actual nparam = %d, user restraining nparam = %d\n"
906                                  "Defaulting to standard tests\n", nparam, restrain_nparam);
907                return check_ecf_params(param, nparam, fitfunc);
908        }
909
910
911        for (i=0; i<nparam; i++) {
912                if (restrain_restraining[i]) {
913                        if (param[i] < restrain_minval[i])
914                                param[i] = restrain_minval[i];
915                        else if (param[i] > restrain_maxval[i])
916                                param[i] = restrain_maxval[i];
917                }
918        }
919
920
921        return 0;
922}
923
924
925/********************************************************************
926
927                                                          MISCELLANY
928
929*********************************************************************/
930
931
932/* From filters */
933int ECF_Find_Float_Max (float data[], int np, float *max_val)
934{
935        int i, maxi;
936        float *ptr=data;
937
938        if (np < 1)        /* daft input */
939                return -1;
940
941        maxi = 0;
942        *max_val = *ptr;
943
944        for (i=1; i<np; i++)
945                if (*++ptr>*max_val) {
946                        *max_val = *ptr;
947                        maxi = i;
948                }
949
950        return maxi;
951}
952
953
954/* Debugging version of printf */
955#ifdef _CVI_
956int dbgprintf(int dbg_level, const char *format, ...)
957{
958        int ret = 0;
959        char msg[1000];
960        va_list va;
961
962        if (ECF_debug >= dbg_level) {
963                va_start(va, format);
964                ret = vsprintf(msg, format, va);  // Here's hoping, no vsnprintf ...
965                MessagePopup("Debug info", msg);
966                va_end(va);
967        }
968
969        return ret;
970}
971#else
972int dbgprintf(int dbg_level, const char *format, ...)
973{
974        int ret = 0;
975        va_list va;
976
977        if (ECF_debug >= dbg_level) {
978                va_start(va, format);
979                ret = vfprintf(stderr, format, va);
980                va_end(va);
981        }
982
983        return ret;
984}
985#endif
986
987#ifdef TESTCHISQ
988int main(int ac, char **av)
989{
990        int i, nu, ret;
991        int debug = 0;
992        float root, chisq;
993        float percents[4] = { 0.50, 0.68, 0.90, 0.95 };
994
995        for (i=0; i<4; i++) {
996                chisq = percents[i];
997                if (debug)
998                        fprintf(stderr, "Finding x for chisq=%.2f:\n", chisq);
999                printf("static float chisq%d[MAXFIT+1] = { 0", (int)(100*chisq + 0.1));
1000                for (nu=1; nu<=20; nu++) {
1001                        ret = GCI_chisq(nu, chisq, &root);
1002                        if (ret != 0)
1003                                fprintf(stderr, "error %d with nu=%d, chisq=%.2f\n",
1004                                                ret, nu, chisq);
1005                        if (debug)
1006                                fprintf(stderr, "nu=%2d  x=%6.2f  chisq(%d,%.3f)=%6.4f\n",
1007                                                nu, root, nu, root, GCI_gammp(0.5*nu, 0.5*root));
1008                        printf(",%s%.2f", nu % 5 == 0 ? "\n                                   " : " ", root);
1009                }
1010                if (debug) fprintf(stderr, "\n");
1011                printf(" };\n");
1012        }
1013}
1014#endif
1015
1016//***************************************** ExportParams ***********************************************/
1017
1018void ECF_ExportParams_start (char path[])
1019{
1020        ecf_exportParams = 1;
1021        if (path) strcpy(ecf_exportParams_path, path);
1022}
1023
1024void ECF_ExportParams_stop (void)
1025{
1026        ecf_exportParams = 0;
1027}
1028
1029static FILE *ecf_exportFileStream;
1030
1031void ecf_ExportParams_OpenFile (void)
1032{
1033        ecf_exportFileStream = fopen(ecf_exportParams_path, "a");
1034}
1035
1036void ecf_ExportParams_CloseFile (void)
1037{
1038        fprintf(ecf_exportFileStream, "\n");
1039        fclose(ecf_exportFileStream);
1040}
1041
1042void ecf_ExportParams (float param[], int nparam, float chisq)
1043{
1044        int i;
1045
1046        for (i=0; i<nparam; i++) fprintf(ecf_exportFileStream, "%g, ", param[i]);
1047        fprintf(ecf_exportFileStream, "%g\n", chisq);
1048}
1049
1050// Emacs settings:
1051// Local variables:
1052// mode: c
1053// c-basic-offset: 4
1054// tab-width: 4
1055// End:
Note: See TracBrowser for help on using the repository browser.