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

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

Work in progress. Working on GCI_marquardt_compute_fn, code is commented out so it compiles.

Tidied up a bit.

Unexpected behavior in EcfUtil.c: Optimization looks good in generated assembly code but runs much slower.

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