/* --------------------------------------------------- This is a general purpose multiparameter fit. This file contains all routines needed to do the whole calculations, except for the Fitfunctions. The parameterfile which must be supplied after the "-t" option has the following format: 1 functionname ; discussed later 2 epsilon maxiter dchi2 ; default: 0.0003 10 1.0 3 first ; first channel to be fitted 4 last ; last channel to be fitted 5 1st [+- error] (x) ; first parameter (const back.) 6 2nd [+- error] (x) ; second parameter (linear back.) 7 3rd [+- error] (x) ; third parameter (quad. back.) 8 4th [+- error] (x) ; fourth parameter (?) . . . The optional "x" or "X" identifies the parameters which are to be fitted. All other parameters are treated as constants. Any number behind the keyword "+-" is ignored on input and is used used as output only. As fit result, a similar table is printed, which can be used as input again. You may also specify an outputfile for the fitted curve. Fit Functions: Fit functions are implemented as external programms which will receive the following parameters: 1. address of array (double) of x-values. (or shared memory ID in UNIX) 2. address of array (double) of results. (or shared memory ID in UNIX) 3. first array element 4. last array element 5. size of array in bytes (used for shmat() under UNIX !) 6. 1st parameter (double) 7. 2nd parameter (double) . . . 26. last parameter (if used) On Amiga Computers the Address represents the real memory location in decimal format (%d) On UNIX Computers this parameter identifies a shared memory ID. See "examplefn.c" On Amiga Computers it is recommended to make the Fit function resident, since it will be called very often. --------------------------------------------------- */ #include #include #include #include #ifdef UNIX #include #include #include #include #endif #ifdef AMIGA #include #endif #ifndef RAND_MAX #define RAND_MAX 2147483647 #endif #define sqr(x) ((x)*(x)) #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a)) #define EPS 0.0003 /* step for differential calculation */ #define MAXITER 10 #define DCHI2 0.01 #define MAXPARM 30 int flag_vv; int *pfree; double *p; double *dfit; double eps = EPS, dchi2 = DCHI2; int maxiter = MAXITER; float *fspc, *ferr, *ftim; double *spc, *err, *tim; double *curv; double *shmx, *shmy; int shmidx, shmidy; int first, last; char *ExtFnName; double fit_chi2; double divide(a,b) double a, b; { double erg; erg = 1.0E10; if( b != 0.0) erg = a / b; return(erg); } /* -------------------------------------------------------------- MATRIX functions -------------------------------------------------------------- */ double invmat(a, b, n, n1) double *a, *b; int n, n1; { /* purpose: to solve an equation system and inverse a matrix parameters: a - on input the coefficient matrix in the eq.sys. a*x=b on output the inverse of the input matrix a b - vector of length n in eq.sys: a*x=b on input on output the roots to the eq.sys. n - number of unknown variables n1 first dimension of the matrix a in the calling program a[n1][n] invmat returns the determinant of a if returned with 0 no solution is given translated from fortran to pascal nov84/Bengt Lindgren translated from pascal to C feb91/Bengt Lindgren */ double *pivot; int *ipivot; int *index; int irow, icolum, i, j, k, l, l1; double amax, swap, t; double determ; if(n>MAXPARM) return(0); /* restricted to max 50*50 matrix */ pivot = (double *) calloc(MAXPARM,sizeof(double)); ipivot = (int *) calloc(MAXPARM,sizeof(int)); index = (int *) calloc(2 * MAXPARM,sizeof(int)); determ = 1.0; for (j = 0; j < n; j++) ipivot[j] = 0; for (i = 0; i < n; i++) { amax = 0; for (j = 0; j < n; j++) if (ipivot[j] != 1) for (k = 0; k < n; k++) { if (ipivot[k] > 1) { free(index); free(ipivot); free(pivot); if(fabs(determ) < 1.0E33) return((double) determ); return((double)1.0E33); } else if ((ipivot[k] < 1) && (fabs(amax) < fabs(a[j*n1+k]))) { irow = j; icolum = k; amax = a[j*n1+k]; } } if (amax == 0) { free(index); free(ipivot); free(pivot); return(0); } ipivot[icolum] = ipivot[icolum] + 1; if (irow != icolum) { determ = -determ; for (l = 0; l < n; l++) { swap = a[irow*n1+l]; a[irow*n1+l] = a[icolum*n1+l]; a[icolum*n1+l] = swap; } swap = b[irow]; b[irow] = b[icolum]; b[icolum] = swap; } index[i] = irow; index[MAXPARM + i] = icolum; pivot[i] = a[icolum*n1+icolum]; determ = determ * (double) pivot[i]; a[icolum*n1+icolum] = 1; for (l = 0; l < n; l++) a[icolum*n1+l] = a[icolum*n1+l] / pivot[i]; b[icolum] = b[icolum] / pivot[i]; for (l1 = 0; l1 < n; l1++) if (l1 != icolum) { t = a[l1*n1+icolum]; a[l1*n1+icolum] = 0; for (l = 0; l < n; l++) a[l1*n1+l] = a[l1*n1+l] - a[icolum*n1+l] * t; b[l1] = b[l1] - b[icolum] * t; } } for (l=n-1; l>=0; l--) { if (index[l] != index[MAXPARM + l]) { irow = index[l]; icolum = index[MAXPARM + l]; for (k = 0; k < n; k++) { swap = a[k*n1+irow]; a[k*n1+irow] = a[k*n1+icolum]; a[k*n1+icolum] = swap; } } } free(index); free(ipivot); free(pivot); return((double)determ); } /* invmat */ help() { printf("fit [Spectrum] -t ParameterTable [-f ResultTable]\n"); printf(" [-o Curve] [-c tica] [-e error] [-v]\n"); printf("parameters:\n"); printf(" -o specify output spectrum name for fit curve\n"); printf(" -t specify file containing parameter table\n"); printf(" -f specify file for the resulting parameter table\n"); printf(" -c generate theory curve, do not fit spectrum, use tica\n"); printf(" -x starting x-value for theory curve\n"); printf(" -e use absolute error for theory curve\n"); printf(" -v verbose: print fit progress\n"); printf("\n"); printf(" The parameterfile which must be supplied after\n"); printf(" the '-t' option has the following format:\n"); printf("\n"); printf(" 1 functionname ; discussed later\n"); printf(" 2 epsilon maxiter dchi2 ; default: 0.0003 10 1.0\n"); printf(" 3 first ; first channel to be fitted\n"); printf(" 4 last ; last channel to be fitted\n"); printf(" 5 1st [+- error] (x) ; first parameter (const back.)\n"); printf(" 6 2nd [+- error] (x) ; second parameter (linear back.)\n"); printf(" 7 3rd [+- error] (x) ; third parameter (quad. back.)\n"); printf(" 8 4th [+- error] (x) ; fourth parameter (?)\n"); printf(" .\n"); printf(" .\n"); printf(" .\n"); printf("\n"); printf(" The optional 'x' or 'X' identifies the parameters\n"); printf(" which are to be fitted. All other parameters are\n"); printf(" treated as constants.\n"); printf(" Any number behind the keyword '+-' is ignored on input\n"); printf(" and is used as output only.\n"); printf(" As fit result, a similar table is printed,\n"); printf(" which can be used as input again.\n"); printf(" You may also specify an outputfile for the fitted\n"); printf(" curve using the '-o' option.\n"); printf("\n"); printf(" Fit Functions:\n"); printf(" Fit functions are implemented as external programms\n"); printf(" which will receive the following parameters:\n"); printf(" 1. address of an array (double) of x-values.\n"); printf(" (or shared memory ID in UNIX)\n"); printf(" 2. address of array (double) of results.\n"); printf(" (or shared memory ID in UNIX)\n"); printf(" 3. first array element\n"); printf(" 4. last array element\n"); printf(" 5. size of array in bytes (used for shmat() under UNIX !)\n"); printf(" 6. 1st parameter (double)\n"); printf(" 7. 2nd parameter (double)\n"); printf(" .\n"); printf(" .\n"); printf(" .\n"); printf(" 26. last parameter (if used)\n"); printf(" On Amiga Computers the Address represents the\n"); printf(" real memory location in decimal format (%d)\n"); printf(" On UNIX Computers this parameter identifies\n"); printf(" a shared memory ID. See 'examplefn.c'\n"); printf(" On Amiga Computers it is recommended to make\n"); printf(" the Fit function resident, since it will be\n"); printf(" called very often.\n"); printf("\n The Fit routine was taken from Bengt Lindgren\n"); exit(0); } fit(first,last) int first,last; { double *fval,*dval; double *b, *h, *df; int *ifit; int i, k, q, iter, freeparm, nfit; double s, dd, w, tmp; fval = (double *) calloc(_MAXSPCLEN,sizeof(double)); b = (double *) calloc(MAXPARM,sizeof(double)); df = (double *) calloc(MAXPARM,sizeof(double)); h = (double *) calloc(sqr(MAXPARM),sizeof(double *)); ifit = (int *) calloc(MAXPARM,sizeof(int)); nfit = 0; /* determine number of parameters to be fitted */ for(i = 0; i < MAXPARM; i++) { if(pfree[i] != 0) ifit[nfit++] = i; } iter = 0; fit_chi2 = 0.0; freeparm = last - first + 1 - nfit; if (freeparm < 1) { free(ifit); free(h); free(df); free(b); free(fval); fprintf(stderr,"too few datapoints to fit and too many param.\n"); return(0); } dval = (double *) calloc((nfit * _MAXSPCLEN),sizeof(double)); if(dval == NULL) { printf("could not allocate enough memory !\n"); exit(-1); } for (i = first; i <= last; i++) { shmx[i] = tim[i]; } do { for (k = 0; k < nfit; k++) { b[k] = 0; for (q = 0; q < nfit; q++) h[(k * MAXPARM) + q] = 0; } call_ext_fn(first,last); for (i = first; i <= last; i++) { fval[i] = shmy[i]; w = divide(1.0,sqr(err[i])); if (iter == 0) fit_chi2 += w * sqr(fval[i] - spc[i]); } for (k = 0; k < nfit; k++) { if(fabs(p[ifit[k]]) > 1) dd = eps * p[ifit[k]]; else dd = eps; tmp = p[ifit[k]]; p[ifit[k]] = tmp + dd; /* assume all p's are of the order 1 */ call_ext_fn(first,last); p[ifit[k]] = tmp; for(i = first; i <= last; i++) dval[(k * _MAXSPCLEN) + i] = shmy[i]; } for (i = first; i <= last; i++) { for (k = 0; k < nfit; k++) { /* calculate derivatives */ if(fabs(p[ifit[k]]) > 1) dd=eps*p[ifit[k]]; else dd=eps; df[k] = (dval[(k * _MAXSPCLEN) + i] - fval[i]) / dd; } for (k = 0; k < nfit; k++) { w = divide(1.0,sqr(err[i])); b[k] += w * df[k] * (spc[i] - fval[i]); for (q = 0; q < nfit; q++) h[(k * MAXPARM) + q] += w * df[k] * df[q]; } } if (iter == 0) { fit_chi2 = fit_chi2 / freeparm; if(flag_vv) printf("iter %d chisq= %10.5f\n",iter,fit_chi2); } dd=invmat(h,b,nfit,MAXPARM); if (dd == 0) { if(flag_vv) printf("*** error determinant =0 ***\n"); break; } for (k = 0; k < nfit; k++) p[ifit[k]] += b[k]; /* b contains the solution */ iter = iter + 1; s = 0.0; for (i = first; i <= last; i++) s = s + sqr(divide(fval[i] - spc[i],err[i])); s = s / (double) freeparm; if(flag_vv) printf("iter %3d chisq= %10.5f\n",iter,s); for (k = 0; k < nfit; k++) { q = ifit[k]; dfit[q] = sqrt(fabs(h[(k * MAXPARM) + k]) * s); if(flag_vv) printf("%2d %8.4f %8.4f\n",q,p[q],dfit[q]); } dd = fabs(s - fit_chi2); if(iter == 1) dd = 100.0; fit_chi2 = s; } while((dd > dchi2) && (iter < maxiter)); free(ifit); free(h); free(df); free(b); free(fval); } /* fit */ /* calculates simulated spc from function */ void simula(tbeg,tinc,chan,df) double tbeg; double tinc; int chan; double df; { int i; double f,rm; rm = (double) RAND_MAX; for (i = 0; i <= chan; i++) { tim[i] = tbeg + (i * tinc); shmx[i] = tim[i]; } call_ext_fn(0,chan); for (i = 0; i <= chan; i++) { f = ((2.0 * ((double) rand()/rm)) - 1.0) * df; err[i] = df; spc[i] = shmy[i] + f; } } /* simula */ main(argc,argv) int argc; char **argv; { int n,max; double x,tica; double df; char *s; FILE *fp; ExtFnName = (char *) malloc(80); s = (char *) malloc(80); p = (double *) calloc(MAXPARM,sizeof(double)); dfit = (double *) calloc(MAXPARM,sizeof(double)); pfree = (int *) calloc(MAXPARM,sizeof(int)); spc = (double *) calloc(_MAXSPCLEN,sizeof(double)); err = (double *) calloc(_MAXSPCLEN,sizeof(double)); tim = (double *) calloc(_MAXSPCLEN,sizeof(double)); fspc = (float *) calloc(_MAXSPCLEN,sizeof(float)); ferr = (float *) calloc(_MAXSPCLEN,sizeof(float)); ftim = (float *) calloc(_MAXSPCLEN,sizeof(float)); curv = (double *) calloc(_MAXSPCLEN,sizeof(double)); alloc_shm(); flag_vv = 0; df = 0.0; if(checkopt(argc,argv,"-t",s)) { read_table(s); } else { fprintf(stderr,"you MUST specify an input table(-t option)\n"); exit(-1); } if(checkopt(argc,argv,"-v",s)) flag_vv = 1; if(checkopt(argc,argv,"-e",s)) df = atosf(s); x = 0.0; if(checkopt(argc,argv,"-x",s)) x = atosf(s); if(checkopt(argc,argv,"-c",s)) { tica = atosf(s); simula(x,tica,last,df); max = last; } else { strcpy(s,argv[1]); fp=fopen(s,"r"); if(fp == NULL) { strcat(s,".spc"); fp = fopen(s,"r"); if(fp == NULL) { printf("spectrum not found %s\n",argv[1]); exit(0); } } fclose(fp); max = readspec(argv[1],fspc,ferr,ftim,s); for(n = 0; n < max; n++) { spc[n] = (double) fspc[n]; err[n] = (double) ferr[n]; tim[n] = (double) ftim[n]; } if(last > max) last = max; fit(first,last-1); tica = tim[2] - tim[1]; simula(tim[0],tica,max,0.0); if((flag_vv != 0)||(checkopt(argc,argv,"-f",s) == 0)) write_table(""); } if(checkopt(argc,argv,"-o",s)) { for(n = 0; n < max; n++) { fspc[n] = (float) spc[n]; ferr[n] = (float) err[n]; ftim[n] = (float) tim[n]; } writespec(s,fspc,ferr,max,2,"fitcurve"); strcat(s,".tim"); writespec(s,ftim,ferr,max,2,"fitcurve"); } if(checkopt(argc,argv,"-f",s)) write_table(s); free_shm(); free(curv); free(ftim); free(ferr); free(fspc); free(tim); free(err); free(spc); free(pfree); free(dfit); free(p); free(s); free(ExtFnName); exit(0); } /* ---------------------------------------------------------------- read parameter table ---------------------------------------------------------------- */ read_table(table) char *table; { int i,j,n; double la,lb; char *s, *z; FILE *fp; s = (char *) malloc(128); z = (char *) malloc(128); fp = fopen(table,"r"); if(fp == NULL) { fprintf(stderr,"could not open parameter table >%s<\n",s); exit(-1); } read_line(s,fp); strcpy(ExtFnName,s); /* read name of fit function */ read_line(s,fp); sscanf(s,"%lf %d %lf",&la,&n,&lb); if(la > 0.0) eps = (double) la; if(n > 0) maxiter = n; if(lb > 0.0) dchi2 = (double) lb; read_line(s,fp); n = atoi(s); first = 0; if(n > 0) first = n; read_line(s,fp); n = atoi(s); last = _MAXSPCLEN; if(n > 0) last = n; for(n = 0; n < MAXPARM; n++) pfree[n] = 0; for(n = 0; n < MAXPARM; n++) { read_line(s,fp); if(feof(fp)) break; i = 0; j = 0; while(s[i] != ' ') z[j++] = s[i++]; z[j] = 0; p[n] = atosf(z); while(s[i] != 0) { if(toupper(s[i]) == 'X') { pfree[n] = 1; if(toupper(s[i + 1]) == '0') pfree[n] = 0; } i = i + 1; } } free(z); free(s); fclose(fp); return(0); } read_line(s,fp) char *s; FILE *fp; { int i,n; char *z; z = (char *) malloc(128); fgets(s,128,fp); while(s[0] == ';') { /* skip comment lines */ if(feof(fp)) { strcpy(s,""); return(0); } fgets(s,128,fp); } i = 0; n = 0; while(s[i] == ' ') i++; while((s[i] != ';') && (s[i] != '\n') && (s[i] != 0)) z[n++] = s[i++]; /* break at comment */ z[n] = 0; strcpy(s,z); free(z); } write_table(table) char *table; { int n; char c; FILE *fp; if(table[0] != 0) { fp = fopen(table,"w"); if(fp == NULL) { fprintf(stderr,"could not open file for output table\n"); fp = stdout; } } else { fp = stdout; } fprintf(fp,"; ----- Fit result -----\n"); fprintf(fp,"; chi^2 = %f\n",fit_chi2); fprintf(fp,"; ----------------------\n"); fprintf(fp,"%s\n",ExtFnName); fprintf(fp,"%f %d %f\n",eps,maxiter,dchi2); fprintf(fp,"%d ; fisrt channel\n%d ; last channel\n",first,last); for(n = 0; n < MAXPARM; n++) { c = ' '; if(pfree[n] != 0) c = 'X'; fprintf(fp,"%f +- %f %c ; parameter %d\n",p[n],dfit[n],c,n); } if(fp != stdout) fclose(fp); return(0); } /* ---------------------------------------------------------------- shared memory allocation ---------------------------------------------------------------- */ alloc_shm() { #ifdef AMIGA shmx = (double *) calloc(_MAXSPCLEN,sizeof(double)); shmy = (double *) calloc(_MAXSPCLEN,sizeof(double)); #endif #ifdef UNIX shmidx = shmget(IPC_PRIVATE,_MAXSPCLEN*sizeof(double),0666); shmidy = shmget(IPC_PRIVATE,_MAXSPCLEN*sizeof(double),0666); if((shmidx < 0) || (shmidy < 0)) { fprintf(stderr,"could not get Shared memory ID, aborting !\n"); exit(-1); } shmx = (double *) shmat(shmidx, 0,0666); shmy = (double *) shmat(shmidy, 0,0666); if((shmx < 0) || (shmy < 0)) { fprintf(stderr,"could not attach Shared memory, aborting !\n"); exit(-1); } signal(SIGFPE,SIG_IGN); /* ignore doubleing Point Exceptions ! */ #endif return(0); } free_shm() { #ifdef UNIX shmdt(shmx); shmctl(shmidx,IPC_RMID,0); shmdt(shmy); shmctl(shmidy,IPC_RMID,0); #endif #ifdef AMIGA free(shmx); free(shmy); #endif return(0); } /* ---------------------------------------------------------------- call external fit function ---------------------------------------------------------------- */ call_ext_fn(first,last) int first,last; { char *s, *fnstr; int n; fnstr = (char *) malloc(2048); s = (char *) malloc(80); #ifdef UNIX sprintf(fnstr,"%s %d %d %d %d %d",ExtFnName,shmidx,shmidy, first,last,_MAXSPCLEN*sizeof(double)); #endif #ifdef AMIGA sprintf(fnstr,"%s %d %d %d %d %d",ExtFnName,shmx,shmy, first,last,_MAXSPCLEN*sizeof(double)); #endif for(n = 0; n < MAXPARM; n++) { sprintf(s," %f",p[n]); strcat(fnstr,s); } #ifdef UNIX system(fnstr); #endif #ifdef AMIGA SysCall(fnstr); #endif free(s); free(fnstr); return(0); } #ifdef AMIGA SysCall(str) char *str; { struct FileHandle *input, *output; input = (struct FileHandle *) Open("nil:",MODE_OLDFILE); output = (struct FileHandle *) Open("nil:",MODE_NEWFILE); Execute(str,input,output); Close(input); Close(output); return(0); } #endif