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

Revision 7710, 28.4 KB checked in by aivar, 8 years ago (diff)

Work in progress.

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