Skip to content

Commit

Permalink
implemented static shift for periodic boundaries as default
Browse files Browse the repository at this point in the history
  • Loading branch information
@lettis
lettis committed May 18, 2017
1 parent d8df689 commit 36edd35
Show file tree
Hide file tree
Showing 4 changed files with 324 additions and 105 deletions.
2 changes: 1 addition & 1 deletion src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ add_library(covariance covariance.cpp file_io.cpp util.cpp matrix.cpp errors.cpp

add_subdirectory(xdrfile)

add_executable (fastpca pca.cpp)
add_executable (fastpca fastpca.cpp)
target_link_libraries (fastpca covariance xdrfile ${Boost_LIBRARIES} ${LAPACK_LIBRARIES})

install(TARGETS fastpca RUNTIME DESTINATION .)
Expand Down
205 changes: 152 additions & 53 deletions src/pca.cpp → src/fastpca.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,15 @@ namespace b_po = boost::program_options;
int main(int argc, char* argv[]) {
bool verbose = false;
bool periodic = false;
bool dynamic_shift = false;
bool use_correlation = false;

b_po::options_description desc (std::string(argv[0]).append(
"\n\n"
"Calculate principle components from large data files.\n"
"Input data should be given as textfiles\n"
"with whitespace separated columns or, alternatively as GROMACS .xtc-files.\n"
"with whitespace separated columns or,"
" alternatively as GROMACS .xtc-files.\n"
"\n"
"options"));

Expand All @@ -60,11 +62,14 @@ int main(int argc, char* argv[]) {
("vec-in,V", b_po::value<std::string>()->default_value(""),
"input (optional): file with already computed eigenvectors")
("stats-in,S", b_po::value<std::string>()->default_value(""),
"input (optional): mean values, sigmas and boundary shifts (shifts only for periodic)."
" Provide this, if you want to project new data onto a previously computed principal space."
" If you do not define the stats of the previous run, means, sigmas and shifts"
" will be re-computed and the resulting projections will not be comparable"
" to the previous ones.")
"input (optional): mean values, sigmas and boundary shifts"
" (shifts only for periodic)."
" Provide this, if you want to project new data onto a previously"
" computed principal space."
" If you do not define the stats of the previous run, means,"
" sigmas and shifts"
" will be re-computed and the resulting projections"
" will not be comparable to the previous ones.")
// output
("proj,p", b_po::value<std::string>()->default_value(""),
"output (optional): file for projected data")
Expand All @@ -73,22 +78,27 @@ int main(int argc, char* argv[]) {
("vec,v", b_po::value<std::string>()->default_value(""),
"output (optional): file for eigenvectors")
("stats,s", b_po::value<std::string>()->default_value(""),
"output (optional): mean values, sigmas and boundary shifts (shifts only for periodic)")
"output (optional): mean values, sigmas and boundary shifts"
" (shifts only for periodic)")
("val,l", b_po::value<std::string>()->default_value(""),
"output (optional): file for eigenvalues")
// parameters
("norm,N", b_po::value(&use_correlation)->zero_tokens(),
"if set, use correlation instead of covariance by normalizing input (default: false)")
"if set, use correlation instead of covariance by"
" normalizing input (default: false)")
("buf,b", b_po::value<std::size_t>()->default_value(2),
"max. allocatable RAM [Gb] (default: 2); WARNING there is some error in the code "
"that leads to twice mem consumption at some point. use only HALF of what you would "
"normally use. will be fixed soon")
"max. allocatable RAM [Gb] (default: 2)")
("periodic,P", b_po::value(&periodic)->zero_tokens(),
"compute covariance and PCA on a torus (i.e. for periodic data like dihedral angles)")
"compute covariance and PCA on a torus"
" (i.e. for periodic data like dihedral angles)")
("dynamic-shift,D", b_po::value(&dynamic_shift)->zero_tokens(),
"use dynamic shifting for periodic projection correction."
" (default: fale, i.e. simply shift to region of lowest density)")
("verbose", b_po::value(&verbose)->zero_tokens(),
"verbose mode (default: not set)")
("nthreads,n", b_po::value<std::size_t>()->default_value(0),
"number of OpenMP threads to use. if set to zero, will use value of OMP_NUM_THREADS (default: 0)");
"number of OpenMP threads to use."
" if set to zero, will use value of OMP_NUM_THREADS (default: 0)");

