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

Revision 7779, 28.6 KB checked in by aivar, 8 years ago (diff)

Go back to using Gaussian method to solve and invert matrices.

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