Some improvements

NobuhiroMoteki · web-flow · commit 92ab42862314 · 2019-12-09T13:59:26.000+09:00
Change the output result: output the aggregate files of all levels at once.
Change the output result: the user can specify the number of aggregates produced by a single execution (using different random seeds).
Improve reliability: Internal tests algorithm for non-overlapping and attached conditions are improved.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -1,7 +1,7 @@
 cmake_minimum_required (VERSION 3.10.0)
 set(CMAKE_CXX_COMPILER "g++")
 set(CMAKE_CXX_FLAGS "-O3 -fopenmp -DNDEBUG -DEIGEN_NO_DEBUG")
-include_directories(/mnt/c/Users/moteki/WSL_Share_Ubuntu1804/C++/eigen_3_3_5)
+include_directories(~/C++/eigen_3_3_5)
 
 project(aggregate_gen)
 
diff --git a/README.md b/README.md
@@ -2,6 +2,11 @@
 ---
 
 ### History:
+#### Oct 23, 2019:
+  * Change the output result: output the aggregate files of all levels at once.
+  * Change the output result: the user can specify the number of aggreagtes produced by a single execution (using different random seeds).
+  * Improve reliability: Internal tests algorithm for non-overlapping and attached conditions are improved.
+
 #### Aug 14, 2018:
   * Bugfix: An algorithm error causing the partial overlapping between some spherules has been fixed.
   * Improve reliability: Internal tests of non-overlapping and attached conditions are added.
@@ -43,12 +48,13 @@
     levels : Number of hierarchical levels in cluster-cluster aggregation (10).
     kf : Fractal prefactor (1.0).
     Df : Fractal dimension (2.5).
