implemented corr-proj with dih-shifts. needs testing.

moldyn · Dec 31, 2015 · 3aa947d · 3aa947d
1 parent 2cbed98
commit 3aa947d
Show file tree

Hide file tree

Showing 5 changed files with 335 additions and 325 deletions.
diff --git a/src/file_io.cpp b/src/file_io.cpp
@@ -121,68 +121,59 @@ namespace FastPCA {
       mem_buf_size /= 4;
       std::vector<double> means;
       std::vector<double> sigmas;
-      if ( (! use_dih_shifts) || use_correlation) {
+      std::vector<double> dih_shifts;
+      std::vector<double> scaled_periodicities;
+      // compute means
+      if ((!use_dih_shifts) || use_correlation) {
         std::tie(std::ignore, std::ignore, means) = FastPCA::Periodic::means(file_in, mem_buf_size);
       }
+      // compute sigmas
       if (use_correlation) {
         sigmas = FastPCA::Periodic::sigmas(file_in, mem_buf_size, means);
-        for (double& s: sigmas) {
-          s = 1.0 / s;
+        scaled_periodicities.resize(sigmas.size(), 2*M_PI);
+        for (std::size_t j=0; j < sigmas.size(); ++j) {
+          // invert sigmas for easier rescaling
+          sigmas[j] = 1.0 / sigmas[j];
+          scaled_periodicities[j] *= sigmas[j];
         }
       }
-
-
-
-
-
-
-//      DataFileReader<double> fh_file_in(file_in, mem_buf_size);
-//      DataFileWriter<double> fh_file_out(file_out);
-
-//TODO: dih-shifts for correlated data
-//        1. compute shifts for correlated data
-//        2. shift by means
-//        3. rescale
-//        4. shift by dih-shifts
-// template: 
-//            verbose && std::cerr << "computing optimal shifts for dihedrals" << std::endl;
-//            dih_shifts = FastPCA::Periodic::dih_shifts(file_input, mem_buf_size);
-//            if (verbose) {
-//              std::cout << "dih shifts:" << std::endl;
-//              for(auto s: dih_shifts) {
-//                std::cout << " " << s;
-//              }
-//              std::cout << std::endl;
-//            }
-
-
-
-
-
-
-// TODO: remove old code: 
-//      bool append_to_file = false;
-//      while ( ! fh_file_in.eof()) {
-//        Matrix<double> m = fh_file_in.next_block();
-//        if (m.n_rows() > 0) {
-//          if (use_dih_shifts) {
-//            // dih shifts: shift minima of dihedral regions to periodic barrier
-//            FastPCA::deg2rad_inplace(m);
-//            FastPCA::Periodic::shift_matrix_columns_inplace(m, dih_shifts);
-//          } else {
-//            // default behaviour: shift by periodic means
-//            FastPCA::deg2rad_inplace(m);
-//            FastPCA::Periodic::shift_matrix_columns_inplace(m, means);
-//          }
-//          if (use_correlation) {
-//            FastPCA::scale_matrix_columns_inplace(m, sigmas);
-//          }
-//          fh_file_out.write(m*eigenvecs, append_to_file);
-//          append_to_file = true;
-//        }
-//      }
+      // compute dih-shifts
+      if (use_dih_shifts) {
+        dih_shifts = FastPCA::Periodic::dih_shifts(file_in, mem_buf_size);
+        if (use_correlation) {
+          for (std::size_t i=0; i < dih_shifts.size(); ++i) {
+            //TODO: check signs
+            dih_shifts[i] = (means[i] - dih_shifts[i]) * sigmas[i];
+          }
+        }
+      }
+      // projections
+      bool append_to_file = false;
+      DataFileReader<double> fh_file_in(file_in, mem_buf_size);
+      DataFileWriter<double> fh_file_out(file_out);
+      read_blockwise(fh_file_in, [&](Matrix<double>& m) {
+        // convert degrees to radians
+        FastPCA::deg2rad_inplace(m);
+        if (use_correlation || (! use_dih_shifts)) {
+          // shift by periodic means
+          FastPCA::Periodic::shift_matrix_columns_inplace(m, means);
+        } else if (use_dih_shifts) {
+          // shift by optimal dih-shifts
+          FastPCA::Periodic::shift_matrix_columns_inplace(m, dih_shifts);
+        }
+        if (use_correlation) {
+          // scale data by sigmas for correlated projections
+          FastPCA::scale_matrix_columns_inplace(m, sigmas);
+          if (use_dih_shifts) {
+          //TODO: check if this is working
+            FastPCA::Periodic::shift_matrix_columns_inplace(m, dih_shifts, scaled_periodicities);
+          }
+        }
+        // output
+        fh_file_out.write(m*eigenvecs, append_to_file);
+        append_to_file = true;
+      });
     }
   } // end namespace FastPCA::Periodic
-
 } // end namespace FastPCA
 
diff --git a/src/file_io.hpp b/src/file_io.hpp
@@ -153,7 +153,6 @@ namespace FastPCA {
                                bool use_correlation,
                                bool use_dih_shifts);
   } // end namespace FastPCA::Periodic
-
 } // end namespace FastPCA
 
 // load template implementations