35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_ISOTOPEWAVELETTRANSFORM_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_ISOTOPEWAVELETTRANSFORM_H
48 #include <boost/math/special_functions/bessel.hpp>
55 #ifdef OPENMS_HAS_CUDA
61 #ifndef CUDA_SAFE_CALL
62 #define CUDA_SAFE_CALL(call) call;
86 template <
typename PeakType>
106 typedef std::multimap<UInt, BoxElement>
Box;
172 (*trans_intens_)[i] = intens;
273 const DoubleReal ampl_cutoff,
const bool check_PPMs);
277 #ifdef OPENMS_HAS_CUDA
283 virtual void finalizeScanCuda();
306 const DoubleReal ampl_cutoff,
const bool check_PPMs);
321 const UInt RT_votes_cutoff,
const Int front_bound = -1,
const Int end_bound = -1);
344 return left_iter->getIntensity() + (right_iter->getIntensity() - left_iter->getIntensity()) / (right_iter->getMZ() - left_iter->getMZ()) * (mz_pos - left_iter->getMZ());
355 return intens_a + (intens_b - intens_a) / (mz_b - mz_a) * (mz_pos - mz_a);
487 const UInt scan_index,
const UInt c,
const bool check_PPMs);
494 const UInt scan_index,
const UInt c,
const bool check_PPMs);
512 if (new_frac_mass - old_frac_mass > 0.5)
517 if (new_frac_mass - old_frac_mass < -0.5)
530 return fabs(mass_a - mass_b) / (0.5 * (mass_a + mass_b)) * 1e6;
552 #ifdef OPENMS_HAS_CUDA
555 UInt largest_array_size_, overall_size_, block_size_, to_load_, to_compute_;
557 void* cuda_device_intens_;
558 void* cuda_device_pos_;
559 void* cuda_device_trans_intens_;
560 void* cuda_device_fwd2_;
561 void* cuda_device_posindices_sorted_;
562 void* cuda_device_trans_intens_sorted_;
563 void* cuda_device_scores_;
564 std::vector<float> cuda_positions_, cuda_intensities_;
565 dim3 dimGrid_, dimBlock_;
572 template <
typename PeakType>
578 template <
typename PeakType>
584 template <
typename PeakType>
590 template <
typename PeakType>
596 template <
typename PeakType>
599 tmp_boxes_ =
new std::vector<std::multimap<DoubleReal, Box> >(1);
603 max_num_peaks_per_pattern_ = 3;
606 #ifdef OPENMS_HAS_CUDA
607 largest_array_size_ = 0;
614 template <
typename PeakType>
617 max_charge_ = max_charge;
618 max_scan_size_ = max_scan_size;
620 intenstype_ = intenstype;
621 tmp_boxes_ =
new std::vector<std::multimap<DoubleReal, Box> >(max_charge);
622 if (max_scan_size <= 0)
631 #ifdef OPENMS_HAS_CUDA
636 max_scan_size_ = max_scan_size * (4 * (max_mz_cutoff_ - 1) - 1);
640 cuda_positions_.reserve(largest_array_size_);
641 cuda_intensities_.reserve(largest_array_size_);
642 indices_.resize(largest_array_size_);
643 for (
UInt q = 0; q < largest_array_size_; ++q)
648 h_data_ = (
float*) malloc(largest_array_size_ *
sizeof(
float));
649 h_pos_ = (
int*) malloc(largest_array_size_ *
sizeof(
int));
655 Int size_estimate((
Int)ceil(max_scan_size_ / (max_mz - min_mz)));
657 psi_.reserve(to_reserve);
658 prod_.reserve(to_reserve);
659 xs_.reserve(to_reserve);
666 Int size_estimate((
Int)ceil(max_scan_size_ / (max_mz - min_mz)));
668 psi_.reserve(to_reserve);
669 prod_.reserve(to_reserve);
670 xs_.reserve(to_reserve);
677 template <
typename PeakType>
680 #ifdef OPENMS_HAS_CUDA
692 template <
typename PeakType>
695 Int spec_size((
Int)c_ref.size());
699 DoubleReal value, T_boundary_left, T_boundary_right, old, c_diff, current, old_pos, my_local_MZ, my_local_lambda, origin, c_mz;
701 for (
Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
704 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_;
709 for (
Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
711 if (current_conv_pos >= spec_size)
713 value += 0.5 * old * min_spacing_;
717 c_mz = c_ref[current_conv_pos].getMZ();
718 c_diff = c_mz + origin;
721 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;
723 value += 0.5 * (current + old) * (c_mz - old_pos);
731 c_trans[my_local_pos].setIntensity(value);
735 template <
typename PeakType>
738 Int spec_size((
Int)c_ref.size());
742 DoubleReal value, T_boundary_left, T_boundary_right, c_diff, current, my_local_MZ, my_local_lambda, origin, c_mz;
744 for (
Int my_local_pos = 0; my_local_pos < spec_size; ++my_local_pos)
753 for (
Int current_conv_pos = std::max(0, my_local_pos - from_max_to_left_); c_diff < T_boundary_right; ++current_conv_pos)
755 if (current_conv_pos >= spec_size)
760 c_mz = c_ref[current_conv_pos].getMZ();
761 c_diff = c_mz + origin;
764 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;
769 c_trans[my_local_pos].setIntensity(value);
773 template <
typename PeakType>
776 data_length_ = (
UInt) c_ref.size();
777 computeMinSpacing(c_ref);
778 Int wavelet_length = 0, quarter_length = 0;
784 for (
UInt i = 0; i < data_length_; ++i)
787 start_iter = c_ref.
MZEnd(c_ref[i].getMZ());
788 end_iter = c_ref.
MZBegin(c_ref[i].getMZ() + c_mz_cutoff);
789 wavelet_length = std::max((
SignedSize) wavelet_length, distance(start_iter, end_iter) + 1);
791 quarter_length = std::max((
SignedSize) quarter_length, distance(end_iter, start_iter) + 1);
798 wavelet_length = (
UInt) ceil(max_mz_cutoff_ / min_spacing_);
803 if (wavelet_length > (
Int) c_ref.size())
805 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;
806 std::cout <<
"Minimal spacing: " << min_spacing_ << std::endl;
807 std::cout <<
"Warning/Error generated at scan with RT " << c_ref.
getRT() <<
"." << std::endl;
811 from_max_to_left_ = max_index;
812 from_max_to_right_ = wavelet_length - 1 - from_max_to_left_;
815 template <
typename PeakType>
818 min_spacing_ = INT_MAX;
819 for (
UInt c_conv_pos = 1; c_conv_pos < c_ref.size(); ++c_conv_pos)
821 min_spacing_ = std::min(min_spacing_, c_ref[c_conv_pos].getMZ() - c_ref[c_conv_pos - 1].getMZ());
825 #ifdef OPENMS_HAS_CUDA
826 template <
typename PeakType>
829 (cudaFree(cuda_device_pos_));
830 (cudaFree(cuda_device_intens_));
831 (cudaFree(cuda_device_trans_intens_));
832 (cudaFree(cuda_device_fwd2_));
833 (cudaFree(cuda_device_trans_intens_sorted_));
834 (cudaFree(cuda_device_posindices_sorted_));
835 (cudaFree(cuda_device_scores_));
838 template <
typename PeakType>
839 int IsotopeWaveletTransform<PeakType>::initializeScanCuda(
const MSSpectrum<PeakType>& scan,
const UInt c)
841 data_length_ = scan.size();
843 std::vector<float> pre_positions(data_length_), pre_intensities(data_length_);
845 min_spacing_ = INT_MAX;
846 pre_positions[0] = scan[0].getMZ();
847 pre_intensities[0] = scan[0].getIntensity();
849 for (
UInt i = 1; i < data_length_; ++i)
851 pre_positions[i] = scan[i].getMZ();
852 c_spacing = pre_positions[i] - pre_positions[i - 1];
853 if (c_spacing < min_spacing_)
855 min_spacing_ = c_spacing;
857 pre_intensities[i] = scan[i].getIntensity();
859 if (min_spacing_ == INT_MAX)
861 std::cout <<
"Scan consits of a single point. Unable to compute transform." << std::endl;
866 UInt wavelet_length = 0, quarter_length = 0;
870 typename MSSpectrum<PeakType>::const_iterator start_iter, end_iter;
871 for (
UInt i = 0; i < data_length_; ++i)
874 start_iter = scan.MZEnd(scan[i].getMZ());
875 end_iter = scan.MZBegin(scan[i].getMZ() + c_mz_cutoff);
876 wavelet_length = std::max((
long int) wavelet_length, distance(start_iter, end_iter) + 1);
878 quarter_length = std::max((
long int) quarter_length, distance(end_iter, start_iter) + 1);
885 wavelet_length = (
UInt) ceil(max_mz_cutoff_ / min_spacing_);
889 if (wavelet_length > data_length_ || wavelet_length == 1)
891 std::cout <<
"Warning: the extremal length of the wavelet is larger (" << wavelet_length <<
") than the number of data points (" << data_length_ <<
"). This might (!) severely affect the transform." << std::endl;
892 std::cout <<
"Minimal spacing: " << min_spacing_ << std::endl;
893 std::cout <<
"Warning/Error generated at scan with RT " << scan.getRT() <<
"." << std::endl;
899 max_index = quarter_length;
905 from_max_to_left_ = max_index;
906 from_max_to_right_ = wavelet_length - 1 - from_max_to_left_;
909 to_load_ = problem_size + from_max_to_left_ + from_max_to_right_;
911 UInt missing_points = problem_size - (data_length_ % problem_size);
912 overall_size_ = wavelet_length - 1 + data_length_ + missing_points;
914 num_elements_ = overall_size_;
915 Int dev_num_elements = 1, tmp = overall_size_ >> 1;
920 dev_num_elements <<= 1;
924 if (num_elements_ > dev_num_elements)
926 dev_num_elements <<= 1;
934 overall_size_ = dev_num_elements;
936 cuda_intensities_.resize(overall_size_, 0); cuda_positions_.resize(overall_size_, 0);
938 float first_pos = pre_positions[0];
939 for (
Int i = 0; i < from_max_to_left_; ++i)
941 cuda_positions_[i] = first_pos - (from_max_to_left_ - i) * min_spacing_;
944 for (
UInt i = 0; i < data_length_; ++i)
946 cuda_positions_[from_max_to_left_ + i] = pre_positions[i];
947 cuda_intensities_[from_max_to_left_ + i] = pre_intensities[i];
950 float last_pos = pre_positions[pre_positions.size() - 1];
951 for (
UInt i = 0; i < missing_points + from_max_to_right_ + dev_num_elements - num_elements_; ++i)
953 cuda_positions_[from_max_to_left_ + data_length_ + i] = last_pos + (i + 1) * min_spacing_;
956 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
957 std::stringstream name; name <<
"cuda_input_" << scan.getRT() <<
".out\0";
958 std::fstream outfile(name.str().c_str(), std::ios::out);
959 for (
size_t i = 0; i < overall_size_; ++i)
960 outfile << cuda_positions_[i] <<
" " << cuda_intensities_[i] << std::endl;
966 to_compute_ = problem_size;
967 dimGrid_ = dim3((data_length_ + missing_points) / problem_size);
969 (cudaMalloc(&cuda_device_posindices_sorted_, overall_size_ *
sizeof(
int)));
970 (cudaMalloc(&cuda_device_pos_, overall_size_ *
sizeof(
float)));
971 (cudaMemcpy(cuda_device_pos_, &(cuda_positions_[0]), overall_size_ *
sizeof(
float), cudaMemcpyHostToDevice));
972 (cudaMalloc(&cuda_device_intens_, overall_size_ *
sizeof(
float)));
973 (cudaMemcpy(cuda_device_intens_, &(cuda_intensities_[0]), overall_size_ *
sizeof(
float), cudaMemcpyHostToDevice));
974 (cudaMalloc(&cuda_device_trans_intens_, overall_size_ *
sizeof(
float)));
975 (cudaMalloc(&cuda_device_fwd2_, overall_size_ *
sizeof(
float)));
976 (cudaMalloc(&cuda_device_trans_intens_sorted_, overall_size_ *
sizeof(
float)));
978 c_sorted_candidate_.resize(overall_size_);
979 scores_.resize(data_length_);
980 zeros_.resize(overall_size_);
981 memset(&zeros_[0], 0., overall_size_ *
sizeof(
float));
983 (cudaMalloc(&cuda_device_scores_, overall_size_ *
sizeof(
float)));
988 template <
typename PeakType>
989 void IsotopeWaveletTransform<PeakType>::getTransformCuda(TransSpectrum& c_trans,
const UInt c)
991 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
992 std::cout <<
"res in vector" << std::endl;
993 std::vector<float> res(overall_size_, 0);
995 (cudaMemcpy(cuda_device_trans_intens_, &zeros_[0], overall_size_ *
sizeof(
float), cudaMemcpyHostToDevice));
997 (cudaMemcpy(cuda_device_fwd2_, &zeros_[0], overall_size_ *
sizeof(
float), cudaMemcpyHostToDevice));
998 getExternalCudaTransforms(dimGrid_, dimBlock_, (
float*)cuda_device_pos_, (
float*)cuda_device_intens_, from_max_to_left_, from_max_to_right_, (
float*)cuda_device_trans_intens_,
999 c + 1, to_load_, to_compute_, data_length_, (
float*)cuda_device_fwd2_, hr_data_);
1001 (cudaMemcpy(cuda_device_trans_intens_sorted_, cuda_device_fwd2_, overall_size_ *
sizeof(
float), cudaMemcpyDeviceToDevice));
1003 (cudaMemcpy(&((*c_trans.trans_intens_)[0]), (
float*)cuda_device_trans_intens_ + from_max_to_left_, data_length_ *
sizeof(
float), cudaMemcpyDeviceToHost));
1004 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1005 (cudaMemcpy(&(res[0]), (
float*)cuda_device_trans_intens_ + from_max_to_left_, data_length_ *
sizeof(
float), cudaMemcpyDeviceToHost));
1006 for (
UInt i = 0; i < data_length_; ++i)
1008 c_trans.setTransIntensity(i, res[i]);
1013 template <
typename PeakType>
1014 int IsotopeWaveletTransform<PeakType>::sortCuda(MSSpectrum<PeakType>& sorted)
1016 (cudaMemcpy(cuda_device_posindices_sorted_, &indices_[0], overall_size_ *
sizeof(
int), cudaMemcpyHostToDevice));
1017 Int gpu_index =
sortOnDevice((
float*)cuda_device_trans_intens_sorted_, (
int*) cuda_device_posindices_sorted_, overall_size_, 0);
1024 (cudaMemcpy(h_data_, (
float*)cuda_device_trans_intens_sorted_ + gpu_index,
sizeof(
float) * (overall_size_ - gpu_index), cudaMemcpyDeviceToHost));
1025 (cudaMemcpy(h_pos_, (
int*)cuda_device_posindices_sorted_ + gpu_index,
sizeof(
int) * (overall_size_ - gpu_index), cudaMemcpyDeviceToHost));
1027 for (
UInt i = 0; i < (overall_size_ - gpu_index); ++i)
1029 sorted[i].setIntensity(h_data_[i]);
1030 sorted[i].setMZ(cuda_positions_[h_pos_[i]]);
1036 template <
typename PeakType>
1037 void IsotopeWaveletTransform<PeakType>::identifyChargeCuda(
const TransSpectrum& candidates,
1038 const UInt scan_index,
const UInt c,
const DoubleReal ampl_cutoff,
const bool check_PPMs)
1040 const MSSpectrum<PeakType>& ref(*candidates.getRefSpectrum());
1041 UInt index, MZ_start, MZ_end;
1042 typename MSSpectrum<PeakType>::iterator iter, bound_iter;
1043 typename MSSpectrum<PeakType>::const_iterator iter_start, iter_end, iter_p, iter2, seed_iter;
1044 DoubleReal mz_cutoff, seed_mz, c_av_intens = 0, c_score = 0, c_sd_intens = 0, threshold = 0, help_mz;
1046 Int gpu_index = sortCuda(c_sorted_candidate_), c_index;
1052 std::vector<UInt> processed(data_length_, 0);
1053 if (ampl_cutoff < 0)
1059 c_av_intens = getAvIntens_(candidates);
1060 c_sd_intens = getSdIntens_(candidates, c_av_intens);
1061 threshold = ampl_cutoff * c_sd_intens + c_av_intens;
1064 Int num_of_scores = overall_size_ - gpu_index;
1066 (cudaMemcpy(cuda_device_scores_, &zeros_[0], num_of_scores *
sizeof(
float), cudaMemcpyHostToDevice));
1068 scoreOnDevice((
int*)cuda_device_posindices_sorted_, (
float*)cuda_device_trans_intens_, (
float*)cuda_device_pos_, (
float*)cuda_device_scores_,
1069 c, num_of_scores, overall_size_, max_num_peaks_per_pattern_, threshold);
1071 (cudaMemcpy(&scores_[0], cuda_device_scores_, num_of_scores *
sizeof(
float), cudaMemcpyDeviceToHost));
1073 std::vector<float>::iterator score_iter;
1074 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1075 std::stringstream stream;
1076 stream <<
"sorted_gpu_" << candidates.getRT() <<
"_" << c + 1 <<
".trans\0";
1077 std::ofstream ofile(stream.str().c_str());
1078 for (c_index = overall_size_ - gpu_index - 1, score_iter = scores_.begin() + num_of_scores - 1; c_index >= 0; --c_index, --score_iter)
1080 ofile << c_sorted_candidate_[c_index].getMZ() <<
"\t" << c_sorted_candidate_[c_index].getIntensity() <<
"\t" << *score_iter << std::endl;
1085 for (c_index = overall_size_ - gpu_index - 1, score_iter = scores_.begin() + num_of_scores - 1; c_index >= 0; --c_index, --score_iter)
1087 seed_mz = c_sorted_candidate_[c_index].getMZ();
1093 index = h_pos_[c_index] - from_max_to_left_;
1094 seed_iter = ref.begin() + index;
1096 if (seed_iter == ref.end() || processed[distance(ref.begin(), seed_iter)] || index <= 0)
1107 iter_end = ref.MZEnd(seed_iter, seed_mz + mz_cutoff / (c + 1.), ref.end());
1109 if (iter_end == ref.end())
1114 MZ_start = distance(ref.begin(), iter_start);
1115 MZ_end = distance(ref.begin(), iter_end);
1117 memset(&(processed[MZ_start]), 1,
sizeof(
UInt) * (MZ_end - MZ_start + 1));
1119 c_score = *score_iter;
1121 if (c_score <= 0 && c_score != -1000)
1129 push2TmpBox_(seed_mz, scan_index, c, c_score, c_sorted_candidate_[c_index].getIntensity(), ref.getRT(), MZ_start, MZ_end);
1133 iter2 = candidates.MZBegin(help_mz);
1135 if (iter2 == candidates.end() || iter2 == candidates.begin())
1143 if (iter2 != candidates.end())
1145 UInt dist = distance(candidates.begin(), iter2);
1146 push2TmpBox_(iter2->getMZ(), scan_index,
c, 0,
1147 getLinearInterpolation((iter2 - 1)->getMZ(), candidates.getTransIntensity(dist - 1), help_mz, iter2->getMZ(), candidates.getTransIntensity(dist)),
1148 candidates.getRT(), MZ_start, MZ_end);
1153 iter2 = candidates.MZBegin(help_mz);
1155 if (iter2 == candidates.end() || iter2 == candidates.begin())
1163 if (iter2 != candidates.end())
1165 UInt dist = distance(candidates.begin(), iter2);
1166 push2TmpBox_(iter2->getMZ(), scan_index,
c, 0,
1167 getLinearInterpolation((iter2 - 1)->getMZ(), candidates.getTransIntensity(dist - 1), help_mz, iter2->getMZ(), candidates.getTransIntensity(dist)),
1168 candidates.getRT(), MZ_start, MZ_end);
1173 clusterSeeds_(candidates, ref, scan_index, c, check_PPMs);
1179 template <
typename PeakType>
1183 Size scan_size(candidates.size());
1186 DoubleReal mz_cutoff, seed_mz, c_av_intens = 0, c_score = 0, c_sd_intens = 0, threshold = 0, help_mz, share, share_pos, bwd, fwd;
1187 UInt MZ_start, MZ_end;
1190 diffed[0].setIntensity(0); diffed[scan_size - 1].setIntensity(0);
1192 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1193 std::stringstream stream;
1194 stream <<
"diffed_" << ref.
getRT() <<
"_" << c + 1 <<
".trans\0";
1195 std::ofstream ofile(stream.str().c_str());
1200 for (
UInt i = 0; i < scan_size - 2; ++i)
1202 share = candidates[i + 1].getIntensity(), share_pos = candidates[i + 1].getMZ();
1203 bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
1204 fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
1206 if (!(bwd >= 0 && fwd <= 0) || share > ref[i + 1].getIntensity())
1208 diffed[i + 1].setIntensity(0);
1211 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1212 ofile << diffed[i + 1].getMZ() <<
"\t" << diffed[i + 1].getIntensity() << std::endl;
1218 for (
UInt i = 0; i < scan_size - 2; ++i)
1220 share = candidates[i + 1].getIntensity(), share_pos = candidates[i + 1].getMZ();
1221 bwd = (share - candidates[i].getIntensity()) / (share_pos - candidates[i].getMZ());
1222 fwd = (candidates[i + 2].getIntensity() - share) / (candidates[i + 2].getMZ() - share_pos);
1224 if (!(bwd >= 0 && fwd <= 0))
1226 diffed[i + 1].setIntensity(0);
1229 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1230 ofile << diffed[i + 1].getMZ() <<
"\t" << diffed[i + 1].getIntensity() << std::endl;
1234 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1243 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1244 std::stringstream stream2;
1245 stream2 <<
"sorted_cpu_" << candidates.
getRT() <<
"_" << c + 1 <<
".trans\0";
1246 std::ofstream ofile2(stream2.str().c_str());
1247 for (iter = c_sorted_candidate_.end() - 1; iter != c_sorted_candidate_.begin(); --iter)
1249 ofile2 << iter->getMZ() <<
"\t" << iter->getIntensity() << std::endl;
1254 std::vector<UInt> processed(scan_size, 0);
1256 if (ampl_cutoff < 0)
1262 c_av_intens = getAvIntens_(candidates);
1263 c_sd_intens = getSdIntens_(candidates, c_av_intens);
1264 threshold = ampl_cutoff * c_sd_intens + c_av_intens;
1267 for (iter = c_sorted_candidate_.end() - 1; iter != c_sorted_candidate_.begin(); --iter)
1269 if (iter->getIntensity() <= 0)
1274 seed_mz = iter->getMZ();
1275 seed_iter = ref.
MZBegin(seed_mz);
1277 if (seed_iter == ref.end() || processed[distance(ref.begin(), seed_iter)])
1288 iter_end = ref.
MZEnd(seed_iter, seed_mz + mz_cutoff / (c + 1.), ref.end());
1289 if (iter_end == ref.end())
1294 MZ_start = distance(ref.begin(), iter_start);
1295 MZ_end = distance(ref.begin(), iter_end);
1297 memset(&(processed[MZ_start]), 1,
sizeof(
UInt) * (MZ_end - MZ_start + 1));
1301 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1302 if (trunc(seed_mz) == 874)
1303 std::cout << seed_mz <<
"\t" << c_score << std::endl;
1306 if (c_score <= 0 && c_score != -1000)
1313 push2TmpBox_(seed_mz, scan_index, c, c_score, iter->getIntensity(), ref.
getRT(), MZ_start, MZ_end);
1316 iter2 = candidates.
MZBegin(help_mz);
1317 if (iter2 == candidates.end() || iter2 == candidates.begin())
1325 if (iter2 != candidates.end())
1327 push2TmpBox_(iter2->getMZ(), scan_index,
c, 0, getLinearInterpolation(iter2 - 1, help_mz, iter2), ref.
getRT(), MZ_start, MZ_end);
1332 iter2 = candidates.
MZBegin(help_mz);
1333 if (iter2 == candidates.end() || iter2 == candidates.begin())
1341 if (iter2 != candidates.end())
1343 push2TmpBox_(iter2->getMZ(), scan_index,
c, 0, getLinearInterpolation(iter2 - 1, help_mz, iter2), ref.
getRT(), MZ_start, MZ_end);
1348 clusterSeeds_(candidates, ref, scan_index, c, check_PPMs);
1351 #if defined(OPENMS_HAS_TBB) && defined(OPENMS_HAS_CUDA)
1352 template <
typename PeakType>
1355 typename std::multimap<DoubleReal, Box>::iterator front_iter, end_iter, best_match, help_iter;
1360 closed_boxes_.insert(*end_iter);
1363 typename std::multimap<DoubleReal, Box>& end_container(this->end_boxes_);
1364 typename std::multimap<DoubleReal, Box>& front_container(later_iwt->
front_boxes_);
1366 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1367 std::cout <<
"FontBox: " << front_container.size() << std::endl;
1368 for (front_iter = front_container.begin(); front_iter != front_container.end(); ++front_iter)
1369 std::cout << front_iter->first <<
"\t" << front_iter->second.size() << std::endl;
1371 std::cout <<
"EndBox: " << end_container.size() << std::endl;
1372 for (front_iter = end_container.begin(); front_iter != end_container.end(); ++front_iter)
1373 std::cout << front_iter->first <<
"\t" << front_iter->second.size() << std::endl;
1376 typename std::multimap<UInt, BoxElement>::iterator biter;
1380 for (front_iter = front_container.begin(); front_iter != front_container.end(); )
1382 best_match = end_container.end(); best_dist = INT_MAX;
1384 for (end_iter = end_container.begin(); end_iter != end_container.end(); ++end_iter)
1387 for (biter = front_iter->second.begin(); biter != front_iter->second.end(); ++biter)
1389 c = std::max(c, biter->second.c);
1392 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1393 std::cout <<
"Trying to match: " << end_iter->first <<
" to front " << front_iter->first << std::endl;
1396 c_dist = fabs(end_iter->first - front_iter->first);
1399 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1400 std::cout <<
"best match " << front_iter->second.begin()->first <<
"\t" << (--(end_iter->second.end()))->first << std::endl;
1402 if ((front_iter->second.begin()->first - (--(end_iter->second.end()))->first) <= RT_interleave + 1)
1405 best_match = end_iter;
1410 if (best_match == end_container.end())
1412 if (front_iter->second.size() >= RT_votes_cutoff)
1414 closed_boxes_.insert(*front_iter);
1421 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1422 std::cout <<
"Merging the boxes: " << front_iter->first <<
"\t" << best_match->first << std::endl;
1425 front_iter->second.insert(best_match->second.begin(), best_match->second.end());
1426 Box replacement(front_iter->second);
1429 DoubleReal c_mz = front_iter->first * (front_iter->second.size() - best_match->second.size()) + best_match->first * best_match->second.size();
1430 c_mz /= ((
DoubleReal) (front_iter->second.size()));
1432 help_iter = front_iter;
1434 std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1435 closed_boxes_.insert(help3);
1437 front_container.erase(front_iter);
1438 end_container.erase(best_match);
1439 front_iter = help_iter;
1444 for (end_iter = end_container.begin(); end_iter != end_container.end(); ++end_iter)
1446 if (end_iter->second.size() >= RT_votes_cutoff)
1448 closed_boxes_.insert(*end_iter);
1457 template <
typename PeakType>
1463 Int signal_size((
Int)candidate.size());
1468 Int p_h_ind = 1, end = 4 * (peak_cutoff - 1) - 1;
1470 std::vector<DoubleReal> positions(end);
1471 for (
Int i = 0; i < end; ++i)
1477 Int start_index = distance(candidate.begin(), candidate.
MZBegin(positions[0])) - 1;
1478 for (
Int v = 1; v <= end; ++v, ++p_h_ind)
1482 if (start_index < signal_size - 1)
1487 while (candidate[start_index].getMZ() < positions[v - 1]);
1489 if (start_index <= 0 || start_index >= signal_size - 1)
1494 c_left_iter2 = candidate.begin() + start_index - 1;
1495 c_right_iter2 = c_left_iter2 + 1;
1497 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());
1499 if (v == (
int)(ceil(end / 2.)))
1505 if (p_h_ind % 2 == 1)
1508 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1509 if (trunc(seed_mz) == 874)
1510 std::cout << -c_val << std::endl;
1516 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1517 if (trunc(seed_mz) == 874)
1518 std::cout << c_val << std::endl;
1523 start_index = distance(candidate.begin(), c_left_iter2);
1526 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1527 std::ofstream ofile_score(
"scores.dat", ios::app);
1528 std::ofstream ofile_check_score(
"check_scores.dat", ios::app);
1529 ofile_score.close();
1530 ofile_check_score.close();
1533 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1534 if (trunc(seed_mz) == 874)
1535 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;
1538 if (c_score - mid_val <= 0)
1543 if (c_score - mid_val <= ampl_cutoff)
1548 if (l_score <= 0 || c_score - l_score - mid_val <= 0)
1556 template <
typename PeakType>
1565 Int p_h_ind = 1, end = 4 * (peak_cutoff - 1) - 1;
1567 std::vector<DoubleReal> positions(end);
1568 for (
Int i = 0; i < end; ++i)
1574 Int start_index = distance(candidate.
begin(), candidate.
MZBegin(positions[0])) - 1;
1575 for (
Int v = 1; v <= end; ++v, ++p_h_ind)
1579 if (start_index < signal_size - 1)
1584 while (candidate.
getMZ(start_index) < positions[v - 1]);
1586 if (start_index <= 0 || start_index >= signal_size - 1)
1591 c_left_iter2 = candidate.
begin() + start_index - 1;
1592 c_right_iter2 = c_left_iter2 + 1;
1595 if (v == (
int)(ceil(end / 2.)))
1601 if (p_h_ind % 2 == 1)
1610 start_index = distance(candidate.
begin(), c_left_iter2);
1613 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1614 std::ofstream ofile_score(
"scores.dat", ios::app);
1615 std::ofstream ofile_check_score(
"check_scores.dat", ios::app);
1616 ofile_score << c_check_point <<
"\t" << c_score << std::endl;
1617 ofile_score.close();
1618 ofile_check_score.close();
1621 if (l_score <= 0 || c_score - l_score - mid_val <= 0 || c_score - mid_val <= ampl_cutoff)
1629 template <
typename PeakType>
1633 typename std::multimap<DoubleReal, Box>::iterator iter;
1634 typename Box::iterator box_iter;
1635 std::vector<BoxElement> final_box;
1636 DoubleReal c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1637 DoubleReal virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1639 typename std::pair<DoubleReal, DoubleReal> c_extend;
1640 for (iter = tmp_boxes_->at(c).begin(); iter != tmp_boxes_->at(c).end(); ++iter)
1643 Box& c_box = iter->second;
1644 av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
1645 virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
1648 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
1650 if (box_iter->second.score == 0)
1655 c_mz = box_iter->second.mz;
1656 virtual_av_intens += box_iter->second.intens;
1657 virtual_av_abs_intens += fabs(box_iter->second.intens);
1658 virtual_av_mz += c_mz * fabs(box_iter->second.intens);
1663 c_mz = box_iter->second.mz;
1664 av_score += box_iter->second.score;
1665 av_intens += box_iter->second.intens;
1666 av_abs_intens += fabs(box_iter->second.intens);
1667 av_mz += c_mz * fabs(box_iter->second.intens);
1674 av_intens = virtual_av_intens / virtual_count;
1676 av_mz = virtual_av_mz / virtual_av_abs_intens;
1682 av_mz /= av_abs_intens;
1686 c_box_element.
mz = av_mz;
1687 c_box_element.
c =
c;
1688 c_box_element.
score = av_score;
1689 c_box_element.
intens = av_intens;
1691 c_box_element.
RT = c_box.begin()->second.RT;
1692 final_box.push_back(c_box_element);
1695 Size num_o_feature = final_box.size();
1696 if (num_o_feature == 0)
1698 tmp_boxes_->at(c).clear();
1703 std::vector<DoubleReal> bwd_diffs(num_o_feature, 0);
1706 for (
Size i = 1; i < num_o_feature; ++i)
1708 bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].mz - final_box[i - 1].mz);
1711 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
1712 std::ofstream ofile_bwd(
"bwd_cpu.dat");
1713 for (
Size i = 0; i < num_o_feature; ++i)
1715 ofile_bwd << final_box[i].mz <<
"\t" << bwd_diffs[i] << std::endl;
1721 for (
Size i = 0; i < num_o_feature - 1; ++i)
1723 while (i < num_o_feature - 2)
1725 if (final_box[i].score > 0 || final_box[i].score == -1000)
1730 if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
1732 checkPositionForPlausibility_(candidate, ref, final_box[i].mz, final_box[i].c, scan_index, check_PPMs, final_box[i].intens, final_box[i].score);
1737 tmp_boxes_->at(c).clear();
1740 template <
typename PeakType>
1744 for (
UInt i = 0; i < scan.size(); ++i)
1746 if (scan[i].getIntensity() >= 0)
1748 av_intens += scan[i].getIntensity();
1751 return av_intens / (
double)scan.size();
1754 template <
typename PeakType>
1758 for (
UInt i = 0; i < scan.size(); ++i)
1760 if (scan[i].getIntensity() >= 0)
1762 intens = scan[i].getIntensity();
1763 res += (intens -
mean) * (intens - mean);
1766 return sqrt(res / (
double)(scan.size() - 1));
1769 template <
typename PeakType>
1772 std::vector<DoubleReal> diffs(scan.size() - 1, 0);
1773 for (
UInt i = 0; i < scan.size() - 1; ++i)
1775 diffs[i] = scan[i + 1].getMZ() - scan[i].getMZ();
1778 sort(diffs.begin(), diffs.end());
1780 for (
UInt i = 0; i < diffs.size() / 2; ++i)
1782 av_MZ_spacing += diffs[i];
1785 return av_MZ_spacing / (diffs.size() / 2);
1788 template <
typename PeakType>
1792 for (
UInt i = 0; i < scan.
size(); ++i)
1802 template <
typename PeakType>
1806 for (
UInt i = 0; i < scan.
size(); ++i)
1811 res += (intens -
mean) * (intens - mean);
1814 return sqrt(res / (
double)(scan.
size() - 1));
1817 template <
typename PeakType>
1823 typename std::multimap<DoubleReal, Box>::iterator upper_iter(open_boxes_.upper_bound(mz));
1824 typename std::multimap<DoubleReal, Box>::iterator lower_iter(open_boxes_.lower_bound(mz));
1826 if (lower_iter != open_boxes_.end())
1829 if (mz != lower_iter->first && lower_iter != open_boxes_.begin())
1835 typename std::multimap<DoubleReal, Box>::iterator insert_iter;
1836 bool create_new_box =
true;
1837 if (lower_iter == open_boxes_.end())
1842 if (!open_boxes_.empty())
1844 if (fabs((--lower_iter)->first - mz) < dist_constraint)
1846 create_new_box =
false;
1847 insert_iter = lower_iter;
1852 create_new_box =
true;
1857 if (upper_iter == open_boxes_.end() && fabs(lower_iter->first - mz) < dist_constraint)
1859 insert_iter = lower_iter;
1860 create_new_box =
false;
1864 create_new_box =
true;
1869 if (upper_iter != open_boxes_.end() && lower_iter != open_boxes_.end())
1874 DoubleReal dist_lower = fabs(lower_iter->first - mz);
1875 DoubleReal dist_upper = fabs(upper_iter->first - mz);
1876 dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1877 dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1879 if (dist_lower >= dist_constraint && dist_upper >= dist_constraint)
1881 create_new_box =
true;
1885 insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1886 create_new_box =
false;
1895 if (create_new_box ==
false)
1897 std::pair<UInt, BoxElement> help2(scan, element);
1898 insert_iter->second.insert(help2);
1901 Box replacement(insert_iter->second);
1905 DoubleReal c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
1906 c_mz /= ((
DoubleReal) insert_iter->second.size());
1909 open_boxes_.erase(insert_iter);
1910 std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
1911 open_boxes_.insert(help3);
1915 std::pair<UInt, BoxElement> help2(scan, element);
1916 std::multimap<UInt, BoxElement> help3;
1917 help3.insert(help2);
1918 std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help4(mz, help3);
1919 open_boxes_.insert(help4);
1923 template <
typename PeakType>
1929 std::multimap<DoubleReal, Box>& tmp_box(tmp_boxes_->at(c));
1930 typename std::multimap<DoubleReal, Box>::iterator upper_iter(tmp_box.upper_bound(mz));
1931 typename std::multimap<DoubleReal, Box>::iterator lower_iter(tmp_box.lower_bound(mz));
1933 if (lower_iter != tmp_box.end())
1936 if (mz != lower_iter->first && lower_iter != tmp_box.begin())
1942 typename std::multimap<DoubleReal, Box>::iterator insert_iter;
1943 bool create_new_box =
true;
1944 if (lower_iter == tmp_box.end())
1949 if (!tmp_box.empty())
1951 if (fabs((--lower_iter)->first - mz) < dist_constraint)
1953 create_new_box =
false;
1954 insert_iter = lower_iter;
1959 create_new_box =
true;
1964 if (upper_iter == tmp_box.end() && fabs(lower_iter->first - mz) < dist_constraint)
1966 insert_iter = lower_iter;
1967 create_new_box =
false;
1971 create_new_box =
true;
1976 if (upper_iter != tmp_box.end() && lower_iter != tmp_box.end())
1979 DoubleReal dist_lower = fabs(lower_iter->first - mz);
1980 DoubleReal dist_upper = fabs(upper_iter->first - mz);
1981 dist_lower = (dist_lower < dist_constraint) ? dist_lower : INT_MAX;
1982 dist_upper = (dist_upper < dist_constraint) ? dist_upper : INT_MAX;
1984 if (dist_lower >= dist_constraint && dist_upper >= dist_constraint)
1986 create_new_box =
true;
1990 insert_iter = (dist_lower < dist_upper) ? lower_iter : upper_iter;
1991 create_new_box =
false;
1999 if (create_new_box ==
false)
2001 std::pair<UInt, BoxElement> help2(scan, element);
2002 insert_iter->second.insert(help2);
2005 Box replacement(insert_iter->second);
2009 DoubleReal c_mz = insert_iter->first * (insert_iter->second.size() - 1) + mz;
2010 c_mz /= ((
DoubleReal) insert_iter->second.size());
2013 tmp_box.erase(insert_iter);
2014 std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help3(c_mz, replacement);
2015 tmp_box.insert(help3);
2019 std::pair<UInt, BoxElement> help2(scan, element);
2020 std::multimap<UInt, BoxElement> help3;
2021 help3.insert(help2);
2023 std::pair<DoubleReal, std::multimap<UInt, BoxElement> > help4(mz, help3);
2024 tmp_box.insert(help4);
2028 template <
typename PeakType>
2030 const UInt RT_votes_cutoff,
const Int front_bound,
const Int end_bound)
2032 typename std::multimap<DoubleReal, Box>::iterator iter, iter2;
2034 if ((
Int)scan_index == end_bound && end_bound != (
Int)map.
size() - 1)
2036 for (iter = open_boxes_.begin(); iter != open_boxes_.end(); ++iter)
2038 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2039 std::cout <<
"LOW THREAD insert in end_box " << iter->first << std::endl;
2040 typename Box::iterator dings;
2041 for (dings = iter->second.begin(); dings != iter->second.end(); ++dings)
2042 std::cout << map[dings->first].getRT() <<
"\t" << dings->second.c + 1 << std::endl;
2044 end_boxes_.insert(*iter);
2046 open_boxes_.
clear();
2050 for (iter = open_boxes_.begin(); iter != open_boxes_.end(); )
2052 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2053 if (front_bound > 0)
2055 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;
2060 UInt lastScan = (--(iter->second.end()))->first;
2061 if (scan_index - lastScan > RT_interleave + 1 || scan_index == map.
size() - 1)
2063 if (iter->second.begin()->first - front_bound <= RT_interleave + 1 && front_bound > 0)
2067 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2068 std::cout <<
"HIGH THREAD insert in front_box " << iter->first << std::endl;
2070 front_boxes_.insert(*iter);
2071 open_boxes_.erase(iter);
2081 if (iter->second.size() >= RT_votes_cutoff)
2085 closed_boxes_.insert(*(--iter));
2087 open_boxes_.erase(iter);
2097 template <
typename PeakType>
2100 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2101 std::cout <<
"**** CHECKING FOR BOX EXTENSIONS ****" << std::endl;
2105 typename Box::const_iterator iter;
2106 std::vector<DoubleReal> elution_profile(box.size());
2108 for (iter = box.begin(); iter != box.end(); ++iter, ++index)
2110 for (
Size i = iter->second.MZ_begin; i != iter->second.MZ_end; ++i)
2112 elution_profile[index] += map[iter->second.RT_index][i].getIntensity();
2114 elution_profile[index] /= iter->second.MZ_end - iter->second.MZ_begin + 1.;
2118 Int max_index = INT_MIN;
2119 for (
Size i = 0; i < elution_profile.size(); ++i)
2121 if (elution_profile[i] > max)
2124 max = elution_profile[i];
2128 Int max_extension = (
Int)(elution_profile.size()) - 2 * max_index;
2131 for (
Size i = 0; i < elution_profile.size(); ++i)
2133 av_elution += elution_profile[i];
2135 av_elution /= (
DoubleReal)elution_profile.size();
2138 for (
Size i = 0; i < elution_profile.size(); ++i)
2140 sd_elution += (av_elution - elution_profile[i]) * (av_elution - elution_profile[i]);
2142 sd_elution /= (
DoubleReal)(elution_profile.size() - 1);
2143 sd_elution = sqrt(sd_elution);
2147 for (iter = box.begin(); iter != box.end(); ++iter, ++index)
2149 av_mz += iter->second.mz;
2150 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2151 std::cout << iter->second.RT <<
"\t" << iter->second.mz <<
"\t" << iter->second.c + 1 << std::endl;
2158 if ((
Int)(box.begin()->second.RT_index) - 1 < 0)
2163 UInt pre_index = box.begin()->second.RT_index - 1;
2167 DoubleReal mz_start = map[pre_index + 1][box.begin()->second.MZ_begin].getMZ();
2168 DoubleReal mz_end = map[pre_index + 1][box.begin()->second.MZ_end].getMZ();
2173 pre_elution += mz_iter->getIntensity();
2178 if (pre_elution <= av_elution - 2 * sd_elution)
2183 Int c_index = max_extension;
2184 Int first_index = box.begin()->second.RT_index;
2185 for (
Int i = 1; i < max_extension; ++i)
2187 c_index = first_index - i;
2194 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2195 std::cout << box.begin()->second.RT <<
"\t" << av_mz <<
"\t" << box.begin()->second.c + 1 <<
"\t" <<
" extending the box " << std::endl;
2198 push2Box_(av_mz, c_index, box.begin()->second.c, box.begin()->second.score, c_iter->getIntensity(),
2199 map[c_index].getRT(), box.
begin()->second.MZ_begin, box.begin()->second.MZ_end);
2203 template <
typename PeakType>
2207 typename std::multimap<DoubleReal, Box>::iterator iter;
2208 typename Box::iterator box_iter;
2209 std::vector<BoxElement> final_box;
2210 DoubleReal c_mz, av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
2211 DoubleReal virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
2213 typename std::pair<DoubleReal, DoubleReal> c_extend;
2214 for (iter = tmp_boxes_->at(c).begin(); iter != tmp_boxes_->at(c).end(); ++iter)
2216 Box& c_box = iter->second;
2217 av_score = 0, av_mz = 0, av_intens = 0, av_abs_intens = 0, count = 0;
2218 virtual_av_mz = 0, virtual_av_intens = 0, virtual_av_abs_intens = 0, virtual_count = 0;
2220 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
2222 if (box_iter->second.score == 0)
2227 c_mz = box_iter->second.mz;
2228 virtual_av_intens += box_iter->second.intens;
2229 virtual_av_abs_intens += fabs(box_iter->second.intens);
2230 virtual_av_mz += c_mz * fabs(box_iter->second.intens);
2235 c_mz = box_iter->second.mz;
2236 av_score += box_iter->second.score;
2237 av_intens += box_iter->second.intens;
2238 av_abs_intens += fabs(box_iter->second.intens);
2239 av_mz += c_mz * fabs(box_iter->second.intens);
2246 av_intens = virtual_av_intens / virtual_count;
2248 av_mz = virtual_av_mz / virtual_av_abs_intens;
2254 av_mz /= av_abs_intens;
2258 c_box_element.
mz = av_mz;
2259 c_box_element.
c =
c;
2260 c_box_element.
score = av_score;
2261 c_box_element.
intens = av_intens;
2263 c_box_element.
RT = c_box.begin()->second.RT;
2265 final_box.push_back(c_box_element);
2268 UInt num_o_feature = final_box.size();
2269 if (num_o_feature == 0)
2271 tmp_boxes_->at(c).clear();
2276 std::vector<DoubleReal> bwd_diffs(num_o_feature, 0);
2279 for (
UInt i = 1; i < num_o_feature; ++i)
2281 bwd_diffs[i] = (final_box[i].intens - final_box[i - 1].intens) / (final_box[i].mz - final_box[i - 1].mz);
2284 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2285 std::ofstream ofile_bwd(
"bwd_gpu.dat");
2286 for (
UInt i = 0; i < num_o_feature; ++i)
2288 ofile_bwd << final_box[i].mz <<
"\t" << bwd_diffs[i] << std::endl;
2293 for (
UInt i = 0; i < num_o_feature - 1; ++i)
2295 while (i < num_o_feature - 2)
2297 if (final_box[i].score > 0 || final_box[i].score == -1000)
2302 if (bwd_diffs[i] > 0 && bwd_diffs[i + 1] < 0)
2304 checkPositionForPlausibility_(candidates, ref, final_box[i].mz, final_box[i].c, scan_index, check_PPMs, final_box[i].intens, final_box[i].score);
2309 tmp_boxes_->at(c).clear();
2312 template <
typename PeakType>
2316 typename std::multimap<DoubleReal, Box>::iterator iter;
2317 typename Box::iterator box_iter;
2319 DoubleReal av_intens = 0, av_ref_intens = 0, av_score = 0, av_mz = 0, av_RT = 0, mz_cutoff, sum_of_ref_intenses_g;
2320 bool restart =
false;
2322 typename std::pair<DoubleReal, DoubleReal> c_extend;
2323 for (iter = closed_boxes_.begin(); iter != closed_boxes_.end(); ++iter)
2325 sum_of_ref_intenses_g = 0;
2326 Box& c_box = iter->second;
2327 std::vector<DoubleReal> charge_votes(max_charge_, 0), charge_binary_votes(max_charge_, 0);
2332 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
2334 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2335 if (trunc(box_iter->second.mz) == 874)
2336 std::cout << box_iter->second.c <<
"\t" << box_iter->second.intens <<
"\t" << box_iter->second.score << std::endl;
2339 if (box_iter->second.score == -1000)
2345 charge_votes[box_iter->second.c] += box_iter->second.intens;
2346 ++charge_binary_votes[box_iter->second.c];
2355 best_charge_index = 0; best_charge_score = 0;
2356 for (
UInt i = 0; i < max_charge_; ++i)
2358 if (charge_votes[i] > best_charge_score)
2360 best_charge_index = i;
2361 best_charge_score = charge_votes[i];
2366 if (charge_binary_votes[best_charge_index] < RT_votes_cutoff && RT_votes_cutoff <= map.
size())
2371 c_charge = best_charge_index + 1;
2373 av_intens = 0, av_ref_intens = 0, av_score = 0, av_mz = 0, av_RT = 0;
2375 std::vector<DPosition<2> > point_set;
2377 for (box_iter = c_box.begin(); box_iter != c_box.end(); ++box_iter)
2379 sum_of_ref_intenses_l = 0;
2380 c_mz = box_iter->second.mz;
2381 c_RT = box_iter->second.RT;
2389 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2390 std::cout <<
"Intenstype: " << intenstype_ << std::endl;
2392 if (intenstype_ ==
"ref")
2397 for (
unsigned int i = 0; i < mz_cutoff; ++i)
2403 while (h_iter != c_spec.begin())
2406 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2407 if (trunc(c_mz) == 874)
2409 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;
2414 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2424 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2425 if (trunc(c_mz) == 874)
2430 sum_of_ref_intenses_l += hc_iter->getIntensity();
2431 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2432 if (trunc(c_mz) == 874)
2434 std::cout << sum_of_ref_intenses_l <<
"********" << std::endl;
2440 if (best_charge_index == box_iter->second.c)
2442 av_score += box_iter->second.score;
2443 av_intens += box_iter->second.intens;
2444 av_ref_intens += box_iter->second.ref_intens;
2445 sum_of_ref_intenses_g += sum_of_ref_intenses_l;
2446 av_mz += c_mz * box_iter->second.intens;
2452 av_ref_intens /= (
DoubleReal)charge_binary_votes[best_charge_index];
2453 av_score /= (
DoubleReal)charge_binary_votes[best_charge_index];
2456 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2457 if (trunc(av_mz) == 874)
2458 std::cout << av_mz <<
"\t" << best_charge_index <<
"\t" << best_charge_score << std::endl;
2465 c_feature.
setConvexHulls(std::vector<ConvexHull2D>(1, c_conv_hull));
2468 if (intenstype_ ==
"corrected")
2471 av_intens /= exp(-2 * lambda) * boost::math::cyl_bessel_i(0, 2 * lambda);
2473 if (intenstype_ ==
"ref")
2475 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2476 if (trunc(c_mz) == 874)
2478 std::cout << sum_of_ref_intenses_g <<
"####" << std::endl;
2482 av_intens = sum_of_ref_intenses_g;
2485 c_feature.
setMZ(av_mz);
2487 c_feature.
setRT(av_RT);
2489 feature_map.push_back(c_feature);
2495 template <
typename PeakType>
2503 iter = candidate.
MZBegin(seed_mz);
2505 if (iter == candidate.begin() || iter == candidate.end())
2510 std::pair<DoubleReal, DoubleReal> reals;
2511 ref_iter = ref.
MZBegin(seed_mz);
2516 reals = checkPPMTheoModel_(ref, iter->getMZ(),
c);
2517 real_mz = reals.first, real_intens = reals.second;
2521 while (h_iter != ref.begin())
2524 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2538 reals = checkPPMTheoModel_(ref, h_iter->getMZ(),
c);
2539 real_mz = reals.first, real_intens = reals.second;
2540 if (real_mz <= 0 || real_intens <= 0)
2544 real_mz = h_iter->getMZ();
2545 real_intens = h_iter->getIntensity();
2550 reals = std::pair<DoubleReal, DoubleReal>(seed_mz, ref_iter->getIntensity());
2551 real_mz = reals.first, real_intens = reals.second;
2553 if (real_mz <= 0 || real_intens <= 0)
2556 while (h_iter != ref.begin())
2559 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2573 real_mz = h_iter->getMZ(), real_intens = h_iter->getIntensity();
2574 if (real_mz <= 0 || real_intens <= 0)
2578 real_mz = h_iter->getMZ();
2579 real_intens = h_iter->getIntensity();
2583 DoubleReal c_score = scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2593 if (real_r_MZ_iter == ref.end())
2599 UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2600 UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2602 if (prev_score == -1000)
2604 push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2608 push2Box_(real_mz, scan_index, c, c_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2613 template <
typename PeakType>
2621 iter = candidate.
MZBegin(seed_mz);
2623 if (iter == candidate.
begin() || iter == candidate.
end())
2628 std::pair<DoubleReal, DoubleReal> reals;
2629 ref_iter = ref.
MZBegin(seed_mz);
2634 reals = checkPPMTheoModel_(ref, iter->getMZ(),
c);
2635 real_mz = reals.first, real_intens = reals.second;
2639 while (h_iter != ref.begin())
2642 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2657 reals = checkPPMTheoModel_(ref, h_iter->getMZ(),
c);
2658 real_mz = reals.first, real_intens = reals.second;
2660 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2661 std::cout <<
"Plausibility check old_mz: " << iter->getMZ() <<
"\t" << real_mz << std::endl;
2664 if (real_mz <= 0 || real_intens <= 0)
2668 real_mz = h_iter->getMZ();
2669 real_intens = h_iter->getIntensity();
2674 reals = std::pair<DoubleReal, DoubleReal>(seed_mz, ref_iter->getIntensity());
2675 real_mz = reals.first, real_intens = reals.second;
2677 if (real_mz <= 0 || real_intens <= 0)
2680 while (h_iter != ref.begin())
2683 if (h_iter->getIntensity() > hc_iter->getIntensity() || (h_iter->getIntensity() == hc_iter->getIntensity() && hc_iter->getIntensity() == 0))
2697 real_mz = h_iter->getMZ(), real_intens = h_iter->getIntensity();
2698 if (real_mz <= 0 || real_intens <= 0)
2702 real_mz = h_iter->getMZ();
2703 real_intens = h_iter->getIntensity();
2707 DoubleReal c_score = scoreThis_(candidate, peak_cutoff, real_mz, c, 0);
2717 if (real_r_MZ_iter == ref.end())
2723 UInt real_mz_begin = distance(ref.begin(), real_l_MZ_iter);
2724 UInt real_mz_end = distance(ref.begin(), real_r_MZ_iter);
2726 if (prev_score == -1000)
2728 push2Box_(real_mz, scan_index, c, prev_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2732 push2Box_(real_mz, scan_index, c, c_score, transintens, ref.
getRT(), real_mz_begin, real_mz_end, real_intens);
2738 template <
typename PeakType>
2742 DoubleReal ppms = getPPMs_(peptideMassRule_(mass), mass);
2745 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2746 std::cout << ::std::setprecision(8) << std::fixed << c_mz <<
"\t =(" <<
"ISO_WAVE" <<
")> " <<
"REJECT \t" << ppms <<
" (rule: " << peptideMassRule_(mass) <<
" got: " << mass <<
")" << std::endl;
2748 return std::pair<DoubleReal, DoubleReal>(-1, -1);
2751 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
2752 std::cout << ::std::setprecision(8) << std::fixed << c_mz <<
"\t =(" <<
"ISO_WAVE" <<
")> " <<
"ACCEPT \t" << ppms <<
" (rule: " << peptideMassRule_(mass) <<
" got: " << mass <<
")" << std::endl;
2754 return std::pair<DoubleReal, DoubleReal>(c_mz, ref.
MZBegin(c_mz)->getIntensity());
static DoubleReal getLambdaL(const DoubleReal m)
Returns the mass-parameter lambda (linear fit).
const int CUDA_INIT_SUCCESS
Definition: IsotopeWaveletConstants.h:92
float getMZ()
Definition: IsotopeWaveletTransform.h:72
A more convenient string class.
Definition: String.h:56
bool intensityPointerComparator(PeakType *a, PeakType *b)
Definition: IsotopeWaveletTransform.h:585
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:55
Size size() const
Definition: MSExperiment.h:117
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
static UInt getMzPeakCutOffAtMonoPos(const DoubleReal mass, const UInt z)
IntensityType getIntensity() const
Definition: Peak2D.h:161
bool myCudaComparator(const cudaHelp &a, const cudaHelp &b)
const double IW_PROTON_MASS
Definition: IsotopeWaveletConstants.h:70
A container for features.
Definition: FeatureMap.h:111
void sortByIntensity(bool reverse=false)
Sorting.
Definition: ConstRefVector.h:753
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
#define NULL
Definition: IsotopeWaveletParallelFor.h:41
const int CUDA_MIN_SORT_SIZE
Definition: IsotopeWaveletConstants.h:108
Iterator begin()
Definition: MSExperiment.h:147
Iterator MZEnd(CoordinateType mz)
Binary search for peak range end (returns the past-the-end iterator)
Definition: MSSpectrum.h:530
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:151
static IsotopeWavelet * init(const DoubleReal max_m, const UInt max_charge)
A 2-dimensional hull representation in [counter]clockwise direction - depending on axis labelling...
Definition: ConvexHull2D.h:75
Iterator MZBegin(CoordinateType mz)
Binary search for peak range begin.
Definition: MSSpectrum.h:506
bool intensityComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:573
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
void setRT(CoordinateType coordinate)
Mutable access to the RT coordinate (index 0)
Definition: Peak2D.h:209
const int CUDA_INIT_FAIL
Definition: IsotopeWaveletConstants.h:93
static DoubleReal mean(IteratorType begin, IteratorType end)
Calculates the mean of a range of values.
Definition: StatisticFunctions.h:69
An LC-MS feature.
Definition: Feature.h:66
int sortOnDevice(float *array, int *pos_indices, int numElements, int padding)
bool positionComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:591
const double IW_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:56
CoordinateType getMZ() const
Returns the m/z coordinate (index 1)
Definition: Peak2D.h:191
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
void setConvexHulls(const std::vector< ConvexHull2D > &hulls)
Set the convex hulls of single mass traces.
const double PEPTIDE_MASS_RULE_BOUND
Definition: IsotopeWaveletConstants.h:52
void setOverallQuality(QualityType q)
Set the overall quality.
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
bool intensityAscendingComparator(const PeakType &a, const PeakType &b)
Definition: IsotopeWaveletTransform.h:579
const int CUDA_EXTENDED_BLOCK_SIZE_MAX
Definition: IsotopeWaveletConstants.h:97
void clear(bool clear_meta_data)
Clears all data and meta data.
Definition: MSExperiment.h:809
float mz
Definition: IsotopeWaveletTransform.h:77
float score
Definition: IsotopeWaveletTransform.h:79
void setCharge(const ChargeType &ch)
Set charge state.
void addPoints(const PointArrayType &points)
const int CUDA_BLOCK_SIZE_MAX
Definition: IsotopeWaveletConstants.h:96
DoubleReal getRT() const
Definition: MSSpectrum.h:215
This vector holds pointer to the elements of another container.
Definition: ConstRefVector.h:72
const double IW_HALF_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:57
float intens
Definition: IsotopeWaveletTransform.h:78
const double PEPTIDE_MASS_RULE_FACTOR
Definition: IsotopeWaveletConstants.h:51
float getIntensity()
Definition: IsotopeWaveletTransform.h:74
void scoreOnDevice(int *sorted_positions_indices, float *trans_intensities, float *pos, float *scores, const int c, const int num_of_scores, const int overall_size, const unsigned int max_peak_cutoff, const float ampl_cutoff)
int Int
Signed integer type.
Definition: Types.h:100
const double IW_QUARTER_NEUTRON_MASS
Definition: IsotopeWaveletConstants.h:58
const unsigned int DEFAULT_NUM_OF_INTERPOLATION_POINTS
Definition: IsotopeWaveletConstants.h:45
void getExternalCudaTransforms(dim3 dimGrid, dim3 dimBlock, float *positions_dev, float *intensities_dev, int from_max_to_left, int from_max_to_right, float *result_dev, const int charge, const int to_load, const int to_compute, const int size, float *fwd2, bool highres=false)
static UInt getNumPeakCutOff(const DoubleReal mass, const UInt z)
static DoubleReal getValueByLambda(const DoubleReal lambda, const DoubleReal 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.
double DoubleReal
Double-precision real type.
Definition: Types.h:118
An internally used class, subsuming several variables.
Definition: IsotopeWaveletTransform.h:69
const double PEPTIDE_MASS_RULE_THEO_PPM_BOUND
Definition: IsotopeWaveletConstants.h:53