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

Revision 7803, 28.0 KB checked in by paulbarber, 8 years ago (diff)

Fixed some memory leaks and some compile warnings.

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