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
Generated by ::viewsrc::