home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / lsqrft15.zip / linear2.c < prev    next >
Text File  |  1993-08-18  |  7KB  |  221 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include "fit.h"
  5.  
  6. /* declarations for data declared externally */
  7. extern int grflag;
  8. extern int wiflag;
  9. extern int veflag;
  10. extern int debug;
  11. extern int MYPLOT;
  12.  
  13. /* linear_fit() does a linear least squares fit to a specified function. */
  14. /* because the functions are defined for non-linear fitting, we have */
  15. /* to jump through some hoops to do the calculation. */
  16.  
  17.  
  18. int linear_fit(double **data,struct data_order order,
  19.                 int num_indep, int ndata, int itmax, double *a,
  20.                 int ma, int *lista, int mfit, double **covar,
  21.                 double *chisq, int (*func)(),
  22.                 char *filename, char *comment){
  23.  
  24. int i, j, k, l, m;  /* indices for loops */
  25. int failed=0;
  26. double **alpha;     /* See Chapter 14 of Numerical Recipes */
  27. double *beta;
  28. double *X;
  29. double *x;
  30. double *dyda;       /* Not really needed, but function might return */
  31.                     /* values for dyda[i], so we have to allocate space */
  32.  
  33. double *aj;         /* The function needs to return the value of each basis */
  34.                     /* function at a given x, X[j](x[i]).  However, function */
  35.                     /* is set up to return value of f(x,a).   */
  36.                     /* To get X[j](x[i]), we send a parameter list where */
  37.                     /* a[i] = 0 if a != j, and a[i] = 1 if a ==j. */
  38.                     /* I call this parameter list aj to distinguish it from a. */
  39.  
  40. int *fita;          /* fita tells the user function  func() which dyda[i]'s to */
  41.                     /* compute.  We don't need any of them for linear fitting */
  42.                     /* so we set all the fita[i]'s to zero */
  43.  
  44. int *dydx_flag;     /* dydx_flag tells the user function which dydx[i]'s to compute. */
  45.                     /* We can't use them for linear fitting, because we are not */
  46.                     /* iterating and don't have a guess for the parameters until */
  47.                     /* the fit is finished.  As a result, our linear fitting does */
  48.                     /* not consider errors in the independent variables.  We set */
  49.                     /* all of the dydx_flag[i]'s to zero, so the function does not */
  50.                     /* have to waste time computing them. */
  51.  
  52. double *dydx;       /* However, since user functions might compute the dydx[i]'s */
  53.                     /* anyway, we need to allocate storage. */
  54.  
  55. double sig2i;       /* sigmayi * sigmay1 */
  56. double *atry;       /* These are the parameter values multiplying the basis functions */
  57.                     /* used in the fit.  If a basis function is not used in the fit, */
  58.                     /* it is multiplied by zero.  atry has dimension mfit and holds */
  59.                     /* only the parameters multiplying basis functions which are fit */
  60. double chisqr;
  61.  
  62.  
  63. alpha=dmatrix(mfit,mfit);
  64. X = dvector(ma);          /* value of basis function at each data point */
  65. beta = dvector(mfit);
  66. dyda=dvector(ma);
  67. fita = ivector(ma);
  68. x = dvector(num_indep);
  69. dydx = dvector(num_indep);
  70. dydx_flag = ivector(num_indep);
  71. aj = dvector(ma);
  72. atry = dvector(mfit);
  73.  
  74. for(j = 0; j < num_indep; j++) dydx_flag[j] = 0;  /* don't calculate dydx[i]'s */
  75. for(j = 0; j < ma; j++) fita[j] = 0;              /* don't calculate dyda[i]'s */
  76.  
  77.  
  78. /* set all alpha[i][j]'s equal to zero */
  79.  
  80. for(i = 0; i < mfit; i++){
  81.     for(j = 0; j < mfit; j++){
  82.         alpha[i][j] = 0;
  83.     }
  84.     beta[i] = 0;
  85. }
  86.  
  87. /* loop summing over data points */
  88. for(i = 0; i < ndata; i++){
  89.  
  90.     /* set up array of independent variable values at i'th data point */
  91.     for(j = 0; j < num_indep; j++){
  92.         x[j] = data[order.x[j]][i];
  93.         if(debug > 1)
  94.             printf("i: %d order.x[%d]: %d x[%d]: %g\n", i, j, order.x[j], j, x[j]);
  95.     }
  96.  
  97.     /* If windowing is on, only consider data points within window */
  98.     if(wiflag == 0 || window(x, num_indep)){
  99.  
  100.         /* calculate value of relevant basis functions at current data point */
  101.         for(j = 0; j < mfit; j++){
  102.  
  103.             /* set up aj[i]'s so function value = basis function value */
  104.             for(l=0; l < ma; l++){
  105.                 if(l != lista[j]) aj[l] = 0;
  106.                 else aj[l] = 1;
  107.             }
  108.             failed = (*func)(x,aj,&X[lista[j]],dyda,ma,0, fita, dydx_flag, dydx, data[order.y][i]);
  109.             if(debug > 1) printf("X[%d]: %g ",lista[j], X[lista[j]]);
  110.             if(failed) break;
  111.         }
  112.         if(debug > 1) printf("\n");
  113.         if(failed) break;
  114.  
  115.         sig2i = data[order.sig][i]*data[order.sig][i];
  116.         if(debug > 1) printf("order.sig: %d data[order.sig][i]: %g sig2i: %g\n",
  117.                                     order.sig, data[order.sig][i], sig2i);
  118.  
  119.         /* add this data point to alpha and beta */
  120.         for(k = 0; k < mfit; k++){
  121.             for(j = 0; j <= k; j++){
  122. /*                if(debug > 1) printf("alpha[%d][%d] += X[%d]*X[%d]/sig2i;\n",
  123.                     j, k, lista[j], lista[k]); */
  124.                 alpha[j][k] += X[lista[j]]*X[lista[k]]/sig2i;
  125.             }
  126.             beta[k] += data[order.y][i]*X[lista[k]]/sig2i;
  127.         }
  128.     }
  129. }
  130.  
  131. if(debug) printf("done looping over data points\n");
  132.  
  133. if(!failed){
  134.  
  135.     if(debug){
  136.         printf("\nalpha: ");
  137.         for(j=0;j<mfit;j++){
  138.             printf("\n %d   ",j);
  139.             for(k=0;k<mfit;k++) printf(" %g ",alpha[j][k]);
  140.         }
  141.     }
  142.  
  143.     /* fill in rest of alpha, alpha is symmetric, and we only  */
  144.     /*calculated upper triangular, we need the whole thing */
  145.     for(k = mfit-1; k >= 0; k--){
  146.         for(j = mfit-1; j > k; j--){
  147.             alpha[j][k] = alpha[k][j];
  148.         }
  149.     }
  150.  
  151.     if(veflag==2){
  152.         printf("\nalpha: ");
  153.         for(j=0;j<mfit;j++){
  154.             printf("\n %d   ",j);
  155.             for(k=0;k<mfit;k++) printf(" %g ",alpha[j][k]);
  156.         }
  157.         printf("\n");
  158.     }
  159.  
  160.     /* solve the matrix equation alpha*atry = beta, where alpha is */
  161.     /* a matrix, and atry and beta are vectors */
  162.     /* covar is the inverse of alpha */
  163.     if(debug) printf("calling solve_for_da()\n");
  164.     solve_for_da(alpha, covar, beta, atry, mfit);
  165.     if(debug) printf("returned from solve_for_da()\n");
  166.  
  167.     if(veflag==2){
  168.         printf("\ncovar: ");
  169.         for(j=0;j<mfit;j++){
  170.             printf("\n %d   ",j);
  171.             for(k=0;k<mfit;k++) printf("\n covar[%d][%d]= %g ",j, k, covar[j][k]);
  172.         }
  173.     printf("\n");
  174.     }
  175.  
  176. if(debug){
  177. for(j=0; j < mfit; j++)printf("\n atry%d = %g", j, atry[j]);
  178. printf("\n");
  179. }
  180.  
  181. for(j = 0; j < ma; j++) a[j] = 0;
  182. for(j = 0; j < mfit; j++) a[lista[j]] = atry[j];
  183.  
  184. /* call calc_yfit to compute chisqr */
  185. failed = calc_yfit(func, data, order, num_indep, a, ndata, ma, &chisqr);
  186.  
  187.  
  188. if(grflag){
  189.     if(MYPLOT && (num_indep == 1)){
  190. #ifdef OS2
  191.         if(debug)printf("in linear_fit(), calling myplot()\n");
  192.         failed=myplot(func,data, order, num_indep, ndata, filename, comment, a, ma);
  193.         if(debug)printf("in linear_fit(), myplot() returned %d\n", failed);
  194. #endif
  195.     }
  196.     else{
  197.         if(debug)printf("in linear_fit(), calling plot()\n");
  198.         failed=plot(func,data, order, num_indep, ndata, filename, comment, a, ma);
  199.         if(debug)printf("in linear_fit(), plot() returned %d\n", failed);
  200.     }
  201. }
  202.  
  203. for(j=0; j < ma; j++)printf("\n a%d = %g", j, a[j]);
  204. *chisq = chisqr;
  205. printf("\n chisqr = %g", *chisq);
  206. printf("\n");
  207. }
  208.  
  209. free_dmatrix(alpha,mfit,mfit);
  210. free(X);
  211. free(beta);
  212. free(dyda);
  213. free(fita);
  214. free(x);
  215. free(dydx);
  216. free(dydx_flag);
  217. free(aj);
  218. free(atry);
  219. return failed;
  220. }
  221.