IsoSpec  1.95
isoSpec++.cpp
1 /*
2  * Copyright (C) 2015-2018 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  int _dimNumber,
53  const int* _isotopeNumbers,
54  const int* _atomCounts,
55  const double* const * _isotopeMasses,
56  const double* const * _isotopeProbabilities
57 ) :
58 disowned(false),
59 dimNumber(_dimNumber),
60 isotopeNumbers(array_copy<int>(_isotopeNumbers, _dimNumber)),
61 atomCounts(array_copy<int>(_atomCounts, _dimNumber)),
62 confSize(_dimNumber * sizeof(int)),
63 allDim(0),
64 marginals(nullptr),
65 modeLProb(0.0)
66 {
67  try{
68  setupMarginals(_isotopeMasses, _isotopeProbabilities);
69  }
70  catch(...)
71  {
72  delete[] isotopeNumbers;
73  delete[] atomCounts;
74  throw;
75  }
76 }
77 
78 Iso::Iso(Iso&& other) :
79 disowned(other.disowned),
80 dimNumber(other.dimNumber),
81 isotopeNumbers(other.isotopeNumbers),
82 atomCounts(other.atomCounts),
83 confSize(other.confSize),
84 allDim(other.allDim),
85 marginals(other.marginals),
86 modeLProb(other.modeLProb)
87 {
88  other.disowned = true;
89 }
90 
91 
92 Iso::Iso(const Iso& other, bool fullcopy) :
93 disowned(fullcopy ? throw std::logic_error("Not implemented") : true),
94 dimNumber(other.dimNumber),
95 isotopeNumbers(fullcopy ? array_copy<int>(other.isotopeNumbers, dimNumber) : other.isotopeNumbers),
96 atomCounts(fullcopy ? array_copy<int>(other.atomCounts, dimNumber) : other.atomCounts),
97 confSize(other.confSize),
98 allDim(other.allDim),
99 marginals(fullcopy ? throw std::logic_error("Not implemented") : other.marginals),
100 modeLProb(other.modeLProb)
101 {}
102 
103 
104 inline void Iso::setupMarginals(const double* const * _isotopeMasses, const double* const * _isotopeProbabilities)
105 {
106  if (marginals == nullptr)
107  {
108  int ii = 0;
109  try
110  {
111  marginals = new Marginal*[dimNumber];
112  while(ii < dimNumber)
113  {
114  allDim += isotopeNumbers[ii];
115  marginals[ii] = new Marginal(
116  _isotopeMasses[ii],
117  _isotopeProbabilities[ii],
118  isotopeNumbers[ii],
119  atomCounts[ii]
120  );
121  modeLProb += marginals[ii]->getModeLProb();
122  ii++;
123  }
124  }
125  catch(...)
126  {
127  ii--;
128  while(ii >= 0)
129  {
130  delete marginals[ii];
131  ii--;
132  }
133  delete[] marginals;
134  marginals = nullptr;
135  throw;
136  }
137  }
138 
139 }
140 
142 {
143  if(!disowned)
144  {
145  if (marginals != nullptr)
146  dealloc_table(marginals, dimNumber);
147  delete[] isotopeNumbers;
148  delete[] atomCounts;
149  }
150 }
151 
152 
154 {
155  double mass = 0.0;
156  for (int ii=0; ii<dimNumber; ii++)
157  mass += marginals[ii]->getLightestConfMass();
158  return mass;
159 }
160 
162 {
163  double mass = 0.0;
164  for (int ii=0; ii<dimNumber; ii++)
165  mass += marginals[ii]->getHeaviestConfMass();
166  return mass;
167 }
168 
169 
170 
171 Iso::Iso(const char* formula) :
172 disowned(false),
173 allDim(0),
174 marginals(nullptr),
175 modeLProb(0.0)
176 {
177  std::vector<const double*> isotope_masses;
178  std::vector<const double*> isotope_probabilities;
179 
180  dimNumber = parse_formula(formula, isotope_masses, isotope_probabilities, &isotopeNumbers, &atomCounts, &confSize);
181 
182  setupMarginals(isotope_masses.data(), isotope_probabilities.data());
183 }
184 
185 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)
186 {
187  // This function is NOT guaranteed to be secure against malicious input. It should be used only for debugging.
188  size_t slen = strlen(formula);
189  // 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
190  // 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
191  // little bit.
192  std::vector<std::pair<const char*, size_t> > elements;
193  std::vector<int> numbers;
194 
195  if(slen == 0)
196  throw invalid_argument("Invalid formula: can't be empty");
197 
198  if(!isdigit(formula[slen-1]))
199  throw invalid_argument("Invalid formula: every element must be followed by a number - write H2O1 and not H2O for water");
200 
201  for(size_t ii=0; ii<slen; ii++)
202  if(!isdigit(formula[ii]) && !isalpha(formula[ii]))
203  throw invalid_argument("Ivalid formula: contains invalid (non-digit, non-alpha) character");
204 
205  size_t position = 0;
206  size_t elem_end = 0;
207  size_t digit_end = 0;
208 
209  while(position < slen)
210  {
211  elem_end = position;
212  while(isalpha(formula[elem_end]))
213  elem_end++;
214  digit_end = elem_end;
215  while(isdigit(formula[digit_end]))
216  digit_end++;
217  elements.emplace_back(&formula[position], elem_end-position);
218  numbers.push_back(atoi(&formula[elem_end]));
219  position = digit_end;
220  }
221 
222  std::vector<int> element_indexes;
223 
224  for (unsigned int i=0; i<elements.size(); i++)
225  {
226  int idx = -1;
227  for(int j=0; j<ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES; j++)
228  {
229  if ((strlen(elem_table_symbol[j]) == elements[i].second) && (strncmp(elements[i].first, elem_table_symbol[j], elements[i].second) == 0))
230  {
231  idx = j;
232  break;
233  }
234  }
235  if(idx < 0)
236  throw invalid_argument("Invalid formula");
237  element_indexes.push_back(idx);
238  }
239 
240  vector<int> _isotope_numbers;
241 
242  for(vector<int>::iterator it = element_indexes.begin(); it != element_indexes.end(); ++it)
243  {
244  int num = 0;
245  int at_idx = *it;
246  int atomicNo = elem_table_atomicNo[at_idx];
247  while(at_idx < ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES && elem_table_atomicNo[at_idx] == atomicNo)
248  {
249  at_idx++;
250  num++;
251  }
252  _isotope_numbers.push_back(num);
253  }
254 
255  for(vector<int>::iterator it = element_indexes.begin(); it != element_indexes.end(); ++it)
256  {
257  isotope_masses.push_back(&elem_table_mass[*it]);
258  isotope_probabilities.push_back(&elem_table_probability[*it]);
259  };
260 
261  const unsigned int dimNumber = elements.size();
262 
263  *isotopeNumbers = array_copy<int>(_isotope_numbers.data(), dimNumber);
264  *atomCounts = array_copy<int>(numbers.data(), dimNumber);
265  *confSize = dimNumber * sizeof(int);
266 
267  return dimNumber;
268 
269 }
270 
271 
272 /*
273  * ----------------------------------------------------------------------------------------------------------
274  */
275 
276 
277 
278 IsoGenerator::IsoGenerator(Iso&& iso, bool alloc_partials) :
279  Iso(std::move(iso)),
280  partialLProbs(alloc_partials ? new double[dimNumber+1] : nullptr),
281  partialMasses(alloc_partials ? new double[dimNumber+1] : nullptr),
282  partialProbs(alloc_partials ? new double[dimNumber+1] : nullptr)
283 {
284  if(alloc_partials)
285  {
286  partialLProbs[dimNumber] = 0.0;
287  partialMasses[dimNumber] = 0.0;
288  partialProbs[dimNumber] = 1.0;
289  }
290 }
291 
292 
294 {
295  if(partialLProbs != nullptr)
296  delete[] partialLProbs;
297  if(partialMasses != nullptr)
298  delete[] partialMasses;
299  if(partialProbs != nullptr)
300  delete[] partialProbs;
301 }
302 
303 
304 
305 /*
306  * ----------------------------------------------------------------------------------------------------------
307  */
308 
309 
310 
311 IsoThresholdGenerator::IsoThresholdGenerator(Iso&& iso, double _threshold, bool _absolute, int tabSize, int hashSize, bool reorder_marginals)
312 : IsoGenerator(std::move(iso)),
313 Lcutoff(_threshold <= 0.0 ? std::numeric_limits<double>::lowest() : (_absolute ? log(_threshold) : log(_threshold) + modeLProb))
314 {
315  counter = new int[dimNumber];
316  maxConfsLPSum = new double[dimNumber-1];
317  marginalResultsUnsorted = new PrecalculatedMarginal*[dimNumber];
318 
319  empty = false;
320 
321  for(int ii=0; ii<dimNumber; ii++)
322  {
323  counter[ii] = 0;
324  marginalResultsUnsorted[ii] = new PrecalculatedMarginal(std::move(*(marginals[ii])),
325  Lcutoff - modeLProb + marginals[ii]->getModeLProb(),
326  true,
327  tabSize,
328  hashSize);
329 
330  if(!marginalResultsUnsorted[ii]->inRange(0))
331  empty = true;
332  }
333 
334  if(reorder_marginals)
335  {
336  OrderMarginalsBySizeDecresing comparator(marginalResultsUnsorted);
337  int* tmpMarginalOrder = new int[dimNumber];
338 
339  for(int ii=0; ii<dimNumber; ii++)
340  tmpMarginalOrder[ii] = ii;
341 
342  std::sort(tmpMarginalOrder, tmpMarginalOrder + dimNumber, comparator);
343  marginalResults = new PrecalculatedMarginal*[dimNumber];
344 
345  for(int ii=0; ii<dimNumber; ii++)
346  marginalResults[ii] = marginalResultsUnsorted[tmpMarginalOrder[ii]];
347 
348  marginalOrder = new int[dimNumber];
349  for(int ii = 0; ii<dimNumber; ii++)
350  marginalOrder[tmpMarginalOrder[ii]] = ii;
351 
352  delete[] tmpMarginalOrder;
353 
354  }
355  else
356  {
357  marginalResults = marginalResultsUnsorted;
358  marginalOrder = nullptr;
359  }
360 
361  lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr();
362 
363  if(dimNumber > 1)
364  maxConfsLPSum[0] = marginalResults[0]->getModeLProb();
365 
366  for(int ii=1; ii<dimNumber-1; ii++)
367  maxConfsLPSum[ii] = maxConfsLPSum[ii-1] + marginalResults[ii]->getModeLProb();
368 
369  lProbs_ptr = lProbs_ptr_start;
370 
371  partialLProbs_second = partialLProbs;
372  partialLProbs_second++;
373 
374  if(!empty)
375  {
376  recalc(dimNumber-1);
377  counter[0]--;
378  lProbs_ptr--;
379  }
380  else
381  {
383  lcfmsv = std::numeric_limits<double>::infinity();
384  }
385 
386 
387 
388 }
389 
391 {
392  for(int ii=0; ii<dimNumber; ii++)
393  {
394  counter[ii] = marginalResults[ii]->get_no_confs()-1;
395  partialLProbs[ii] = -std::numeric_limits<double>::infinity();
396  }
397  partialLProbs[dimNumber] = -std::numeric_limits<double>::infinity();
398  lProbs_ptr = lProbs_ptr_start + marginalResults[0]->get_no_confs()-1;
399 }
400 
402 {
403  // Smarter algorithm forthcoming in 2.0
404  size_t ret = 0;
406  ret++;
407  reset();
408  return ret;
409 }
410 
412 {
413  if(empty)
414  {
416  return;
417  }
418 
419  partialLProbs[dimNumber] = 0.0;
420 
421  memset(counter, 0, sizeof(int)*dimNumber);
422  recalc(dimNumber-1);
423  counter[0]--;
424 
425  lProbs_ptr = lProbs_ptr_start - 1;
426 }
427 
428 /*
429  * ------------------------------------------------------------------------------------------------------------------------
430  */
431 
432 IsoOrderedGenerator::IsoOrderedGenerator(Iso&& iso, int _tabSize, int _hashSize) :
433 IsoGenerator(std::move(iso), false), allocator(dimNumber, _tabSize)
434 {
435  partialLProbs = &currentLProb;
436  partialMasses = &currentMass;
437  partialProbs = &currentProb;
438 
439  marginalResults = new MarginalTrek*[dimNumber];
440 
441  for(int i = 0; i<dimNumber; i++)
442  marginalResults[i] = new MarginalTrek(std::move(*(marginals[i])), _tabSize, _hashSize);
443 
444  logProbs = new const vector<double>*[dimNumber];
445  masses = new const vector<double>*[dimNumber];
446  marginalConfs = new const vector<int*>*[dimNumber];
447 
448  for(int i = 0; i<dimNumber; i++)
449  {
450  masses[i] = &marginalResults[i]->conf_masses();
451  logProbs[i] = &marginalResults[i]->conf_lprobs();
452  marginalConfs[i] = &marginalResults[i]->confs();
453  }
454 
455  topConf = allocator.newConf();
456  memset(
457  reinterpret_cast<char*>(topConf) + sizeof(double),
458  0,
459  sizeof(int)*dimNumber
460  );
461 
462  *(reinterpret_cast<double*>(topConf)) =
463  combinedSum(
464  getConf(topConf),
465  logProbs,
466  dimNumber
467  );
468 
469  pq.push(topConf);
470 
471 }
472 
473 
475 {
476  dealloc_table<MarginalTrek*>(marginalResults, dimNumber);
477  delete[] logProbs;
478  delete[] masses;
479  delete[] marginalConfs;
480  partialLProbs = nullptr;
481  partialMasses = nullptr;
482  partialProbs = nullptr;
483 }
484 
485 
487 {
488  if(pq.size() < 1)
489  return false;
490 
491 
492  topConf = pq.top();
493  pq.pop();
494 
495  int* topConfIsoCounts = getConf(topConf);
496 
497  currentLProb = *(reinterpret_cast<double*>(topConf));
498  currentMass = combinedSum( topConfIsoCounts, masses, dimNumber );
499  currentProb = exp(currentLProb);
500 
501  ccount = -1;
502  for(int j = 0; j < dimNumber; ++j)
503  {
504  if(marginalResults[j]->probeConfigurationIdx(topConfIsoCounts[j] + 1))
505  {
506  if(ccount == -1)
507  {
508  topConfIsoCounts[j]++;
509  *(reinterpret_cast<double*>(topConf)) = combinedSum(topConfIsoCounts, logProbs, dimNumber);
510  pq.push(topConf);
511  topConfIsoCounts[j]--;
512  ccount = j;
513  }
514  else
515  {
516  void* acceptedCandidate = allocator.newConf();
517  int* acceptedCandidateIsoCounts = getConf(acceptedCandidate);
518  memcpy(acceptedCandidateIsoCounts, topConfIsoCounts, confSize);
519 
520  acceptedCandidateIsoCounts[j]++;
521 
522  *(reinterpret_cast<double*>(acceptedCandidate)) = combinedSum(acceptedCandidateIsoCounts, logProbs, dimNumber);
523 
524  pq.push(acceptedCandidate);
525  }
526  }
527  if(topConfIsoCounts[j] > 0)
528  break;
529  }
530  if(ccount >=0)
531  topConfIsoCounts[ccount]++;
532 
533  return true;
534 }
535 
536 
537 
538 /*
539  * ---------------------------------------------------------------------------------------------------
540  */
541 
542 
543 
544 
545 #if !ISOSPEC_BUILDING_R
546 
547 void printConfigurations(
548  const std::tuple<double*,double*,int*,int>& results,
549  int dimNumber,
550  int* isotopeNumbers
551 ){
552  int m = 0;
553 
554  for(int i=0; i<std::get<3>(results); i++){
555 
556  std::cout << "Mass = " << std::get<0>(results)[i] <<
557  "\tand log-prob = " << std::get<1>(results)[i] <<
558  "\tand prob = " << exp(std::get<1>(results)[i]) <<
559  "\tand configuration =\t";
560 
561 
562  for(int j=0; j<dimNumber; j++){
563  for(int k=0; k<isotopeNumbers[j]; k++ )
564  {
565  std::cout << std::get<2>(results)[m] << " ";
566  m++;
567  }
568  std::cout << '\t';
569  }
570 
571 
572  std::cout << std::endl;
573  }
574 }
575 
576 #endif /* !ISOSPEC_BUILDING_R */
577 
578 
579 
580 IsoLayeredGenerator::IsoLayeredGenerator( Iso&& iso,
581  double _targetCoverage,
582  double _percentageToExpand,
583  int _tabSize,
584  int _hashSize,
585  bool trim
586 ) : IsoGenerator(std::move(iso)),
587 allocator(dimNumber, _tabSize),
588 candidate(new int[dimNumber]),
589 targetCoverage(_targetCoverage >= 1.0 ? 10000.0 : _targetCoverage), // If the user wants the entire spectrum,
590  // give it to him - and make sure we don't terminate
591  // early because of rounding errors
592 percentageToExpand(_percentageToExpand),
593 do_trim(trim),
594 layers(0),
595 generator_position(-1)
596 {
597  marginalResults = new MarginalTrek*[dimNumber];
598 
599  for(int i = 0; i<dimNumber; i++)
600  marginalResults[i] = new MarginalTrek(std::move(*(marginals[i])), _tabSize, _hashSize);
601 
602  logProbs = new const vector<double>*[dimNumber];
603  masses = new const vector<double>*[dimNumber];
604  marginalConfs = new const vector<int*>*[dimNumber];
605 
606  for(int i = 0; i<dimNumber; i++)
607  {
608  masses[i] = &marginalResults[i]->conf_masses();
609  logProbs[i] = &marginalResults[i]->conf_lprobs();
610  marginalConfs[i] = &marginalResults[i]->confs();
611  }
612 
613  void* topConf = allocator.newConf();
614  memset(reinterpret_cast<char*>(topConf) + sizeof(double), 0, sizeof(int)*dimNumber);
615  *(reinterpret_cast<double*>(topConf)) = combinedSum(getConf(topConf), logProbs, dimNumber);
616 
617  current = new std::vector<void*>();
618  next = new std::vector<void*>();
619 
620  current->push_back(topConf);
621 
622  lprobThr = (*reinterpret_cast<double*>(topConf));
623 
624  if(targetCoverage > 0.0)
625  while(advanceToNextLayer()) {};
626 }
627 
628 
629 IsoLayeredGenerator::~IsoLayeredGenerator()
630 {
631  if(current != nullptr)
632  delete current;
633  if(next != nullptr)
634  delete next;
635  delete[] logProbs;
636  delete[] masses;
637  delete[] marginalConfs;
638  delete[] candidate;
639  dealloc_table(marginalResults, dimNumber);
640 }
641 
642 bool IsoLayeredGenerator::advanceToNextLayer()
643 {
644  layers += 1;
645  double maxFringeLprob = -std::numeric_limits<double>::infinity();
646 
647  if(current == nullptr)
648  return false;
649  int accepted_in_this_layer = 0;
650  Summator prob_in_this_layer(totalProb);
651 
652  void* topConf;
653 
654  while(current->size() > 0)
655  {
656  topConf = current->back();
657  current->pop_back();
658 
659  double top_lprob = getLProb(topConf);
660 
661  if(top_lprob >= lprobThr)
662  {
663  newaccepted.push_back(topConf);
664  accepted_in_this_layer++;
665  prob_in_this_layer.add(exp(top_lprob));
666  }
667  else
668  {
669  next->push_back(topConf);
670  continue;
671  }
672 
673  int* topConfIsoCounts = getConf(topConf);
674 
675  for(int j = 0; j < dimNumber; ++j)
676  {
677  // candidate cannot refer to a position that is
678  // out of range of the stored marginal distribution.
679  if(marginalResults[j]->probeConfigurationIdx(topConfIsoCounts[j] + 1))
680  {
681  memcpy(candidate, topConfIsoCounts, confSize);
682  candidate[j]++;
683 
684  void* acceptedCandidate = allocator.newConf();
685  int* acceptedCandidateIsoCounts = getConf(acceptedCandidate);
686  memcpy( acceptedCandidateIsoCounts, candidate, confSize);
687 
688  double newConfProb = combinedSum(
689  candidate,
690  logProbs,
691  dimNumber
692  );
693 
694 
695 
696  *(reinterpret_cast<double*>(acceptedCandidate)) = newConfProb;
697 
698  if(newConfProb >= lprobThr)
699  current->push_back(acceptedCandidate);
700  else
701  {
702  next->push_back(acceptedCandidate);
703  if(newConfProb > maxFringeLprob)
704  maxFringeLprob = top_lprob;
705  }
706  }
707  if(topConfIsoCounts[j] > 0)
708  break;
709  }
710  }
711 
712  if(next == nullptr || next->size() < 1)
713  return false;
714  else
715  {
716  if(prob_in_this_layer.get() < targetCoverage)
717  {
718  std::vector<void*>* nnew = current;
719  nnew->clear();
720  current = next;
721  next = nnew;
722  int howmany = floor(current->size()*percentageToExpand);
723  lprobThr = getLProb(quickselect(current->data(), howmany, 0, current->size()));
724  totalProb = prob_in_this_layer;
725  }
726  else
727  {
728  delete next;
729  next = nullptr;
730  delete current;
731  current = nullptr;
732  int start = 0;
733  int end = accepted_in_this_layer - 1;
734  void* swapspace;
735 
736  if(do_trim)
737  {
738  void** lastLayer = &(newaccepted.data()[newaccepted.size()-accepted_in_this_layer]);
739 
740  Summator qsprob(totalProb);
741  while(totalProb.get() < targetCoverage)
742  {
743  if(start == end)
744  break;
745 
746  // Partition part
747 
748  int len = end - start;
749 #if ISOSPEC_BUILDING_R
750  int pivot = len/2 + start; // We're very definitely NOT switching to R to use a RNG, and if R sees us use C RNG it complains...
751 #else
752  int pivot = rand() % len + start;
753 #endif
754  void* pval = lastLayer[pivot];
755  double pprob = getLProb(pval);
756  mswap(lastLayer[pivot], lastLayer[end-1]);
757  int loweridx = start;
758  for(int i=start; i<end-1; i++)
759  {
760  if(getLProb(lastLayer[i]) > pprob)
761  {
762  mswap(lastLayer[i], lastLayer[loweridx]);
763  loweridx++;
764  }
765  }
766  mswap(lastLayer[end-1], lastLayer[loweridx]);
767 
768  // Selection part
769 
770  Summator leftProb(qsprob);
771  for(int i=start; i<=loweridx; i++)
772  {
773  leftProb.add(exp(getLProb(lastLayer[i])));
774  }
775  if(leftProb.get() < targetCoverage)
776  {
777  start = loweridx+1;
778  qsprob = leftProb;
779  }
780  else
781  end = loweridx;
782  }
783  int accend = newaccepted.size()-accepted_in_this_layer+start+1;
784 
785  totalProb = qsprob;
786  newaccepted.resize(accend);
787  return true;
788  }
789  else // No trimming
790  {
791  totalProb = prob_in_this_layer;
792  return true;
793  }
794  }
795  }
796  return true;
797 
798 }
799 
801 {
802  generator_position++;
803  if(generator_position < newaccepted.size())
804  {
805  partialLProbs[0] = getLProb(newaccepted[generator_position]);
806  partialMasses[0] = combinedSum(getConf(newaccepted[generator_position]), masses, dimNumber);
807  partialProbs[0] = exp(partialLProbs[0]);
808  return true;
809  }
810  else
811  return false;
812 }
813 
814 
815 
816 
817 } // namespace IsoSpec
818 
IsoSpec::Iso::allDim
int allDim
Definition: isoSpec++.h:71
IsoSpec::IsoGenerator::partialProbs
double * partialProbs
Definition: isoSpec++.h:134
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:153
IsoSpec::PrecalculatedMarginal
Precalculated Marginal class.
Definition: marginalTrek++.h:213
IsoSpec::Marginal
The marginal distribution class (a subisotopologue).
Definition: marginalTrek++.h:45
IsoSpec::OrderMarginalsBySizeDecresing
Definition: operators.h:130
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:274
IsoSpec::IsoLayeredGenerator::advanceToNextConfiguration
bool advanceToNextConfiguration() override final
Advance to the next, not yet visited, most probable isotopologue.
Definition: isoSpec++.cpp:800
IsoSpec::IsoThresholdGenerator::count_confs
size_t count_confs()
Definition: isoSpec++.cpp:401
IsoSpec::Iso::getHeaviestPeakMass
double getHeaviestPeakMass() const
Get the mass of the heaviest peak in the isotopic distribution.
Definition: isoSpec++.cpp:161
IsoSpec::IsoThresholdGenerator::advanceToNextConfiguration
ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
Advance to the next, not yet visited, most probable isotopologue.
Definition: isoSpec++.h:296
IsoSpec
Definition: allocator.cpp:21
IsoSpec::Marginal::getModeLProb
double getModeLProb() const
Get the log-probability of the mode subisotopologue.
Definition: marginalTrek++.h:110
IsoSpec::IsoOrderedGenerator::~IsoOrderedGenerator
virtual ~IsoOrderedGenerator()
Destructor.
Definition: isoSpec++.cpp:474
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:115
IsoSpec::Iso
The Iso class for the calculation of the isotopic distribution.
Definition: isoSpec++.h:52
IsoSpec::quickselect
void * quickselect(void **array, int n, int start, int end)
Quickly select the n'th positional statistic, including the weights.
Definition: misc.cpp:28
IsoSpec::IsoOrderedGenerator::IsoOrderedGenerator
IsoOrderedGenerator(Iso &&iso, int _tabSize=1000, int _hashSize=1000)
The move-contstructor.
Definition: isoSpec++.cpp:432
IsoSpec::Iso::~Iso
virtual ~Iso()
Destructor.
Definition: isoSpec++.cpp:141
IsoSpec::Iso::isotopeNumbers
int * isotopeNumbers
Definition: isoSpec++.h:68
IsoSpec::IsoGenerator::partialMasses
double * partialMasses
Definition: isoSpec++.h:133
IsoSpec::Iso::modeLProb
double modeLProb
Definition: isoSpec++.h:73
IsoSpec::IsoThresholdGenerator::reset
void reset()
Definition: isoSpec++.cpp:411
IsoSpec::IsoOrderedGenerator::advanceToNextConfiguration
bool advanceToNextConfiguration() override final
Advance to the next, not yet visited, most probable isotopologue.
Definition: isoSpec++.cpp:486
IsoSpec::PrecalculatedMarginal::get_no_confs
unsigned int get_no_confs() const
Get the number of precomputed subisotopologues.
Definition: marginalTrek++.h:294
IsoSpec::IsoGenerator::IsoGenerator
IsoGenerator(Iso &&iso, bool alloc_partials=true)
Move constructor.
Definition: isoSpec++.cpp:278
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:311
IsoSpec::IsoGenerator::partialLProbs
double * partialLProbs
Definition: isoSpec++.h:132
IsoSpec::IsoGenerator::~IsoGenerator
virtual ~IsoGenerator()
Destructor.
Definition: isoSpec++.cpp:293
IsoSpec::Iso::Iso
Iso(int _dimNumber, const int *_isotopeNumbers, const int *_atomCounts, const double *const *_isotopeMasses, const double *const *_isotopeProbabilities)
General constructror.
Definition: isoSpec++.cpp:51
IsoSpec::IsoThresholdGenerator::terminate_search
void terminate_search()
Block the subsequent search of isotopologues.
Definition: isoSpec++.cpp:390
IsoSpec::IsoGenerator
The generator of isotopologues.
Definition: isoSpec++.h:129
IsoSpec::MarginalTrek
The marginal distribution class (a subisotopologue).
Definition: marginalTrek++.h:141
IsoSpec::Iso::atomCounts
int * atomCounts
Definition: isoSpec++.h:69
IsoSpec::Iso::disowned
bool disowned
Definition: isoSpec++.h:65