From 96c4f65c555af8e4c6a66298ba8c947ea5251c08 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Sun, 24 Oct 2021 22:11:33 +0200 Subject: [PATCH 01/11] neue Matrixstruktur in MRM Scoring --- .../OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h | 27 +- .../OpenMS/OPENSWATHALGO/ALGO/Scoring.h | 4 +- .../OPENSWATHALGO/DATASTRUCTURES/Matrix.cpp | 41 ++ .../OPENSWATHALGO/DATASTRUCTURES/Matrix.h | 403 ++++++++++++++++++ .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 372 ++++++++-------- .../source/OPENSWATHALGO/ALGO/Scoring.cpp | 27 +- .../openswathalgo/MRMScoring_test.cpp | 116 +++-- .../openswathalgo/Scoring_test.cpp | 18 +- 8 files changed, 716 insertions(+), 292 deletions(-) create mode 100644 src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.cpp create mode 100644 src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.h diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h index 23f941af2e5..37f2ac528f5 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h @@ -38,6 +38,7 @@ #include // for isnan #include +#include #include #include "OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h" @@ -80,7 +81,7 @@ namespace OpenSwath /// Cross Correlation array typedef OpenSwath::Scoring::XCorrArrayType XCorrArrayType; /// Cross Correlation matrix - typedef std::vector > XCorrMatrixType; + typedef OpenMS::Matrix XCorrMatrixType; typedef std::string String; @@ -113,7 +114,7 @@ namespace OpenSwath /** @name Scores */ //@{ /// Initialize the scoring object and building the cross-correlation matrix - void initializeXCorrMatrix(const std::vector< std::vector< double > >& data); + void initializeXCorrMatrix(const std::vector< std::vector< double > >& data); //andere matrix? /// Initialize the scoring object and building the cross-correlation matrix void initializeXCorrMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids); @@ -220,23 +221,23 @@ namespace OpenSwath std::vector& signal_noise_estimators); /// non-mutable access to the MI matrix - const std::vector< std::vector > & getMIMatrix() const; + const OpenMS::Matrix & getMIMatrix() const; //@} /// non-mutable access to the MI contrast matrix - const std::vector< std::vector > & getMIContrastMatrix() const; + const OpenMS::Matrix & getMIContrastMatrix() const; //@} /// non-mutable access to the MI precursor contrast matrix - const std::vector< std::vector > & getMIPrecursorContrastMatrix() const; + const OpenMS::Matrix & getMIPrecursorContrastMatrix() const; //@} /// non-mutable access to the MI precursor combined matrix - const std::vector< std::vector > & getMIPrecursorCombinedMatrix() const; + const OpenMS::Matrix & getMIPrecursorCombinedMatrix() const; //@} /// Initialize the scoring object and building the MI matrix - void initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids); + void initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids); //ref? /// Initialize the scoring object and building the MI matrix of chromatograms of set1 (e.g. identification transitions) vs set2 (e.g. detection transitions) void initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids_set1, std::vector native_ids_set2); @@ -282,20 +283,20 @@ namespace OpenSwath //@} /// the precomputed mutual information matrix - std::vector< std::vector > mi_matrix_; - + //std::vector< std::vector > mi_matrix_; + OpenMS::Matrix mi_matrix_; /// the precomputed contrast mutual information matrix - std::vector< std::vector > mi_contrast_matrix_; + OpenMS::Matrix mi_contrast_matrix_; /// the precomputed mutual information matrix of the MS1 trace - std::vector< std::vector > mi_precursor_matrix_; + OpenMS::Matrix mi_precursor_matrix_; /// the precomputed contrast mutual information matrix against the MS1 trace - std::vector< std::vector > mi_precursor_contrast_matrix_; + OpenMS::Matrix mi_precursor_contrast_matrix_; //@} /// the precomputed contrast mutual information matrix with the MS1 trace - std::vector< std::vector > mi_precursor_combined_matrix_; + OpenMS::Matrix mi_precursor_combined_matrix_; //@} }; diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h index 1b03beba28d..4abfc7aa123 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h @@ -117,11 +117,11 @@ namespace OpenSwath /// Calculate crosscorrelation on std::vector data (which is first normalized) /// NOTE: this replaces calcxcorr OPENSWATHALGO_DLLAPI XCorrArrayType normalizedCrossCorrelation(std::vector& data1, - std::vector& data2, const int& maxdelay, const int& lag); + std::vector& data2, int maxdelay, int lag); /// Calculate crosscorrelation on std::vector data without normalization OPENSWATHALGO_DLLAPI XCorrArrayType calculateCrossCorrelation(const std::vector& data1, - const std::vector& data2, const int& maxdelay, const int& lag); + const std::vector& data2, int maxdelay, int lag); /// Find best peak in an cross-correlation (highest apex) OPENSWATHALGO_DLLAPI XCorrArrayType::const_iterator xcorrArrayGetMaxPeak(const XCorrArrayType & array); diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.cpp b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.cpp new file mode 100644 index 00000000000..7c94e0b1319 --- /dev/null +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.cpp @@ -0,0 +1,41 @@ +// -------------------------------------------------------------------------- +// OpenMS -- Open-Source Mass Spectrometry +// -------------------------------------------------------------------------- +// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen, +// ETH Zurich, and Freie Universitaet Berlin 2002-2021. +// +// This software is released under a three-clause BSD license: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// * Neither the name of any author or any participating institution +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// For a full list of authors, refer to the file AUTHORS. +// -------------------------------------------------------------------------- +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING +// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// -------------------------------------------------------------------------- +// $Maintainer: Timo Sachsenberg $ +// $Authors: $ +// -------------------------------------------------------------------------- + +#include + +namespace OpenMS +{ + Matrix default_matrix_int; + Matrix default_matrix_double; +} diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.h new file mode 100644 index 00000000000..619f7493502 --- /dev/null +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/DATASTRUCTURES/Matrix.h @@ -0,0 +1,403 @@ +// -------------------------------------------------------------------------- +// OpenMS -- Open-Source Mass Spectrometry +// -------------------------------------------------------------------------- +// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen, +// ETH Zurich, and Freie Universitaet Berlin 2002-2021. +// +// This software is released under a three-clause BSD license: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// * Neither the name of any author or any participating institution +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// For a full list of authors, refer to the file AUTHORS. +// -------------------------------------------------------------------------- +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING +// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// -------------------------------------------------------------------------- +// $Maintainer: Timo Sachsenberg $ +// $Authors: $ +// -------------------------------------------------------------------------- + +#pragma once + + +#include // pow() +#include +#include +#include + +namespace OpenMS +{ + + /** + @brief A two-dimensional matrix. Similar to std::vector, but uses a binary + operator(,) for element access. + + Think of it as a random access container. You can also generate gray + scale images. This data structure is not designed to be used for linear algebra, + but rather a simple two-dimensional array. + + The following member functions of the base class std::vector + can also be used: + +
    +
  • begin
    • +
  • end
    • +
  • rbegin
    • +
  • rend
    • +
  • front
    • +
  • back
    • +
  • assign
    • +
  • empty
    • +
  • size
    • +
  • capacity
    • +
  • max_size
    • +
