From 3aa947db7d5ba7f8bce86d0ecb402ddda0629e8c Mon Sep 17 00:00:00 2001 From: Florian Sittel Date: Thu, 31 Dec 2015 17:06:13 +0100 Subject: [PATCH] implemented corr-proj with dih-shifts. needs testing. --- src/file_io.cpp | 101 +++++----- src/file_io.hpp | 1 - src/file_io.hxx | 516 ++++++++++++++++++++++++------------------------ src/util.hpp | 23 ++- src/util.hxx | 19 +- 5 files changed, 335 insertions(+), 325 deletions(-) diff --git a/src/file_io.cpp b/src/file_io.cpp index 2196715..53a4833 100644 --- a/src/file_io.cpp +++ b/src/file_io.cpp @@ -121,68 +121,59 @@ namespace FastPCA { mem_buf_size /= 4; std::vector means; std::vector sigmas; - if ( (! use_dih_shifts) || use_correlation) { + std::vector dih_shifts; + std::vector 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 fh_file_in(file_in, mem_buf_size); -// DataFileWriter 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 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 fh_file_in(file_in, mem_buf_size); + DataFileWriter fh_file_out(file_out); + read_blockwise(fh_file_in, [&](Matrix& 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 index eaa9a36..ecd3cbf 100644 --- 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 diff --git a/src/file_io.hxx b/src/file_io.hxx index 69983e6..9ecf406 100644 --- a/src/file_io.hxx +++ b/src/file_io.hxx @@ -30,297 +30,297 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace FastPCA { -namespace { // begin local namespace -namespace UseCSV { - template - Matrix read_chunk(LTS::AsciiParser& parser, std::size_t n_lines) { - std::size_t n_read_lines = n_lines; - std::vector m = std::move(parser.next_n_lines_continuous(n_read_lines, - LTS::AsciiParser::COL_MAJOR)); - return Matrix(m, n_read_lines, parser.n_cols()); - } - - template - void write_vector(std::ofstream& fh, std::vector v, bool in_line) { - if (fh.good()) { - //int n_digits = std::numeric_limits::digits10; - int n_digits = 6; - fh.precision(n_digits); - fh.setf(std::ios::scientific); - for (auto it=v.begin(); it != v.end(); ++it) { - fh.width(n_digits+10); // 8 extra chars for spacing, +/- and scientific formatting + namespace { // begin local namespace + namespace UseCSV { + template + Matrix read_chunk(LTS::AsciiParser& parser, std::size_t n_lines) { + std::size_t n_read_lines = n_lines; + std::vector m = std::move(parser.next_n_lines_continuous(n_read_lines, + LTS::AsciiParser::COL_MAJOR)); + return Matrix(m, n_read_lines, parser.n_cols()); + } + + template + void write_vector(std::ofstream& fh, std::vector v, bool in_line) { + if (fh.good()) { + //int n_digits = std::numeric_limits::digits10; + int n_digits = 6; + fh.precision(n_digits); + fh.setf(std::ios::scientific); + for (auto it=v.begin(); it != v.end(); ++it) { + fh.width(n_digits+10); // 8 extra chars for spacing, +/- and scientific formatting + if (in_line) { + fh << " "; + } + fh << (double) *it; + if ( ! in_line) { + fh << std::endl; + } + } if (in_line) { - fh << " "; + fh << std::endl; } - fh << (double) *it; - if ( ! in_line) { + } else { + throw FileNotFoundError(); + } + } + + template + void write_vector(std::string file_out, std::vector v, bool append) { + std::ofstream fh; + if (append) { + fh.open(file_out, std::ios::out | std::ios::app); + } else { + fh.open(file_out, std::ios::out); + } + write_vector(fh, v); + } + + template + void write_matrix(std::string file_out, Matrix m, bool append) { + std::ofstream fh; + if (append) { + fh.open(file_out, std::ios::out | std::ios::app); + } else { + fh.open(file_out, std::ios::out); + } + if (fh.is_open()) { + //int n_digits = std::numeric_limits::digits10; + int n_digits = 6; + fh.precision(n_digits); + fh.setf(std::ios::scientific); + for (std::size_t i=0; i < m.n_rows(); ++i) { + for (std::size_t j=0; j < m.n_cols(); ++j) { + fh.width(n_digits+8); // 8 extra chars for spacing, +/- and scientific formatting + fh << m(i, j); + } fh << std::endl; } + } else { + throw FileNotFoundError(); } - if (in_line) { - fh << std::endl; + } + } // end namespace UseCSV + + namespace UseXTC { + RvecArray::RvecArray(std::size_t n_atoms) + : n_atoms(n_atoms) { + if (this->n_atoms > 0) { + this->values = static_cast(calloc(this->n_atoms, sizeof(this->values[0]))); } - } else { - throw FileNotFoundError(); } - } - - template - void write_vector(std::string file_out, std::vector v, bool append) { - std::ofstream fh; - if (append) { - fh.open(file_out, std::ios::out | std::ios::app); - } else { - fh.open(file_out, std::ios::out); + + RvecArray::~RvecArray() { + free(this->values); } - write_vector(fh, v); - } - - template - void write_matrix(std::string file_out, Matrix m, bool append) { - std::ofstream fh; - if (append) { - fh.open(file_out, std::ios::out | std::ios::app); - } else { - fh.open(file_out, std::ios::out); + + std::size_t n_atoms(std::string file_in) { + int n_atoms=0; + read_xtc_natoms(file_in.c_str(), &n_atoms); + return (std::size_t) n_atoms; } - if (fh.is_open()) { - //int n_digits = std::numeric_limits::digits10; - int n_digits = 6; - fh.precision(n_digits); - fh.setf(std::ios::scientific); - for (std::size_t i=0; i < m.n_rows(); ++i) { - for (std::size_t j=0; j < m.n_cols(); ++j) { - fh.width(n_digits+8); // 8 extra chars for spacing, +/- and scientific formatting - fh << m(i, j); + + template + Matrix read_chunk(std::shared_ptr file, std::size_t n_cols, std::size_t n_lines) { + int n_fake_atoms = n_cols/3; + int return_code = 0; + int n_step = 0; + float time_step = 1.0; + matrix box_vec; + RvecArray coords(n_fake_atoms); + float prec = 1000.0; + Matrix m(n_lines, n_cols); + + std::size_t n_lines_0 = n_lines; + for (; n_lines > 0; --n_lines) { + // read next frame + return_code = read_xtc(file.get(), n_fake_atoms, &n_step, + &time_step, box_vec, coords.values, &prec); + if (return_code == exdrOK) { + // add next frame to input matrix + std::vector r(n_cols); + for (int i = 0; i < n_fake_atoms; ++i) { + for (int j = 0; j < 3; ++j) { + m(n_lines_0-n_lines, i*3+j) = (T) coords.values[i][j]; + } + } + } else if (return_code == exdrENDOFFILE) { + //EOF: if no line has been read, just return an empty matrix. + // (i.e. return 'm' as initialized) + break; + } else if (return_code == exdrFILENOTFOUND) { + throw std::runtime_error("xtc file not found"); + break; + } else { + throw std::runtime_error("strange error while reading xtc file"); + break; } - fh << std::endl; } - } else { - throw FileNotFoundError(); + if (n_lines != 0) { + // resize m if needed + m = std::move(m.limited_rows(n_lines_0-n_lines)); + } + return m; } - } -} // end namespace UseCSV - -namespace UseXTC { - RvecArray::RvecArray(std::size_t n_atoms) - : n_atoms(n_atoms) { - if (this->n_atoms > 0) { - this->values = static_cast(calloc(this->n_atoms, sizeof(this->values[0]))); + + template + void write_vector(std::string file_out, std::vector v, bool append) { + // write vector into column, keeping two additional columns empty. + // this is needed since xtc files only can store cartesian coordinates, + // i.e. packets of three numbers. + // vectors will be written in first 'x' field with values in different + // rows (aka time-frames) to get it column-oriented as first column, + // if xtc-data is dumped (e.g. with gmxdump). + Matrix m(v.size(), 3); + for (std::size_t i=0; i < v.size(); ++i) { + m(i,0) = v[i]; + m(i,1) = 0.0; + m(i,2) = 0.0; + } + write_matrix(file_out, m, append); } - } - - RvecArray::~RvecArray() { - free(this->values); - } - - std::size_t n_atoms(std::string file_in) { - int n_atoms=0; - read_xtc_natoms(file_in.c_str(), &n_atoms); - return (std::size_t) n_atoms; - } - - template - Matrix read_chunk(std::shared_ptr file, std::size_t n_cols, std::size_t n_lines) { - int n_fake_atoms = n_cols/3; - int return_code = 0; - int n_step = 0; - float time_step = 1.0; - matrix box_vec; - RvecArray coords(n_fake_atoms); - float prec = 1000.0; - Matrix m(n_lines, n_cols); - - std::size_t n_lines_0 = n_lines; - for (; n_lines > 0; --n_lines) { - // read next frame - return_code = read_xtc(file.get(), n_fake_atoms, &n_step, - &time_step, box_vec, coords.values, &prec); - if (return_code == exdrOK) { - // add next frame to input matrix - std::vector r(n_cols); - for (int i = 0; i < n_fake_atoms; ++i) { - for (int j = 0; j < 3; ++j) { - m(n_lines_0-n_lines, i*3+j) = (T) coords.values[i][j]; + + template + void write_matrix(std::string file_out, Matrix m, bool append) { + std::shared_ptr file(xdrfile_open(file_out.c_str(), append ? "a" : "w"), + // deleter function called when pointer goes out of scope + [=](XDRFILE* f) { + xdrfile_close(f); + }); + // faked box matrix for xtc-file + float fake_box_matrix[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}; + // xtc-files are only readable/writable in 3D cartesian coordinates + // for atoms. if the number of matrix columns is not divisable by three, + // add zero-valued columns to fill. + int n_fake_atoms = m.n_cols() / 3; + if (m.n_cols() % 3 == 0) { + ++n_fake_atoms; + } + std::size_t n_fake_cols = 3*n_fake_atoms; + std::size_t n_cols = m.n_cols(); + RvecArray buf(3*n_fake_atoms); + for (std::size_t r=0; r < m.n_rows(); ++r) { + for (std::size_t c=0; c < n_fake_cols; ++c) { + if (c < n_cols) { + buf.values[c/3][c%3] = m(r,c); + } else { + buf.values[c/3][c%3] = 0.0; } } - } else if (return_code == exdrENDOFFILE) { - //EOF: if no line has been read, just return an empty matrix. - // (i.e. return 'm' as initialized) - break; - } else if (return_code == exdrFILENOTFOUND) { - throw std::runtime_error("xtc file not found"); - break; - } else { - throw std::runtime_error("strange error while reading xtc file"); - break; + write_xtc(file.get(), n_fake_atoms, r, 1.0, fake_box_matrix, buf.values, 1000.0); } } - if (n_lines != 0) { - // resize m if needed - m = std::move(m.limited_rows(n_lines_0-n_lines)); + } // end UseXTC namespace + } // end local namespace + + + template + DataFileReader::DataFileReader(std::string filename, std::size_t buf_size_bytes) + : _filename(filename), + _buf_size_bytes(buf_size_bytes), + _n_cols(0), + _eof(false) { + this->_ftype = filename_suffix(filename); + + if (this->_ftype == XTC) { + // use XTC file + //this->_fh_xtc = UseXTC::XTCFile(filename, "r"); + this->_fh_xtc = std::shared_ptr(xdrfile_open(filename.c_str(), "r"), + // deleter function called when pointer goes out of scope + [](XDRFILE* f) { + xdrfile_close(f); + }); + this->_n_cols = 3*UseXTC::n_atoms(filename); + } else { + this->_parser.open(filename); + this->_n_cols = this->_parser.n_cols(); } - return m; } - + template - void write_vector(std::string file_out, std::vector v, bool append) { - // write vector into column, keeping two additional columns empty. - // this is needed since xtc files only can store cartesian coordinates, - // i.e. packets of three numbers. - // vectors will be written in first 'x' field with values in different - // rows (aka time-frames) to get it column-oriented as first column, - // if xtc-data is dumped (e.g. with gmxdump). - Matrix m(v.size(), 3); - for (std::size_t i=0; i < v.size(); ++i) { - m(i,0) = v[i]; - m(i,1) = 0.0; - m(i,2) = 0.0; - } - write_matrix(file_out, m, append); + DataFileReader::DataFileReader(std::string filename) + : DataFileReader(filename, FastPCA::gigabytes_to_bytes(1)) { } template - void write_matrix(std::string file_out, Matrix m, bool append) { - std::shared_ptr file(xdrfile_open(file_out.c_str(), append ? "a" : "w"), - // deleter function called when pointer goes out of scope - [=](XDRFILE* f) { - xdrfile_close(f); - }); - // faked box matrix for xtc-file - float fake_box_matrix[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}}; - // xtc-files are only readable/writable in 3D cartesian coordinates - // for atoms. if the number of matrix columns is not divisable by three, - // add zero-valued columns to fill. - int n_fake_atoms = m.n_cols() / 3; - if (m.n_cols() % 3 == 0) { - ++n_fake_atoms; + Matrix DataFileReader::next_block(std::size_t n_rows) { + std::size_t n_possible_rows; + if (n_rows == 0){ + n_possible_rows = this->_buf_size_bytes / (sizeof(T) * this->_n_cols); + } else { + n_possible_rows = n_rows; } - std::size_t n_fake_cols = 3*n_fake_atoms; - std::size_t n_cols = m.n_cols(); - RvecArray buf(3*n_fake_atoms); - for (std::size_t r=0; r < m.n_rows(); ++r) { - for (std::size_t c=0; c < n_fake_cols; ++c) { - if (c < n_cols) { - buf.values[c/3][c%3] = m(r,c); - } else { - buf.values[c/3][c%3] = 0.0; + if (this->_ftype == XTC) { + if (this->_n_cols == 0) { + throw FileFormatError(); + } + Matrix m = UseXTC::read_chunk(this->_fh_xtc, this->_n_cols, n_possible_rows); + if (m.n_rows() == 0) { + this->_eof = true; + return Matrix(); + } else { + return m; + } + } else { + // assume ASCII file + if (this->_parser.eof()) { + this->_eof = true; + return Matrix(); + } else { + if (this->_n_cols == 0) { + throw FileFormatError(); } + //return UseCSV::read_chunk(this->_fh, this->_n_cols, n_possible_rows); + return UseCSV::read_chunk(this->_parser, n_possible_rows); } - write_xtc(file.get(), n_fake_atoms, r, 1.0, fake_box_matrix, buf.values, 1000.0); } } -} // end UseXTC namespace -} // end local namespace - - -template -DataFileReader::DataFileReader(std::string filename, std::size_t buf_size_bytes) - : _filename(filename), - _buf_size_bytes(buf_size_bytes), - _n_cols(0), - _eof(false) { - this->_ftype = filename_suffix(filename); - - if (this->_ftype == XTC) { - // use XTC file - //this->_fh_xtc = UseXTC::XTCFile(filename, "r"); - this->_fh_xtc = std::shared_ptr(xdrfile_open(filename.c_str(), "r"), - // deleter function called when pointer goes out of scope - [](XDRFILE* f) { - xdrfile_close(f); - }); - this->_n_cols = 3*UseXTC::n_atoms(filename); - } else { - this->_parser.open(filename); - this->_n_cols = this->_parser.n_cols(); + + template + std::size_t DataFileReader::n_cols() { + return this->_n_cols; } -} - -template -DataFileReader::DataFileReader(std::string filename) - : DataFileReader(filename, gigabytes_to_bytes(1)) { -} - -template -Matrix DataFileReader::next_block(std::size_t n_rows) { - std::size_t n_possible_rows; - if (n_rows == 0){ - n_possible_rows = this->_buf_size_bytes / (sizeof(T) * this->_n_cols); - } else { - n_possible_rows = n_rows; + + template + bool DataFileReader::eof() const { + return this->_eof; } - if (this->_ftype == XTC) { - if (this->_n_cols == 0) { - throw FileFormatError(); - } - Matrix m = UseXTC::read_chunk(this->_fh_xtc, this->_n_cols, n_possible_rows); - if (m.n_rows() == 0) { - this->_eof = true; - return Matrix(); - } else { - return m; - } - } else { - // assume ASCII file - if (this->_parser.eof()) { - this->_eof = true; - return Matrix(); - } else { - if (this->_n_cols == 0) { - throw FileFormatError(); - } - //return UseCSV::read_chunk(this->_fh, this->_n_cols, n_possible_rows); - return UseCSV::read_chunk(this->_parser, n_possible_rows); - } + + //// + + template + DataFileWriter::DataFileWriter(std::string filename) + : _filename(filename) { + this->_ftype = filename_suffix(filename); } -} - -template -std::size_t DataFileReader::n_cols() { - return this->_n_cols; -} - -template -bool DataFileReader::eof() const { - return this->_eof; -} - -//// - -template -DataFileWriter::DataFileWriter(std::string filename) - : _filename(filename) { - this->_ftype = filename_suffix(filename); -} - -template -void DataFileWriter::write(std::vector v, bool append) { - switch (this->_ftype) { - case XTC: - UseXTC::write_vector(this->_filename, v, append); - break; - case UNKNOWN: - UseCSV::write_vector(this->_filename, v, append); - break; - default: - throw std::runtime_error("unknown file type to write vector"); + + template + void DataFileWriter::write(std::vector v, bool append) { + switch (this->_ftype) { + case XTC: + UseXTC::write_vector(this->_filename, v, append); + break; + case UNKNOWN: + UseCSV::write_vector(this->_filename, v, append); + break; + default: + throw std::runtime_error("unknown file type to write vector"); + } } -} - -template -void DataFileWriter::write(Matrix m, bool append) { - switch (this->_ftype) { - case XTC: - UseXTC::write_matrix(this->_filename, m, append); - break; - case UNKNOWN: - UseCSV::write_matrix(this->_filename, m, append); - break; - default: - throw std::runtime_error("unknown file type given to write matrix"); + + template + void DataFileWriter::write(Matrix m, bool append) { + switch (this->_ftype) { + case XTC: + UseXTC::write_matrix(this->_filename, m, append); + break; + case UNKNOWN: + UseCSV::write_matrix(this->_filename, m, append); + break; + default: + throw std::runtime_error("unknown file type given to write matrix"); + } } -} } // end namespace FastPCA diff --git a/src/util.hpp b/src/util.hpp index 882d2c7..608c1e3 100644 --- a/src/util.hpp +++ b/src/util.hpp @@ -26,6 +26,16 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "matrix.hpp" + +// resolve circular dependency +// (used in file_io) +namespace FastPCA { + constexpr std::size_t + gigabytes_to_bytes(std::size_t gb) { + return gb * 1073741824; + } +} + #include "file_io.hpp" #include @@ -33,11 +43,6 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace FastPCA { - constexpr std::size_t - gigabytes_to_bytes(std::size_t gb) { - return gb * 1073741824; - } - bool is_comment_or_empty(std::string line); @@ -95,10 +100,16 @@ namespace FastPCA { * periodic interval inside boundaries of * [-periodicity/2, +periodicity/2]. */ - tamplate + template constexpr T normalized(T var, T periodicity); + template + void + shift_matrix_columns_inplace(Matrix& m + , std::vector shifts + , std::vector periodicities); + template void shift_matrix_columns_inplace(Matrix& m, std::vector shifts); diff --git a/src/util.hxx b/src/util.hxx index 8fea00f..ce2c71d 100644 --- a/src/util.hxx +++ b/src/util.hxx @@ -130,33 +130,42 @@ namespace FastPCA { } namespace Periodic { - tamplate + template constexpr T normalized(T var, T periodicity) { - return fmod(var + periodicity/2.0) - periodicity/2.0; + return fmod(var+0.5*periodicity, periodicity) - 0.5*periodicity; } template void shift_matrix_columns_inplace(Matrix& m - , std::vector shifts) { + , std::vector shifts + , std::vector periodicities) { std::size_t i,j; const std::size_t n_rows = m.n_rows(); const std::size_t n_cols = m.n_cols(); #pragma omp parallel for default(none)\ private(i,j)\ firstprivate(n_rows,n_cols)\ - shared(m,shifts) + shared(m,shifts, periodicities) for (j=0; j < n_cols; ++j) { for (i=0; i < n_rows; ++i) { m(i,j) = m(i,j) - shifts[j]; // periodic boundary checks // TODO: test and remove old code //m(i,j) = atan2(sin(m(i,j)), cos(m(i,j))); - m(i,j) = normalized(m(i,j), 2*M_PI); + m(i,j) = normalized(m(i,j), periodicities[j]); } } } + + template + void + shift_matrix_columns_inplace(Matrix& m + , std::vector shifts) { + std::vector p(shifts.size(), 2*M_PI); + shift_matrix_columns_inplace(m, shifts, p); + } } // FastPCA::Periodic } // namespace FastPCA