#include /* The 2010 version of the ECF library. This has basically been completely rewritten to avoid license issues. Also, this takes account of the fact that we may be handling Poisson noise. This file contains utility code for various functions needed by the ECF library, and which are not particularly specific to the single or global fitting routines. */ #include #include #include #include #ifdef _CVI_ #include #endif #include "EcfInternal.h" /* For #defines */ int ECF_debug = 0; /******************************************************************** GAUSS-JORDAN AND COVAR-SORTING ROUTINES *********************************************************************/ /* Linear equation solution by Gauss-Jordan elimination. a[0..n-1][0..n-1] is the input matrix, b[0..n] is input containing the single right-hand side vector. On output, a is replaced by its matrix inverse and b is replaced by the corresponding set of solution vectors. */ #define SWAP(a,b) { temp=(a); (a)=(b); (b)=temp; } /* * Based on LinearEquationSolving from linearEquations.c by Henry Guennadi Levkin. */ int GCI_gauss_jordan(float **a, int n, float *b) { float max; float tmp; int i, j, k, m; // base row of matrix for (k = 0; k < n - 1; ++k) { // search for line with max element max = fabs(a[k][k]); m = k; for (i = k + 1; i < n; ++i) { if (max < fabs(a[i][k])) // row i col k { max = a[i][k]; m = i; } } // permutation of base line (index k) and max element line (index m) if (m != k) { for (i = k; i < n; ++i) { tmp = a[k][i]; a[k][i] = a[m][i]; a[m][i] = tmp; } tmp = b[k]; b[k] = b[m]; b[m] = tmp; } if (0.0 == a[k][k]) { return -2; // singular matrix } // triangulation of matrix with coefficients for (j = k + 1; j < n; ++j) // current row of matrix { tmp = -a[j][k] / a[k][k]; for (i = k; i < n; ++i) { a[j][i] += tmp * a[k][i]; } b[j] += tmp * b[k]; // free member recalculation } } for (k = n - 1; k >= 0; --k) { for (i = k + 1; i < n; ++i) { b[k] -= a[k][i] * b[i]; } b[k] /= a[k][k]; } return 0; } //============================================================================== // return 1 if system not solving // nDim - system dimension // pfMatr - matrix with coefficients // pfVect - vector with free members // pfSolution - vector with system solution // pfMatr becames trianglular after function call // pfVect changes after function call // // Developer: Henry Guennadi Levkin // //============================================================================== int LinearEquationsSolving(int nDim, double* pfMatr, double* pfVect, double* pfSolution) { double fMaxElem; double fAcc; int i , j, k, m; for(k=0; k<(nDim-1); k++) // base row of matrix { // search of line with max element fMaxElem = fabs( pfMatr[k*nDim + k] ); m = k; for(i=k+1; i=0; k--) { pfSolution[k] = pfVect[k]; for(i=(k+1); i=0; j--) { if (paramfree[j]) { for (i=0; i 0) { xa = x / param[2]; /* xa = x/param[2] */ lxa = log(xa); /* lxa = log(x/param[2]) */ xah = exp(lxa / param[3]); /* xah = exp(log(x/param[2])/param[3]) = (x/param[2])^(1/param[3]) */ dy_dparam[1] = ex = exp(-xah); /* ex = exp(-(x/param[2])^(1/param[3])) */ ex *= param[1]; /* ex = param[1] * exp(-(x/param[2])^(1/param[3])) */ *y = ex; /* y is now correct */ ex *= xah / param[3]; /* ex = param[1] * exp(...) * (x/param[2])^(1/param[3]) * 1/param[3] */ dy_dparam[2] = ex / param[2]; dy_dparam[3] = ex * lxa / param[3]; } else if (x > -1e-10) { // Almost zero *y = param[1]; dy_dparam[1] = 1; dy_dparam[2] = dy_dparam[3] = 0; } else { /* Ouch: x<0 */ fprintf(stderr, "Can't have x < 0 in stretched exponential!!\n"); } } /* This is actually essentially the same as the other version; because of the power (x/tau)^(1/h), we cannot use the same tricks as for the multiexponential case, unfortunately. */ int stretchedexp_array(float xincr, float param[], float *y, float **dy_dparam, int nx, int nparam) { int i; float ex, lxa, xah, a2inv, a3inv; double xa, xaincr; if (xincr == 0) return -1; xa=0; xaincr = xincr / param[2]; a2inv = 1/param[2]; a3inv = 1/param[3]; /* When x=0 */ y[0] = param[1]; dy_dparam[0][1] = 1; dy_dparam[0][2] = dy_dparam[0][3] = 0; for (i=1; i MAX_Z) return -21; if (param[1] < MIN_A || param[1] > MAX_A) return -22; if (param[2] < MIN_TAU || param[2] > MAX_TAU) return -23; break; case 5: if (param[0] < MIN_Z || param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3])) || param[0] > MAX_Z) return -21; if (param[1] < MIN_A || param[1] > MAX_A) return -22; if (param[2] < MIN_TAU || param[2] > MAX_TAU) return -23; if (param[3] < MIN_A || param[3] > MAX_A) return -24; if (param[4] < MIN_TAU || param[4] > MAX_TAU) return -25; break; case 7: if (param[0] < MIN_Z || param[0] < -MIN_Z_FACTOR * (fabs(param[1]) + fabs(param[3]) + fabs(param[5])) || param[0] > MAX_Z) return -21; if (param[1] < MIN_A || param[1] > MAX_A) return -22; if (param[2] < MIN_TAU || param[2] > MAX_TAU) return -23; if (param[3] < MIN_A || param[3] > MAX_A) return -24; if (param[4] < MIN_TAU || param[4] > MAX_TAU) return -25; if (param[5] < MIN_A || param[5] > MAX_A) return -26; if (param[6] < MIN_TAU || param[6] > MAX_TAU) return -27; break; } } else if (fitfunc == GCI_stretchedexp) { if (param[0] < MIN_Z || param[0] < -MIN_Z_FACTOR * fabs(param[1]) || param[0] > MAX_Z) return -21; if (param[1] < MIN_A || param[1] > MAX_A) return -22; if (param[2] < MIN_TAU || param[2] > MAX_TAU) return -23; if (param[3] < MIN_H || param[3] > MAX_H) return -24; } return 0; } /* For the user-specified version, we have some global variables to store the settings between invocations. */ static int restrain_nparam=0; /* How many parameters have we set up? */ static int restrain_restraining[MAXFIT]; /* Do we check parameter i? */ static float restrain_minval[MAXFIT]; /* Minimum acceptable parameter value */ static float restrain_maxval[MAXFIT]; /* Maximum acceptable parameter value */ int GCI_set_restrain_limits(int nparam, int restrain[], float minval[], float maxval[]) { int i; if (nparam < 0 || nparam > MAXFIT) return -1; /* We're going to be doing something, so clear the memory */ restrain_nparam = 0; for (i=0; i maxval[i]) return -2; restrain_minval[i] = minval[i]; restrain_maxval[i] = maxval[i]; } else restrain_restraining[i] = 0; } restrain_nparam = nparam; return 0; } /* original version int check_ecf_user_params (float param[], int nparam, void (*fitfunc)(float, float [], float *, float [], int)) { int i; if (restrain_nparam != nparam) { dbgprintf(0, "Using user-defined parameter restraining with " "wrong number of parameters:\n" "actual nparam = %d, user restraining nparam = %d\n" "Defaulting to standard tests\n", nparam, restrain_nparam); return check_ecf_params(param, nparam, fitfunc); } for (i=0; i restrain_maxval[i])) return -21; } } return 0; } */ // new version from J Gilbey 31.03.03 int check_ecf_user_params (float param[], int nparam, void (*fitfunc)(float, float [], float *, float [], int)) { int i; if (restrain_nparam != nparam) { dbgprintf(0, "Using user-defined parameter restraining with " "wrong number of parameters:\n" "actual nparam = %d, user restraining nparam = %d\n" "Defaulting to standard tests\n", nparam, restrain_nparam); return check_ecf_params(param, nparam, fitfunc); } for (i=0; i restrain_maxval[i]) param[i] = restrain_maxval[i]; } } return 0; } /******************************************************************** MISCELLANY *********************************************************************/ /* From filters */ int ECF_Find_Float_Max (float data[], int np, float *max_val) { int i, maxi; float *ptr=data; if (np < 1) /* daft input */ return -1; maxi = 0; *max_val = *ptr; for (i=1; i*max_val) { *max_val = *ptr; maxi = i; } return maxi; } /* Debugging version of printf */ #ifdef _CVI_ int dbgprintf(int dbg_level, const char *format, ...) { int ret = 0; char msg[1000]; va_list va; if (ECF_debug >= dbg_level) { va_start(va, format); ret = vsprintf(msg, format, va); // Here's hoping, no vsnprintf ... MessagePopup("Debug info", msg); va_end(va); } return ret; } #else int dbgprintf(int dbg_level, const char *format, ...) { int ret = 0; va_list va; if (ECF_debug >= dbg_level) { va_start(va, format); ret = vfprintf(stderr, format, va); va_end(va); } return ret; } #endif #ifdef TESTCHISQ int main(int ac, char **av) { int i, nu, ret; int debug = 0; float root, chisq; float percents[4] = { 0.50, 0.68, 0.90, 0.95 }; for (i=0; i<4; i++) { chisq = percents[i]; if (debug) fprintf(stderr, "Finding x for chisq=%.2f:\n", chisq); printf("static float chisq%d[MAXFIT+1] = { 0", (int)(100*chisq + 0.1)); for (nu=1; nu<=20; nu++) { ret = GCI_chisq(nu, chisq, &root); if (ret != 0) fprintf(stderr, "error %d with nu=%d, chisq=%.2f\n", ret, nu, chisq); if (debug) fprintf(stderr, "nu=%2d x=%6.2f chisq(%d,%.3f)=%6.4f\n", nu, root, nu, root, GCI_gammp(0.5*nu, 0.5*root)); printf(",%s%.2f", nu % 5 == 0 ? "\n " : " ", root); } if (debug) fprintf(stderr, "\n"); printf(" };\n"); } } #endif //***************************************** ExportParams ***********************************************/ void ECF_ExportParams_start (char path[]) { ecf_exportParams = 1; if (path) strcpy(ecf_exportParams_path, path); } void ECF_ExportParams_stop (void) { ecf_exportParams = 0; } static FILE *ecf_exportFileStream; void ecf_ExportParams_OpenFile (void) { ecf_exportFileStream = fopen(ecf_exportParams_path, "a"); } void ecf_ExportParams_CloseFile (void) { fprintf(ecf_exportFileStream, "\n"); fclose(ecf_exportFileStream); } void ecf_ExportParams (float param[], int nparam, float chisq) { int i; for (i=0; i