+ + @ingroup Datastructures + */ + template + class Matrix : + protected std::vector + { + protected: + typedef std::vector Base; + + public: + + ///@name STL compliance type definitions + //@{ + typedef Base container_type; + + typedef typename Base::difference_type difference_type; + typedef typename Base::size_type size_type; + + typedef typename Base::const_iterator const_iterator; + typedef typename Base::const_reverse_iterator const_reverse_iterator; + typedef typename Base::iterator iterator; + typedef typename Base::reverse_iterator reverse_iterator; + + typedef typename Base::const_reference const_reference; + typedef typename Base::pointer pointer; + typedef typename Base::reference reference; + typedef typename Base::value_type value_type; + + typedef typename Base::allocator_type allocator_type; + //@} + + ///@name OpenMS compliance type definitions + //@{ + typedef Base ContainerType; + typedef difference_type DifferenceType; + typedef size_type SizeType; + + typedef const_iterator ConstIterator; + typedef const_reverse_iterator ConstReverseIterator; + typedef iterator Iterator; + typedef reverse_iterator ReverseIterator; + + typedef const_reference ConstReference; + typedef pointer Pointer; + typedef reference Reference; + typedef value_type ValueType; + + typedef allocator_type AllocatorType; + //@} + + ///@name Constructors, assignment, and destructor + //@{ + Matrix() : + Base(), + rows_(0), + cols_(0) + {} + + Matrix(const SizeType rows, const SizeType cols, ValueType value = ValueType()) : + Base(rows * cols, value), + rows_(rows), + cols_(cols) + {} + + Matrix(const Matrix& source) : + Base(source), + rows_(source.rows_), + cols_(source.cols_) + {} + + Matrix& operator=(const Matrix& rhs) + { + Base::operator=(rhs); + rows_ = rhs.rows_; + cols_ = rhs.cols_; + return *this; + } + + ~Matrix() {} + //@} + + ///@name Accessors + //@{ + const_reference operator()(size_type const i, size_type const j) const + { + return getValue(i, j); + } + + reference operator()(size_type const i, size_type const j) + { + return getValue(i, j); + } + + const_reference getValue(size_type const i, size_type const j) const + { + return Base::operator[](index(i, j)); + } + + reference getValue(size_type const i, size_type const j) + { + return Base::operator[](index(i, j)); + } + + void setValue(size_type const i, size_type const j, value_type value) + { + Base::operator[](index(i, j)) = value; + } + + /// Return the i-th row of the matrix as a vector. + container_type row(size_type const i) const + { +#ifdef OPENMS_DEBUG + if (i >= rows_) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, i, rows_); +#endif + container_type values(cols_); + for (size_type j = 0; j < cols_; j++) + { + values[j] = Base::operator[](index(i, j)); + } + return values; + } + + /// Return the i-th column of the matrix as a vector. + container_type col(size_type const i) const + { +#ifdef OPENMS_DEBUG + if (i >= cols_) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, i, cols_); +#endif + container_type values(rows_); + for (size_type j = 0; j < rows_; j++) + { + values[j] = Base::operator[](index(j, i)); + } + return values; + } + + //@} + + + /** + @name Pure access declarations + + These make begin(), end() etc. from container_type accessible. + */ + //@{ + public: + + using Base::begin; + using Base::end; + using Base::rbegin; + using Base::rend; + + using Base::front; + using Base::back; + using Base::assign; + + using Base::empty; + using Base::size; + + using Base::capacity; + using Base::max_size; + + //@} + + void clear() + { + Base::clear(); + rows_ = 0; + cols_ = 0; + } + + void resize(size_type i, size_type j, value_type value = value_type()) + { + rows_ = i; + cols_ = j; + Base::resize(rows_ * cols_, value); + } + + void resize(std::pair const& size_pair, value_type value = value_type()) + { + rows_ = size_pair.first; + cols_ = size_pair.second; + Base::resize(rows_ * cols_, value); + } + + /// Number of rows + SizeType rows() const + { + return rows_; + } + + /// Number of columns + SizeType cols() const + { + return cols_; + } + + std::pair sizePair() const + { + return std::pair(rows_, cols_); + } + + /** + @brief Calculate the index into the underlying vector from row and column. + Note that Matrix uses the (row,column) lexicographic ordering for indexing. + */ + SizeType const index(SizeType row, SizeType col) const + { +#ifdef OPENMS_DEBUG + if (row >= rows_) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, row, rows_); + if (col >= cols_) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, col, cols_); +#endif + return row * cols_ + col; + } + + /** + @brief Calculate the row and column from an index into the underlying vector. + Note that Matrix uses the (row,column) lexicographic ordering for indexing. + */ + std::pair const indexPair(std::size_t index) const + { +#ifdef OPENMS_DEBUG + if (index >= size()) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, index, size() - 1); +#endif + return std::pair(index / cols_, index % cols_); + } + + /** + @brief Calculate the column from an index into the underlying vector. + Note that Matrix uses the (row,column) lexicographic ordering for indexing. + */ + SizeType colIndex(SizeType index) const + { +#ifdef OPENMS_DEBUG + if (index >= size()) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, index, size() - 1); +#endif + return index % cols_; + } + + /** + @brief Calculate the row from an index into the underlying vector. + Note that Matrix uses the (row,column) lexicographic ordering for indexing. + */ + SizeType rowIndex(SizeType index) const + { +#ifdef OPENMS_DEBUG + if (index >= size()) throw Exception::IndexOverflow(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, index, size() - 1); +#endif + + return index / cols_; + } + + /** + @brief Equality comparator. + + If matrices have different row or column numbers, throws a precondition exception. + */ + bool operator==(Matrix const& rhs) const + { + OPENMS_PRECONDITION(cols_ == rhs.cols_, + "Matrices have different row sizes."); + OPENMS_PRECONDITION(rows_ == rhs.rows_, + "Matrices have different column sizes."); + return static_cast::Base const&>(*this) == static_cast::Base const&>(rhs); + } + + /** + @brief Less-than comparator. Comparison is done lexicographically: first by row, then by column. + + If matrices have different row or column numbers, throws a precondition exception. + */ + bool operator<(Matrix const& rhs) const + { + OPENMS_PRECONDITION(cols_ == rhs.cols_, + "Matrices have different row sizes."); + OPENMS_PRECONDITION(rows_ == rhs.rows_, + "Matrices have different column sizes."); + return static_cast::Base const&>(*this) < static_cast::Base const&>(rhs); + } + + /// set matrix to 2D arrays values + template + void setMatrix(const ValueType matrix[ROWS][COLS]) + { + resize(ROWS, COLS); + for (SizeType i = 0; i < this->rows_; ++i) + { + for (SizeType j = 0; j < this->cols_; ++j) + { + setValue(i, j, matrix[i][j]); + } + } + } + + const Base& asVector() + { + return *this; + } + + protected: + + ///@name Data members + //@{ + /// Number of rows (height of a column) + SizeType rows_; + /// Number of columns (width of a row) + SizeType cols_; + //@} + + }; // class Matrix + + /** + @brief Print the contents to a stream. + + @relatesalso Matrix + */ + template + std::ostream& operator<<(std::ostream& os, const Matrix& matrix) + { + typedef typename Matrix::size_type size_type; + for (size_type i = 0; i < matrix.rows(); ++i) + { + for (size_type j = 0; j < matrix.cols(); ++j) + { + os << std::setprecision(6) << std::setw(6) << matrix(i, j) << ' '; + } + os << std::endl; + } + return os; + } + +} // namespace OpenMS + diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 6b9ee6278cd..848daf4d8a6 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -34,6 +34,7 @@ #include #include +#include #include //#define MRMSCORING_TESTING #include @@ -51,16 +52,15 @@ namespace OpenSwath void MRMScoring::initializeXCorrMatrix(const std::vector< std::vector< double > >& data) { - xcorr_matrix_.resize(data.size()); + xcorr_matrix_.resize(data.size(), data.size()); for (std::size_t i = 0; i < data.size(); i++) { - xcorr_matrix_[i].resize(data.size()); + std::vector< double > tmp1(data[i]); for (std::size_t j = i; j < data.size(); j++) { // compute normalized cross correlation - std::vector< double > tmp1(data[i]); std::vector< double > tmp2(data[j]); - xcorr_matrix_[i][j] = Scoring::normalizedCrossCorrelation(tmp1, tmp2, boost::numeric_cast(data[i].size()), 1); + xcorr_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(tmp1, tmp2, boost::numeric_cast(data[i].size()), 1)); } } } @@ -83,12 +83,11 @@ namespace OpenSwath void MRMScoring::initializeXCorrMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids) { std::vector intensityi, intensityj; - xcorr_matrix_.resize(native_ids.size()); + xcorr_matrix_.resize(native_ids.size(), native_ids.size()); for (std::size_t i = 0; i < native_ids.size(); i++) { String native_id = native_ids[i]; FeatureType fi = mrmfeature->getFeature(native_id); - xcorr_matrix_[i].resize(native_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < native_ids.size(); j++) @@ -98,7 +97,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation - xcorr_matrix_[i][j] = Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1); + xcorr_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); } } } @@ -106,12 +105,11 @@ namespace OpenSwath void MRMScoring::initializeXCorrContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids_set1, const std::vector& native_ids_set2) { std::vector intensityi, intensityj; - xcorr_contrast_matrix_.resize(native_ids_set1.size()); + xcorr_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); for (std::size_t i = 0; i < native_ids_set1.size(); i++) { String native_id = native_ids_set1[i]; FeatureType fi = mrmfeature->getFeature(native_id); - xcorr_contrast_matrix_[i].resize(native_ids_set2.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids_set2.size(); j++) @@ -121,7 +119,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation - xcorr_contrast_matrix_[i][j] = Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1); + xcorr_contrast_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); } } } @@ -129,12 +127,11 @@ namespace OpenSwath void MRMScoring::initializeXCorrPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids) { std::vector intensityi, intensityj; - xcorr_precursor_matrix_.resize(precursor_ids.size()); + xcorr_precursor_matrix_.resize(precursor_ids.size(), precursor_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { String precursor_id = precursor_ids[i]; FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - xcorr_precursor_matrix_[i].resize(precursor_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < precursor_ids.size(); j++) @@ -144,7 +141,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation - xcorr_precursor_matrix_[i][j] = Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1); + xcorr_precursor_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); } } } @@ -152,12 +149,11 @@ namespace OpenSwath void MRMScoring::initializeXCorrPrecursorContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids, const std::vector& native_ids) { std::vector intensityi, intensityj; - xcorr_precursor_contrast_matrix_.resize(precursor_ids.size()); + xcorr_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { String precursor_id = precursor_ids[i]; FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - xcorr_precursor_contrast_matrix_[i].resize(native_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids.size(); j++) @@ -167,23 +163,22 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation - xcorr_precursor_contrast_matrix_[i][j] = Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1); + xcorr_precursor_contrast_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); } } } void MRMScoring::initializeXCorrPrecursorContrastMatrix(const std::vector< std::vector< double > >& data_precursor, const std::vector< std::vector< double > >& data_fragments) { - xcorr_precursor_contrast_matrix_.resize(data_precursor.size()); + xcorr_precursor_contrast_matrix_.resize(data_precursor.size(), data_fragments.size()); for (std::size_t i = 0; i < data_precursor.size(); i++) - { - xcorr_precursor_contrast_matrix_[i].resize(data_fragments.size()); + { + std::vector< double > tmp1(data_precursor[i]); for (std::size_t j = 0; j < data_fragments.size(); j++) { // compute normalized cross correlation - std::vector< double > tmp1(data_precursor[i]); std::vector< double > tmp2(data_fragments[j]); - xcorr_precursor_contrast_matrix_[i][j] = Scoring::normalizedCrossCorrelation(tmp1, tmp2, boost::numeric_cast(tmp1.size()), 1); + xcorr_precursor_contrast_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(tmp1, tmp2, boost::numeric_cast(tmp1.size()), 1)); #ifdef MRMSCORING_TESTING std::cout << " fill xcorr_precursor_contrast_matrix_ "<< tmp1.size() << " / " << tmp2.size() << " : " << xcorr_precursor_contrast_matrix_[i][j].data.size() << std::endl; #endif @@ -209,11 +204,11 @@ namespace OpenSwath features.push_back(fj); } - xcorr_precursor_combined_matrix_.resize(features.size()); + xcorr_precursor_combined_matrix_.resize(features.size(), features.size()); for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; - xcorr_precursor_combined_matrix_[i].resize(features.size()); + //xcorr_precursor_combined_matrix_[i].resize(features.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < features.size(); j++) @@ -222,7 +217,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation - xcorr_precursor_combined_matrix_[i][j] = Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1); + xcorr_precursor_combined_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); } } } @@ -235,15 +230,15 @@ namespace OpenSwath // return $deltascore_mean + $deltascore_stdev double MRMScoring::calcXcorrCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_matrix_.size(); i++) + for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) { - for (std::size_t j = i; j < xcorr_matrix_.size(); j++) + for (std::size_t j = i; j < xcorr_matrix_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first)); + deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->first)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first) << std::endl; #endif @@ -262,16 +257,15 @@ namespace OpenSwath double MRMScoring::calcXcorrCoelutionWeightedScore( const std::vector& normalized_library_intensity) { - OPENSWATH_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); #ifdef MRMSCORING_TESTING double weights = 0; #endif - std::vector deltas; - for (std::size_t i = 0; i < xcorr_matrix_.size(); i++) + double deltas{0}; + for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) { - deltas.push_back( - std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->first) + deltas += (std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, i))->first) * normalized_library_intensity[i] * normalized_library_intensity[i]); #ifdef MRMSCORING_TESTING @@ -279,11 +273,10 @@ namespace OpenSwath normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl; weights += normalized_library_intensity[i] * normalized_library_intensity[i]; #endif - for (std::size_t j = i + 1; j < xcorr_matrix_.size(); j++) + for (std::size_t j = i + 1; j < xcorr_matrix_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back( - std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first) + deltas += (std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->first) * normalized_library_intensity[i] * normalized_library_intensity[j] * 2); #ifdef MRMSCORING_TESTING @@ -299,24 +292,21 @@ namespace OpenSwath std::cout << " all weights sum " << weights << std::endl; #endif - return std::accumulate(deltas.begin(), deltas.end(), 0.0); + return deltas; } double MRMScoring::calcXcorrContrastCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.size() > 0 && xcorr_contrast_matrix_[0].size() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_contrast_matrix_.size(); i++) + for (auto e : xcorr_contrast_matrix_) { - for (std::size_t j = 0; j < xcorr_contrast_matrix_[0].size(); j++) - { - // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first)); + // first is the X value (RT), should be an int + deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(e)->first)); #ifdef MRMSCORING_TESTING - std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first) << std::endl; + std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first) << std::endl; #endif - } } OpenSwath::mean_and_stddev msc; @@ -330,21 +320,21 @@ namespace OpenSwath std::vector MRMScoring::calcSeparateXcorrContrastCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.size() > 0 && xcorr_contrast_matrix_[0].size() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_contrast_matrix_.size(); i++) + for (std::size_t i = 0; i < xcorr_contrast_matrix_.rows(); i++) { double deltas_id = 0; - for (std::size_t j = 0; j < xcorr_contrast_matrix_[0].size(); j++) + for (std::size_t j = 0; j < xcorr_contrast_matrix_.cols(); j++) { // first is the X value (RT), should be an int - deltas_id += std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first); + deltas_id += std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_.getValue(i,j))->first); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first) << std::endl; #endif } - deltas.push_back(deltas_id / xcorr_contrast_matrix_[0].size()); + deltas.push_back(deltas_id / xcorr_contrast_matrix_.cols()); } return deltas; @@ -352,15 +342,15 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_precursor_matrix_.size(); i++) + for (std::size_t i = 0; i < xcorr_precursor_matrix_.rows(); i++) { - for (std::size_t j = i; j < xcorr_precursor_matrix_.size(); j++) + for (std::size_t j = i; j < xcorr_precursor_matrix_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_[i][j])->first)); + deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_.getValue(i, j))->first)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_[i][j])->first) << std::endl; #endif @@ -378,19 +368,16 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorContrastCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_.size() > 0 && xcorr_precursor_contrast_matrix_[0].size() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_.rows() > 0 && xcorr_precursor_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_precursor_contrast_matrix_.size(); i++) + for (auto e : xcorr_precursor_contrast_matrix_) { - for (std::size_t j = 0; j < xcorr_precursor_contrast_matrix_[0].size(); j++) - { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_[i][j])->first)); + deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(e)->first)); #ifdef MRMSCORING_TESTING - std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_[i][j])->first) << std::endl; + std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_[i][j])->first) << std::endl; #endif - } } OpenSwath::mean_and_stddev msc; @@ -404,21 +391,20 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorCombinedCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_precursor_combined_matrix_.size(); i++) + for (std::size_t i = 0; i < xcorr_precursor_combined_matrix_.rows(); i++) { - for (std::size_t j = i; j < xcorr_precursor_combined_matrix_.size(); j++) + for (std::size_t j = i; j < xcorr_precursor_combined_matrix_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_[i][j])->first)); + deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_.getValue(i, j))->first)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_[i][j])->first) << std::endl; #endif } } - OpenSwath::mean_and_stddev msc; msc = std::for_each(deltas.begin(), deltas.end(), msc); double deltas_mean = msc.mean(); @@ -436,45 +422,53 @@ namespace OpenSwath /// double MRMScoring::calcXcorrShapeScore() { - OPENSWATH_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); - std::vector intensities; - for (std::size_t i = 0; i < xcorr_matrix_.size(); i++) + size_t element_number{0}; + double intensities{0}; + for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) { - for (std::size_t j = i; j < xcorr_matrix_.size(); j++) + for (std::size_t j = i; j < xcorr_matrix_.rows(); j++) { // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second); + intensities +=Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->second; + element_number++; } - } + }/* OpenSwath::mean_and_stddev msc; msc = std::for_each(intensities.begin(), intensities.end(), msc); return msc.mean(); + double intensities{0}; + for(auto e : xcorr_matrix_) + { + intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; + }*/ + // xcorr_matrix_ is a triangle matrix + //size_t element_number= xcorr_matrix_.rows()*xcorr_matrix_.rows()/2 + (xcorr_matrix_.rows()+1)/2; + return intensities / element_number; } double MRMScoring::calcXcorrShapeWeightedScore( const std::vector& normalized_library_intensity) { - OPENSWATH_PRECONDITION(xcorr_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); // TODO (hroest) : check implementation // see _calc_weighted_xcorr_shape_score in MRM_pgroup.pm // -- they only multiply up the intensity once - std::vector intensities; - for (std::size_t i = 0; i < xcorr_matrix_.size(); i++) + double intensities{0}; + for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) { - intensities.push_back( - Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->second + intensities += (Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, i))->second * normalized_library_intensity[i] * normalized_library_intensity[i]); #ifdef MRMSCORING_TESTING std::cout << "_xcorr_weighted " << i << " " << i << " " << Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->second << " weight " << normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl; #endif - for (std::size_t j = i + 1; j < xcorr_matrix_.size(); j++) + for (std::size_t j = i + 1; j < xcorr_matrix_.rows(); j++) { - intensities.push_back( - Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second + intensities += (Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->second * normalized_library_intensity[i] * normalized_library_intensity[j] * 2); #ifdef MRMSCORING_TESTING @@ -483,41 +477,35 @@ namespace OpenSwath #endif } } - return std::accumulate(intensities.begin(), intensities.end(), 0.0); + return intensities; } double MRMScoring::calcXcorrContrastShapeScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.size() > 0 && xcorr_contrast_matrix_[0].size() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); - std::vector intensities; - for (std::size_t i = 0; i < xcorr_contrast_matrix_.size(); i++) + double intensities{0}; + for(auto e : xcorr_contrast_matrix_) { - for (std::size_t j = 0; j < xcorr_contrast_matrix_[0].size(); j++) - { - // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->second); - } + intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + return intensities; } std::vector MRMScoring::calcSeparateXcorrContrastShapeScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.size() > 0 && xcorr_contrast_matrix_[0].size() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector intensities; - for (std::size_t i = 0; i < xcorr_contrast_matrix_.size(); i++) + for (std::size_t i = 0; i < xcorr_contrast_matrix_.rows(); i++) { double intensities_id = 0; - for (std::size_t j = 0; j < xcorr_contrast_matrix_[0].size(); j++) + for (std::size_t j = 0; j < xcorr_contrast_matrix_.cols(); j++) { // second is the Y value (intensity) - intensities_id += Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->second; + intensities_id += Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_.getValue(i, j))->second; } - intensities.push_back(intensities_id / xcorr_contrast_matrix_[0].size()); + intensities.push_back(intensities_id / xcorr_contrast_matrix_.cols()); } return intensities; @@ -525,56 +513,56 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorShapeScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); - + OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + /* std::vector intensities; - for (std::size_t i = 0; i < xcorr_precursor_matrix_.size(); i++) + for (std::size_t i = 0; i < xcorr_precursor_matrix_.rows(); i++) { - for (std::size_t j = i; j < xcorr_precursor_matrix_.size(); j++) + for (std::size_t j = i; j < xcorr_precursor_matrix_.rows(); j++) { // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_[i][j])->second); + intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_.getValue(i, j))->second); } } OpenSwath::mean_and_stddev msc; msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + return msc.mean();*/ + double intensities{0}; + for(auto e : xcorr_precursor_matrix_) + { + intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; + } + //xcorr_precursor_matrix_ is a triangle matrix + size_t element_number = xcorr_precursor_matrix_.rows()*xcorr_precursor_matrix_.rows()/2 + (xcorr_precursor_matrix_.rows()+1)/2; + return intensities / element_number; } double MRMScoring::calcXcorrPrecursorContrastShapeScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_.size() > 0 && xcorr_precursor_contrast_matrix_[0].size() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_.rows() > 0 && xcorr_precursor_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); - std::vector intensities; - for (std::size_t i = 0; i < xcorr_precursor_contrast_matrix_.size(); i++) + + double intensities{0}; + for(auto e : xcorr_precursor_contrast_matrix_) { - for (std::size_t j = 0; j < xcorr_precursor_contrast_matrix_[0].size(); j++) - { - // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_[i][j])->second); - } + intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + return intensities / xcorr_precursor_contrast_matrix_.size(); } double MRMScoring::calcXcorrPrecursorCombinedShapeScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_.size() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); - std::vector intensities; - for (std::size_t i = 0; i < xcorr_precursor_combined_matrix_.size(); i++) + + double intensities{0}; + for(auto e : xcorr_precursor_combined_matrix_) { - for (std::size_t j = i; j < xcorr_precursor_combined_matrix_.size(); j++) - { - // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_[i][j])->second); - } + intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + //xc_precursor-combined_matrix_ is a triangle matrix + size_t element_number = xcorr_precursor_combined_matrix_.rows()*xcorr_precursor_combined_matrix_.rows()/2 + (xcorr_precursor_combined_matrix_.rows()+1)/2; + return intensities / element_number; } void MRMScoring::calcLibraryScore(OpenSwath::IMRMFeature* mrmfeature, const std::vector& transitions, @@ -688,22 +676,22 @@ namespace OpenSwath return sn_scores; } - const std::vector< std::vector >& MRMScoring::getMIMatrix() const + const OpenMS::Matrix & MRMScoring::getMIMatrix() const { return mi_matrix_; } - const std::vector< std::vector >& MRMScoring::getMIContrastMatrix() const + const OpenMS::Matrix & MRMScoring::getMIContrastMatrix() const { return mi_contrast_matrix_; } - const std::vector< std::vector >& MRMScoring::getMIPrecursorContrastMatrix() const + const OpenMS::Matrix & MRMScoring::getMIPrecursorContrastMatrix() const { return mi_precursor_contrast_matrix_; } - const std::vector< std::vector >& MRMScoring::getMIPrecursorCombinedMatrix() const + const OpenMS::Matrix & MRMScoring::getMIPrecursorCombinedMatrix() const { return mi_precursor_combined_matrix_; } @@ -711,12 +699,12 @@ namespace OpenSwath void MRMScoring::initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids) { std::vector intensityi, intensityj; - mi_matrix_.resize(native_ids.size()); + mi_matrix_.resize(native_ids.size(),native_ids.size()); for (std::size_t i = 0; i < native_ids.size(); i++) { String native_id = native_ids[i]; FeatureType fi = mrmfeature->getFeature(native_id); - mi_matrix_[i].resize(native_ids.size()); + intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < native_ids.size(); j++) @@ -726,7 +714,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information - mi_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_matrix_.setValue(i,j,Scoring::rankedMutualInformation(intensityi, intensityj)); } } } @@ -734,12 +722,12 @@ namespace OpenSwath void MRMScoring::initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids_set1, std::vector native_ids_set2) { std::vector intensityi, intensityj; - mi_contrast_matrix_.resize(native_ids_set1.size()); + mi_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); for (std::size_t i = 0; i < native_ids_set1.size(); i++) { String native_id = native_ids_set1[i]; FeatureType fi = mrmfeature->getFeature(native_id); - mi_contrast_matrix_[i].resize(native_ids_set2.size()); + //mi_contrast_matrix_[i].resize(native_ids_set2.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids_set2.size(); j++) @@ -749,7 +737,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information - mi_contrast_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj)); } } } @@ -757,12 +745,12 @@ namespace OpenSwath void MRMScoring::initializeMIPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector precursor_ids) { std::vector intensityi, intensityj; - mi_precursor_matrix_.resize(precursor_ids.size()); + mi_precursor_matrix_.resize(precursor_ids.size(),precursor_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { String precursor_id = precursor_ids[i]; FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - mi_precursor_matrix_[i].resize(precursor_ids.size()); + //mi_precursor_matrix_[i].resize(precursor_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < precursor_ids.size(); j++) @@ -772,7 +760,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information - mi_precursor_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_precursor_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj)); } } } @@ -780,12 +768,12 @@ namespace OpenSwath void MRMScoring::initializeMIPrecursorContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids, const std::vector& native_ids) { std::vector intensityi, intensityj; - mi_precursor_contrast_matrix_.resize(precursor_ids.size()); + mi_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { String precursor_id = precursor_ids[i]; FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - mi_precursor_contrast_matrix_[i].resize(native_ids.size()); + //mi_precursor_contrast_matrix_[i].resize(native_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids.size(); j++) @@ -795,7 +783,7 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information - mi_precursor_contrast_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_precursor_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj)); } } } @@ -818,11 +806,11 @@ namespace OpenSwath features.push_back(fj); } - mi_precursor_combined_matrix_.resize(features.size()); + mi_precursor_combined_matrix_.resize(features.size(), features.size()); for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; - mi_precursor_combined_matrix_[i].resize(features.size()); + //mi_precursor_combined_matrix_[i].resize(features.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < features.size(); j++) @@ -831,124 +819,108 @@ namespace OpenSwath intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information - mi_precursor_combined_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_precursor_combined_matrix_.setValue(i ,j, Scoring::rankedMutualInformation(intensityi, intensityj)); } } } double MRMScoring::calcMIScore() { - OPENSWATH_PRECONDITION(mi_matrix_.size() > 1, "Expect mutual information matrix of at least 2x2"); + OPENSWATH_PRECONDITION(mi_matrix_.rows() > 1, "Expect mutual information matrix of at least 2x2"); - std::vector mi_scores; - for (std::size_t i = 0; i < mi_matrix_.size(); i++) + double mi_scores{0}; + for(auto e : mi_matrix_) { - for (std::size_t j = i; j < mi_matrix_.size(); j++) - { - mi_scores.push_back(mi_matrix_[i][j]); - } + mi_scores += e; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(mi_scores.begin(), mi_scores.end(), msc); - return msc.mean(); + //mi_matrix_ is a triangular matrix + size_t element_number = mi_matrix_.rows() * mi_matrix_.rows() / 2 + (mi_matrix_.rows() + 1) / 2; + return mi_scores / element_number; } double MRMScoring::calcMIWeightedScore( const std::vector& normalized_library_intensity) { - OPENSWATH_PRECONDITION(mi_matrix_.size() > 1, "Expect mutual information matrix of at least 2x2"); + OPENSWATH_PRECONDITION(mi_matrix_.rows() > 1, "Expect mutual information matrix of at least 2x2"); - std::vector mi_scores; - for (std::size_t i = 0; i < mi_matrix_.size(); i++) + double mi_scores{0}; + for (std::size_t i = 0; i < mi_matrix_.rows(); i++) { - mi_scores.push_back(mi_matrix_[i][i] + mi_scores += mi_matrix_.getValue(i,i) * normalized_library_intensity[i] - * normalized_library_intensity[i]); + * normalized_library_intensity[i]; #ifdef MRMSCORING_TESTING std::cout << "_mi_weighted " << i << " " << i << " " << mi_matrix_[i][i] << " weight " << normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl; #endif - for (std::size_t j = i + 1; j < mi_matrix_.size(); j++) + for (std::size_t j = i + 1; j < mi_matrix_.rows(); j++) { - mi_scores.push_back(mi_matrix_[i][j] + mi_scores += mi_matrix_.getValue(i,j) * normalized_library_intensity[i] - * normalized_library_intensity[j] * 2); + * normalized_library_intensity[j] * 2; #ifdef MRMSCORING_TESTING std::cout << "_mi_weighted " << i << " " << j << " " << mi_matrix_[i][j] << " weight " << normalized_library_intensity[i] * normalized_library_intensity[j] * 2 << std::endl; #endif } } - return std::accumulate(mi_scores.begin(), mi_scores.end(), 0.0); + return mi_scores; } double MRMScoring::calcMIPrecursorScore() { - OPENSWATH_PRECONDITION(mi_precursor_matrix_.size() > 1, "Expect mutual information matrix of at least 2x2"); + OPENSWATH_PRECONDITION(mi_precursor_matrix_.rows() > 1, "Expect mutual information matrix of at least 2x2"); - std::vector mi_scores; - for (std::size_t i = 0; i < mi_precursor_matrix_.size(); i++) + double mi_scores{0}; + for(auto e : mi_precursor_matrix_) { - for (std::size_t j = i; j < mi_precursor_matrix_.size(); j++) - { - mi_scores.push_back(mi_precursor_matrix_[i][j]); - } + mi_scores += e; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(mi_scores.begin(), mi_scores.end(), msc); - return msc.mean(); + //mi_precursor_matrix_ is a triangular matrix + size_t element_number = mi_precursor_matrix_.rows()*mi_precursor_matrix_.rows()/2 + (mi_precursor_matrix_.rows()+1)/2; + return mi_scores / element_number; } double MRMScoring::calcMIPrecursorContrastScore() { - OPENSWATH_PRECONDITION(mi_precursor_contrast_matrix_.size() > 0 && mi_precursor_contrast_matrix_[0].size() > 1, "Expect mutual information matrix of at least 1x2"); + OPENSWATH_PRECONDITION(mi_precursor_contrast_matrix_.rows() > 0 && mi_precursor_contrast_matrix_.cols() > 1, "Expect mutual information matrix of at least 1x2"); - std::vector mi_scores; - for (std::size_t i = 0; i < mi_precursor_contrast_matrix_.size(); i++) + double mi_scores{0}; + for(auto e : mi_precursor_contrast_matrix_) { - for (std::size_t j = 0; j < mi_precursor_contrast_matrix_[0].size(); j++) - { - mi_scores.push_back(mi_precursor_contrast_matrix_[i][j]); - } + mi_scores += e; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(mi_scores.begin(), mi_scores.end(), msc); - return msc.mean(); + return mi_scores / mi_precursor_contrast_matrix_.size(); } double MRMScoring::calcMIPrecursorCombinedScore() { - OPENSWATH_PRECONDITION(mi_precursor_combined_matrix_.size() > 1, "Expect mutual information matrix of at least 2x2"); + OPENSWATH_PRECONDITION(mi_precursor_combined_matrix_.rows() > 1, "Expect mutual information matrix of at least 2x2"); - std::vector mi_scores; - for (std::size_t i = 0; i < mi_precursor_combined_matrix_.size(); i++) + double mi_scores{0}; + + for(auto e: mi_precursor_combined_matrix_) { - for (std::size_t j = 0; j < mi_precursor_combined_matrix_[0].size(); j++) - { - mi_scores.push_back(mi_precursor_combined_matrix_[i][j]); - } + mi_scores += e; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(mi_scores.begin(), mi_scores.end(), msc); - return msc.mean(); + return mi_scores / mi_precursor_combined_matrix_.size(); } std::vector MRMScoring::calcSeparateMIContrastScore() { - OPENSWATH_PRECONDITION(mi_contrast_matrix_.size() > 0 && mi_contrast_matrix_[0].size() > 1, "Expect mutual information matrix of at least 1x2"); + OPENSWATH_PRECONDITION(mi_contrast_matrix_.rows() > 0 && mi_contrast_matrix_.cols() > 1, "Expect mutual information matrix of at least 1x2"); std::vector mi_scores; - for (std::size_t i = 0; i < mi_contrast_matrix_.size(); i++) + mi_scores.resize(mi_contrast_matrix_.rows()); + for (std::size_t i = 0; i < mi_contrast_matrix_.rows(); i++) { double mi_scores_id = 0; - for (std::size_t j = 0; j < mi_contrast_matrix_[0].size(); j++) + for (std::size_t j = 0; j < mi_contrast_matrix_.cols(); j++) { - mi_scores_id += mi_contrast_matrix_[i][j]; + mi_scores_id += mi_contrast_matrix_.getValue(i,j); } - mi_scores.push_back(mi_scores_id / mi_contrast_matrix_[0].size()); + mi_scores[i] = mi_scores_id / mi_contrast_matrix_.cols(); } - return mi_scores; } - } diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp index d150306b99c..03c07d18c54 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp @@ -159,7 +159,7 @@ namespace OpenSwath::Scoring } XCorrArrayType normalizedCrossCorrelation(std::vector& data1, - std::vector& data2, const int& maxdelay, const int& lag = 1) + std::vector& data2, int maxdelay, int lag = 1) //const ref entfernt { OPENSWATH_PRECONDITION(data1.size() != 0 && data1.size() == data2.size(), "Both data vectors need to have the same length"); @@ -175,7 +175,7 @@ namespace OpenSwath::Scoring } XCorrArrayType calculateCrossCorrelation(const std::vector& data1, - const std::vector& data2, const int& maxdelay, const int& lag) + const std::vector& data2, int maxdelay, int lag) //const ref entfernt { OPENSWATH_PRECONDITION(data1.size() != 0 && data1.size() == data2.size(), "Both data vectors need to have the same length"); @@ -270,6 +270,29 @@ namespace OpenSwath::Scoring std::vector computeRank(const std::vector& v_temp) { + /*std::vector ranks{}; + ranks.resize(v_temp.size()); + std::iota(ranks.begin(), ranks.end(), 0); + std::sort(ranks.begin(), ranks.end(), + [&v_temp](unsigned int i, unsigned int j) { return v_temp[i] < v_temp[j]; }); + unsigned int tmp = 0; + for (unsigned int i = 1; i < ranks.size(); ++i) + { + if (v_temp[tmp] == v_temp[ranks[i]]) + { + tmp = ranks[i]; + ranks[i] = ranks[i-1]; + } + else + { + tmp = ranks[i]; + ranks[i] = i; + } + } + + return ranks; + }*/ + std::vector > v_sort(v_temp.size()); for (unsigned int i = 0; i < v_sort.size(); ++i) { diff --git a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp index a2dceee4a8a..2922ea1cc3e 100644 --- a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp +++ b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp @@ -70,23 +70,22 @@ void fill_mock_objects(MockMRMFeature * imrmfeature, std::vector& n native_ids.push_back("group1"); native_ids.push_back("group2"); + /* static const double arr1[] = { 5.97543668746948, 4.2749171257019, 3.3301842212677, 4.08597040176392, 5.50307035446167, 5.24326848983765, 8.40812492370605, 2.83419919013977, 6.94378805160522, 7.69957494735718, 4.08597040176392 }; + std::vector intensity1 (arr1, arr1 + sizeof(arr1) / sizeof(arr1[0]) ); - static const double arr2[] = - { - 15.8951349258423, 41.5446395874023, 76.0746307373047, 109.069435119629, 111.90364074707, 169.79216003418, - 121.043930053711, 63.0136985778809, 44.6150207519531, 21.4926776885986, 7.93575811386108 - }; - std::vector intensity2 (arr2, arr2 + sizeof(arr2) / sizeof(arr2[0]) ); - static const double arr3[] = - { - 0.0, 110.0, 200.0, 270.0, 320.0, 350.0, 360.0, 350.0, 320.0, 270.0, 200.0 - }; - std::vector ms1intensity (arr3, arr3 + sizeof(arr3) / sizeof(arr3[0]) ); + */ + std::vector intensity1 {5.97543668746948, 4.2749171257019, 3.3301842212677, 4.08597040176392, 5.50307035446167, 5.24326848983765, + 8.40812492370605, 2.83419919013977, 6.94378805160522, 7.69957494735718, 4.08597040176392}; + + std::vector intensity2 {15.8951349258423, 41.5446395874023, 76.0746307373047, 109.069435119629, 111.90364074707, 169.79216003418, + 121.043930053711, 63.0136985778809, 44.6150207519531, 21.4926776885986, 7.93575811386108}; + + std::vector ms1intensity {0.0, 110.0, 200.0, 270.0, 320.0, 350.0, 360.0, 350.0, 320.0, 270.0, 200.0}; boost::shared_ptr f1_ptr = boost::shared_ptr(new MockFeature()); boost::shared_ptr f2_ptr = boost::shared_ptr(new MockFeature()); @@ -113,33 +112,18 @@ void fill_mock_objects2(MockMRMFeature * imrmfeature, std::vector& native_ids.push_back("group1"); native_ids.push_back("group2"); - static const double arr1[] = - { - 5.97543668746948, 4.2749171257019, 3.3301842212677, 4.08597040176392, 5.50307035446167, 5.24326848983765, - 8.40812492370605, 2.83419919013977, 6.94378805160522, 7.69957494735718, 4.08597040176392 - }; - std::vector intensity1 (arr1, arr1 + sizeof(arr1) / sizeof(arr1[0]) ); - static const double arr2[] = - { - 15.8951349258423, 41.5446395874023, 76.0746307373047, 109.069435119629, 111.90364074707, 169.79216003418, - 121.043930053711, 63.0136985778809, 44.6150207519531, 21.4926776885986, 7.93575811386108 - }; - std::vector intensity2 (arr2, arr2 + sizeof(arr2) / sizeof(arr2[0]) ); - static const double arr3[] = - { - 0.0, 110.0, 200.0, 270.0, 320.0, 350.0, 360.0, 350.0, 320.0, 270.0, 200.0 - }; - std::vector ms1intensity1 (arr3, arr3 + sizeof(arr3) / sizeof(arr3[0]) ); - static const double arr4[] = - { - 10.0, 115.0, 180.0, 280.0, 330.0, 340.0, 390.0, 320.0, 300.0, 250.0, 100.0 - }; - std::vector ms1intensity2 (arr4, arr4 + sizeof(arr4) / sizeof(arr4[0]) ); - static const double arr5[] = - { - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 - }; - std::vector ms1intensity3 (arr5, arr5 + sizeof(arr5) / sizeof(arr5[0]) ); + + std::vector intensity1 {5.97543668746948, 4.2749171257019, 3.3301842212677, 4.08597040176392, 5.50307035446167, 5.24326848983765, + 8.40812492370605, 2.83419919013977, 6.94378805160522, 7.69957494735718, 4.08597040176392}; + + std::vector intensity2 {15.8951349258423, 41.5446395874023, 76.0746307373047, 109.069435119629, 111.90364074707, 169.79216003418, + 121.043930053711, 63.0136985778809, 44.6150207519531, 21.4926776885986, 7.93575811386108}; + + std::vector ms1intensity1 {0.0, 110.0, 200.0, 270.0, 320.0, 350.0, 360.0, 350.0, 320.0, 270.0, 200.0}; + + std::vector ms1intensity2 {0.0, 115.0, 180.0, 280.0, 330.0, 340.0, 390.0, 320.0, 300.0, 250.0, 100.0}; + + std::vector ms1intensity3 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; boost::shared_ptr f1_ptr = boost::shared_ptr(new MockFeature()); boost::shared_ptr f2_ptr = boost::shared_ptr(new MockFeature()); @@ -251,13 +235,13 @@ BOOST_AUTO_TEST_CASE(initializeXCorrMatrix) mrmscore.initializeXCorrMatrix(imrmfeature, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getXCorrMatrix().size(), 2) - TEST_EQUAL(mrmscore.getXCorrMatrix()[0].size(), 2) - TEST_EQUAL(mrmscore.getXCorrMatrix()[0][0].data.size(), 23) + TEST_EQUAL(mrmscore.getXCorrMatrix().rows(), 2) + TEST_EQUAL(mrmscore.getXCorrMatrix().cols(), 2) + TEST_EQUAL(mrmscore.getXCorrMatrix().getValue(0, 0).data.size(), 23) // test auto-correlation = xcorrmatrix_0_0 const OpenSwath::Scoring::XCorrArrayType auto_correlation = - mrmscore.getXCorrMatrix()[0][0]; + mrmscore.getXCorrMatrix().getValue(0, 0); TEST_EQUAL(auto_correlation.data[11].first, 0) TEST_EQUAL(auto_correlation.data[12].first, 1) @@ -273,7 +257,7 @@ BOOST_AUTO_TEST_CASE(initializeXCorrMatrix) // test cross-correlation = xcorrmatrix_0_1 const OpenSwath::Scoring::XCorrArrayType cross_correlation = - mrmscore.getXCorrMatrix()[0][1]; + mrmscore.getXCorrMatrix().getValue(0, 1); TEST_REAL_SIMILAR(cross_correlation.data[13].second, -0.31165141) // find(2)->second, TEST_REAL_SIMILAR(cross_correlation.data[12].second, -0.35036919) // find(1)->second, @@ -298,8 +282,8 @@ BOOST_AUTO_TEST_CASE(initializeXCorrPrecursorContrastMatrix) mrmscore.initializeXCorrPrecursorContrastMatrix(imrmfeature, precursor_ids, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getXCorrPrecursorContrastMatrix().size(), 3) - TEST_EQUAL(mrmscore.getXCorrPrecursorContrastMatrix()[0].size(), 2) + TEST_EQUAL(mrmscore.getXCorrPrecursorContrastMatrix().rows(), 3) + TEST_EQUAL(mrmscore.getXCorrPrecursorContrastMatrix().cols(), 2) } END_SECTION @@ -316,8 +300,8 @@ BOOST_AUTO_TEST_CASE(initializeXCorrPrecursorCombinedMatrix) mrmscore.initializeXCorrPrecursorCombinedMatrix(imrmfeature, precursor_ids, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getXCorrPrecursorCombinedMatrix().size(), 5) - TEST_EQUAL(mrmscore.getXCorrPrecursorCombinedMatrix()[0].size(), 5) + TEST_EQUAL(mrmscore.getXCorrPrecursorCombinedMatrix().rows(), 5) + TEST_EQUAL(mrmscore.getXCorrPrecursorCombinedMatrix().cols(), 5) } END_SECTION @@ -333,13 +317,13 @@ BOOST_AUTO_TEST_CASE(initializeXCorrContrastMatrix) mrmscore.initializeXCorrContrastMatrix(imrmfeature, native_ids, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getXCorrContrastMatrix().size(), 2) - TEST_EQUAL(mrmscore.getXCorrContrastMatrix()[0].size(), 2) - TEST_EQUAL(mrmscore.getXCorrContrastMatrix()[0][0].data.size(), 23) + TEST_EQUAL(mrmscore.getXCorrContrastMatrix().rows(), 2) + TEST_EQUAL(mrmscore.getXCorrContrastMatrix().cols(), 2) + TEST_EQUAL(mrmscore.getXCorrContrastMatrix().getValue(0, 0).data.size(), 23) // test auto-correlation = xcorrmatrix_0_0 const OpenSwath::Scoring::XCorrArrayType auto_correlation = - mrmscore.getXCorrContrastMatrix()[0][0]; + mrmscore.getXCorrContrastMatrix().getValue(0, 0); TEST_REAL_SIMILAR(auto_correlation.data[11].second, 1) // find(0)->second, TEST_REAL_SIMILAR(auto_correlation.data[12].second, -0.227352707759245) // find(1)->second, TEST_REAL_SIMILAR(auto_correlation.data[10].second, -0.227352707759245) // find(-1)->second, @@ -348,7 +332,7 @@ BOOST_AUTO_TEST_CASE(initializeXCorrContrastMatrix) // // test cross-correlation = xcorrmatrix_0_1 const OpenSwath::Scoring::XCorrArrayType cross_correlation = - mrmscore.getXCorrContrastMatrix()[0][1]; + mrmscore.getXCorrContrastMatrix().getValue(0, 1); TEST_REAL_SIMILAR(cross_correlation.data[13].second, -0.31165141) // find(2)->second, TEST_REAL_SIMILAR(cross_correlation.data[12].second, -0.35036919) // find(1)->second, TEST_REAL_SIMILAR(cross_correlation.data[11].second, 0.03129565) // find(0)->second, @@ -682,13 +666,13 @@ mean(m4) mrmscore.initializeMIMatrix(imrmfeature, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getMIMatrix().size(), 2) - TEST_EQUAL(mrmscore.getMIMatrix()[0].size(), 2) + TEST_EQUAL(mrmscore.getMIMatrix().rows(), 2) + TEST_EQUAL(mrmscore.getMIMatrix().cols(), 2) - TEST_REAL_SIMILAR(mrmscore.getMIMatrix()[0][0], 3.2776) - TEST_REAL_SIMILAR(mrmscore.getMIMatrix()[0][1], 3.2776) - TEST_REAL_SIMILAR(mrmscore.getMIMatrix()[1][1], 3.4594) - TEST_REAL_SIMILAR(mrmscore.getMIMatrix()[1][0], 0) // value not initialized for lower diagonal half of matrix + TEST_REAL_SIMILAR(mrmscore.getMIMatrix().getValue(0, 0), 3.2776) + TEST_REAL_SIMILAR(mrmscore.getMIMatrix().getValue(0, 1), 3.2776) + TEST_REAL_SIMILAR(mrmscore.getMIMatrix().getValue(1, 1), 3.4594) + TEST_REAL_SIMILAR(mrmscore.getMIMatrix().getValue(1, 0), 0) // value not initialized for lower diagonal half of matrix } END_SECTION @@ -705,8 +689,8 @@ BOOST_AUTO_TEST_CASE(initializeMIPrecursorContrastMatrix) mrmscore.initializeMIPrecursorContrastMatrix(imrmfeature, precursor_ids, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getMIPrecursorContrastMatrix().size(), 3) - TEST_EQUAL(mrmscore.getMIPrecursorContrastMatrix()[0].size(), 2) + TEST_EQUAL(mrmscore.getMIPrecursorContrastMatrix().rows(), 3) + TEST_EQUAL(mrmscore.getMIPrecursorContrastMatrix().cols(), 2) } END_SECTION @@ -723,8 +707,8 @@ BOOST_AUTO_TEST_CASE(initializeMIPrecursorCombinedMatrix) mrmscore.initializeMIPrecursorCombinedMatrix(imrmfeature, precursor_ids, native_ids); delete imrmfeature; - TEST_EQUAL(mrmscore.getMIPrecursorCombinedMatrix().size(), 5) - TEST_EQUAL(mrmscore.getMIPrecursorCombinedMatrix()[0].size(), 5) + TEST_EQUAL(mrmscore.getMIPrecursorCombinedMatrix().rows(), 5) + TEST_EQUAL(mrmscore.getMIPrecursorCombinedMatrix().cols(), 5) } END_SECTION @@ -740,10 +724,10 @@ BOOST_AUTO_TEST_CASE(initializeMIContrastMatrix) mrmscore.initializeMIContrastMatrix(imrmfeature, native_ids, native_ids); delete imrmfeature; - TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix()[0][0], 3.2776) - TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix()[0][1], 3.2776) - TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix()[1][1], 3.4594) - TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix()[1][0], 3.2776) + TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix().getValue(0, 0), 3.2776) + TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix().getValue(0, 1), 3.2776) + TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix().getValue(1, 1), 3.4594) + TEST_REAL_SIMILAR(mrmscore.getMIContrastMatrix().getValue(1, 0), 3.2776) } END_SECTION diff --git a/src/tests/class_tests/openswathalgo/Scoring_test.cpp b/src/tests/class_tests/openswathalgo/Scoring_test.cpp index ea55f53ad3d..f5acb0b5266 100644 --- a/src/tests/class_tests/openswathalgo/Scoring_test.cpp +++ b/src/tests/class_tests/openswathalgo/Scoring_test.cpp @@ -388,18 +388,18 @@ y = [5.97543668746948 4.2749171257019 3.3301842212677 4.08597040176392 5.5030703 */ - static const double arr1[] = - { - 5.97543668746948, 4.2749171257019, 3.3301842212677, 4.08597040176392, 5.50307035446167, 5.24326848983765, - 8.40812492370605, 2.83419919013977, 6.94378805160522, 7.69957494735718, 4.08597040176392 - }; - static const double arr2[] = + + + /*static const double arr2[] = { 15.8951349258423, 41.5446395874023, 76.0746307373047, 109.069435119629, 111.90364074707, 169.79216003418, 121.043930053711, 63.0136985778809, 44.6150207519531, 21.4926776885986, 7.93575811386108 - }; - std::vector data1 (arr1, arr1 + sizeof(arr1) / sizeof(arr1[0]) ); - std::vector data2 (arr2, arr2 + sizeof(arr2) / sizeof(arr2[0]) ); + };*/ + std::vector data1 {5.97543668746948, 4.2749171257019, 3.3301842212677, 4.08597040176392, 5.50307035446167, 5.24326848983765, + 8.40812492370605, 2.83419919013977, 6.94378805160522, 7.69957494735718, 4.08597040176392}; + + std::vector data2 {15.8951349258423, 41.5446395874023, 76.0746307373047, 109.069435119629, 111.90364074707, 169.79216003418, + 121.043930053711, 63.0136985778809, 44.6150207519531, 21.4926776885986, 7.93575811386108}; auto result = Scoring::computeRank(data1); From 9f9c7719dda5a92c48af39efd5b92feb02f34073 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Mon, 25 Oct 2021 19:11:00 +0200 Subject: [PATCH 02/11] const und ref --- .../include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h | 6 +++--- .../include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h | 4 ++-- .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h index 37f2ac528f5..a471af8d8f8 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h @@ -237,13 +237,13 @@ namespace OpenSwath //@} /// Initialize the scoring object and building the MI matrix - void initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids); //ref? + void initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids); /// Initialize the scoring object and building the MI matrix of chromatograms of set1 (e.g. identification transitions) vs set2 (e.g. detection transitions) - void initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids_set1, std::vector native_ids_set2); + void initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids_set1, const std::vector& native_ids_set2); /// Initialize the scoring object and building the MI matrix - void initializeMIPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector precursor_ids); + void initializeMIPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids); /// Initialize the mutual information vector against the MS1 trace void initializeMIPrecursorContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids, const std::vector& native_ids); diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h index 4abfc7aa123..6e9ed2c9b2a 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h @@ -117,11 +117,11 @@ namespace OpenSwath /// Calculate crosscorrelation on std::vector data (which is first normalized) /// NOTE: this replaces calcxcorr OPENSWATHALGO_DLLAPI XCorrArrayType normalizedCrossCorrelation(std::vector& data1, - std::vector& data2, int maxdelay, int lag); + std::vector& data2, const int maxdelay, const int lag); /// Calculate crosscorrelation on std::vector data without normalization OPENSWATHALGO_DLLAPI XCorrArrayType calculateCrossCorrelation(const std::vector& data1, - const std::vector& data2, int maxdelay, int lag); + const std::vector& data2, const int maxdelay, const int lag); /// Find best peak in an cross-correlation (highest apex) OPENSWATHALGO_DLLAPI XCorrArrayType::const_iterator xcorrArrayGetMaxPeak(const XCorrArrayType & array); diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 848daf4d8a6..54b35feeaad 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -560,7 +560,7 @@ namespace OpenSwath { intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; } - //xc_precursor-combined_matrix_ is a triangle matrix + //xcorr_precursor-combined_matrix_ is a triangle matrix size_t element_number = xcorr_precursor_combined_matrix_.rows()*xcorr_precursor_combined_matrix_.rows()/2 + (xcorr_precursor_combined_matrix_.rows()+1)/2; return intensities / element_number; } @@ -696,7 +696,7 @@ namespace OpenSwath return mi_precursor_combined_matrix_; } - void MRMScoring::initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids) + void MRMScoring::initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids) { std::vector intensityi, intensityj; mi_matrix_.resize(native_ids.size(),native_ids.size()); @@ -719,7 +719,7 @@ namespace OpenSwath } } - void MRMScoring::initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector native_ids_set1, std::vector native_ids_set2) + void MRMScoring::initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids_set1, const std::vector& native_ids_set2) { std::vector intensityi, intensityj; mi_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); @@ -742,7 +742,7 @@ namespace OpenSwath } } - void MRMScoring::initializeMIPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, std::vector precursor_ids) + void MRMScoring::initializeMIPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids) { std::vector intensityi, intensityj; mi_precursor_matrix_.resize(precursor_ids.size(),precursor_ids.size()); From 73e9cceec5c186d2f88fa0d1a9889f2a44487f1e Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Tue, 26 Oct 2021 10:29:41 +0200 Subject: [PATCH 03/11] fix xcorr_precursor_combined & contrast --- src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 7 +++++-- src/tests/class_tests/openswathalgo/MRMScoring_test.cpp | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 54b35feeaad..13c0d1d67a3 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -556,9 +556,12 @@ namespace OpenSwath double intensities{0}; - for(auto e : xcorr_precursor_combined_matrix_) + for(size_t i = 0; i < xcorr_precursor_combined_matrix_.rows(); i++) { - intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; + for(size_t j = i; j < xcorr_precursor_combined_matrix_.cols(); j++) + { + intensities += Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_.getValue(i,j))->second; + } } //xcorr_precursor-combined_matrix_ is a triangle matrix size_t element_number = xcorr_precursor_combined_matrix_.rows()*xcorr_precursor_combined_matrix_.rows()/2 + (xcorr_precursor_combined_matrix_.rows()+1)/2; diff --git a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp index 2922ea1cc3e..4993bab6b89 100644 --- a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp +++ b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp @@ -121,7 +121,7 @@ void fill_mock_objects2(MockMRMFeature * imrmfeature, std::vector& std::vector ms1intensity1 {0.0, 110.0, 200.0, 270.0, 320.0, 350.0, 360.0, 350.0, 320.0, 270.0, 200.0}; - std::vector ms1intensity2 {0.0, 115.0, 180.0, 280.0, 330.0, 340.0, 390.0, 320.0, 300.0, 250.0, 100.0}; + std::vector ms1intensity2 {10.0, 115.0, 180.0, 280.0, 330.0, 340.0, 390.0, 320.0, 300.0, 250.0, 100.0}; std::vector ms1intensity3 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; From a13267384c4b785b803f73be1dc3982fb52e80c0 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Tue, 26 Oct 2021 13:27:55 +0200 Subject: [PATCH 04/11] doppelungen in initialize entfernt --- .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 98 +++++-------------- 1 file changed, 27 insertions(+), 71 deletions(-) diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 13c0d1d67a3..9b905cc3cc2 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -86,14 +86,12 @@ namespace OpenSwath xcorr_matrix_.resize(native_ids.size(), native_ids.size()); for (std::size_t i = 0; i < native_ids.size(); i++) { - String native_id = native_ids[i]; - FeatureType fi = mrmfeature->getFeature(native_id); + FeatureType fi = mrmfeature->getFeature(native_ids[i]); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < native_ids.size(); j++) { - String native_id2 = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id2); + FeatureType fj = mrmfeature->getFeature(native_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation @@ -108,14 +106,12 @@ namespace OpenSwath xcorr_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); for (std::size_t i = 0; i < native_ids_set1.size(); i++) { - String native_id = native_ids_set1[i]; - FeatureType fi = mrmfeature->getFeature(native_id); + FeatureType fi = mrmfeature->getFeature(native_ids_set1[i]); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids_set2.size(); j++) { - String native_id2 = native_ids_set2[j]; - FeatureType fj = mrmfeature->getFeature(native_id2); + FeatureType fj = mrmfeature->getFeature(native_ids_set2[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation @@ -130,14 +126,12 @@ namespace OpenSwath xcorr_precursor_matrix_.resize(precursor_ids.size(), precursor_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < precursor_ids.size(); j++) { - String precursor_id2 = precursor_ids[j]; - FeatureType fj = mrmfeature->getPrecursorFeature(precursor_id2); + FeatureType fj = mrmfeature->getPrecursorFeature(precursor_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation @@ -151,15 +145,13 @@ namespace OpenSwath std::vector intensityi, intensityj; xcorr_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) - { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); + { + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids.size(); j++) { - String native_id = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id); + FeatureType fj = mrmfeature->getFeature(native_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute normalized cross correlation @@ -192,15 +184,13 @@ namespace OpenSwath std::vector features; for (std::size_t i = 0; i < precursor_ids.size(); i++) - { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); + { + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); features.push_back(fi); } for (std::size_t j = 0; j < native_ids.size(); j++) { - String native_id = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id); + FeatureType fj = mrmfeature->getFeature(native_ids[j]); features.push_back(fj); } @@ -208,7 +198,6 @@ namespace OpenSwath for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; - //xcorr_precursor_combined_matrix_[i].resize(features.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < features.size(); j++) @@ -434,17 +423,7 @@ namespace OpenSwath intensities +=Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->second; element_number++; } - }/* - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); - double intensities{0}; - for(auto e : xcorr_matrix_) - { - intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; - }*/ - // xcorr_matrix_ is a triangle matrix - //size_t element_number= xcorr_matrix_.rows()*xcorr_matrix_.rows()/2 + (xcorr_matrix_.rows()+1)/2; + } return intensities / element_number; } @@ -514,19 +493,7 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorShapeScore() { OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); - /* - std::vector intensities; - for (std::size_t i = 0; i < xcorr_precursor_matrix_.rows(); i++) - { - for (std::size_t j = i; j < xcorr_precursor_matrix_.rows(); j++) - { - // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_.getValue(i, j))->second); - } - } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean();*/ + double intensities{0}; for(auto e : xcorr_precursor_matrix_) { @@ -705,15 +672,13 @@ namespace OpenSwath mi_matrix_.resize(native_ids.size(),native_ids.size()); for (std::size_t i = 0; i < native_ids.size(); i++) { - String native_id = native_ids[i]; - FeatureType fi = mrmfeature->getFeature(native_id); + FeatureType fi = mrmfeature->getFeature(native_ids[i]); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < native_ids.size(); j++) { - String native_id2 = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id2); + FeatureType fj = mrmfeature->getFeature(native_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information @@ -727,16 +692,14 @@ namespace OpenSwath std::vector intensityi, intensityj; mi_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); for (std::size_t i = 0; i < native_ids_set1.size(); i++) - { - String native_id = native_ids_set1[i]; - FeatureType fi = mrmfeature->getFeature(native_id); + { + FeatureType fi = mrmfeature->getFeature(native_ids_set1[i]); //mi_contrast_matrix_[i].resize(native_ids_set2.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids_set2.size(); j++) { - String native_id2 = native_ids_set2[j]; - FeatureType fj = mrmfeature->getFeature(native_id2); + FeatureType fj = mrmfeature->getFeature(native_ids_set2[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information @@ -751,15 +714,12 @@ namespace OpenSwath mi_precursor_matrix_.resize(precursor_ids.size(),precursor_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - //mi_precursor_matrix_[i].resize(precursor_ids.size()); + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = i; j < precursor_ids.size(); j++) { - String precursor_id2 = precursor_ids[j]; - FeatureType fj = mrmfeature->getPrecursorFeature(precursor_id2); + FeatureType fj = mrmfeature->getPrecursorFeature(precursor_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information @@ -773,16 +733,14 @@ namespace OpenSwath std::vector intensityi, intensityj; mi_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) - { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); + { + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); //mi_precursor_contrast_matrix_[i].resize(native_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < native_ids.size(); j++) { - String native_id = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id); + FeatureType fj = mrmfeature->getFeature(native_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); // compute ranked mutual information @@ -797,15 +755,13 @@ namespace OpenSwath std::vector features; for (std::size_t i = 0; i < precursor_ids.size(); i++) - { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); + { + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); features.push_back(fi); } for (std::size_t j = 0; j < native_ids.size(); j++) { - String native_id = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id); + FeatureType fj = mrmfeature->getFeature(native_ids[j]); features.push_back(fj); } From f3c2b09915cf73f2fb05bd94d6f605cfacef0a03 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Thu, 28 Oct 2021 13:03:21 +0200 Subject: [PATCH 05/11] neue Matritzen, Scoring wird direkt bei der initialisierung berechnet --- .../OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h | 16 +- .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 157 +++++++++++------- .../openswathalgo/MRMScoring_test.cpp | 85 ++++++++++ 3 files changed, 197 insertions(+), 61 deletions(-) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h index a471af8d8f8..11ce8d88aa9 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h @@ -267,21 +267,35 @@ namespace OpenSwath /// the precomputed cross correlation matrix XCorrMatrixType xcorr_matrix_; + /// contains max Peaks from xcorr_matrix_ + OpenMS::Matrix xcorr_matrix_max_peak_; + OpenMS::Matrix xcorr_matrix_max_peak_sec_; /// the precomputed contrast cross correlation XCorrMatrixType xcorr_contrast_matrix_; //@} + /// contains max Peaks from xcorr_contrast_matrix_ + OpenMS::Matrix xcorr_contrast_matrix_max_peak_; + OpenMS::Matrix xcorr_contrast_matrix_max_peak_sec_; /// the precomputed cross correlation matrix of the MS1 trace XCorrMatrixType xcorr_precursor_matrix_; + /// contains max Peaks from xcorr_precursor_matrix_ + OpenMS::Matrix xcorr_precursor_matrix_max_peak_; + OpenMS::Matrix xcorr_precursor_matrix_max_peak_sec_; /// the precomputed cross correlation against the MS1 trace XCorrMatrixType xcorr_precursor_contrast_matrix_; //@} + /// contains max Peaks from xcorr_precursor_contrast_matrix_ + OpenMS::Matrix xcorr_precursor_contrast_matrix_max_peak_; + OpenMS::Matrix xcorr_precursor_contrast_matrix_max_peak_sec_; /// the precomputed cross correlation with the MS1 trace XCorrMatrixType xcorr_precursor_combined_matrix_; //@} - + /// contains max Peaks from xcorr_precursor_combined_matrix_; + OpenMS::Matrix xcorr_precursor_combined_matrix_max_peak_; + OpenMS::Matrix xcorr_precursor_combined_matrix_max_peak_sec_; /// the precomputed mutual information matrix //std::vector< std::vector > mi_matrix_; OpenMS::Matrix mi_matrix_; diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 9b905cc3cc2..800765b91e2 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -53,6 +53,9 @@ namespace OpenSwath void MRMScoring::initializeXCorrMatrix(const std::vector< std::vector< double > >& data) { xcorr_matrix_.resize(data.size(), data.size()); + xcorr_matrix_max_peak_.resize(data.size(), data.size()); + xcorr_matrix_max_peak_sec_.resize(data.size(), data.size()); + for (std::size_t i = 0; i < data.size(); i++) { std::vector< double > tmp1(data[i]); @@ -61,6 +64,9 @@ namespace OpenSwath // compute normalized cross correlation std::vector< double > tmp2(data[j]); xcorr_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(tmp1, tmp2, boost::numeric_cast(data[i].size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j)); + xcorr_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_matrix_max_peak_sec_.setValue(i, j, x->second); } } } @@ -84,6 +90,8 @@ namespace OpenSwath { std::vector intensityi, intensityj; xcorr_matrix_.resize(native_ids.size(), native_ids.size()); + xcorr_matrix_max_peak_.resize(native_ids.size(), native_ids.size()); + xcorr_matrix_max_peak_sec_.resize(native_ids.size(), native_ids.size()); for (std::size_t i = 0; i < native_ids.size(); i++) { FeatureType fi = mrmfeature->getFeature(native_ids[i]); @@ -96,6 +104,9 @@ namespace OpenSwath fj->getIntensity(intensityj); // compute normalized cross correlation xcorr_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j)); + xcorr_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_matrix_max_peak_sec_.setValue(i, j, x->second); } } } @@ -104,6 +115,8 @@ namespace OpenSwath { std::vector intensityi, intensityj; xcorr_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); + xcorr_contrast_matrix_max_peak_.resize(native_ids_set1.size(), native_ids_set2.size()); + xcorr_contrast_matrix_max_peak_sec_.resize(native_ids_set1.size(), native_ids_set2.size()); for (std::size_t i = 0; i < native_ids_set1.size(); i++) { FeatureType fi = mrmfeature->getFeature(native_ids_set1[i]); @@ -116,6 +129,9 @@ namespace OpenSwath fj->getIntensity(intensityj); // compute normalized cross correlation xcorr_contrast_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_.getValue(i, j)); + xcorr_contrast_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_contrast_matrix_max_peak_sec_.setValue(i, j, x->second); } } } @@ -124,6 +140,8 @@ namespace OpenSwath { std::vector intensityi, intensityj; xcorr_precursor_matrix_.resize(precursor_ids.size(), precursor_ids.size()); + xcorr_precursor_matrix_max_peak_.resize(precursor_ids.size(), precursor_ids.size()); + xcorr_precursor_matrix_max_peak_sec_.resize(precursor_ids.size(), precursor_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); @@ -136,6 +154,9 @@ namespace OpenSwath fj->getIntensity(intensityj); // compute normalized cross correlation xcorr_precursor_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_.getValue(i, j)); + xcorr_precursor_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_precursor_matrix_max_peak_sec_.setValue(i, j, x->second); } } } @@ -144,6 +165,8 @@ namespace OpenSwath { std::vector intensityi, intensityj; xcorr_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size()); + xcorr_precursor_contrast_matrix_max_peak_.resize(precursor_ids.size(), native_ids.size()); + xcorr_precursor_contrast_matrix_max_peak_sec_.resize(precursor_ids.size(), native_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); @@ -156,6 +179,9 @@ namespace OpenSwath fj->getIntensity(intensityj); // compute normalized cross correlation xcorr_precursor_contrast_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_.getValue(i, j)); + xcorr_precursor_contrast_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_precursor_contrast_matrix_max_peak_sec_.setValue(i, j, x->second); } } } @@ -163,6 +189,8 @@ namespace OpenSwath void MRMScoring::initializeXCorrPrecursorContrastMatrix(const std::vector< std::vector< double > >& data_precursor, const std::vector< std::vector< double > >& data_fragments) { xcorr_precursor_contrast_matrix_.resize(data_precursor.size(), data_fragments.size()); + xcorr_precursor_contrast_matrix_max_peak_.resize(data_precursor.size(), data_fragments.size()); + xcorr_precursor_contrast_matrix_max_peak_sec_.resize(data_precursor.size(), data_fragments.size()); for (std::size_t i = 0; i < data_precursor.size(); i++) { std::vector< double > tmp1(data_precursor[i]); @@ -171,6 +199,9 @@ namespace OpenSwath // compute normalized cross correlation std::vector< double > tmp2(data_fragments[j]); xcorr_precursor_contrast_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(tmp1, tmp2, boost::numeric_cast(tmp1.size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_.getValue(i, j)); + xcorr_precursor_contrast_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_precursor_contrast_matrix_max_peak_sec_.setValue(i, j, x->second); #ifdef MRMSCORING_TESTING std::cout << " fill xcorr_precursor_contrast_matrix_ "<< tmp1.size() << " / " << tmp2.size() << " : " << xcorr_precursor_contrast_matrix_[i][j].data.size() << std::endl; #endif @@ -195,6 +226,8 @@ namespace OpenSwath } xcorr_precursor_combined_matrix_.resize(features.size(), features.size()); + xcorr_precursor_combined_matrix_max_peak_.resize(features.size(), features.size()); + xcorr_precursor_combined_matrix_max_peak_sec_.resize(features.size(), features.size()); for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; @@ -207,6 +240,9 @@ namespace OpenSwath fj->getIntensity(intensityj); // compute normalized cross correlation xcorr_precursor_combined_matrix_.setValue(i, j, Scoring::normalizedCrossCorrelation(intensityi, intensityj, boost::numeric_cast(intensityi.size()), 1)); + auto x = Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_.getValue(i, j)); + xcorr_precursor_combined_matrix_max_peak_.setValue(i, j, std::abs(x->first)); + xcorr_precursor_combined_matrix_max_peak_sec_.setValue(i, j, x->second); } } } @@ -219,15 +255,16 @@ namespace OpenSwath // return $deltascore_mean + $deltascore_stdev double MRMScoring::calcXcorrCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_max_peak_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_matrix_max_peak_.rows(); i++) { - for (std::size_t j = i; j < xcorr_matrix_.rows(); j++) + for (std::size_t j = i; j < xcorr_matrix_max_peak_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->first)); + //deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->first)); + deltas.push_back(xcorr_matrix_max_peak_.getValue(i,j)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->first) << std::endl; #endif @@ -246,15 +283,15 @@ namespace OpenSwath double MRMScoring::calcXcorrCoelutionWeightedScore( const std::vector& normalized_library_intensity) { - OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_max_peak_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); #ifdef MRMSCORING_TESTING double weights = 0; #endif double deltas{0}; - for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_matrix_max_peak_.rows(); i++) { - deltas += (std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, i))->first) + deltas += (xcorr_matrix_max_peak_.getValue(i, i)//std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, i))->first) * normalized_library_intensity[i] * normalized_library_intensity[i]); #ifdef MRMSCORING_TESTING @@ -262,10 +299,10 @@ namespace OpenSwath normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl; weights += normalized_library_intensity[i] * normalized_library_intensity[i]; #endif - for (std::size_t j = i + 1; j < xcorr_matrix_.rows(); j++) + for (std::size_t j = i + 1; j < xcorr_matrix_max_peak_.rows(); j++) { // first is the X value (RT), should be an int - deltas += (std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->first) + deltas += (xcorr_matrix_max_peak_.getValue(i, j)//std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->first) * normalized_library_intensity[i] * normalized_library_intensity[j] * 2); #ifdef MRMSCORING_TESTING @@ -286,13 +323,13 @@ namespace OpenSwath double MRMScoring::calcXcorrContrastCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_max_peak_.rows() > 0 && xcorr_contrast_matrix_max_peak_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector deltas; - for (auto e : xcorr_contrast_matrix_) + for (auto e : xcorr_contrast_matrix_max_peak_) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(e)->first)); + deltas.push_back(e); //(std::abs(Scoring::xcorrArrayGetMaxPeak(e)->first)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first) << std::endl; #endif @@ -309,21 +346,21 @@ namespace OpenSwath std::vector MRMScoring::calcSeparateXcorrContrastCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_max_peak_.rows() > 0 && xcorr_contrast_matrix_max_peak_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); - std::vector deltas; - for (std::size_t i = 0; i < xcorr_contrast_matrix_.rows(); i++) + std::vector deltas; + for (std::size_t i = 0; i < xcorr_contrast_matrix_max_peak_.rows(); i++) { double deltas_id = 0; - for (std::size_t j = 0; j < xcorr_contrast_matrix_.cols(); j++) + for (std::size_t j = 0; j < xcorr_contrast_matrix_max_peak_.cols(); j++) { // first is the X value (RT), should be an int - deltas_id += std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_.getValue(i,j))->first); + deltas_id += xcorr_contrast_matrix_max_peak_.getValue(i, j); #ifdef MRMSCORING_TESTING - std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first) << std::endl; + std::cout << "&&_xcoel append " << xcorr_contrast_matrix_max_peak_getValue(i, j) << std::endl; #endif } - deltas.push_back(deltas_id / xcorr_contrast_matrix_.cols()); + deltas.push_back(deltas_id / xcorr_contrast_matrix_max_peak_.cols()); } return deltas; @@ -331,15 +368,15 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_matrix_max_peak_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_precursor_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_precursor_matrix_max_peak_.rows(); i++) { - for (std::size_t j = i; j < xcorr_precursor_matrix_.rows(); j++) + for (std::size_t j = i; j < xcorr_precursor_matrix_max_peak_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_.getValue(i, j))->first)); + deltas.push_back(xcorr_precursor_matrix_max_peak_.getValue(i, j)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_[i][j])->first) << std::endl; #endif @@ -357,13 +394,13 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorContrastCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_.rows() > 0 && xcorr_precursor_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_max_peak_.rows() > 0 && xcorr_precursor_contrast_matrix_max_peak_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector deltas; - for (auto e : xcorr_precursor_contrast_matrix_) + for (auto e : xcorr_precursor_contrast_matrix_max_peak_) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(e)->first)); + deltas.push_back(e); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_contrast_matrix_[i][j])->first) << std::endl; #endif @@ -380,15 +417,15 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorCombinedCoelutionScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_max_peak_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); std::vector deltas; - for (std::size_t i = 0; i < xcorr_precursor_combined_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_precursor_combined_matrix_max_peak_.rows(); i++) { - for (std::size_t j = i; j < xcorr_precursor_combined_matrix_.rows(); j++) + for (std::size_t j = i; j < xcorr_precursor_combined_matrix_max_peak_.rows(); j++) { // first is the X value (RT), should be an int - deltas.push_back(std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_.getValue(i, j))->first)); + deltas.push_back(xcorr_precursor_combined_matrix_max_peak_.getValue(i, j)); #ifdef MRMSCORING_TESTING std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_[i][j])->first) << std::endl; #endif @@ -411,16 +448,16 @@ namespace OpenSwath /// double MRMScoring::calcXcorrShapeScore() { - OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_max_peak_sec_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); size_t element_number{0}; double intensities{0}; - for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_matrix_max_peak_sec_.rows(); i++) { - for (std::size_t j = i; j < xcorr_matrix_.rows(); j++) + for (std::size_t j = i; j < xcorr_matrix_max_peak_sec_.rows(); j++) { // second is the Y value (intensity) - intensities +=Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->second; + intensities += xcorr_matrix_max_peak_sec_.getValue(i, j); element_number++; } } @@ -430,24 +467,24 @@ namespace OpenSwath double MRMScoring::calcXcorrShapeWeightedScore( const std::vector& normalized_library_intensity) { - OPENSWATH_PRECONDITION(xcorr_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_matrix_max_peak_sec_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); // TODO (hroest) : check implementation // see _calc_weighted_xcorr_shape_score in MRM_pgroup.pm // -- they only multiply up the intensity once double intensities{0}; - for (std::size_t i = 0; i < xcorr_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_matrix_max_peak_sec_.rows(); i++) { - intensities += (Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, i))->second + intensities += (xcorr_matrix_max_peak_sec_.getValue(i, i) * normalized_library_intensity[i] * normalized_library_intensity[i]); #ifdef MRMSCORING_TESTING std::cout << "_xcorr_weighted " << i << " " << i << " " << Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][i])->second << " weight " << normalized_library_intensity[i] * normalized_library_intensity[i] << std::endl; #endif - for (std::size_t j = i + 1; j < xcorr_matrix_.rows(); j++) + for (std::size_t j = i + 1; j < xcorr_matrix_max_peak_sec_.rows(); j++) { - intensities += (Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_.getValue(i, j))->second + intensities += (xcorr_matrix_max_peak_sec_.getValue(i, j) * normalized_library_intensity[i] * normalized_library_intensity[j] * 2); #ifdef MRMSCORING_TESTING @@ -461,30 +498,30 @@ namespace OpenSwath double MRMScoring::calcXcorrContrastShapeScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_max_peak_sec_.rows() > 0 && xcorr_contrast_matrix_max_peak_sec_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); double intensities{0}; - for(auto e : xcorr_contrast_matrix_) + for(auto e : xcorr_contrast_matrix_max_peak_sec_) { - intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; + intensities += e; } return intensities; } std::vector MRMScoring::calcSeparateXcorrContrastShapeScore() { - OPENSWATH_PRECONDITION(xcorr_contrast_matrix_.rows() > 0 && xcorr_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_contrast_matrix_max_peak_sec_.rows() > 0 && xcorr_contrast_matrix_max_peak_sec_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); std::vector intensities; - for (std::size_t i = 0; i < xcorr_contrast_matrix_.rows(); i++) + for (std::size_t i = 0; i < xcorr_contrast_matrix_max_peak_sec_.rows(); i++) { double intensities_id = 0; - for (std::size_t j = 0; j < xcorr_contrast_matrix_.cols(); j++) + for (std::size_t j = 0; j < xcorr_contrast_matrix_max_peak_sec_.cols(); j++) { // second is the Y value (intensity) - intensities_id += Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_.getValue(i, j))->second; + intensities_id += xcorr_contrast_matrix_max_peak_sec_.getValue(i,j); } - intensities.push_back(intensities_id / xcorr_contrast_matrix_.cols()); + intensities.push_back(intensities_id / xcorr_contrast_matrix_max_peak_sec_.cols()); } return intensities; @@ -492,46 +529,46 @@ namespace OpenSwath double MRMScoring::calcXcorrPrecursorShapeScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_matrix_max_peak_sec_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); double intensities{0}; - for(auto e : xcorr_precursor_matrix_) + for(auto e : xcorr_precursor_matrix_max_peak_sec_) { - intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; + intensities += e; } //xcorr_precursor_matrix_ is a triangle matrix - size_t element_number = xcorr_precursor_matrix_.rows()*xcorr_precursor_matrix_.rows()/2 + (xcorr_precursor_matrix_.rows()+1)/2; + size_t element_number = xcorr_precursor_matrix_max_peak_sec_.rows()*xcorr_precursor_matrix_.rows()/2 + (xcorr_precursor_matrix_max_peak_sec_.rows()+1)/2; return intensities / element_number; } double MRMScoring::calcXcorrPrecursorContrastShapeScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_.rows() > 0 && xcorr_precursor_contrast_matrix_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_contrast_matrix_max_peak_sec_.rows() > 0 && xcorr_precursor_contrast_matrix_max_peak_sec_.cols() > 1, "Expect cross-correlation matrix of at least 1x2"); double intensities{0}; - for(auto e : xcorr_precursor_contrast_matrix_) + for(auto e : xcorr_precursor_contrast_matrix_max_peak_sec_) { - intensities += Scoring::xcorrArrayGetMaxPeak(e)->second; + intensities += e; } - return intensities / xcorr_precursor_contrast_matrix_.size(); + return intensities / xcorr_precursor_contrast_matrix_max_peak_sec_.size(); } double MRMScoring::calcXcorrPrecursorCombinedShapeScore() { - OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); + OPENSWATH_PRECONDITION(xcorr_precursor_combined_matrix_max_peak_sec_.rows() > 1, "Expect cross-correlation matrix of at least 2x2"); double intensities{0}; - for(size_t i = 0; i < xcorr_precursor_combined_matrix_.rows(); i++) + for(size_t i = 0; i < xcorr_precursor_combined_matrix_max_peak_sec_.rows(); i++) { - for(size_t j = i; j < xcorr_precursor_combined_matrix_.cols(); j++) + for(size_t j = i; j < xcorr_precursor_combined_matrix_max_peak_sec_.cols(); j++) { - intensities += Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_.getValue(i,j))->second; + intensities += xcorr_precursor_combined_matrix_max_peak_sec_.getValue(i, j); } } //xcorr_precursor-combined_matrix_ is a triangle matrix - size_t element_number = xcorr_precursor_combined_matrix_.rows()*xcorr_precursor_combined_matrix_.rows()/2 + (xcorr_precursor_combined_matrix_.rows()+1)/2; + size_t element_number = xcorr_precursor_combined_matrix_max_peak_sec_.rows()*xcorr_precursor_combined_matrix_max_peak_sec_.rows()/2 + (xcorr_precursor_combined_matrix_max_peak_sec_.rows()+1)/2; return intensities / element_number; } diff --git a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp index 4993bab6b89..3a1fc169b8e 100644 --- a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp +++ b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp @@ -284,6 +284,24 @@ BOOST_AUTO_TEST_CASE(initializeXCorrPrecursorContrastMatrix) TEST_EQUAL(mrmscore.getXCorrPrecursorContrastMatrix().rows(), 3) TEST_EQUAL(mrmscore.getXCorrPrecursorContrastMatrix().cols(), 2) + + std::vector sum_matrix; + for(auto e : mrmscore.getXCorrPrecursorContrastMatrix()) + { + double sum{0}; + for(size_t i = 0; i < e.data.size(); ++i) + { + sum += abs(e.data[i].second); + } + sum_matrix.push_back(sum); + } + + TEST_REAL_SIMILAR(sum_matrix[0], 3.40949220) + TEST_REAL_SIMILAR(sum_matrix[1], 6.19794611) + TEST_REAL_SIMILAR(sum_matrix[2], 3.68912454) + TEST_REAL_SIMILAR(sum_matrix[3], 6.60757921) + TEST_REAL_SIMILAR(sum_matrix[4], 0) + TEST_REAL_SIMILAR(sum_matrix[5], 0) } END_SECTION @@ -302,9 +320,63 @@ BOOST_AUTO_TEST_CASE(initializeXCorrPrecursorCombinedMatrix) TEST_EQUAL(mrmscore.getXCorrPrecursorCombinedMatrix().rows(), 5) TEST_EQUAL(mrmscore.getXCorrPrecursorCombinedMatrix().cols(), 5) + + std::vector sum_matrix; + for(auto e : mrmscore.getXCorrPrecursorCombinedMatrix()) + { + double sum{0}; + for(size_t i = 0; i < e.data.size(); ++i) + { + sum += abs(e.data[i].second); + } + sum_matrix.push_back(sum); + } + + TEST_REAL_SIMILAR(sum_matrix[0], 5.86440677) + TEST_REAL_SIMILAR(sum_matrix[1], 6.05410398) + TEST_REAL_SIMILAR(sum_matrix[2], 0) + TEST_REAL_SIMILAR(sum_matrix[3], 3.40949220) + TEST_REAL_SIMILAR(sum_matrix[4], 6.19794611) + TEST_REAL_SIMILAR(sum_matrix[5], 6.05410398) + TEST_REAL_SIMILAR(sum_matrix[6], 6.30751744) + TEST_REAL_SIMILAR(sum_matrix[7], 0) + TEST_REAL_SIMILAR(sum_matrix[8], 3.68912454) + TEST_REAL_SIMILAR(sum_matrix[9], 6.60757921) + TEST_REAL_SIMILAR(sum_matrix[10], 0) + TEST_REAL_SIMILAR(sum_matrix[11], 0) + TEST_REAL_SIMILAR(sum_matrix[12], 0) + TEST_REAL_SIMILAR(sum_matrix[13], 0) + TEST_REAL_SIMILAR(sum_matrix[14], 0) + TEST_REAL_SIMILAR(sum_matrix[15], 3.40949220) + TEST_REAL_SIMILAR(sum_matrix[16], 3.68912454) + TEST_REAL_SIMILAR(sum_matrix[17], 0) + TEST_REAL_SIMILAR(sum_matrix[18], 3.13711983) + TEST_REAL_SIMILAR(sum_matrix[19], 3.57832717) + TEST_REAL_SIMILAR(sum_matrix[20], 6.19794611) + TEST_REAL_SIMILAR(sum_matrix[21], 6.60757921) + TEST_REAL_SIMILAR(sum_matrix[22], 0) + TEST_REAL_SIMILAR(sum_matrix[23], 3.57832717) + TEST_REAL_SIMILAR(sum_matrix[25], 0) } END_SECTION +/*BOOST_AUTO_TEST_CASE(initializeXCorrPrecursorMatrix) +{ + MockMRMFeature * imrmfeature = new MockMRMFeature(); + MRMScoring mrmscore; + + std::vector precursor_ids; + fill_mock_objects(imrmfeature, precursor_ids); + + //initialize the XCorr vector + mrmscore.initializeXCorrPrecursorMatrix(imrmfeature, precursor_ids); + delete imrmfeature; + + TEST_EQUAL(mrmscore.getXCorrPrecursorMatrix().rows(), 3) + TEST_EQUAL(mrmscore.getXCorrPrecurso().cols(), 2) +} +END_SECTION*/ + BOOST_AUTO_TEST_CASE(initializeXCorrContrastMatrix) { MockMRMFeature * imrmfeature = new MockMRMFeature(); @@ -691,6 +763,12 @@ BOOST_AUTO_TEST_CASE(initializeMIPrecursorContrastMatrix) TEST_EQUAL(mrmscore.getMIPrecursorContrastMatrix().rows(), 3) TEST_EQUAL(mrmscore.getMIPrecursorContrastMatrix().cols(), 2) + double sum{0}; + for(auto e : mrmscore.getMIPrecursorContrastMatrix()) + { + sum += e; + } + TEST_REAL_SIMILAR(sum, 12.01954465) } END_SECTION @@ -709,6 +787,13 @@ BOOST_AUTO_TEST_CASE(initializeMIPrecursorCombinedMatrix) TEST_EQUAL(mrmscore.getMIPrecursorCombinedMatrix().rows(), 5) TEST_EQUAL(mrmscore.getMIPrecursorCombinedMatrix().cols(), 5) + + double sum{0}; + for(auto e : mrmscore.getMIPrecursorCombinedMatrix()) + { + sum += e; + } + TEST_REAL_SIMILAR(sum, 48.98726953) } END_SECTION From 20f4f50438387d37bf296f30ef60593620b6c9c5 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Thu, 4 Nov 2021 09:15:51 +0100 Subject: [PATCH 06/11] delete some comments --- .../OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h | 4 ++-- .../source/OPENSWATHALGO/ALGO/Scoring.cpp | 23 ------------------- 2 files changed, 2 insertions(+), 25 deletions(-) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h index 6a10d97f7ab..31d239a8504 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h @@ -112,7 +112,7 @@ namespace OpenSwath /** @name Scores */ //@{ /// Initialize the scoring object and building the cross-correlation matrix - void initializeXCorrMatrix(const std::vector< std::vector< double > >& data); //andere matrix? + void initializeXCorrMatrix(const std::vector< std::vector< double > >& data); /// Initialize the scoring object and building the cross-correlation matrix void initializeXCorrMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids); @@ -295,7 +295,7 @@ namespace OpenSwath OpenMS::Matrix xcorr_precursor_combined_matrix_max_peak_; OpenMS::Matrix xcorr_precursor_combined_matrix_max_peak_sec_; /// the precomputed mutual information matrix - //std::vector< std::vector > mi_matrix_; + OpenMS::Matrix mi_matrix_; /// the precomputed contrast mutual information matrix OpenMS::Matrix mi_contrast_matrix_; diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp index 03c07d18c54..4401d51b066 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp @@ -270,29 +270,6 @@ namespace OpenSwath::Scoring std::vector computeRank(const std::vector& v_temp) { - /*std::vector ranks{}; - ranks.resize(v_temp.size()); - std::iota(ranks.begin(), ranks.end(), 0); - std::sort(ranks.begin(), ranks.end(), - [&v_temp](unsigned int i, unsigned int j) { return v_temp[i] < v_temp[j]; }); - unsigned int tmp = 0; - for (unsigned int i = 1; i < ranks.size(); ++i) - { - if (v_temp[tmp] == v_temp[ranks[i]]) - { - tmp = ranks[i]; - ranks[i] = ranks[i-1]; - } - else - { - tmp = ranks[i]; - ranks[i] = i; - } - } - - return ranks; - }*/ - std::vector > v_sort(v_temp.size()); for (unsigned int i = 0; i < v_sort.size(); ++i) { From e88b2315a9866ecf873f2c4689033ce7f3461e15 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Thu, 4 Nov 2021 11:41:15 +0100 Subject: [PATCH 07/11] new branch --- src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 64a6084fdd8..1609b4522d7 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -808,7 +808,6 @@ namespace OpenSwath for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; - //mi_precursor_combined_matrix_[i].resize(features.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < features.size(); j++) From d25139f690f214900fbf0dbc580ba19607c5bdb4 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Thu, 4 Nov 2021 17:34:04 +0100 Subject: [PATCH 08/11] new computeRanks, only in MIcombinedMatrix --- .../OpenMS/OPENSWATHALGO/ALGO/Scoring.h | 4 ++- .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 8 +++-- .../source/OPENSWATHALGO/ALGO/Scoring.cpp | 33 +++++++++++++++++++ 3 files changed, 42 insertions(+), 3 deletions(-) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h index 49610e0d03c..695d8618988 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h @@ -131,10 +131,12 @@ namespace OpenSwath OPENSWATHALGO_DLLAPI void normalize_sum(double x[], unsigned int n); // Compute rank of vector elements - OPENSWATHALGO_DLLAPI std::vector computeRank(const std::vector& w); + OPENSWATHALGO_DLLAPI std::vector computeRank(const std::vector& v_temp); + OPENSWATHALGO_DLLAPI void computeRank2(const std::vector& v, std::vector& ranks); // Estimate rank-transformed mutual information between two vectors of data points OPENSWATHALGO_DLLAPI double rankedMutualInformation(std::vector& data1, std::vector& data2); + OPENSWATHALGO_DLLAPI double rankedMutualInformation2(std::vector& data1, std::vector& data2); //@} diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 1609b4522d7..f3f1340e2b7 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -44,6 +44,7 @@ #include // for isnan #include + namespace OpenSwath { @@ -803,22 +804,25 @@ namespace OpenSwath FeatureType fj = mrmfeature->getFeature(native_ids[j]); features.push_back(fj); } - + std::vector rank_vec1, rank_vec2; mi_precursor_combined_matrix_.resize(features.size(), features.size()); for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank2(intensityi, rank_vec1); for (std::size_t j = 0; j < features.size(); j++) { FeatureType fj = features[j]; intensityj.clear(); fj->getIntensity(intensityj); + Scoring::computeRank2(intensityj, rank_vec2); // compute ranked mutual information - mi_precursor_combined_matrix_.setValue(i ,j, Scoring::rankedMutualInformation(intensityi, intensityj)); + mi_precursor_combined_matrix_.setValue(i , j, Scoring::rankedMutualInformation2(rank_vec1, rank_vec2)); } } + } double MRMScoring::calcMIScore() diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp index 4401d51b066..8463c935e96 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp @@ -291,6 +291,28 @@ namespace OpenSwath::Scoring } return result; } + void computeRank2(const std::vector& v_temp, std::vector& ranks) + { + std::vector > v_sort(v_temp.size()); + + for (unsigned int i = 0; i < v_sort.size(); ++i) { + v_sort[i] = std::make_pair(v_temp[i], i); + } + + std::sort(v_sort.begin(), v_sort.end()); + + std::pair rank; + + ranks.resize(v_temp.size()); + for (unsigned int i = 0; i < v_sort.size(); ++i) + { + if (v_sort[i].first != rank.first) + { + rank = std::make_pair(v_sort[i].first, i); + } + ranks[v_sort[i].second] = rank.second; + } + } double rankedMutualInformation(std::vector& data1, std::vector& data2) { @@ -305,6 +327,17 @@ namespace OpenSwath::Scoring double result = calcMutualInformation(arr_int_data1, arr_int_data2, int_data1.size()); + return result; + } + double rankedMutualInformation2(std::vector& data1, std::vector& data2) + { + OPENSWATH_PRECONDITION(data1.size() != 0 && data1.size() == data2.size(), "Both data vectors need to have the same length"); + + unsigned int* arr_int_data1 = &data1[0]; + unsigned int* arr_int_data2 = &data2[0]; + + double result = calcMutualInformation(arr_int_data1, arr_int_data2, data1.size()); + return result; } } //namespace OpenMS // namespace Scoring From abc6bddee25999564bd29105f1704e0b846c5ec7 Mon Sep 17 00:00:00 2001 From: "khue.nm" Date: Thu, 4 Nov 2021 19:34:25 +0100 Subject: [PATCH 09/11] Moved computeRank outside of inner loop --- .../OpenMS/OPENSWATHALGO/ALGO/Scoring.h | 6 ++++ .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 6 +++- .../source/OPENSWATHALGO/ALGO/Scoring.cpp | 35 +++++++++++++++++++ 3 files changed, 46 insertions(+), 1 deletion(-) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h index 49610e0d03c..13f1ddff926 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h @@ -133,9 +133,15 @@ namespace OpenSwath // Compute rank of vector elements OPENSWATHALGO_DLLAPI std::vector computeRank(const std::vector& w); + /// Compute rank of vector elements + OPENSWATHALGO_DLLAPI void computeRank_inplace(const std::vector& w, std::vector& rank_vector); + // Estimate rank-transformed mutual information between two vectors of data points OPENSWATHALGO_DLLAPI double rankedMutualInformation(std::vector& data1, std::vector& data2); + // Estimate rank-transformed mutual information between two vectors of rank + OPENSWATHALGO_DLLAPI double preCalcRankedMutualInformation(std::vector& rank_vector1, std::vector& rank_vector2); + //@} } diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 1609b4522d7..894a8083c44 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -804,19 +804,23 @@ namespace OpenSwath features.push_back(fj); } + std::vector rank_vector1, rank_vector2; + mi_precursor_combined_matrix_.resize(features.size(), features.size()); for (std::size_t i = 0; i < features.size(); i++) { FeatureType fi = features[i]; intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank_inplace(intensityi, rank_vector1); for (std::size_t j = 0; j < features.size(); j++) { FeatureType fj = features[j]; intensityj.clear(); fj->getIntensity(intensityj); + Scoring::computeRank_inplace(intensityj, rank_vector2); // compute ranked mutual information - mi_precursor_combined_matrix_.setValue(i ,j, Scoring::rankedMutualInformation(intensityi, intensityj)); + mi_precursor_combined_matrix_.setValue(i ,j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp index 4401d51b066..188fb072bec 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp @@ -292,6 +292,29 @@ namespace OpenSwath::Scoring return result; } + void computeRank_inplace(const std::vector& v_temp, std::vector& result) + { + std::vector > v_sort(v_temp.size()); + + for (unsigned int i = 0; i < v_sort.size(); ++i) { + v_sort[i] = std::make_pair(v_temp[i], i); + } + + std::sort(v_sort.begin(), v_sort.end()); + + std::pair rank; + result.resize(v_temp.size()); + + for (unsigned int i = 0; i < v_sort.size(); ++i) + { + if (v_sort[i].first != rank.first) + { + rank = std::make_pair(v_sort[i].first, i); + } + result[v_sort[i].second] = rank.second; + } + } + double rankedMutualInformation(std::vector& data1, std::vector& data2) { OPENSWATH_PRECONDITION(data1.size() != 0 && data1.size() == data2.size(), "Both data vectors need to have the same length"); @@ -307,4 +330,16 @@ namespace OpenSwath::Scoring return result; } + + double preCalcRankedMutualInformation(std::vector& rank_vector1, std::vector& rank_vector2) + { + OPENSWATH_PRECONDITION(rank_vector1.size() != 0 && rank_vector1.size() == rank_vector1.size(), "Both data vectors need to have the same length"); + + unsigned int* arr_int_data1 = &rank_vector1[0]; + unsigned int* arr_int_data2 = &rank_vector2[0]; + + double result = calcMutualInformation(arr_int_data1, arr_int_data2, rank_vector1.size()); + + return result; + } } //namespace OpenMS // namespace Scoring From d3ee09708dc7e190c24e85672f53879d4511a772 Mon Sep 17 00:00:00 2001 From: Vincent Musch Date: Fri, 5 Nov 2021 15:57:44 +0100 Subject: [PATCH 10/11] little fix --- src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp index 8463c935e96..9e087e5f236 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp @@ -231,7 +231,8 @@ namespace OpenSwath::Scoring // sigma_1 * sigma_2 * n denominator = sqrt(sqsum1 * sqsum2); } - + //avoids division in the for loop + denominator = 1/denominator; XCorrArrayType result; result.data.reserve( (size_t)std::ceil((2*maxdelay + 1) / lag)); int cnt = 0; @@ -257,12 +258,12 @@ namespace OpenSwath::Scoring if (denominator > 0) { - result.data.push_back(std::make_pair(delay, sxy/denominator)); + result.data.emplace_back(delay, sxy*denominator); } else { // e.g. if all datapoints are zero - result.data.push_back(std::make_pair(delay, 0)); + result.data.emplace_back(delay, 0); } } return result; From ce8937eba1029dda85ca6c3457ae3b21dfbdf660 Mon Sep 17 00:00:00 2001 From: "khue.nm" Date: Wed, 10 Nov 2021 11:04:31 +0100 Subject: [PATCH 11/11] precompute rank for all intialize MI matrix funcs --- .../OpenMS/OPENSWATHALGO/ALGO/Scoring.h | 5 +-- .../source/OPENSWATHALGO/ALGO/MRMScoring.cpp | 31 ++++++++++++++----- .../source/OPENSWATHALGO/ALGO/Scoring.cpp | 31 +++---------------- 3 files changed, 28 insertions(+), 39 deletions(-) diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h index 13f1ddff926..6a3c5a66ccc 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h @@ -131,10 +131,7 @@ namespace OpenSwath OPENSWATHALGO_DLLAPI void normalize_sum(double x[], unsigned int n); // Compute rank of vector elements - OPENSWATHALGO_DLLAPI std::vector computeRank(const std::vector& w); - - /// Compute rank of vector elements - OPENSWATHALGO_DLLAPI void computeRank_inplace(const std::vector& w, std::vector& rank_vector); + OPENSWATHALGO_DLLAPI void computeRank(const std::vector& v_temp, std::vector& result); // Estimate rank-transformed mutual information between two vectors of data points OPENSWATHALGO_DLLAPI double rankedMutualInformation(std::vector& data1, std::vector& data2); diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp index 894a8083c44..b0bfe4a4aae 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/MRMScoring.cpp @@ -708,6 +708,7 @@ namespace OpenSwath void MRMScoring::initializeMIMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids) { std::vector intensityi, intensityj; + std::vector rank_vector1, rank_vector2; mi_matrix_.resize(native_ids.size(),native_ids.size()); for (std::size_t i = 0; i < native_ids.size(); i++) { @@ -715,13 +716,15 @@ namespace OpenSwath intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank(intensityi, rank_vector1); for (std::size_t j = i; j < native_ids.size(); j++) { FeatureType fj = mrmfeature->getFeature(native_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_matrix_.setValue(i,j,Scoring::rankedMutualInformation(intensityi, intensityj)); + mi_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } @@ -729,6 +732,7 @@ namespace OpenSwath void MRMScoring::initializeMIContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& native_ids_set1, const std::vector& native_ids_set2) { std::vector intensityi, intensityj; + std::vector rank_vector1, rank_vector2; mi_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size()); for (std::size_t i = 0; i < native_ids_set1.size(); i++) { @@ -736,13 +740,16 @@ namespace OpenSwath //mi_contrast_matrix_[i].resize(native_ids_set2.size()); intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank(intensityi, rank_vector1); for (std::size_t j = 0; j < native_ids_set2.size(); j++) { FeatureType fj = mrmfeature->getFeature(native_ids_set2[j]); intensityj.clear(); fj->getIntensity(intensityj); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj)); + // is this symmetric? + mi_contrast_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } @@ -750,19 +757,22 @@ namespace OpenSwath void MRMScoring::initializeMIPrecursorMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids) { std::vector intensityi, intensityj; + std::vector rank_vector1, rank_vector2; mi_precursor_matrix_.resize(precursor_ids.size(),precursor_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank(intensityi, rank_vector1); for (std::size_t j = i; j < precursor_ids.size(); j++) { FeatureType fj = mrmfeature->getPrecursorFeature(precursor_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_precursor_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj)); + mi_precursor_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } @@ -770,6 +780,7 @@ namespace OpenSwath void MRMScoring::initializeMIPrecursorContrastMatrix(OpenSwath::IMRMFeature* mrmfeature, const std::vector& precursor_ids, const std::vector& native_ids) { std::vector intensityi, intensityj; + std::vector rank_vector1, rank_vector2; mi_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size()); for (std::size_t i = 0; i < precursor_ids.size(); i++) { @@ -777,13 +788,15 @@ namespace OpenSwath //mi_precursor_contrast_matrix_[i].resize(native_ids.size()); intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank(intensityi, rank_vector1); for (std::size_t j = 0; j < native_ids.size(); j++) { FeatureType fj = mrmfeature->getFeature(native_ids[j]); intensityj.clear(); fj->getIntensity(intensityj); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_precursor_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj)); + mi_precursor_contrast_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } @@ -812,15 +825,17 @@ namespace OpenSwath FeatureType fi = features[i]; intensityi.clear(); fi->getIntensity(intensityi); - Scoring::computeRank_inplace(intensityi, rank_vector1); - for (std::size_t j = 0; j < features.size(); j++) + Scoring::computeRank(intensityi, rank_vector1); + for (std::size_t j = i; j < features.size(); j++) { FeatureType fj = features[j]; intensityj.clear(); fj->getIntensity(intensityj); - Scoring::computeRank_inplace(intensityj, rank_vector2); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_precursor_combined_matrix_.setValue(i ,j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); + double curr_mi_score = Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2); + mi_precursor_combined_matrix_.setValue(i ,j, curr_mi_score); + mi_precursor_combined_matrix_.setValue(j ,i, curr_mi_score); } } } diff --git a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp index 188fb072bec..98755e476fa 100644 --- a/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp +++ b/src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp @@ -268,31 +268,7 @@ namespace OpenSwath::Scoring return result; } - std::vector computeRank(const std::vector& v_temp) - { - std::vector > v_sort(v_temp.size()); - - for (unsigned int i = 0; i < v_sort.size(); ++i) { - v_sort[i] = std::make_pair(v_temp[i], i); - } - - std::sort(v_sort.begin(), v_sort.end()); - - std::pair rank; - std::vector result(v_temp.size()); - - for (unsigned int i = 0; i < v_sort.size(); ++i) - { - if (v_sort[i].first != rank.first) - { - rank = std::make_pair(v_sort[i].first, i); - } - result[v_sort[i].second] = rank.second; - } - return result; - } - - void computeRank_inplace(const std::vector& v_temp, std::vector& result) + void computeRank(const std::vector& v_temp, std::vector& result) { std::vector > v_sort(v_temp.size()); @@ -320,8 +296,9 @@ namespace OpenSwath::Scoring OPENSWATH_PRECONDITION(data1.size() != 0 && data1.size() == data2.size(), "Both data vectors need to have the same length"); // rank the data - std::vector int_data1 = computeRank(data1); - std::vector int_data2 = computeRank(data2); + std::vector int_data1, int_data2; + computeRank(data1,int_data1); + computeRank(data2,int_data2); unsigned int* arr_int_data1 = &int_data1[0]; unsigned int* arr_int_data2 = &int_data2[0];