b_po::variables_map args;
try {
Expand Down Expand Up @@ -140,94 +150,183 @@ int main(int argc, char* argv[]) {
//// input
if (input_stats_file_given) {
// -S ?
verbose && std::cerr << "loading stats (shifts, means, sigmas) from file" << std::endl;
FastPCA::DataFileReader<double> fh_shifts(args["stats-in"].as<std::string>());
verbose
&& std::cerr << "loading stats (shifts, means, sigmas) from file"
<< std::endl;
FastPCA::DataFileReader<double>
fh_shifts(args["stats-in"].as<std::string>());
stats = fh_shifts.next_block();
} else {
verbose && std::cerr << "computing stats (means, sigmas, ...)" << std::endl;
verbose
&& std::cerr << "computing stats (means, sigmas, ...)"
<< std::endl;
if (periodic) {
stats = FastPCA::Periodic::stats(file_input, mem_buf_size);
verbose
&& dynamic_shift
&& std::cerr << "using dynamic shifting";
verbose
&& ( ! dynamic_shift)
&& std::cerr << "using static shifting";
verbose && std::cerr <<" for periodic barrier correction"
<< std::endl;
stats = FastPCA::Periodic::stats(file_input
, mem_buf_size
, dynamic_shift);
} else {
stats = FastPCA::stats(file_input, mem_buf_size);
stats = FastPCA::stats(file_input
, mem_buf_size);
}
}
if (input_eigenvec_file_given && input_covmat_file_given) {
// -V and -C ?
std::cerr << "error: Providing the covariance matrix and the eigenvectors does not make sense." << std::endl
<< " If we have the vectors, we do not need the covariance matrix," << std::endl
<< " therefore the options '-C' and '-V' are mutually exclusive." << std::endl;
std::cerr << "error: Providing the covariance matrix and"
" the eigenvectors does not make sense."
<< std::endl
<< " If we have the vectors, we do not need"
" the covariance matrix,"
<< std::endl
<< " therefore the options '-C' and '-V'"
" are mutually exclusive."
<< std::endl;
return EXIT_FAILURE;
} else if (input_eigenvec_file_given) {
// -V ?
verbose && std::cerr << "loading eigenvectors from file" << std::endl;
FastPCA::DataFileReader<double> fh_vecs(args["vec-in"].as<std::string>());
verbose && std::cerr << "loading eigenvectors from file"
<< std::endl;
FastPCA::DataFileReader<double>
fh_vecs(args["vec-in"].as<std::string>());
vecs = fh_vecs.next_block();
if (projection_file_given && (! input_stats_file_given)) {
std::cerr << "warning: You specified a projection output and gave" << std::endl
<< " pre-computed eigenvectors." << std::endl
<< " Projection is, however, centered to the mean." << std::endl
<< " Without stats (means, sigma, [dihedral shifts])," << std::endl
<< " you will in general get non-comparable results" << std::endl
<< " to previous projections!" << std::endl;
std::cerr << "warning: You specified a projection output and gave"
<< std::endl
<< " pre-computed eigenvectors."
<< std::endl
<< " Projection is, however,"
" centered to the mean."
<< std::endl
<< " Without stats (means, sigma,"
" [dihedral shifts]),"
<< std::endl
<< " you will in general get"
" non-comparable results"
<< std::endl
<< " to previous projections!"
<< std::endl;
}
} else {
if (input_covmat_file_given) {
// -C ?
verbose && std::cerr << "loading covariance matrix from file" << std::endl;
FastPCA::DataFileReader<double> cov_in(args["cov-in"].as<std::string>());
cov = FastPCA::SymmetricMatrix<double>(cov_in.next_block(cov_in.n_cols()));
verbose
&& std::cerr << "loading covariance matrix from file"
<< std::endl;
FastPCA::DataFileReader<double>
cov_in(args["cov-in"].as<std::string>());
cov = FastPCA::SymmetricMatrix<double>(
cov_in.next_block(cov_in.n_cols()));
} else {
if (periodic) {
verbose && ( ! use_correlation) && std::cerr << "constructing covariance matrix for periodic data" << std::endl;
verbose && use_correlation && std::cerr << "constructing correlation matrix for periodic data" << std::endl;
cov = FastPCA::Periodic::covariance_matrix(file_input, mem_buf_size, use_correlation, stats);
verbose
&& ( ! use_correlation)
&& std::cerr << "constructing covariance matrix"
" for periodic data"
<< std::endl;
verbose
&& use_correlation
&& std::cerr << "constructing correlation matrix"
" for periodic data"
<< std::endl;
cov = FastPCA::Periodic::covariance_matrix(file_input
, mem_buf_size
, use_correlation
, stats);
} else {
verbose && ( ! use_correlation) && std::cerr << "constructing covariance matrix" << std::endl;
verbose && use_correlation && std::cerr << "constructing correlation matrix" << std::endl;
cov = FastPCA::covariance_matrix(file_input, mem_buf_size, use_correlation, stats);
verbose
&& ( ! use_correlation)
&& std::cerr << "constructing covariance matrix"
<< std::endl;
verbose
&& use_correlation
&& std::cerr << "constructing correlation matrix"
<< std::endl;
cov = FastPCA::covariance_matrix(file_input
, mem_buf_size
, use_correlation
, stats);
}
}
vecs = cov.eigenvectors();
}
//// output
// -c ?
if (covmat_file_given) {
verbose && ( ! use_correlation) && std::cerr << "writing covariance matrix" << std::endl;
verbose && use_correlation && std::cerr << "writing correlation matrix" << std::endl;
verbose
&& ( ! use_correlation)
&& std::cerr << "writing covariance matrix"
<< std::endl;
verbose
&& use_correlation
&& std::cerr << "writing correlation matrix"
<< std::endl;
std::string covmat_file = args["cov"].as<std::string>();
FastPCA::DataFileWriter<double>(covmat_file).write(FastPCA::Matrix<double>(cov));
FastPCA::DataFileWriter<double>(covmat_file)
.write(FastPCA::Matrix<double>(cov));
}
// -v ?
if (eigenvec_file_given) {
verbose && std::cerr << "solving eigensystem/writing eigenvectors matrix" << std::endl;
verbose
&& std::cerr << "solving eigensystem/writing eigenvectors matrix"
<< std::endl;
std::string eigenvec_file = args["vec"].as<std::string>();
FastPCA::DataFileWriter<double>(eigenvec_file).write(cov.eigenvectors());
FastPCA::DataFileWriter<double>(eigenvec_file)
.write(cov.eigenvectors());
}
// -l ?
if (eigenval_file_given) {
verbose && std::cerr << "solving eigensystem/writing eigenvalues matrix" << std::endl;
verbose
&& std::cerr << "solving eigensystem/writing eigenvalues matrix"
<< std::endl;
std::string eigenval_file = args["val"].as<std::string>();
FastPCA::DataFileWriter<double>(eigenval_file).write(cov.eigenvalues());
FastPCA::DataFileWriter<double>(eigenval_file)
.write(cov.eigenvalues());
}
// -p ?
if (projection_file_given) {
std::string projection_file = args["proj"].as<std::string>();
if (periodic) {
verbose && std::cerr << "computing projections for periodic data" << std::endl;
FastPCA::Periodic::calculate_projections(file_input, projection_file, vecs, mem_buf_size, use_correlation, stats);
verbose
&& std::cerr << "computing projections for periodic data"
<< std::endl;
FastPCA::Periodic::calculate_projections(file_input
, projection_file
, vecs
, mem_buf_size
, use_correlation
, stats);
} else {
verbose && std::cerr << "computing projections" << std::endl;
FastPCA::calculate_projections(file_input, projection_file, vecs, mem_buf_size, use_correlation, stats);
verbose
&& std::cerr << "computing projections"
<< std::endl;
FastPCA::calculate_projections(file_input
, projection_file
, vecs
, mem_buf_size
, use_correlation
, stats);
}
}
// -s ?
if (stats_file_given) {
verbose && std::cerr << "writing stats" << std::endl;
FastPCA::DataFileWriter<double>(args["stats"].as<std::string>()).write(stats);
verbose
&& std::cerr << "writing stats"
<< std::endl;
FastPCA::DataFileWriter<double>(args["stats"].as<std::string>())
.write(stats);
}
} else {
std::cerr << "please specify at least one output file!" << std::endl;
std::cerr << desc << std::endl;
std::cerr << "please specify at least one output file!"
<< std::endl
<< desc
<< std::endl;
return EXIT_FAILURE;
}
} else {
Expand Down
Loading

0 comments on commit 36edd35

Please sign in to comment.