Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages

FeatureFinderAlgorithmIsotopeWavelet.h

Go to the documentation of this file.
00001 // --------------------------------------------------------------------------
00002 //                   OpenMS -- Open-Source Mass Spectrometry
00003 // --------------------------------------------------------------------------
00004 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
00005 // ETH Zurich, and Freie Universitaet Berlin 2002-2012.
00006 //
00007 // This software is released under a three-clause BSD license:
00008 //  * Redistributions of source code must retain the above copyright
00009 //    notice, this list of conditions and the following disclaimer.
00010 //  * Redistributions in binary form must reproduce the above copyright
00011 //    notice, this list of conditions and the following disclaimer in the
00012 //    documentation and/or other materials provided with the distribution.
00013 //  * Neither the name of any author or any participating institution
00014 //    may be used to endorse or promote products derived from this software
00015 //    without specific prior written permission.
00016 // For a full list of authors, refer to the file AUTHORS.
00017 // --------------------------------------------------------------------------
00018 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00019 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00020 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00021 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
00022 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00023 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00024 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
00025 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
00026 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
00027 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
00028 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00029 //
00030 // --------------------------------------------------------------------------
00031 // $Maintainer: Rene Hussong $
00032 // $Authors: $
00033 // --------------------------------------------------------------------------
00034 
00035 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMISOTOPEWAVELET_H
00036 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_FEATUREFINDERALGORITHMISOTOPEWAVELET_H
00037 
00038 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/IsotopeWaveletTransform.h>
00039 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeatureFinderAlgorithm.h>
00040 #include <OpenMS/CONCEPT/ProgressLogger.h>
00041 #include <OpenMS/KERNEL/FeatureMap.h>
00042 #include <OpenMS/FORMAT/MzDataFile.h>
00043 #include <iostream>
00044 #include <time.h>
00045 #include <algorithm>
00046 
00047 #ifdef OPENMS_HAS_TBB
00048 #include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/IsotopeWaveletParallelFor.h>
00049 #include <tbb/task_scheduler_init.h>
00050 #include <tbb/pipeline.h>
00051 #include <tbb/parallel_for.h>
00052 #endif
00053 
00054 
00055 namespace OpenMS
00056 {
00075   template <typename PeakType, typename FeatureType>
00076   class FeatureFinderAlgorithmIsotopeWavelet :
00077     public FeatureFinderAlgorithm<PeakType, FeatureType>
00078   {
00079 #ifdef OPENMS_HAS_TBB
00080     friend class IsotopeWaveletParallelFor<PeakType, FeatureType>;
00081 #endif
00082 
00083 public:
00084 
00085     typedef FeatureFinderAlgorithm<PeakType, FeatureType> Base;
00086 
00087 
00089     FeatureFinderAlgorithmIsotopeWavelet()
00090     {
00091       this->defaults_.setValue("max_charge", 3, "The maximal charge state to be considered.");
00092       this->defaults_.setMinInt("max_charge", 1);
00093 
00094       this->defaults_.setValue("intensity_threshold", -1., "The final threshold t' is build upon the formula: t' = av+t*sd, "
00095                                                            "where t is the intensity_threshold, av the average intensity within the wavelet transformed signal "
00096                                                            "and sd the standard deviation of the transform. "
00097                                                            "If you set intensity_threshold=-1, t' will be zero.\n"
00098                                                            "As the 'optimal' value for this parameter is highly data dependent, we would recommend to start "
00099                                                            "with -1, which will also extract features with very low signal-to-noise ratio. Subsequently, one "
00100                                                            "might increase the threshold to find an optimized trade-off between false positives and true positives. "
00101                                                            "Depending on the dynamic range of your spectra, suitable value ranges include: -1, [0:10], and if your data "
00102                                                            "features even very high intensity values, t can also adopt values up to around 30. "
00103                                                            "Please note that this parameter is not of an integer type, s.t. you can also use t:=0.1, e.g.");
00104       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"));
00105       this->defaults_.setValidStrings("intensity_type", StringList::create("ref,trans,corrected"));
00106 
00107       this->defaults_.setValue("check_ppm", "false", "Enables/disables a ppm test vs. the averagine model, i.e. "
00108                                                      "potential peptide masses are checked for plausibility. In addition, a heuristic correcting potential mass shifts induced by the wavelet is applied.", StringList::create("advanced"));
00109       this->defaults_.setValidStrings("check_ppm", StringList::create("true,false"));
00110 
00111       this->defaults_.setValue("hr_data", "false", "Must be true in case of high-resolution data, i.e. "
00112                                                    "for spectra featuring large m/z-gaps (present in FTICR and Orbitrap data, e.g.). Please check "
00113                                                    "a single MS scan out of your recording, if you are unsure.");
00114       this->defaults_.setValidStrings("hr_data", StringList::create("true,false"));
00115 
00116 #if (defined(OPENMS_HAS_CUDA) || defined(OPENMS_HAS_TBB))
00117       this->defaults_.setValue("parallel:use_gpus", "-1", "A comma-separated list of IDs corresponding to the GPU devices to use.\n"
00118                                                           "'-1' disables parallelization (CUDA/TBB) at all.\n");
00119 #endif
00120 
00121       this->defaults_.setValue("sweep_line:rt_votes_cutoff", 5, "Defines the minimum number of "
00122                                                                 "subsequent scans where a pattern must occur to be considered as a feature.", StringList::create("advanced"));
00123       this->defaults_.setMinInt("sweep_line:rt_votes_cutoff", 0);
00124       this->defaults_.setValue("sweep_line:rt_interleave", 1, "Defines the maximum number of "
00125                                                               "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"));
00126       this->defaults_.setMinInt("sweep_line:rt_interleave", 0);
00127 
00128       this->defaultsToParam_();
00129     }
00130 
00132     virtual ~FeatureFinderAlgorithmIsotopeWavelet()
00133     {
00134     }
00135 
00136     typename IsotopeWaveletTransform<PeakType>::TransSpectrum * prepareHRDataCuda(const UInt i, IsotopeWaveletTransform<PeakType> * iwt)
00137     {
00138       MSSpectrum<PeakType> * new_spec(createHRData(i));
00139       typename IsotopeWaveletTransform<PeakType>::TransSpectrum * c_trans = new typename IsotopeWaveletTransform<PeakType>::TransSpectrum(new_spec);
00140       iwt->initializeScanCuda(*new_spec);
00141 
00142       return c_trans;
00143     }
00144 
00145     MSSpectrum<PeakType> * createHRData(const UInt i)
00146     {
00147       MSSpectrum<PeakType> spec((*this->map_)[i]);
00148 
00149       const MSSpectrum<PeakType> & specr((*this->map_)[i]);
00150 
00151       for (UInt j = 0; j < spec.size() - 1; ++j)
00152       {
00153         spec[j].setMZ(-1 * (specr[j + 1].getMZ() - specr[j].getMZ()));
00154         spec[j].setIntensity((specr[j].getIntensity() + specr[j + 1].getIntensity()));
00155       }
00156       spec[spec.size() - 1].setMZ(-1); spec[spec.size() - 1].setIntensity(-1);
00157 
00158       ConstRefVector<MSSpectrum<PeakType> > c_sorted_spec(spec.begin(), spec.end());
00159       //Sort in ascending order according to the intensities present in the transform
00160       c_sorted_spec.sortByPosition();
00161 
00162 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00163       std::ofstream ofilex("spacings.trans");
00164       for (UInt j = 0; j < spec.size() - 1; ++j)
00165       {
00166         ofilex << ::std::setprecision(12) << std::fixed << spec[j].getMZ() << "\t" << spec[j].getIntensity() << std::endl;
00167       }
00168       ofilex.close();
00169 #endif
00170 
00171       UInt pos = 0;
00172       while (c_sorted_spec[pos].getIntensity() <= 0)
00173       {
00174         if (++pos >= c_sorted_spec.size())
00175         {
00176           std::cout << "Detected empty scan or a scan that cannot be interpolated with zeros in HR mode. " << std::endl;
00177           std::cout << "Please check scan # " << i << " of your data set." << std::endl;
00178           exit(-1);
00179         }
00180       }
00181       DoubleReal bound = -1 * c_sorted_spec[pos].getMZ();
00182 
00183       if (bound > (1. / max_charge_) / 2.)
00184       {
00185         //that might be case for simulated spectra,
00186         //which might show a very artificial spacing
00187         bound = (1. / max_charge_) / 2. / 4.;
00188       }
00189 
00190       MSSpectrum<PeakType> * new_spec = new MSSpectrum<PeakType>;
00191       new_spec->reserve(200000);
00192       new_spec->setRT(((*this->map_)[i]).getRT());
00193       PeakType p; p.setMZ(specr[0].getMZ()); p.setIntensity(specr[0].getIntensity());
00194       new_spec->push_back(p);
00195 
00196       UInt count;
00197       for (UInt j = 0; j < spec.size() - 1; ++j)
00198       {
00199         count = 0;
00200         while (-spec[j].getMZ() - count * bound > bound)
00201         {
00202           ++count;
00203           p.setMZ(specr[j].getMZ() + count * bound); p.setIntensity(0);
00204           new_spec->push_back(p);
00205         }
00206         p.setMZ(specr[j + 1].getMZ()); p.setIntensity(specr[j + 1].getIntensity());
00207         new_spec->push_back(p);
00208       }
00209 
00210 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00211       std::ofstream ofiley("new_spec.trans");
00212       for (UInt j = 0; j < new_spec->size(); ++j)
00213       {
00214         ofiley << ::std::setprecision(12) << std::fixed << (*new_spec)[j].getMZ() << "\t" << (*new_spec)[j].getIntensity() << std::endl;
00215       }
00216       ofiley.close();
00217 #endif
00218 
00219       return new_spec;
00220     }
00221 
00223     void run()
00224     {
00225       DoubleReal max_mz = this->map_->getMax()[1];
00226       DoubleReal min_mz = this->map_->getMin()[1];
00227 
00228       Size max_size = 0;
00229 #ifdef OPENMS_HAS_CUDA
00230       if (use_cuda_)               //some preprocessing necessary for the GPU computation
00231       {
00232         for (UInt i = 0; i < this->map_->size(); ++i)
00233         {
00234           max_size = std::max(max_size, (*this->map_)[i].size());
00235         }
00236       }
00237 #endif
00238 
00239       //Check for useless RT_votes_cutoff_ parameter
00240       if (RT_votes_cutoff_ > this->map_->size())
00241       {
00242         real_RT_votes_cutoff_ = 0;
00243       }
00244       else
00245       {
00246         real_RT_votes_cutoff_ = RT_votes_cutoff_;
00247       }
00248 
00249       this->ff_->setLogType(ProgressLogger::CMD);
00250       progress_counter_ = 0;
00251       this->ff_->startProgress(0, 2 * this->map_->size() * max_charge_, "analyzing spectra");
00252 
00253 #if defined(OPENMS_HAS_TBB) && defined(OPENMS_HAS_CUDA)
00254       if (use_tbb_)
00255       {
00256         UInt num_gpus = this->gpu_ids_.size();
00257         tbb::task_scheduler_init init(num_gpus);
00258         std::vector<IsotopeWaveletTransform<PeakType> *> iwts(num_gpus);
00259 
00260         for (UInt t = 0; t < num_gpus; ++t)
00261         {
00262           iwts[t] = new IsotopeWaveletTransform<PeakType>(min_mz, max_mz, max_charge_, max_size, true, hr_data_, intensity_type_);
00263         }
00264 
00265         static tbb::affinity_partitioner ap;
00266         //The parallel execution over all available GPU devices
00267         tbb::parallel_for(tbb::blocked_range<size_t>(0, num_gpus, 1), IsotopeWaveletParallelFor<PeakType, FeatureType>(iwts, this), ap);
00268 
00269 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00270         std::cout << "Merging."; std::cout.flush();
00271 #endif
00272 
00273         for (UInt t = 1; t < num_gpus; ++t)
00274         {
00275           iwts[0]->mergeFeatures(iwts[t], RT_interleave_, RT_votes_cutoff_);
00276         }
00277 
00278 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00279         std::cout << "Final mapping."; std::cout.flush();
00280 #endif
00281         *this->features_ = iwts[0]->mapSeeds2Features(*this->map_, real_RT_votes_cutoff_);
00282 
00283         for (UInt t = 0; t < num_gpus; ++t)
00284         {
00285           delete (iwts[t]);
00286         }
00287       }
00288 #else
00289       if (use_tbb_)
00290       {
00291         std::cerr << "Error: You requested computation via TBB, but OpenMS has not been configured for TBB usage." << std::endl;
00292         std::cerr << "Error: You need to rebuild OpenMS using the configure flag \"--enable-tbb-release\" or \"--enable-tbb-debug\"." << std::endl;
00293         std::cerr << "Error: Please note that the multithreaded FeatureFinder needs necessarily the CUDA library, which must be enabled with \"--enable-cuda\"." << std::endl;
00294       }
00295 #endif
00296 
00297       if (!use_tbb_)
00298       {
00299         IsotopeWaveletTransform<PeakType> * iwt = new IsotopeWaveletTransform<PeakType>(min_mz, max_mz, max_charge_, max_size, use_cuda_, hr_data_, intensity_type_);
00300 #ifdef OPENMS_HAS_CUDA
00301         if (use_cuda_)
00302         {
00303           cudaSetDevice(gpu_ids_[0]);
00304           std::cout << "Using device with ID: " << gpu_ids_[0] << std::endl;
00305           cudaDeviceProp props;
00306           cudaGetDeviceProperties(&props, gpu_ids_[0]);
00307           std::cout << "This device is named: " << props.name << std::endl;
00308         }
00309 #endif
00310 
00311         for (UInt i = 0; i < this->map_->size(); ++i)
00312         {
00313           const MSSpectrum<PeakType> & c_ref((*this->map_)[i]);
00314 
00315 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00316           std::cout << ::std::fixed << ::std::setprecision(6) << "Spectrum " << i + 1 << " (" << (*this->map_)[i].getRT() << ") of " << this->map_->size() << " ... ";
00317           std::cout.flush();
00318 #endif
00319 
00320           if (c_ref.size() <= 1)               //unable to do transform anything
00321           {
00322 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00323             std::cout << "scan empty or consisting of a single data point. Skipping." << std::endl;
00324 #endif
00325             this->ff_->setProgress(progress_counter_ += 2);
00326             continue;
00327           }
00328 
00329           if (!use_cuda_)
00330           {
00331             if (!hr_data_)                 //Lowres data
00332             {
00333               iwt->initializeScan((*this->map_)[i]);
00334               for (UInt c = 0; c < max_charge_; ++c)
00335               {
00336                 MSSpectrum<PeakType> c_trans(c_ref);
00337 
00338                 iwt->getTransform(c_trans, c_ref, c);
00339 
00340 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00341                 std::stringstream stream;
00342                 stream << "cpu_lowres_" << c_ref.getRT() << "_" << c + 1 << ".trans\0";
00343                 std::ofstream ofile(stream.str().c_str());
00344                 for (UInt k = 0; k < c_ref.size(); ++k)
00345                 {
00346                   ofile << ::std::setprecision(8) << std::fixed << c_trans[k].getMZ() << "\t" << c_trans[k].getIntensity() << "\t" << c_ref[k].getIntensity() << std::endl;
00347                 }
00348                 ofile.close();
00349 #endif
00350 
00351 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00352                 std::cout << "transform O.K. ... "; std::cout.flush();
00353 #endif
00354                 this->ff_->setProgress(++progress_counter_);
00355 
00356                 iwt->identifyCharge(c_trans, c_ref, i, c, intensity_threshold_, check_PPMs_);
00357 
00358 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00359                 std::cout << "charge recognition O.K. ... "; std::cout.flush();
00360 #endif
00361                 this->ff_->setProgress(++progress_counter_);
00362               }
00363             }
00364             else                 //Highres data
00365             {
00366               MSSpectrum<PeakType> * new_spec(NULL);
00367               for (UInt c = 0; c < max_charge_; ++c)
00368               {
00369                 new_spec = createHRData(i);
00370                 iwt->initializeScan(*new_spec, c);
00371                 MSSpectrum<PeakType> c_trans(*new_spec);
00372 
00373                 iwt->getTransformHighRes(c_trans, *new_spec, c);
00374 
00375 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00376                 std::stringstream stream;
00377                 stream << "cpu_highres_" << new_spec->getRT() << "_" << c + 1 << ".trans\0";
00378                 std::ofstream ofile(stream.str().c_str());
00379                 for (UInt k = 0; k < new_spec->size(); ++k)
00380                 {
00381                   ofile << ::std::setprecision(8) << std::fixed << c_trans[k].getMZ() << "\t" << c_trans[k].getIntensity() << "\t" << (*new_spec)[k].getIntensity() << std::endl;
00382                 }
00383                 ofile.close();
00384 #endif
00385 
00386 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00387                 std::cout << "transform O.K. ... "; std::cout.flush();
00388 #endif
00389                 this->ff_->setProgress(++progress_counter_);
00390 
00391                 iwt->identifyCharge(c_trans, *new_spec, i, c, intensity_threshold_, check_PPMs_);
00392 
00393 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00394                 std::cout << "charge recognition O.K. ... "; std::cout.flush();
00395 #endif
00396                 this->ff_->setProgress(++progress_counter_);
00397 
00398                 delete (new_spec); new_spec = NULL;
00399               }
00400             }
00401           }
00402           else               //use CUDA but not TBB
00403           {
00404 #ifdef OPENMS_HAS_CUDA
00405             bool success = true;
00406             typename IsotopeWaveletTransform<PeakType>::TransSpectrum * c_trans(NULL); MSSpectrum<PeakType> * new_spec(NULL);
00407             if (!hr_data_)                     //LowRes data
00408             {
00409               c_trans = new typename IsotopeWaveletTransform<PeakType>::TransSpectrum(&(*this->map_)[i]);
00410               success = iwt->initializeScanCuda((*this->map_)[i]) == Constants::CUDA_INIT_SUCCESS;
00411 
00412               if (success)
00413               {
00414                 for (UInt c = 0; c < max_charge_; ++c)
00415                 {
00416                   iwt->getTransformCuda(*c_trans, c);
00417 
00418 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00419                   std::stringstream stream;
00420                   stream << "gpu_lowres_" << ((*this->map_)[i]).getRT() << "_" << c + 1 << ".trans\0";
00421                   std::ofstream ofile(stream.str().c_str());
00422                   for (UInt k = 0; k < c_trans->size(); ++k)
00423                   {
00424                     ofile << ::std::setprecision(8) << std::fixed << c_trans->getMZ(k) << "\t" <<  c_trans->getTransIntensity(k) << "\t" << c_trans->getRefIntensity(k) << std::endl;
00425                   }
00426                   ofile.close();
00427 #endif
00428 
00429 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00430                   std::cout << "cuda transform for charge " << c + 1 << "  O.K. ... "; std::cout.flush();
00431 #endif
00432                   this->ff_->setProgress(++progress_counter_);
00433 
00434                   iwt->identifyChargeCuda(*c_trans, i, c, intensity_threshold_, check_PPMs_);
00435 
00436 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00437                   std::cout << "cuda charge recognition for charge " << c + 1 << " O.K." << std::endl;
00438 #endif
00439                   this->ff_->setProgress(++progress_counter_);
00440                 }
00441                 iwt->finalizeScanCuda();
00442               }
00443               else
00444               {
00445                 std::cout << "Warning/Error generated at scan " << i << " (" << ((*this->map_)[i]).getRT() << ")." << std::endl;
00446               }
00447             }
00448             else                     //HighRes data
00449             {
00450               c_trans = prepareHRDataCuda(i, iwt);
00451               for (UInt c = 0; c < max_charge_; ++c)
00452               {
00453                 iwt->getTransformCuda(*c_trans, c);
00454 
00455 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00456                 std::stringstream stream;
00457                 stream << "gpu_highres_" << ((*this->map_)[i]).getRT() << "_" << c + 1 << ".trans\0";
00458                 std::ofstream ofile(stream.str().c_str());
00459                 for (UInt k = 0; k < c_trans->size(); ++k)
00460                 {
00461                   ofile << ::std::setprecision(8) << std::fixed << c_trans->getMZ(k) << "\t" <<  c_trans->getTransIntensity(k)  << "\t" << c_trans->getRefIntensity(k) << std::endl;
00462                 }
00463                 ofile.close();
00464 #endif
00465 
00466 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00467                 std::cout << "cuda transform for charge " << c + 1 << "  O.K. ... "; std::cout.flush();
00468 #endif
00469                 this->ff_->setProgress(++progress_counter_);
00470 
00471                 iwt->identifyChargeCuda(*c_trans, i, c, intensity_threshold_, check_PPMs_);
00472 
00473 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00474                 std::cout << "cuda charge recognition for charge " << c + 1 << " O.K." << std::endl;
00475 #endif
00476                 this->ff_->setProgress(++progress_counter_);
00477               }
00478               c_trans->destroy();
00479               iwt->finalizeScanCuda();
00480             }
00481 
00482             delete (new_spec); new_spec = NULL;
00483             delete (c_trans); c_trans = NULL;
00484 
00485 #else
00486             std::cerr << "Error: You requested computation on GPU, but OpenMS has not been configured for CUDA usage." << std::endl;
00487             std::cerr << "Error: You need to rebuild OpenMS using the configure flag \"--enable-cuda\"." << std::endl;
00488 #endif
00489           }
00490 
00491           iwt->updateBoxStates(*this->map_, i, RT_interleave_, real_RT_votes_cutoff_);
00492 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00493           std::cout << "updated box states." << std::endl;
00494 #endif
00495 
00496           std::cout.flush();
00497         }
00498 
00499         this->ff_->endProgress();
00500 
00501         //Forces to empty OpenBoxes_ and to synchronize ClosedBoxes_
00502         iwt->updateBoxStates(*this->map_, INT_MAX, RT_interleave_, real_RT_votes_cutoff_);
00503 
00504 #ifdef OPENMS_DEBUG_ISOTOPE_WAVELET
00505         std::cout << "Final mapping."; std::cout.flush();
00506 #endif
00507         *this->features_ = iwt->mapSeeds2Features(*this->map_, real_RT_votes_cutoff_);
00508 
00509         delete (iwt);
00510       }
00511       else
00512       {
00513 #ifndef OPENMS_HAS_TBB
00514         std::cerr << "Error: You requested multi-threaded computation via threading building blocks, but OpenMS has not been configured for TBB usage." << std::endl;
00515         std::cerr << "Error: You need to rebuild OpenMS with -DENABLE_TBB=ON." << std::endl;
00516 #endif
00517       }
00518     }
00519 
00520     static const String getProductName()
00521     {
00522       return "isotope_wavelet";
00523     }
00524 
00525     static FeatureFinderAlgorithm<PeakType, FeatureType> * create()
00526     {
00527       return new FeatureFinderAlgorithmIsotopeWavelet();
00528     }
00529 
00530 protected:
00531 
00533     struct BoxElement
00534     {
00535       DoubleReal mz;
00536       UInt c;       
00537       DoubleReal score;
00538       DoubleReal intens;
00539       DoubleReal RT;       
00540     };
00541 
00542     typedef std::map<UInt, BoxElement> Box;     
00543 
00544     UInt max_charge_;     
00545     DoubleReal intensity_threshold_;     
00546     UInt RT_votes_cutoff_, real_RT_votes_cutoff_, RT_interleave_;     
00547     String use_gpus_, intensity_type_;
00548     bool use_tbb_, use_cuda_, check_PPMs_, hr_data_;
00549     std::vector<UInt> gpu_ids_;     
00550 
00551 #if defined(OPENMS_HAS_TBB) && defined(OPENMS_HAS_CUDA)
00552     tbb::atomic<int> progress_counter_;
00553     Int device_num_, gpu_to_exclude_;
00554 #else
00555     Int progress_counter_;
00556 #endif
00557 
00558     void updateMembers_()
00559     {
00560       max_charge_ = this->param_.getValue("max_charge");
00561       intensity_threshold_ = this->param_.getValue("intensity_threshold");
00562       RT_votes_cutoff_ = this->param_.getValue("sweep_line:rt_votes_cutoff");
00563       RT_interleave_ = this->param_.getValue("sweep_line:rt_interleave");
00564       IsotopeWavelet::setMaxCharge(max_charge_);
00565       check_PPMs_ = ((String)(this->param_.getValue("check_ppm")) == "true");
00566       hr_data_ = ((String)(this->param_.getValue("hr_data")) == "true");
00567       intensity_type_ = ((String)(this->param_.getValue("intensity_type")));
00568 #if defined(OPENMS_HAS_CUDA) || defined(OPENMS_HAS_TBB)
00569       use_gpus_ = this->param_.getValue("parallel:use_gpus");
00570       std::vector<String> tokens;
00571       if (!use_gpus_.split(',', tokens))
00572       {
00573         tokens.push_back(use_gpus_);
00574       }
00575       //Attention: updateMembers_ can be called several times!
00576       gpu_ids_.clear();
00577       if (tokens[0].trim().toInt() == -1)           //no parallelization
00578       {
00579         use_cuda_ = false;
00580         use_tbb_ = false;
00581         return;
00582       }
00583       gpu_ids_.push_back(tokens[0].trim().toInt());
00584       use_cuda_ = true;
00585       use_tbb_ = false;
00586       for (UInt i = 1; i < tokens.size(); ++i)
00587       {
00588         if (tokens[i].trim().toInt() == (Int) gpu_ids_[i - 1])
00589         {
00590           continue;
00591         }
00592         gpu_ids_.push_back(tokens[i].trim().toInt());
00593         use_tbb_ = true;
00594       }
00595 #else
00596       use_cuda_ = false;
00597       use_tbb_ = false;
00598 #endif
00599     }
00600 
00601   };
00602 
00603 } //namespace
00604 
00605 #endif

OpenMS / TOPP release 1.10.0 Documentation generated on Thu Mar 7 2013 09:42:39 using doxygen 1.7.1