| Diploma Thesis Percolation Simulation C++ Sourcecode Documentation |
00001
00002 // divide_and_conquer.h
00003 // divides spheres in parts, solves the partial problems
00004 // and combines the cluster
00005 //
00006 // Andreas Krueger 14.7.2000
00007 // v2.0 , last changes: 26.8.2000
00008 //
00009
00010 namespace counters {
00011 using analyze::edgespheres;
00012 using analyze::analyze_clusters;
00013 using analyze::erasehisto;
00014 using analyze::makehistogram;
00015
00016
00017 // divides list original->lowerlist+upperlist,
00018 // depending on value 'border' in coordinate 'direction'
00019 void divide_all_spheres (sphere *array,
00020 NUMLIST & original,
00021 NUMLIST & lowerlist,
00022 NUMLIST & upperlist,
00023 COORDFLOAT border,
00024 int direction){
00025
00026 NUMLIST::iterator iter1;
00027
00028 for (iter1 = original.begin(); iter1 != original.end(); ++iter1){
00029 if ( lower(array[*iter1].c, border, direction))
00030 lowerlist.insert(lowerlist.end(), *iter1);
00031 else upperlist.insert(upperlist.end(), *iter1);
00032 }
00033 }
00034
00035
00036
00037 #ifndef TIMEMEASUREMENT
00038
00039 NUMBER combine (sphere* spheres,
00040 NUMLIST &sph1,
00041 NUMLIST &sph2,
00042 cluson* clst,
00043 NUMBER firstclno,
00044 NUMBER c1,
00045 NUMBER c2,
00046 COORDFLOAT border,
00047 int direction)
00048 {
00049 if (sph1.empty() || sph2.empty()) return c2; // if one of the parts is empty,
00050 // do nothing
00051
00052 sphere s1=spheres[*sph1.begin()]; sphere s2=spheres[*sph2.begin()];
00053 sphere temp=s1; // for a faster overlap-check
00054
00055 if (s1.clno > s2.clno ) {
00056 errorout("combine is not possible, \nbecause clusno's in first list greater than in second list");
00057 cout <<s1.clno <<"\t"<< s2.clno<<endl;
00058 exit (1);
00059 }
00060 if (lower(s2.c-s1.c, 0,direction)) {
00061 errorout("combine is not possible, \nbecause coordinates of first list are higher than of second list");
00062 cout <<"(in direction "<<direction<<" ): s1.c="<<s1.c<<" s2.c="<<s2.c<<endl;
00063 exit(1);
00064 }
00065
00066 REAL radius1=find_biggestradius(spheres,sph1);
00067 REAL radius2=find_biggestradius(spheres,sph2);
00068 REAL radius=( (radius1 > radius2) ? radius1 : radius2);
00069
00070 NUMLIST e1, e2;
00071 edgespheres(spheres, sph1, e1,border, 2*radius, direction);
00072 edgespheres(spheres, sph2, e2,border, 2*radius, -direction);
00073 // cout <<"number of edgespheres: "<< e1.size()<<"\t"<<e2.size()<<endl;
00074
00075 // now find all overlapping edgespheres -> clstoverlap[ ]
00076 // = list of directly overlapping clno's (both ways "forwards" and "backwards"!)
00077 // !!! one clno might be mentioned several times !!!
00078
00079 clusterneighbours *clstoverlap = new clusterneighbours[c2-firstclno];
00080 // for the clusternumbers [firstclno;...c1.............;c2-1]
00081 // use the (smaller) array [0;...........c1-firstclno...;c2-firstclno-1]
00082
00083 NUMLIST::iterator esph1, esph2;
00084
00085 // idea: sort first perpendicular to cutline,
00086 // then compare only those within a certain distance
00087 // but at the moment, this checks _all_ edgespheres:
00088 for (esph1=e1.begin();esph1!=e1.end();esph1++){
00089 for (esph2=e2.begin();esph2!=e2.end();esph2++){
00090 s1= spheres[*esph1];
00091 s2= spheres[*esph2];
00092 if (overlap2 (s1,s2,temp)) { // overlapping spheres -> overlapping clusters
00093 clstoverlap[s1.clno-firstclno].clusno.push_back ( s2.clno); // forwards
00094 clstoverlap[s2.clno-firstclno].clusno.push_back ( s1.clno); // backwards
00095 }
00096 }
00097 }
00098
00099 e1.clear(); e2.clear();
00100
00101 // now follow all cluster-connections in clstoverlap[]
00102 // delete the multiple entries
00103 // and collect all clno's of overlapping clusters
00104 // in the list with the smallest clusternumber
00105
00106 NUMLIST::iterator nneighb ;
00107 NUMBER clno;
00108 for (clno=firstclno;clno<=c1-1;clno++){
00109 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno; // alias (!)
00110
00111 if (! neighbours.empty() ){ // no overlap
00112 if ( *(neighbours.begin()) != -1) // already counted
00113 {
00114 for (nneighb=neighbours.begin();nneighb!=neighbours.end();nneighb++){
00115
00116 while ( (*nneighb==clno) || (*(clstoverlap[*nneighb-firstclno].clusno.begin())==-1 )) {
00117 nneighb=neighbours.erase(nneighb);
00118 if (nneighb==neighbours.end()) break;
00119 }
00120 if (nneighb==neighbours.end()) break;
00121
00122 neighbours.splice(neighbours.end(), clstoverlap[*nneighb-firstclno].clusno);
00123 clstoverlap[*nneighb-firstclno].clusno.push_front(-1);
00124 }
00125 }
00126 else neighbours.clear();
00127 }
00128 }
00129
00130
00131 /*
00132 // the following is not necessary if the tests below are switched off
00133 // because the combination is done in [firstclno;c1-1] only
00134 //
00135 // clear the backwards connections (that were set to -1 in the above process)
00136 for (clno=c1;clno<=c2-1;clno++){
00137 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno;
00138 if ( neighbours.size() != 0 && *(neighbours.begin()) == -1 ) {
00139 neighbours.clear();
00140 }
00141 }
00142
00143 // only a test (time-consuming!):
00144 // test if all self-references and "-1" have been erased
00145 NUMBER s;
00146 for (clno=firstclno;clno<=c2-1;clno++){
00147 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno;
00148 if ( neighbours.size()!=0 && (*(neighbours.begin())) == -1) {
00149 error("ALARM: -1 left");
00150 cout << "clno="<<clno<<endl;
00151 cout << "firstclno="<<firstclno<<" c1="<<c1<<" c2="<<c2<<endl;
00152 exit(1);
00153 }
00154 s=neighbours.size();
00155 neighbours.remove(clno);
00156 s=s-neighbours.size();
00157 if (s !=0) {
00158 error("still self-references in clusterlist");
00159 cout << " clno="<<clno<<" # of self references=" <<s<<endl;
00160 exit(1);
00161 }
00162 }
00163
00164 // only a test (time consuming!)
00165 // test if all the clusters have been combined
00166 // by looking for non-empty backwards connection lists
00167 for (clno = c1 ; clno<= c2-1;clno++){
00168 if (clstoverlap[clno-firstclno].clusno.size()!= 0) {
00169 error ("not found every connection during combining");
00170 cout <<"clno="<<clno;
00171 cout <<" overlap["<<clno<<"].clusno.size = "<<clstoverlap[clno-firstclno].clusno.size()<<endl;
00172 NUMLIST::iterator cl;
00173 for (cl=clstoverlap[clno-firstclno].clusno.begin();cl!=clstoverlap[clno-firstclno].clusno.end();cl++){
00174 cout <<*cl<< " ";
00175 }
00176 cout <<endl;
00177 exit(1);
00178 }
00179 }
00180
00181
00182 // show the resulting cluster combination
00183 cout <<"\ncombined clusters are"<<endl;
00184 for (clno=firstclno;clno<=c1-1;clno++){
00185
00186 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno;
00187 if (! neighbours.empty() ){
00188 cout <<"{"<< clno<<"+ ";
00189 for (nneighb=neighbours.begin();nneighb!=neighbours.end();nneighb++){
00190 cout <<*nneighb<< " "<<flush;
00191 }
00192 cout <<"}\t"<<flush;
00193 }
00194 }
00195 */
00196
00197 // now combine the clusters of spheres according to clstoverlap
00198 // and set the right clusternumber to each sphere
00199
00200 // idea: perhaps renumbering doesn't have to be done here, but below
00201
00202 for (clno=firstclno;clno<=c1-1;clno++){
00203 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno;
00204 if (! neighbours.empty() ){
00205 for (nneighb=neighbours.begin();nneighb!=neighbours.end();nneighb++){
00206 set_clusternumber (spheres, clst[*nneighb].sphlist,clno);
00207 clst[clno].sphlist.splice(clst[clno].sphlist.end(), clst[*nneighb].sphlist);
00208 }
00209 }
00210 }
00211
00212 delete[] clstoverlap;
00213
00214
00215 // clean clustertable of empty clusters
00216 // and set the right clno to each sphere
00217
00218 NUMBER highest=c2-1;
00219
00220 for (clno=firstclno;clno<=c2-1;clno++){
00221 if (! clst[clno].sphlist.empty()) highest=clno;
00222 else{
00223 // find first non-empty
00224 NUMBER j;
00225 for (j=clno+1;j<=c2-1;j++) if (! clst[j].sphlist.empty()) break;
00226 if (j==c2) {clno=c2;break;} // the rest is also empty
00227
00228 else{
00229 swapcluson(clst[clno],clst[j]);
00230 set_clusternumber (spheres, clst[clno].sphlist,clno);
00231 if (!clst[j].sphlist.empty()) {
00232 errorout ("error in getting rid of empty clusters (after swap)");
00233 cout <<"should be empty: clno j="<<j<<endl;
00234 cout <<"was swapped: clno ="<<clno<<endl;
00235 exit (1);
00236 }
00237 highest=clno;
00238 }
00239 }
00240 }
00241 c2=highest+1; // first empty clusternumber
00242
00243 return (c2);
00244 }
00245
00246
00247
00248 #else // if TIMEMEASUREMENT wanted
00249
00250
00251
00252
00253 // global stop-clocks for TIMEMEASUREMENT
00254 class t_combine{
00255 public:
00256 clock_t checks;
00257 clock_t findedges;
00258 clock_t createobjects;
00259 clock_t edgesphere_ovrlap;
00260 clock_t ovrlpclst_relabel;
00261 clock_t clustersplice;
00262 clock_t cleanclustertable;
00263 t_combine(){checks=0; findedges=0; createobjects=0;
00264 edgesphere_ovrlap=0;ovrlpclst_relabel=0;
00265 clustersplice=0; cleanclustertable=0;};
00266 clock_t sum(){return checks+findedges+createobjects+edgesphere_ovrlap+
00267 ovrlpclst_relabel+clustersplice+cleanclustertable;};
00268 };
00269 t_combine tcom;
00270
00271 NUMBER combine (sphere* spheres,
00272 NUMLIST &sph1,
00273 NUMLIST &sph2,
00274 cluson* clst,
00275 NUMBER firstclno,
00276 NUMBER c1,
00277 NUMBER c2,
00278 COORDFLOAT border,
00279 int direction)
00280 {
00281 clock_t ttemp;
00282 ttemp=clock();
00283
00284 if (sph1.empty() || sph2.empty()) return c2; // if one of the parts is empty,
00285 // do nothing
00286
00287 sphere s1=spheres[*sph1.begin()]; sphere s2=spheres[*sph2.begin()];
00288 sphere temp=s1; // for a faster overlap-check
00289
00290 if (s1.clno > s2.clno ) {
00291 errorout("combine is not possible, \nbecause clusno's in first list greater than in second list");
00292 cout <<s1.clno <<"\t"<< s2.clno<<endl;
00293 exit (1);
00294 }
00295 if (lower(s2.c-s1.c, 0,direction)) {
00296 errorout("combine is not possible, \nbecause coordinates of first list are higher than of second list");
00297 cout <<"(in direction "<<direction<<" ): s1.c="<<s1.c<<" s2.c="<<s2.c<<endl;
00298 exit(1);
00299 }
00300
00301 REAL radius1=find_biggestradius(spheres,sph1);
00302 REAL radius2=find_biggestradius(spheres,sph2);
00303 REAL radius=( (radius1 > radius2) ? radius1 : radius2);
00304
00305 tcom.checks+=clock()-ttemp;
00306 ttemp=clock();
00307
00308 NUMLIST e1, e2;
00309 edgespheres(spheres, sph1, e1,border, 2*radius, direction);
00310 edgespheres(spheres, sph2, e2,border, 2*radius, -direction);
00311 // cout <<"number of edgespheres: "<< e1.size()<<"\t"<<e2.size()<<endl;
00312
00313 // now find all overlapping edgespheres -> clstoverlap[ ]
00314 // = list of directly overlapping clno's (both ways "forwards" and "backwards"!)
00315 // !!! one clno might be mentioned several times !!!
00316
00317 tcom.findedges+=clock()-ttemp;
00318 ttemp=clock();
00319
00320 clusterneighbours *clstoverlap = new clusterneighbours[c2-firstclno];
00321 // for the clusternumbers [firstclno;...c1.............;c2-1]
00322 // use the (smaller) array [0;...........c1-firstclno...;c2-firstclno-1]
00323
00324 NUMLIST::iterator esph1, esph2;
00325
00326 #ifndef BOXING_ON
00327
00328 for (esph1=e1.begin();esph1!=e1.end();esph1++){
00329 for (esph2=e2.begin();esph2!=e2.end();esph2++){
00330 s1= spheres[*esph1];
00331 s2= spheres[*esph2];
00332 if (overlap2 (s1,s2,temp)) { // overlapping spheres -> overlapping clusters
00333 clstoverlap[s1.clno-firstclno].clusno.push_back ( s2.clno); // forwards
00334 clstoverlap[s2.clno-firstclno].clusno.push_back ( s1.clno); // backwards
00335 }
00336 }
00337 }
00338
00339 e1.clear(); e2.clear();
00340
00341 #else // if BOXING_ON
00342
00343 // this has changed after boxing (START)
00344 int dim=getdim(spheres[e1.front()]);
00345
00346 NUMBER nslots=GRIDSIZE/(4*radius); // heuristic: shouldn't probably be > .../(2*radius)
00347 nslots=nslots<1?1:nslots;
00348 NUMBER maxslots=(NUMBER) grid::maximal_slots_that_can_be_numbered(dim);
00349 //cout<<" |nsl:"<<nslots<<", maxsl:"<<maxslots<<"| ";
00350 nslots=nslots>maxslots?maxslots:nslots;
00351 NUMBER nboxes=pow(nslots,dim);
00352
00353 std::vector<NUMLIST > leftbox, rightbox;
00354 leftbox.resize(nboxes); rightbox.resize(nboxes);
00355
00356 tcom.createobjects+=clock()-ttemp;
00357 ttemp=clock();
00358
00359 NUMBER sn1=grid::spheres_into_boxes(spheres, e1, 0, GRIDSIZE, nslots, 2*radius, leftbox);
00360 NUMBER sn2=grid::spheres_into_boxes(spheres, e2, 0, GRIDSIZE, nslots, 2*radius, rightbox);
00361 if (sn1<0 || sn2<0) {
00362 cerr<<"problems with boxing the spheres into the grid\n";
00363 cerr<<"sn1="<<sn1<<" sn2="<<sn2<<" please write to cpp__at__AndreasKrueger__dot__de\n"<<endl;
00364 waitanykey();
00365 exit(1);
00366 }
00367 //cout<<" |sn1="<<sn1<<" sn2="<<sn2<<"| "<<flush;
00368
00369 e1.clear(); e2.clear();
00370
00371 NUMBER bn;
00372 for (bn=0;bn<nboxes;bn++){
00373
00374 for (esph1=leftbox[bn].begin();esph1!=leftbox[bn].end();esph1++){
00375 for (esph2=rightbox[bn].begin();esph2!=rightbox[bn].end();esph2++){
00376 s1= spheres[*esph1];
00377 s2= spheres[*esph2];
00378 if (overlap2 (s1,s2,temp)) { // overlapping spheres -> overlapping clusters
00379 clstoverlap[s1.clno-firstclno].clusno.push_back ( s2.clno); // forwards
00380 clstoverlap[s2.clno-firstclno].clusno.push_back ( s1.clno); // backwards
00381 }
00382 }
00383 }
00384 }
00385
00386 // this has changed after boxing (END)
00387
00388 #endif // if BOXING_ON
00389
00390
00391 tcom.edgesphere_ovrlap+=clock()-ttemp;
00392
00393
00394 // now follow all cluster-connections in clstoverlap[]
00395 // delete the multiple entries
00396 // and collect all clno's of overlapping clusters
00397 // in the list with the smallest clusternumber
00398
00399 ttemp=clock();
00400 NUMLIST::iterator nneighb ;
00401 NUMBER clno;
00402 for (clno=firstclno;clno<=c1-1;clno++){
00403 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno; // alias (!)
00404
00405 if (! neighbours.empty() ){ // no overlap
00406 if ( *(neighbours.begin()) != -1) // already counted
00407 {
00408 for (nneighb=neighbours.begin();nneighb!=neighbours.end();nneighb++){
00409
00410 while ( (*nneighb==clno) || (*(clstoverlap[*nneighb-firstclno].clusno.begin())==-1 )) {
00411 nneighb=neighbours.erase(nneighb);
00412 if (nneighb==neighbours.end()) break;
00413 }
00414 if (nneighb==neighbours.end()) break;
00415
00416 neighbours.splice(neighbours.end(), clstoverlap[*nneighb-firstclno].clusno);
00417 clstoverlap[*nneighb-firstclno].clusno.push_front(-1);
00418 }
00419 }
00420 else neighbours.clear();
00421 }
00422 }
00423
00424 tcom.ovrlpclst_relabel+=clock()-ttemp;
00425 ttemp=clock();
00426
00427 // now combine the clusters of spheres according to clstoverlap
00428 // and set the right clusternumber to each sphere
00429
00430 // idea: perhaps renumbering doesn't have to be done here, but below
00431
00432 for (clno=firstclno;clno<=c1-1;clno++){
00433 NUMLIST &neighbours=clstoverlap[clno-firstclno].clusno;
00434 if (! neighbours.empty() ){
00435 for (nneighb=neighbours.begin();nneighb!=neighbours.end();nneighb++){
00436 set_clusternumber (spheres, clst[*nneighb].sphlist,clno);
00437 clst[clno].sphlist.splice(clst[clno].sphlist.end(), clst[*nneighb].sphlist);
00438 }
00439 }
00440 }
00441
00442 tcom.clustersplice+=clock()-ttemp;
00443
00444 ttemp=clock();
00445 delete[] clstoverlap;
00446 tcom.createobjects+=clock()-ttemp;
00447
00448 // clean clustertable of empty clusters
00449 // and set the right clno to each sphere
00450
00451 ttemp=clock();
00452
00453 NUMBER highest=c2-1;
00454
00455 for (clno=firstclno;clno<=c2-1;clno++){
00456 if (! clst[clno].sphlist.empty()) highest=clno;
00457 else{
00458 // find first non-empty
00459 NUMBER j;
00460 for (j=clno+1;j<=c2-1;j++) if (! clst[j].sphlist.empty()) break;
00461 if (j==c2) {clno=c2;break;} // the rest is also empty
00462
00463 else{
00464 swapcluson(clst[clno],clst[j]);
00465 set_clusternumber (spheres, clst[clno].sphlist,clno);
00466 if (!clst[j].sphlist.empty()) {
00467 errorout ("error in getting rid of empty clusters (after swap)");
00468 cout <<"should be empty: clno j="<<j<<endl;
00469 cout <<"was swapped: clno ="<<clno<<endl;
00470 exit (1);
00471 }
00472 highest=clno;
00473 }
00474 }
00475 }
00476
00477 tcom.cleanclustertable+=clock()-ttemp;
00478 ttemp=clock();
00479
00480 c2=highest+1; // first empty clusternumber
00481 return (c2);
00482 }
00483
00484 #endif // TIMEMEASUREMENT
00485
00486
00487 // the first tests: divide once, twice, three times...
00488
00489 // only dividing once, without analysis of histogram
00490 // used for divide_two_times & three_times
00491 NUMBER divide_once_count_and_combine (sphere* spheres,
00492 sphere* temp,
00493 NUMLIST &all,
00494 cluson *clst,
00495 NUMBER firstclno,
00496 COORDFLOAT border,
00497 int direction)
00498 {
00499
00500 NUMBER c1,c2;
00501 NUMLIST sph1,sph2;
00502 divide_all_spheres(spheres, all, sph1,sph2,border,direction);
00503
00504 c1=list_to_given_array_and_find_cluster(spheres, temp, sph1, firstclno);
00505 c2=list_to_given_array_and_find_cluster(spheres, temp, sph2, c1);
00506 cout <<"found "<<c1-firstclno<<" and "<<c2-c1<<" clusters ";
00507 cout <<"(c1="<<c1<<", c2="<<c2<<")";
00508
00509 make_clusterlist_array (spheres, all, clst);
00510
00511 c2=combine(spheres,sph1,sph2,clst,firstclno,c1,c2,border,direction);
00512 return (c2);
00513 }
00514
00515 // divide once, count and analyze histogram
00516 NUMBER divide_once_count_and_analyze(sphere* spheres,
00517 sphere* temp,
00518 NUMLIST &all,
00519 NUMBER N,
00520 NUMBER *histogram,
00521 NUMBER &biggestcl,
00522 REAL &averagecl)
00523 {
00524
00525 cluson *clusters=new cluson[N+1];
00526 NUMBER nextclno;
00527
00528 nextclno=divide_once_count_and_combine (spheres, temp, all, clusters,1,GRIDSIZE/2, 1);
00529 cout <<" -> combined to: "<<nextclno-1<<" clusters"<<endl;
00530
00531 NUMBER biggestcl_no=analyze_clusters (clusters, 1, nextclno-1,histogram,N,biggestcl,averagecl);
00532
00533 delete[] clusters;
00534
00535 return (biggestcl_no);
00536 }
00537
00538
00539 NUMBER divide_two_times_count_and_analyze(sphere* spheres,
00540 sphere* temp,
00541 NUMLIST &all,
00542 NUMBER N,
00543 NUMBER *histogram,
00544 NUMBER &biggestcl,
00545 REAL &averagecl)
00546 {
00547
00548 cluson *clusters=new cluson[N+1];
00549
00550 NUMBER nextclno1, nextclno2;
00551 NUMBER cl2;
00552 NUMLIST sph1,sph2;
00553
00554 COORDFLOAT border1=GRIDSIZE/2;
00555
00556 divide_all_spheres(spheres, all,sph1,sph2, border1,1);
00557
00558 //cout <<"A) count clusters in left part..."<<endl;
00559 nextclno1=divide_once_count_and_combine (spheres, temp, sph1, clusters,1,border1, 2);
00560 cout <<" -> combined to: "<<nextclno1-1<<" clusters"<<endl;
00561
00562 //cout <<"B) count clusters in right part..."<<endl;
00563 nextclno2=divide_once_count_and_combine (spheres, temp, sph2, clusters,nextclno1,border1, 2);
00564 cout <<" -> combined to: "<<nextclno2-nextclno1<<" clusters"<<endl;
00565
00566 cout <<"Now combine the two results...";
00567 cl2=combine (spheres,sph1,sph2,clusters,1,nextclno1,nextclno2,border1,1);
00568 cout <<" -> combined to: "<<cl2-1<<" clusters"<<endl;
00569
00570 NUMBER biggestcl_no=analyze_clusters (clusters, 1, cl2-1,histogram,N,biggestcl,averagecl);
00571
00572 delete[] clusters;
00573
00574 return (biggestcl_no);
00575 }
00576
00577 NUMBER divide_three_times_count_and_analyze(sphere* spheres,
00578 sphere* temp,
00579 NUMLIST &all,
00580 NUMBER N,
00581 NUMBER *histogram,
00582 NUMBER &biggestcl,
00583 REAL &averagecl)
00584 {
00585 cluson *clusters=new cluson[N+1];
00586
00587 NUMBER nextclno1, nextclno2;
00588 NUMBER cl1, cl2;
00589 NUMLIST sph1,sph2,sph3,sph4,sph5,sph6;
00590
00591 COORDFLOAT border1=GRIDSIZE/2;
00592 COORDFLOAT border2=GRIDSIZE/4;
00593
00594 divide_all_spheres(spheres, all, sph1,sph2,border1,1);
00595
00596 divide_all_spheres(spheres, sph1,sph3,sph4,border1,2);
00597 nextclno1=divide_once_count_and_combine (spheres, temp, sph3, clusters,1,border1-border2, 1); cout<<endl;
00598 nextclno2=divide_once_count_and_combine (spheres, temp, sph4, clusters,nextclno1,border1-border2, 1); cout<<endl;
00599 cl1=combine (spheres,sph3,sph4,clusters,1,nextclno1,nextclno2,border1,2);
00600 cout <<" -> combined to: "<<cl1-1<<" clusters"<<endl;
00601
00602 divide_all_spheres(spheres, sph2,sph5,sph6,border1,2);
00603 nextclno1=divide_once_count_and_combine (spheres, temp, sph5, clusters,cl1,border1+border2, 1); cout <<endl;
00604 nextclno2=divide_once_count_and_combine (spheres, temp, sph6, clusters,nextclno1,border1+border2, 1); cout<<endl;
00605 cl2=combine (spheres,sph5,sph6,clusters,cl1,nextclno1,nextclno2,border1,2);
00606 cout <<" -> combined to: "<<cl2-cl1<<" clusters"<<endl;
00607
00608 cout <<"Now these two combined to: ";
00609 cl2=combine(spheres,sph1,sph2,clusters,1,cl1,cl2,border1,1);
00610 cout<<cl2-1<<" clusters"<<endl;
00611
00612 NUMBER biggestcl_no=analyze_clusters (clusters, 1, cl2-1,histogram,N,biggestcl,averagecl);
00613
00614 delete[] clusters;
00615
00616 return (biggestcl_no);
00617 }
00618
00619
00620
00621 // the divide d-times algo:
00622 // this algo divides the space into smaller square cells
00623 //
00624 // with the same number of cuts in each dimensional direction
00625 // -> cuts = dividings * dimension
00626 // #cells = 2^cuts = 2^( dividings * dimension)
00627 // advantage: easy to code
00628 // disadvantage: very big steps
00629 // 4dim: 2 dividings -> 2^(4*2) = 256 boxes
00630 // 3 dividings -> 2^(4*3) = 4096 boxes
00631 //
00632 //
00633 // !!! NOT USED ANYMORE - only for testing !!!
00634
00635 COUNTER maximal_dividings_in_one_coordinate(COORDFLOAT spacelength, REAL radius) {
00636 REAL boxlength=2*radius;
00637 REAL maximal_boxes_in_one_direction=spacelength/boxlength;
00638 REAL maximal_dividings=log(maximal_boxes_in_one_direction) / log (2);
00639 return ( (int) maximal_dividings );
00640 }
00641
00642 NUMBER partcount (sphere* spheres,
00643 NUMLIST &all,
00644 COORDFLOAT l,
00645 COORDFLOAT r,
00646 COORDFLOAT lm,
00647 COORDFLOAT rm,
00648 COUNTER divstep,
00649 COUNTER divmax,
00650 int direction,
00651 int dimension,
00652 cluson *clusters,
00653 NUMBER firstclno)
00654 {
00655 NUMBER nextclno;
00656
00657 if (divmax==0)
00658 {
00659 nextclno=list_to_temp_array_and_find_cluster(spheres, all, firstclno);
00660 make_clusterlist_array (spheres, all, clusters);
00661 return nextclno;
00662 }
00663
00664 if ( direction==dimension ) {
00665 if ( divstep==divmax ) {
00666 NUMLIST s1, s2;
00667 NUMBER ncl1, ncl2;
00668 divide_all_spheres (spheres, all, s1,s2,(l+r)/2,direction);
00669 //cout <<"\ndir="<<direction<<"\tborders l,r=\t"<<l<<"\t"<<r;
00670 //cout <<"\t# of spheres="<<s1.size()<<" and "<<s2.size()<<endl;
00671 ncl1=list_to_temp_array_and_find_cluster(spheres, s1, firstclno);
00672 ncl2=list_to_temp_array_and_find_cluster(spheres, s2, ncl1);
00673 make_clusterlist_array (spheres, all, clusters);
00674 nextclno=combine(spheres, s1,s2,clusters,firstclno,ncl1,ncl2,(l+r)/2,direction);
00675 }
00676 else {
00677 NUMLIST s1, s2;
00678 NUMBER ncl1, ncl2;
00679 divide_all_spheres (spheres, all, s1,s2,(l+r)/2,direction);
00680 //cout <<"diV ";
00681 ncl1=partcount(spheres, s1, l,(l+r)/2, lm,rm, divstep+1, divmax, direction,dimension, clusters, firstclno);
00682 ncl2=partcount(spheres, s2, (l+r)/2,r, lm,rm, divstep+1, divmax, direction,dimension, clusters, ncl1);
00683 nextclno=combine(spheres, s1,s2,clusters,firstclno,ncl1,ncl2,(l+r)/2,direction);
00684 }
00685 }
00686
00687 else {
00688 if (divstep <= divmax) {
00689 NUMLIST s1, s2;
00690 NUMBER ncl1, ncl2;
00691 divide_all_spheres (spheres, all, s1,s2,(l+r)/2,direction);
00692 //cout <<"diV ";
00693 ncl1=partcount(spheres, s1, l,(l+r)/2, lm,rm, divstep+1, divmax, direction,dimension, clusters, firstclno);
00694 ncl2=partcount(spheres, s2, (l+r)/2,r, lm,rm, divstep+1, divmax, direction,dimension, clusters, ncl1);
00695 nextclno=combine(spheres, s1,s2,clusters,firstclno,ncl1,ncl2,(l+r)/2,direction);
00696 }
00697
00698 else {
00699 //cout <<"diM ";
00700 nextclno=partcount(spheres, all, lm,rm, lm,rm, 1, divmax, direction+1,dimension, clusters, firstclno);
00701 }
00702 }
00703
00704 return nextclno;
00705 }
00706
00707
00708 NUMBER divide_d_times_count_and_analyze(sphere* spheres,
00709 NUMLIST &all,
00710 NUMBER N,
00711 NUMBER d,
00712 NUMBER *histogram,
00713 NUMBER &biggestcl,
00714 REAL &averagecl)
00715 {
00716 cluson *clusters=new cluson[N+1];
00717
00718 int dimension=getdim( spheres[*all.begin()].c );
00719
00720 NUMBER nextcluster=partcount(spheres, all, 0, GRIDSIZE, (COORDFLOAT)0, GRIDSIZE, 1, d, 1,dimension, clusters, 1);
00721
00722 NUMBER biggestcl_no=analyze_clusters (clusters, 1, nextcluster-1,histogram,N,biggestcl,averagecl);
00723
00724 delete[] clusters;
00725
00726 return (biggestcl_no);
00727 }
00728
00729
00730
00731
00732 COUNTER choose_optimal_cuts (int dim, NUMBER N, REAL R){
00733
00734 REAL rc=give_radius(ff_critical_guessed(N, dim), N, GRIDSIZE, dim);
00735 REAL rs=give_radius(saturating_fillingfactor(dim, N), N, GRIDSIZE, dim);
00736 // cout <<"rc="<<rc<<" rs="<<rs<<" this R="<<R<<endl;
00737 REAL lowest_N=39.0625;
00738 NUMBER Nclose=rounded(lowest_N * rounded(pow(2, 1+rounded(log(N/(2*lowest_N))/log(2)))));
00739 if (Nclose>480000) Nclose=480000;
00740 // cout <<"N="<<N<<" Nclose="<<Nclose<<endl;
00741
00742 COUNTER cuts_c=choose_optimal_cuts_ffc(dim, Nclose);
00743 COUNTER cuts_s=choose_optimal_cuts_ffs(dim, Nclose);
00744 if (cuts_s > cuts_c) cuts_s=cuts_c;
00745
00746 if (dim==1) return cuts_c; // because then rs==rc
00747
00748 REAL cuts=cuts_c + (cuts_s-cuts_c)* (R-rc)/(rs-rc);
00749
00750 // cout <<"cuts_c="<<cuts_c<<" cuts_s="<<cuts_s<<" cuts_R="<<cuts<<" return "<<endl;
00751 cuts=rounded(cuts);
00752 return (COUNTER)(cuts>0?cuts:0);
00753 }
00754
00755 void test_cuts_function(){
00756 cout<<"type in dim, N"<<endl;
00757 int dim; COUNTER N;
00758 cin>>dim>>N;
00759 REAL R=0;
00760 while (R!=-1){
00761 cout <<"\ntype in R (-1 for end)";
00762 cin>>R;
00763 cout <<choose_optimal_cuts(dim, N, R);
00764 }
00765 }
00766
00767
00768
00769
00770
00771 // the divide-by-c-cuts counters:
00772 // this algo divides the space into smaller cells,
00773 // if necessary non-square, rectangular
00774 //
00775 // with a total number of cuts in all dimensional directions
00776 // -> cuts = dividings.dim1 + dividings.dim2 + ...+ dividings.dimd
00777 // #cells = 2^cuts
00778 // advantage: smaller steps
00779 // 4dim: 8 cuts -> 2^8 = 256 boxes (cuts=2+2+2+2)
00780 // 9 cuts -> 2^9 = 512 boxes (cuts=3+2+2+2)
00781
00782 // !!! this one is used NOW !!!
00783
00784 #ifndef TIMEMEASUREMENT
00785
00786 NUMBER dividecount (sphere* spheres, NUMLIST &all,
00787 COORDFLOAT l, COORDFLOAT r,
00788 COORDFLOAT lm,COORDFLOAT rm,
00789 COUNTER divstep, COUNTER divcounter,
00790 int direction, int dimension,
00791 cluson *clusters, NUMBER firstclno)
00792 {
00793 NUMBER nextclno;
00794
00795 if (divstep < divcounter/direction ) {
00796 NUMLIST s1, s2;
00797 NUMBER ncl1, ncl2;
00798 divide_all_spheres (spheres, all, s1,s2,(l+r)/2,direction);
00799 ncl1=dividecount(spheres, s1, // RECURSION!
00800 l,(l+r)/2, lm,rm,
00801 divstep+1, divcounter,
00802 direction,dimension,
00803 clusters, firstclno);
00804 ncl2=dividecount(spheres, s2, // RECURSION!
00805 (l+r)/2,r, lm,rm,
00806 divstep+1, divcounter,
00807 direction,dimension,
00808 clusters, ncl1);
00809 nextclno=combine(spheres, s1,s2,clusters,firstclno,ncl1,ncl2,(l+r)/2,direction);
00810 }
00811
00812 else {
00813 if (direction==1) {
00814 // the core clusterCounter
00815 // (naive recursion or iteration with full network)
00816 nextclno=list_to_temp_array_and_find_cluster(spheres, all, firstclno);
00817 make_clusterlist_array (spheres, all, clusters);
00818 }
00819 else {
00820 // continue in next-lower direction (divcounter-divstep cuts left to do)
00821 nextclno=dividecount(spheres, all,
00822 lm,rm, lm,rm,
00823 0, divcounter-divstep,
00824 direction-1,dimension,
00825 clusters, firstclno);
00826 }
00827 }
00828 return nextclno;
00829 }
00830
00831 #else // if defined TIMEMEASUREMENT
00832
00833 class timedividecount{
00834 public:
00835 clock_t cellcount;
00836 clock_t clusterlist;
00837 clock_t divide;
00838 clock_t combine;
00839 clock_t dividerecursion;
00840 timedividecount(){cellcount=0; clusterlist=0; divide=0; combine=0;dividerecursion=0;};
00841 };
00842
00843 timedividecount tdiv;
00844
00845
00846 NUMBER dividecount (sphere* spheres, NUMLIST &all,
00847 COORDFLOAT l, COORDFLOAT r,
00848 COORDFLOAT lm,COORDFLOAT rm,
00849 COUNTER divstep, COUNTER divcounter,
00850 int direction, int dimension,
00851 cluson *clusters, NUMBER firstclno)
00852 {
00853 clock_t ttemp;
00854 ttemp=clock();
00855
00856 NUMBER nextclno;
00857 if (divstep < divcounter/direction ) {
00858 NUMLIST s1, s2;
00859 NUMBER ncl1, ncl2;
00860
00861 ttemp=clock();
00862 divide_all_spheres (spheres, all, s1,s2,(l+r)/2,direction);
00863 tdiv.divide+=clock()-ttemp;
00864
00865 ncl1=dividecount(spheres, s1, // RECURSION!
00866 l,(l+r)/2, lm,rm,
00867 divstep+1, divcounter,
00868 direction,dimension,
00869 clusters, firstclno);
00870 ncl2=dividecount(spheres, s2, // RECURSION!
00871 (l+r)/2,r, lm,rm,
00872 divstep+1, divcounter,
00873 direction,dimension,
00874 clusters, ncl1);
00875
00876 ttemp=clock();
00877 nextclno=combine(spheres, s1,s2,clusters,firstclno,ncl1,ncl2,(l+r)/2,direction);
00878 tdiv.combine+=clock()-ttemp;
00879 }
00880
00881 else {
00882 if (direction==1) {
00883 // the core clusterCounter
00884 // (naive recursion or iteration with full network)
00885
00886 ttemp=clock();
00887 nextclno=list_to_temp_array_and_find_cluster(spheres, all, firstclno);
00888 tdiv.cellcount+=clock()-ttemp;
00889 ttemp=clock();
00890 make_clusterlist_array (spheres, all, clusters);
00891 tdiv.clusterlist+=clock()-ttemp;
00892 }
00893 else {
00894 // continue in next-lower direction (divcounter-divstep cuts left to do)
00895 nextclno=dividecount(spheres, all,
00896 lm,rm, lm,rm,
00897 0, divcounter-divstep,
00898 direction-1,dimension,
00899 clusters, firstclno);
00900 }
00901 }
00902 return nextclno;
00903 }
00904
00905 #endif // defined TIMEMEASUREMENT
00906
00907
00908 // This intermediate step returns an array of clusons=spherelist
00909 // used for visualization
00910 NUMBER count_analyze_and_return_clusters(NUMBER c,
00911 sphere* spheres, NUMLIST &all,
00912 cluson *clusters,
00913 NUMBER *histogram,
00914 NUMBER &biggestcl, REAL &averagecl)
00915 {
00916 // if (c>=1) {
00917 int dimension=getdim( spheres[*all.begin()].c );
00918
00919 NUMBER nextcluster=dividecount(spheres, all,
00920 0,GRIDSIZE,0,GRIDSIZE,
00921 0, c,
00922 dimension,dimension,
00923 clusters, 1);
00924
00925 erasehisto(histogram,1,all.size());
00926 NUMBER biggestcl_no;
00927 NUMBER totalnumber= makehistogram ( clusters,
00928 1,nextcluster-1,
00929 histogram,
00930 biggestcl_no,
00931 biggestcl,
00932 averagecl);
00933
00934 if (totalnumber != (NUMBER)all.size()) {
00935 errorout("not all the spheres are in clusters!");
00936 cout << "number of spheres: N="<<all.size();
00937 cout <<" counted="<<totalnumber<<endl;
00938 exit(1);
00939 }
00940 else ; // cout <<"N="<<all.size()<<" counted="<<totalnumber<<endl;
00941
00942 // }
00943 // else exit(1);
00944
00945 return (nextcluster-1); // returns the number of clusters (M0)
00946
00947 }
00948
00949
00950 // the same as count_analyze_and_return_clusters, but throws away the cluson-array
00951 NUMBER divide_by_c_cuts_count_and_analyze(sphere* spheres, NUMLIST &all,
00952 NUMBER c,
00953 NUMBER *histogram,
00954 NUMBER &biggestcl,
00955 REAL &averagecl)
00956 {
00957 cluson *clusters=new cluson[all.size()+1];
00958
00959 NUMBER number_of_clusters= count_analyze_and_return_clusters(c,
00960 spheres, all,
00961 clusters, histogram,
00962 biggestcl, averagecl);
00963
00964 delete[] clusters;
00965
00966 return number_of_clusters;
00967 }
00968
00969
00970 // this one is used NOW (in find_ffc, into_file8, etc.)
00971 // it returns an array of clusons=spherelists
00972
00973 #ifndef TIMEMEASUREMENT
00974
00975 NUMBER divide_count_and_analyze(sphere* spheres, NUMLIST &all,
00976 NUMBER c,
00977 cluson *clusters,
00978 NUMBER &biggestcl,
00979 REAL &mean_clsz)
00980 {
00981 NUMBER *histogram=new NUMBER[all.size()+1];
00982 if (histogram==NULL) exit_out_of_memory("divide_count_and_analyze(...): NUMBER *histogram");
00983 erasehisto(histogram,1,all.size());
00984
00985 int dimension=getdim( spheres[*all.begin()].c );
00986
00987 NUMBER numberof_cl=dividecount(spheres, all,
00988 0,GRIDSIZE,(COORDFLOAT)0,GRIDSIZE,
00989 0, c,
00990 dimension,dimension,
00991 clusters, 1);
00992
00993 NUMBER biggestcl_no;
00994 NUMBER totalnumber= makehistogram ( clusters,
00995 1,numberof_cl-1,
00996 histogram,
00997 biggestcl_no,
00998 biggestcl,
00999 mean_clsz);
01000 if (totalnumber != (NUMBER)all.size()) {
01001 errorout("not all the spheres are in clusters!");
01002 cout << "number of spheres: N="<<all.size();
01003 cout <<" counted="<<totalnumber<<endl;
01004 waitanykey();
01005 }
01006
01007 delete[] histogram;
01008
01009 return numberof_cl-1;
01010 }
01011
01012 #else // if TIMEMEASUREMENT wanted
01013
01014 NUMBER divide_count_and_analyze(sphere* spheres, NUMLIST &all,
01015 NUMBER c,
01016 cluson *clusters,
01017 NUMBER &biggestcl,
01018 REAL &mean_clsz)
01019 {
01020 NUMBER *histogram=new NUMBER[all.size()+1];
01021 if (histogram==NULL) exit_out_of_memory("divide_count_and_analyze(...): NUMBER *histogram");
01022 erasehisto(histogram,1,all.size());
01023
01024 int dimension=getdim( spheres[*all.begin()].c );
01025
01026 clock_t ttemp=clock();
01027 NUMBER numberof_cl=dividecount(spheres, all,
01028 0,GRIDSIZE,(COORDFLOAT)0,GRIDSIZE,
01029 0, c,
01030 dimension,dimension,
01031 clusters, 1);
01032 tdiv.dividerecursion+=clock()-ttemp;
01033
01034 NUMBER biggestcl_no;
01035 NUMBER totalnumber= makehistogram ( clusters,
01036 1,numberof_cl-1,
01037 histogram,
01038 biggestcl_no,
01039 biggestcl,
01040 mean_clsz);
01041 if (totalnumber != (NUMBER)all.size()) {
01042 errorout("not all the spheres are in clusters!");
01043 cout << "number of spheres: N="<<all.size();
01044 cout <<" counted="<<totalnumber<<endl;
01045 waitanykey();
01046 }
01047
01048 delete[] histogram;
01049
01050 return numberof_cl-1;
01051 }
01052
01053 #endif // TIMEMEASUREMENT
01054
01055
01056 } // end of namespace counters
01057
01058
| Diploma Thesis Sourcecode
Documentation check out the text and the executable binaries |