IsoSpec  2.1.2
marginalTrek++.h
1 /*
2  * Copyright (C) 2015-2020 Mateusz Łącki and Michał Startek.
3  *
4  * This file is part of IsoSpec.
5  *
6  * IsoSpec is free software: you can redistribute it and/or modify
7  * it under the terms of the Simplified ("2-clause") BSD licence.
8  *
9  * IsoSpec is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12  *
13  * You should have received a copy of the Simplified BSD Licence
14  * along with IsoSpec. If not, see <https://opensource.org/licenses/BSD-2-Clause>.
15  */
16 
17 #pragma once
18 
19 #include <queue>
20 #include <algorithm>
21 #include <vector>
22 #include <functional>
23 #include <utility>
24 #include "conf.h"
25 #include "allocator.h"
26 #include "operators.h"
27 #include "summator.h"
28 #include "pod_vector.h"
29 
30 
31 namespace IsoSpec
32 {
33 
35 
42 class Marginal
43 {
44  private:
45  bool disowned;
46  protected:
47  const unsigned int isotopeNo;
48  const unsigned int atomCnt;
49  const double* const atom_lProbs;
50  const double* const atom_masses;
51  const double loggamma_nominator;
52  Conf mode_conf;
53  double mode_lprob;
56  public:
58 
65  Marginal(
66  const double* _masses, // masses size = logProbs size = isotopeNo
67  const double* _probs,
68  int _isotopeNo,
69  int _atomCnt
70  );
71 
72  // Get rid of the C++ generated assignment constructor.
73  Marginal& operator= (const Marginal& other) = delete;
74 
76  Marginal(const Marginal& other);
77 
79  Marginal(Marginal&& other);
80 
82  virtual ~Marginal();
83 
85 
88  inline int get_isotopeNo() const { return isotopeNo; }
89 
90  inline const double* get_lProbs() const { return atom_lProbs; }
91 
93 
96  double getLightestConfMass() const;
97 
99 
102  double getHeaviestConfMass() const;
103 
105 
109  double getMonoisotopicConfMass() const;
110 
112 
115  inline double getModeMass() { ensureModeConf(); return calc_mass(mode_conf, atom_masses, isotopeNo); }
116 
118 
121  inline double getModeLProb() { ensureModeConf(); return mode_lprob; }
122 
124  inline double fastGetModeLProb() { return mode_lprob; }
125 
127 
130 // inline double getModeProb() const { return exp(getModeLProb()); }
131 
133  Conf computeModeConf() const;
134 
136 
139  inline double getSmallestLProb() const { return atomCnt * *std::min_element(atom_lProbs, atom_lProbs+isotopeNo); }
140 
142 
145  double getAtomAverageMass() const;
146 
148 
151  inline double getTheoreticalAverageMass() const { return getAtomAverageMass() * atomCnt; }
152 
154 
158  protected:
159  ISOSPEC_FORCE_INLINE double unnormalized_logProb(Conf conf) const { double ret = 0.0; for(size_t ii = 0; ii < isotopeNo; ii++) ret += minuslogFactorial(conf[ii]) + conf[ii] * atom_lProbs[ii]; return ret; }
160  ISOSPEC_FORCE_INLINE double logProb(Conf conf) const { return loggamma_nominator + unnormalized_logProb(conf); }
161  public:
163  double variance() const;
164 
166  double getLogSizeEstimate(double logEllipsoidRadius) const;
167 
168  inline void ensureModeConf() { if (mode_conf == nullptr) setupMode(); }
169  private:
170  void setupMode();
171 };
172 
173 
175 class MarginalTrek : public Marginal
176 {
177  private:
178  int current_count;
179  const ConfOrderMarginal orderMarginal;
180  std::priority_queue<ProbAndConfPtr, pod_vector<ProbAndConfPtr> > pq;
181  Allocator<int> allocator;
182  pod_vector<double> _conf_lprobs;
183  pod_vector<double> _conf_masses;
184  pod_vector<int*> _confs;
185 
187  bool add_next_conf();
188 
189  public:
191 
195  MarginalTrek(
196  Marginal&& m,
197  int tabSize = 1000,
198  int hashSize = 1000
199  ); // NOLINT(runtime/explicit) - Constructor deliberately left usable as a conversion.
200 
201  MarginalTrek(const MarginalTrek& other) = delete;
202  MarginalTrek& operator=(const MarginalTrek& other) = delete;
203 
205 
211  inline bool probeConfigurationIdx(int idx)
212  {
213  while(current_count <= idx)
214  if(!add_next_conf())
215  return false;
216  return true;
217  }
218 
220 
223  inline double getModeLProb() const { return mode_lprob; }
224 
225 
226  inline const pod_vector<double>& conf_lprobs() const { return _conf_lprobs; }
227  inline const pod_vector<double>& conf_masses() const { return _conf_masses; }
228  inline const pod_vector<Conf>& confs() const { return _confs; }
229 
230 
231  virtual ~MarginalTrek();
232 };
233 
234 
236 
245 {
246  protected:
247  pod_vector<Conf> configurations;
248  Conf* confs;
249  unsigned int no_confs;
250  double* masses;
251  pod_vector<double> lProbs;
252  double* probs;
253  Allocator<int> allocator;
254  public:
256 
264  Marginal&& m,
265  double lCutOff,
266  bool sort = true,
267  int tabSize = 1000,
268  int hashSize = 1000
269  );
270 
271  PrecalculatedMarginal(const PrecalculatedMarginal& other) = delete;
272  PrecalculatedMarginal& operator=(const PrecalculatedMarginal& other) = delete;
273 
275  virtual ~PrecalculatedMarginal();
276 
278 
281  inline bool inRange(unsigned int idx) const { return idx < no_confs; }
282 
284 
288  inline const double& get_lProb(int idx) const { return lProbs[idx]; }
289 
291 
295  inline const double& get_prob(int idx) const { return probs[idx]; }
296 
298 
302  inline const double& get_mass(int idx) const { return masses[idx]; }
303 
305 
308  inline const double* get_lProbs_ptr() const { return lProbs.data(); }
309 
311 
314  inline const double* get_masses_ptr() const { return masses; }
315 
316 
318 
322  inline const Conf& get_conf(int idx) const { return confs[idx]; }
323 
325 
328  inline unsigned int get_no_confs() const { return no_confs; }
329 
331 
334  inline double getModeLProb() const { return mode_lprob; }
335 };
336 
337 
338 
340 
343 class LayeredMarginal : public Marginal
344 {
345  private:
346  double current_threshold;
347  pod_vector<Conf> configurations;
348  pod_vector<Conf> fringe;
349  pod_vector<double> fringe_unn_lprobs;
350  Allocator<int> allocator;
351  const ConfEqual equalizer;
352  const KeyHasher keyHasher;
353  pod_vector<double> lProbs;
354  pod_vector<double> probs;
355  pod_vector<double> masses;
356  double* guarded_lProbs;
357 
358  public:
360 
364  LayeredMarginal(Marginal&& m, int tabSize = 1000, int hashSize = 1000); // NOLINT(runtime/explicit) - constructor deliberately left usable as a conversion
365 
366  LayeredMarginal(const LayeredMarginal& other) = delete;
367  LayeredMarginal& operator=(const LayeredMarginal& other) = delete;
368 
370 
374  bool extend(double new_threshold, bool do_sort = true);
375 
377  inline double get_lProb(int idx) const { return guarded_lProbs[idx]; } // access to idx == -1 is valid and gives a guardian of +inf
378 
380  inline double get_prob(int idx) const { return probs[idx]; }
381 
383  inline double get_mass(int idx) const { return masses[idx]; }
384 
386  inline const double* get_lProbs_ptr() const { return lProbs.data()+1; }
387 
389  inline const Conf& get_conf(int idx) const { return configurations[idx]; }
390 
392  inline unsigned int get_no_confs() const { return configurations.size(); }
393 
395  double get_min_mass() const;
396 
398  double get_max_mass() const;
399 
401 
404  inline double getModeLProb() const { return mode_lprob; }
405 };
406 
407 
408 
409 } // namespace IsoSpec
IsoSpec::Marginal::getMonoisotopicConfMass
double getMonoisotopicConfMass() const
Get the mass of the monoisotopic subisotopologue.
Definition: marginalTrek++.cpp:268
IsoSpec::PrecalculatedMarginal::PrecalculatedMarginal
PrecalculatedMarginal(Marginal &&m, double lCutOff, bool sort=true, int tabSize=1000, int hashSize=1000)
The move constructor (disowns the Marginal).
Definition: marginalTrek++.cpp:405
IsoSpec::Marginal::variance
double variance() const
Calculate the variance of the theoretical distribution describing the subisotopologue.
Definition: marginalTrek++.cpp:289
IsoSpec::Marginal::Marginal
Marginal(const double *_masses, const double *_probs, int _isotopeNo, int _atomCnt)
Class constructor.
Definition: marginalTrek++.cpp:165
pod_vector< double >
IsoSpec::PrecalculatedMarginal::get_masses_ptr
const double * get_masses_ptr() const
Get the table of the masses of subisotopologues.
Definition: marginalTrek++.h:314
IsoSpec::PrecalculatedMarginal
Precalculated Marginal class.
Definition: marginalTrek++.h:244
IsoSpec::Marginal
The marginal distribution class (a subisotopologue).
Definition: marginalTrek++.h:42
IsoSpec::LayeredMarginal::get_no_confs
unsigned int get_no_confs() const
Get the number of precomputed subisotopologues, see details in PrecalculatedMarginal::get_no_confs.
Definition: marginalTrek++.h:392
IsoSpec::Marginal::unnormalized_logProb
ISOSPEC_FORCE_INLINE double unnormalized_logProb(Conf conf) const
Calculate the log-probability of a given subisotopologue.
Definition: marginalTrek++.h:159
IsoSpec::PrecalculatedMarginal::get_lProbs_ptr
const double * get_lProbs_ptr() const
Get the table of the log-probabilities of subisotopologues.
Definition: marginalTrek++.h:308
IsoSpec::PrecalculatedMarginal::get_mass
const double & get_mass(int idx) const
Get the mass of the idx-th subisotopologue.
Definition: marginalTrek++.h:302
IsoSpec::LayeredMarginal::LayeredMarginal
LayeredMarginal(Marginal &&m, int tabSize=1000, int hashSize=1000)
Move constructor: specializes the Marginal class.
Definition: marginalTrek++.cpp:518
IsoSpec::Allocator< int >
IsoSpec::Marginal::fastGetModeLProb
double fastGetModeLProb()
Get the log-probability of the mode subisotopologue. Results undefined if ensureModeConf() wasn't cal...
Definition: marginalTrek++.h:124
IsoSpec::Marginal::getHeaviestConfMass
double getHeaviestConfMass() const
Get the mass of the heaviest subisotopologue.
Definition: marginalTrek++.cpp:259
IsoSpec::Marginal::getSmallestLProb
double getSmallestLProb() const
The the log-probability of the lightest subisotopologue.
Definition: marginalTrek++.h:139
IsoSpec::MarginalTrek::getModeLProb
double getModeLProb() const
Get the log-probability of the mode subisotopologue.
Definition: marginalTrek++.h:223
IsoSpec
Definition: allocator.cpp:20
IsoSpec::PrecalculatedMarginal::get_conf
const Conf & get_conf(int idx) const
Get the counts of isotopes that define the subisotopologue.
Definition: marginalTrek++.h:322
IsoSpec::PrecalculatedMarginal::get_lProb
const double & get_lProb(int idx) const
Get the log-probability of the idx-th subisotopologue.
Definition: marginalTrek++.h:288
IsoSpec::PrecalculatedMarginal::~PrecalculatedMarginal
virtual ~PrecalculatedMarginal()
Destructor.
Definition: marginalTrek++.cpp:504
IsoSpec::LayeredMarginal::get_mass
double get_mass(int idx) const
get the mass of the idx-th subisotopologue, see details in PrecalculatedMarginal::get_mass.
Definition: marginalTrek++.h:383
IsoSpec::Marginal::getLightestConfMass
double getLightestConfMass() const
Get the mass of the lightest subisotopologue.
Definition: marginalTrek++.cpp:250
IsoSpec::Marginal::mode_lprob
double mode_lprob
Definition: marginalTrek++.h:53
IsoSpec::ConfOrderMarginal
Definition: operators.h:87
IsoSpec::KeyHasher
Definition: operators.h:28
IsoSpec::Marginal::mode_conf
Conf mode_conf
Definition: marginalTrek++.h:52
IsoSpec::Marginal::atom_masses
const double *const atom_masses
Definition: marginalTrek++.h:50
IsoSpec::Marginal::computeModeConf
Conf computeModeConf() const
The the probability of the mode subisotopologue.
Definition: marginalTrek++.cpp:235
IsoSpec::LayeredMarginal::getModeLProb
double getModeLProb() const
Get the log-probability of the mode subisotopologue.
Definition: marginalTrek++.h:404
IsoSpec::Marginal::atomCnt
const unsigned int atomCnt
Definition: marginalTrek++.h:48
IsoSpec::ConfEqual
Definition: operators.h:53
IsoSpec::LayeredMarginal::get_conf
const Conf & get_conf(int idx) const
get the counts of isotopes that define the subisotopologue, see details in PrecalculatedMarginal::get...
Definition: marginalTrek++.h:389
IsoSpec::PrecalculatedMarginal::get_no_confs
unsigned int get_no_confs() const
Get the number of precomputed subisotopologues.
Definition: marginalTrek++.h:328
IsoSpec::Marginal::getAtomAverageMass
double getAtomAverageMass() const
The average mass of a single atom.
Definition: marginalTrek++.cpp:281
IsoSpec::Marginal::getLogSizeEstimate
double getLogSizeEstimate(double logEllipsoidRadius) const
Return estimated logarithm of size of the marginal at a given ellipsoid radius.
Definition: marginalTrek++.cpp:301
IsoSpec::Marginal::getTheoreticalAverageMass
double getTheoreticalAverageMass() const
The theoretical average mass of the molecule.
Definition: marginalTrek++.h:151
IsoSpec::LayeredMarginal::get_min_mass
double get_min_mass() const
Get the minimal mass in current layer.
Definition: marginalTrek++.cpp:641
IsoSpec::LayeredMarginal
LayeredMarginal class.
Definition: marginalTrek++.h:343
IsoSpec::LayeredMarginal::get_lProbs_ptr
const double * get_lProbs_ptr() const
get the pointer to lProbs array. Accessing index -1 is legal and returns a guardian of -inf....
Definition: marginalTrek++.h:386
IsoSpec::MarginalTrek::probeConfigurationIdx
bool probeConfigurationIdx(int idx)
Check if the table of computed subisotopologues does not have to be extended.
Definition: marginalTrek++.h:211
IsoSpec::Marginal::get_isotopeNo
int get_isotopeNo() const
Get the number of isotopes of the investigated element.
Definition: marginalTrek++.h:88
IsoSpec::Marginal::atom_lProbs
const double *const atom_lProbs
Definition: marginalTrek++.h:49
IsoSpec::Marginal::~Marginal
virtual ~Marginal()
Destructor.
Definition: marginalTrek++.cpp:224
IsoSpec::LayeredMarginal::get_prob
double get_prob(int idx) const
get the probability of the idx-th subisotopologue, see details in PrecalculatedMarginal::get_eProb.
Definition: marginalTrek++.h:380
IsoSpec::Marginal::getModeMass
double getModeMass()
The the mass of the mode subisotopologue.
Definition: marginalTrek++.h:115
IsoSpec::LayeredMarginal::extend
bool extend(double new_threshold, bool do_sort=true)
Extend the set of computed subisotopologues to those above the new threshold.
Definition: marginalTrek++.cpp:529
IsoSpec::MarginalTrek::MarginalTrek
MarginalTrek(Marginal &&m, int tabSize=1000, int hashSize=1000)
Move constructor: specializes the Marginal class.
Definition: marginalTrek++.cpp:323
IsoSpec::PrecalculatedMarginal::inRange
bool inRange(unsigned int idx) const
Is there a subisotopologue with a given number?
Definition: marginalTrek++.h:281
IsoSpec::MarginalTrek
The marginal distribution class (a subisotopologue).
Definition: marginalTrek++.h:175
IsoSpec::LayeredMarginal::get_lProb
double get_lProb(int idx) const
get the log-probability of the idx-th subisotopologue, see details in PrecalculatedMarginal::get_lPro...
Definition: marginalTrek++.h:377
IsoSpec::PrecalculatedMarginal::get_prob
const double & get_prob(int idx) const
Get the probability of the idx-th subisotopologue.
Definition: marginalTrek++.h:295
IsoSpec::Marginal::isotopeNo
const unsigned int isotopeNo
Definition: marginalTrek++.h:47
IsoSpec::Marginal::loggamma_nominator
const double loggamma_nominator
Definition: marginalTrek++.h:51
IsoSpec::LayeredMarginal::get_max_mass
double get_max_mass() const
Get the maximal mass in current layer.
Definition: marginalTrek++.cpp:651
IsoSpec::Marginal::getModeLProb
double getModeLProb()
Get the log-probability of the mode subisotopologue.
Definition: marginalTrek++.h:121
IsoSpec::PrecalculatedMarginal::getModeLProb
double getModeLProb() const
Get the log-probability of the mode subisotopologue.
Definition: marginalTrek++.h:334