home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / dimslib / dims_analyze.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-26  |  2.8 KB  |  140 lines
  1. /* Patrick Worfolk
  2.    January 24, 1990
  3.    Computes slope of the best least squares fit
  4.    to the logs of a set of data.
  5.    As per Numerical Recipes.
  6.    */
  7.  
  8. #include <stdio.h>
  9. #include <math.h>
  10.  
  11. static double sqrarg;
  12. #define SQR(a) (sqrarg=(a),sqrarg*sqrarg)
  13.  
  14. /* Compute least squares fit to the line y=a+bx and
  15.    compute standard deviations and chi square stats
  16.    */
  17. void dims_fit(x,y,n,a,b,siga,sigb,chi2)
  18. double *x,*y,*a,*b,*siga,*sigb,*chi2;
  19. int n;
  20. {
  21.     int i;
  22.     double t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
  23.  
  24.     *b=0.0;
  25.     for (i=0;i<n;i++) {
  26.         sx += x[i];
  27.         sy += y[i];
  28.     }
  29.     ss=n;
  30.     sxoss=sx/ss;
  31.     for (i=0;i<n;i++) {
  32.         t=x[i]-sxoss;
  33.         st2 += t*t;
  34.         *b += t*y[i];
  35.     }
  36.     *b /= st2;
  37.     *a=(sy-sx*(*b))/ss;
  38.     *siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
  39.     *sigb=sqrt(1.0/st2);
  40.     *chi2=0.0;
  41.     for (i=0;i<n;i++)
  42.         *chi2 += SQR(y[i]-(*a)-(*b)*x[i]);
  43.     sigdat=sqrt((*chi2)/(n-2));
  44.     *siga *= sigdat;
  45.     *sigb *= sigdat;
  46. }
  47. /* Analyze the dimension data by computing the best
  48.    linear fit to the chosen section of curve.
  49.    */
  50. void dims_analyze()
  51. {
  52.     double *logx,*logy,a,b,siga,sigb,chi2,log2(),*dvector(),*x,*y;
  53.     int i,m,n,dim,hi,lo;
  54.     FILE *fp,*fopen();
  55.     char s[80];
  56.     extern int dims_type_option,dims_lsf_start,dims_lsf_end;
  57.  
  58.  
  59.     /* Read input file to array */
  60.     sprintf(s,"kaos.dims.tmpdat");
  61.     fp = fopen(s,"r");
  62.  
  63.     dim = 2;
  64.     x = dvector(0,dim);
  65.     y = dvector(0,dim);
  66.     n = 0;
  67.     while ( fscanf(fp,"%lf %lf",x+n,y+n) != EOF){
  68.         n++;
  69.         if(n>=dim){
  70.             dim *= 2;
  71.             x = (double *) realloc(x,(unsigned)sizeof(double) * (dim+1));
  72.             y = (double *) realloc(y,(unsigned)sizeof(double) * (dim+1));
  73.         }
  74.     }
  75.     fclose(fp);
  76.  
  77.     lo=dims_lsf_start;
  78.     hi=dims_lsf_end;
  79.  
  80.     if (hi>n) {
  81.         system_mess_proc(1,"Cannot compute using all points.");
  82.         hi=n;
  83.     }
  84.     if (lo<1) {
  85.         system_mess_proc(1,"Cannot start before first point.");
  86.         lo=1;
  87.     }
  88.     if (lo>=hi) {
  89.         system_mess_proc(1,"There are no valid data points");
  90.         return;
  91.     }
  92.  
  93.     m = hi-lo;
  94.     /* Allocate memory */
  95.     logx = dvector(0,m);
  96.     if (!logx) {
  97.         system_mess_proc(1,"(Dim analyze) Insufficient memory");
  98.         return;
  99.     }
  100.     logy = dvector(0,m);
  101.     if (!logy) {
  102.         system_mess_proc(1,"(Dim analyze) Insufficient memory");
  103.         return;
  104.     }
  105.  
  106.     /* copy to log log data */
  107.     for (i=lo; i<=hi; i++) {
  108.         logx[i-lo] = log2(x[i-1]);
  109.         logy[i-lo] = log2(y[i-1]);
  110.     }
  111.  
  112.     /* send output */
  113.     if (dims_type_option==0) {
  114.         system_mess_proc(0,"Fractal dimension analysis:");
  115.     }
  116.     else if (dims_type_option==2) {
  117.         system_mess_proc(0,"Information dimension analysis:");
  118.     }
  119.  
  120.     /* send data to be analyzed */
  121.     if (lo+1==hi) {
  122.         sprintf(s,"Slope is: %le",(logy[1]-logy[0])/(logx[1]-logx[0]));
  123.     }
  124.     else {
  125.         dims_fit(logx,logy,m+1,&a,&b,&siga,&sigb,&chi2);
  126.         sprintf(s,"Best fit slope = %le",b);
  127.         system_mess_proc(0,s);
  128.         sprintf(s,"sd = %le, chi square = %le",sigb,chi2);
  129.     }
  130.     system_mess_proc(0,s);
  131.  
  132.  
  133.     /* free memory */
  134.     (void) free_dvector(logx,0,m);
  135.     (void) free_dvector(logy,0,m);
  136.     (void) free_dvector(x,0,dim);
  137.     (void) free_dvector(y,0,dim);
  138.  
  139. }
  140.