Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Swath4 #122

Open
wants to merge 12 commits into
base: develop
Choose a base branch
from
Open

Swath4 #122

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
e88b231
new branch
Nov 4, 2021
d25139f
new computeRanks, only in MIcombinedMatrix
Nov 4, 2021
d3ee097
little fix
Nov 5, 2021
dc65a58
aktualisiern
Nov 12, 2021
88fb258
Merge branch 'develop' of https://github.com/OpenMS/OpenMS into swath4
Nov 12, 2021
2aa9ff7
new computeRank
Nov 12, 2021
3e59247
delete old function
Nov 12, 2021
5330cc6
little fixes
Nov 16, 2021
6c90292
fix
Nov 18, 2021
16902f0
new MutualInformation fkt
Nov 25, 2021
de64c00
.h definition
Nov 26, 2021
d4299dd
tests
Nov 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions src/openms/include/OpenMS/ANALYSIS/OPENSWATH/MRMTransitionGroupPicker.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,12 +445,13 @@ namespace OpenMS
double transition_total_mi = 0;
if (compute_total_mi_)
{
std::vector<unsigned int> rank_vec1, rank_vec2;
std::vector<double> chrom_vect_id, chrom_vect_det;
for (typename SpectrumT::const_iterator it = chromatogram.begin(); it != chromatogram.end(); it++)
{
chrom_vect_id.push_back(it->getIntensity());
}

OpenSwath::Scoring::computeRank(chrom_vect_id, rank_vec2);
// compute baseline mutual information
int transition_total_mi_norm = 0;
for (Size m = 0; m < transition_group.getTransitions().size(); m++)
Expand All @@ -463,7 +464,8 @@ namespace OpenMS
{
chrom_vect_det.push_back(it->getIntensity());
}
transition_total_mi += OpenSwath::Scoring::rankedMutualInformation(chrom_vect_det, chrom_vect_id);
OpenSwath::Scoring::computeRank(chrom_vect_det, rank_vec1);
transition_total_mi += OpenSwath::Scoring::rankedMutualInformation(rank_vec1, rank_vec2);
transition_total_mi_norm++;
}
}
Expand Down
26 changes: 20 additions & 6 deletions src/openms/source/ANALYSIS/OPENSWATH/MRMScoring.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -709,19 +709,22 @@ namespace OpenSwath
{
std::vector<double> intensityi, intensityj;
mi_matrix_.resize(native_ids.size(),native_ids.size());
std::vector<unsigned int> rank_vec1, rank_vec2;
for (std::size_t i = 0; i < native_ids.size(); i++)
{
FeatureType fi = mrmfeature->getFeature(native_ids[i]);

intensityi.clear();
fi->getIntensity(intensityi);
Scoring::computeRank(intensityi, rank_vec1);
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_vec2);
// compute ranked mutual information
mi_matrix_.setValue(i,j,Scoring::rankedMutualInformation(intensityi, intensityj));
mi_matrix_.setValue(i,j,Scoring::rankedMutualInformation(rank_vec1, rank_vec2));
}
}
}
Expand All @@ -730,19 +733,22 @@ namespace OpenSwath
{
std::vector<double> intensityi, intensityj;
mi_contrast_matrix_.resize(native_ids_set1.size(), native_ids_set2.size());
std::vector<unsigned int> rank_vec1, rank_vec2;
for (std::size_t i = 0; i < native_ids_set1.size(); i++)
{
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_vec1);
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_vec2);
// compute ranked mutual information
mi_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj));
mi_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(rank_vec1, rank_vec2));
}
}
}
Expand All @@ -751,18 +757,21 @@ namespace OpenSwath
{
std::vector<double> intensityi, intensityj;
mi_precursor_matrix_.resize(precursor_ids.size(),precursor_ids.size());
std::vector<unsigned int> rank_vec1, rank_vec2;
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_vec1);
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_vec2);
// compute ranked mutual information
mi_precursor_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj));
mi_precursor_matrix_.setValue(i, j, Scoring::rankedMutualInformation(rank_vec1, rank_vec2));
}
}
}
Expand All @@ -771,19 +780,22 @@ namespace OpenSwath
{
std::vector<double> intensityi, intensityj;
mi_precursor_contrast_matrix_.resize(precursor_ids.size(), native_ids.size());
std::vector<unsigned int> rank_vec1, rank_vec2;
for (std::size_t i = 0; i < precursor_ids.size(); i++)
{
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_vec1);
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_vec2);
// compute ranked mutual information
mi_precursor_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(intensityi, intensityj));
mi_precursor_contrast_matrix_.setValue(i, j, Scoring::rankedMutualInformation(rank_vec1, rank_vec2));
}
}
}
Expand All @@ -803,20 +815,22 @@ namespace OpenSwath
FeatureType fj = mrmfeature->getFeature(native_ids[j]);
features.push_back(fj);
}

