20 #include <unordered_map>
25 #include "dirtyAllocator.h"
27 #include "operators.h"
28 #include "marginalTrek++.h"
31 #if ISOSPEC_BUILDING_R
41 unsigned int parse_formula(
const char* formula,
42 std::vector<const double*>& isotope_masses,
43 std::vector<const double*>& isotope_probabilities,
46 unsigned int* confSize);
53 class ISOSPEC_EXPORT_SYMBOL
Iso {
63 void setupMarginals(
const double*
const * _isotopeMasses,
64 const double*
const * _isotopeProbabilities);
89 const int* _isotopeNumbers,
90 const int* _atomCounts,
91 const double*
const * _isotopeMasses,
92 const double*
const * _isotopeProbabilities
96 Iso(
const char* formula);
99 inline Iso(
const std::string& formula) :
Iso(formula.c_str()) {};
109 Iso(
const Iso& other,
bool fullcopy);
115 double getLightestPeakMass()
const;
118 double getHeaviestPeakMass()
const;
125 double getMonoisotopicPeakMass()
const;
131 double getUnlikeliestPeakLProb()
const;
134 double getModeMass()
const;
137 double getTheoreticalAverageMass()
const;
146 void addElement(
int atomCount,
int noIsotopes,
const double* isotopeMasses,
const double* isotopeProbabilities);
166 virtual bool advanceToNextConfiguration() = 0;
168 ISOSPEC_FORCE_INLINE
virtual bool advanceToNextConfigurationWithinLayer() {
return advanceToNextConfiguration(); };
169 ISOSPEC_FORCE_INLINE
virtual bool nextLayer(
double) {
return false; };
175 virtual double lprob()
const {
return partialLProbs[0]; };
181 virtual double mass()
const {
return partialMasses[0]; };
187 virtual double prob()
const {
return partialProbs[0]; };
190 virtual void get_conf_signature(
int* space)
const = 0;
211 std::priority_queue<void*,std::vector<void*>,
ConfOrder> pq;
214 const std::vector<double>** logProbs;
215 const std::vector<double>** masses;
216 const std::vector<int*>** marginalConfs;
223 bool advanceToNextConfiguration()
override final;
232 int* c = getConf(topConf);
237 for(
int ii=0; ii<dimNumber; ii++)
239 memcpy(space, marginalResults[ii]->confs()[c[ii]], isotopeNumbers[ii]*
sizeof(
int));
240 space += isotopeNumbers[ii];
268 double* maxConfsLPSum;
269 const double Lcutoff;
274 const double* lProbs_ptr;
275 const double* lProbs_ptr_start;
276 double* partialLProbs_second;
277 double partialLProbs_second_val, lcfmsv;
283 counter[0] = lProbs_ptr - lProbs_ptr_start;
284 if(marginalOrder !=
nullptr)
285 for(
int ii=0; ii<dimNumber; ii++)
287 int jj = marginalOrder[ii];
288 memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[jj]), isotopeNumbers[ii]*
sizeof(
int));
289 space += isotopeNumbers[ii];
292 for(
int ii=0; ii<dimNumber; ii++)
294 memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[ii]), isotopeNumbers[ii]*
sizeof(
int));
295 space += isotopeNumbers[ii];
309 IsoThresholdGenerator(
Iso&& iso,
double _threshold,
bool _absolute=
true,
int _tabSize=1000,
int _hashSize=1000,
bool reorder_marginals =
true);
319 if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
327 lProbs_ptr = lProbs_ptr_start;
329 int * cntr_ptr = counter;
331 while(idx<dimNumber-1)
339 partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->
get_lProb(counter[idx]);
340 if(partialLProbs[idx] + maxConfsLPSum[idx-1] >= Lcutoff)
342 partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->
get_mass(counter[idx]);
343 partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->
get_prob(counter[idx]);
354 ISOSPEC_FORCE_INLINE
double lprob() const override final {
return partialLProbs_second_val + (*(lProbs_ptr)); };
355 ISOSPEC_FORCE_INLINE
double mass() const override final {
return partialMasses[1] + marginalResults[0]->
get_mass(lProbs_ptr - lProbs_ptr_start); };
356 ISOSPEC_FORCE_INLINE
double prob() const override final {
return partialProbs[1] * marginalResults[0]->
get_prob(lProbs_ptr - lProbs_ptr_start); };
359 void terminate_search();
371 size_t count_confs();
375 ISOSPEC_FORCE_INLINE
void recalc(
int idx)
377 for(; idx > 0; idx--)
379 partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->
get_lProb(counter[idx]);
380 partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->
get_mass(counter[idx]);
381 partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->
get_prob(counter[idx]);
383 partialLProbs_second_val = *partialLProbs_second;
384 partialLProbs[0] = *partialLProbs_second + marginalResults[0]->
get_lProb(counter[0]);
385 lcfmsv = Lcutoff - partialLProbs_second_val;
398 double* maxConfsLPSum;
399 double currentLThreshold, lastLThreshold;
404 const double* lProbs_ptr;
405 const double* lProbs_ptr_start;
406 const double** resetPositions;
407 double* partialLProbs_second;
408 double partialLProbs_second_val, lcfmsv, last_lcfmsv;
413 counter[0] = lProbs_ptr - lProbs_ptr_start;
414 if(marginalOrder !=
nullptr)
415 for(
int ii=0; ii<dimNumber; ii++)
417 int jj = marginalOrder[ii];
418 memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[jj]), isotopeNumbers[ii]*
sizeof(
int));
419 space += isotopeNumbers[ii];
422 for(
int ii=0; ii<dimNumber; ii++)
424 memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[ii]), isotopeNumbers[ii]*
sizeof(
int));
425 space += isotopeNumbers[ii];
430 inline double get_currentLThreshold()
const {
return currentLThreshold; };
432 IsoLayeredGenerator(Iso&& iso,
int _tabSize=1000,
int _hashSize=1000,
bool reorder_marginals =
true,
double t_prob_hint = 0.99);
434 ~IsoLayeredGenerator();
440 if(advanceToNextConfigurationWithinLayer())
442 }
while(nextLayer(-2.0));
446 ISOSPEC_FORCE_INLINE
bool advanceToNextConfigurationWithinLayer() override final
451 if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
458 ISOSPEC_FORCE_INLINE
double lprob() const override final {
return partialLProbs_second_val + (*(lProbs_ptr)); };
459 ISOSPEC_FORCE_INLINE
double mass() const override final {
return partialMasses[1] + marginalResults[0]->
get_mass(lProbs_ptr - lProbs_ptr_start); };
460 ISOSPEC_FORCE_INLINE
double prob() const override final {
return partialProbs[1] * marginalResults[0]->
get_prob(lProbs_ptr - lProbs_ptr_start); };
463 void terminate_search();
467 ISOSPEC_FORCE_INLINE
void recalc(
int idx)
469 for(; idx > 0; idx--)
471 partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->
get_lProb(counter[idx]);
472 partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->
get_mass(counter[idx]);
473 partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->
get_prob(counter[idx]);
475 partialLProbs_second_val = *partialLProbs_second;
476 partialLProbs[0] = partialLProbs_second_val + marginalResults[0]->
get_lProb(counter[0]);
477 lcfmsv = currentLThreshold - partialLProbs_second_val;
478 last_lcfmsv = lastLThreshold - partialLProbs_second_val;
481 virtual bool nextLayer(
double offset)
override final;
491 #if !ISOSPEC_BUILDING_R
493 void printConfigurations(
494 const std::tuple<double*,double*,int*,int>& results,