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
32namespace IsoSpec
33{
34
35class FixedEnvelope;
36
37double AbyssalWassersteinDistanceGrad(FixedEnvelope* const* envelopes, const double* scales, double* ret_gradient, size_t N, double abyss_depth_exp, double abyss_depth_the);
38
39
40class 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