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