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

Revision 7706, 27.3 KB checked in by aivar, 9 years ago (diff)

Changed GCI_gauss_jordan to be GCI_solve and GCI_invert. The current default CGI_solve uses Gaussian Elimination. Added inversion method. Added Lower/Upper Decomposition code for both solve and invert, but it is commented out.

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