std::vector<unsigned int> 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::computeRank(intensityi, rank_vec1);
for (std::size_t j = 0; j < features.size(); j++)
{
FeatureType fj = features[j];
intensityj.clear();
fj->getIntensity(intensityj);
Scoring::computeRank(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::rankedMutualInformation(rank_vec1, rank_vec2));
}
}
}
Expand Down
21 changes: 19 additions & 2 deletions src/openswathalgo/include/OpenMS/OPENSWATHALGO/ALGO/Scoring.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,17 @@ namespace OpenSwath
iterator end() {return data.end();}
const_iterator end() const {return data.end();}
};

struct jpstate
{
std::vector<double> jointProbabilityVector;
int numJointStates;
std::vector<double> firstProbabilityVector;
int numFirstStates;
std::vector<double> secondProbabilityVector;
int numSecondStates;
std::vector<unsigned int> jointPositionVector;
};
//@}

/** @name Helper functions */
Expand Down Expand Up @@ -131,10 +142,16 @@ namespace OpenSwath
OPENSWATHALGO_DLLAPI void normalize_sum(double x[], unsigned int n);

// Compute rank of vector elements
OPENSWATHALGO_DLLAPI std::vector<unsigned int> computeRank(const std::vector<double>& w);
//OPENSWATHALGO_DLLAPI std::vector<unsigned int> computeRank(const std::vector<double>& v_temp);
OPENSWATHALGO_DLLAPI void computeRank(const std::vector<double>& v, std::vector<unsigned int>& ranks);

// Estimate rank-transformed mutual information between two vectors of data points
OPENSWATHALGO_DLLAPI double rankedMutualInformation(std::vector<double>& data1, std::vector<double>& data2);
//OPENSWATHALGO_DLLAPI double rankedMutualInformation(std::vector<double>& data1, std::vector<double>& data2);
OPENSWATHALGO_DLLAPI double rankedMutualInformation(std::vector<unsigned int>& data1, std::vector<unsigned int>& data2);

OPENSWATHALGO_DLLAPI unsigned int maxElem(const std::vector<unsigned int>& arr);
OPENSWATHALGO_DLLAPI jpstate calcJointProbability(const std::vector<unsigned int>& firstVector,const std::vector<unsigned int>& secondVector,const int& vectorLength);
OPENSWATHALGO_DLLAPI double mutualInformation(jpstate& state,const std::vector<unsigned int>& firstVector,const std::vector<unsigned int>& secondVector);

//@}

Expand Down
144 changes: 116 additions & 28 deletions src/openswathalgo/source/OPENSWATHALGO/ALGO/Scoring.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
// $Authors: Hannes Roest $
// --------------------------------------------------------------------------

#include <iostream>
#include <OpenMS/OPENSWATHALGO/ALGO/Scoring.h>
#include <OpenMS/OPENSWATHALGO/Macros.h>
#include <cmath>
Expand All @@ -55,7 +56,7 @@ namespace OpenSwath::Scoring
return;
}
auto inverse_sum = 1 / sumx; // precompute inverse since division is expensive!
for (int i = 0; i < n; ++i)
for (unsigned int i = 0; i < n; ++i)
{
x[i] *= inverse_sum;
}
Expand Down Expand Up @@ -153,30 +154,33 @@ namespace OpenSwath::Scoring
{
stdev = 1; // all data is equal
}
stdev = 1/stdev;
for (std::size_t i = 0; i < data.size(); i++)
{
data[i] = (data[i] - mean) / stdev;
data[i] = (data[i] - mean) * stdev;
}
}

XCorrArrayType normalizedCrossCorrelation(std::vector<double>& data1,
std::vector<double>& data2, int maxdelay, int lag = 1)
std::vector<double>& 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");

// normalize the data
standardize_data(data1);
standardize_data(data2);
XCorrArrayType result = calculateCrossCorrelation(data1, data2, maxdelay, lag);
for (XCorrArrayType::iterator it = result.begin(); it != result.end(); ++it)

double d = 1.0 / data1.size();
for(auto& e : result)
{
it->second = it->second / data1.size();
e.second *= d;
}
return result;
}

XCorrArrayType calculateCrossCorrelation(const std::vector<double>& data1,
const std::vector<double>& data2, int maxdelay, int lag)
const std::vector<double>& 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");

Expand Down Expand Up @@ -232,7 +236,8 @@ namespace OpenSwath::Scoring
// sigma_1 * sigma_2 * n
denominator = sqrt(sqsum1 * sqsum2);
}
denominator = 1/denominator; // inverse denominator for faster calculation
//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;
Expand Down Expand Up @@ -269,43 +274,126 @@ namespace OpenSwath::Scoring
return result;
}