+    num_agg : Number of aggregates generated using different random seeds.
 
     If you skip the [parameter_list], the code uses the default values in the parenthesis. Total number of spherules is *N* = num_sph_SA*2<sup>levels</sup>.
 ```
 
 ### Output data
-  The result is written in a text file "agg_N???_kf???_Df???.out", where the three strings ??? indicate the total number of spherules *N*, the fractal prefactor *kf*, and the fractal dimension *Df*.
+  The result is written in a text file "agg???_N???_kf???_Df???.out", where the three strings ??? indicate the aggregate id number (0 ~ num_agg), total number of spherules *N*, the fractal prefactor *kf*, and the fractal dimension *Df*, respectively.
 
   In this file, each row shows the cartesian coordinate "x y z" of each spherules relative to the centroid of cluster, where the
   length is scaled such that spherule radius is 1.0.
diff --git a/aggregate_gen_CCA.cpp b/aggregate_gen_CCA.cpp
@@ -12,34 +12,39 @@ using namespace std;
 
 
 //return pos_sph matrix
-MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& levels, const float& kf, const float& Df, const float& tol){
+vector<MatrixXd> aggregate_gen_CCA(const double& a, const int& num_sph_SA, const int& levels, const double& kf, const double& Df, const double& tol){
+
+  vector<MatrixXd> pos_sph_agg_output; // vector of aggregate geometry from level 0 to "levels".
   int iter1max {1000};
   int iter2max {10000};
 
   int num_SA_agg= pow(2,levels);
-  int totnum_sph= num_SA_agg*num_sph_SA;
-  cout << "totnum_sph= " << totnum_sph << endl;
+  int final_totnum_sph= num_SA_agg*num_sph_SA;
+  cout << "final totnum_sph= " << final_totnum_sph << endl;
 
-  vector<MatrixXf> pos_sph_SAagg;
+  vector<MatrixXd> pos_sph_SAagg;
   for(int i= 0; i<num_SA_agg; ++i ) pos_sph_SAagg.push_back(aggregate_gen_SA(a,num_sph_SA,kf,Df,tol));
 
-  vector<MatrixXf> pos_sph_agg_old;
+  pos_sph_agg_output.push_back(pos_sph_SAagg[0]); // level 0 aggregate
+
+
+  vector<MatrixXd> pos_sph_agg_old;
 
   // start time
   chrono::time_point<chrono::steady_clock> CCA_start_time= chrono::steady_clock::now();
 
   for(int k= 1; k <= levels; ++k){
       cout << "Generating aggregates in " << k << "th level ... please wait ..." << endl;
-      vector<MatrixXf> pos_sph_agg;
+      vector<MatrixXd> pos_sph_agg;
 
       for(int L= 1; L <= pow(2,levels-k); ++L){
           //cout << "generating " << L << " th " << "aggregate in level" << k << endl;
           int N1= num_sph_SA*pow(2,k-1);
           int N2= N1;
 
           // aggregation process in k-th level
-          MatrixXf pos_sph_agg1(3,N1);
-          MatrixXf pos_sph_agg2(3,N2);
+          MatrixXd pos_sph_agg1(3,N1);
+          MatrixXd pos_sph_agg2(3,N2);
 
           if(k==1){
             pos_sph_agg1= pos_sph_SAagg[2*L-2];
@@ -50,27 +55,27 @@ MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& lev
           }
 
           // random rotation of aggregate 1 around the centroid (origin)
-          float theta,phi; // polar angles
+          double theta,phi; // polar angles
           theta= acos(2.0f*urf(re)-1.0f);
           phi= 2*M_PI*urf(re);
 
-          Vector3f nvec; // unit vector of rotation axis
+          Vector3d nvec; // unit vector of rotation axis
           nvec << sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta); //axis of rotation randomly choosen over unit sphere
 
-          float psi; // angle of rotation
+          double psi; // angle of rotation
           psi= 2*M_PI*urf(re); // randomly choosen rotation angle
 
-          Matrix3f rot_mat; // 3D rotation matrix
-          rot_mat= AngleAxisf(psi,nvec);
+          Matrix3d rot_mat; // 3D rotation matrix
+          rot_mat= AngleAxisd(psi,nvec);
 
           for(int i= 0; i < pos_sph_agg1.cols(); ++i) pos_sph_agg1.col(i)=rot_mat*pos_sph_agg1.col(i); // random rotation of agg1
 
-          float R1_2; // square of gyration radius of agg1
-          float R2_2; // square of gyration radius of agg2
+          double R1_2; // square of gyration radius of agg1
+          double R2_2; // square of gyration radius of agg2
           R1_2= a*a + pos_sph_agg1.colwise().squaredNorm().mean();
           R2_2= a*a + pos_sph_agg2.colwise().squaredNorm().mean();
 
-          float rN;
+          double rN;
           rN= a*a*(N1+N2)*(N1+N2)/(N1*N2)*pow((N1+N2)/kf,2/Df)-(N1+N2)/N2*R1_2-(N1+N2)/N1*R2_2;
           rN= sqrt(rN);
 
@@ -81,26 +86,22 @@ MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& lev
 
               pos_sph_agg1.colwise() -= pos_sph_agg1.rowwise().mean(); // translate the centroid of agg1 to origin
 
-              float phi;
+              double phi;
               phi= 2*M_PI*urf(re);
-              float u; // uniform random number [-1.0,1.0)
-              u= 2*urf(re)-1.0f;
+              double u; // uniform random number [-1.0,1.0)
+              u= 2*urf(re)-1.0;
 
-              Vector3f cen_pos_agg1; // centroid of agg1
-              cen_pos_agg1 << sqrt(1.0f-u*u)*cos(phi), sqrt(1.0f-u*u)*sin(phi), u;
+              Vector3d cen_pos_agg1; // centroid of agg1
+              cen_pos_agg1 << sqrt(1.0-u*u)*cos(phi), sqrt(1.0-u*u)*sin(phi), u;
               cen_pos_agg1 *= rN; // cen_pos_agg1 is set to a uniform random cartesian position (x,y,z) over the spherical surface of radius rN
               pos_sph_agg1.colwise() += cen_pos_agg1;
 
-
-              float theta;
-              theta= acos(2.0f*urf(re)-1.0f);
+              double theta;
+              theta= acos(2.0*urf(re)-1.0);
               phi= 2*M_PI*urf(re);
-              Vector3f nvec; // unit vector of rotation axis
+              Vector3d nvec; // unit vector of rotation axis
               nvec << sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta); //axis of rotation randomly choosen over unit sphere
 
-
-
-
               // random rotation of aggregate 2 around the centroid (origin)
               #pragma omp parallel for
               for(int iter2= 0; iter2 < iter2max; ++iter2){
@@ -109,31 +110,31 @@ MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& lev
                   found_local= found;
                   if(found) continue;
 
-                  MatrixXf pos_sph_agg1_thread(3,N1); // thread local copy of pos_sph_agg1
-                  MatrixXf pos_sph_agg2_thread(3,N2); // thread local copy of pos_sph_agg2
+                  MatrixXd pos_sph_agg1_thread(3,N1); // thread local copy of pos_sph_agg1
+                  MatrixXd pos_sph_agg2_thread(3,N2); // thread local copy of pos_sph_agg2
                   pos_sph_agg1_thread= pos_sph_agg1;
                   pos_sph_agg2_thread= pos_sph_agg2;
 
-                  float psi; // angle of rotation
+                  double psi; // angle of rotation
 
                   #pragma omp critical
                   {
                      psi= 2*M_PI*urf(re);
                   }
 
-                  Matrix3f rot_mat; // 3D rotation matrix
-                  rot_mat= AngleAxisf(psi,nvec);
+                  Matrix3d rot_mat; // 3D rotation matrix
+                  rot_mat= AngleAxisd(psi,nvec);
                   for(int i= 0; i < N2; ++i) pos_sph_agg2_thread.col(i)=rot_mat*pos_sph_agg2_thread.col(i); // random rotation of agg2
                   //  cout << "iter2 pos_sph_agg2 " << endl << pos_sph_agg2 << endl;
 
                   bool attached {false};
                   for(int i= 0; i < N2; ++i){
-                      RowVectorXf dist_pair(N1);
+                      RowVectorXd dist_pair(N1);
                       dist_pair=(pos_sph_agg1_thread.colwise()-pos_sph_agg2_thread.col(i)).colwise().norm();
                     //  cout << "dist_pair" << endl << dist_pair << endl;
-                      float mindist= dist_pair.minCoeff();
+                      double mindist= dist_pair.minCoeff();
                       //cout << "2*a, mindist " << 2*a << ' ' << mindist << endl;
-                      bool overrapped {mindist <= (2*a-tol)};
+                      bool overrapped {mindist <= 2*a};
                       if(overrapped) {
                         //exit; // 20180813 removed
                         attached= false; // 20180813 added
@@ -146,7 +147,7 @@ MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& lev
                   // construction of new aggregate if attached condition is satisfied
                   if(attached){
                       //cout << "attached " << attached << endl;
-                      MatrixXf pos_sph_agg_tmp(3,N1+N2);
+                      MatrixXd pos_sph_agg_tmp(3,N1+N2);
                       pos_sph_agg_tmp.block(0,0,3,N1)= pos_sph_agg1_thread;
                       pos_sph_agg_tmp.block(0,N1,3,N2)= pos_sph_agg2_thread;
                       pos_sph_agg_tmp.colwise() -= pos_sph_agg_tmp.rowwise().mean(); // translate the centroid of new agg to origin
@@ -171,6 +172,7 @@ MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& lev
       }// L loop
 
       pos_sph_agg_old= move(pos_sph_agg);
+      pos_sph_agg_output.push_back(pos_sph_agg_old[0]);
 
   }// k loop
 
@@ -179,37 +181,43 @@ MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& lev
   cout << "CCA computation time: "
   << chrono::duration_cast<chrono::milliseconds>(CCA_end_time - CCA_start_time).count() << " milliseconds" << endl;
 
-  MatrixXf pos_sph_agg_final;
-  pos_sph_agg_final = pos_sph_agg_old[0];
-
-  assert(totnum_sph == pos_sph_agg_final.cols());
-
-  // final check of the attached condition
-  double buf_factor= 5.0; //buffer factor of torelance
-  bool non_overlapped;
-  ArrayXf min_diff_from2a(totnum_sph);
-  for(int i= 0; i< totnum_sph; ++i){
-      ArrayXf distance_from_this_sph(totnum_sph);
-      distance_from_this_sph= (pos_sph_agg_final.colwise()-pos_sph_agg_final.col(i)).colwise().norm().array();
-      non_overlapped= {(distance_from_this_sph-(2*a-buf_factor*tol) > 0)(i) == false};
-      for(int j= 0; j< totnum_sph; ++j){
-          if(j != i) {
-            non_overlapped= non_overlapped && ((distance_from_this_sph-(2*a-buf_factor*tol) > 0)(j) == true);
+  // geometry check for output
+  bool non_overlapped_all {true};
+  bool attached_at_least_1pair_all {true};
+
+  for(int k= 0; k <= levels; ++k){
+      MatrixXd pos_sph_agg_k_out {pos_sph_agg_output[k]};
+      int num_sph = num_sph_SA*pow(2,k);
+      assert(num_sph == pos_sph_agg_k_out.cols());
+      // final check of the attached condition
+      bool non_overlapped;
+      ArrayXd min_diff_from2a(num_sph);
+      for(int i= 0; i< num_sph; ++i){
+          ArrayXd distance_from_this_sph(num_sph);
+          distance_from_this_sph= (pos_sph_agg_k_out.colwise()-pos_sph_agg_k_out.col(i)).colwise().norm().array();
+          non_overlapped= {(distance_from_this_sph-2*a >= 0)(i) == false};
+          for(int j= 0; j< num_sph; ++j){
+              if(j != i) {
+                non_overlapped= non_overlapped && ((distance_from_this_sph-2*a > 0)(j) == true);
+              }
           }
+         //cout << i << ", " << (distance_from_this_sph-(2*a-buf_factor*tol) >= 0).transpose() << endl ;
+          min_diff_from2a(i)=((pos_sph_agg_k_out.colwise()-pos_sph_agg_k_out.col(i)).colwise().norm().array()-2*a).array().abs().minCoeff();
       }
-     //cout << i << ", " << (distance_from_this_sph-(2*a-buf_factor*tol) >= 0).transpose() << endl ;
-      min_diff_from2a(i)=((pos_sph_agg_final.colwise()-pos_sph_agg_final.col(i)).colwise().norm().array()-2*a).array().abs().minCoeff();
+      bool attached_at_least_1pair {min_diff_from2a.maxCoeff() < tol};
+      string non_overlapped_test= {non_overlapped ? "Success" : "Failed"};
+      string attached_at_least_1pair_test= {attached_at_least_1pair ? "Success" : "Failed"};
+      cout << k << "th level, non_overlapped_check         : " << non_overlapped_test << endl;
+      cout << k << "th level, attached_at_least_1pair_check: " << attached_at_least_1pair_test << endl;
+
+      non_overlapped_all = non_overlapped_all && non_overlapped;
+      attached_at_least_1pair_all = attached_at_least_1pair_all && attached_at_least_1pair;
   }
-  bool attached_at_least_1pair {min_diff_from2a.maxCoeff() < buf_factor*tol};
-
-  string non_overlapped_test= {non_overlapped ? "Success" : "Failed"};
-  string attached_at_least_1pair_test= {attached_at_least_1pair ? "Success" : "Failed"};
-  cout << "non_overlapped_check: " << non_overlapped_test << endl;
-  cout << "attached_at_least_1pair_check: " << attached_at_least_1pair_test << endl;
-  if(non_overlapped && attached_at_least_1pair){
-    return pos_sph_agg_final;
+
+  if(non_overlapped_all && attached_at_least_1pair_all){
+    return pos_sph_agg_output;
   }else{
-    cout << "Geometry check failed! Please try other (kf, Df) pair." << endl;
+    cout << "Geometry check failed! Please change tol." << endl;
     exit(1);
   }
 
diff --git a/aggregate_gen_CCA.hpp b/aggregate_gen_CCA.hpp
@@ -1,5 +1,5 @@
 
 extern std::default_random_engine re;
-extern std::uniform_real_distribution<float> urf;
+extern std::uniform_real_distribution<double> urf;
 
-Eigen::MatrixXf aggregate_gen_CCA(const float& a, const int& num_sph_SA, const int& levels, const float& kf, const float& Df, const float& tol);
+std::vector<Eigen::MatrixXd> aggregate_gen_CCA(const double& a, const int& num_sph_SA, const int& levels, const double& kf, const double& Df, const double& tol);
diff --git a/aggregate_gen_SA.cpp b/aggregate_gen_SA.cpp
@@ -6,16 +6,16 @@ using namespace Eigen;
 using namespace std;
 
 //return pos_sph matrix
-MatrixXf aggregate_gen_SA(const float& a, const int& num_sph, const float& kf, const float& Df, const float& tol){
-    MatrixXf pos_sph= MatrixXf::Zero(3,num_sph);
+MatrixXd aggregate_gen_SA(const double& a, const int& num_sph, const double& kf, const double& Df, const double& tol){
+    MatrixXd pos_sph= MatrixXd::Zero(3,num_sph);
     pos_sph.col(0) << -a,0,0;
     pos_sph.col(1) << a,0,0;
 
-    Vector3f pos_cand(3);
-    float rN;
-    float phi;
-    float u;
-    ArrayXf dist_sph;
+    Vector3d pos_cand(3);
+    double rN;
+    double phi;
+    double u;
+    ArrayXd dist_sph;
 
     for(int N=3; N <= num_sph; ++N){
         rN= N*N*a*a/(N-1)*pow(N/kf,2/Df)-N*a*a/(N-1)-N*a*a*pow((N-1)/kf,2/Df);
@@ -32,7 +32,7 @@ MatrixXf aggregate_gen_SA(const float& a, const int& num_sph, const float& kf, c
                 dist_sph(i)=(pos_sph.col(i)-pos_cand).norm();
             }
 
-            if(dist_sph.minCoeff() < 2*a+tol && dist_sph.minCoeff() > 2*a-tol){
+            if(dist_sph.minCoeff() < 2*a+tol && dist_sph.minCoeff() > 2*a){
                 pos_sph.col(N-1)= pos_cand;
                 pos_sph.colwise() -= pos_sph.rowwise().mean(); // translate the centroid to origin
                 break;
@@ -41,7 +41,7 @@ MatrixXf aggregate_gen_SA(const float& a, const int& num_sph, const float& kf, c
     }
 
     // final check of the attached condition
-    ArrayXf mindist_c(num_sph);
+    ArrayXd mindist_c(num_sph);
     for(int i= 0; i< num_sph; ++i){
         mindist_c(i)=((pos_sph.colwise()-pos_sph.col(i)).colwise().norm().array()-2*a).array().abs().minCoeff();
     }
diff --git a/aggregate_gen_SA.hpp b/aggregate_gen_SA.hpp
@@ -1,5 +1,5 @@
 
 extern std::default_random_engine re;
-extern std::uniform_real_distribution<float> urf;
+extern std::uniform_real_distribution<double> urf;
 
-Eigen::MatrixXf aggregate_gen_SA(const float& a, const int& num_sph, const float& kf, const float& Df, const float& tol);
+Eigen::MatrixXd aggregate_gen_SA(const double& a, const int& num_sph, const double& kf, const double& Df, const double& tol);
diff --git a/aggregate_gen_main.cpp b/aggregate_gen_main.cpp
