35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMISOTOPEWAVELET_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMISOTOPEWAVELET_H
49 #include <tbb/task_scheduler_init.h>
50 #include <tbb/pipeline.h>
51 #include <tbb/parallel_for.h>
75 template <
typename PeakType,
typename FeatureType>
91 this->
defaults_.
setValue(
"max_charge", 3,
"The maximal charge state to be considered.");
94 this->
defaults_.
setValue(
"intensity_threshold", -1.,
"The final threshold t' is build upon the formula: t' = av+t*sd, "
95 "where t is the intensity_threshold, av the average intensity within the wavelet transformed signal "
96 "and sd the standard deviation of the transform. "
97 "If you set intensity_threshold=-1, t' will be zero.\n"
98 "As the 'optimal' value for this parameter is highly data dependent, we would recommend to start "
99 "with -1, which will also extract features with very low signal-to-noise ratio. Subsequently, one "
100 "might increase the threshold to find an optimized trade-off between false positives and true positives. "
101 "Depending on the dynamic range of your spectra, suitable value ranges include: -1, [0:10], and if your data "
102 "features even very high intensity values, t can also adopt values up to around 30. "
103 "Please note that this parameter is not of an integer type, s.t. you can also use t:=0.1, e.g.");
104 this->
defaults_.
setValue(
"intensity_type",
"ref",
"Determines the intensity type returned for the identified features. 'ref' (default) returns the sum of the intensities of each isotopic peak within an isotope pattern. 'trans' refers to the intensity of the monoisotopic peak within the wavelet transform. 'corrected' refers also to the transformed intensity with an attempt to remove the effects of the convolution. While the latter ones might be preferable for qualitative analyses, 'ref' might be the best option to obtain quantitative results. Please note that intensity values might be spoiled (in particular for the option 'ref'), as soon as patterns overlap (see also the explanations given in the class documentation of FeatureFinderAlgorihtmIsotopeWavelet).",
StringList::create(
"advanced"));
107 this->
defaults_.
setValue(
"check_ppm",
"false",
"Enables/disables a ppm test vs. the averagine model, i.e. "
108 "potential peptide masses are checked for plausibility. In addition, a heuristic correcting potential mass shifts induced by the wavelet is applied.",
StringList::create(
"advanced"));
111 this->
defaults_.
setValue(
"hr_data",
"false",
"Must be true in case of high-resolution data, i.e. "
112 "for spectra featuring large m/z-gaps (present in FTICR and Orbitrap data, e.g.). Please check "
113 "a single MS scan out of your recording, if you are unsure.");
116 #if (defined(OPENMS_HAS_CUDA) || defined(OPENMS_HAS_TBB))
117 this->
defaults_.
setValue(
"parallel:use_gpus",
"-1",
"A comma-separated list of IDs corresponding to the GPU devices to use.\n"
118 "'-1' disables parallelization (CUDA/TBB) at all.\n");
121 this->
defaults_.
setValue(
"sweep_line:rt_votes_cutoff", 5,
"Defines the minimum number of "
122 "subsequent scans where a pattern must occur to be considered as a feature.",
StringList::create(
"advanced"));
124 this->
defaults_.
setValue(
"sweep_line:rt_interleave", 1,
"Defines the maximum number of "
125 "scans (w.r.t. rt_votes_cutoff) where an expected pattern is missing. There is usually no reason to change the default value.",
StringList::create(
"advanced"));
140 iwt->initializeScanCuda(*new_spec);
151 for (
UInt j = 0; j < spec.size() - 1; ++j)
153 spec[j].setMZ(-1 * (specr[j + 1].getMZ() - specr[j].getMZ()));
154 spec[j].setIntensity((specr[j].getIntensity() + specr[j + 1].getIntensity()));
156 spec[spec.size() - 1].setMZ(-1); spec[spec.size() - 1].setIntensity(-1);
162 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
163 std::ofstream ofilex(
"spacings.trans");
164 for (
UInt j = 0; j < spec.size() - 1; ++j)
166 ofilex << ::std::setprecision(12) << std::fixed << spec[j].getMZ() <<
"\t" << spec[j].getIntensity() << std::endl;
172 while (c_sorted_spec[pos].getIntensity() <= 0)
174 if (++pos >= c_sorted_spec.size())
176 std::cout <<
"Detected empty scan or a scan that cannot be interpolated with zeros in HR mode. " << std::endl;
177 std::cout <<
"Please check scan # " << i <<
" of your data set." << std::endl;
181 DoubleReal bound = -1 * c_sorted_spec[pos].getMZ();
191 new_spec->reserve(200000);
192 new_spec->
setRT(((*this->
map_)[i]).getRT());
194 new_spec->push_back(p);
197 for (
UInt j = 0; j < spec.size() - 1; ++j)
200 while (-spec[j].getMZ() - count * bound > bound)
204 new_spec->push_back(p);
207 new_spec->push_back(p);
210 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
211 std::ofstream ofiley(
"new_spec.trans");
212 for (
UInt j = 0; j < new_spec->size(); ++j)
214 ofiley << ::std::setprecision(12) << std::fixed << (*new_spec)[j].getMZ() <<
"\t" << (*new_spec)[j].getIntensity() << std::endl;
229 #ifdef OPENMS_HAS_CUDA
234 max_size = std::max(max_size, (*this->
map_)[i].size());
253 #if defined(OPENMS_HAS_TBB) && defined(OPENMS_HAS_CUDA)
257 tbb::task_scheduler_init init(num_gpus);
258 std::vector<IsotopeWaveletTransform<PeakType> *> iwts(num_gpus);
260 for (
UInt t = 0; t < num_gpus; ++t)
265 static tbb::affinity_partitioner ap;
269 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
270 std::cout <<
"Merging."; std::cout.flush();
273 for (
UInt t = 1; t < num_gpus; ++t)
278 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
279 std::cout <<
"Final mapping."; std::cout.flush();
283 for (
UInt t = 0; t < num_gpus; ++t)
291 std::cerr <<
"Error: You requested computation via TBB, but OpenMS has not been configured for TBB usage." << std::endl;
292 std::cerr <<
"Error: You need to rebuild OpenMS using the configure flag \"--enable-tbb-release\" or \"--enable-tbb-debug\"." << std::endl;
293 std::cerr <<
"Error: Please note that the multithreaded FeatureFinder needs necessarily the CUDA library, which must be enabled with \"--enable-cuda\"." << std::endl;
300 #ifdef OPENMS_HAS_CUDA
304 std::cout <<
"Using device with ID: " <<
gpu_ids_[0] << std::endl;
305 cudaDeviceProp props;
306 cudaGetDeviceProperties(&props,
gpu_ids_[0]);
307 std::cout <<
"This device is named: " << props.name << std::endl;
315 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
316 std::cout << ::std::fixed << ::std::setprecision(6) <<
"Spectrum " << i + 1 <<
" (" << (*this->
map_)[i].getRT() <<
") of " << this->
map_->
size() <<
" ... ";
320 if (c_ref.size() <= 1)
322 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
323 std::cout <<
"scan empty or consisting of a single data point. Skipping." << std::endl;
340 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
341 std::stringstream stream;
342 stream <<
"cpu_lowres_" << c_ref.
getRT() <<
"_" <<
c + 1 <<
".trans\0";
343 std::ofstream ofile(stream.str().c_str());
344 for (
UInt k = 0;
k < c_ref.size(); ++
k)
346 ofile << ::std::setprecision(8) << std::fixed << c_trans[
k].getMZ() <<
"\t" << c_trans[
k].getIntensity() <<
"\t" << c_ref[
k].getIntensity() << std::endl;
351 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
352 std::cout <<
"transform O.K. ... "; std::cout.flush();
358 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
359 std::cout <<
"charge recognition O.K. ... "; std::cout.flush();
375 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
376 std::stringstream stream;
377 stream <<
"cpu_highres_" << new_spec->
getRT() <<
"_" <<
c + 1 <<
".trans\0";
378 std::ofstream ofile(stream.str().c_str());
379 for (
UInt k = 0;
k < new_spec->size(); ++
k)
381 ofile << ::std::setprecision(8) << std::fixed << c_trans[
k].getMZ() <<
"\t" << c_trans[
k].getIntensity() <<
"\t" << (*new_spec)[
k].getIntensity() << std::endl;
386 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
387 std::cout <<
"transform O.K. ... "; std::cout.flush();
393 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
394 std::cout <<
"charge recognition O.K. ... "; std::cout.flush();
398 delete (new_spec); new_spec =
NULL;
404 #ifdef OPENMS_HAS_CUDA
416 iwt->getTransformCuda(*c_trans,
c);
418 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
419 std::stringstream stream;
420 stream <<
"gpu_lowres_" << ((*this->
map_)[i]).getRT() <<
"_" <<
c + 1 <<
".trans\0";
421 std::ofstream ofile(stream.str().c_str());
429 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
430 std::cout <<
"cuda transform for charge " <<
c + 1 <<
" O.K. ... "; std::cout.flush();
436 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
437 std::cout <<
"cuda charge recognition for charge " << c + 1 <<
" O.K." << std::endl;
441 iwt->finalizeScanCuda();
445 std::cout <<
"Warning/Error generated at scan " << i <<
" (" << ((*this->
map_)[i]).getRT() <<
")." << std::endl;
453 iwt->getTransformCuda(*c_trans,
c);
455 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
456 std::stringstream stream;
457 stream <<
"gpu_highres_" << ((*this->
map_)[i]).getRT() <<
"_" <<
c + 1 <<
".trans\0";
458 std::ofstream ofile(stream.str().c_str());
466 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
467 std::cout <<
"cuda transform for charge " <<
c + 1 <<
" O.K. ... "; std::cout.flush();
473 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
474 std::cout <<
"cuda charge recognition for charge " << c + 1 <<
" O.K." << std::endl;
479 iwt->finalizeScanCuda();
482 delete (new_spec); new_spec =
NULL;
483 delete (c_trans); c_trans =
NULL;
486 std::cerr <<
"Error: You requested computation on GPU, but OpenMS has not been configured for CUDA usage." << std::endl;
487 std::cerr <<
"Error: You need to rebuild OpenMS using the configure flag \"--enable-cuda\"." << std::endl;
492 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
493 std::cout <<
"updated box states." << std::endl;
504 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
505 std::cout <<
"Final mapping."; std::cout.flush();
513 #ifndef OPENMS_HAS_TBB
514 std::cerr <<
"Error: You requested multi-threaded computation via threading building blocks, but OpenMS has not been configured for TBB usage." << std::endl;
515 std::cerr <<
"Error: You need to rebuild OpenMS with -DENABLE_TBB=ON." << std::endl;
522 return "isotope_wavelet";
542 typedef std::map<UInt, BoxElement>
Box;
551 #if defined(OPENMS_HAS_TBB) && defined(OPENMS_HAS_CUDA)
553 Int device_num_, gpu_to_exclude_;
568 #if defined(OPENMS_HAS_CUDA) || defined(OPENMS_HAS_TBB)
570 std::vector<String> tokens;
577 if (tokens[0].trim().toInt() == -1)
583 gpu_ids_.push_back(tokens[0].trim().toInt());
586 for (
UInt i = 1; i < tokens.size(); ++i)
588 if (tokens[i].trim().toInt() == (
Int)
gpu_ids_[i - 1])
592 gpu_ids_.push_back(tokens[i].trim().toInt());
const int CUDA_INIT_SUCCESS
Definition: IsotopeWaveletConstants.h:92
FeatureMapType * features_
Output data pointer.
Definition: FeatureFinderAlgorithm.h:144
Param defaults_
Container for default parameters. This member should be filled in the constructor of derived classes!...
Definition: DefaultParamHandler.h:155
void setValue(const String &key, const DataValue &value, const String &description="", const StringList &tags=StringList())
Sets a value.
A more convenient string class.
Definition: String.h:56
FeatureFinder * ff_
Pointer to the calling FeatureFinder that is used to access the feature flags.
Definition: FeatureFinderAlgorithm.h:147
Size size() const
Definition: MSExperiment.h:117
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
String use_gpus_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:547
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:148
bool use_tbb_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:548
const MapType * map_
Input data pointer.
Definition: FeatureFinderAlgorithm.h:141
FeatureFinderAlgorithm< PeakType, FeatureType > Base
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:85
static void setMaxCharge(const UInt max_charge)
Sets the max_charge parameter.
Definition: IsotopeWavelet.h:116
DoubleReal score
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:537
#define NULL
Definition: IsotopeWaveletParallelFor.h:41
std::map< UInt, BoxElement > Box
Key: RT (index), value: BoxElement.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:542
Abstract base class for FeatureFinder algorithms.
Definition: FeatureFinderAlgorithm.h:74
UInt RT_interleave_
The number of subsequent scans a pattern must cover in order to be considered as signal.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:546
static const String getProductName()
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:520
UInt c
Note, this is not the charge (it is charge-1!!!)
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:536
DoubleReal mz
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:535
const PositionType & getMax() const
Returns the maximum position.
Definition: RangeManager.h:115
DoubleReal intens
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:538
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
Internally used data structure for the sweep line algorithm.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:533
void endProgress() const
Ends the progress display.
const DataValue & getValue(const String &key) const
Returns a value of a parameter.
virtual ~FeatureFinderAlgorithmIsotopeWavelet()
Destructor.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:132
bool use_cuda_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:548
const PositionType & getMin() const
Returns the minimum position.
Definition: RangeManager.h:109
DoubleReal intensity_threshold_
The only parameter of the isotope wavelet.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:545
Implements the isotope wavelet feature finder.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:76
static StringList create(const String &list, const char splitter= ',')
Returns a list that is created by splitting the given (comma-separated) string (String are not trimme...
DoubleReal RT
The elution time (not the scan index)
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:539
MSSpectrum< PeakType > * createHRData(const UInt i)
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:145
void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:558
String intensity_type_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:547
void setValidStrings(const String &key, const std::vector< String > &strings)
Sets the valid strings for the parameter key.
IsotopeWaveletTransform< PeakType >::TransSpectrum * prepareHRDataCuda(const UInt i, IsotopeWaveletTransform< PeakType > *iwt)
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:136
Command line progress.
Definition: ProgressLogger.h:68
void setMinInt(const String &key, Int min)
Sets the minimum value for the integer or integer list parameter key.
void setRT(DoubleReal rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:221
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void startProgress(SignedSize begin, SignedSize end, const String &label) const
Initializes the progress display.
bool hr_data_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:548
A class for distributing the data over several GPUs using Intel Threading Building Blocks...
Definition: IsotopeWaveletParallelFor.h:52
Int progress_counter_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:555
UInt max_charge_
The maximal charge state we will consider.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:544
void setProgress(SignedSize value) const
Sets the current progress.
FeatureFinderAlgorithmIsotopeWavelet()
Default Constructor.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:89
DoubleReal getRT() const
Definition: MSSpectrum.h:215
This vector holds pointer to the elements of another container.
Definition: ConstRefVector.h:72
void sortByPosition()
Lexicographically sorts the elements by their position.
Definition: ConstRefVector.h:766
UInt RT_votes_cutoff_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:546
int Int
Signed integer type.
Definition: Types.h:100
bool split(const char splitter, std::vector< String > &substrings, bool quote_protect=false) const
Splits a string into substrings using splitter as delimiter.
static FeatureFinderAlgorithm< PeakType, FeatureType > * create()
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:525
void run()
The working horse of this class.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:223
void setLogType(LogType type) const
Sets the progress log that should be used. The default type is NONE!
void defaultsToParam_()
Updates the parameters after the defaults have been set in the constructor.
bool check_PPMs_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:548
UInt real_RT_votes_cutoff_
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:546
std::vector< UInt > gpu_ids_
A list of all GPU devices that can be used.
Definition: FeatureFinderAlgorithmIsotopeWavelet.h:549