Changeset 7b67a3 for molecuilder/src/analyzer.cpp

Timestamp:

Aug 18, 2008, 8:35:11 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

e6971b

Parents:

d24e8f0

Message:

Splitting MatrrixContainer::ParseMatrix in ParseMatrix

ParseMatrix just parses only a given matrix into a given part of the MatrixContainer::matrix array).
ParseFragmentMatrix calls ParseMatrix for a sequence of numbered files containing the matrices.
ForceMatrix::ParseFragmentMatrix overrides the routine from MatrixContainer in order to handle the parsing of the last and one ma
trix in a special way
EnergyMatrix::ParseFragmentMatrix overrides the routine from MatrixContainer in order to handle the parsing of the last and one matrix in a special way
Realloc() was fixed: ReAlloc would break if given oldptr is NULL, but now we simply switch to doing a malloc instead and only admonishing the wrong call to ReAlloc instead of Malloc
analyser.cpp and joiner.cpp have been adapted to these new calling schemes

File:

• molecuilder/src/analyzer.cpp (modified) (7 diffs)

molecuilder/src/analyzer.cpp

@@ -29 +29 @@
   ForceMatrix Shielding;
   ForceMatrix ShieldingPAS;
-  EnergyMatrix Time;
+  MatrixContainer Time;
   EnergyMatrix EnergyFragments;
   EnergyMatrix HcorrectionFragments;
@@ -51 +51 @@
   bool Hcorrected = true;
   double norm;
-  int counter;
   
   cout << "ANOVA Analyzer" << endl;
@@ -83 +82 @@
   
   // ------------- Parse through all Fragment subdirs --------
-  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
-  Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
-  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
-  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
+  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
+  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
+  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
+  if (!Time.ParseMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
+    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
@@ -112 +111 @@
   
   // ---------- Parse fragment files created by 'joiner' into an array -------------
-  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
-  if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
-  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
+  if (!EnergyFragments.ParseMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
+  if (Hcorrected) HcorrectionFragments.ParseMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
+  if (!ForceFragments.ParseMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
   if (periode != NULL) { // also look for PAS values
-    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
-    if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
+    if (!ShieldingFragments.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
+    if (!ShieldingPASFragments.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
   }
 
@@ -161 +160 @@
   
   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
+  Time.SetLastMatrix(0., 0);
+  for (int BondOrder=KeySet.Order;BondOrder--;)
+    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
+      for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
+        for(int k=Time.ColumnCounter;k--;) {
+          Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
+        }
   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
-    for(int k=0;k<Time.ColumnCounter;k++) {
+    for(int k=0;k<Time.ColumnCounter;k++)
       output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
-    }
     output << endl;
   }
   output << endl;
   output.close();
-  for(int k=0;k<Time.ColumnCounter;k++)
-    Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
 
   // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
@@ -179 +182 @@
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
-  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
   if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
-  if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
-  for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
-    for(int k=Time.ColumnCounter;k--;) {
-      Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
-    }
-  counter = 0;
+  Time.SetLastMatrix(0., 0);
   output << "#Order\tFrag.No.\t" << Time.Header << endl;
-  output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
@@ -195 +191 @@
           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
         }
-    counter += KeySet.FragmentsPerOrder[BondOrder];
-    output << BondOrder+1 << "\t" << counter;
-    output2 << BondOrder+1 << "\t" << counter;
-    for(int k=0;k<Time.ColumnCounter;k++) {
-      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
-      if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
-        output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
-      else
-        output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
-    }
-    output << endl;
-    output2 << endl;
+    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
+    for(int k=0;k<Time.ColumnCounter;k++)
+      output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][1][k];
+    output << endl;
   }
   output.close();
-  output2.close();
-
 
   // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM