diff --git a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h index 073b1b3e65f..31d239a8504 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/MRMScoring.h @@ -36,6 +36,7 @@ #include +#include #include #include "OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h" @@ -78,7 +79,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; @@ -218,29 +219,29 @@ 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, 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); @@ -264,36 +265,50 @@ 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_; /// 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 e65d58242c1..6a3c5a66ccc 100644 --- a/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h +++ b/src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h @@ -115,11 +115,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, 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, const int& maxdelay, const 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); @@ -131,11 +131,14 @@ 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 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); + // 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/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 a5ec9c024ad..b0bfe4a4aae 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 @@ -53,16 +54,21 @@ namespace OpenSwath void MRMScoring::initializeXCorrMatrix(const std::vector< std::vector< double > >& data) { - xcorr_matrix_.resize(data.size()); + 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++) { - 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)); + 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); } } } @@ -85,22 +91,24 @@ 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()); + 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++) { - String native_id = native_ids[i]; - FeatureType fi = mrmfeature->getFeature(native_id); - xcorr_matrix_[i].resize(native_ids.size()); + 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 - 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)); + 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); } } } @@ -108,22 +116,24 @@ 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()); + 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++) { - String native_id = native_ids_set1[i]; - FeatureType fi = mrmfeature->getFeature(native_id); - xcorr_contrast_matrix_[i].resize(native_ids_set2.size()); + 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 - 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)); + 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); } } } @@ -131,22 +141,24 @@ 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()); + 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++) { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - xcorr_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 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)); + 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); } } } @@ -154,38 +166,44 @@ 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()); + 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++) - { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - xcorr_precursor_contrast_matrix_[i].resize(native_ids.size()); + { + 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 - 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)); + 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); } } } 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()); + 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++) - { - 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)); + 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 @@ -199,23 +217,22 @@ 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); } - xcorr_precursor_combined_matrix_.resize(features.size()); + 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]; - xcorr_precursor_combined_matrix_[i].resize(features.size()); intensityi.clear(); fi->getIntensity(intensityi); for (std::size_t j = 0; j < features.size(); j++) @@ -224,7 +241,10 @@ 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)); + 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); } } } @@ -237,15 +257,16 @@ 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_max_peak_.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_max_peak_.rows(); i++) { - for (std::size_t j = i; j < xcorr_matrix_.size(); 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_[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 @@ -264,16 +285,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_max_peak_.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_max_peak_.rows(); i++) { - deltas.push_back( - std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[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 @@ -281,11 +301,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_max_peak_.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 += (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 @@ -301,24 +320,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_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_.size(); i++) + for (auto e : xcorr_contrast_matrix_max_peak_) { - 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(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; + std::cout << "&&_xcoel append " << std::abs(Scoring::xcorrArrayGetMaxPeak(xcorr_contrast_matrix_[i][j])->first) << std::endl; #endif - } } OpenSwath::mean_and_stddev msc; @@ -332,21 +348,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_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_.size(); 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_[0].size(); 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_[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_[0].size()); + deltas.push_back(deltas_id / xcorr_contrast_matrix_max_peak_.cols()); } return deltas; @@ -354,15 +370,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_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_.size(); 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_.size(); 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_[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 @@ -380,19 +396,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_max_peak_.rows() > 0 && xcorr_precursor_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_precursor_contrast_matrix_.size(); i++) + for (auto e : xcorr_precursor_contrast_matrix_max_peak_) { - 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(e); #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; @@ -406,21 +419,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_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_.size(); 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_.size(); 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_[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 } } - OpenSwath::mean_and_stddev msc; msc = std::for_each(deltas.begin(), deltas.end(), msc); double deltas_mean = msc.mean(); @@ -438,45 +450,43 @@ namespace OpenSwath /// double MRMScoring::calcXcorrShapeScore() { - OPENSWATH_PRECONDITION(xcorr_matrix_.size() > 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"); - 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_max_peak_sec_.rows(); i++) { - for (std::size_t j = i; j < xcorr_matrix_.size(); j++) + for (std::size_t j = i; j < xcorr_matrix_max_peak_sec_.rows(); j++) { // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second); + intensities += xcorr_matrix_max_peak_sec_.getValue(i, j); + element_number++; } } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + 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_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 - 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_max_peak_sec_.rows(); i++) { - intensities.push_back( - Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[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_.size(); j++) + for (std::size_t j = i + 1; j < xcorr_matrix_max_peak_sec_.rows(); j++) { - intensities.push_back( - Scoring::xcorrArrayGetMaxPeak(xcorr_matrix_[i][j])->second + intensities += (xcorr_matrix_max_peak_sec_.getValue(i, j) * normalized_library_intensity[i] * normalized_library_intensity[j] * 2); #ifdef MRMSCORING_TESTING @@ -485,41 +495,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_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_.size(); i++) + double intensities{0}; + for(auto e : xcorr_contrast_matrix_max_peak_sec_) { - 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 += e; } - 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_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_.size(); 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_[0].size(); 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_[i][j])->second; + intensities_id += xcorr_contrast_matrix_max_peak_sec_.getValue(i,j); } - intensities.push_back(intensities_id / xcorr_contrast_matrix_[0].size()); + intensities.push_back(intensities_id / xcorr_contrast_matrix_max_peak_sec_.cols()); } return intensities; @@ -527,56 +531,47 @@ 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_max_peak_sec_.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++) + double intensities{0}; + for(auto e : xcorr_precursor_matrix_max_peak_sec_) { - for (std::size_t j = i; j < xcorr_precursor_matrix_.size(); j++) - { - // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_matrix_[i][j])->second); - } + intensities += e; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + //xcorr_precursor_matrix_ is a triangle matrix + 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_.size() > 0 && xcorr_precursor_contrast_matrix_[0].size() > 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"); - 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_max_peak_sec_) { - 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 += e; } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + return intensities / xcorr_precursor_contrast_matrix_max_peak_sec_.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_max_peak_sec_.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(size_t i = 0; i < xcorr_precursor_combined_matrix_max_peak_sec_.rows(); i++) { - for (std::size_t j = i; j < xcorr_precursor_combined_matrix_.size(); j++) + for(size_t j = i; j < xcorr_precursor_combined_matrix_max_peak_sec_.cols(); j++) { - // second is the Y value (intensity) - intensities.push_back(Scoring::xcorrArrayGetMaxPeak(xcorr_precursor_combined_matrix_[i][j])->second); + intensities += xcorr_precursor_combined_matrix_max_peak_sec_.getValue(i, j); } } - OpenSwath::mean_and_stddev msc; - msc = std::for_each(intensities.begin(), intensities.end(), msc); - return msc.mean(); + //xcorr_precursor-combined_matrix_ is a triangle matrix + 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; } void MRMScoring::calcLibraryScore(OpenSwath::IMRMFeature* mrmfeature, const std::vector& transitions, @@ -690,91 +685,94 @@ 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_; } - 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()); + 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++) { - String native_id = native_ids[i]; - FeatureType fi = mrmfeature->getFeature(native_id); - mi_matrix_[i].resize(native_ids.size()); + FeatureType fi = mrmfeature->getFeature(native_ids[i]); + intensityi.clear(); fi->getIntensity(intensityi); + Scoring::computeRank(intensityi, rank_vector1); 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); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } - 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()); + 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++) - { - String native_id = native_ids_set1[i]; - FeatureType fi = mrmfeature->getFeature(native_id); - mi_contrast_matrix_[i].resize(native_ids_set2.size()); + { + FeatureType fi = mrmfeature->getFeature(native_ids_set1[i]); + //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++) { - 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); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_contrast_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + // is this symmetric? + mi_contrast_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } - 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()); + 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++) { - 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); + Scoring::computeRank(intensityi, rank_vector1); 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); + Scoring::computeRank(intensityj, rank_vector2); // compute ranked mutual information - mi_precursor_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_precursor_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } @@ -782,22 +780,23 @@ 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()); + 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++) - { - String precursor_id = precursor_ids[i]; - FeatureType fi = mrmfeature->getPrecursorFeature(precursor_id); - mi_precursor_contrast_matrix_[i].resize(native_ids.size()); + { + FeatureType fi = mrmfeature->getPrecursorFeature(precursor_ids[i]); + //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++) { - String native_id = native_ids[j]; - FeatureType fj = mrmfeature->getFeature(native_id); + 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_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + mi_precursor_contrast_matrix_.setValue(i, j, Scoring::preCalcRankedMutualInformation(rank_vector1, rank_vector2)); } } } @@ -808,149 +807,136 @@ 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); } - mi_precursor_combined_matrix_.resize(features.size()); + 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]; - mi_precursor_combined_matrix_[i].resize(features.size()); intensityi.clear(); fi->getIntensity(intensityi); - 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(intensityj, rank_vector2); // compute ranked mutual information - mi_precursor_combined_matrix_[i][j] = Scoring::rankedMutualInformation(intensityi, intensityj); + 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); } } } 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..98755e476fa 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"); @@ -268,7 +268,7 @@ namespace OpenSwath::Scoring return result; } - std::vector computeRank(const std::vector& v_temp) + void computeRank(const std::vector& v_temp, std::vector& result) { std::vector > v_sort(v_temp.size()); @@ -279,7 +279,7 @@ namespace OpenSwath::Scoring std::sort(v_sort.begin(), v_sort.end()); std::pair rank; - std::vector result(v_temp.size()); + result.resize(v_temp.size()); for (unsigned int i = 0; i < v_sort.size(); ++i) { @@ -289,7 +289,6 @@ namespace OpenSwath::Scoring } result[v_sort[i].second] = rank.second; } - return result; } double rankedMutualInformation(std::vector& data1, std::vector& data2) @@ -297,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]; @@ -307,4 +307,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 diff --git a/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp b/src/tests/class_tests/openswathalgo/MRMScoring_test.cpp index a2dceee4a8a..3a1fc169b8e 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 {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}; 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,26 @@ 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) + + 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 @@ -316,11 +318,65 @@ 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) + + 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(); @@ -333,13 +389,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 +404,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 +738,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 +761,14 @@ 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) + double sum{0}; + for(auto e : mrmscore.getMIPrecursorContrastMatrix()) + { + sum += e; + } + TEST_REAL_SIMILAR(sum, 12.01954465) } END_SECTION @@ -723,8 +785,15 @@ 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) + + double sum{0}; + for(auto e : mrmscore.getMIPrecursorCombinedMatrix()) + { + sum += e; + } + TEST_REAL_SIMILAR(sum, 48.98726953) } END_SECTION @@ -740,10 +809,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);