Following of the second_generation for the

Working Betti numbers
PRISBteam · May 8, 2023 · 7495635 · 7495635
1 parent e1c6956
commit 7495635
Show file tree

Hide file tree

Showing 10 changed files with 412 additions and 158 deletions.
diff --git a/config/characterisation.ini b/config/characterisation.ini
@@ -13,13 +13,13 @@ S_mean = 0
 S_skew = 0
 
 [edges_lab]
-el_active = 1
-config_entropy = 1
-S_mean = 1
-S_skew = 1
-j_fractions = 1
-d_fractions = 1
-analytical = 1
+el_active = 0
+config_entropy = 0
+S_mean = 0
+S_skew = 0
+j_fractions = 0
+d_fractions = 0
+analytical = 0
 
 [nodes_lab]
 nl_active = 0
@@ -28,7 +28,8 @@ S_mean = 0
 S_skew = 0
 
 [spectra_lab]
-laplacians = 0
+calc_steps_numb = 50
+laplacians = 1
 laplacians_spectra = 0
 laplacians_betti = 0
 

diff --git a/config/main.ini b/config/main.ini
@@ -11,7 +11,7 @@ PCC_Writer=ON
 
 [general]
 dim=3
-source_dir=/Users/user/Dropbox/OFFICE/NEPER/Kinetic_PCC_examples_/Voronoi_1k/
-;source_dir=/Users/user/Dropbox/OFFICE/EGitHub/_Kinetic_PCC_examples/Siying500Complex/
-output_dir=/Users/user/Dropbox/OFFICE/NEPER/resultsFebruary2023/v1k_processing/
-;output_dir=/Users/user/Dropbox/OFFICE/NEPER/resultsMay2023/500cells_processing/
+; source_dir=/Users/user/Dropbox/OFFICE/NEPER/Kinetic_PCC_examples_/Voronoi_1k/
+source_dir=/Users/user/Dropbox/OFFICE/EGitHub/_Kinetic_PCC_examples/Siying500Complex/
+;output_dir=/Users/user/Dropbox/OFFICE/NEPER/resultsFebruary2023/v1k_processing/
+output_dir=/Users/user/Dropbox/OFFICE/NEPER/resultsMay2023/500cells_processing/
diff --git a/config/processing.ini b/config/processing.ini
@@ -14,7 +14,7 @@ pmax_fraction3 = 0.0
 
 [faces]
 face_types_number = 1
-pf_mode = S
+pf_mode = R
 ; "S" - read from source, "R" or "F" -> Smax, "D" -> Smin (max S_deviatoric (skew) part) or "L"
 source = /Users/user/Dropbox/OFFICE/NEPER/resultsFebruary2023/v1k_processing/Random099/s_cells_sequence.txt
 pf_index = 0

diff --git a/config/writer.ini b/config/writer.ini
@@ -1,25 +1,27 @@
 ; writer specifications
 
 [sequences]
-isSequencesOutput = 1
-isDesignvectorsOutput = 1
+isSequencesOutput = 0
+isDesignvectorsOutput = 0
 
 [entropic_polyhedrons]
 isPolyhedronFractions = 0
 
-
 [entropic_faces]
 isFaceFractions = 0
 isConfEntropy = 0
 
 [entropic_edges]
-isConfEntropy = 1
-isFractions = 1
-isDegreeFractions = 1
+isConfEntropy = 0
+isFractions = 0
+isDegreeFractions = 0
 
 [entropic_nodes]
 isNodeFractions = 0
 
 [entropic_analytical]