std::vector<unsigned int> computeRank(const std::vector<double>& v_temp)
void computeRank(const std::vector<double>& v_temp, std::vector<unsigned int>& ranks_out)
{
std::vector<unsigned int> 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]; });
ranks_out.clear();
ranks_out.resize(v_temp.size());
double x = 0;
unsigned int y = 0;
for(unsigned int i = 0; i < ranks.size();++i)
{
if(v_temp[ranks[i]] != x)
{
x = v_temp[ranks[i]];
y = i;
}
ranks_out[ranks[i]] = y;
}
}

unsigned int maxElem(const std::vector<unsigned int>& arr)
{
unsigned int max = arr[0];
for(auto e : arr)
{
if(e > max) max = e;
}
return max+1;
}



jpstate calcJointProbability(const std::vector<unsigned int>& firstVector,const std::vector<unsigned int>& secondVector,const int& vectorLength)
{
std::vector<std::pair<float, unsigned int> > v_sort(v_temp.size());
jpstate state;
double length = 1.0 / vectorLength;
unsigned int firstNumStates = maxElem(firstVector);
unsigned int secondNumStates = maxElem(secondVector);
unsigned int jointNumStates = firstNumStates * secondNumStates;

std::vector<unsigned int> firstStateCounts(firstNumStates, 0);
std::vector<unsigned int> secondStateCounts(secondNumStates, 0);
std::vector<unsigned int> jointStateCounts(jointNumStates, 0);
std::vector<unsigned int> jointPosition(firstNumStates, 0);

std::vector<double> firstStateProbs(firstNumStates, 0.0);
std::vector<double> secondStateProbs(secondNumStates, 0.0);
std::vector<double> jointStateProbs(jointNumStates, 0.0);

for(int i = 0; i < vectorLength; i++)
{
firstStateCounts[firstVector[i]] += 1;
secondStateCounts[secondVector[i]] += 1;
jointPosition[i] = secondVector[i] * firstNumStates + firstVector[i];
jointStateCounts[jointPosition[i]] += 1;
}

for (unsigned int i = 0; i < firstNumStates; i++) {
firstStateProbs[i] = firstStateCounts[i] * length;
}

for (unsigned int i = 0; i < secondNumStates; i++) {
secondStateProbs[i] = secondStateCounts[i] * length;
}

for (unsigned int i = 0; i < v_sort.size(); ++i) {
v_sort[i] = std::make_pair(v_temp[i], i);
for (unsigned int i = 0; i < jointNumStates; i++) {
jointStateProbs[i] = jointStateCounts[i] * length;
}

std::sort(v_sort.begin(), v_sort.end());
state.jointPositionVector = jointPosition;
state.jointProbabilityVector = jointStateProbs;
state.numJointStates = jointNumStates;
state.firstProbabilityVector = firstStateProbs;
state.numFirstStates = firstNumStates;
state.secondProbabilityVector = secondStateProbs;
state.numSecondStates = secondNumStates;

std::pair<double, unsigned int> rank;
std::vector<unsigned int> result(v_temp.size());
return state;
}

double mutualInformation(jpstate& state,const std::vector<unsigned int>& firstVector,const std::vector<unsigned int>& secondVector)
{
double mutualInformation = 0.0;
//int firstIndex,secondIndex;

for (unsigned int i = 0; i < v_sort.size(); ++i)
/*
** I(X;Y) = \sum_x \sum_y p(x,y) * \log (p(x,y)/p(x)p(y))
*/
for (unsigned int i = 0; i < firstVector.size(); i++)
{
if (v_sort[i].first != rank.first)

int j = state.jointPositionVector[i];
if(state.jointProbabilityVector[j] != 0)
{
rank = std::make_pair(v_sort[i].first, i);
/*double division is probably more stable than multiplying two small numbers together
** mutualInformation += state.jointProbabilityVector[i] * log(state.jointProbabilityVector[i] / (state.firstProbabilityVector[firstIndex] * state.secondProbabilityVector[secondIndex]));
*/
mutualInformation += state.jointProbabilityVector[j] *
log2(state.jointProbabilityVector[j] / state.firstProbabilityVector[firstVector[i]] /
state.secondProbabilityVector[secondVector[i]]);
state.jointProbabilityVector[j] = 0;
}
result[v_sort[i].second] = rank.second;
}
return result;
return mutualInformation;
}

double rankedMutualInformation(std::vector<double>& data1, std::vector<double>& data2)
double rankedMutualInformation(std::vector<unsigned int>& data1, std::vector<unsigned int>& data2)
{
OPENSWATH_PRECONDITION(data1.size() != 0 && data1.size() == data2.size(), "Both data vectors need to have the same length");

// rank the data
std::vector<unsigned int> int_data1 = computeRank(data1);
std::vector<unsigned int> int_data2 = computeRank(data2);

unsigned int* arr_int_data1 = &int_data1[0];
unsigned int* arr_int_data2 = &int_data2[0];

double result = calcMutualInformation(arr_int_data1, arr_int_data2, int_data1.size());
jpstate state = calcJointProbability(data1, data2, data1.size());

double result = mutualInformation(state, data1, data2);
/*
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
Loading