IsoSpec  2.2.1
fixedEnvelopes.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 <cstdlib>
20 #include <algorithm>
21 #include <vector>
22 #include <utility>
23 
24 #include "isoSpec++.h"
25 
26 #ifdef DEBUG
27 #define ISOSPEC_INIT_TABLE_SIZE 16
28 #else
29 #define ISOSPEC_INIT_TABLE_SIZE 1024
30 #endif
31 
32 namespace IsoSpec
33 {
34 
35 class FixedEnvelope;
36 
37 double AbyssalWassersteinDistanceGrad(FixedEnvelope* const* envelopes, const double* scales, double* ret_gradient, size_t N, double abyss_depth_exp, double abyss_depth_the);
38 
39 
40 class ISOSPEC_EXPORT_SYMBOL FixedEnvelope {
41  protected:
42  double* _masses;
43  double* _probs;
44  int* _confs;
45  size_t _confs_no;
46  int allDim;
47  bool sorted_by_mass;
48  bool sorted_by_prob;
49  double total_prob;
50  size_t current_size;
51  double* tmasses;
52  double* tprobs;
53  int* tconfs;
54  int allDimSizeofInt;
55 
56  public:
57  ISOSPEC_FORCE_INLINE FixedEnvelope() : _masses(nullptr),
58  _probs(nullptr),
59  _confs(nullptr),
60  _confs_no(0),
61  allDim(0),
62  sorted_by_mass(false),
63  sorted_by_prob(false),
64  total_prob(0.0),
65  current_size(0),
66  allDimSizeofInt(0)
67  // Deliberately not initializing tmasses, tprobs, tconfs
68  {};
69 
70  FixedEnvelope(const FixedEnvelope& other);
72 
73  FixedEnvelope(double* masses, double* probs, size_t confs_no, bool masses_sorted = false, bool probs_sorted = false, double _total_prob = NAN);
74 
75  virtual ~FixedEnvelope()
76  {
77  free(_masses);
78  free(_probs);
79  free(_confs);
80  };
81 
82  FixedEnvelope operator+(const FixedEnvelope& other) const;
83  FixedEnvelope operator*(const FixedEnvelope& other) const;
84 
85  inline size_t confs_no() const { return _confs_no; }
86  inline int getAllDim() const { return allDim; }
87 
88  inline const double* masses() const { return _masses; }
89  inline const double* probs() const { return _probs; }
90  inline const int* confs() const { return _confs; }
91 
92  inline double* release_masses() { double* ret = _masses; _masses = nullptr; return ret; }
93  inline double* release_probs() { double* ret = _probs; _probs = nullptr; return ret; }
94  inline int* release_confs() { int* ret = _confs; _confs = nullptr; return ret; }
95  inline void release_everything() { _confs = nullptr; _probs = _masses = nullptr; }
96 
97 
98  inline double mass(size_t i) const { return _masses[i]; }
99  inline double prob(size_t i) const { return _probs[i]; }
100  inline const int* conf(size_t i) const { return _confs + i*allDim; }
101 
102  void sort_by_mass();
103  void sort_by_prob();
104 
105  double get_total_prob();
106  void scale(double factor);
107  void normalize();
108  void shift_mass(double shift);
109  void resample(size_t ionic_current, double beta_bias = 1.0);
110 
111  double empiric_average_mass();
112  double empiric_variance();
113  double empiric_stddev() { return sqrt(empiric_variance()); }
114 
115  double WassersteinDistance(FixedEnvelope& other);
116  double OrientedWassersteinDistance(FixedEnvelope& other);
117  double AbyssalWassersteinDistance(FixedEnvelope& other, double abyss_depth, double other_scale = 1.0);
118  std::tuple<double, double, double> WassersteinMatch(FixedEnvelope& other, double flow_distance, double other_scale = 1.0);
119 
120 
121  static FixedEnvelope LinearCombination(const std::vector<const FixedEnvelope*>& spectra, const std::vector<double>& intensities);
122  static FixedEnvelope LinearCombination(const FixedEnvelope* const * spectra, const double* intensities, size_t size);
123 
124 
125  FixedEnvelope bin(double bin_width = 1.0, double middle = 0.0);
126 
127  private:
128  void sort_by(double* order);
129 
130 
131  protected:
132  template<typename T, bool tgetConfs> ISOSPEC_FORCE_INLINE void store_conf(const T& generator)
133  {
134  *tmasses = generator.mass(); tmasses++;
135  *tprobs = generator.prob(); tprobs++;
136  constexpr_if(tgetConfs) { generator.get_conf_signature(tconfs); tconfs += allDim; }
137  }
138 
139  ISOSPEC_FORCE_INLINE void store_conf(double _mass, double _prob)
140  {
141  if(_confs_no == current_size)
142  reallocate_memory<false>(current_size*2);
143 
144  *tprobs = _prob;
145  *tmasses = _mass;
146  tprobs++;
147  tmasses++;
148 
149  _confs_no++;
150  }
151 
152  template<bool tgetConfs> ISOSPEC_FORCE_INLINE void swap(size_t idx1, size_t idx2, ISOSPEC_MAYBE_UNUSED int* conf_swapspace)
153  {
154  std::swap<double>(this->_probs[idx1], this->_probs[idx2]);
155  std::swap<double>(this->_masses[idx1], this->_masses[idx2]);
156  constexpr_if(tgetConfs)
157  {
158  int* c1 = this->_confs + (idx1*this->allDim);
159  int* c2 = this->_confs + (idx2*this->allDim);
160  memcpy(conf_swapspace, c1, this->allDimSizeofInt);
161  memcpy(c1, c2, this->allDimSizeofInt);
162  memcpy(c2, conf_swapspace, this->allDimSizeofInt);
163  }
164  }
165 
166  template<bool tgetConfs> void reallocate_memory(size_t new_size);
167  void slow_reallocate_memory(size_t new_size);
168 
169  public:
170  template<bool tgetConfs> void threshold_init(Iso&& iso, double threshold, bool absolute);
171 
172  template<bool tgetConfs, typename GenType = IsoLayeredGenerator> void addConfILG(const GenType& generator)
173  {
174  if(this->_confs_no == this->current_size)
175  this->template reallocate_memory<tgetConfs>(this->current_size*2);
176 
177  this->template store_conf<GenType, tgetConfs>(generator);
178  this->_confs_no++;
179  }
180 
181  template<bool tgetConfs> void total_prob_init(Iso&& iso, double target_prob, bool trim);
182 
183  static FixedEnvelope FromThreshold(Iso&& iso, double threshold, bool absolute, bool tgetConfs = false)
184  {
185  FixedEnvelope ret;
186 
187  if(tgetConfs)
188  ret.threshold_init<true>(std::move(iso), threshold, absolute);
189  else
190  ret.threshold_init<false>(std::move(iso), threshold, absolute);
191  return ret;
192  }
193 
194  inline static FixedEnvelope FromThreshold(const Iso& iso, double _threshold, bool _absolute, bool tgetConfs = false)
195  {
196  return FromThreshold(Iso(iso, false), _threshold, _absolute, tgetConfs);
197  }
198 
199  static FixedEnvelope FromTotalProb(Iso&& iso, double target_total_prob, bool optimize, bool tgetConfs = false)
200  {
201  FixedEnvelope ret;
202 
203  if(tgetConfs)
204  ret.total_prob_init<true>(std::move(iso), target_total_prob, optimize);
205  else
206  ret.total_prob_init<false>(std::move(iso), target_total_prob, optimize);
207 
208  return ret;
209  }
210 
211  inline static FixedEnvelope FromTotalProb(const Iso& iso, double _target_total_prob, bool _optimize, bool tgetConfs = false)
212  {
213  return FromTotalProb(Iso(iso, false), _target_total_prob, _optimize, tgetConfs);
214  }
215 
216  template<bool tgetConfs> void stochastic_init(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias);
217 
218  inline static FixedEnvelope FromStochastic(Iso&& iso, size_t _no_molecules, double _precision = 0.9999, double _beta_bias = 5.0, bool tgetConfs = false)
219  {
220  FixedEnvelope ret;
221 
222  if(tgetConfs)
223  ret.stochastic_init<true>(std::move(iso), _no_molecules, _precision, _beta_bias);
224  else
225  ret.stochastic_init<false>(std::move(iso), _no_molecules, _precision, _beta_bias);
226 
227  return ret;
228  }
229 
230  static FixedEnvelope FromStochastic(const Iso& iso, size_t _no_molecules, double _precision = 0.9999, double _beta_bias = 5.0, bool tgetConfs = false)
231  {
232  return FromStochastic(Iso(iso, false), _no_molecules, _precision, _beta_bias, tgetConfs);
233  }
234 
235  static FixedEnvelope Binned(Iso&& iso, double target_total_prob, double bin_width, double bin_middle = 0.0);
236  static FixedEnvelope Binned(const Iso& iso, double target_total_prob, double bin_width, double bin_middle = 0.0)
237  {
238  return Binned(Iso(iso, false), target_total_prob, bin_width, bin_middle);
239  }
240 
241  friend double AbyssalWassersteinDistanceGrad(FixedEnvelope* const* envelopes, const double* scales, double* ret_gradient, size_t N, double abyss_depth_exp, double abyss_depth_the);
242 };
243 
244 } // namespace IsoSpec
The Iso class for the calculation of the isotopic distribution.
Definition: isoSpec++.h:49