-isEdgeFractions = 1
-isEdgeConfEntropies = 1
+isEdgeFractions = 0
+isEdgeConfEntropies = 0
+
+[component_analysis]
+isBetti = 1
diff --git a/src/lib/PCC_Characterisation/LaplaciansLab.h b/src/lib/PCC_Characterisation/LaplaciansLab.h
diff --git a/src/lib/PCC_Characterisation/PCC_StructureCharacterisation.h b/src/lib/PCC_Characterisation/PCC_StructureCharacterisation.h
@@ -10,8 +10,8 @@
 ProcessedComplex PCC_StructureCharacterisation(CellsDesign &new_cells_design) {
     ProcessedComplex PCC_characteristics; // module output
 /// Read simulation configuration from file :: the number of special face types and calculating parameters. Then Output of the current configuration to the screen
-    std::vector<int> charlabs_polyhedrons, charlabs_faces, charlabs_edges, charlabs_nodes;
-    std::vector<double> ConfigVector = config_reader_characterisation(source_path, charlabs_polyhedrons, charlabs_faces, charlabs_edges, charlabs_nodes, Out_logfile_stream); // vector<double> (!) from ini_readers.h library
+    std::vector<int> charlabs_polyhedrons, charlabs_faces, charlabs_edges, charlabs_nodes, charlabs_laplacians;
+    std::vector<double> ConfigVector = config_reader_characterisation(source_path, charlabs_polyhedrons, charlabs_faces, charlabs_edges, charlabs_nodes, charlabs_laplacians, Out_logfile_stream); // vector<double> (!) from ini_readers.h library
 
 /// # 1 # Creation of the State Vectors
 // Faces
@@ -26,49 +26,119 @@ ProcessedComplex PCC_StructureCharacterisation(CellsDesign &new_cells_design) {
 
     } // end for(auto fseq...)
 
-    for (int i = 0; i < 4; ++i) { // for all types of cells
+    for (int i = 0; i < 4; ++i) { /// for all types of cells
+
+        /// # 2 # Measures: fractions and entropies
+        if(i == 1 && charlabs_edges.size() > 0) { // Edges lab
+        if (charlabs_edges.at(0) == 1) {
+                vector<int> EdgeTypes_char(CellNumbs.at(1 + (dim - 3)),
+                                           0); // vector<int> in the form [ 0 2 3 3 2 1 ...] with the TJs type ID as its values
+                std::vector<double> j_edge_fractions(dim + 1, 0), d_edge_fractions(dim,
+                                                                                   0); // fractions of (1) edges of different types and (2) edges of different degrees
+
+                /// for each state:
+                tuple<double, double> conf_entropy_t;
+                for (auto current_seq: PCC_characteristics.face_process_seq) {
+                    EdgeTypes_char.clear();
+                    std::fill(j_edge_fractions.begin(), j_edge_fractions.end(), 0);
+                    std::fill(d_edge_fractions.begin(), d_edge_fractions.end(), 0);
+
+                    EdgeTypes_char = Edge_types_byFaces(CellNumbs, current_seq, j_edge_fractions, d_edge_fractions);
+
+                    conf_entropy_t = Configuration_Entropy_tuple(j_edge_fractions); // conf entropy
+
+                    PCC_characteristics.e_entropy_mean_vector.push_back(std::get<0>(
+                            conf_entropy_t)); // std::tuple<double, double> Configuration_Entropy_tuple(std::vector<double> const &j_fractions) based on Edges fraction vector
+                    PCC_characteristics.e_entropy_skrew_vector.push_back(std::get<1>(conf_entropy_t)); // tuple
+                    //// MUST BE IMPROVED (!) :
+                    PCC_characteristics.e_entropy_full_vector.push_back(
+                            std::get<0>(conf_entropy_t) + std::get<1>(conf_entropy_t));
+
+                    PCC_characteristics.je_fractions_vector.push_back(j_edge_fractions);
+                    PCC_characteristics.de_fractions_vector.push_back(d_edge_fractions);
+                } // for (auto current_seq : PCC_characteristics.face_process_seq)
+
+                if (charlabs_edges.at(4) == 1) { // analytical
+                    /// Measures: Analytical solutions for edge fractions and configuration entropies
+                    std::vector<tuple<double, double>> AnalyticalRandEntropies, AnalyticalCrystEntropies;
+                    PCC_characteristics.j_analytical_rand_vector = TJsAnalytics_random(CellNumbs.at(2 + (dim - 3)),
+                                                                                       AnalyticalRandEntropies);
+                    PCC_characteristics.AnRandEntropies_vector = AnalyticalRandEntropies;
+
+                    PCC_characteristics.j_analytical_cryst_vector = TJsAnalytics_crystallography(
+                            CellNumbs.at(2 + (dim - 3)), AnalyticalCrystEntropies);
+                    PCC_characteristics.AnCrystEntropies_vector = AnalyticalCrystEntropies;
+
+                    // analytical edge degree fractions
+                    PCC_characteristics.d_analytical_rand_vector = TJDsAnalytics_random(CellNumbs.at(2 + (dim - 3)));
+                    PCC_characteristics.d_analytical_cryst_vector = TJDsAnalytics_crystallography(
+                            CellNumbs.at(2 + (dim - 3)));
+                } // enf if(charlabs_edges.at(4) == 1)
+
+            } //  end of if(charlabs_edges.at(0) == 1)
+            } //  end of if(i == 1 && charlabs_edges.at(0) == 1) { // Edges lab
 
-    /// # 2 # Measures: fractions and entropies
-        if(i == 1 && charlabs_edges.at(0) == 1) { // Edges lab
-            vector<int> EdgeTypes_char(CellNumbs.at(1 + (dim - 3)), 0); // vector<int> in the form [ 0 2 3 3 2 1 ...] with the TJs type ID as its values
-            std::vector<double> j_edge_fractions(dim+1,0), d_edge_fractions(dim,0); // fractions of (1) edges of different types and (2) edges of different degrees
-
-            /// for each state:
-            tuple<double,double> conf_entropy_t;
-            for (auto current_seq : PCC_characteristics.face_process_seq) {
-                EdgeTypes_char.clear();
-                std::fill(j_edge_fractions.begin(), j_edge_fractions.end(),0);
-                std::fill(d_edge_fractions.begin(), d_edge_fractions.end(),0);
-
-                EdgeTypes_char = Edge_types_byFaces(CellNumbs, current_seq, j_edge_fractions, d_edge_fractions);
-
-                conf_entropy_t = Configuration_Entropy_tuple(j_edge_fractions); // conf entropy
-
-                PCC_characteristics.e_entropy_mean_vector.push_back(std::get<0>(conf_entropy_t)); // std::tuple<double, double> Configuration_Entropy_tuple(std::vector<double> const &j_fractions) based on Edges fraction vector
-                PCC_characteristics.e_entropy_skrew_vector.push_back(std::get<1>(conf_entropy_t)); // tuple
-                //// MUST BE IMPROVED (!) :
-                PCC_characteristics.e_entropy_full_vector.push_back( std::get<0>(conf_entropy_t) + std::get<1>(conf_entropy_t));
-
-                PCC_characteristics.je_fractions_vector.push_back(j_edge_fractions);
-                PCC_characteristics.de_fractions_vector.push_back(d_edge_fractions);
-            }
-
-            if(charlabs_edges.at(4) == 1) { // analytical
-            /// Measures: Analytical solutions for edge fractions and configuration entropies
-                std::vector<tuple<double, double>> AnalyticalRandEntropies, AnalyticalCrystEntropies;
-                PCC_characteristics.j_analytical_rand_vector = TJsAnalytics_random(CellNumbs.at(2 + (dim - 3)), AnalyticalRandEntropies);
-                PCC_characteristics.AnRandEntropies_vector = AnalyticalRandEntropies;
-
-                PCC_characteristics.j_analytical_cryst_vector = TJsAnalytics_crystallography(CellNumbs.at(2 + (dim - 3)), AnalyticalCrystEntropies);
-                PCC_characteristics.AnCrystEntropies_vector = AnalyticalCrystEntropies;
-
-                // analytical edge degree fractions
-                PCC_characteristics.d_analytical_rand_vector = TJDsAnalytics_random(CellNumbs.at(2 + (dim - 3)));
-                PCC_characteristics.d_analytical_cryst_vector = TJDsAnalytics_crystallography(CellNumbs.at(2 + (dim - 3)));
-            } // enf if(charlabs_edges.at(4) == 1)
-
-            } //         if(i == 1 && charlabs_edges.at(0) == 1) { // Edges lab
     } // end for (int i = 0; i < 4; ++i)
 
+    /// # 3 # Laplacians of specail cells
+    if(charlabs_laplacians.at(1) == 1) {
+// AN - Nodes (Vertices or 0-Cells) sparse adjacency matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values
+// AE - Edges (Triple Junctions or 1-Cells) sparse adjacency matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values
+// AF - Faces (Grain Boundaries or 2-Cells) sparse adjacency matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values
+// AG - Grains (Volumes or 3-Cells) sparse adjacency matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values
+// MEN - Edges - Nodes (1-Cells to 0-Cells) sparse incidence matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values
+// MFE - Faces - Edges (2-Cells to 1-Cells) sparse incidence matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values
+// MGF - Grains - Faces (3-Cells to 2-Cells) sparse incidence matrix  in the form of three column {i, j, value}, where i and j are the indices of elements with non-zero values // odir - source path
+
+/// Adjacency sparse matrix for nodes /// Adjacency sparse matrix for edges /// Adjacency sparse matrix for faces /// Adjacency sparse matrix for grains
+//    SpMat ANS(CellNumbs.at(0), CellNumbs.at(0)), AES(CellNumbs.at(1 + (dim - 3)), CellNumbs.at(1 + (dim - 3))), AFS(CellNumbs.at(2 + (dim - 3)), CellNumbs.at(2 + (dim - 3))), AGS(CellNumbs.at(3 + (dim - 3)), CellNumbs.at(3 + (dim - 3)));
+//    ANS = SMatrixReader(paths.at(0), (CellNumbs.at(0)), (CellNumbs.at(0))); //all Nodes
+//    ANS = 0.5 * (ANS + SparseMatrix<double>(ANS.transpose())); // Full matrix instead of triagonal
+//    AES = SMatrixReader(paths.at(1 + (dim - 3)), (CellNumbs.at(1 + (dim - 3))), (CellNumbs.at(1 + (dim - 3)))); //all Edges
+//    AES = 0.5 * (AES + SparseMatrix<double>(AES.transpose())); // Full matrix instead of triagonal
+//    AFS = SMatrixReader(paths.at(2 + (dim - 3)), (CellNumbs.at(2 + (dim - 3))), (CellNumbs.at(2 + (dim - 3)))); //all Faces
+//    AFS = 0.5 * (AFS + SparseMatrix<double>(AFS.transpose())); // Full matrix instead of triagonal
+//    AGS = SMatrixReader(paths.at(3 + (dim - 3)), (CellNumbs.at(3 + (dim - 3))), (CellNumbs.at(3 + (dim - 3)))); //all Volumes
+//    AGS = 0.5 * (AGS + SparseMatrix<double>(AGS.transpose())); // Full matrix instead of triagonal
+
+/// Incidence sparse matrix for Edges and Nodes /// Incidence sparse matrix for Faces and Edges /// Incidence sparse matrix for Grains and Faces
+        SpMat ENS(CellNumbs.at(0), CellNumbs.at(1)), FES(CellNumbs.at(1 + (dim - 3)), CellNumbs.at(2 + (dim - 3))), GFS(
+                CellNumbs.at(2 + (dim - 3)), CellNumbs.at(3 + (dim - 3)));
+
+        ENS = SMatrixReader(paths.at(4 + (dim - 3)), (CellNumbs.at(0)), (CellNumbs.at(1))); //all Nodes-Edges
+        FES = SMatrixReader(paths.at(5 + (dim - 3)), (CellNumbs.at(1 + (dim - 3))),
+                            (CellNumbs.at(2 + (dim - 3)))); //all Edges-Faces
+        GFS = SMatrixReader(paths.at(6 + (dim - 3)), (CellNumbs.at(2 + (dim - 3))),
+                            (CellNumbs.at(3 + (dim - 3)))); //all Faces-Grains
+
+        double number_of_steps = (double) charlabs_laplacians.at(0); // reading from Characterisation.ini file (!)
+        unsigned int d_seq = std::floor((double) PCC_characteristics.face_process_seq.size() / number_of_steps);
+        std::vector<unsigned int> new_current_seq;
+        std::vector<double> new_Betti_numbers;
+
+        std::vector<unsigned int> Betti_calc_time;
+
+#pragma omp parallel for // parallel execution by OpenMP
+        for (unsigned int i = 0; i <= number_of_steps; ++i) {
+            new_current_seq = PCC_characteristics.face_process_seq.at(i*d_seq);
+            cout << "---------------------------------------------------------------------------" << endl;
+            cout << "[CHAR] Current special faces fraction: " << new_current_seq.size()/ (double) CellNumbs.at(2 + (dim - 3)) << endl;
+            Out_logfile_stream << "---------------------------------------------------------------------------" << endl;
+            Out_logfile_stream << "[CHAR] Current special faces fraction: " << new_current_seq.size()/ (double) CellNumbs.at(2 + (dim - 3)) << endl;
+            // Elapsing time for 3 Betti numbers
+            clock_t t;
+            Betti_calc_time.push_back(clock());
+
+            new_Betti_numbers = OperatorsBetti(new_current_seq, ENS, FES, GFS);
+            cout << "[CHAR] p & Betti numbers:  " << new_current_seq.size()/ (double) CellNumbs.at(2) << "  " << new_Betti_numbers.at(0) << "  " << new_Betti_numbers.at(1) << "  " << new_Betti_numbers.at(2) << endl;
+            Out_logfile_stream << "[CHAR] p & Betti numbers:  " << new_current_seq.size()/ (double) CellNumbs.at(2) << "  " << new_Betti_numbers.at(0) << "  " << new_Betti_numbers.at(1) << "  " << new_Betti_numbers.at(2) << endl;
+            //
+            Betti_calc_time.back() = clock() - Betti_calc_time.back();
+            cout << "[CHAR] Calculation time [s] of the Betti numbers is equal to   " << Betti_calc_time.back()/ pow(10.0,6.0)<< endl;
+            Out_logfile_stream << "[CHAR] Calculation time [s] of the Betti numbers is equal to   " << Betti_calc_time.back()/ pow(10.0,6.0)<< endl;
+
+            PCC_characteristics.Betti_vector.push_back({new_current_seq.size()/(double) CellNumbs.at(2), new_Betti_numbers.at(0), new_Betti_numbers.at(1), new_Betti_numbers.at(2)});
+        }
+    } // end of  for (int i = 0; i < number_of_steps, ++i)
     return PCC_characteristics;
 }
diff --git a/src/lib/PCC_Objects.h b/src/lib/PCC_Objects.h
@@ -695,18 +695,22 @@ class ProcessedComplex { // Essential for Characterisation module
         pcc_design = pcc_design;
     }
 
+    // Sequences of special k-cells
     std::vector<vector<unsigned int>> face_process_seq;
     std::vector<vector<int>> face_process_state;
 
+    // Entropic analysis
     std::vector<double> e_entropy_mean_vector, e_entropy_skrew_vector, e_entropy_full_vector;
-
     std::vector<vector<double>> je_fractions_vector, de_fractions_vector;
 
-    // analytical solutions
+    // Analytical solutions
     std::vector<vector<double>> j_analytical_rand_vector, d_analytical_rand_vector;
     std::vector<vector<double>> j_analytical_cryst_vector, d_analytical_cryst_vector;
     std::vector<tuple<double, double>> AnRandEntropies_vector, AnCrystEntropies_vector;
 
+    // Laplacian lab
+    std::vector<std::vector<double>> Betti_vector;
+
 }; /// end of class ProcessedComplex
 
 #endif

