home *** CD-ROM | disk | FTP | other *** search
/ Dr. CD ROM (Annual Premium Edition) / premium.zip / premium / IBMOS2_1 / LSQRFT13.ZIP / mrqfit2.c < prev    next >
Text File  |  1993-07-10  |  7KB  |  282 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include "fit.h"
  5.  
  6. /* calculates alpha, beta, and chisq */
  7. /* see Numerical Recipes in C for definition of these */
  8. /* the Levenberg Marquardt algorithm requires us to slove */
  9. /* the equation alpha*da = beta (where alpha is a matrix */
  10. /* and da and beta are vectors) for da */
  11. /* this function does this by gauss-jordan */
  12. /* elimination */
  13.  
  14. /* declarations for data declared externally */
  15. extern int grflag;
  16. extern int gnuopen;
  17. extern FILE *pstream;
  18. extern double *xmin, *xmax;
  19. extern int wiflag;
  20. extern int veflag;
  21. extern int debug;
  22. extern char *GNUPLOT;
  23.  
  24. int window(double *x, int num_indep);
  25.  
  26. int mrqfit(data,order,num_indep,ndata,itmax,a,ma,lista,mfit,
  27.                 covar,chisq,func,filename,comment)
  28. double **data;   
  29. struct data_order order;
  30. int num_indep;
  31. int ndata;
  32. int itmax;
  33. double *a;
  34. int ma, *lista, mfit;
  35. double **covar, *chisq;
  36. int (*func)();
  37. char *filename;
  38. char *comment;
  39. {
  40. int i, j, k;  /* indices for loops */
  41. int failed;
  42. double **alpha;
  43. /* alamda determines how much we change the fitting parameters */
  44. double alamda;
  45. char cmd[20]="";
  46. double *beta;
  47. double *atry; /* parameters to try and see if we reduce chisqr */
  48. double **alpha_try;
  49. double ochisq; /* best chisqr so far */
  50. double *da;
  51. FILE *stream;
  52.  
  53. /* allocate space for arrays */
  54. atry = dvector(ma);
  55. alpha_try=dmatrix(ma,ma);
  56. beta = dvector(ma);
  57. da = dvector(ma);
  58. alpha=dmatrix(ma,ma);
  59.  
  60. stream=NULL;
  61.  
  62. /* calculate alpha, beta, and chisq for current value of parameters */
  63. if(debug) printf("in mrqfit, calling alpha_beta_chisqr()\n");
  64. failed=alpha_beta_chisq(data,order,num_indep,ndata,a,ma,lista,mfit,
  65.                             alpha,beta,chisq,func);
  66. if(debug) printf("in mrqfit, alpha_beta_chisqr() returned %d\n", failed);
  67. if(failed){
  68.     free_dmatrix(alpha,ma,ma);
  69.     free(atry);
  70.     free_dmatrix(alpha_try,ma,ma);
  71.     free(beta);
  72.     free(da);
  73.     return 0;
  74. }
  75.  
  76. /* start alamda small */
  77. alamda = 1e-3;
  78. i = 0;
  79.  
  80. /* loop until user tells us to quit */
  81. while( strcmp(cmd,"q") != 0){
  82.  
  83. /* loop specified number of iterations */
  84. while( i <= itmax ){
  85.     ochisq = *chisq;
  86.     for(j=0;j<mfit;j++){
  87.         for(k=0;k<mfit;k++){
  88.             alpha_try[j][k] = alpha[j][k];
  89.         }
  90.         alpha_try[j][j] = alpha[j][j] + alamda;
  91.     }
  92.     if(debug)printf("in mrqfit, calling solve_for_da()\n");
  93.     solve_for_da(alpha_try, covar, beta, da, mfit);
  94.     if(debug)printf("in mrqfit, returned from solve_for_da()\n");
  95.  
  96.     for(j=0;j<ma;j++) atry[j] = a[j];
  97.     for(j = 0; j < mfit; j++) atry[lista[j]] += da[j];
  98.  
  99. if(veflag==2){
  100.         printf("\nalpha: ");
  101.         for(j=0;j<mfit;j++){
  102.             printf("\n %d   ",lista[j]);
  103.             for(k=0;k<mfit;k++) printf(" %g ",alpha_try[j][k]);
  104.         }
  105.         printf("\ncovar: ");
  106.         for(j=0;j<mfit;j++){
  107.             printf("\n %d   ",lista[j]);
  108.             for(k=0;k<mfit;k++) printf("\n covar[%d][%d]= %g ",lista[j], lista[k], covar[j][k]);
  109.         }
  110.         printf("\natry:");
  111.         for(j=0;j<mfit;j++) printf("\n atry%d = %g",lista[j], atry[lista[j]]);
  112.         printf("\n");
  113.     }
  114.  
  115.     if(debug)printf("in mrqfit, calling alpha_beta_chisqr()\n");
  116.     failed=alpha_beta_chisq(data,order,num_indep,ndata,atry,ma,
  117.                             lista,mfit,alpha,beta,chisq,func);
  118.     if(debug)printf("in mrqfit, alpha_beta_chisqr() returned %d\n", failed);
  119.  
  120. if(failed){
  121.     free_dmatrix(alpha,ma,ma);
  122.     free(atry);
  123.     free_dmatrix(alpha_try,ma,ma);
  124.     free(beta);
  125.     free(da);
  126.     return 0;
  127. }
  128.  
  129.     if(*chisq >= ochisq){
  130.         alamda *= 10;
  131.         if(alamda > 1e15) alamda /= 3e14;
  132.         *chisq = ochisq;
  133.     }
  134.     else{
  135.         if(debug){
  136.             if(stream == NULL){
  137.                 printf("opening lamda.sts\n");
  138.                 stream = fopen("lamda.sts","at");
  139.             }
  140.             if(stream != NULL){
  141.                 fprintf(stream,"%g\n",alamda);
  142.                 printf("alamda: %g\n",alamda);
  143.             }
  144.         }
  145.         alamda *= 0.1;
  146.         if(alamda < 1e-15) alamda *= 3e14;
  147.         for(j=0;j < mfit; j++)a[lista[j]] = atry[lista[j]];
  148.     }
  149.  
  150.     if(veflag){
  151.         for(j=0; j < ma; j++)printf("\n a%d = %g", j, a[j]);
  152.         printf("\n chisqr = %g", *chisq);
  153.         printf("\n alamda = %g", alamda);
  154.         printf("\n");
  155.     }
  156.     else printf("iteration: %d chisqr: %g\n", i, *chisq);
  157.     i++;
  158. }
  159.  
  160. if(grflag){
  161.     if(debug)printf("in mrqfit, calling plot()\n");
  162.     failed=plot(func,data, order, num_indep, ndata, filename, comment, a, ma);
  163.     if(debug)printf("in mrqfit, plot() returned %d\n", failed);
  164. }
  165.  
  166. printf("q quit, g graphs, n does not graph, anything else continues ");
  167. gets(cmd);
  168. i = 0;
  169. if(strcmp(cmd,"g") == 0) grflag = 1;
  170. if(strcmp(cmd,"n") == 0) grflag = 0;
  171.  
  172. }
  173. free_dmatrix(alpha,ma,ma);
  174. free(atry);
  175. free_dmatrix(alpha_try,ma,ma);
  176. free(beta);
  177. free(da);
  178. if(stream != NULL) fclose(stream);
  179.  
  180. return 0;
  181. }
  182.  
  183.  
  184. int alpha_beta_chisq(data, order,num_indep,ndata,a,ma,
  185.                                 lista,mfit,alpha,beta,chisq,funcs)
  186. double **data;
  187. struct data_order order;
  188. int num_indep;
  189. int ndata;
  190. double *a;
  191. int ma, *lista,mfit;
  192. double **alpha, *beta, *chisq;
  193. int (*funcs)(double *,double *,double *, double *,int, int, 
  194.                     int *, int *, double *, double);
  195. {
  196. int i,j,k;
  197. int failed;
  198. double *x, *y, *sig,*xsig;
  199. double ymod, *dyda, sig2i, wt,dy;
  200. int *fita, *dydx_flag;
  201. double *dydx;
  202. double sigi;
  203.  
  204. dyda=dvector(ma);
  205. fita = ivector(ma);
  206. x = dvector(num_indep);
  207. dydx_flag = ivector(num_indep);
  208. dydx = dvector(num_indep);
  209.  
  210. y = data[order.y];
  211. sig = data[order.sig];
  212.  
  213. for(i = 0; i < num_indep; i++){
  214.     if(order.xsig[i] >= 0){
  215.         xsig = data[order.xsig[i]];
  216.         dydx_flag[i] = 1;
  217.     }
  218. }
  219.  
  220. for(i = 0; i < ma; i++){
  221.     fita[i] = 0;
  222.     for(j = 0; j < mfit; j++) if(lista[j] == i) fita[i] = 1;
  223. }
  224.  
  225. *chisq = 0;
  226. for(j=0; j < mfit;j++){
  227.     beta[j] = 0;
  228.     for(k=0;k<mfit;k++) alpha[j][k] = 0;
  229. }
  230.  
  231. for(i=0;i<ndata;i++){
  232.     for(j = 0; j < num_indep; j++){
  233.         x[j] = data[order.x[j]][i];
  234.     }
  235.  
  236.     if(wiflag == 0 || window(x, num_indep)){
  237.  
  238.         failed = (*funcs)(x,a,&ymod,dyda,ma,1, fita, dydx_flag, dydx, y[i]);
  239.         if(failed){
  240.         free(dyda);
  241.         free(fita);
  242.         free(x);
  243.         free(dydx_flag);
  244.         free(dydx);
  245.         return 1;
  246.         }
  247.         sig2i = (sig[i]*sig[i]);
  248.         for(j = 0; j < num_indep; j++)
  249.             if(order.xsig[j] >= 0) sig2i += (dydx[j])*(dydx[j])*xsig[i]*xsig[i];
  250.         if(debug==2) printf("i: %d x[0]: %g ymod: %g y: %g sig2i: %g\n", i, x[0], ymod, y[i], sig2i);
  251.         dy = y[i] - ymod;
  252.         for(j=0;j<mfit;j++){
  253.             wt=dyda[lista[j]]/sig2i;
  254.             for(k=0;k<=j;k++)alpha[j][k] += wt*dyda[lista[k]];
  255.             beta[j] += dy*wt;
  256.         }
  257.         *chisq += dy*dy/sig2i;
  258.     }
  259. }
  260.  
  261. for(j=1;j<mfit;j++)
  262.     for(k=0;k<=j-1;k++) alpha[k][j]=alpha[j][k];
  263.  
  264. free(dyda);
  265. free(fita);
  266. free(x);
  267. free(dydx_flag);
  268. free(dydx);
  269. return 0;
  270. }
  271.  
  272. int window(double *x, int num_indep){
  273.  
  274. int i;
  275.  
  276. for(i = 0; i < num_indep; i++)
  277.     if(x[i] < xmin[i] || x[i] > xmax[i]) return 0;
  278.  
  279. return 1;
  280. }
  281.  
  282.