Search:

Return to previous page

Contents of file 'boxcount/boxcount.c':



    1   /*7:*/
    2   #line 758 "./boxcount.w"
    3   
    4   /*8:*/
    5   #line 789 "./boxcount.w"
    6   
    7   #include <math.h> 
    8   #include <time.h> 
    9   #include <stdio.h> 
   10   #include <stdlib.h> 
   11   #include <string.h> 
   12   #include <ctype.h> 
   13   
   14   /*:8*/
   15   #line 759 "./boxcount.w"
   16   
   17   /*9:*/
   18   #line 810 "./boxcount.w"
   19   
   20   #define VERSION "1.6"
   21   #define COPYRIGHT "Copyright (C) 2006-2011, Fredrik Jonsson"
   22   #define SUCCESS (0)
   23   #define FAILURE (1)
   24   #define NCHMAX (256)
   25   #define APPEND 1
   26   #define OVERWRITE 2
   27   
   28   /*:9*/
   29   #line 760 "./boxcount.w"
   30   
   31   /*10:*/
   32   #line 825 "./boxcount.w"
   33   
   34   extern char*optarg;
   35   char*progname;
   36   
   37   /*:10*/
   38   #line 761 "./boxcount.w"
   39   
   40   /*23:*/
   41   #line 1306 "./boxcount.w"
   42   
   43   /*24:*/
   44   #line 1320 "./boxcount.w"
   45   
   46   void moment(double data[],int nnmin,int nnmax,double*ave,double*adev,
   47   double*sdev,double*var,double*skew,double*curt){
   48   int j;
   49   double ep= 0.0,s,p;
   50   if(nnmax-nnmin<=1){
   51   fprintf(stderr,"%s: Error in routine moment()! ",progname);
   52   fprintf(stderr,"(nnmax-nnmin>1 is a requirement)\n");
   53   exit(FAILURE);
   54   }
   55   s= 0.0;
   56   for(j= nnmin;j<=nnmax;j++)s+= data[j];
   57   *ave= s/(nnmax-nnmin+1);
   58   *adev= (*var)= (*skew)= (*curt)= 0.0;
   59   for(j= nnmin;j<=nnmax;j++){
   60   *adev+= fabs(s= data[j]-(*ave));
   61   ep+= s;
   62   ep+= s;
   63   *var+= (p= s*s);
   64   *skew+= (p*= s);
   65   *curt+= (p*= s);
   66   }
   67   *adev/= (nnmax-nnmin+1);
   68   *var= (*var-ep*ep/(nnmax-nnmin+1))/(nnmax-nnmin);
   69   *sdev= sqrt(*var);
   70   if(*var){
   71   *skew/= ((nnmax-nnmin+1)*(*var)*(*sdev));
   72   *curt= (*curt)/((nnmax-nnmin+1)*(*var)*(*var))-3.0;
   73   }else{
   74   fprintf(stderr,"%s: Error in routine moment()! ",progname);
   75   fprintf(stderr,"No skew/kurtosis for zero variance\n");
   76   exit(FAILURE);
   77   }
   78   }
   79   
   80   /*:24*/
   81   #line 1307 "./boxcount.w"
   82   
   83   /*25:*/
   84   #line 1361 "./boxcount.w"
   85   
   86   long int num_coordinate_pairs(FILE*trajectory_file){
   87   double tmp;
   88   int tmpch;
   89   long int mm= 0;
   90   fseek(trajectory_file,0L,SEEK_SET);
   91   while((tmpch= getc(trajectory_file))!=EOF){
   92   ungetc(tmpch,trajectory_file);
   93   fscanf(trajectory_file,"%lf",&tmp);
   94   fscanf(trajectory_file,"%lf",&tmp);
   95   mm++;
   96   tmpch= getc(trajectory_file);
   97   while((tmpch!=EOF)&&(!isdigit(tmpch)))tmpch= getc(trajectory_file);
   98   if(tmpch!=EOF)ungetc(tmpch,trajectory_file);
   99   }
  100   fseek(trajectory_file,0L,SEEK_SET);
  101   return(mm);
  102   }
  103   
  104   
  105   /*:25*/
  106   #line 1308 "./boxcount.w"
  107   
  108   /*26:*/
  109   #line 1390 "./boxcount.w"
  110   
  111   /*27:*/
  112   #line 1399 "./boxcount.w"
  113   
  114   short pathcharacter(int ch){
  115   return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
  116   (ch=='-')||(ch=='+'));
  117   }
  118   
  119   /*:27*/
  120   #line 1391 "./boxcount.w"
  121   
  122   /*28:*/
  123   #line 1411 "./boxcount.w"
  124   
  125   char*strip_away_path(char filename[]){
  126   int j,k= 0;
  127   while(pathcharacter(filename[k]))k++;
  128   j= (--k);
  129   while(isalnum((int)(filename[j])))j--;
  130   j++;
  131   return(&filename[j]);
  132   }
  133   
  134   /*:28*/
  135   #line 1392 "./boxcount.w"
  136   
  137   
  138   /*:26*/
  139   #line 1309 "./boxcount.w"
  140   
  141   /*29:*/
  142   #line 1426 "./boxcount.w"
  143   
  144   /*30:*/
  145   #line 1437 "./boxcount.w"
  146   
  147   double*dvector(long nl,long nh){
  148   double*v;
  149   v= (double*)malloc((size_t)((nh-nl+2)*sizeof(double)));
  150   if(!v){
  151   fprintf(stderr,"Error: Allocation failure in dvector()\n");
  152   exit(FAILURE);
  153   }
  154   return v-nl+1;
  155   }
  156   
  157   /*:30*/
  158   #line 1427 "./boxcount.w"
  159   
  160   /*33:*/
  161   #line 1473 "./boxcount.w"
  162   
  163   void free_ivector(int*v,long nl,long nh){
  164   free((char*)(v+nl-1));
  165   }
  166   
  167   /*:33*//*35:*/
  168   #line 1495 "./boxcount.w"
  169   
  170   void free_livector(long int*v,long nl,long nh){
  171   free((char*)(v+nl-1));
  172   }
  173   
  174   /*:35*/
  175   #line 1428 "./boxcount.w"
  176   
  177   /*32:*/
  178   #line 1459 "./boxcount.w"
  179   
  180   int*ivector(long nl,long nh){
  181   int*v;
  182   v= (int*)malloc((size_t)((nh-nl+2)*sizeof(int)));
  183   if(!v){
  184   fprintf(stderr,"Error: Allocation failure in ivector()\n");
  185   exit(FAILURE);
  186   }
  187   return v-nl+1;
  188   }
  189   
  190   /*:32*//*34:*/
  191   #line 1481 "./boxcount.w"
  192   
  193   long int*livector(long nl,long nh){
  194   long int*v;
  195   v= (long int*)malloc((size_t)((nh-nl+2)*sizeof(long int)));
  196   if(!v){
  197   fprintf(stderr,"Error: Allocation failure in livector()\n");
  198   exit(FAILURE);
  199   }
  200   return v-nl+1;
  201   }
  202   
  203   /*:34*/
  204   #line 1429 "./boxcount.w"
  205   
  206   /*31:*/
  207   #line 1451 "./boxcount.w"
  208   
  209   void free_dvector(double*v,long nl,long nh){
  210   free((char*)(v+nl-1));
  211   }
  212   
  213   /*:31*/
  214   #line 1430 "./boxcount.w"
  215   
  216   /*36:*/
  217   #line 1504 "./boxcount.w"
  218   
  219   short int**simatrix(long nrl,long nrh,long ncl,long nch)
  220   {
  221   long i,nrow= nrh-nrl+1,ncol= nch-ncl+1;
  222   short int**m;
  223   m= (short int**)malloc((size_t)((nrow+1)*sizeof(short int*)));
  224   if(!m){
  225   fprintf(stderr,"%s: Allocation failure 1 in simatrix()\n",progname);
  226   exit(FAILURE);
  227   }
  228   m+= 1;
  229   m-= nrl;
  230   m[nrl]= (short int*)malloc((size_t)((nrow*ncol+1)*sizeof(short int)));
  231   if(!m[nrl]){
  232   fprintf(stderr,"%s: Allocation failure 2 in simatrix()\n",progname);
  233   exit(FAILURE);
  234   }
  235   m[nrl]+= 1;
  236   m[nrl]-= ncl;
  237   for(i= nrl+1;i<=nrh;i++)m[i]= m[i-1]+ncol;
  238   return m;
  239   }
  240   
  241   /*:36*/
  242   #line 1431 "./boxcount.w"
  243   
  244   /*37:*/
  245   #line 1530 "./boxcount.w"
  246   
  247   void free_simatrix(short int**m,long nrl,long nrh,long ncl,long nch){
  248   free((char*)(m[nrl]+ncl-1));
  249   free((char*)(m+nrl-1));
  250   }
  251   
  252   /*:37*/
  253   #line 1432 "./boxcount.w"
  254   
  255   
  256   /*:29*/
  257   #line 1310 "./boxcount.w"
  258   
  259   /*38:*/
  260   #line 1538 "./boxcount.w"
  261   
  262   void showsomehelp(void){
  263   fprintf(stderr,"Usage: %s [options]\n",progname);
  264   fprintf(stderr,"Options:\n");
  265   fprintf(stderr," -h, --help\n");
  266   fprintf(stderr,"    Display this help message and exit clean.\n");
  267   fprintf(stderr," -v, --verbose\n");
  268   fprintf(stderr,"    Toggle verbose mode on/off.\n");
  269   fprintf(stderr," -o, --outputfile <str>\n");
  270   fprintf(stderr,"    Specifies the base name of the output files where the program\n");
  271   fprintf(stderr,"    is to save the calculated data. If the --outputfile or -o\n");
  272   fprintf(stderr,"    option is followed by 'append' the estimate for the fractal\n");
  273   fprintf(stderr,"    dimension will be appended to the file named <str>.dat, which\n");
  274   fprintf(stderr,"    will be created if it does not exist. If the following word\n");
  275   fprintf(stderr,"    instead is `overwrite' the file will instead be overwritten.\n");
  276   fprintf(stderr," -i, --trajectoryfile\n");
  277   fprintf(stderr,"    Specifies the input trajectory of the graph whose fractal\n");
  278   fprintf(stderr,"    dimension is to be estimated.\n");
  279   fprintf(stderr," -w, --computationwindow llx <num> lly <num> urx <num> ury <num>\n");
  280   fprintf(stderr,"    This option explicitly specifies the domain over which the\n");
  281   fprintf(stderr,"    box-counting algorithm will be applied, in terms of the\n");
  282   fprintf(stderr,"    lower-left and upper-right corners (llx,lly) and (urx,ury),\n");
  283   fprintf(stderr,"    respectively. By specifying this option, any automatic\n");
  284   fprintf(stderr,"    calculation of computational window will be neglected.\n");
  285   fprintf(stderr," -m, --boxmaps\n");
  286   fprintf(stderr,"    If this option is present, the program will generate\n");
  287   fprintf(stderr,"    MetaPost code for maps of the distribution of boxes.\n");
  288   fprintf(stderr,"    In doing so, also the input trajectory is included in\n");
  289   fprintf(stderr,"    the graphs. The convention for the naming of the output\n");
  290   fprintf(stderr,"    map files is that they are saved as <str>.<N>.dat,\n");
  291   fprintf(stderr,"    where <str> is the base filename as specified using the\n");
  292   fprintf(stderr,"    -o or --outputfile option, <N> is the automatically appended\n");
  293   fprintf(stderr,"    current level of resolution refinement in the box-counting,\n");
  294   fprintf(stderr,"    and where '.dat' is the file suffix as automatically appended\n");
  295   fprintf(stderr,"    by the program.\n");
  296   fprintf(stderr," --graphsize <width in mm> <height in mm>\n");
  297   fprintf(stderr,"    If the -m or --boxmaps option is present at the command line,\n");
  298   fprintf(stderr,"    then the --graphsize option will override the default graph\n");
  299   fprintf(stderr,"    size of the generated boxmaps. (Default graph size is 80 mm\n");
  300   fprintf(stderr,"    width and 56 mm height.)\n");
  301   fprintf(stderr," -Nmin, --minlevel <num>\n");
  302   fprintf(stderr,"    Specifies the minimum level Nmin of grid refinement that \n");
  303   fprintf(stderr,"    will be used in the evaluation of the estimate of the fractal\n");
  304   fprintf(stderr,"    dimension. With this option specified, the coarsest level\n");
  305   fprintf(stderr,"    used in the box-counting will be a [(2^Nmin)x(2^Nmin)]-grid\n");
  306   fprintf(stderr,"    of boxes.\n");
  307   fprintf(stderr," -Nmax, --maxlevel <num>\n");
  308   fprintf(stderr,"    This option is similar to the --minlevel or -Nmin options,\n");
  309   fprintf(stderr,"    with the difference that it instead specifies the maximum\n");
  310   fprintf(stderr,"    level Nmax of grid refinement that will be used in the\n");
  311   fprintf(stderr,"    evaluation of the estimate of the fractal dimension.\n");
  312   }
  313   
  314   /*:38*/
  315   #line 1311 "./boxcount.w"
  316   
  317   /*39:*/
  318   #line 1647 "./boxcount.w"
  319   
  320   short int lines_intersect(double p1x,double p1y,double q1x,double q1y,
  321   double p2x,double p2y,double q2x,double q2y){
  322   double a,b,c,d,e,f,det,tmp1,tmp2;
  323   short int intersect;
  324   a= q1x-p1x;
  325   b= p2x-q2x;
  326   c= q1y-p1y;
  327   d= p2y-q2y;
  328   e= p2x-p1x;
  329   f= p2y-p1y;
  330   det= a*d-b*c;
  331   tmp1= e*d-b*f;
  332   tmp2= a*f-e*c;
  333   intersect= 0;
  334   if(det> 0){
  335   if(((0.0<=tmp1)&&(tmp1<=det))&&((0.0<=tmp2)&&(tmp2<=det)))intersect= 1;
  336   }else if(det<0){
  337   if(((det<=tmp1)&&(tmp1<=0.0))&&((det<=tmp2)&&(tmp2<=0.0)))intersect= 1;
  338   }
  339   return(intersect);
  340   }
  341   
  342   /*:39*/
  343   #line 1312 "./boxcount.w"
  344   
  345   /*40:*/
  346   #line 1692 "./boxcount.w"
  347   
  348   short int box_intersection(double px,double py,double qx,double qy,
  349   double llx,double lly,double urx,double ury){
  350   if(lines_intersect(px,py,qx,qy,llx,lly,urx,lly))
  351   return(1);
  352   if(lines_intersect(px,py,qx,qy,urx,lly,urx,ury))
  353   return(1);
  354   if(lines_intersect(px,py,qx,qy,urx,ury,llx,ury))
  355   return(1);
  356   if(lines_intersect(px,py,qx,qy,llx,ury,llx,lly))
  357   return(1);
  358   return(0);
  359   }
  360   
  361   /*:40*/
  362   #line 1313 "./boxcount.w"
  363   
  364   /*41:*/
  365   #line 1714 "./boxcount.w"
  366   
  367   /*42:*/
  368   #line 1762 "./boxcount.w"
  369   
  370   long unsigned int get_num_covering_boxes_with_boxmaps(
  371   double*x_traj,double*y_traj,
  372   long int mm,int resolution,double global_llx,double global_lly,
  373   double global_urx,double global_ury,char boxgraph_filename[],
  374   double boxgraph_width,double boxgraph_height,
  375   char trajectory_filename[]){
  376   short int**box;
  377   long unsigned int i,j,m,nn,istart,jstart,istop,jstop,iincr,jincr,num_boxes;
  378   double*x_box,*y_box;
  379   
  380   double px,py,qx,qy;
  381   FILE*boxgraph_file;
  382   
  383   /*43:*/
  384   #line 1805 "./boxcount.w"
  385   
  386   {
  387   nn= 1;
  388   for(i= 1;i<=resolution;i++)nn= 2*nn;
  389   box= simatrix(1,nn,1,nn);
  390   for(i= 1;i<=nn;i++)for(j= 1;j<=nn;j++)box[i][j]= 0;
  391   x_box= dvector(1,nn+1);
  392   y_box= dvector(1,nn+1);
  393   for(m= 1;m<=nn+1;m++){
  394   x_box[m]= global_llx+((double)(m-1))*(global_urx-global_llx)/((double)(nn));
  395   y_box[m]= global_lly+((double)(m-1))*(global_ury-global_lly)/((double)(nn));
  396   }
  397   }
  398   
  399   /*:43*/
  400   #line 1776 "./boxcount.w"
  401   
  402   /*44:*/
  403   #line 1822 "./boxcount.w"
  404   
  405   {
  406   istart= 0;
  407   jstart= 0;
  408   for(m= 1;m<=nn;m++){
  409   if((x_box[m]<=x_traj[1])&&(x_traj[1]<=x_box[m+1]))istart= m;
  410   if((y_box[m]<=y_traj[1])&&(y_traj[1]<=y_box[m+1]))jstart= m;
  411   }
  412   if((istart==0)||(jstart==0)){
  413   fprintf(stderr,"%s: Error! Cannot find box indices of 1st coordinate!\n",
  414   progname);
  415   fprintf(stderr,"%s: Please check data for input trajetory.\n",progname);
  416   exit(FAILURE);
  417   }
  418   }
  419   
  420   
  421   /*:44*/
  422   #line 1777 "./boxcount.w"
  423   
  424   for(m= 1;m<=mm-1;m++){
  425   /*45:*/
  426   #line 1842 "./boxcount.w"
  427   {
  428   px= x_traj[m];
  429   py= y_traj[m];
  430   qx= x_traj[m+1];
  431   qy= y_traj[m+1];
  432   istop= istart;
  433   jstop= jstart;
  434   if(px<qx){
  435   iincr= 1;
  436   while(x_box[istop+1]<qx)istop++;
  437   }else{
  438   iincr= -1;
  439   while(x_box[istop]> qx)istop--;
  440   }
  441   if(py<qy){
  442   jincr= 1;
  443   while(y_box[jstop+1]<qy)jstop++;
  444   }else{
  445   jincr= -1;
  446   while(y_box[jstop]> qy)jstop--;
  447   }
  448   if(0==1){
  449   fprintf(stdout,"%s: Endpoint box indices: (i=%ld,j=%ld)\n",
  450   progname,istop,jstop);
  451   }
  452   }
  453   
  454   /*:45*/
  455   #line 1779 "./boxcount.w"
  456   
  457   /*46:*/
  458   #line 1878 "./boxcount.w"
  459   
  460   {
  461   for(i= istart;i!=(istop+iincr);i+= iincr){
  462   for(j= jstart;j!=(jstop+jincr);j+= jincr){
  463   if(box_intersection(px,py,qx,qy,
  464   x_box[i],y_box[j],x_box[i+1],y_box[j+1])){
  465   box[i][j]= 1;
  466   }
  467   }
  468   }
  469   istart= istop;
  470   jstart= jstop;
  471   }
  472   
  473   /*:46*/
  474   #line 1780 "./boxcount.w"
  475   
  476   }
  477   /*47:*/
  478   #line 1898 "./boxcount.w"
  479   {
  480   if(!strcmp(boxgraph_filename,"")){
  481   boxgraph_file= NULL;
  482   }else{
  483   if((boxgraph_file= fopen(boxgraph_filename,"w"))==NULL){
  484   fprintf(stderr,"%s: Could not open %s for box graphs!\n",
  485   progname,boxgraph_filename);
  486   exit(FAILURE);
  487   }
  488   fseek(boxgraph_file,0L,SEEK_SET);
  489   fprintf(boxgraph_file,"input graph;\n");
  490   fprintf(boxgraph_file,"def box(expr i,j)=\n");
  491   fprintf(boxgraph_file,"begingroup\n");
  492   fprintf(boxgraph_file,"llx:=%2.8f+(i-1)*%2.8f;\n",
  493   global_llx,(global_urx-global_llx)/(nn));
  494   fprintf(boxgraph_file,"lly:=%2.8f+(j-1)*%2.8f;\n",
  495   global_lly,(global_ury-global_lly)/(nn));
  496   fprintf(boxgraph_file,"urx:=%2.8f+(i)*%2.8f;\n",
  497   global_llx,(global_urx-global_llx)/(nn));
  498   fprintf(boxgraph_file,"ury:=%2.8f+(j)*%2.8f;\n",
  499   global_lly,(global_ury-global_lly)/(nn));
  500   fprintf(boxgraph_file,"gdraw (llx,lly)--(urx,lly);\n");
  501   fprintf(boxgraph_file,"gdraw (urx,lly)--(urx,ury);\n");
  502   fprintf(boxgraph_file,"gdraw (urx,ury)--(llx,ury);\n");
  503   fprintf(boxgraph_file,"gdraw (llx,ury)--(llx,lly);\n");
  504   fprintf(boxgraph_file,"endgroup\n");
  505   fprintf(boxgraph_file,"enddef;\n");
  506   fprintf(boxgraph_file,"beginfig(1);\n");
  507   fprintf(boxgraph_file,"w:=%2.4fmm; h:=%2.4fmm;\n",
  508   boxgraph_width,boxgraph_height);
  509   fprintf(boxgraph_file,"draw begingraph(w,h);\n");
  510   fprintf(boxgraph_file,"pickup pencircle scaled .3pt;\n");
  511   fprintf(boxgraph_file,"setrange(%2.6f,%2.6f,%2.6f,%2.6f);\n",
  512   global_llx,global_lly,global_urx,global_ury);
  513   }
  514   }
  515   
  516   /*:47*/
  517   #line 1782 "./boxcount.w"
  518   
  519   /*48:*/
  520   #line 1964 "./boxcount.w"
  521   {
  522   if(boxgraph_file!=NULL){
  523   fprintf(boxgraph_file,"pickup pencircle scaled .5pt;\n");
  524   i= 0;
  525   for(m= 1;m<=mm;m++){
  526   if(i==0){
  527   if(m==1){
  528   fprintf(boxgraph_file,"gdraw (%2.4f,%2.4f)",
  529   x_traj[m],y_traj[m]);
  530   }else if(2<mm){
  531   fprintf(boxgraph_file,"gdraw (%2.4f,%2.4f)--(%2.4f,%2.4f)",
  532   x_traj[m-1],y_traj[m-1],x_traj[m],y_traj[m]);
  533   }
  534   }else{
  535   fprintf(boxgraph_file,"--(%2.4f,%2.4f)",x_traj[m],y_traj[m]);
  536   }
  537   i++;
  538   if((i==5)||(m==mm)){
  539   fprintf(boxgraph_file,";\n");
  540   i= 0;
  541   }
  542   }
  543   }
  544   }
  545   
  546   /*:48*/
  547   #line 1783 "./boxcount.w"
  548   
  549   /*49:*/
  550   #line 1999 "./boxcount.w"
  551   
  552   {
  553   num_boxes= 0;
  554   for(i= 1;i<=nn;i++){
  555   for(j= 1;j<=nn;j++){
  556   if(box[i][j]==1){
  557   num_boxes++;
  558   if(boxgraph_file!=NULL){
  559   fprintf(boxgraph_file,"box(%ld,%ld);\n",i,j);
  560   }
  561   }
  562   }
  563   }
  564   free_simatrix(box,1,nn,1,nn);
  565   }
  566   
  567   /*:49*/
  568   #line 1784 "./boxcount.w"
  569   
  570   /*50:*/
  571   #line 2020 "./boxcount.w"
  572   {
  573   if(boxgraph_file!=NULL){
  574   fprintf(boxgraph_file,"autogrid(itick bot,itick lft);\n");
  575   fprintf(boxgraph_file,"glabel.bot(btex $x$ etex,OUT);\n");
  576   fprintf(boxgraph_file,"glabel.lft(btex $y$ etex,OUT);\n");
  577   fprintf(boxgraph_file,"endgraph;\n");
  578   fprintf(boxgraph_file,"endfig;\n");
  579   fprintf(boxgraph_file,"end\n");
  580   }
  581   }
  582   
  583   /*:50*/
  584   #line 1785 "./boxcount.w"
  585   
  586   /*51:*/
  587   #line 2033 "./boxcount.w"
  588   {
  589   if(boxgraph_file!=NULL){
  590   fprintf(stdout,
  591   "%s: MetaPost output box distribution graph written to %s.\n",
  592   progname,boxgraph_filename);
  593   fclose(boxgraph_file);
  594   }
  595   }
  596   
  597   /*:51*/
  598   #line 1786 "./boxcount.w"
  599   
  600   return(num_boxes);
  601   }
  602   
  603   /*:42*/
  604   #line 1715 "./boxcount.w"
  605   
  606   /*52:*/
  607   #line 2078 "./boxcount.w"
  608   
  609   long unsigned int get_num_covering_boxes(double*x_traj,double*y_traj,long int mm,
  610   int i,double global_llx,double global_lly,
  611   double global_urx,double global_ury){
  612   return(get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i,
  613   global_llx,global_lly,global_urx,global_ury,"",0.0,0.0,""));
  614   }
  615   
  616   /*:52*/
  617   #line 1716 "./boxcount.w"
  618   
  619   
  620   /*:41*/
  621   #line 1314 "./boxcount.w"
  622   
  623   
  624   
  625   /*:23*/
  626   #line 762 "./boxcount.w"
  627   
  628   
  629   int main(int argc,char*argv[])
  630   {
  631   /*11:*/
  632   #line 904 "./boxcount.w"
  633   
  634   short verbose,user_set_compwin,output_mode,make_boxmaps;
  635   long int*num_boxes;
  636   time_t initime;
  637   time_t now;
  638   long mm,nn,nnmin,nnmax;
  639   double global_llx,global_lly,global_urx,global_ury;
  640   double*x_traj,*y_traj,*x,*y;
  641   FILE*trajectory_file,*frac_estimate_file,*boxmap_file;
  642   char trajectory_filename[NCHMAX],output_filename[NCHMAX],frac_estimate_filename[NCHMAX];
  643   char boxmap_filename[NCHMAX];
  644   double boxgraph_width,boxgraph_height;
  645   int no_arg;
  646   long i;
  647   double*fracdimen_estimates,ave,adev,sdev,var,skew,curt;
  648   
  649   /*:11*/
  650   #line 766 "./boxcount.w"
  651   
  652   /*12:*/
  653   #line 922 "./boxcount.w"
  654   
  655   verbose= 0;
  656   user_set_compwin= 0;
  657   output_mode= OVERWRITE;
  658   make_boxmaps= 0;
  659   nnmin= 0;
  660   nnmax= 10;
  661   strcpy(output_filename,"out");
  662   strcpy(trajectory_filename,"");
  663   boxgraph_width= 80.0;
  664   boxgraph_height= 56.0;
  665   trajectory_file= NULL;
  666   frac_estimate_file= NULL;
  667   boxmap_file= NULL;
  668   initime= time(NULL);
  669   
  670   /*:12*/
  671   #line 767 "./boxcount.w"
  672   
  673   /*13:*/
  674   #line 944 "./boxcount.w"
  675   
  676   {
  677   progname= strip_away_path(argv[0]);
  678   fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
  679   no_arg= argc;
  680   while(--argc){
  681   if(!strcmp(argv[no_arg-argc],"-o")||
  682   !strcmp(argv[no_arg-argc],"--outputfile")){
  683   --argc;
  684   if(!strcmp(argv[no_arg-argc],"append")||
  685   !strcmp(argv[no_arg-argc],"a")){
  686   output_mode= APPEND;
  687   }else if(!strcmp(argv[no_arg-argc],"overwrite")||
  688   !strcmp(argv[no_arg-argc],"o")){
  689   output_mode= OVERWRITE;
  690   }else{
  691   fprintf(stderr,"%s: Error in '-o' or '--outputfile' option!",
  692   progname);
  693   exit(FAILURE);
  694   }
  695   --argc;
  696   strcpy(output_filename,argv[no_arg-argc]);
  697   }else if(!strcmp(argv[no_arg-argc],"-w")||
  698   !strcmp(argv[no_arg-argc],"--computationwindow")){
  699   user_set_compwin= 1;
  700   --argc;
  701   if(!strcmp(argv[no_arg-argc],"llx")){
  702   --argc;
  703   if(!sscanf(argv[no_arg-argc],"%lf",&global_llx)){
  704   fprintf(stderr,"%s: Error in parsing lower-left x-value.\n",
  705   progname);
  706   exit(FAILURE);
  707   }
  708   }else{
  709   fprintf(stderr,"%s: Error in computation window option\n",
  710   progname);
  711   fprintf(stderr,"%s: Expecting 'llx'\n",progname);
  712   exit(FAILURE);
  713   }
  714   --argc;
  715   if(!strcmp(argv[no_arg-argc],"lly")){
  716   --argc;
  717   if(!sscanf(argv[no_arg-argc],"%lf",&global_lly)){
  718   fprintf(stderr,"%s: Error in parsing lower-left y-value.\n",
  719   progname);
  720   exit(FAILURE);
  721   }
  722   }else{
  723   fprintf(stderr,"%s: Error in computation window option\n",
  724   progname);
  725   fprintf(stderr,"%s: Expecting 'lly'\n",progname);
  726   exit(FAILURE);
  727   }
  728   --argc;
  729   if(!strcmp(argv[no_arg-argc],"urx")){
  730   --argc;
  731   if(!sscanf(argv[no_arg-argc],"%lf",&global_urx)){
  732   fprintf(stderr,"%s: Error in parsing lower-left x-value.\n",
  733   progname);
  734   exit(FAILURE);
  735   }
  736   }else{
  737   fprintf(stderr,"%s: Error in computation window option\n",
  738   progname);
  739   fprintf(stderr,"%s: Expecting 'urx'\n",progname);
  740   exit(FAILURE);
  741   }
  742   --argc;
  743   if(!strcmp(argv[no_arg-argc],"ury")){
  744   --argc;
  745   if(!sscanf(argv[no_arg-argc],"%lf",&global_ury)){
  746   fprintf(stderr,"%s: Error in parsing lower-left y-value.\n",
  747   progname);
  748   exit(FAILURE);
  749   }
  750   }else{
  751   fprintf(stderr,"%s: Error in computation window option\n",
  752   progname);
  753   fprintf(stderr,"%s: Expecting 'ury'\n",progname);
  754   exit(FAILURE);
  755   }
  756   }else if(!strcmp(argv[no_arg-argc],"-i")||
  757   !strcmp(argv[no_arg-argc],"--trajectoryfile")){
  758   --argc;
  759   strcpy(trajectory_filename,argv[no_arg-argc]);
  760   }else if(!strcmp(argv[no_arg-argc],"-v")||
  761   !strcmp(argv[no_arg-argc],"--verbose")){
  762   verbose= (verbose?0:1);
  763   }else if(!strcmp(argv[no_arg-argc],"-h")||
  764   !strcmp(argv[no_arg-argc],"--help")){
  765   showsomehelp();
  766   exit(SUCCESS);
  767   }else if(!strcmp(argv[no_arg-argc],"-m")||
  768   !strcmp(argv[no_arg-argc],"--boxmaps")){
  769   make_boxmaps= 1;
  770   }else if(!strcmp(argv[no_arg-argc],"--graphsize")){
  771   --argc;
  772   if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_width)){
  773   fprintf(stderr,"%s: Error in width of '--graphsize' option.\n",
  774   progname);
  775   exit(FAILURE);
  776   }
  777   --argc;
  778   if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_height)){
  779   fprintf(stderr,"%s: Error in height of '--graphsize' option.\n",
  780   progname);
  781   exit(FAILURE);
  782   }
  783   }else if(!strcmp(argv[no_arg-argc],"--minlevel")||
  784   !strcmp(argv[no_arg-argc],"-Nmin")){
  785   --argc;
  786   if(!sscanf(argv[no_arg-argc],"%ld",&nnmin)){
  787   fprintf(stderr,"%s: Error in '--minlevel' or '-Nmin' option.\n",
  788   progname);
  789   exit(FAILURE);
  790   }
  791   }else if(!strcmp(argv[no_arg-argc],"--maxlevel")||
  792   !strcmp(argv[no_arg-argc],"-Nmax")){
  793   --argc;
  794   if(!sscanf(argv[no_arg-argc],"%ld",&nnmax)){
  795   fprintf(stderr,"%s: Error in '--maxlevel' or '-Nmax' option.\n",
  796   progname);
  797   exit(FAILURE);
  798   }
  799   }else{
  800   fprintf(stderr,"%s: Specified option '%s' invalid!\n",
  801   progname,argv[no_arg-argc]);
  802   showsomehelp();
  803   exit(FAILURE);
  804   }
  805   }
  806   }
  807   
  808   /*:13*/
  809   #line 768 "./boxcount.w"
  810   
  811   /*14:*/
  812   #line 1079 "./boxcount.w"
  813   
  814   {
  815   fprintf(stdout,"%s: Program execution started %s",progname,ctime(&initime));
  816   }
  817   
  818   /*:14*/
  819   #line 769 "./boxcount.w"
  820   
  821   /*15:*/
  822   #line 1089 "./boxcount.w"
  823   {
  824   if(!strcmp(trajectory_filename,"")){
  825   fprintf(stderr,"%s: No input trajectory specified!\n",progname);
  826   fprintf(stderr,"%s: Please use the '--trajectoryfile' option.\n",
  827   progname);
  828   fprintf(stderr,"%s: Use '--help' option to display a help message.\n",
  829   progname);
  830   exit(FAILURE);
  831   }
  832   if((trajectory_file= fopen(trajectory_filename,"r"))==NULL){
  833   fprintf(stderr,"%s: Could not open %s for loading trajectory!\n",
  834   progname,trajectory_filename);
  835   exit(FAILURE);
  836   }
  837   mm= num_coordinate_pairs(trajectory_file);
  838   x_traj= dvector(1,mm);
  839   y_traj= dvector(1,mm);
  840   for(i= 1;i<=mm;i++){
  841   fscanf(trajectory_file,"%lf",&x_traj[i]);
  842   fscanf(trajectory_file,"%lf",&y_traj[i]);
  843   }
  844   fclose(trajectory_file);
  845   if(verbose){
  846   fprintf(stdout,"%s: Loaded %ld coordinate pairs from file '%s'.\n",
  847   progname,mm,trajectory_filename);
  848   }
  849   }
  850   
  851   /*:15*/
  852   #line 770 "./boxcount.w"
  853   
  854   /*16:*/
  855   #line 1119 "./boxcount.w"
  856   {
  857   if(!strcmp(output_filename,"")){
  858   fprintf(stderr,"%s: No output base name specified!\n",progname);
  859   fprintf(stderr,"%s: Please use the '--outputfile' option.\n",
  860   progname);
  861   exit(FAILURE);
  862   }
  863   sprintf(frac_estimate_filename,"%s.dat",output_filename);
  864   if(output_mode==APPEND){
  865   if((frac_estimate_file= fopen(frac_estimate_filename,"a"))==NULL){
  866   fprintf(stderr,"%s: Could not open %s for output!\n",
  867   progname,frac_estimate_filename);
  868   exit(FAILURE);
  869   }
  870   fseek(frac_estimate_file,0L,SEEK_END);
  871   
  872   }else if(output_mode==OVERWRITE){
  873   if((frac_estimate_file= fopen(frac_estimate_filename,"w"))==NULL){
  874   fprintf(stderr,"%s: Could not open %s for loading trajectory!\n",
  875   progname,frac_estimate_filename);
  876   exit(FAILURE);
  877   }
  878   fseek(frac_estimate_file,0L,SEEK_SET);
  879   
  880   }else{
  881   fprintf(stderr,"%s: Error: Output mode (output_mode) undefined!",progname);
  882   exit(FAILURE);
  883   }
  884   }
  885   
  886   /*:16*/
  887   #line 771 "./boxcount.w"
  888   
  889   /*17:*/
  890   #line 1158 "./boxcount.w"
  891   
  892   {
  893   if(!user_set_compwin){
  894   global_llx= x_traj[1];
  895   global_lly= y_traj[1];
  896   global_urx= global_llx;
  897   global_ury= global_lly;
  898   for(i= 1;i<=mm;i++){
  899   if(x_traj[i]> global_urx)global_urx= x_traj[i];
  900   if(y_traj[i]> global_ury)global_ury= y_traj[i];
  901   if(x_traj[i]<global_llx)global_llx= x_traj[i];
  902   if(y_traj[i]<global_lly)global_lly= y_traj[i];
  903   }
  904   if(verbose){
  905   fprintf(stdout,"%s: Global box-counting window:\n",progname);
  906   fprintf(stdout,"%s:   (llx,lly)=(%2.8f,%2.8f)\n",progname,
  907   global_llx,global_lly);
  908   fprintf(stdout,"%s:   (urx,ury)=(%2.8f,%2.8f)\n",progname,
  909   global_urx,global_ury);
  910   }
  911   }
  912   }
  913   
  914   /*:17*/
  915   #line 772 "./boxcount.w"
  916   
  917   /*18:*/
  918   #line 1195 "./boxcount.w"
  919   {
  920   num_boxes= livector(1,nnmax);
  921   nn= 1;
  922   for(i= 1;i<=nnmin-1;i++)nn= 2*nn;
  923   for(i= nnmin;i<=nnmax;i++){
  924   nn= 2*nn;
  925   if(make_boxmaps){
  926   
  927   sprintf(boxmap_filename,"%s-%02ld.mp",output_filename,i);
  928   num_boxes[i]= get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i,
  929   global_llx,global_lly,global_urx,global_ury,boxmap_filename,
  930   boxgraph_width,boxgraph_height,trajectory_filename);
  931   }else{
  932   
  933   num_boxes[i]= get_num_covering_boxes(x_traj,y_traj,mm,i,
  934   global_llx,global_lly,global_urx,global_ury);
  935   }
  936   if(verbose){
  937   fprintf(stdout,"%s: N=%ld (%ldx%ld-grid of boxes): ",progname,i,nn,nn);
  938   fprintf(stdout,"Trajectory covered by %ld boxes\n",num_boxes[i]);
  939   }
  940   }
  941   }
  942   
  943   /*:18*/
  944   #line 773 "./boxcount.w"
  945   
  946   /*19:*/
  947   #line 1245 "./boxcount.w"
  948   {
  949   x= dvector(nnmin,nnmax);
  950   y= dvector(nnmin,nnmax);
  951   fracdimen_estimates= dvector(nnmin,nnmax);
  952   nn= 1;
  953   for(i= 1;i<=nnmax;i++){
  954   nn= 2*nn;
  955   if(nnmin<=i)x[i]= log((double)nn);
  956   }
  957   for(i= nnmin;i<=nnmax;i++)y[i]= log(num_boxes[i]);
  958   for(i= nnmin;i<=nnmax;i++)fracdimen_estimates[i]= y[i]/x[i];
  959   if(verbose){
  960   for(i= nnmin;i<=nnmax;i++){
  961   fprintf(stdout,"%s: N=%ld, fractal dimension estimate=%f\n",progname,
  962   i,fracdimen_estimates[i]);
  963   }
  964   }
  965   moment(fracdimen_estimates,nnmin,nnmax,&ave,&adev,&sdev,&var,&skew,&curt);
  966   if(verbose){
  967   fprintf(stdout,"%s: Estimate of fractal dimension: %f\n",progname,ave);
  968   fprintf(stdout,"%s: Average deviation: %f\n",progname,adev);
  969   fprintf(stdout,"%s: Standard deviation: %f\n",progname,sdev);
  970   fprintf(stdout,"%s: Skewness: %f\n",progname,skew);
  971   }
  972   free_livector(num_boxes,1,nnmax);
  973   free_dvector(fracdimen_estimates,nnmin,nnmax);
  974   free_dvector(x,nnmin,nnmax);
  975   free_dvector(y,nnmin,nnmax);
  976   }
  977   
  978   /*:19*/
  979   #line 774 "./boxcount.w"
  980   
  981   /*20:*/
  982   #line 1277 "./boxcount.w"
  983   {
  984   if(output_mode==APPEND){
  985   fseek(frac_estimate_file,0L,SEEK_END);
  986   }else if(output_mode==OVERWRITE){
  987   fseek(frac_estimate_file,0L,SEEK_SET);
  988   }
  989   fprintf(frac_estimate_file,"%f %f %f\n",ave,sdev,skew);
  990   }
  991   
  992   
  993   /*:20*/
  994   #line 775 "./boxcount.w"
  995   
  996   /*21:*/
  997   #line 1289 "./boxcount.w"
  998   {
  999   fclose(frac_estimate_file);
 1000   }
 1001   
 1002   /*:21*/
 1003   #line 776 "./boxcount.w"
 1004   
 1005   /*22:*/
 1006   #line 1295 "./boxcount.w"
 1007   {
 1008   now= time(NULL);
 1009   if(verbose)
 1010   fprintf(stdout,"%s: Total execution time: %d s\n",progname,
 1011   ((int)difftime(now,initime)));
 1012   fprintf(stdout,"%s: Program execution closed %s",progname,ctime(&now));
 1013   }
 1014   
 1015   /*:22*/
 1016   #line 777 "./boxcount.w"
 1017   
 1018   return(SUCCESS);
 1019   }
 1020   
 1021   /*:7*/
 1022   

Return to previous page

Generated by ::viewsrc::

Last modified Wednesday 15 Feb 2023