48 #include <boost/math/special_functions/bessel.hpp> 57 #pragma clang diagnostic push 58 #pragma clang diagnostic ignored "-Wfloat-equal" 59 #pragma clang diagnostic ignored "-Wconversion" 60 #pragma clang diagnostic ignored "-Wshorten-64-to-32" 68 template <
typename PeakType>
88 typedef std::multimap<UInt, BoxElement>
Box;
102 reference_(nullptr), trans_intens_(nullptr)
108 reference_(reference)
110 trans_intens_ =
new std::vector<float>(reference_->size(), 0.0);
116 delete (trans_intens_);
121 delete (trans_intens_);
122 trans_intens_ = NULL;
130 return reference_->getRT();
136 return (*reference_)[i].getMZ();
142 return (*reference_)[i].getIntensity();
148 return (*trans_intens_)[i];
154 (*trans_intens_)[i] =
intens;
160 return trans_intens_->size();
177 inline typename MSSpectrum::const_iterator
MZBegin(
const double mz)
const 179 return reference_->MZBegin(mz);
184 inline typename MSSpectrum::const_iterator
MZEnd(
const double mz)
const 186 return reference_->MZEnd(mz);
191 inline typename MSSpectrum::const_iterator
end()
const 193 return reference_->end();
198 inline typename MSSpectrum::const_iterator
begin()
const 200 return reference_->begin();
255 const double ampl_cutoff,
const bool check_PPMs);
266 const UInt RT_votes_cutoff,
const Int front_bound = -1,
const Int end_bound = -1);
284 inline double getLinearInterpolation(
const typename MSSpectrum::const_iterator& left_iter,
const double mz_pos,
const typename MSSpectrum::const_iterator& right_iter)
286 return left_iter->getIntensity() + (right_iter->getIntensity() - left_iter->getIntensity()) / (right_iter->getMZ() - left_iter->getMZ()) * (mz_pos - left_iter->getMZ());
295 inline double getLinearInterpolation(
const double mz_a,
const double intens_a,
const double mz_pos,
const double mz_b,
const double intens_b)
297 return intens_a + (intens_b - intens_a) / (mz_b - mz_a) * (mz_pos - mz_a);
340 const double seed_mz,
const UInt c,
const double ampl_cutoff);
349 const double seed_mz,
const UInt c,
const double ampl_cutoff);
360 const UInt c,
const UInt scan_index,
const bool check_PPMs,
const double transintens,
const double prev_score);
370 const UInt c,
const UInt scan_index,
const bool check_PPMs,
const double transintens,
const double prev_score);
414 const double intens,
const double rt,
const UInt MZ_begin,
const UInt MZ_end);
429 const UInt scan_index,
const UInt c,
const bool check_PPMs);
436 const UInt scan_index,
const UInt c,
const bool check_PPMs);
450 double old_frac_mass = c_mass - (
Int)(c_mass);
452 double new_frac_mass = new_mass - (
Int)(new_mass);
454 if (new_frac_mass - old_frac_mass > 0.5)
459 if (new_frac_mass - old_frac_mass < -0.5)
470 inline double getPPMs_(
const double mass_a,
const double mass_b)
const 472 return fabs(mass_a - mass_b) / (0.5 * (mass_a + mass_b)) * 1e6;
494 template <
typename PeakType>
500 template <
typename PeakType>
506 template <
typename PeakType>
512 template <
typename PeakType>
518 template <
typename PeakType>
521 tmp_boxes_ =
new std::vector<std::multimap<double, Box> >(1);
530 template <
typename PeakType>
537 tmp_boxes_ =
new std::vector<std::multimap<double, Box> >(max_charge);
538 if (max_scan_size <= 0)
549 psi_.reserve(to_reserve);
550 prod_.reserve(to_reserve);
551 xs_.reserve(to_reserve);
556 template <
typename PeakType>
562 template <
typename PeakType>
565 Int spec_size((
Int)c_ref.size());
569 double value, T_boundary_left, T_boundary_right, old, c_diff, current, old_pos, my_local_MZ, my_local_lambda, origin, c_mz;
571 for (
Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
574 old = 0; old_pos = (my_local_pos -
from_max_to_left_ - 1 >= 0) ? c_ref[my_local_pos - from_max_to_left_ - 1].getMZ() : c_ref[0].getMZ() -
min_spacing_;
579 for (
Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
581 if (current_conv_pos >= spec_size)
587 c_mz = c_ref[current_conv_pos].getMZ();
588 c_diff = c_mz + origin;
591 current = c_diff > T_boundary_left && c_diff <= T_boundary_right ?
IsotopeWavelet::getValueByLambda(my_local_lambda, c_diff * charge + 1.) * c_ref[current_conv_pos].getIntensity() : 0;
593 value += 0.5 * (current + old) * (c_mz - old_pos);
601 c_trans[my_local_pos].setIntensity(value);
605 template <
typename PeakType>
608 Int spec_size((
Int)c_ref.size());
612 double value, T_boundary_left, T_boundary_right, c_diff, current, my_local_MZ, my_local_lambda, origin, c_mz;
614 for (
Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
623 for (
Int current_conv_pos = std::max(0, my_local_pos -
from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
625 if (current_conv_pos >= spec_size)
630 c_mz = c_ref[current_conv_pos].getMZ();
631 c_diff = c_mz + origin;
634 current = c_diff > T_boundary_left && c_diff <= T_boundary_right ?
IsotopeWavelet::getValueByLambda(my_local_lambda, c_diff * charge + 1.) * c_ref[current_conv_pos].getIntensity() : 0;
639 c_trans[my_local_pos].setIntensity(value);
643 template <
typename PeakType>
648 Int wavelet_length = 0, quarter_length = 0;
653 typename MSSpectrum::const_iterator start_iter, end_iter;
657 start_iter = c_ref.
MZEnd(c_ref[i].getMZ());
658 end_iter = c_ref.
MZBegin(c_ref[i].getMZ() + c_mz_cutoff);
659 wavelet_length = std::max((
SignedSize) wavelet_length, distance(start_iter, end_iter) + 1);
661 quarter_length = std::max((
SignedSize) quarter_length, distance(end_iter, start_iter) + 1);
673 if (wavelet_length > (
Int) c_ref.size())
675 std::cout <<
"Warning: the extremal length of the wavelet is larger (" << wavelet_length <<
") than the number of data points (" << c_ref.size() <<
"). This might (!) severely affect the transform." << std::endl;
676 std::cout <<
"Minimal spacing: " <<
min_spacing_ << std::endl;
677 std::cout <<
"Warning/Error generated at scan with RT " << c_ref.
getRT() <<
"." << std::endl;
685 template <
typename PeakType>
689 for (
UInt c_conv_pos = 1; c_conv_pos < c_ref.size(); ++c_conv_pos)
695 template <
typename PeakType>
697 const MSSpectrum& ref,
const UInt scan_index,
const UInt c,
const double ampl_cutoff,
const bool check_PPMs)
699 Size scan_size(candidates.size());
701 typename MSSpectrum::const_iterator iter_start, iter_end, iter_p, seed_iter, iter2;
702 double mz_cutoff, seed_mz, c_av_intens = 0, c_score = 0, c_sd_intens = 0, threshold = 0, help_mz, share, share_pos, bwd, fwd;
706 diffed[0].setIntensity(0); diffed[scan_size - 1].setIntensity(0);
708 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 709 std::stringstream stream;
710 stream <<
"diffed_" << ref.
getRT() <<
"_" << c + 1 <<
".trans\0";
711 std::ofstream ofile(stream.str().c_str());
716 for (
UInt i = 0; i < scan_size - 2; ++i)
718 share = candidates[i + 1].getIntensity(); share_pos = candidates[i + 1].getMZ();
719 bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
720 fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
722 if (!(bwd >= 0 && fwd <= 0) || share > ref[i + 1].getIntensity())
724 diffed[i + 1].setIntensity(0);
727 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 728 ofile << diffed[i + 1].getMZ() <<
"\t" << diffed[i + 1].getIntensity() << std::endl;
734 for (
UInt i = 0; i < scan_size - 2; ++i)
736 share = candidates[i + 1].getIntensity(); share_pos = candidates[i + 1].getMZ();
737 bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
738 fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
740 if (!(bwd >= 0 && fwd <= 0))
742 diffed[i + 1].setIntensity(0);
745 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 746 ofile << diffed[i + 1].getMZ() <<
"\t" << diffed[i + 1].getIntensity() << std::endl;
750 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 759 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 760 std::stringstream stream2;
761 stream2 <<
"sorted_cpu_" << candidates.
getRT() <<
"_" << c + 1 <<
".trans\0";
762 std::ofstream ofile2(stream2.str().c_str());
763 for (iter = c_sorted_candidate.end() - 1; iter != c_sorted_candidate.begin(); --iter)
765 ofile2 << iter->getMZ() <<
"\t" << iter->getIntensity() << std::endl;
770 std::vector<UInt> processed(scan_size, 0);
780 threshold = ampl_cutoff * c_sd_intens + c_av_intens;
783 for (iter = c_sorted_candidate.end() - 1; iter != c_sorted_candidate.begin(); --iter)
785 if (iter->getIntensity() <= 0)
790 seed_mz = iter->getMZ();
791 seed_iter = ref.
MZBegin(seed_mz);
793 if (seed_iter == ref.end() || processed[distance(ref.begin(), seed_iter)])
804 iter_end = ref.
MZEnd(seed_iter, seed_mz + mz_cutoff / (c + 1.), ref.end());
805 if (iter_end == ref.end())
810 MZ_start = distance(ref.begin(), iter_start);
811 MZ_end = distance(ref.begin(), iter_end);
813 memset(&(processed[MZ_start]), 1,
sizeof(
UInt) * (MZ_end - MZ_start + 1));
817 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 818 if (trunc(seed_mz) == 874)
819 std::cout << seed_mz <<
"\t" << c_score << std::endl;
822 if (c_score <= 0 && c_score != -1000)
832 iter2 = candidates.
MZBegin(help_mz);
833 if (iter2 == candidates.end() || iter2 == candidates.begin())
841 if (iter2 != candidates.end())
848 iter2 = candidates.
MZBegin(help_mz);
849 if (iter2 == candidates.end() || iter2 == candidates.begin())
857 if (iter2 != candidates.end())
867 template <
typename PeakType>
869 const UInt peak_cutoff,
const double seed_mz,
const UInt c,
const double ampl_cutoff)
871 double c_score = 0, c_val;
872 typename MSSpectrum::const_iterator c_left_iter2, c_right_iter2;
873 Int signal_size((
Int)candidate.size());
878 Int p_h_ind = 1, end = 4 * (peak_cutoff - 1) - 1;
880 std::vector<double> positions(end);
881 for (
Int i = 0; i < end; ++i)
886 double l_score = 0, mid_val = 0;
887 Int start_index = distance(candidate.begin(), candidate.
MZBegin(positions[0])) - 1;
888 for (
Int v = 1; v <= end; ++v, ++p_h_ind)
892 if (start_index < signal_size - 1)
897 while (candidate[start_index].getMZ() < positions[v - 1]);
899 if (start_index <= 0 || start_index >= signal_size - 1)
904 c_left_iter2 = candidate.begin() + start_index - 1;
905 c_right_iter2 = c_left_iter2 + 1;
907 c_val = c_left_iter2->getIntensity() + (c_right_iter2->getIntensity() - c_left_iter2->getIntensity()) / (c_right_iter2->getMZ() - c_left_iter2->getMZ()) * (positions[v - 1] - c_left_iter2->getMZ());
909 if (v == (
int)(ceil(end / 2.)))
915 if (p_h_ind % 2 == 1)
918 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 919 if (trunc(seed_mz) == 874)
920 std::cout << -c_val << std::endl;
926 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 927 if (trunc(seed_mz) == 874)
928 std::cout << c_val << std::endl;
933 start_index = distance(candidate.begin(), c_left_iter2);
936 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 937 std::ofstream ofile_score(
"scores.dat", ios::app);
938 std::ofstream ofile_check_score(
"check_scores.dat", ios::app);
940 ofile_check_score.close();
943 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 944 if (trunc(seed_mz) == 874)
945 std::cout <<
"final_score: " << seed_mz <<
"\t" << c_score <<
"\t l_score: " << l_score <<
"\t" << c_score - l_score - mid_val <<
"\t" << c_score - mid_val <<
"\t" << ampl_cutoff << std::endl;
948 if (c_score - mid_val <= 0)
953 if (c_score - mid_val <= ampl_cutoff)
958 if (l_score <= 0 || c_score - l_score - mid_val <= 0)
966 template <
typename PeakType>
968 const UInt peak_cutoff,
const double seed_mz,
const UInt c,
const double ampl_cutoff)
970 double c_score = 0, c_val;
971 typename MSSpectrum::const_iterator c_left_iter2, c_right_iter2;
975 Int p_h_ind = 1, end = 4 * (peak_cutoff - 1) - 1;
977 std::vector<double> positions(end);
978 for (
Int i = 0; i < end; ++i)
983 double l_score = 0, mid_val = 0;
984 Int start_index = distance(candidate.
begin(), candidate.
MZBegin(positions[0])) - 1;
985 for (
Int v = 1; v <= end; ++v, ++p_h_ind)
989 if (start_index < signal_size - 1)
994 while (candidate.
getMZ(start_index) < positions[v - 1]);
996 if (start_index <= 0 || start_index >= signal_size - 1)
1001 c_left_iter2 = candidate.
begin() + start_index - 1;
1002 c_right_iter2 = c_left_iter2 + 1;
1005 if (v == (
int)(ceil(end / 2.)))
1011 if (p_h_ind % 2 == 1)
1020 start_index = distance(candidate.
begin(), c_left_iter2);
1023 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1024 std::ofstream ofile_score(
"scores.dat", ios::app);
1025 std::ofstream ofile_check_score(
"check_scores.dat", ios::app);
1026 ofile_score << c_check_point <<
"\t" << c_score << std::endl;
1027 ofile_score.close();
1028 ofile_check_score.close();
1031 if (l_score <= 0 || c_score - l_score - mid_val <= 0 || c_score - mid_val <= ampl_cutoff)
1039 template <
typename PeakType>
1043 typename std::multimap<double, Box>::iterator iter;
1044 typename Box::iterator box_iter;
1045 std::vector<BoxElement> final_box;
1046 double c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1047 double virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1049 typename std::pair<double, double> c_extend;
1053 Box& c_box = iter->second;
1054 av_score = 0; av_mz = 0; av_intens = 0; av_abs_intens = 0; count = 0;
1055 virtual_av_mz = 0; virtual_av_intens = 0; virtual_av_abs_intens = 0; virtual_count = 0;
1058 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1060 if (box_iter->second.score == 0)
1065 c_mz = box_iter->second.mz;
1066 virtual_av_intens += box_iter->second.intens;
1067 virtual_av_abs_intens += fabs(box_iter->second.intens);
1068 virtual_av_mz += c_mz * fabs(box_iter->second.intens);
1073 c_mz = box_iter->second.mz;
1074 av_score += box_iter->second.score;
1075 av_intens += box_iter->second.intens;
1076 av_abs_intens += fabs(box_iter->second.intens);
1077 av_mz += c_mz * fabs(box_iter->second.intens);
1084 av_intens = virtual_av_intens / virtual_count;
1086 av_mz = virtual_av_mz / virtual_av_abs_intens;
1092 av_mz /= av_abs_intens;
1096 c_box_element.
mz = av_mz;
1097 c_box_element.
c =
c;
1098 c_box_element.
score = av_score;
1099 c_box_element.
intens = av_intens;
1101 c_box_element.
RT = c_box.begin()->second.RT;
1102 final_box.push_back(c_box_element);
1105 Size num_o_feature = final_box.size();
1106 if (num_o_feature == 0)
1113 std::vector<double> bwd_diffs(num_o_feature, 0);
1116 for (
Size i = 1; i < num_o_feature; ++i)
1118 bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].
mz - final_box[i - 1].
mz);
1121 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1122 std::ofstream ofile_bwd(
"bwd_cpu.dat");
1123 for (
Size i = 0; i < num_o_feature; ++i)
1125 ofile_bwd << final_box[i].mz <<
"\t" << bwd_diffs[i] << std::endl;
1131 for (
Size i = 0; i < num_o_feature - 1; ++i)
1133 while (i < num_o_feature - 2)
1135 if (final_box[i].
score > 0 || final_box[i].
score == -1000)
1140 if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
1150 template <
typename PeakType>
1153 double av_intens = 0;
1154 for (
UInt i = 0; i < scan.size(); ++i)
1156 if (scan[i].getIntensity() >= 0)
1158 av_intens += scan[i].getIntensity();
1161 return av_intens / (
double)scan.size();
1164 template <
typename PeakType>
1168 for (
UInt i = 0; i < scan.size(); ++i)
1170 if (scan[i].getIntensity() >= 0)
1172 intens = scan[i].getIntensity();
1176 return sqrt(res / (
double)(scan.size() - 1));
1179 template <
typename PeakType>
1182 std::vector<double> diffs(scan.size() - 1, 0);
1183 for (
UInt i = 0; i < scan.size() - 1; ++i)
1185 diffs[i] = scan[i + 1].getMZ() - scan[i].getMZ();
1188 sort(diffs.begin(), diffs.end());
1189 double av_MZ_spacing = 0;
1190 for (
UInt i = 0; i < diffs.size() / 2; ++i)
1192 av_MZ_spacing += diffs[i];
1195 return av_MZ_spacing / (diffs.size() / 2);
1198 template <
typename PeakType>
1201 double av_intens = 0;
1202 for (
UInt i = 0; i < scan.
size(); ++i)
1212 template <
typename PeakType>
1216 for (
UInt i = 0; i < scan.
size(); ++i)
1224 return sqrt(res / (
double)(scan.
size() - 1));
1227 template <
typename PeakType>
1233 typename std::multimap<double, Box>::iterator upper_iter(
open_boxes_.upper_bound(mz));
1234 typename std::multimap<double, Box>::iterator lower_iter(
open_boxes_.lower_bound(mz));
1239 if (mz != lower_iter->first && lower_iter !=
open_boxes_.begin())
1245 typename std::multimap<double, Box>::iterator insert_iter;
1246 bool create_new_box =
true;
1254 if (fabs((--lower_iter)->first - mz) < dist_constraint)
1256 create_new_box =
false;
1257 insert_iter = lower_iter;
1262 create_new_box =
true;
1267 if (upper_iter ==
open_boxes_.end() && fabs(lower_iter->first - mz) < dist_constraint)
1269 insert_iter = lower_iter;
1270 create_new_box =
false;
1274 create_new_box =
true;
1284 double dist_lower = fabs(lower_iter->first - mz);
1285 double dist_upper = fabs(upper_iter->first - mz);
1286 dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1287 dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1289 if (dist_lower >= dist_constraint && dist_upper >= dist_constraint)
1291 create_new_box =
true;
1295 insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1296 create_new_box =
false;
1305 if (create_new_box ==
false)
1307 std::pair<UInt, BoxElement> help2(scan, element);
1308 insert_iter->second.insert(help2);
1311 Box replacement(insert_iter->second);
1315 double c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
1316 c_mz /= ((
double) insert_iter->second.size());
1320 std::pair<double, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1325 std::pair<UInt, BoxElement> help2(scan, element);
1326 std::multimap<UInt, BoxElement> help3;
1327 help3.insert(help2);
1328 std::pair<double, std::multimap<UInt, BoxElement> > help4(mz, help3);
1333 template <
typename PeakType>
1339 std::multimap<double, Box>& tmp_box(
tmp_boxes_->at(c));
1340 typename std::multimap<double, Box>::iterator upper_iter(tmp_box.upper_bound(mz));
1341 typename std::multimap<double, Box>::iterator lower_iter(tmp_box.lower_bound(mz));
1343 if (lower_iter != tmp_box.end())
1346 if (mz != lower_iter->first && lower_iter != tmp_box.begin())
1352 typename std::multimap<double, Box>::iterator insert_iter;
1353 bool create_new_box =
true;
1354 if (lower_iter == tmp_box.end())
1359 if (!tmp_box.empty())
1361 if (fabs((--lower_iter)->first - mz) < dist_constraint)
1363 create_new_box =
false;
1364 insert_iter = lower_iter;
1369 create_new_box =
true;
1374 if (upper_iter == tmp_box.end() && fabs(lower_iter->first - mz) < dist_constraint)
1376 insert_iter = lower_iter;
1377 create_new_box =
false;
1381 create_new_box =
true;
1386 if (upper_iter != tmp_box.end() && lower_iter != tmp_box.end())
1389 double dist_lower = fabs(lower_iter->first - mz);
1390 double dist_upper = fabs(upper_iter->first - mz);
1391 dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1392 dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1394 if (dist_lower >= dist_constraint && dist_upper >= dist_constraint)
1396 create_new_box =
true;
1400 insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1401 create_new_box =
false;
1409 if (create_new_box ==
false)
1411 std::pair<UInt, BoxElement> help2(scan, element);
1412 insert_iter->second.insert(help2);
1415 Box replacement(insert_iter->second);
1419 double c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
1420 c_mz /= ((
double) insert_iter->second.size());
1423 tmp_box.erase(insert_iter);
1424 std::pair<double, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1425 tmp_box.insert(help3);
1429 std::pair<UInt, BoxElement> help2(scan, element);
1430 std::multimap<UInt, BoxElement> help3;
1431 help3.insert(help2);
1433 std::pair<double, std::multimap<UInt, BoxElement> > help4(mz, help3);
1434 tmp_box.insert(help4);
1438 template <
typename PeakType>
1440 const UInt RT_votes_cutoff,
const Int front_bound,
const Int end_bound)
1442 typename std::multimap<double, Box>::iterator iter, iter2;
1444 if ((
Int)scan_index == end_bound && end_bound != (
Int)map.
size() - 1)
1448 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1449 std::cout <<
"LOW THREAD insert in end_box " << iter->first << std::endl;
1450 typename Box::iterator dings;
1451 for (dings = iter->second.begin(); dings != iter->second.end(); ++dings)
1452 std::cout << map[dings->first].getRT() <<
"\t" << dings->second.c + 1 << std::endl;
1462 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1463 if (front_bound > 0)
1465 std::cout <<
"HIGH THREAD open box. " << iter->first <<
"\t current scan index " << scan_index <<
"\t" << ((iter->second.begin()))->first <<
"\t of last scan " << map.
size() - 1 <<
"\t" << front_bound << std::endl;
1470 UInt lastScan = (--(iter->second.end()))->first;
1471 if (scan_index - lastScan > RT_interleave + 1 || scan_index == map.
size() - 1)
1473 if (iter->second.begin()->first - front_bound <= RT_interleave + 1 && front_bound > 0)
1477 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1478 std::cout <<
"HIGH THREAD insert in front_box " << iter->first << std::endl;
1491 if (iter->second.size() >= RT_votes_cutoff)
1507 template <
typename PeakType>
1510 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1511 std::cout <<
"**** CHECKING FOR BOX EXTENSIONS ****" << std::endl;
1515 typename Box::const_iterator iter;
1516 std::vector<double> elution_profile(box.size());
1518 for (iter = box.begin(); iter != box.end(); ++iter, ++index)
1520 for (
Size i = iter->second.MZ_begin; i != iter->second.MZ_end; ++i)
1522 elution_profile[index] += map[iter->second.RT_index][i].getIntensity();
1524 elution_profile[index] /= iter->second.MZ_end - iter->second.MZ_begin + 1.;
1528 Int max_index = INT_MIN;
1529 for (
Size i = 0; i < elution_profile.size(); ++i)
1531 if (elution_profile[i] > max)
1534 max = elution_profile[i];
1538 Int max_extension = (
Int)(elution_profile.size()) - 2 * max_index;
1540 double av_elution = 0;
1541 for (
Size i = 0; i < elution_profile.size(); ++i)
1543 av_elution += elution_profile[i];
1545 av_elution /= (
double)elution_profile.size();
1547 double sd_elution = 0;
1548 for (
Size i = 0; i < elution_profile.size(); ++i)
1550 sd_elution += (av_elution - elution_profile[i]) * (av_elution - elution_profile[i]);
1552 sd_elution /= (
double)(elution_profile.size() - 1);
1553 sd_elution = sqrt(sd_elution);
1557 for (iter = box.begin(); iter != box.end(); ++iter, ++index)
1559 av_mz += iter->second.mz;
1560 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1561 std::cout << iter->second.RT <<
"\t" << iter->second.mz <<
"\t" << iter->second.c + 1 << std::endl;
1564 av_mz /= (
double)box.size();
1568 if ((
Int)(box.begin()->second.RT_index) - 1 < 0)
1573 UInt pre_index = box.begin()->second.RT_index - 1;
1574 typename MSSpectrum::const_iterator c_iter = map[pre_index].MZBegin(av_mz);
1575 double pre_elution = 0;
1577 double mz_start = map[pre_index + 1][box.begin()->second.MZ_begin].getMZ();
1578 double mz_end = map[pre_index + 1][box.begin()->second.MZ_end].getMZ();
1580 typename MSSpectrum::const_iterator mz_start_iter = map[pre_index].MZBegin(mz_start), mz_end_iter = map[pre_index].MZBegin(mz_end);
1581 for (
typename MSSpectrum::const_iterator mz_iter = mz_start_iter; mz_iter != mz_end_iter; ++mz_iter)
1583 pre_elution += mz_iter->getIntensity();
1588 if (pre_elution <= av_elution - 2 * sd_elution)
1593 Int c_index = max_extension;
1594 Int first_index = box.begin()->second.RT_index;
1595 for (
Int i = 1; i < max_extension; ++i)
1597 c_index = first_index - i;
1604 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1605 std::cout << box.begin()->second.RT <<
"\t" << av_mz <<
"\t" << box.begin()->second.c + 1 <<
"\t" <<
" extending the box " << std::endl;
1608 push2Box_(av_mz, c_index, box.begin()->second.c, box.begin()->second.score, c_iter->getIntensity(),
1609 map[c_index].getRT(), box.
begin()->second.MZ_begin, box.begin()->second.MZ_end);
1613 template <
typename PeakType>
1617 typename std::multimap<double, Box>::iterator iter;
1618 typename Box::iterator box_iter;
1619 std::vector<BoxElement> final_box;
1620 double c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1621 double virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1623 typename std::pair<double, double> c_extend;
1626 Box& c_box = iter->second;
1627 av_score = 0; av_mz = 0; av_intens = 0; av_abs_intens = 0; count = 0;
1628 virtual_av_mz = 0; virtual_av_intens = 0; virtual_av_abs_intens = 0; virtual_count = 0;
1630 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1632 if (box_iter->second.score == 0)
1637 c_mz = box_iter->second.mz;
1638 virtual_av_intens += box_iter->second.intens;
1639 virtual_av_abs_intens += fabs(box_iter->second.intens);
1640 virtual_av_mz += c_mz * fabs(box_iter->second.intens);
1645 c_mz = box_iter->second.mz;
1646 av_score += box_iter->second.score;
1647 av_intens += box_iter->second.intens;
1648 av_abs_intens += fabs(box_iter->second.intens);
1649 av_mz += c_mz * fabs(box_iter->second.intens);
1656 av_intens = virtual_av_intens / virtual_count;
1658 av_mz = virtual_av_mz / virtual_av_abs_intens;
1664 av_mz /= av_abs_intens;
1668 c_box_element.
mz = av_mz;
1669 c_box_element.
c =
c;
1670 c_box_element.
score = av_score;
1671 c_box_element.
intens = av_intens;
1673 c_box_element.
RT = c_box.begin()->second.RT;
1675 final_box.push_back(c_box_element);
1678 UInt num_o_feature = final_box.size();
1679 if (num_o_feature == 0)
1686 std::vector<double> bwd_diffs(num_o_feature, 0);
1689 for (
UInt i = 1; i < num_o_feature; ++i)
1691 bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].
mz - final_box[i - 1].
mz);
1694 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1695 std::ofstream ofile_bwd(
"bwd_gpu.dat");
1696 for (
UInt i = 0; i < num_o_feature; ++i)
1698 ofile_bwd << final_box[i].mz <<
"\t" << bwd_diffs[i] << std::endl;
1703 for (
UInt i = 0; i < num_o_feature - 1; ++i)
1705 while (i < num_o_feature - 2)
1707 if (final_box[i].
score > 0 || final_box[i].
score == -1000)
1712 if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
1722 template <
typename PeakType>
1726 typename std::multimap<double, Box>::iterator iter;
1727 typename Box::iterator box_iter;
1728 UInt best_charge_index;
double best_charge_score, c_mz, c_RT;
UInt c_charge;
1729 double av_intens = 0, av_ref_intens = 0, av_score = 0, av_mz = 0, av_RT = 0, mz_cutoff, sum_of_ref_intenses_g;
1730 bool restart =
false;
1732 typename std::pair<double, double> c_extend;
1735 sum_of_ref_intenses_g = 0;
1736 Box& c_box = iter->second;
1742 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1744 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1745 if (trunc(box_iter->second.mz) == 874)
1746 std::cout << box_iter->second.c <<
"\t" << box_iter->second.intens <<
"\t" << box_iter->second.score << std::endl;
1749 if (box_iter->second.score == -1000)
1755 charge_votes[box_iter->second.c] += box_iter->second.intens;
1756 ++charge_binary_votes[box_iter->second.c];
1765 best_charge_index = 0; best_charge_score = 0;
1768 if (charge_votes[i] > best_charge_score)
1770 best_charge_index = i;
1771 best_charge_score = charge_votes[i];
1776 if (charge_binary_votes[best_charge_index] < RT_votes_cutoff && RT_votes_cutoff <= map.
size())
1781 c_charge = best_charge_index + 1;
1783 av_intens = 0; av_ref_intens = 0; av_score = 0; av_mz = 0; av_RT = 0;
1785 std::vector<DPosition<2> > point_set;
1786 double sum_of_ref_intenses_l;
1787 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1789 sum_of_ref_intenses_l = 0;
1790 c_mz = box_iter->second.mz;
1791 c_RT = box_iter->second.RT;
1797 point_set.push_back(
DPosition<2>(c_RT, c_mz + mz_cutoff / (
double)c_charge));
1799 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1800 std::cout <<
"Intenstype: " <<
intenstype_ << std::endl;
1805 const MSSpectrum& c_spec(map[box_iter->second.RT_index]);
1807 for (
unsigned int i = 0; i < mz_cutoff; ++i)
1813 while (h_iter != c_spec.begin())
1816 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1817 if (trunc(c_mz) == 874)
1819 std::cout <<
"cmz: " << c_mz + i *
Constants::IW_NEUTRON_MASS / c_charge <<
"\t" << hc_iter->getMZ() <<
"\t" << hc_iter->getIntensity() <<
"\t" << h_iter->getMZ() <<
"\t" << h_iter->getIntensity() << std::endl;
1824 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
1834 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1835 if (trunc(c_mz) == 874)
1840 sum_of_ref_intenses_l += hc_iter->getIntensity();
1841 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1842 if (trunc(c_mz) == 874)
1844 std::cout << sum_of_ref_intenses_l <<
"********" << std::endl;
1850 if (best_charge_index == box_iter->second.c)
1852 av_score += box_iter->second.score;
1853 av_intens += box_iter->second.intens;
1854 av_ref_intens += box_iter->second.ref_intens;
1855 sum_of_ref_intenses_g += sum_of_ref_intenses_l;
1856 av_mz += c_mz * box_iter->second.intens;
1862 av_ref_intens /= (
double)charge_binary_votes[best_charge_index];
1863 av_score /= (
double)charge_binary_votes[best_charge_index];
1864 av_RT /= (
double)c_box.size();
1866 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1867 if (trunc(av_mz) == 874)
1868 std::cout << av_mz <<
"\t" << best_charge_index <<
"\t" << best_charge_score << std::endl;
1875 c_feature.
setConvexHulls(std::vector<ConvexHull2D>(1, c_conv_hull));
1881 av_intens /= exp(-2 * lambda) * boost::math::cyl_bessel_i(0, 2 * lambda);
1885 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 1886 if (trunc(c_mz) == 874)
1888 std::cout << sum_of_ref_intenses_g <<
"####" << std::endl;
1892 av_intens = sum_of_ref_intenses_g;
1895 c_feature.
setMZ(av_mz);
1897 c_feature.
setRT(av_RT);
1899 feature_map.push_back(c_feature);
1905 template <
typename PeakType>
1907 const MSSpectrum& ref,
const double seed_mz,
const UInt c,
const UInt scan_index,
const bool check_PPMs,
const double transintens,
const double prev_score)
1909 typename MSSpectrum::const_iterator iter, ref_iter;
1913 iter = candidate.
MZBegin(seed_mz);
1915 if (iter == candidate.begin() || iter == candidate.end())
1920 std::pair<double, double> reals;
1921 ref_iter = ref.
MZBegin(seed_mz);
1923 double real_mz, real_intens;
1927 real_mz = reals.first; real_intens = reals.second;
1930 typename MSSpectrum::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
1931 while (h_iter != ref.begin())
1934 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
1949 real_mz = reals.first; real_intens = reals.second;
1950 if (real_mz <= 0 || real_intens <= 0)
1954 real_mz = h_iter->getMZ();
1955 real_intens = h_iter->getIntensity();
1960 reals = std::pair<double, double>(seed_mz, ref_iter->getIntensity());
1961 real_mz = reals.first; real_intens = reals.second;
1963 if (real_mz <= 0 || real_intens <= 0)
1965 typename MSSpectrum::const_iterator h_iter = ref_iter, hc_iter = ref_iter;
1966 while (h_iter != ref.begin())
1969 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
1983 real_mz = h_iter->getMZ(); real_intens = h_iter->getIntensity();
1984 if (real_mz <= 0 || real_intens <= 0)
1988 real_mz = h_iter->getMZ();
1989 real_intens = h_iter->getIntensity();
1993 double c_score =
scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2002 typename MSSpectrum::const_iterator real_r_MZ_iter = ref.
MZBegin(real_l_MZ_iter, real_mz + mz_cutoff / (c + 1.), ref.end());
2003 if (real_r_MZ_iter == ref.end())
2009 UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2010 UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2012 if (prev_score == -1000)
2014 push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2018 push2Box_(real_mz, scan_index, c, c_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2023 template <
typename PeakType>
2025 const MSSpectrum& ref,
const double seed_mz,
const UInt c,
const UInt scan_index,
const bool check_PPMs,
const double transintens,
const double prev_score)
2027 typename MSSpectrum::const_iterator iter, ref_iter;
2031 iter = candidate.
MZBegin(seed_mz);
2033 if (iter == candidate.
begin() || iter == candidate.
end())
2038 std::pair<double, double> reals;
2039 ref_iter = ref.
MZBegin(seed_mz);
2041 double real_mz, real_intens;
2045 real_mz = reals.first; real_intens = reals.second;
2048 auto h_iter = ref_iter, hc_iter = ref_iter;
2049 while (h_iter != ref.begin())
2052 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2068 real_mz = reals.first; real_intens = reals.second;
2070 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 2071 std::cout <<
"Plausibility check old_mz: " << iter->getMZ() <<
"\t" << real_mz << std::endl;
2074 if (real_mz <= 0 || real_intens <= 0)
2078 real_mz = h_iter->getMZ();
2079 real_intens = h_iter->getIntensity();
2084 reals = std::pair<double, double>(seed_mz, ref_iter->getIntensity());
2085 real_mz = reals.first; real_intens = reals.second;
2087 if (real_mz <= 0 || real_intens <= 0)
2089 auto h_iter = ref_iter, hc_iter = ref_iter;
2090 while (h_iter != ref.begin())
2093 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2107 real_mz = h_iter->getMZ(); real_intens = h_iter->getIntensity();
2108 if (real_mz <= 0 || real_intens <= 0)
2112 real_mz = h_iter->getMZ();
2113 real_intens = h_iter->getIntensity();
2117 double c_score =
scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2126 typename MSSpectrum::const_iterator real_r_MZ_iter = ref.
MZBegin(real_l_MZ_iter, real_mz + mz_cutoff / (c + 1.), ref.end());
2127 if (real_r_MZ_iter == ref.end())
2133 UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2134 UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2136 if (prev_score == -1000)
2138 push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2142 push2Box_(real_mz, scan_index, c, c_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2148 template <
typename PeakType>
2155 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 2156 std::cout << ::std::setprecision(8) << std::fixed << c_mz <<
"\t =(" <<
"ISO_WAVE" <<
")> " <<
"REJECT \t" << ppms <<
" (rule: " <<
peptideMassRule_(mass) <<
" got: " << mass <<
")" << std::endl;
2158 return std::pair<double, double>(-1, -1);
2161 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET 2162 std::cout << ::std::setprecision(8) << std::fixed << c_mz <<
"\t =(" <<
"ISO_WAVE" <<
")> " <<
"ACCEPT \t" << ppms <<
" (rule: " <<
peptideMassRule_(mass) <<
" got: " << mass <<
")" << std::endl;
2164 return std::pair<double, double>(c_mz, ref.
MZBegin(c_mz)->getIntensity());
2169 #pragma clang diagnostic pop double getSdIntens_(const TransSpectrum &scan, const double mean)
Computes the standard deviation (neglecting negative values) of the (transformed) intensities of scan...
Definition: IsotopeWaveletTransform.h:1213
bool hr_data_
Definition: IsotopeWaveletTransform.h:485
double getPPMs_(const double mass_a, const double mass_b) const
Returns the parts-per-million deviation of the masses.
Definition: IsotopeWaveletTransform.h:470
IntensityType getIntensity() const
Definition: Peak2D.h:166
Size size() const
Definition: IsotopeWaveletTransform.h:158
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:196
UInt MZ_end
Index.
Definition: IsotopeWaveletTransform.h:85
Iterator MZEnd(CoordinateType mz)
Binary search for peak range end (returns the past-the-end iterator)
const double PEPTIDE_MASS_RULE_FACTOR
Definition: IsotopeWaveletConstants.h:49
int Int
Signed integer type.
Definition: Types.h:102
std::vector< float > zeros_
Definition: IsotopeWaveletTransform.h:491
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:214
virtual ~IsotopeWaveletTransform()
Destructor.
Definition: IsotopeWaveletTransform.h:557
std::vector< double > c_spacings_
Definition: IsotopeWaveletTransform.h:480
const double IW_HALF_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:55
virtual void identifyCharge(const MSSpectrum &candidates, const MSSpectrum &ref, const UInt scan_index, const UInt c, const double ampl_cutoff, const bool check_PPMs)
Given an isotope wavelet transformed spectrum candidates, this function assigns to every significant ...
Definition: IsotopeWaveletTransform.h:696
void sortByIntensity(bool reverse=false)
Sorting.
Definition: ConstRefVector.h:741
void setCharge(const ChargeType &ch)
Set charge state.
std::multimap< double, Box > open_boxes_
Definition: IsotopeWaveletTransform.h:476
double getAvMZSpacing_(const MSSpectrum &scan)
Computes the average MZ spacing of scan.
Definition: IsotopeWaveletTransform.h:1180
void setTransIntensity(const UInt i, const double intens)
Definition: IsotopeWaveletTransform.h:152
IsotopeWaveletTransform()
Default Constructor.
Definition: IsotopeWaveletTransform.h:519
double getMinSpacing() const
Definition: IsotopeWaveletTransform.h:312
virtual double scoreThis_(const TransSpectrum &candidate, const UInt peak_cutoff, const double seed_mz, const UInt c, const double ampl_cutoff)
Given a candidate for an isotopic pattern, this function computes the corresponding score...
Definition: IsotopeWaveletTransform.h:967
UInt MZ_begin
Index.
Definition: IsotopeWaveletTransform.h:84
Size getMaxScanSize() const
Definition: IsotopeWaveletTransform.h:317
double ref_intens
Definition: IsotopeWaveletTransform.h:81
double getTransIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:146
std::vector< double > prod_
Definition: IsotopeWaveletTransform.h:480
Int from_max_to_right_
Definition: IsotopeWaveletTransform.h:487
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
double getAvIntens_(const TransSpectrum &scan)
Computes the average (transformed) intensity (neglecting negative values) of scan.
Definition: IsotopeWaveletTransform.h:1199
Size size() const
Definition: MSExperiment.h:127
std::vector< float > scores_
Definition: IsotopeWaveletTransform.h:491
FeatureMap mapSeeds2Features(const PeakMap &map, const UInt RT_votes_cutoff)
Filters the candidates further more and maps the internally used data structures to the OpenMS framew...
Definition: IsotopeWaveletTransform.h:1723
Int from_max_to_left_
Definition: IsotopeWaveletTransform.h:487
MSSpectrum::const_iterator MZEnd(const double mz) const
Definition: IsotopeWaveletTransform.h:184
An LC-MS feature.
Definition: Feature.h:70
UInt data_length_
Definition: IsotopeWaveletTransform.h:484
virtual void initializeScan(const MSSpectrum &c_ref, const UInt c=0)
Definition: IsotopeWaveletTransform.h:644
double mz
The monoisotopic position.
Definition: IsotopeWaveletTransform.h:77
MSSpectrum::const_iterator MZBegin(const double mz) const
Definition: IsotopeWaveletTransform.h:177
std::vector< double > c_mzs_
Definition: IsotopeWaveletTransform.h:480
const double PEPTIDE_MASS_RULE_BOUND
Definition: IsotopeWaveletConstants.h:50
MSSpectrum::const_iterator begin() const
Definition: IsotopeWaveletTransform.h:198
std::multimap< double, Box > front_boxes_
Definition: IsotopeWaveletTransform.h:476
double max_mz_cutoff_
Definition: IsotopeWaveletTransform.h:490
static double getValueByLambda(const double lambda, const double tz1)
Returns the value of the isotope wavelet at position t via a fast table lookup. Usually, you do not need to call this function. Please use.
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:77
Size max_scan_size_
Definition: IsotopeWaveletTransform.h:483
Iterator begin()
Definition: MSExperiment.h:157
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
double getRT() const
Definition: IsotopeWaveletTransform.h:128
static UInt getNumPeakCutOff(const double mass, const UInt z)
Mutable iterator for the ConstRefVector.
Definition: ConstRefVector.h:236
std::multimap< double, Box > end_boxes_
Definition: IsotopeWaveletTransform.h:476
void clusterSeeds_(const TransSpectrum &candidates, const MSSpectrum &ref, const UInt scan_index, const UInt c, const bool check_PPMs)
Clusters the seeds stored by push2TmpBox_.
Definition: IsotopeWaveletTransform.h:1614
double getMZ(const UInt i) const
Definition: IsotopeWaveletTransform.h:134
double getSigma() const
Definition: IsotopeWaveletTransform.h:300
const double IW_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:54
double score
The associated score.
Definition: IsotopeWaveletTransform.h:79
const double IW_PROTON_MASS
Definition: IsotopeWaveletConstants.h:68
TransSpectrum()
Definition: IsotopeWaveletTransform.h:101
virtual void push2Box_(const double mz, const UInt scan, UInt c, const double score, const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end, const double ref_intens)
Inserts a potential isotopic pattern into an open box or - if no such box exists - creates a new one...
Definition: IsotopeWaveletTransform.h:1228
bool intensityPointerComparator(PeakType *a, PeakType *b)
Definition: IsotopeWaveletTransform.h:507
double min_spacing_
Definition: IsotopeWaveletTransform.h:490
const MSSpectrum * getRefSpectrum()
Definition: IsotopeWaveletTransform.h:164
std::vector< double > interpol_ys_
Definition: IsotopeWaveletTransform.h:481
virtual ~TransSpectrum()
Definition: IsotopeWaveletTransform.h:114
A class implementing the isotope wavelet transform. If you just want to find features using the isoto...
Definition: IsotopeWaveletTransform.h:69
UInt max_charge_
Definition: IsotopeWaveletTransform.h:484
double av_MZ_spacing_
Definition: IsotopeWaveletTransform.h:479
void setSigma(const double sigma)
Definition: IsotopeWaveletTransform.h:305
String intenstype_
Definition: IsotopeWaveletTransform.h:486
double getLinearInterpolation(const double mz_a, const double intens_a, const double mz_pos, const double mz_b, const double intens_b)
Computes a linear (intensity) interpolation.
Definition: IsotopeWaveletTransform.h:295
std::vector< int > indices_
Definition: IsotopeWaveletTransform.h:488
void setOverallQuality(QualityType q)
Set the overall quality.
double intens
The transformed intensity at the monoisotopic mass.
Definition: IsotopeWaveletTransform.h:80
A container for features.
Definition: FeatureMap.h:95
Internally (only by GPUs) used data structure . It allows efficient data exchange between CPU and GPU...
Definition: IsotopeWaveletTransform.h:94
UInt RT_index
The elution time (map) index.
Definition: IsotopeWaveletTransform.h:83
virtual std::multimap< double, Box > getClosedBoxes()
Returns the closed boxes.
Definition: IsotopeWaveletTransform.h:276
void addPoints(const PointArrayType &points)
void sampleTheCMarrWavelet_(const MSSpectrum &scan, const Int wavelet_length, const Int mz_index, const UInt charge)
Internally used data structure.
Definition: IsotopeWaveletTransform.h:75
unsigned int UInt
Unsigned integer type.
Definition: Types.h:94
const double PEPTIDE_MASS_RULE_THEO_PPM_BOUND
Definition: IsotopeWaveletConstants.h:51
double sigma_
Definition: IsotopeWaveletTransform.h:479
virtual void getTransformHighRes(MSSpectrum &c_trans, const MSSpectrum &c_ref, const UInt c)
Computes the isotope wavelet transform of charge state c.
Definition: IsotopeWaveletTransform.h:606
std::vector< double > psi_
Definition: IsotopeWaveletTransform.h:480
virtual void computeMinSpacing(const MSSpectrum &c_ref)
Definition: IsotopeWaveletTransform.h:686
std::vector< double > interpol_xs_
Definition: IsotopeWaveletTransform.h:481
UInt max_num_peaks_per_pattern_
Definition: IsotopeWaveletTransform.h:484
virtual void getTransform(MSSpectrum &c_trans, const MSSpectrum &c_ref, const UInt c)
Computes the isotope wavelet transform of charge state c.
Definition: IsotopeWaveletTransform.h:563
const MSSpectrum * getRefSpectrum() const
Definition: IsotopeWaveletTransform.h:170
This vector holds pointer to the elements of another container.
Definition: ConstRefVector.h:70
std::vector< float > * trans_intens_
The intensities of the transform.
Definition: IsotopeWaveletTransform.h:206
The representation of a 1D spectrum.
Definition: MSSpectrum.h:67
const MSSpectrum * reference_
The reference spectrum.
Definition: IsotopeWaveletTransform.h:205
UInt c
Note, this is not the charge (it is charge-1!!!)
Definition: IsotopeWaveletTransform.h:78
virtual bool checkPositionForPlausibility_(const TransSpectrum &candidate, const MSSpectrum &ref, const double seed_mz, const UInt c, const UInt scan_index, const bool check_PPMs, const double transintens, const double prev_score)
A ugly but necessary function to handle "off-by-1-Dalton predictions" due to idiosyncrasies of the da...
Definition: IsotopeWaveletTransform.h:2024
void extendBox_(const PeakMap &map, const Box &box)
A currently still necessary function that extends the box box in order to capture also signals whose ...
Definition: IsotopeWaveletTransform.h:1508
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:172
const double IW_QUARTER_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:56
TransSpectrum(const MSSpectrum *reference)
Definition: IsotopeWaveletTransform.h:107
bool positionComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:513
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:72
void updateBoxStates(const PeakMap &map, const Size scan_index, const UInt RT_interleave, const UInt RT_votes_cutoff, const Int front_bound=-1, const Int end_bound=-1)
A function keeping track of currently open and closed sweep line boxes. This function is used by the ...
Definition: IsotopeWaveletTransform.h:1439
std::vector< double > xs_
Definition: IsotopeWaveletTransform.h:480
std::multimap< double, Box > closed_boxes_
Definition: IsotopeWaveletTransform.h:476
virtual void push2TmpBox_(const double mz, const UInt scan, UInt charge, const double score, const double intens, const double rt, const UInt MZ_begin, const UInt MZ_end)
Essentially the same function as.
Definition: IsotopeWaveletTransform.h:1334
void setConvexHulls(const std::vector< ConvexHull2D > &hulls)
Set the convex hulls of single mass traces.
virtual std::pair< double, double > checkPPMTheoModel_(const MSSpectrum &ref, const double c_mz, const UInt c)
Definition: IsotopeWaveletTransform.h:2149
static double getLambdaL(const double m)
Returns the mass-parameter lambda (linear fit).
A more convenient string class.
Definition: String.h:58
virtual void destroy()
Definition: IsotopeWaveletTransform.h:119
double peptideMassRule_(const double c_mass) const
Returns the monoisotopic mass (with corresponding decimal values) we would expect at c_mass...
Definition: IsotopeWaveletTransform.h:447
double RT
The elution time (not the scan index)
Definition: IsotopeWaveletTransform.h:82
std::multimap< UInt, BoxElement > Box
Key: RT index, value: BoxElement.
Definition: IsotopeWaveletTransform.h:88
static UInt getMzPeakCutOffAtMonoPos(const double mass, const UInt z)
static double mean(IteratorType begin, IteratorType end)
Calculates the mean of a range of values.
Definition: StatisticFunctions.h:133
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:54
MSSpectrum::const_iterator end() const
Definition: IsotopeWaveletTransform.h:191
Iterator MZBegin(CoordinateType mz)
Binary search for peak range begin.
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:202
bool intensityComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:495
std::vector< std::multimap< double, Box > > * tmp_boxes_
Definition: IsotopeWaveletTransform.h:477
bool intensityAscendingComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:501
static IsotopeWavelet * init(const double max_m, const UInt max_charge)
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
const unsigned int DEFAULT_NUM_OF_INTERPOLATION_POINTS
Definition: IsotopeWaveletConstants.h:43
double getLinearInterpolation(const typename MSSpectrum::const_iterator &left_iter, const double mz_pos, const typename MSSpectrum::const_iterator &right_iter)
Computes a linear (intensity) interpolation.
Definition: IsotopeWaveletTransform.h:284
double getRefIntensity(const UInt i) const
Definition: IsotopeWaveletTransform.h:140