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

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

Added GPL license to all source code.

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