| Diploma Thesis Percolation Simulation C++ Sourcecode Documentation |
00001 // starters.h
00002 // (if necessary) asks for parameters
00003 // and starts the counters, etc.
00004 //
00005 // (c) 2000 Andreas Krueger
00006 // Version 1.0
00007 // last change: 31.10.2000
00008
00009 namespace starters {
00010
00011 using statistics::find_table_minmax;
00012 using statistics::neighbourhood_average;
00013 using statistics::muvarskewkurt;
00014 using statistics::calculate_moments;
00015 using statistics::manipdata;
00016 using statistics::averaging;
00017 using results::one_result;
00018 using results::all_results;
00019 using results::set_N_and_dim;
00020 using counters::naive_count_and_analyze;
00021 using counters::sort_naivecount_and_analyze;
00022 using counters::naive_listcount_and_analyze;
00023 using counters::copy_l_count_and_analyze;
00024 using counters::divide_all_spheres;
00025 using counters::choose_optimal_cuts;
00026 using counters::divide_by_c_cuts_count_and_analyze;
00027 using counters::throw_dnc_count_analyze_step; // used NOW
00028 using analyze::edgespheres;
00029 using analyze::show_twohisto;
00030 using analyze::sameresult;
00031 using analyze::occuring_clusternumbers;
00032 using datafiles::set_and_show_filename;
00033 using datafiles::write_intro;
00034 using datafiles::write_one_step;
00035 using datafiles::closefile;
00036 using datafiles::setfilename;
00037 using datafiles::write_intro7;
00038 using datafiles::write_results7;
00039 using datafiles::save_results8;
00040 using frontend::introduction;
00041 using frontend::ask_parameter;
00042 using frontend::ask_for_ff_range;
00043 using frontend::ask_for_parameters;
00044 using ff::ffnoc_guessed;
00045
00046
00047 void into_screen(int dimension) {
00048
00049 srand( (unsigned)time( NULL ) ); // set the random generator
00050
00051 NUMBER N; // actual number of spheres
00052 REAL R; // radius of spheres
00053 introduction (N,R, dimension); // ask for actual number of spheres and radius
00054 cout <<"N="<<N<<" R="<<R<<endl;
00055
00056 sphere *scatter = new sphere[N+1];
00057 set_dim(scatter,0,N,dimension);
00058 sphere *copy = new sphere[N+1];
00059 set_dim(copy,0,N,dimension);
00060 NUMBER *histogram=new NUMBER[(N+1)]; // initialize histogram
00061
00062 NUMBER biggestcl=0;
00063 REAL averagecl=0;
00064 clock_t time=0;
00065
00066 REAL sum_biggestcl=0;
00067 REAL sum_averagecl=0;
00068 REAL sum_time=0;
00069
00070 int loops;
00071 cout <<"number of loops: ";
00072 cin>> loops;
00073
00074 set_radius(scatter,1,N,R);
00075 int i;
00076 for (i=1;i<=loops;i++){
00077 throw_spheres(scatter, 1, N); // make randomly thrown spheres from 1 to N in array scatter
00078 set_clusternumber(scatter,1,N,0);
00079
00080 time=sort_naivecount_and_analyze(scatter, copy, 1, N, histogram, biggestcl, averagecl);
00081
00082 cout <<"loop "<<i<<" time: "<< ms(time) <<" ms";
00083 cout << " average size :"<<averagecl;
00084 cout <<" biggest: "<<biggestcl<<endl;
00085
00086 sum_time+=time;
00087 sum_biggestcl+=biggestcl;
00088 sum_averagecl+=averagecl;
00089 }
00090
00091 sum_time/=loops;
00092 sum_biggestcl/=loops;
00093 sum_averagecl/=loops;
00094
00095 waitanykey();
00096
00097 cout <<loops<<" loops took an averaged time of each time: "<< ms((clock_t)sum_time) <<" ms";
00098 cout << "\n averaged average size :"<<sum_averagecl;
00099 cout <<"\n averaged biggest: "<<sum_biggestcl<<endl;
00100
00101 waitanykey();
00102
00103 cout << "\nby the way...you wanted: "<<N<<" spheres with radius "<<R;
00104 cout << "\nNow try the same number with different radii... ";
00105
00106 cout << endl<<endl;
00107
00108
00109
00110 delete[] scatter;
00111
00112 delete[] copy;
00113 delete[] histogram;
00114 }
00115
00116
00117 void test_list (int dimension=2) {
00118
00119 NUMBER N;
00120 cout <<"Please type in NUMBER of spheres (a low number, eg. 15 !): ";
00121 cin >> N;
00122
00123 sphere *scatter = new sphere[N+1];
00124 set_dim(scatter,0,N,dimension);
00125 sphere *copy1 = new sphere[N+1];
00126 set_dim(copy1,0,N,dimension);
00127 sphere *copy2 = new sphere[N+1];
00128 set_dim(copy2,0,N,dimension);
00129
00130 throw_spheres(scatter, 1, N); // randomly throw spheres from 1 to N in array scatter
00131 set_clustersize(scatter,1,N,0);
00132
00133 cout <<"\nPlease type in the radius of one sphere (eg. 1000): ";
00134 REAL R; cin >> R;
00135
00136 set_radius(scatter,1,N,R);
00137
00138 NUMLIST all;
00139 increasing_integers (all, N); // fill the list "all" with the numbers of all spheres
00140
00141 NUMLIST sph1;
00142 NUMLIST sph2;
00143
00144
00145 NUMBER biggestcl1, biggestcl2;
00146 REAL averagecl1, averagecl2 ;
00147 NUMBER *histogram1=new NUMBER[(N+1)];
00148 NUMBER *histogram2=new NUMBER[(N+1)];
00149
00150
00151 divide_all_spheres(scatter, all,sph1,sph2,(COORDFLOAT)(GRIDSIZE/2),1);
00152
00153 NUMLIST::iterator iter1;
00154 NUMLIST::iterator iter2;
00155
00156 set_clusternumber_l(scatter,sph1,0);
00157 cout << "\n-> "<<sph1.size()<<" with first component lower than "<<GRIDSIZE/2<<endl;
00158 for (iter1 = sph1.begin(); iter1 != sph1.end(); iter1++){
00159 cout << *iter1 << "\t";
00160 cout << scatter[*iter1] <<endl;
00161 }
00162
00163 REAL radius=find_biggestradius(scatter, all);
00164
00165 NUMLIST e1;
00166 edgespheres(scatter, sph1, e1, GRIDSIZE/2,2*radius, +1);
00167
00168 cout << "Egdespheres (close enough to the "<<GRIDSIZE/2<<" border) are: ";
00169 for (iter2=e1.begin();iter2 !=e1.end();iter2++){
00170 cout << *iter2 <<" ";
00171 }
00172 cout << endl;
00173
00174 set_clusternumber_l(scatter,sph2,0);
00175 naive_listcount_and_analyze(scatter, sph2,N,histogram1, biggestcl1, averagecl1);
00176
00177 cout << "\n-> "<<sph2.size()<<" are in the rest (clustercounted only in this part!):"<< endl;
00178 for (iter2 = sph2.begin(); iter2 != sph2.end(); iter2++){
00179 cout << *iter2 ;
00180 cout << "\t"<< scatter[*iter2] <<endl;
00181 }
00182
00183 NUMLIST e2;
00184 edgespheres(scatter, sph2, e2, GRIDSIZE/2, 2*radius, -1);
00185
00186 cout << "Egdespheres (close enough to the "<<GRIDSIZE/2<<" border) are: ";
00187 for (iter2=e2.begin();iter2 !=e2.end();iter2++){
00188 cout << *iter2 <<" ";
00189 }
00190 cout << endl;
00191
00192 copy_array(scatter, copy1,1,N); // backup the spherearray for testing later
00193 set_clusternumber_l(scatter,sph2,0);
00194 copy_l_count_and_analyze(scatter, copy2,sph2,histogram2, biggestcl2, averagecl2);
00195
00196 if ( ! same_array(copy1,scatter,1,N)) // did the copy_l_count_a.... do the right thing?
00197 errorout ("not the same spheres anymore");
00198
00199 delete[] scatter, copy1,copy2, histogram1, histogram2;
00200 }
00201
00202
00203 void listcount (int dimension=2) {
00204
00205 NUMBER N;
00206 cout <<"Please type in NUMBER of spheres: ";
00207 cin >> N;
00208
00209 sphere *scatter = new sphere[N+1];
00210 set_dim(scatter,0,N,dimension);
00211 sphere *copy1 = new sphere[N+1];
00212 set_dim(copy1,0,N,dimension);
00213 sphere *copy2 = new sphere[N+1];
00214 set_dim(copy2,0,N,dimension);
00215
00216 throw_spheres(scatter, 1, N); // randomly throw spheres from 1 to N in array scatter
00217 set_clustersize(scatter,1,N,0);
00218
00219 double rr=R_critical(N,GRIDSIZE, dimension);
00220 if (rr!=-1) {
00221 cout << "\nThe theoretically critical radius in ";
00222 cout <<dimension<<"dim is "<<rr;
00223 }
00224 else {
00225 rr=R_critical_guessed(N,GRIDSIZE,dimension);
00226 cout << "The critical radius in ";
00227 cout <<dimension<<"dim is approximately "<<rr;
00228 }
00229 cout <<"\nPlease type in the radius of one sphere: ";
00230 REAL R; cin >> R;
00231
00232 set_radius(scatter,1,N,R);
00233
00234 NUMLIST all;
00235 increasing_integers (all, N); // fill the list "all" with the numbers of all spheres
00236
00237 NUMBER biggestcl1, biggestcl2;
00238 REAL averagecl1, averagecl2 ;
00239 NUMBER *histogram1=new NUMBER[(N+1)];
00240 NUMBER *histogram2=new NUMBER[(N+1)];
00241
00242 cout<<"\nfind clusters of spheres given in a list..."<<flush;
00243 set_clusternumber(scatter,1,N,0);
00244 clock_t time1=naive_listcount_and_analyze(scatter, all,N,histogram1, biggestcl1, averagecl1);
00245 cout<<"\nfind clusters of spheres given in an array..."<<flush;
00246 set_clusternumber(scatter,1,N,0);
00247 clock_t time2=naive_count_and_analyze(scatter,1,N,histogram2, biggestcl2, averagecl2);
00248
00249 cout<<"\nfind clusters of spheres given in a list but copy to an array first..."<<flush;
00250 copy_array(scatter, copy1,1,N); // backup the spherearray for testing later
00251 set_clusternumber(scatter,1,N,0);
00252 clock_t time3=copy_l_count_and_analyze(scatter, copy2,all,histogram1, biggestcl1, averagecl1);
00253
00254 cout <<"size\t\tlist-\tnormal-count"<<endl;
00255
00256 show_twohisto(histogram1,histogram2, biggestcl1);
00257
00258 cout <<"\n ****** same resulting histograms like in normalcount? "; // did the analyze do the right thing?
00259 cout <<(sameresult(histogram1, histogram2,1,N) ? "yes":"no")<<" ******\n";
00260 if ((biggestcl1!=biggestcl2)||(averagecl1!=averagecl2))
00261 errorout("different results in biggest or average cluster");
00262 if ( ! same_array(copy1,scatter,1,N)) // did the copy_l_count_a.... do the right thing?
00263 errorout ("not the same spheres anymore");
00264
00265 cout <<"\n time.listcount= "<<ms(time1)<<"\tt.normalcount= "<<ms(time2)<<"\tt.listcopied_normalcount= "<<ms(time3)<<endl;
00266
00267
00268 delete[] scatter, copy1,copy2, histogram1, histogram2;
00269 }
00270
00271
00272
00273
00274 void into_file(int dimension) { // radius { loop {count} }
00275 // generates 1 file line per line (1 line per loop)
00276
00277 clock_t totaltime=clock();
00278 srand( (unsigned)time( NULL ) ); // set the random generator
00279 NUMBER N; // actual number of spheres
00280 int loops, FILE_STEPS; // averaging loops per radius and radius-steps
00281
00282 ask_parameter(N, loops, FILE_STEPS, FILE_FROM, FILE_TO);
00283
00284 sphere *scatter = new sphere[N+1];
00285 set_dim(scatter,0,N,dimension);
00286 sphere *copy = new sphere[N+1];
00287 set_dim(copy,0,N,dimension);
00288 NUMBER *histogram=new NUMBER[(N+1)];
00289
00290 char filename[80];
00291 REAL rfrom, rto, rstep;
00292 set_and_show_filename(filename, FILE_FROM, FILE_TO, FILE_STEPS, loops, GRIDSIZE, dimension, N, FILEHEAD, rfrom, rto, rstep);
00293 cout <<"\n\nIf you want to stop, press <q>"<<endl;
00294 write_intro(filename,N,GRIDSIZE, loops, dimension);
00295
00296 for (double radius =rfrom;radius<=rto;radius=radius+rstep) {
00297 cout <<"*******************************\n r= ";
00298 cout <<radius<<"\t f-factor= "<<flush;
00299 REAL fillingfactor=give_fillingfactor(radius,N,GRIDSIZE, dimension);
00300 cout <<fillingfactor<<endl;
00301
00302 NUMBER biggestcl=0, mostfreqcl=0;
00303 REAL averagecl=0;
00304 clock_t time=0;
00305
00306 REAL sum_biggestcl=0;
00307 REAL sum_averagecl=0;
00308 clock_t sum_time=0;
00309
00310 set_radius(scatter, 1,N,radius);
00311
00312 for (int i=1;i<=loops;i++){
00313 throw_spheres(scatter, 1, N); // randomly throw spheres from 1 to N in array scatter
00314 set_clusternumber(scatter,1,N,0); // no sphere is found yet
00315
00316 time=sort_naivecount_and_analyze(scatter, copy, 1, N, histogram, biggestcl, averagecl);
00317
00318 cout <<"loop "<<i<<" time: "<< ms(time)<<" ms";
00319 cout <<" average size :"<<averagecl;
00320 cout <<" biggest: "<<biggestcl;
00321 cout <<" mst frequent: "<<mostfreqcl<<endl;
00322
00323 sum_time+=time;
00324 sum_biggestcl+=biggestcl;
00325 sum_averagecl+=averagecl;
00326 } // jump to next loop
00327
00328 sum_biggestcl/=loops;
00329 sum_averagecl/=loops;
00330
00331 cout <<loops<<" loops took a time of: "<< ms(sum_time) <<" ms";
00332 cout << "\n averaged average size :"<<sum_averagecl;
00333 cout <<"\n averaged biggest: "<<sum_biggestcl;
00334
00335 write_one_step(filename, radius, fillingfactor, sum_biggestcl, sum_averagecl, sum_time);
00336
00337 //if (_kbhit()) { if (getch()=='q') { closefile( filename , totaltime);exit(0);}}
00338
00339 // jump to next radius
00340 }
00341
00342 closefile( filename, totaltime);
00343
00344
00345 delete[] scatter;
00346 delete[] copy;
00347 delete[] histogram;
00348 }
00349
00350
00351 void into_file2(int dimension) { // loop { all radii {count} }
00352 // # of files = 'loops' (1 file per loop)
00353
00354
00355 srand( (unsigned)time( NULL ) ); // set the random generator
00356 NUMBER N; // actual number of spheres
00357 int loops, FILE_STEPS; // averaging loops per radius and radius-steps
00358
00359 ask_parameter(N, loops, FILE_STEPS, FILE_FROM, FILE_TO);
00360
00361 sphere *scatter = new sphere[N+1];
00362 set_dim(scatter,0,N,dimension);
00363 sphere *copy = new sphere[N+1];
00364 set_dim(copy,0,N,dimension);
00365 NUMBER *histogram=new NUMBER[(N+1)];
00366
00367 char filename[80];
00368 REAL rfrom, rto, rstep;
00369 set_and_show_filename(filename, FILE_FROM, FILE_TO, FILE_STEPS, loops, GRIDSIZE, dimension, N, FILEHEAD,
00370 rfrom, rto, rstep); // radius intervall is set here
00371
00372 NUMBER biggestcl;
00373 REAL averagecl;
00374 REAL fillingfactor;
00375
00376 REAL *sum_biggestcl=new REAL[FILE_STEPS+1]; // these results get better
00377 REAL *sum_averagecl=new REAL[FILE_STEPS+1]; // in each loop (arithmetic mean)
00378 clock_t *sum_time=new clock_t[FILE_STEPS+1]; //
00379
00380 int index;
00381 for (index=0;index<=FILE_STEPS;index++){
00382 sum_biggestcl[index]=0;
00383 sum_averagecl[index]=0;
00384 sum_time[index]=0;
00385 }
00386
00387 clock_t totaltime=clock();
00388
00389 for (int i=1;i<=loops;i++){
00390 cout <<"** Loop "<<i<<"/"<<loops<<" **"<<endl;
00391
00392 index=0;
00393 REAL radius;
00394 for (radius =rfrom;radius<=rto;radius=radius+rstep) {
00395
00396 cout <<" radiusstep= "<<index<<"\t"<<flush;
00397
00398 throw_spheres(scatter, 1, N); // randomly throw spheres from 1 to N in array scatter
00399 set_radius(scatter, 1,N,radius);
00400
00401 biggestcl=0; averagecl=0;
00402
00403 set_clusternumber(scatter,1,N,0); // no sphere is found yet
00404 sum_time[index]+=sort_naivecount_and_analyze(scatter, copy, 1, N, histogram, biggestcl, averagecl);
00405
00406 sum_biggestcl[index]=(sum_biggestcl[index] *(i-1)+biggestcl) /i;
00407 sum_averagecl[index]=(sum_averagecl[index] *(i-1)+averagecl) /i;
00408
00409 index++;
00410
00411 } // next radius
00412
00413
00414 cout <<"\n "<<i<<" loops took total time of: "<< ms(clock()-totaltime)/1000<<" seconds\n"<<endl;
00415
00416 setfilename(filename,FILEHEAD,N,dimension,i, FILE_FROM, FILE_TO, FILE_STEPS);
00417 write_intro(filename,N,GRIDSIZE, i, dimension);
00418
00419 index=0;
00420 for (radius =rfrom;radius<=rto;radius=radius+rstep) {
00421
00422 fillingfactor=give_fillingfactor(radius,N,GRIDSIZE, dimension);
00423 write_one_step(filename, radius, fillingfactor,
00424 sum_biggestcl[index],
00425 sum_averagecl[index],
00426 sum_time[index]);
00427 index++;
00428 }
00429
00430 closefile( filename, totaltime);
00431
00432
00433 } // next loop
00434
00435
00436 delete[] sum_biggestcl;
00437 delete[] sum_averagecl;
00438 delete[] sum_time;
00439
00440 delete[] scatter;
00441 delete[] copy;
00442 delete[] histogram;
00443 }
00444
00445
00446 void test_statistical(){
00447 cout <<"\nHi!\nTest statistics with a list of random numbers.";
00448 cout <<"\nHow many:"<<flush;
00449 int N=0;
00450 while (N<1) cin>>N;
00451 std::list<float> L;
00452 for (int i=0;i<N;i++) L.push_back(rand()/(float)RAND_MAX);
00453 muvarskewkurt<double> A;
00454 int dummy=42;
00455 double dummy2=1.0; // for calculate_moments result type
00456 A=statistics::calculate_moments<double,float>(&manipdata<float>::linear, dummy, L, dummy2);
00457 cout <<"The distribution has these statistical properties:\n";
00458 cout <<"(average mu, variance, skewness, kurtosis)\n" <<A<<endl;
00459 cout <<"Old, simple routine:\n";
00460 double mu=0, sdev=0;
00461 averaging(L, mu, sdev);
00462 cout <<"mu="<<mu<<ASCII_PLUSMINUS<<sdev<<endl;
00463 }
00464
00465 void test_occuring_clusternumbers(){
00466 NUMBER N=200;
00467 int dimension=2;
00468 REAL radius=0.03;
00469 COUNTER cuts=0;
00470 NUMBER numberofcl; NUMBER biggestcl; REAL averagecl;
00471 NUMBER *histogram=new NUMBER[(N+1)];
00472 sphere *spheres = new sphere[N+1]; set_dim(spheres,0,N,dimension);
00473 NUMLIST all;
00474
00475 increasing_integers (all, N); // fill the list "all" with the numbers of all spheres
00476 throw_spheres(spheres, 1, N); // randomly throw spheres from 1 to N in array spheres
00477 set_radius(spheres, 1,N,radius);
00478
00479 set_clusternumber(spheres,1,N,0); // no sphere is found yet (initialization for counter)
00480 numberofcl=divide_by_c_cuts_count_and_analyze(spheres, all, cuts,histogram, biggestcl, averagecl);
00481
00482 NUMLIST clusternumbers;
00483 occuring_clusternumbers(spheres, all, clusternumbers);
00484
00485 cout <<"\n# of clusters: "<<numberofcl<<" biggestcl: "<<biggestcl<<" meancl:"<<averagecl<<"\noccuring clusternumbers:\n";
00486 NUMLIST::iterator clno;
00487 for (clno=clusternumbers.begin();clno!=clusternumbers.end();clno++){
00488 cout<<*clno<<"\t";
00489 }
00490 cout <<"\n";
00491 }
00492
00493
00494 void into_file7_pt1(int dimension, NUMBER N,
00495 COUNTER min_loops, COUNTER additional_loops,
00496 COUNTER ff_steps, REAL ff_from, REAL ff_to)
00497 {
00498 clock_t starttime=clock();
00499
00500 char filename[80];
00501 REAL rfrom, rto, rstep;
00502 rfrom=give_radius(ff_from, N, GRIDSIZE, dimension);
00503 rto=give_radius(ff_to, N, GRIDSIZE, dimension);
00504 rstep=(rto-rfrom)/ff_steps;
00505 COUNTER cuts=choose_optimal_cuts(dimension, N, (rto+rfrom)/2);
00506
00507 // initialize and zero all measurement containers
00508 sphere *scatter = new sphere[N+1];
00509 set_dim(scatter,0,N,dimension);
00510 NUMLIST all;
00511 increasing_integers (all, N); // fill the list "all" with the numbers of all spheres
00512 NUMBER *histogram=new NUMBER[(N+1)];
00513 NUMBER numberofcl; NUMBER biggestcl; REAL averagecl;
00514
00515 REAL *R=new REAL[ff_steps+1]; // radius
00516 REAL *ffactor=new REAL[ff_steps+1]; // fillingfactor
00517
00518 NUMBER *this_biggestcl=new NUMBER[ff_steps+1]; // these two are the results
00519 REAL *this_averagecl=new REAL[ff_steps+1]; // of the last run
00520 REAL *sum_biggestcl=new REAL[ff_steps+1]; // these four results get better
00521 REAL *sum_square_biggestcl=new REAL[ff_steps+1]; // in each loop (arithmetic mean)
00522 REAL *sum_averagecl=new REAL[ff_steps+1]; //
00523 REAL *sum_square_averagecl=new REAL[ff_steps+1]; //
00524 REAL *fluctuation_biggestcl=new REAL[ff_steps+1]; // fluctuation / susceptibility of B.cl and Av.cl
00525 REAL *fluctuation_averagecl=new REAL[ff_steps+1]; // = <B^2> - <B>^2
00526
00527
00528 NUMBER *this_numberofcl=new NUMBER[ff_steps+1]; // |
00529 REAL *sum_numberofcl=new REAL[ff_steps+1]; // | the same as above for
00530 REAL *sum_square_numberofcl=new REAL[ff_steps+1]; // | M0 - the total number of clusters
00531 REAL *fluctuation_numberofcl=new REAL[ff_steps+1]; // |
00532
00533 COUNTER *loop_no=new COUNTER[ff_steps+1]; // this is the number of averaging loops for this radius
00534 clock_t *sum_time=new clock_t[ff_steps+1]; // this contains the totaltime for this radius
00535
00536 COUNTER index;
00537 for (index=0;index<=ff_steps;index++){
00538 sum_numberofcl[index]=0;
00539 sum_square_numberofcl[index]=0;
00540 sum_biggestcl[index]=0;
00541 sum_square_biggestcl[index]=0;
00542 sum_averagecl[index]=0;
00543 sum_square_averagecl[index]=0;
00544 sum_time[index]=0;
00545 loop_no[index]=0;
00546 }
00547
00548 index=0;
00549 REAL radius;
00550 for (radius =rfrom;radius<=rto;radius=radius+rstep) {
00551 R[index]=radius;
00552 ffactor[index]=give_fillingfactor(radius,N,GRIDSIZE, dimension);
00553 index++;
00554 }
00555
00556 REAL min_fluctuation=0;
00557 REAL ffactor_min_fluct;
00558 REAL max_fluctuation=0;
00559 REAL ffactor_max_fluct;
00560
00561 // here starts every loop
00562 // for the first averaging loops (minimum)
00563 COUNTER loop;
00564 for (loop=1;loop<=(min_loops);loop++){
00565 cout <<"** Loop "<<loop<<"/"<<min_loops<<"(min) ; ( N="<<N<<" dim="<<dimension<<" cuts="<<cuts<<") **"<<endl;
00566
00567 // here starts every radiusstep
00568 index=0;
00569 for (radius =rfrom;radius<=rto;radius=radius+rstep) {
00570
00571 cout <<" radiusstep= "<<index<<" (f-factor="<<ffactor[index]<<")\t"<<flush;
00572
00573 // throw new spheres and count the clusters
00574 throw_spheres(scatter, 1, N); // randomly throw spheres from 1 to N in array scatter
00575 set_radius(scatter, 1,N,radius);
00576 set_clusternumber(scatter,1,N,0); // no sphere is found yet (initialization for counter)
00577 clock_t onerun=clock();
00578 numberofcl=divide_by_c_cuts_count_and_analyze(scatter, all, cuts,histogram, biggestcl, averagecl);
00579 sum_time[index]+=(clock()-onerun);
00580
00581 // store and calculate the results (biggest cluster size & average cluster size)
00582 this_numberofcl[index]=numberofcl;
00583 this_biggestcl[index]=biggestcl;
00584 this_averagecl[index]=averagecl;
00585
00586 sum_numberofcl[index]=(sum_numberofcl[index] *((REAL)(loop_no[index]))+numberofcl) / (loop_no[index]+1);
00587 sum_biggestcl[index]=(sum_biggestcl[index] *((REAL)(loop_no[index]))+biggestcl) / (loop_no[index]+1);
00588 sum_averagecl[index]=(sum_averagecl[index] *((REAL)(loop_no[index]))+averagecl) / (loop_no[index]+1);
00589
00590 sum_square_numberofcl[index]=(sum_square_numberofcl[index] *((REAL)(loop_no[index]))+ pow(numberofcl,2) ) / (loop_no[index]+1);
00591 sum_square_biggestcl[index]=(sum_square_biggestcl[index] *((REAL)(loop_no[index]))+ pow(biggestcl,2) ) / (loop_no[index]+1);
00592 sum_square_averagecl[index]=(sum_square_averagecl[index] *((REAL)(loop_no[index]))+ pow(averagecl,2) ) / (loop_no[index]+1);
00593
00594 fluctuation_numberofcl[index]=sum_square_numberofcl[index] - pow (sum_numberofcl[index],2);
00595 fluctuation_biggestcl[index]=sum_square_biggestcl[index] - pow (sum_biggestcl[index],2);
00596 fluctuation_averagecl[index]=sum_square_averagecl[index] - pow (sum_averagecl[index],2);
00597 loop_no[index]++;
00598
00599 index++;
00600 } // next radius
00601
00602 find_table_minmax(0, ff_steps, // find the smallest and biggest fluctuation
00603 ffactor, fluctuation_biggestcl,
00604 min_fluctuation, ffactor_min_fluct,
00605 max_fluctuation, ffactor_max_fluct);
00606
00607 cout <<"\n "<<loop<<" loops took total time of: "<< ms(clock()-starttime)/1000<<" seconds"<<endl;
00608 cout <<" Biggest fluctuation at: ffactor="<<ffactor_max_fluct;
00609 cout <<"\n Now saving all results...\n"<<endl;
00610
00611 // save to network harddisc:
00612 setfilename(filename,FILEHEAD1,N,dimension,loop, ff_from, ff_to, ff_steps);
00613 write_intro7(filename,N,GRIDSIZE, loop, dimension, cuts);
00614 write_results7 (filename, ff_steps,
00615 R, ffactor,
00616 sum_numberofcl, sum_biggestcl, sum_averagecl,
00617 fluctuation_numberofcl, fluctuation_biggestcl, fluctuation_averagecl,
00618 this_numberofcl, this_biggestcl, this_averagecl,
00619 loop_no,
00620 sum_time,
00621 clock()-starttime,
00622 min_fluctuation, ffactor_min_fluct,
00623 max_fluctuation, ffactor_max_fluct);
00624 // local copy for security reasons:
00625 setfilename(filename,FILEHEAD2,N,dimension,loop, ff_from, ff_to, ff_steps);
00626 write_intro7(filename,N,GRIDSIZE, loop, dimension, cuts);
00627 write_results7 (filename, ff_steps,
00628 R, ffactor,
00629 sum_numberofcl, sum_biggestcl, sum_averagecl,
00630 fluctuation_numberofcl, fluctuation_biggestcl, fluctuation_averagecl,
00631 this_numberofcl, this_biggestcl, this_averagecl,
00632 loop_no,
00633 sum_time,
00634 clock()-starttime,
00635 min_fluctuation, ffactor_min_fluct,
00636 max_fluctuation, ffactor_max_fluct);
00637
00638 } // next loop
00639
00640
00641 // now continue doing the same, but simulation
00642 // only for the radii that are fluctuating most
00643 // (these are the "additional loops")
00644
00645 cout <<"\n\n-------------------------------------------------------------------"<<endl;
00646 cout <<"Additional Loops following:"<<endl;
00647 cout <<"Now continue only for the fillingfactors/radii, that fluctuate most.";
00648 cout <<"\nSo the results will be better around the critical value!"<<endl;
00649 cout <<"-------------------------------------------------------------------"<<endl;
00650
00651 REAL fluct_threshold; // if the (one loop before) fluctuation is higher, do another simulation for this radius
00652 REAL estimate_fluctuation;
00653 COUNTER neighbourhood;
00654
00655 for (loop=min_loops+1;loop<=(min_loops + additional_loops);loop++){
00656 cout <<"** Loop "<<loop<<"/"<<min_loops+additional_loops<<"(max) ; ( N="<<N<<" dim="<<dimension<<" cuts="<<cuts<<") **"<<endl;
00657
00658 fluct_threshold= max_fluctuation * ((REAL)loop / (REAL)(min_loops+additional_loops) ) ;
00659
00660 // here starts every radiusstep
00661 index=0;
00662 for (radius =rfrom;radius<=rto;radius=radius+rstep) {
00663
00664 //TODO: find a good value for this ( depending on rough / medium /fine radius-stepsize)
00665 neighbourhood=1; // average over 1*2+1=3 points
00666 estimate_fluctuation=neighbourhood_average (fluctuation_biggestcl, index, ff_steps,neighbourhood);
00667
00668 if (estimate_fluctuation >= fluct_threshold)
00669 {
00670 cout <<" radiusstep= "<<index<<" (f-factor="<<ffactor[index]<<")\t"<<flush;
00671
00672 // throw new spheres and count the clusters
00673 throw_spheres(scatter, 1, N); // randomly throw spheres from 1 to N in array scatter
00674 set_radius(scatter, 1,N,radius);
00675 set_clusternumber(scatter,1,N,0); // no sphere is found yet (initialization for counter)
00676 clock_t onerun=clock();
00677 numberofcl=divide_by_c_cuts_count_and_analyze(scatter, all, cuts,histogram, biggestcl, averagecl);
00678 sum_time[index]+=(clock()-onerun);
00679
00680 // store and calculate the results (biggest cluster size & average cluster size)
00681
00682 this_numberofcl[index]=numberofcl;
00683 this_biggestcl[index]=biggestcl;
00684 this_averagecl[index]=averagecl;
00685
00686 // repaired mistake (!) : not divide & multiply by loop, but by loop_no[index]
00687
00688 sum_numberofcl[index]=(sum_numberofcl[index] *((REAL)(loop_no[index]))+numberofcl) / (loop_no[index]+1);
00689 sum_biggestcl[index]=(sum_biggestcl[index] *((REAL)(loop_no[index]))+biggestcl) / (loop_no[index]+1);
00690 sum_averagecl[index]=(sum_averagecl[index] *((REAL)(loop_no[index]))+averagecl) / (loop_no[index]+1);
00691
00692 sum_square_numberofcl[index]=(sum_square_numberofcl[index] *((REAL)(loop_no[index]))+ pow(numberofcl,2) ) / (loop_no[index]+1);
00693 sum_square_biggestcl[index]=(sum_square_biggestcl[index] *((REAL)(loop_no[index]))+ pow(biggestcl,2) ) / (loop_no[index]+1);
00694 sum_square_averagecl[index]=(sum_square_averagecl[index] *((REAL)(loop_no[index]))+ pow(averagecl,2) ) / (loop_no[index]+1);
00695
00696 fluctuation_numberofcl[index]=sum_square_numberofcl[index] - pow (sum_numberofcl[index],2);
00697 fluctuation_biggestcl[index]=sum_square_biggestcl[index] - pow (sum_biggestcl[index],2);
00698 fluctuation_averagecl[index]=sum_square_averagecl[index] - pow (sum_averagecl[index],2);
00699
00700 loop_no[index]++;
00701 }
00702
00703 index++;
00704 } // next radius
00705
00706 find_table_minmax(0, ff_steps, // find the smallest and biggest fluctuation
00707 ffactor, fluctuation_biggestcl,
00708 min_fluctuation, ffactor_min_fluct,
00709 max_fluctuation, ffactor_max_fluct);
00710
00711 cout <<"\n "<<loop<<" loops took total time of: "<< ms(clock()-starttime)/1000<<" seconds"<<endl;
00712 cout <<" Biggest fluctuation at: ffactor="<<ffactor_max_fluct;
00713 cout <<"\n Now saving all results...\n"<<endl;
00714
00715 // save to network harddisc:
00716 setfilename(filename,FILEHEAD1,N,dimension,loop, ff_from, ff_to, ff_steps);
00717 write_intro7(filename,N,GRIDSIZE, loop, dimension, cuts);
00718 write_results7 (filename, ff_steps,
00719 R, ffactor,
00720 sum_numberofcl, sum_biggestcl, sum_averagecl,
00721 fluctuation_numberofcl, fluctuation_biggestcl, fluctuation_averagecl,
00722 this_numberofcl, this_biggestcl, this_averagecl,
00723 loop_no,
00724 sum_time,
00725 clock()-starttime,
00726 min_fluctuation, ffactor_min_fluct,
00727 max_fluctuation, ffactor_max_fluct);
00728 // local copy for security reasons:
00729 setfilename(filename,FILEHEAD2,N,dimension,loop, ff_from, ff_to, ff_steps);
00730 write_intro7(filename,N,GRIDSIZE, loop, dimension, cuts);
00731 write_results7 (filename, ff_steps,
00732 R, ffactor,
00733 sum_numberofcl, sum_biggestcl, sum_averagecl,
00734 fluctuation_numberofcl, fluctuation_biggestcl, fluctuation_averagecl,
00735 this_numberofcl, this_biggestcl, this_averagecl,
00736 loop_no,
00737 sum_time,
00738 clock()-starttime,
00739 min_fluctuation, ffactor_min_fluct,
00740 max_fluctuation, ffactor_max_fluct);
00741 } // next loop
00742
00743
00744
00745 delete[] R, ffactor;
00746 delete[] this_numberofcl, this_averagecl, this_biggestcl;
00747 delete[] sum_numberofcl, sum_square_numberofcl ;
00748 delete[] sum_biggestcl, sum_square_biggestcl ;
00749 delete[] sum_averagecl, sum_square_averagecl ;
00750 delete[] fluctuation_numberofcl, fluctuation_biggestcl, fluctuation_averagecl;
00751 delete[] sum_time;
00752 delete[] loop_no;
00753
00754 delete[] scatter;
00755 delete[] histogram;
00756 }
00757
00758 // END of file7
00759 //////////////////////////////
00760
00761
00762 void find_member_minmax(COUNTER tablestart, COUNTER tableend,
00763 REAL * ffactor, all_results *res,
00764 muvarskewkurt<REAL> (all_results::* parameter),
00765 measure<REAL> (muvarskewkurt<REAL>::* moment),
00766 REAL &min_x, REAL &min_y, REAL &max_x, REAL &max_y)
00767 {
00768 REAL value=((res[tablestart].*parameter).*moment).av();
00769 min_y=value;
00770 min_x=ffactor[tablestart];
00771 max_y=value;
00772 max_x=ffactor[tablestart];
00773 for (COUNTER index=tablestart;index<=tableend;index++) {
00774 value=((res[index].*parameter).*moment).av();
00775 if ( value < min_y) {
00776 min_y=value;
00777 min_x=ffactor[index];
00778 }
00779 if ( value > max_y) {
00780 max_y=value;
00781 max_x=ffactor[index];
00782 }
00783 }
00784 }
00785
00786
00787 void into_file8_pt1(int dim, NUMBER N,
00788 COUNTER min_loops, COUNTER additional_loops,
00789 COUNTER ff_steps, REAL ff_from, REAL ff_to)
00790 {
00791
00792 // TODO: clock() -> 1000*clock()/CLOCKS_PER_SEC
00793 //
00794 // TODO: fluctuation threshold auch für mean_clustersize!
00795
00796 clock_t starttime=clock();
00797
00798 // initialize sphere container
00799 sphere *spheres = new sphere[N+1];
00800 if (spheres==NULL) exit_out_of_memory("into_file8_pt1(...): sphere *spheres");
00801 set_dim(spheres,0,N,dim);
00802 NUMLIST all; increasing_integers (all, N); // fill the list "all" with the numbers of all spheres
00803
00804 all_results *res = new all_results[ff_steps+1]; // result container-table
00805 if (res==NULL) exit_out_of_memory("into_file8_pt1(...): all_results *res");
00806 set_N_and_dim(res,ff_steps,N,dim);
00807
00808 // TODO: put R and ffactor into all_results::all_results(N,dim)
00809 // !!! but will it be needed elsewhere???
00810 // then: make_?_table -> set_N_and_dim..._and_R_and_ff(...,ff_from,ff_to_,ff_steps)
00811 REAL *R=new REAL[ff_steps+1]; // radius-table
00812 if (R==NULL) exit_out_of_memory("into_file8_pt1(...): REAL *R");
00813 REAL *ffactor=new REAL[ff_steps+1]; // fillingfactor-table
00814 if (R==NULL) exit_out_of_memory("into_file8_pt1(...): REAL *ffactor");
00815
00816 // TODO: make_linear_table -> R
00817 // TODO: make_ff_table -> ffactor
00818 REAL rfrom, rto, rstep;
00819 rfrom=give_radius(ff_from, N, GRIDSIZE, dim);
00820 rto =give_radius(ff_to, N, GRIDSIZE, dim);
00821 rstep=(rto-rfrom)/ff_steps;
00822
00823 COUNTER cuts=choose_optimal_cuts(dim, N, (rto+rfrom)/2);
00824
00825 REAL radius=rfrom;
00826 COUNTER index;
00827 for (index=0;index<=ff_steps;index++) {
00828 R[index]=radius;
00829 ffactor[index]=give_fillingfactor(radius,N,GRIDSIZE, dim);
00830 radius+=rstep;
00831 }
00832
00833 // temporary variables
00834 REAL min_fluctuation=0; REAL ffactor_min_fluct;
00835 REAL max_fluctuation=0; REAL ffactor_max_fluct;
00836 clock_t clcount_time, spcl_search_time;
00837 one_result step;
00838
00839 // here starts every loop
00840 // for the first averaging loops (minimum)
00841 COUNTER loop;
00842 for (loop=1;loop<=(min_loops);loop++){
00843 cout <<"** Loop "<<loop<<"/"<<min_loops<<"(min) ; ( N="<<N;
00844 cout <<" dim="<<dim<<" cuts="<<cuts<<") **"<<endl;
00845
00846 radius=rfrom;
00847 for (index=0;index<=ff_steps;index++) {
00848
00849 cout <<" radiusstep= "<<index<<" (f-factor=";
00850 cout <<give_fillingfactor(radius,N,GRIDSIZE, dim)<<")\t"<<flush;
00851
00852 // here is the cluster-counter
00853 step=throw_dnc_count_analyze_step(spheres, all, radius, clcount_time, spcl_search_time);
00854
00855 res[index].add_result(step);
00856 res[index].compute_averages();
00857 res[index].sum_time += clcount_time;
00858 res[index].sum_spsearch_time += spcl_search_time;
00859
00860 radius+=rstep;
00861 } // next radius
00862
00863 find_member_minmax(0, ff_steps, ffactor, res,
00864 &all_results::spcl_cum, // of the spanningcluster cumulant
00865 &muvarskewkurt<REAL>::variance, // the variance
00866 ffactor_min_fluct, min_fluctuation,
00867 ffactor_max_fluct, max_fluctuation);
00868 cout <<"\n "<<loop<<" loops took total time of: "<< ms(clock()-starttime)*1000<<" seconds"<<endl;
00869 cout <<" Biggest fluctuation at: ffactor="<<ffactor_max_fluct;
00870 cout <<"\n Now saving all results...\n"<<endl;
00871 save_results8(dim,N,loop,ff_from,ff_to,ff_steps,cuts,
00872 R,ffactor,res,
00873 starttime,
00874 min_fluctuation,ffactor_min_fluct,max_fluctuation, ffactor_max_fluct);
00875 } // next loop
00876
00877
00878 // now continue doing the same, but simulation
00879 // only for the radii that are fluctuating most
00880 // (these are the "additional loops")
00881
00882 cout <<"\n\n-------------------------------------------------------------------"<<endl;
00883 cout <<"Additional Loops following:"<<endl;
00884 cout <<"Now continue only for the fillingfactors/radii, that fluctuate most.";
00885 cout <<"\nSo the results will be better around the critical value!"<<endl;
00886 cout <<"-------------------------------------------------------------------"<<endl;
00887
00888 REAL fluct_threshold; // if the (one loop before) fluctuation is higher, do another simulation for this radius
00889 REAL fluctuation;
00890
00891 for (loop=min_loops+1;loop<=(min_loops + additional_loops);loop++){
00892 cout <<"** Loop "<<loop<<"/"<<min_loops+additional_loops<<"(add) ; ( N="<<N;
00893 cout <<" dim="<<dim<<" cuts="<<cuts<<") **"<<endl;
00894
00895 fluct_threshold= max_fluctuation * ((REAL)loop / (REAL)(min_loops+additional_loops) ) ;
00896
00897 radius=rfrom;
00898 for (index=0;index<=ff_steps;index++) {
00899 fluctuation=res[index].spcl_cum.variance.av();
00900
00901 if (fluctuation >= fluct_threshold)
00902 {
00903 cout <<" radiusstep= "<<index<<" (f-factor=";
00904 cout <<give_fillingfactor(radius,N,GRIDSIZE, dim)<<")\t"<<flush;
00905
00906 // here is the cluster-counter
00907 step=throw_dnc_count_analyze_step(spheres, all, radius, clcount_time, spcl_search_time);
00908
00909 res[index].add_result(step);
00910 res[index].compute_averages();
00911 res[index].sum_time += clcount_time;
00912 res[index].sum_spsearch_time += spcl_search_time;
00913 }
00914 else res[index].no_result(); // don't print this in resultfile
00915
00916 radius+=rstep;
00917 } // next radius
00918
00919 find_member_minmax(0, ff_steps, ffactor, res,
00920 &all_results::spcl_cum, // of the spanningcluster cumulant
00921 &muvarskewkurt<REAL>::variance, // the variance
00922 ffactor_min_fluct, min_fluctuation,
00923 ffactor_max_fluct, max_fluctuation);
00924 cout <<"\n "<<loop<<" loops took total time of: "<< ms(clock()-starttime)/1000<<" seconds\n";
00925 cout <<" Biggest fluctuation at: ffactor="<<ffactor_max_fluct;
00926 cout <<"\n Now saving all results...\n"<<endl;
00927 save_results8(dim,N,loop,ff_from,ff_to,ff_steps,cuts,
00928 R,ffactor,res,
00929 starttime,
00930 min_fluctuation,ffactor_min_fluct,max_fluctuation, ffactor_max_fluct);
00931 } // next loop
00932
00933 delete[] R, ffactor;
00934 delete[] res;
00935 delete[] spheres;
00936 }
00937
00938
00939 void into_file7_starter(int dimension){
00940 srand( SEED4 ); // set the random generator to seed 4
00941 NUMBER N; // actual number of spheres
00942 COUNTER min_loops, additional_loops, ff_steps; // averaging loops per radius and radius-steps
00943 REAL ff_from, ff_to;
00944 ask_for_ff_range(ff_from, ff_to);
00945 ask_for_parameters(N, min_loops, additional_loops, ff_steps, ff_from, ff_to);
00946 into_file7_pt1(dimension, N, min_loops, additional_loops, ff_steps, ff_from, ff_to);
00947 }
00948
00949 void into_file8_starter(int dimension){
00950 srand( SEED4 ); // set the random generator to seed 4
00951 NUMBER N; // actual number of spheres
00952 COUNTER min_loops, additional_loops, ff_steps; // averaging loops per radius and radius-steps
00953 REAL ff_from, ff_to;
00954 ask_for_ff_range(ff_from, ff_to);
00955 ask_for_parameters(N, min_loops, additional_loops, ff_steps, ff_from, ff_to);
00956
00957 into_file8_pt1(dimension, N, min_loops, additional_loops, ff_steps, ff_from, ff_to);
00958 }
00959
00960
00961
00962
00963
00964 void unitspheredimensions(int first, int last) {
00965 FILE *outFile;
00966 outFile = fopen( FILENAME1, "w" );
00967
00968 cout <<"volume of unitsphere in n dimensions. Write results to "<<FILENAME1<<endl;
00969 waitanykey();
00970
00971 for (int i=first; i<=last ;i++) {
00972 cout <<i<<"\t"<<unitsphere(i)<<endl;
00973 fprintf(outFile,FORMAT,i); fprintf(outFile,"%s"," ");
00974 fprintf(outFile,FORMATEXACT,unitsphere(i)); fprintf(outFile,"%c",'\n');
00975 }
00976 fclose( outFile ) ;
00977 }
00978
00979 void test_saturating_fillingfactor(){
00980 cout <<"\ndim\tN\tffs_guessed(dim,N)";
00981 NUMBER N;
00982 for(int dim=1;dim<=11;dim++){
00983 for (REAL floatingpointN=39.0625;floatingpointN<=10000;floatingpointN*=2){
00984 N=(NUMBER)floatingpointN;
00985 cout<<"\n"<<dim<<"\t"<<N<<"\t"<<saturating_fillingfactor(dim,N)<<flush;
00986 }
00987 }
00988 cout<<"\n"<<endl;
00989 }
00990
00991
00992 void max_ff(int startdim, int enddim){
00993 cout <<"\nIf there are two touching spheres in opposite corners (which have a radius that is half of the diagonal of the space), we get the fillingfactor that forces percolation (a 'maximal fillingfactor'):\n"<<endl;
00994
00995 cout <<"Write results to "<<FILENAME2<<endl;
00996 FILE *outFile;
00997 outFile = fopen( FILENAME2, "w" );
00998 if (!outFile) cerr<<"Couldn't open file "<<FILENAME2<<"\n";
00999
01000 int dim;
01001 for(dim=startdim;dim<=enddim;dim++){
01002 cout << "dim= "<<setw(2)<<dim<<": ff.max= ";
01003 cout << percolating_fillingfactor_for_two_spheres (dim)<<"\t";
01004 cout << percolating_fillingfactor_for_two_spheres_verbose (dim, 42)<<"\n";
01005 fprintf(outFile,FORMAT,dim);
01006 fprintf(outFile,"%s","\t");
01007 fprintf(outFile,FORMATEXACT,percolating_fillingfactor_for_two_spheres(dim));
01008 fprintf(outFile,"%c",'\n');
01009 }
01010 fclose(outFile);
01011 }
01012
01013 void test_ffnoc1percent(string filename){
01014 cout <<"ffnoc1% fitfunction\nsaved to "<<filename<<"\n\n";
01015 FILE* outFile= fopen( filename.c_str(), "w" );
01016 int dim; NUMBER N; double NN; double ffnoc1percent;
01017 for (dim=1;dim<=11;dim++){
01018 for (NN=156.25; NN<=320000;NN*=2){
01019 N=(NUMBER)NN;
01020 ffnoc1percent=fillingfactor_noc1percent(dim,N);
01021 fprintf(outFile,"%d\t",dim);
01022 fprintf(outFile,"%d\t",N);
01023 fprintf(outFile,"%f\n",ffnoc1percent);
01024 cout<<dim<<"\t"<<N<<"\t"<<ffnoc1percent<<"\n";
01025 }
01026 }
01027 fclose(outFile);
01028 }
01029
01030 } // end of namespace starters
| Diploma Thesis Sourcecode
Documentation check out the text and the executable binaries |