diff --git a/src/lib/PCC_Writer/PCC_Writer.h b/src/lib/PCC_Writer/PCC_Writer.h
@@ -16,9 +16,14 @@ config_reader_writer(source_path, writer_specifications, Out_logfile_stream); //
 int output_counter = 0; // special counter for output numeration
 
     /// Output special and ordinary face sequences to the output directory specified in config.txt
-    PCC_CellSequences_Writer(new_cells_design, output_counter);
-    PCC_Entropic_Writer(pcc_processed, output_counter);
-    PCC_AnalyticalFractions_Writer(pcc_processed, output_counter); // analytical solutions
+    if (writer_specifications.at(0) == 1)
+        PCC_CellSequences_Writer(new_cells_design, output_counter);
+    if (writer_specifications.at(2) == 1)
+        PCC_Entropic_Writer(pcc_processed, output_counter);
+    if (writer_specifications.at(5) == 1)
+        PCC_AnalyticalFractions_Writer(pcc_processed, output_counter); // analytical solutions
+    if(writer_specifications.at(7) == 1)
+        PCC_Laplacians_Writer(pcc_processed, output_counter);
 
 /*
        char* stype = const_cast<char*>(P_type.c_str()); // Processing type variable

diff --git a/src/lib/PCC_Writer/Writer_Functions.h b/src/lib/PCC_Writer/Writer_Functions.h
@@ -175,6 +175,19 @@ void PCC_AnalyticalFractions_Writer(ProcessedComplex &pcc_processed, int &output
 
 } // END of PCC_AnalyticalFractions_Writer( )
 
+void PCC_Laplacians_Writer(ProcessedComplex &pcc_processed, int &output_counter) {
+    ofstream OutBettiFile; // Betti numbers output
+    string betti_odir = output_dir + "Betti numbers.txt"s; // output directory
+
+    for (auto BV : pcc_processed.Betti_vector)
+        OutBettiFile << BV.at(0) << "  " << BV.at(1) << "  " << BV.at(2) << "  " << BV.at(3) << endl;
+
+    OutBettiFile.close();
+
+    cout << "(" << output_counter++ << ")  " << "Betti numbers of the special cell structure evolution has been successfully written in " << betti_odir << endl;
+    Out_logfile_stream << "(" << output_counter << ")  " << "Betti numbers of the special cell structure evolution has been successfully written in " << betti_odir << endl;
+} // END of PCC_Laplacians_Writer()
+
 /// * NEW HEAP * ///
 
 /**