OpenMS
SignalToNoiseEstimatorMeanIterative.h
Go to the documentation of this file.
1 // Copyright (c) 2002-present, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
2 // SPDX-License-Identifier: BSD-3-Clause
3 //
4 // --------------------------------------------------------------------------
5 // $Maintainer: Chris Bielow $
6 // $Authors: $
7 // --------------------------------------------------------------------------
8 //
9 
10 #pragma once
11 
15 #include <vector>
16 #include <algorithm> //for std::max_element
17 
18 namespace OpenMS
19 {
43  template <typename Container = MSSpectrum>
45  public SignalToNoiseEstimator<Container>
46  {
47 
48 public:
49 
52 
56 
59 
61 
62 
65  {
66  //set the name for DefaultParamHandler error messages
67  this->setName("SignalToNoiseEstimatorMeanIterative");
68 
69  defaults_.setValue("max_intensity", -1, "maximal intensity considered for histogram construction. By default, it will be calculated automatically (see auto_mode)." \
70  " Only provide this parameter if you know what you are doing (and change 'auto_mode' to '-1')!" \
71  " All intensities EQUAL/ABOVE 'max_intensity' will not be added to the histogram." \
72  " If you choose 'max_intensity' too small, the noise estimate might be too small as well." \
73  " If chosen too big, the bins become quite large (which you could counter by increasing 'bin_count', which increases runtime).", {"advanced"});
74  defaults_.setMinInt("max_intensity", -1);
75 
76  defaults_.setValue("auto_max_stdev_factor", 3.0, "parameter for 'max_intensity' estimation (if 'auto_mode' == 0): mean + 'auto_max_stdev_factor' * stdev", {"advanced"});
77  defaults_.setMinFloat("auto_max_stdev_factor", 0.0);
78  defaults_.setMaxFloat("auto_max_stdev_factor", 999.0);
79 
80 
81  defaults_.setValue("auto_max_percentile", 95, "parameter for 'max_intensity' estimation (if 'auto_mode' == 1): auto_max_percentile th percentile", {"advanced"});
82  defaults_.setMinInt("auto_max_percentile", 0);
83  defaults_.setMaxInt("auto_max_percentile", 100);
84 
85  defaults_.setValue("auto_mode", 0, "method to use to determine maximal intensity: -1 --> use 'max_intensity'; 0 --> 'auto_max_stdev_factor' method (default); 1 --> 'auto_max_percentile' method", {"advanced"});
86  defaults_.setMinInt("auto_mode", -1);
87  defaults_.setMaxInt("auto_mode", 1);
88 
89  defaults_.setValue("win_len", 200.0, "window length in Thomson");
90  defaults_.setMinFloat("win_len", 1.0);
91 
92  defaults_.setValue("bin_count", 30, "number of bins for intensity values");
93  defaults_.setMinInt("bin_count", 3);
94 
95  defaults_.setValue("stdev_mp", 3.0, "multiplier for stdev", {"advanced"});
96  defaults_.setMinFloat("stdev_mp", 0.01);
97  defaults_.setMaxFloat("stdev_mp", 999.0);
98 
99  defaults_.setValue("min_required_elements", 10, "minimum number of elements required in a window (otherwise it is considered sparse)");
100  defaults_.setMinInt("min_required_elements", 1);
101 
102  defaults_.setValue("noise_for_empty_window", std::pow(10.0, 20), "noise value used for sparse windows", {"advanced"});
103 
105  }
106 
109  SignalToNoiseEstimator<Container>(source)
110  {
111  updateMembers_();
112  }
113 
119  {
120  if (&source == this) return *this;
121 
123  updateMembers_();
124  return *this;
125  }
126 
128 
129 
132  {}
133 
134 
135 protected:
136 
137 
142  void computeSTN_(const Container& c) override
143  {
144  //first element in the scan
145  PeakIterator scan_first_ = c.begin();
146  //last element in the scan
147  PeakIterator scan_last_ = c.end();
148 
149  // reset counter for sparse windows
150  double sparse_window_percent = 0;
151 
152  // reset the results
153  stn_estimates_.clear();
154  stn_estimates_.resize(c.size());
155 
156  // maximal range of histogram needs to be calculated first
157  if (auto_mode_ == AUTOMAXBYSTDEV)
158  {
159  // use MEAN+auto_max_intensity_*STDEV as threshold
160  GaussianEstimate gauss_global = SignalToNoiseEstimator<Container>::estimate_(scan_first_, scan_last_);
161  max_intensity_ = gauss_global.mean + std::sqrt(gauss_global.variance) * auto_max_stdev_Factor_;
162  }
163  else if (auto_mode_ == AUTOMAXBYPERCENT)
164  {
165  // get value at "auto_max_percentile_"th percentile
166  // we use a histogram approach here as well.
167  if ((auto_max_percentile_ < 0) || (auto_max_percentile_ > 100))
168  {
170  throw Exception::InvalidValue(__FILE__,
171  __LINE__,
172  OPENMS_PRETTY_FUNCTION,
173  "auto_mode is on AUTOMAXBYPERCENT! auto_max_percentile is not in [0,100]. Use setAutoMaxPercentile(<value>) to change it!",
174  s);
175  }
176 
177  std::vector<int> histogram_auto(100, 0);
178 
179  // find maximum of current scan
180  auto maxIt = std::max_element(c.begin(), c.end() ,[](const PeakType& a, const PeakType& b){ return a.getIntensity() > b.getIntensity();});
181  typename PeakType::IntensityType maxInt = maxIt->getIntensity();
182 
183  double bin_size = maxInt / 100;
184 
185  // fill histogram
186  for(auto& run : c)
187  {
188  ++histogram_auto[(int) (((run).getIntensity() - 1) / bin_size)];
189  }
190 
191  // add up element counts in histogram until ?th percentile is reached
192  int elements_below_percentile = (int) (auto_max_percentile_ * c.size() / 100);
193  int elements_seen = 0;
194  int i = -1;
195  PeakIterator run = scan_first_;
196 
197  while (run != scan_last_ && elements_seen < elements_below_percentile)
198  {
199  ++i;
200  elements_seen += histogram_auto[i];
201  ++run;
202  }
203 
204  max_intensity_ = (((double)i) + 0.5) * bin_size;
205  }
206  else //if (auto_mode_ == MANUAL)
207  {
208  if (max_intensity_ <= 0)
209  {
211  throw Exception::InvalidValue(__FILE__,
212  __LINE__,
213  OPENMS_PRETTY_FUNCTION,
214  "auto_mode is on MANUAL! max_intensity is <=0. Needs to be positive! Use setMaxIntensity(<value>) or enable auto_mode!",
215  s);
216  }
217  }
218 
219  if (max_intensity_ < 0)
220  {
221  std::cerr << "TODO SignalToNoiseEstimatorMedian: the max_intensity_ value should be positive! " << max_intensity_ << std::endl;
222  return;
223  }
224 
225  PeakIterator window_pos_center = scan_first_;
226  PeakIterator window_pos_borderleft = scan_first_;
227  PeakIterator window_pos_borderright = scan_first_;
228 
229  double window_half_size = win_len_ / 2;
230  double bin_size = std::max(1.0, max_intensity_ / bin_count_); // at least size of 1 for intensity bins
231 
232  std::vector<int> histogram(bin_count_, 0);
233  std::vector<double> bin_value(bin_count_, 0);
234  // calculate average intensity that is represented by a bin
235  for (int bin = 0; bin < bin_count_; bin++)
236  {
237  histogram[bin] = 0;
238  bin_value[bin] = (bin + 0.5) * bin_size;
239  }
240  // index of last valid bin during iteration
241  int hist_rightmost_bin;
242  // bin in which a datapoint would fall
243  int to_bin;
244  // mean & stdev of the histogram
245  double hist_mean;
246  double hist_stdev;
247 
248  // tracks elements in current window, which may vary because of unevenly spaced data
249  int elements_in_window = 0;
250  int window_count = 0;
251 
252  double noise; // noise value of a datapoint
253 
255  SignalToNoiseEstimator<Container>::startProgress(0, c.size(), "noise estimation of data");
256 
257  // MAIN LOOP
258  while (window_pos_center != scan_last_)
259  {
260  // erase all elements from histogram that will leave the window on the LEFT side
261  while ((*window_pos_borderleft).getMZ() < (*window_pos_center).getMZ() - window_half_size)
262  {
263  //std::cout << "S: " << (*window_pos_borderleft).getMZ() << " " << ( (*window_pos_center).getMZ() - window_half_size ) << "\n";
264  to_bin = (int) ((std::max((*window_pos_borderleft).getIntensity(), 0.0f)) / bin_size);
265  if (to_bin < bin_count_)
266  {
267  --histogram[to_bin];
268  --elements_in_window;
269  }
270  ++window_pos_borderleft;
271  }
272 
273  //std::printf("S1: %E %E\n", (*window_pos_borderright).getMZ(), (*window_pos_center).getMZ() + window_half_size);
274 
275 
276  // add all elements to histogram that will enter the window on the RIGHT side
277  while ((window_pos_borderright != scan_last_)
278  && ((*window_pos_borderright).getMZ() < (*window_pos_center).getMZ() + window_half_size))
279  {
280  //std::printf("Sb: %E %E %E\n", (*window_pos_borderright).getMZ(), (*window_pos_center).getMZ() + window_half_size, (*window_pos_borderright).getMZ() - ((*window_pos_center).getMZ() + window_half_size));
281 
282  to_bin = (int) ((std::max((*window_pos_borderright).getIntensity(), 0.0f)) / bin_size);
283  if (to_bin < bin_count_)
284  {
285  ++histogram[to_bin];
286  ++elements_in_window;
287  }
288  ++window_pos_borderright;
289  }
290 
291  if (elements_in_window < min_required_elements_)
292  {
293  noise = noise_for_empty_window_;
294  ++sparse_window_percent;
295  }
296  else
297  {
298 
299  hist_rightmost_bin = bin_count_;
300 
301  // do iteration on histogram and find threshold
302  for (int i = 0; i < 3; ++i)
303  {
304  // mean
305  hist_mean = 0;
306  for (int bin = 0; bin < hist_rightmost_bin; ++bin)
307  {
308  //std::cout << "V: " << bin << " " << hist_mean << " " << histogram[bin] << " " << elements_in_window << " " << bin_value[bin] << "\n";
309  // immediate division is numerically more stable
310  hist_mean += histogram[bin] / (double) elements_in_window * bin_value[bin];
311  }
312  //hist_mean = hist_mean / elements_in_window;
313 
314  // stdev
315  hist_stdev = 0;
316  for (int bin = 0; bin < hist_rightmost_bin; ++bin)
317  {
318  double tmp(bin_value[bin] - hist_mean);
319  hist_stdev += histogram[bin] / (double) elements_in_window * tmp * tmp;
320  }
321  hist_stdev = std::sqrt(hist_stdev);
322 
323  //determine new threshold (i.e. the rightmost bin we consider)
324  int estimate = (int) ((hist_mean + hist_stdev * stdev_ - 1) / bin_size + 1);
325  //std::cout << "E: " << hist_mean << " " << hist_stdev << " " << stdev_ << " " << bin_size<< " " << estimate << "\n";
326  hist_rightmost_bin = std::min(estimate, bin_count_);
327  }
328 
329  // just avoid division by 0
330  noise = std::max(1.0, hist_mean);
331  }
332 
333  // store result
334  stn_estimates_[window_count] = (*window_pos_center).getIntensity() / noise;
335 
336 
337 
338  // advance the window center by one datapoint
339  ++window_pos_center;
340  ++window_count;
341  // update progress
343 
344  } // end while
345 
347 
348  sparse_window_percent = sparse_window_percent * 100 / window_count;
349  // warn if percentage of sparse windows is above 20%
350  if (sparse_window_percent > 20)
351  {
352  std::cerr << "WARNING in SignalToNoiseEstimatorMeanIterative: "
353  << sparse_window_percent
354  << "% of all windows were sparse. You should consider increasing 'win_len' or increasing 'min_required_elements'"
355  << " You should also check the MaximalIntensity value (or the parameters for its heuristic estimation)"
356  << " If it is too low, then too many high intensity peaks will be discarded, which leads to a sparse window!"
357  << std::endl;
358  }
359 
360  return;
361 
362  } // end of shiftWindow_
363 
365  void updateMembers_() override
366  {
367  max_intensity_ = (double)param_.getValue("max_intensity");
368  auto_max_stdev_Factor_ = (double)param_.getValue("auto_max_stdev_factor");
369  auto_max_percentile_ = param_.getValue("auto_max_percentile");
370  auto_mode_ = param_.getValue("auto_mode");
371  win_len_ = (double)param_.getValue("win_len");
372  bin_count_ = param_.getValue("bin_count");
373  stdev_ = (double)param_.getValue("stdev_mp");
374  min_required_elements_ = param_.getValue("min_required_elements");
375  noise_for_empty_window_ = (double)param_.getValue("noise_for_empty_window");
376  stn_estimates_.clear();
377  }
378 
388  double win_len_;
392  double stdev_;
398 
399 
400 
401 
402  };
403 
404 } // namespace OpenMS
405 
void defaultsToParam_()
Updates the parameters after the defaults have been set in the constructor.
Param param_
Container for current parameters.
Definition: DefaultParamHandler.h:139
Param defaults_
Container for default parameters. This member should be filled in the constructor of derived classes!
Definition: DefaultParamHandler.h:146
void setName(const String &name)
Mutable access to the name.
Invalid value exception.
Definition: Exception.h:305
void setMaxFloat(const std::string &key, double max)
Sets the maximum value for the floating point or floating point list parameter key.
void setMaxInt(const std::string &key, int max)
Sets the maximum value for the integer or integer list parameter key.
const ParamValue & getValue(const std::string &key) const
Returns a value of a parameter.
void setMinInt(const std::string &key, int min)
Sets the minimum value for the integer or integer list parameter key.
void setValue(const std::string &key, const ParamValue &value, const std::string &description="", const std::vector< std::string > &tags=std::vector< std::string >())
Sets a value.
void setMinFloat(const std::string &key, double min)
Sets the minimum value for the floating point or floating point list parameter key.
float IntensityType
Intensity type.
Definition: Peak2D.h:36
void setProgress(SignedSize value) const
Sets the current progress.
void startProgress(SignedSize begin, SignedSize end, const String &label) const
Initializes the progress display.
void endProgress(UInt64 bytes_processed=0) const
Estimates the signal/noise (S/N) ratio of each data point in a scan based on an iterative scheme whic...
Definition: SignalToNoiseEstimatorMeanIterative.h:46
SignalToNoiseEstimator< Container >::PeakIterator PeakIterator
Definition: SignalToNoiseEstimatorMeanIterative.h:57
SignalToNoiseEstimatorMeanIterative()
default constructor
Definition: SignalToNoiseEstimatorMeanIterative.h:64
double win_len_
range of data points which belong to a window in Thomson
Definition: SignalToNoiseEstimatorMeanIterative.h:388
double stdev_
multiplier for the stdev of intensities
Definition: SignalToNoiseEstimatorMeanIterative.h:392
double noise_for_empty_window_
Definition: SignalToNoiseEstimatorMeanIterative.h:397
~SignalToNoiseEstimatorMeanIterative() override
Destructor.
Definition: SignalToNoiseEstimatorMeanIterative.h:131
SignalToNoiseEstimatorMeanIterative(const SignalToNoiseEstimatorMeanIterative &source)
Copy Constructor.
Definition: SignalToNoiseEstimatorMeanIterative.h:108
double max_intensity_
maximal intensity considered during binning (values above get discarded)
Definition: SignalToNoiseEstimatorMeanIterative.h:380
double auto_max_percentile_
parameter for initial automatic estimation of "max_intensity_" percentile or a stdev
Definition: SignalToNoiseEstimatorMeanIterative.h:384
void computeSTN_(const Container &c) override
Definition: SignalToNoiseEstimatorMeanIterative.h:142
void updateMembers_() override
overridden function from DefaultParamHandler to keep members up to date, when a parameter is changed
Definition: SignalToNoiseEstimatorMeanIterative.h:365
int min_required_elements_
minimal number of elements a window needs to cover to be used
Definition: SignalToNoiseEstimatorMeanIterative.h:394
SignalToNoiseEstimator< Container >::PeakType PeakType
Definition: SignalToNoiseEstimatorMeanIterative.h:58
SignalToNoiseEstimator< Container >::GaussianEstimate GaussianEstimate
Definition: SignalToNoiseEstimatorMeanIterative.h:60
int auto_mode_
determines which method shall be used for estimating "max_intensity_". valid are MANUAL=-1,...
Definition: SignalToNoiseEstimatorMeanIterative.h:386
IntensityThresholdCalculation
method to use for estimating the maximal intensity that is used for histogram calculation
Definition: SignalToNoiseEstimatorMeanIterative.h:51
@ MANUAL
Definition: SignalToNoiseEstimatorMeanIterative.h:51
@ AUTOMAXBYSTDEV
Definition: SignalToNoiseEstimatorMeanIterative.h:51
@ AUTOMAXBYPERCENT
Definition: SignalToNoiseEstimatorMeanIterative.h:51
SignalToNoiseEstimatorMeanIterative & operator=(const SignalToNoiseEstimatorMeanIterative &source)
Definition: SignalToNoiseEstimatorMeanIterative.h:118
int bin_count_
number of bins in the histogram
Definition: SignalToNoiseEstimatorMeanIterative.h:390
double auto_max_stdev_Factor_
parameter for initial automatic estimation of "max_intensity_": a stdev multiplier
Definition: SignalToNoiseEstimatorMeanIterative.h:382
This class represents the abstract base class of a signal to noise estimator.
Definition: SignalToNoiseEstimator.h:33
double variance
variance of estimated Gaussian
Definition: SignalToNoiseEstimator.h:108
SignalToNoiseEstimator & operator=(const SignalToNoiseEstimator &source)
Assignment operator.
Definition: SignalToNoiseEstimator.h:60
PeakIterator::value_type PeakType
Definition: SignalToNoiseEstimator.h:40
GaussianEstimate estimate_(const PeakIterator &scan_first_, const PeakIterator &scan_last_) const
calculate mean & stdev of intensities of a spectrum
Definition: SignalToNoiseEstimator.h:113
double mean
mean of estimated Gaussian
Definition: SignalToNoiseEstimator.h:107
std::vector< double > stn_estimates_
stores the noise estimate for each peak
Definition: SignalToNoiseEstimator.h:146
Container::const_iterator PeakIterator
Definition: SignalToNoiseEstimator.h:39
protected struct to store parameters my, sigma for a Gaussian distribution
Definition: SignalToNoiseEstimator.h:106
A more convenient string class.
Definition: String.h:34
const double c
Definition: Constants.h:188
Main OpenMS namespace.
Definition: openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h:19