IsoSpec  1.95
isoSpec++.cpp
1 /*
2  * Copyright (C) 2015-2019 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 
18 #include <cmath>
19 #include <algorithm>
20 #include <vector>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <tuple>
24 #include <unordered_map>
25 #include <queue>
26 #include <utility>
27 #include <iostream>
28 #include <iomanip>
29 #include <cctype>
30 #include <stdexcept>
31 #include <string>
32 #include <limits>
33 #include <assert.h>
34 #include <ctype.h>
35 #include "platform.h"
36 #include "conf.h"
37 #include "dirtyAllocator.h"
38 #include "operators.h"
39 #include "summator.h"
40 #include "marginalTrek++.h"
41 #include "isoSpec++.h"
42 #include "misc.h"
43 #include "element_tables.h"
44 
45 
46 using namespace std;
47 
48 namespace IsoSpec
49 {
50 
51 Iso::Iso() :
52 disowned(false),
53 dimNumber(0),
54 isotopeNumbers(new int[0]),
55 atomCounts(new int[0]),
56 confSize(0),
57 allDim(0),
58 marginals(new Marginal*[0]),
59 modeLProb(0.0)
60 {}
61 
62 
63 Iso::Iso(
64  int _dimNumber,
65  const int* _isotopeNumbers,
66  const int* _atomCounts,
67  const double* const * _isotopeMasses,
68  const double* const * _isotopeProbabilities
69 ) :
70 disowned(false),
71 dimNumber(_dimNumber),
72 isotopeNumbers(array_copy<int>(_isotopeNumbers, _dimNumber)),
73 atomCounts(array_copy<int>(_atomCounts, _dimNumber)),
74 confSize(_dimNumber * sizeof(int)),
75 allDim(0),
76 marginals(nullptr),
77 modeLProb(0.0)
78 {
79  try{
80  setupMarginals(_isotopeMasses, _isotopeProbabilities);
81  }
82  catch(...)
83  {
84  delete[] isotopeNumbers;
85  delete[] atomCounts;
86  throw;
87  }
88 }
89 
90 Iso::Iso(Iso&& other) :
91 disowned(other.disowned),
92 dimNumber(other.dimNumber),
93 isotopeNumbers(other.isotopeNumbers),
94 atomCounts(other.atomCounts),
95 confSize(other.confSize),
96 allDim(other.allDim),
97 marginals(other.marginals),
98 modeLProb(other.modeLProb)
99 {
100  other.disowned = true;
101 }
102 
103 
104 Iso::Iso(const Iso& other, bool fullcopy) :
105 disowned(fullcopy ? throw std::logic_error("Not implemented") : true),
106 dimNumber(other.dimNumber),
107 isotopeNumbers(fullcopy ? array_copy<int>(other.isotopeNumbers, dimNumber) : other.isotopeNumbers),
108 atomCounts(fullcopy ? array_copy<int>(other.atomCounts, dimNumber) : other.atomCounts),
109 confSize(other.confSize),
110 allDim(other.allDim),
111 marginals(fullcopy ? throw std::logic_error("Not implemented") : other.marginals),
112 modeLProb(other.modeLProb)
113 {}
114 
115 
116 inline void Iso::setupMarginals(const double* const * _isotopeMasses, const double* const * _isotopeProbabilities)
117 {
118  if (marginals == nullptr)
119  {
120  int ii = 0;
121  try
122  {
123  marginals = new Marginal*[dimNumber];
124  while(ii < dimNumber)
125  {
126  allDim += isotopeNumbers[ii];
127  marginals[ii] = new Marginal(
128  _isotopeMasses[ii],
129  _isotopeProbabilities[ii],
130  isotopeNumbers[ii],
131  atomCounts[ii]
132  );
133  modeLProb += marginals[ii]->getModeLProb();
134  ii++;
135  }
136  }
137  catch(...)
138  {
139  ii--;
140  while(ii >= 0)
141  {
142  delete marginals[ii];
143  ii--;
144  }
145  delete[] marginals;
146  marginals = nullptr;
147  throw;
148  }
149  }
150 
151 }
152 
154 {
155  if(!disowned)
156  {
157  if (marginals != nullptr)
158  dealloc_table(marginals, dimNumber);
159  delete[] isotopeNumbers;
160  delete[] atomCounts;
161  }
162 }
163 
164 
166 {
167  double mass = 0.0;
168  for (int ii=0; ii<dimNumber; ii++)
169  mass += marginals[ii]->getLightestConfMass();
170  return mass;
171 }
172 
174 {
175  double mass = 0.0;
176  for (int ii=0; ii<dimNumber; ii++)
177  mass += marginals[ii]->getHeaviestConfMass();
178  return mass;
179 }
180 
182 {
183  double mass = 0.0;
184  for (int ii=0; ii<dimNumber; ii++)
185  mass += marginals[ii]->getMonoisotopicConfMass();
186  return mass;
187 }
188 
190 {
191  double lprob = 0.0;
192  for (int ii=0; ii<dimNumber; ii++)
193  lprob += marginals[ii]->getSmallestLProb();
194  return lprob;
195 }
196 
197 double Iso::getModeMass() const
198 {
199  double mass = 0.0;
200  for (int ii=0; ii<dimNumber; ii++)
201  mass += marginals[ii]->getModeMass();
202  return mass;
203 }
204 
206 {
207  double mass = 0.0;
208  for (int ii=0; ii<dimNumber; ii++)
209  mass += marginals[ii]->getTheoreticalAverageMass();
210  return mass;
211 }
212 
213 
214 Iso::Iso(const char* formula) :
215 disowned(false),
216 allDim(0),
217 marginals(nullptr),
218 modeLProb(0.0)
219 {
220  std::vector<const double*> isotope_masses;
221  std::vector<const double*> isotope_probabilities;
222 
223  dimNumber = parse_formula(formula, isotope_masses, isotope_probabilities, &isotopeNumbers, &atomCounts, &confSize);
224 
225  setupMarginals(isotope_masses.data(), isotope_probabilities.data());
226 }
227 
228 
229 void Iso::addElement(int atomCount, int noIsotopes, const double* isotopeMasses, const double* isotopeProbabilities)
230 {
231  Marginal* m = new Marginal(isotopeMasses, isotopeProbabilities, noIsotopes, atomCount);
232  modeLProb += m->getModeLProb();
233  realloc_append<int>(&isotopeNumbers, noIsotopes, dimNumber);
234  realloc_append<int>(&atomCounts, atomCount, dimNumber);
235  realloc_append<Marginal*>(&marginals, m, dimNumber);
236  dimNumber++;
237  confSize += sizeof(int);
238  allDim += noIsotopes;
239 
240 }
241 
242 unsigned int parse_formula(const char* formula, std::vector<const double*>& isotope_masses, std::vector<const double*>& isotope_probabilities, int** isotopeNumbers, int** atomCounts, unsigned int* confSize)
243 {
244  // This function is NOT guaranteed to be secure against malicious input. It should be used only for debugging.
245  size_t slen = strlen(formula);
246  // Yes, it would be more elegant to use std::string here, but it's the only promiment place where it would be used in IsoSpec, and avoiding it here
247  // means we can run the whole thing through Clang's memory sanitizer without the need for instrumented libc++/libstdc++. That's worth messing with char pointers a
248  // little bit.
249  std::vector<std::pair<const char*, size_t> > elements;
250  std::vector<int> numbers;
251 
252  if(slen == 0)
253  throw invalid_argument("Invalid formula: can't be empty");
254 
255  if(!isdigit(formula[slen-1]))
256  throw invalid_argument("Invalid formula: every element must be followed by a number - write H2O1 and not H2O for water");
257 
258  for(size_t ii=0; ii<slen; ii++)
259  if(!isdigit(formula[ii]) && !isalpha(formula[ii]))
260  throw invalid_argument("Invalid formula: contains invalid (non-digit, non-alpha) character");
261 
262  size_t position = 0;
263  size_t elem_end = 0;
264  size_t digit_end = 0;
265 
266  while(position < slen)
267  {
268  elem_end = position;
269  while(isalpha(formula[elem_end]))
270  elem_end++;
271  digit_end = elem_end;
272  while(isdigit(formula[digit_end]))
273  digit_end++;
274  elements.emplace_back(&formula[position], elem_end-position);
275  numbers.push_back(atoi(&formula[elem_end]));
276  position = digit_end;
277  }
278 
279  std::vector<int> element_indexes;
280 
281  for (unsigned int i=0; i<elements.size(); i++)
282  {
283  int idx = -1;
284  for(int j=0; j<ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES; j++)
285  {
286  if ((strlen(elem_table_symbol[j]) == elements[i].second) && (strncmp(elements[i].first, elem_table_symbol[j], elements[i].second) == 0))
287  {
288  idx = j;
289  break;
290  }
291  }
292  if(idx < 0)
293  throw invalid_argument("Invalid formula");
294  element_indexes.push_back(idx);
295  }
296 
297  vector<int> _isotope_numbers;
298 
299  for(vector<int>::iterator it = element_indexes.begin(); it != element_indexes.end(); ++it)
300  {
301  int num = 0;
302  int at_idx = *it;
303  int atomicNo = elem_table_atomicNo[at_idx];
304  while(at_idx < ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES && elem_table_atomicNo[at_idx] == atomicNo)
305  {
306  at_idx++;
307  num++;
308  }
309  _isotope_numbers.push_back(num);
310  }
311 
312  for(vector<int>::iterator it = element_indexes.begin(); it != element_indexes.end(); ++it)
313  {
314  isotope_masses.push_back(&elem_table_mass[*it]);
315  isotope_probabilities.push_back(&elem_table_probability[*it]);
316  };
317 
318  const unsigned int dimNumber = elements.size();
319 
320  *isotopeNumbers = array_copy<int>(_isotope_numbers.data(), dimNumber);
321  *atomCounts = array_copy<int>(numbers.data(), dimNumber);
322  *confSize = dimNumber * sizeof(int);
323 
324  return dimNumber;
325 
326 }
327 
328 
329 /*
330  * ----------------------------------------------------------------------------------------------------------
331  */
332 
333 
334 
335 IsoGenerator::IsoGenerator(Iso&& iso, bool alloc_partials) :
336  Iso(std::move(iso)),
337  partialLProbs(alloc_partials ? new double[dimNumber+1] : nullptr),
338  partialMasses(alloc_partials ? new double[dimNumber+1] : nullptr),
339  partialProbs(alloc_partials ? new double[dimNumber+1] : nullptr)
340 {
341  if(alloc_partials)
342  {
343  partialLProbs[dimNumber] = 0.0;
344  partialMasses[dimNumber] = 0.0;
345  partialProbs[dimNumber] = 1.0;
346  }
347 }
348 
349 
351 {
352  if(partialLProbs != nullptr)
353  delete[] partialLProbs;
354  if(partialMasses != nullptr)
355  delete[] partialMasses;
356  if(partialProbs != nullptr)
357  delete[] partialProbs;
358 }
359 
360 
361 
362 /*
363  * ----------------------------------------------------------------------------------------------------------
364  */
365 
366 
367 
368 IsoThresholdGenerator::IsoThresholdGenerator(Iso&& iso, double _threshold, bool _absolute, int tabSize, int hashSize, bool reorder_marginals)
369 : IsoGenerator(std::move(iso)),
370 Lcutoff(_threshold <= 0.0 ? std::numeric_limits<double>::lowest() : (_absolute ? log(_threshold) : log(_threshold) + modeLProb))
371 {
372  counter = new int[dimNumber];
373  maxConfsLPSum = new double[dimNumber-1];
374  marginalResultsUnsorted = new PrecalculatedMarginal*[dimNumber];
375 
376  empty = false;
377 
378  for(int ii=0; ii<dimNumber; ii++)
379  {
380  counter[ii] = 0;
381  marginalResultsUnsorted[ii] = new PrecalculatedMarginal(std::move(*(marginals[ii])),
382  Lcutoff - modeLProb + marginals[ii]->getModeLProb(),
383  true,
384  tabSize,
385  hashSize);
386 
387  if(!marginalResultsUnsorted[ii]->inRange(0))
388  empty = true;
389  }
390 
391  if(reorder_marginals && dimNumber > 1)
392  {
393  OrderMarginalsBySizeDecresing comparator(marginalResultsUnsorted);
394  int* tmpMarginalOrder = new int[dimNumber];
395 
396  for(int ii=0; ii<dimNumber; ii++)
397  tmpMarginalOrder[ii] = ii;
398 
399  std::sort(tmpMarginalOrder, tmpMarginalOrder + dimNumber, comparator);
400  marginalResults = new PrecalculatedMarginal*[dimNumber];
401 
402  for(int ii=0; ii<dimNumber; ii++)
403  marginalResults[ii] = marginalResultsUnsorted[tmpMarginalOrder[ii]];
404 
405  marginalOrder = new int[dimNumber];
406  for(int ii = 0; ii<dimNumber; ii++)
407  marginalOrder[tmpMarginalOrder[ii]] = ii;
408 
409  delete[] tmpMarginalOrder;
410 
411  }
412  else
413  {
414  marginalResults = marginalResultsUnsorted;
415  marginalOrder = nullptr;
416  }
417 
418  lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr();
419 
420  if(dimNumber > 1)
421  maxConfsLPSum[0] = marginalResults[0]->getModeLProb();
422 
423  for(int ii=1; ii<dimNumber-1; ii++)
424  maxConfsLPSum[ii] = maxConfsLPSum[ii-1] + marginalResults[ii]->getModeLProb();
425 
426  lProbs_ptr = lProbs_ptr_start;
427 
428  partialLProbs_second = partialLProbs;
429  partialLProbs_second++;
430 
431  if(!empty)
432  {
433  recalc(dimNumber-1);
434  counter[0]--;
435  lProbs_ptr--;
436  }
437  else
438  {
440  lcfmsv = std::numeric_limits<double>::infinity();
441  }
442 
443 
444 
445 }
446 
448 {
449  for(int ii=0; ii<dimNumber; ii++)
450  {
451  counter[ii] = marginalResults[ii]->get_no_confs()-1;
452  partialLProbs[ii] = -std::numeric_limits<double>::infinity();
453  }
454  partialLProbs[dimNumber] = -std::numeric_limits<double>::infinity();
455  lProbs_ptr = lProbs_ptr_start + marginalResults[0]->get_no_confs()-1;
456 }
457 
459 {
460  // Smarter algorithm forthcoming in 2.0
461  size_t ret = 0;
463  ret++;
464  reset();
465  return ret;
466 }
467 
469 {
470  if(empty)
471  {
473  return;
474  }
475 
476  partialLProbs[dimNumber] = 0.0;
477 
478  memset(counter, 0, sizeof(int)*dimNumber);
479  recalc(dimNumber-1);
480  counter[0]--;
481 
482  lProbs_ptr = lProbs_ptr_start - 1;
483 }
484 
485 IsoThresholdGenerator::~IsoThresholdGenerator()
486 {
487  delete[] counter;
488  delete[] maxConfsLPSum;
489  if (marginalResultsUnsorted != marginalResults)
490  delete[] marginalResultsUnsorted;
491  dealloc_table(marginalResults, dimNumber);
492  if(marginalOrder != nullptr)
493  delete[] marginalOrder;
494 }
495 
496 
497 /*
498  * ------------------------------------------------------------------------------------------------------------------------
499  */
500 
501 
502 IsoLayeredGenerator::IsoLayeredGenerator(Iso&& iso, int tabSize, int hashSize, bool reorder_marginals, double t_prob_hint)
503 : IsoGenerator(std::move(iso))
504 {
505  counter = new int[dimNumber];
506  maxConfsLPSum = new double[dimNumber-1];
507  currentLThreshold = nextafter(getModeLProb(), -std::numeric_limits<double>::infinity());
508  lastLThreshold = std::numeric_limits<double>::min();
509  marginalResultsUnsorted = new LayeredMarginal*[dimNumber];
510  resetPositions = new const double*[dimNumber];
511 
512  for(int ii=0; ii<dimNumber; ii++)
513  {
514  counter[ii] = 0;
515  marginalResultsUnsorted[ii] = new LayeredMarginal(std::move(*(marginals[ii])), tabSize, hashSize);
516  }
517 
518  if(reorder_marginals && dimNumber > 1)
519  {
520  double* marginal_priorities = new double[dimNumber];
521 
522  /*
523  * We shall now use Gaussian approximations of the marginal multinomial distributions to estimate
524  * how many configurations we shall need to visit from each marginal. This should be approximately
525  * proportional to the volume of the optimal P-ellipsoid of the marginal, which, in turn is defined
526  * by the quantile function of the chi-square distribution plus some modifications.
527  *
528  * We're dropping the constant factor and the (monotonic) exp() transform - these will be used as keys
529  * for sorting, so only the ordering is important.
530  */
531 
532  double K = allDim - dimNumber;
533 
534  double log_R2 = log(InverseChiSquareCDF2(K, t_prob_hint));
535 
536  for(int ii = 0; ii < dimNumber; ii++)
537  {
538  const int i = marginalResultsUnsorted[ii]->get_isotopeNo();
539  if(i <= 1)
540  marginal_priorities[ii] = 0.0;
541  else
542  {
543  double k = static_cast<double>(i - 1);
544  const int n = atomCounts[ii];
545 
546  double sum_lprobs = 0.0;
547  for(int jj = 0; jj < i; jj++)
548  sum_lprobs += marginalResultsUnsorted[ii]->get_lProbs()[jj];
549 
550  double sum_rademacher = 0.0;
551  for(int jj = 1; jj < i; jj++)
552  sum_rademacher += log1p((static_cast<double>(jj)) / static_cast<double>(n));
553 
554  marginal_priorities[ii] = -(sum_lprobs/2.0 + sum_rademacher - lgamma((k+2.0)/2.0) + k/2.0 * (log_R2 + log2pluslogpi + log(n)));
555  }
556  }
557 
558  int* tmpMarginalOrder = new int[dimNumber];
559 
560  for(int ii=0; ii<dimNumber; ii++)
561  tmpMarginalOrder[ii] = ii;
562 
563  TableOrder<double> TO(marginal_priorities);
564 
565  std::sort(tmpMarginalOrder, tmpMarginalOrder + dimNumber, TO);
566  marginalResults = new LayeredMarginal*[dimNumber];
567 
568  for(int ii=0; ii<dimNumber; ii++)
569  marginalResults[ii] = marginalResultsUnsorted[tmpMarginalOrder[ii]];
570 
571  marginalOrder = new int[dimNumber];
572  for(int ii = 0; ii<dimNumber; ii++)
573  marginalOrder[tmpMarginalOrder[ii]] = ii;
574 
575  delete[] tmpMarginalOrder;
576  delete[] marginal_priorities;
577  }
578  else
579  {
580  marginalResults = marginalResultsUnsorted;
581  marginalOrder = nullptr;
582  }
583 
584  lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr();
585 
586  if(dimNumber > 1)
587  maxConfsLPSum[0] = marginalResults[0]->getModeLProb();
588 
589  for(int ii=1; ii<dimNumber-1; ii++)
590  maxConfsLPSum[ii] = maxConfsLPSum[ii-1] + marginalResults[ii]->getModeLProb();
591 
592  lProbs_ptr = lProbs_ptr_start;
593 
594  partialLProbs_second = partialLProbs;
595  partialLProbs_second++;
596 
597  counter[0]--;
598  lProbs_ptr--;
599  lastLThreshold = 10.0;
600  nextLayer(-0.00001);
601 }
602 
603 bool IsoLayeredGenerator::nextLayer(double offset)
604 {
605  size_t first_mrg_size = marginalResults[0]->get_no_confs();
606 
607  if(lastLThreshold < getUnlikeliestPeakLProb())
608  return false;
609 
610  lastLThreshold = currentLThreshold;
611  currentLThreshold += offset;
612 
613  for(int ii=0; ii<dimNumber; ii++)
614  {
615  marginalResults[ii]->extend(currentLThreshold - modeLProb + marginalResults[ii]->getModeLProb());
616  counter[ii] = 0;
617  }
618 
619  lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr(); // vector relocation might have happened
620 
621  lProbs_ptr = lProbs_ptr_start + first_mrg_size - 1;
622 
623  for(int ii=0; ii<dimNumber; ii++)
624  resetPositions[ii] = lProbs_ptr;
625 
626  recalc(dimNumber-1);
627 
628  return true;
629 }
630 
631 bool IsoLayeredGenerator::carry()
632 {
633  // If we reached this point, a carry is needed
634 
635  int idx = 0;
636 
637  int * cntr_ptr = counter;
638 
639  while(idx<dimNumber-1)
640  {
641  *cntr_ptr = 0;
642  idx++;
643  cntr_ptr++;
644  (*cntr_ptr)++;
645  partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
646  if(partialLProbs[idx] + maxConfsLPSum[idx-1] >= currentLThreshold)
647  {
648  partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
649  partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
650  recalc(idx-1);
651  lProbs_ptr = resetPositions[idx];
652 
653  while(*lProbs_ptr <= last_lcfmsv)
654  lProbs_ptr--;
655 
656  for(int ii=0; ii<idx; ii++)
657  resetPositions[ii] = lProbs_ptr;
658 
659  return true;
660  }
661  }
662 
663  return false;
664 }
665 
666 
668 {
669  for(int ii=0; ii<dimNumber; ii++)
670  {
671  counter[ii] = marginalResults[ii]->get_no_confs()-1;
672  partialLProbs[ii] = -std::numeric_limits<double>::infinity();
673  }
674  partialLProbs[dimNumber] = -std::numeric_limits<double>::infinity();
675  lProbs_ptr = lProbs_ptr_start + marginalResults[0]->get_no_confs()-1;
676 }
677 
678 IsoLayeredGenerator::~IsoLayeredGenerator()
679 {
680  delete[] counter;
681  delete[] maxConfsLPSum;
682  delete[] resetPositions;
683  if (marginalResultsUnsorted != marginalResults)
684  delete[] marginalResultsUnsorted;
685  dealloc_table(marginalResults, dimNumber);
686  if(marginalOrder != nullptr)
687  delete[] marginalOrder;
688 }
689 
690 
691 /*
692  * ------------------------------------------------------------------------------------------------------------------------
693  */
694 
695 
696 IsoOrderedGenerator::IsoOrderedGenerator(Iso&& iso, int _tabSize, int _hashSize) :
697 IsoGenerator(std::move(iso), false), allocator(dimNumber, _tabSize)
698 {
699  partialLProbs = &currentLProb;
700  partialMasses = &currentMass;
701  partialProbs = &currentProb;
702 
703  marginalResults = new MarginalTrek*[dimNumber];
704 
705  for(int i = 0; i<dimNumber; i++)
706  marginalResults[i] = new MarginalTrek(std::move(*(marginals[i])), _tabSize, _hashSize);
707 
708  logProbs = new const vector<double>*[dimNumber];
709  masses = new const vector<double>*[dimNumber];
710  marginalConfs = new const vector<int*>*[dimNumber];
711 
712  for(int i = 0; i<dimNumber; i++)
713  {
714  masses[i] = &marginalResults[i]->conf_masses();
715  logProbs[i] = &marginalResults[i]->conf_lprobs();
716  marginalConfs[i] = &marginalResults[i]->confs();
717  }
718 
719  topConf = allocator.newConf();
720  memset(
721  reinterpret_cast<char*>(topConf) + sizeof(double),
722  0,
723  sizeof(int)*dimNumber
724  );
725 
726  *(reinterpret_cast<double*>(topConf)) =
727  combinedSum(
728  getConf(topConf),
729  logProbs,
730  dimNumber
731  );
732 
733  pq.push(topConf);
734 
735 }
736 
737 
739 {
740  dealloc_table<MarginalTrek*>(marginalResults, dimNumber);
741  delete[] logProbs;
742  delete[] masses;
743  delete[] marginalConfs;
744  partialLProbs = nullptr;
745  partialMasses = nullptr;
746  partialProbs = nullptr;
747 }
748 
749 
751 {
752  if(pq.size() < 1)
753  return false;
754 
755 
756  topConf = pq.top();
757  pq.pop();
758 
759  int* topConfIsoCounts = getConf(topConf);
760 
761  currentLProb = *(reinterpret_cast<double*>(topConf));
762  currentMass = combinedSum( topConfIsoCounts, masses, dimNumber );
763  currentProb = exp(currentLProb);
764 
765  ccount = -1;
766  for(int j = 0; j < dimNumber; ++j)
767  {
768  if(marginalResults[j]->probeConfigurationIdx(topConfIsoCounts[j] + 1))
769  {
770  if(ccount == -1)
771  {
772  topConfIsoCounts[j]++;
773  *(reinterpret_cast<double*>(topConf)) = combinedSum(topConfIsoCounts, logProbs, dimNumber);
774  pq.push(topConf);
775  topConfIsoCounts[j]--;
776  ccount = j;
777  }
778  else
779  {
780  void* acceptedCandidate = allocator.newConf();
781  int* acceptedCandidateIsoCounts = getConf(acceptedCandidate);
782  memcpy(acceptedCandidateIsoCounts, topConfIsoCounts, confSize);
783 
784  acceptedCandidateIsoCounts[j]++;
785 
786  *(reinterpret_cast<double*>(acceptedCandidate)) = combinedSum(acceptedCandidateIsoCounts, logProbs, dimNumber);
787 
788  pq.push(acceptedCandidate);
789  }
790  }
791  if(topConfIsoCounts[j] > 0)
792  break;
793  }
794  if(ccount >=0)
795  topConfIsoCounts[ccount]++;
796 
797  return true;
798 }
799 
800 
801 
802 /*
803  * ---------------------------------------------------------------------------------------------------
804  */
805 
806 
807 
808 
809 #if !ISOSPEC_BUILDING_R
810 
811 void printConfigurations(
812  const std::tuple<double*,double*,int*,int>& results,
813  int dimNumber,
814  int* isotopeNumbers
815 ){
816  int m = 0;
817 
818  for(int i=0; i<std::get<3>(results); i++){
819 
820  std::cout << "Mass = " << std::get<0>(results)[i] <<
821  "\tand log-prob = " << std::get<1>(results)[i] <<
822  "\tand prob = " << exp(std::get<1>(results)[i]) <<
823  "\tand configuration =\t";
824 
825 
826  for(int j=0; j<dimNumber; j++){
827  for(int k=0; k<isotopeNumbers[j]; k++ )
828  {
829  std::cout << std::get<2>(results)[m] << " ";
830  m++;
831  }
832  std::cout << '\t';
833  }
834 
835 
836  std::cout << std::endl;
837  }
838 }
839 
840 #endif /* !ISOSPEC_BUILDING_R */
841 
842 
843 
844 
845 } // namespace IsoSpec
846 
IsoSpec::Iso::allDim
int allDim
Definition: isoSpec++.h:71
IsoSpec::IsoGenerator::partialProbs
double * partialProbs
Definition: isoSpec++.h:159
IsoSpec::Iso::confSize
unsigned int confSize
Definition: isoSpec++.h:70
IsoSpec::Iso::getLightestPeakMass
double getLightestPeakMass() const
Get the mass of the lightest peak in the isotopic distribution.
Definition: isoSpec++.cpp:165
IsoSpec::PrecalculatedMarginal
Precalculated Marginal class.
Definition: marginalTrek++.h:226
IsoSpec::Marginal
The marginal distribution class (a subisotopologue).
Definition: marginalTrek++.h:44
IsoSpec::OrderMarginalsBySizeDecresing
Definition: operators.h:130
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:363
IsoSpec::LayeredMarginal::extend
bool extend(double new_threshold)
Extend the set of computed subisotopologues to those above the new threshold.
Definition: marginalTrek++.cpp:492
IsoSpec::IsoLayeredGenerator::recalc
ISOSPEC_FORCE_INLINE void recalc(int idx)
Recalculate the current partial log-probabilities, masses, and probabilities.
Definition: isoSpec++.h:467
IsoSpec::Iso::dimNumber
int dimNumber
Definition: isoSpec++.h:67
IsoSpec::PrecalculatedMarginal::get_lProbs_ptr
const double * get_lProbs_ptr() const
Get the table of the log-probabilities of subisotopologues.
Definition: marginalTrek++.h:287
IsoSpec::IsoThresholdGenerator::count_confs
size_t count_confs()
Definition: isoSpec++.cpp:458
IsoSpec::Iso::getHeaviestPeakMass
double getHeaviestPeakMass() const
Get the mass of the heaviest peak in the isotopic distribution.
Definition: isoSpec++.cpp:173
IsoSpec::IsoThresholdGenerator::advanceToNextConfiguration
ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
Advance to the next, not yet visited, most probable isotopologue.
Definition: isoSpec++.h:315
IsoSpec
Definition: allocator.cpp:21
IsoSpec::Marginal::getModeLProb
double getModeLProb() const
Get the log-probability of the mode subisotopologue.
Definition: marginalTrek++.h:118
IsoSpec::IsoOrderedGenerator::~IsoOrderedGenerator
virtual ~IsoOrderedGenerator()
Destructor.
Definition: isoSpec++.cpp:738
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:354
IsoSpec::Iso::getModeLProb
double getModeLProb() const
Get the log-probability of the mode-configuration (if there are many modes, they share this value).
Definition: isoSpec++.h:128
IsoSpec::Iso
The Iso class for the calculation of the isotopic distribution.
Definition: isoSpec++.h:53
IsoSpec::Iso::getUnlikeliestPeakLProb
double getUnlikeliestPeakLProb() const
Get the logprobability of the least probable subisotopologue.
Definition: isoSpec++.cpp:189
IsoSpec::IsoOrderedGenerator::IsoOrderedGenerator
IsoOrderedGenerator(Iso &&iso, int _tabSize=1000, int _hashSize=1000)
The move-contstructor.
Definition: isoSpec++.cpp:696
IsoSpec::Iso::~Iso
virtual ~Iso()
Destructor.
Definition: isoSpec++.cpp:153
IsoSpec::Iso::isotopeNumbers
int * isotopeNumbers
Definition: isoSpec++.h:68
IsoSpec::Iso::getMonoisotopicPeakMass
double getMonoisotopicPeakMass() const
Definition: isoSpec++.cpp:181
IsoSpec::IsoGenerator::partialMasses
double * partialMasses
Definition: isoSpec++.h:158
IsoSpec::Iso::modeLProb
double modeLProb
Definition: isoSpec++.h:73
IsoSpec::IsoThresholdGenerator::reset
void reset()
Definition: isoSpec++.cpp:468
IsoSpec::IsoOrderedGenerator::advanceToNextConfiguration
bool advanceToNextConfiguration() override final
Advance to the next, not yet visited, most probable isotopologue.
Definition: isoSpec++.cpp:750
IsoSpec::PrecalculatedMarginal::get_no_confs
unsigned int get_no_confs() const
Get the number of precomputed subisotopologues.
Definition: marginalTrek++.h:307
IsoSpec::IsoGenerator::IsoGenerator
IsoGenerator(Iso &&iso, bool alloc_partials=true)
Move constructor.
Definition: isoSpec++.cpp:335
IsoSpec::Iso::marginals
Marginal ** marginals
Definition: isoSpec++.h:72
IsoSpec::IsoThresholdGenerator::IsoThresholdGenerator
IsoThresholdGenerator(Iso &&iso, double _threshold, bool _absolute=true, int _tabSize=1000, int _hashSize=1000, bool reorder_marginals=true)
The move-constructor.
Definition: isoSpec++.cpp:368
IsoSpec::IsoGenerator::partialLProbs
double * partialLProbs
Definition: isoSpec++.h:157
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:357
IsoSpec::IsoGenerator::~IsoGenerator
virtual ~IsoGenerator()
Destructor.
Definition: isoSpec++.cpp:350
IsoSpec::Iso::addElement
void addElement(int atomCount, int noIsotopes, const double *isotopeMasses, const double *isotopeProbabilities)
Add an element to the molecule. Note: this method can only be used BEFORE Iso is used to construct an...
Definition: isoSpec++.cpp:229
IsoSpec::Marginal::get_isotopeNo
int get_isotopeNo() const
Get the number of isotopes of the investigated element.
Definition: marginalTrek++.h:91
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:351
IsoSpec::IsoThresholdGenerator::terminate_search
void terminate_search()
Block the subsequent search of isotopologues.
Definition: isoSpec++.cpp:447
IsoSpec::Iso::getModeMass
double getModeMass() const
Get the mass of the mode-configuration (if there are many modes, it is undefined which one will be se...
Definition: isoSpec++.cpp:197
IsoSpec::IsoGenerator
The generator of isotopologues.
Definition: isoSpec++.h:154
IsoSpec::MarginalTrek
The marginal distribution class (a subisotopologue).
Definition: marginalTrek++.h:154
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:348
IsoSpec::Iso::getTheoreticalAverageMass
double getTheoreticalAverageMass() const
Get the theoretical average mass of the molecule.
Definition: isoSpec++.cpp:205
IsoSpec::Iso::atomCounts
int * atomCounts
Definition: isoSpec++.h:69
IsoSpec::IsoLayeredGenerator::terminate_search
void terminate_search()
Block the subsequent search of isotopologues.
Definition: isoSpec++.cpp:667