35 #ifndef OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
36 #define OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
60 #ifndef OPENMS_SYSTEM_STOPWATCH_H
63 #include <boost/math/special_functions/acosh.hpp>
64 #include <gsl/gsl_vector.h>
65 #include <gsl/gsl_multifit_nlin.h>
66 #include <gsl/gsl_blas.h>
109 tolerance_mz_ = tolerance_mz;
110 param_.setValue(
"2d:tolerance_mz", tolerance_mz);
118 max_peak_distance_ = max_peak_distance;
119 param_.setValue(
"2d:max_peak_distance", max_peak_distance);
128 param_.setValue(
"delta_abs_error", eps_abs);
137 param_.setValue(
"delta_rel_error", eps_rel);
145 max_iteration_ = max_iteration;
146 param_.setValue(
"iterations", max_iteration);
154 penalties_ = penalties;
155 param_.setValue(
"penalties:position", penalties.
pos);
156 param_.setValue(
"penalties:height", penalties.
height);
157 param_.setValue(
"penalties:left_width", penalties.
lWidth);
158 param_.setValue(
"penalties:right_width", penalties.
rWidth);
177 template <
typename InputSpectrumIterator,
typename OutputPeakType>
178 void optimize(InputSpectrumIterator first,
179 InputSpectrumIterator last,
187 std::vector<std::pair<SignedSize, SignedSize> >
signal2D;
199 std::multimap<DoubleReal, IsotopeCluster>
iso_map_;
202 std::multimap<DoubleReal, IsotopeCluster>::const_iterator
curr_region_;
236 static Int residual2D_(
const gsl_vector* x,
void* params, gsl_vector* f);
239 static Int jacobian2D_(
const gsl_vector* x,
void* params, gsl_matrix* J);
241 static Int evaluate2D_(
const gsl_vector* x,
void* params, gsl_vector* f, gsl_matrix* J);
248 std::vector<DoubleReal>::iterator searchInScan_(std::vector<DoubleReal>::iterator scan_begin,
249 std::vector<DoubleReal>::iterator scan_end,
253 template <
typename InputSpectrumIterator,
typename OutputPeakType>
254 void optimizeRegions_(InputSpectrumIterator& first,
255 InputSpectrumIterator& last,
259 template <
typename InputSpectrumIterator,
typename OutputPeakType>
260 void optimizeRegionsScanwise_(InputSpectrumIterator& first,
261 InputSpectrumIterator& last,
266 template <
typename InputSpectrumIterator,
typename OutputPeakType>
268 InputSpectrumIterator& first,
269 InputSpectrumIterator& last,
275 void findMatchingPeaks_(std::multimap<DoubleReal, IsotopeCluster>::iterator& it,
281 void updateMembers_();
285 template <
typename InputSpectrumIterator,
typename OutputPeakType>
290 if ((
UInt)distance(first, last) != ms_exp.
size())
292 throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
"Error in Two2Optimization: Raw and peak map do not have the same number of spectra");
300 for (
Size i = 0; i < ms_exp.
size(); ++i)
303 if (ms_exp[i].getFloatDataArrays().
size() < 6)
305 throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
"Error in Two2Optimization: Not enough meta data arrays present (1:area, 5:shape, 3:left width, 4:right width)");
307 bool area = ms_exp[i].getFloatDataArrays()[1].getName() ==
"maximumIntensity";
308 bool wleft = ms_exp[i].getFloatDataArrays()[3].getName() ==
"leftWidth";
309 bool wright = ms_exp[i].getFloatDataArrays()[4].getName() ==
"rightWidth";
310 bool shape = ms_exp[i].getFloatDataArrays()[5].getName() ==
"peakShape";
312 if (!area || !wleft || !wright || !shape)
314 throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__,
"Error in Two2Optimization: One or several meta data arrays missing (1:intensity, 5:shape, 3:left width, 4:right width)");
318 typedef typename InputSpectrumIterator::value_type InputSpectrumType;
319 typedef typename InputSpectrumType::value_type
PeakType;
326 std::cout <<
"empty experiment" << std::endl;
330 std::vector<DoubleReal> iso_last_scan;
331 std::vector<DoubleReal> iso_curr_scan;
332 std::vector<std::multimap<DoubleReal, IsotopeCluster>::iterator> clusters_last_scan;
333 std::vector<std::multimap<DoubleReal, IsotopeCluster>::iterator> clusters_curr_scan;
334 std::multimap<DoubleReal, IsotopeCluster>::iterator cluster_iter;
335 DoubleReal current_rt = ms_exp_it->getRT(), last_rt = 0;
341 UInt current_charge = 0;
345 for (
UInt curr_scan = 0; ms_exp_it + curr_scan != ms_exp_it_end; ++curr_scan)
347 Size nr_peaks_in_scan = (ms_exp_it + curr_scan)->size();
348 if (nr_peaks_in_scan == 0)
352 current_rt = (ms_exp_it + curr_scan)->getRT();
356 iso_last_scan = iso_curr_scan;
357 iso_curr_scan.clear();
358 clusters_last_scan = clusters_curr_scan;
359 clusters_curr_scan.clear();
362 std::cout <<
"Next scan with rt: " << current_rt << std::endl;
363 std::cout <<
"Next scan, rt = " << current_rt <<
" last_rt: " << last_rt << std::endl;
364 std::cout <<
"---------------------------------------------------------------------------" << std::endl;
370 ((lower_bound(first, last, s,
typename SpectrumType::RTLess()) - 1)->getRT() == last_rt))
374 for (
UInt curr_peak = 0; curr_peak < (ms_exp_it + curr_scan)->size() - 1; ++curr_peak)
378 DoubleReal curr_mz = (peak_it + curr_peak)->getMZ();
379 DoubleReal dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - curr_mz;
384 std::cout <<
"Isotopic pattern found ! " << std::endl;
385 std::cout <<
"We are at: " << (peak_it + curr_peak)->getMZ() <<
" " << curr_mz << std::endl;
387 if (!iso_last_scan.empty())
389 std::sort(iso_last_scan.begin(), iso_last_scan.end());
391 std::vector<DoubleReal>::iterator it =
392 searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);
396 if (delta_mz > tolerance_mz)
398 mz_in_hash = curr_mz;
406 new_cluster.
scans.push_back(curr_scan);
407 cluster_iter =
iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
415 cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];
418 if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
419 == cluster_iter->second.scans.end())
421 cluster_iter->second.scans.push_back(curr_scan);
438 mz_in_hash = curr_mz;
443 new_cluster.
scans.push_back(curr_scan);
444 cluster_iter =
iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
454 cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
456 iso_curr_scan.push_back(mz_in_hash);
457 clusters_curr_scan.push_back(cluster_iter);
460 cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
461 iso_curr_scan.push_back((peak_it + curr_peak)->getMZ());
462 clusters_curr_scan.push_back(cluster_iter);
465 if ((curr_peak + 1) >= nr_peaks_in_scan)
467 dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - (peak_it + curr_peak)->getMZ();
472 && curr_peak < (nr_peaks_in_scan - 1))
474 cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak + 1));
475 iso_curr_scan.push_back((peak_it + curr_peak + 1)->getMZ());
476 clusters_curr_scan.push_back(cluster_iter);
479 if (curr_peak >= nr_peaks_in_scan - 1)
481 dist2nextpeak = (peak_it + curr_peak + 1)->getMZ() - (peak_it + curr_peak)->getMZ();
491 if (!iso_last_scan.empty())
493 std::sort(iso_last_scan.begin(), iso_last_scan.end());
495 std::vector<DoubleReal>::iterator it =
496 searchInScan_(iso_last_scan.begin(), iso_last_scan.end(), curr_mz);
500 if (delta_mz > tolerance_mz)
502 mz_in_hash = curr_mz;
510 new_cluster.
scans.push_back(curr_scan);
511 cluster_iter =
iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
519 cluster_iter = clusters_last_scan[distance(iso_last_scan.begin(), it)];
522 if (find(cluster_iter->second.scans.begin(), cluster_iter->second.scans.end(), curr_scan)
523 == cluster_iter->second.scans.end())
525 cluster_iter->second.scans.push_back(curr_scan);
542 mz_in_hash = curr_mz;
547 new_cluster.
scans.push_back(curr_scan);
548 cluster_iter =
iso_map_.insert(std::pair<DoubleReal, IsotopeCluster>(mz_in_hash, new_cluster));
558 cluster_iter->second.peaks.insert(std::pair<UInt, UInt>(curr_scan, curr_peak));
560 iso_curr_scan.push_back(mz_in_hash);
561 clusters_curr_scan.push_back(cluster_iter);
569 last_rt = current_rt;
573 std::cout <<
iso_map_.size() <<
" isotopic clusters were found ! " << std::endl;
583 template <
typename InputSpectrumIterator,
typename OutputPeakType>
585 InputSpectrumIterator& last,
590 for (std::multimap<DoubleReal, IsotopeCluster>::iterator it =
iso_map_.begin();
595 std::cout <<
"element: " << counter << std::endl;
596 std::cout <<
"mz: " << it->first << std::endl <<
"rts: ";
598 std::cout << std::endl <<
"peaks: ";
599 IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
600 for (; iter != it->second.peaks.end(); ++iter)
601 std::cout << ms_exp[iter->first].getRT() <<
" " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
604 std::cout << std::endl << std::endl;
629 gsl_vector* start_value = gsl_vector_alloc(nr_parameters);
630 gsl_vector_set_zero(start_value);
633 std::map<Int, std::vector<PeakIndex> >::iterator m_peaks_it = d.
matching_peaks.begin();
634 DoubleReal av_mz = 0, av_lw = 0, av_rw = 0, avr_height = 0, height;
635 Int peak_counter = 0;
636 Int diff_peak_counter = 0;
640 av_mz = 0, av_lw = 0, av_rw = 0, avr_height = 0;
641 std::vector<PeakIndex>::iterator iter_iter = (m_peaks_it)->second.begin();
642 for (; iter_iter != m_peaks_it->second.end(); ++iter_iter)
644 height = ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[1][(iter_iter)->peak];
645 avr_height += height;
646 av_mz += (iter_iter)->getPeak(ms_exp).getMZ() * height;
647 av_lw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[3][(iter_iter)->peak] * height;
648 av_rw += ms_exp[(iter_iter)->spectrum].getFloatDataArrays()[4][(iter_iter)->peak] * height;
649 gsl_vector_set(start_value, peak_counter, height);
652 gsl_vector_set(start_value, d.
total_nr_peaks + 3 * diff_peak_counter, av_mz / avr_height);
653 gsl_vector_set(start_value, d.
total_nr_peaks + 3 * diff_peak_counter + 1, av_lw / avr_height);
654 gsl_vector_set(start_value, d.
total_nr_peaks + 3 * diff_peak_counter + 2, av_rw / avr_height);
659 std::cout <<
"----------------------------\n\nstart_value: " << std::endl;
660 for (
Size k = 0;
k < start_value->size; ++
k)
662 std::cout << gsl_vector_get(start_value,
k) << std::endl;
665 Int num_positions = 0;
670 std::cout << d.
signal2D[i + 1].second <<
" - " << d.
signal2D[i].second <<
" +1 " << std::endl;
675 std::cout <<
"num_positions : " << num_positions << std::endl;
679 gsl_multifit_function_fdf fit_function;
684 fit_function.n = std::max(num_positions + 1, (
Int)(nr_parameters));
685 fit_function.p = nr_parameters;
686 fit_function.params = &d;
688 std::cout <<
"fit_function.n " << fit_function.n
689 <<
"\tfit_function.p " << fit_function.p << std::endl;
691 const gsl_multifit_fdfsolver_type* type = gsl_multifit_fdfsolver_lmsder;
693 gsl_multifit_fdfsolver* fit = gsl_multifit_fdfsolver_alloc(type,
694 std::max(num_positions + 1, (
Int)(nr_parameters)),
697 gsl_multifit_fdfsolver_set(fit, &fit_function, start_value);
703 std::cout <<
"Before optimization: ||f|| = " << gsl_blas_dnrm2(fit->f) << std::endl;
712 status = gsl_multifit_fdfsolver_iterate(fit);
714 std::cout <<
"Iteration " << iteration <<
"; Status " << gsl_strerror(status) <<
"; " << std::endl;
715 std::cout <<
"||f|| = " << gsl_blas_dnrm2(fit->f) << std::endl;
716 std::cout <<
"Number of parms: " << nr_parameters << std::endl;
717 std::cout <<
"Delta: " << gsl_blas_dnrm2(fit->dx) << std::endl;
721 if (status != GSL_CONTINUE)
728 std::cout <<
"Finished! No. of iterations" << iteration << std::endl;
729 std::cout <<
"Delta: " << gsl_blas_dnrm2(fit->dx) << std::endl;
731 std::cout <<
"After optimization: || f || = " << gsl_blas_dnrm2(fit->f) << std::endl;
732 std::cout <<
"chisq/dof = " << pow(chi, 2.0) / (num_positions - nr_parameters);
735 std::cout <<
"----------------------------------------------\n\nnachher" << std::endl;
736 for (
Size k = 0;
k < fit->x->size; ++
k)
738 std::cout << gsl_vector_get(fit->x,
k) << std::endl;
742 std::map<Int, std::vector<PeakIndex> >::iterator itv
747 for (
Size j = 0; j < itv->second.size(); ++j)
751 std::cout <<
"pos: " << itv->second[j].getPeak(ms_exp).getMZ() <<
"\nint: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[1][itv->second[j].peak]
752 <<
"\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
753 <<
"\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] <<
"\n";
757 ms_exp[itv->second[j].spectrum][itv->second[j].peak].setMZ(mz);
758 DoubleReal height = (gsl_vector_get(fit->x, peak_idx));
759 ms_exp[itv->second[j].spectrum].getFloatDataArrays()[1][itv->second[j].peak] = height;
761 ms_exp[itv->second[j].spectrum].getFloatDataArrays()[3][itv->second[j].peak] = left_width;
763 ms_exp[itv->second[j].spectrum].getFloatDataArrays()[4][itv->second[j].peak] = right_width;
767 DoubleReal x_left_endpoint = mz - 1 / left_width* sqrt(height / 1 - 1);
768 DoubleReal x_rigth_endpoint = mz + 1 / right_width* sqrt(height / 1 - 1);
769 DoubleReal area_left = -height / left_width* atan(left_width * (x_left_endpoint - mz));
770 DoubleReal area_right = -height / right_width* atan(right_width * (mz - x_rigth_endpoint));
771 ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
775 DoubleReal x_left_endpoint = mz - 1 / left_width* boost::math::acosh(sqrt(height / 0.001));
776 DoubleReal x_rigth_endpoint = mz + 1 / right_width* boost::math::acosh(sqrt(height / 0.001));
777 DoubleReal area_left = -height / left_width * (sinh(left_width * (mz - x_left_endpoint)) / cosh(left_width * (mz - x_left_endpoint)));
778 DoubleReal area_right = -height / right_width * (sinh(right_width * (mz - x_rigth_endpoint)) / cosh(right_width * (mz - x_rigth_endpoint)));
779 ms_exp[itv->second[j].spectrum][itv->second[j].peak].setIntensity(area_left + area_right);
784 std::cout <<
"pos: " << itv->second[j].getPeak(ms_exp).getMZ() <<
"\nint: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[1][itv->second[j].peak]
785 <<
"\nlw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[3][itv->second[j].peak]
786 <<
"\nrw: " << itv->second[j].getSpectrum(ms_exp).getFloatDataArrays()[4][itv->second[j].peak] <<
"\n";
800 gsl_multifit_fdfsolver_free(fit);
801 gsl_vector_free(start_value);
807 template <
typename InputSpectrumIterator,
typename OutputPeakType>
809 InputSpectrumIterator& last,
822 if (dv.isEmpty() || dv.toString() ==
"")
825 penalties.
pos = (
float)dv;
828 if (dv.isEmpty() || dv.toString() ==
"")
831 penalties.
lWidth = (
float)dv;
834 if (dv.isEmpty() || dv.toString() ==
"")
837 penalties.
rWidth = (
float)dv;
839 std::cout << penalties.
pos <<
" "
840 << penalties.
rWidth <<
" "
841 << penalties.
lWidth << std::endl;
854 if (dv.isEmpty() || dv.toString() ==
"")
857 max_iteration = (
UInt)dv;
861 if (dv.isEmpty() || dv.toString() ==
"")
868 if (dv.isEmpty() || dv.toString() ==
"")
873 std::vector<PeakShape> peak_shapes;
877 for (std::multimap<DoubleReal, IsotopeCluster>::iterator it =
iso_map_.begin();
883 std::cerr <<
"element: " << counter << std::endl;
884 std::cerr <<
"mz: " << it->first << std::endl <<
"rts: ";
885 for (
Size i = 0; i < it->second.scans.size(); ++i)
886 std::cerr << it->second.scans[i] <<
"\n";
887 std::cerr << std::endl <<
"peaks: ";
888 IsotopeCluster::IndexSet::const_iterator iter = it->second.peaks.begin();
889 for (; iter != it->second.peaks.end(); ++iter)
890 std::cerr << ms_exp[iter->first].getRT() <<
" " << (ms_exp[iter->first][iter->second]).getMZ() << std::endl;
892 std::cerr << std::endl << std::endl;
913 data.
signal.reserve(size);
917 data.
positions.push_back(ms_it->getMZ());
918 data.
signal.push_back(ms_it->getIntensity());
924 pair.first = d.
iso_map_iter->second.peaks.begin()->first + idx;
926 IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(d.
iso_map_iter->second.peaks.begin(),
933 IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(d.
iso_map_iter->second.peaks.begin(),
937 while (set_iter != set_iter2)
939 const Size peak_index = set_iter->second;
942 spec[peak_index].getMZ(),
945 spec[peak_index].getIntensity(),
946 std::vector<Peak1D>::iterator(),
947 std::vector<Peak1D>::iterator(),
949 peak_shapes.push_back(shape);
957 OptimizePick opt(penalties, max_iteration, eps_abs, eps_rel);
959 std::cout <<
"vorher\n";
961 for (
Size p = 0; p < peak_shapes.size(); ++p)
963 std::cout << peak_shapes[p].mz_position <<
"\t" << peak_shapes[p].height
964 <<
"\t" << peak_shapes[p].left_width <<
"\t" << peak_shapes[p].right_width << std::endl;
969 std::cout <<
"nachher\n";
970 for (
Size p = 0; p < peak_shapes.size(); ++p)
972 std::cout << peak_shapes[p].mz_position <<
"\t" << peak_shapes[p].height
973 <<
"\t" << peak_shapes[p].left_width <<
"\t" << peak_shapes[p].right_width << std::endl;
977 pair.first = d.
iso_map_iter->second.peaks.begin()->first + idx;
979 set_iter = lower_bound(d.
iso_map_iter->second.peaks.begin(),
983 while (p < peak_shapes.size())
986 spec[set_iter->second].setMZ(peak_shapes[p].mz_position);
998 spec[set_iter->second].setIntensity(area_left + area_right);
1007 spec[set_iter->second].setIntensity(area_left + area_right);
1013 peak_shapes.clear();
1020 template <
typename InputSpectrumIterator,
typename OutputPeakType>
1022 InputSpectrumIterator& first,
1023 InputSpectrumIterator& last,
1029 typedef typename InputSpectrumIterator::value_type InputExperimentType;
1030 typedef typename InputExperimentType::value_type InputPeakType;
1031 typedef std::multimap<DoubleReal, IsotopeCluster> MapType;
1039 MapType::iterator iso_map_iter =
iso_map_.begin();
1040 for (
Size i = 0; i < iso_map_idx; ++i)
1044 std::cout <<
"rt begin: " << exp[iso_map_iter->second.scans[0]].getRT()
1045 <<
"\trt end: " << exp[iso_map_iter->second.scans[iso_map_iter->second.scans.size() - 1]].getRT()
1046 <<
" \t" << iso_map_iter->second.scans.
size() <<
" scans"
1051 for (
Size i = 0; i < iso_map_iter->second.scans.size(); ++i)
1056 rt = exp[iso_map_iter->second.scans[i]].getRT();
1062 std::cout << exp_it->getRT() <<
" vs " << iter->getRT() << std::endl;
1066 pair.first = iso_map_iter->second.peaks.
begin()->first + i;
1068 IsotopeCluster::IndexSet::const_iterator set_iter = lower_bound(iso_map_iter->second.peaks.begin(),
1069 iso_map_iter->second.peaks.end(),
1073 first_peak_mz = (exp_it->
begin() + set_iter->second)->getMZ() - 1;
1077 IsotopeCluster::IndexSet::const_iterator set_iter2 = lower_bound(iso_map_iter->second.peaks.begin(),
1078 iso_map_iter->second.peaks.end(),
1081 if (i == iso_map_iter->second.scans.size() - 1)
1083 set_iter2 = iso_map_iter->second.peaks.end();
1086 else if (set_iter2 != iso_map_iter->second.peaks.begin())
1089 last_peak_mz = (exp_it->
begin() + set_iter2->second)->getMZ() + 1;
1092 peak.setPosition(first_peak_mz);
1094 = lower_bound(iter->begin(), iter->end(), peak,
typename InputPeakType::PositionLess());
1095 if (raw_data_iter != iter->
begin())
1099 DoubleReal intensity = raw_data_iter->getIntensity();
1101 while (raw_data_iter != iter->
begin() && (raw_data_iter - 1)->getIntensity() < intensity &&
1102 (raw_data_iter - 1)->getIntensity() > noise_level)
1105 intensity = raw_data_iter->getIntensity();
1109 left.first = distance(first, iter);
1110 left.second = raw_data_iter - iter->
begin();
1112 std::cout <<
"left: " << iter->getRT() <<
"\t" << raw_data_iter->getMZ() << std::endl;
1115 peak.setPosition(last_peak_mz + 1);
1117 = upper_bound(iter->begin(), iter->end(), peak,
typename InputPeakType::PositionLess());
1118 if (raw_data_iter == iter->
end())
1120 intensity = raw_data_iter->getIntensity();
1122 while (raw_data_iter + 1 != iter->
end() && (raw_data_iter + 1)->getIntensity() < intensity)
1125 intensity = raw_data_iter->getIntensity();
1126 if ((raw_data_iter + 1 != iter->
end()) && (raw_data_iter + 1)->getIntensity() > noise_level)
1129 right.first = left.first;
1130 right.second = raw_data_iter - iter->
begin();
1132 std::cout <<
"right: " << iter->getRT() <<
"\t" << raw_data_iter->getMZ() << std::endl;
1140 std::cout << first_peak_mz <<
"\t" << last_peak_mz << std::endl;
1146 #endif //OPENMS_TRANSFORMATIONS_RAW2PEAK_TWODOPTIMIZATION_H
std::vector< double > positions
Positions and intensity values of the raw data.
Definition: OptimizePick.h:123
Internal representation of a peak shape (used by the PeakPickerCWT)
Definition: PeakShape.h:51
void getRegionEndpoints_(MSExperiment< OutputPeakType > &exp, InputSpectrumIterator &first, InputSpectrumIterator &last, Size iso_map_idx, DoubleReal noise_level, TwoDOptimization::Data &d)
Get the indices of the first and last raw data point of this region.
Definition: TwoDOptimization.h:1021
Definition: OptimizePick.h:120
std::map< Int, std::vector< PeakIndex > > matching_peaks_
Indices of peaks in the adjacent scans matching peaks in the scan with no. ref_scan.
Definition: TwoDOptimization.h:212
Size size() const
Definition: MSExperiment.h:117
void setMZTolerance(DoubleReal tolerance_mz)
Mutable access to the matching epsilon.
Definition: TwoDOptimization.h:107
std::pair< Size, Size > IndexPair
An index e.g. in an MSExperiment.
Definition: IsotopeCluster.h:48
std::vector< SpectrumType >::iterator Iterator
Mutable iterator.
Definition: MSExperiment.h:101
const OptimizationFunctions::PenaltyFactorsIntensity & getPenalties() const
Non-mutable access to the minimal number of adjacent scans.
Definition: TwoDOptimization.h:150
DoubleReal getMZTolerance() const
Non-mutable access to the matching epsilon.
Definition: TwoDOptimization.h:105
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
Comparison of mz_positions.
Definition: PeakShape.h:154
DoubleReal height
Maximum intensity of the peak shape.
Definition: PeakShape.h:131
DoubleReal eps_rel_
Convergence Parameter: Maximal relative error.
Definition: TwoDOptimization.h:219
Comparator for the retention time.
Definition: MSSpectrum.h:94
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
static Int residual2D_(const gsl_vector *x, void *params, gsl_vector *f)
Function computing estimated signal and its deviation to the experimental signal*/.
Peak2D PeakType
Definition: MassTrace.h:49
Helper struct (contains the size of an area and a raw data container)
Definition: TwoDOptimization.h:185
std::vector< DoubleReal >::iterator searchInScan_(std::vector< DoubleReal >::iterator scan_begin, std::vector< DoubleReal >::iterator scan_end, DoubleReal current_mz)
Iterator begin()
Definition: MSExperiment.h:147
virtual ~TwoDOptimization()
Destructor.
Definition: TwoDOptimization.h:98
std::multimap< DoubleReal, IsotopeCluster > iso_map_
stores the retention time of each isotopic cluster
Definition: TwoDOptimization.h:199
ContainerType::iterator Iterator
Mutable iterator.
Definition: MSSpectrum.h:123
double lWidth
Penalty factor for the peak shape's left width parameter.
Definition: OptimizePick.h:91
DoubleReal tolerance_mz_
threshold for the difference in the peak position of two matching peaks
Definition: TwoDOptimization.h:208
void setMaxRelError(DoubleReal eps_rel)
Mutable access to the maximal relative error.
Definition: TwoDOptimization.h:134
Class to hold strings, numeric values, lists of strings and lists of numeric values.
Definition: DataValue.h:57
void optimizeRegions_(InputSpectrumIterator &first, InputSpectrumIterator &last, MSExperiment< OutputPeakType > &ms_exp)
Definition: TwoDOptimization.h:584
std::multimap< DoubleReal, IsotopeCluster >::iterator iso_map_iter
Definition: TwoDOptimization.h:188
OptimizationFunctions::PenaltyFactorsIntensity penalties
Definition: TwoDOptimization.h:193
DoubleReal eps_abs_
Convergence Parameter: Maximal absolute error.
Definition: TwoDOptimization.h:216
DoubleReal getMaxPeakDistance() const
Non-mutable access to the maximal peak distance in a cluster.
Definition: TwoDOptimization.h:114
Iterator end()
Definition: MSExperiment.h:157
DoubleReal getMaxRelError() const
Non-mutable access to the maximal relative error.
Definition: TwoDOptimization.h:132
ConstIterator RTBegin(CoordinateType rt) const
Fast search for spectrum range begin.
Definition: MSExperiment.h:349
void setMaxAbsError(DoubleReal eps_abs)
Mutable access to the maximal absolute error.
Definition: TwoDOptimization.h:125
Base::iterator iterator
Definition: MSExperiment.h:114
void findMatchingPeaks_(std::multimap< DoubleReal, IsotopeCluster >::iterator &it, MSExperiment<> &ms_exp)
Identify matching peak in a peak cluster.
static Int jacobian2D_(const gsl_vector *x, void *params, gsl_matrix *J)
Function computing the Jacobian */.
Type
Peak shape type (asymmetric lorentzian or asymmetric hyperbolic secans squared).
Definition: PeakShape.h:68
void setPenalties(OptimizationFunctions::PenaltyFactorsIntensity &penalties)
Mutable access to the minimal number of adjacent scans.
Definition: TwoDOptimization.h:152
A method or algorithm argument contains illegal values.
Definition: Exception.h:634
ChargedIndexSet peaks
peaks in this cluster
Definition: IsotopeCluster.h:72
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
void optimize(std::vector< PeakShape > &peaks, Data &data)
Start the optimization of the peak shapes peaks. The original peak shapes will be subsituted by the o...
bool real_2D_
Optimization considering all scans of a cluster or optimization of each scan separately.
Definition: TwoDOptimization.h:225
DoubleReal max_peak_distance_
upper bound for distance between two peaks belonging to the same region
Definition: TwoDOptimization.h:205
Definition: PeakShape.h:70
Size total_nr_peaks
Definition: TwoDOptimization.h:189
std::vector< DoubleReal > signal
Definition: TwoDOptimization.h:195
void setMaxPeakDistance(DoubleReal max_peak_distance)
Mutable access to the maximal peak distance in a cluster.
Definition: TwoDOptimization.h:116
OptimizationFunctions::PenaltyFactorsIntensity penalties_
Penalty factors for some parameters in the optimization.
Definition: TwoDOptimization.h:229
std::vector< std::pair< SignedSize, SignedSize > > signal2D
Definition: TwoDOptimization.h:187
DoubleReal right_width
Right width parameter.
Definition: PeakShape.h:137
std::vector< DoubleReal > positions
Definition: TwoDOptimization.h:194
MSExperiment< Peak1D >::ConstIterator raw_data_first
Definition: TwoDOptimization.h:192
double pos
Penalty factor for the peak shape's position.
Definition: OptimizePick.h:89
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
std::multimap< DoubleReal, IsotopeCluster >::const_iterator curr_region_
Pointer to the current region.
Definition: TwoDOptimization.h:202
MSExperiment picked_peaks
Definition: TwoDOptimization.h:191
double rWidth
Penalty factor for the peak shape's right width parameter.
Definition: OptimizePick.h:93
void setRT(DoubleReal rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:221
std::vector< Size > scans
the scans of this cluster
Definition: IsotopeCluster.h:75
DoubleReal left_width
Left width parameter.
Definition: PeakShape.h:135
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
bool empty() const
Definition: MSExperiment.h:127
void setMaxIterations(UInt max_iteration)
Mutable access to the maximal number of iterations.
Definition: TwoDOptimization.h:143
Stores information about an isotopic cluster (i.e. potential peptide charge variants) ...
Definition: IsotopeCluster.h:45
Class for the penalty factors used during the optimization.
Definition: OptimizePick.h:71
Int charge
charge estimate (convention: zero means "no charge estimate")
Definition: IsotopeCluster.h:62
UInt max_iteration_
Convergence Parameter: Maximal number of iterations.
Definition: TwoDOptimization.h:222
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
std::vector< double > signal
Definition: OptimizePick.h:124
UInt getMaxIterations() const
Non-mutable access to the maximal number of iterations.
Definition: TwoDOptimization.h:141
DoubleReal getMaxAbsError() const
Non-mutable access to the maximal absolute error.
Definition: TwoDOptimization.h:123
std::map< Int, std::vector< PeakIndex > > matching_peaks
Definition: TwoDOptimization.h:190
Class for comparison of std::pair using first ONLY e.g. for use with std::sort.
Definition: ComparatorUtils.h:326
const FloatDataArrays & getFloatDataArrays() const
Returns a const reference to the float meta data arrays.
Definition: MSSpectrum.h:270
void optimizeRegionsScanwise_(InputSpectrumIterator &first, InputSpectrumIterator &last, MSExperiment< OutputPeakType > &ms_exp)
Definition: TwoDOptimization.h:808
DoubleReal height
Definition: OptimizePeakDeconvolution.h:85
int Int
Signed integer type.
Definition: Types.h:100
void optimize(InputSpectrumIterator first, InputSpectrumIterator last, MSExperiment< OutputPeakType > &ms_exp, bool real2D=true)
Find two dimensional peak clusters and optimize their peak parameters.
Definition: TwoDOptimization.h:286
static Int evaluate2D_(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
Function that calls residual2D and jacobian2D*/.
This class provides the non-linear optimization of the peak parameters.
Definition: OptimizePick.h:116
This class provides the two-dimensional optimization of the picked peak parameters.
Definition: TwoDOptimization.h:86
Class for the penalty factors used during the optimization.
Definition: OptimizePeakDeconvolution.h:66
double DoubleReal
Double-precision real type.
Definition: Types.h:118
DoubleReal mz_position
Centroid position.
Definition: PeakShape.h:133