29 #include "marginalTrek++.h"
31 #include "allocator.h"
32 #include "operators.h"
34 #include "element_tables.h"
58 for(
int i = 0; i < isotopeNo; ++i)
59 res[i] =
static_cast<int>( atomCnt * exp(lprobs[i]) ) + 1;
64 for(
int i = 0; i < isotopeNo; ++i) s += res[i];
66 int diff = atomCnt - s;
78 int coordDiff = res[i] - diff;
86 diff = abs(coordDiff);
94 double LP = unnormalized_logProb(res, lprobs, isotopeNo);
100 for(
int ii = 0; ii < isotopeNo; ii++)
101 for(
int jj = 0; jj < isotopeNo; jj++)
102 if(ii != jj && res[ii] > 0)
106 NLP = unnormalized_logProb(res, lprobs, isotopeNo);
107 if(NLP > LP || (NLP == LP && ii > jj))
129 for(
int ii = 0; ii < isoNo; ii++)
130 if(probs[ii] <= 0.0 || probs[ii] > 1.0)
131 throw std::invalid_argument(
"All isotope probabilities p must fulfill: 0.0 < p <= 1.0");
133 double* ret =
new double[isoNo];
136 for(
int i = 0; i < isoNo; i++)
138 ret[i] = log(probs[i]);
139 for(
int j = 0; j < ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES; j++)
140 if(elem_table_probability[j] == probs[i])
142 ret[i] = elem_table_log_probability[j];
149 double get_loggamma_nominator(
int x)
152 double ret = lgamma(x+1);
156 int verify_atom_cnt(
int atomCnt)
158 #if !ISOSPEC_BUILDING_OPENMS
159 if(ISOSPEC_G_FACT_TABLE_SIZE-1 <= atomCnt)
160 throw std::length_error(
"Subisotopologue too large, size limit (that is, the maximum number of atoms of a single element in a molecule) is: " + std::to_string(ISOSPEC_G_FACT_TABLE_SIZE-1));
166 const double* _masses,
167 const double* _probs,
172 isotopeNo(_isotopeNo),
173 atomCnt(verify_atom_cnt(_atomCnt)),
175 atom_masses(array_copy<double>(_masses, _isotopeNo)),
176 loggamma_nominator(get_loggamma_nominator(_atomCnt)),
183 isotopeNo(other.isotopeNo),
184 atomCnt(other.atomCnt),
185 atom_lProbs(array_copy<double>(other.atom_lProbs, isotopeNo)),
186 atom_masses(array_copy<double>(other.atom_masses, isotopeNo)),
187 loggamma_nominator(other.loggamma_nominator)
204 disowned(other.disowned),
205 isotopeNo(other.isotopeNo),
206 atomCnt(other.atomCnt),
207 atom_lProbs(other.atom_lProbs),
208 atom_masses(other.atom_masses),
209 loggamma_nominator(other.loggamma_nominator)
211 other.disowned =
true;
212 if(other.mode_conf ==
nullptr)
242 void Marginal::setupMode()
244 ISOSPEC_IMPOSSIBLE(
mode_conf !=
nullptr);
252 double ret_mass = std::numeric_limits<double>::infinity();
253 for(
unsigned int ii = 0; ii <
isotopeNo; ii++)
261 double ret_mass = 0.0;
262 for(
unsigned int ii = 0; ii <
isotopeNo; ii++)
270 double found_prob = -std::numeric_limits<double>::infinity();
271 double found_mass = 0.0;
272 for(
unsigned int ii = 0; ii <
isotopeNo; ii++)
284 for(
unsigned int ii = 0; ii <
isotopeNo; ii++)
304 return -std::numeric_limits<double>::infinity();
306 const double i =
static_cast<double>(
isotopeNo);
307 const double k = i - 1.0;
308 const double n =
static_cast<double>(
atomCnt);
310 double sum_lprobs = 0.0;
311 for(
int jj = 0; jj < i; jj++)
314 double log_V_simplex = k * log(n) - lgamma(i);
315 double log_N_simplex = lgamma(n+i) - lgamma(n+1.0) - lgamma(i);
316 double log_V_ellipsoid = (k * (log(n) + logpi + logEllipsoidRadius) + sum_lprobs) * 0.5 - lgamma((i+1)*0.5);
318 return log_N_simplex + log_V_ellipsoid - log_V_simplex;
330 orderMarginal(atom_lProbs, isotopeNo),
332 allocator(isotopeNo, tabSize)
334 int* initialConf = allocator.makeCopy(
mode_conf);
344 bool MarginalTrek::add_next_conf()
350 if(pq.size() < 1)
return false;
353 Conf topConf = pq.top().second;
358 _confs.push_back(topConf);
361 _conf_lprobs.push_back(logprob);
363 for(
unsigned int j = 0; j <
isotopeNo; ++j )
370 for(
unsigned int i = 0; i <
isotopeNo; ++i )
377 Conf acceptedCandidate = allocator.makeCopy(topConf);
379 ++acceptedCandidate[i];
380 --acceptedCandidate[j];
384 pq.push({new_prob, acceptedCandidate});
399 MarginalTrek::~MarginalTrek()
411 allocator(isotopeNo, tabSize)
413 Conf currentConf = allocator.makeCopy(
mode_conf);
414 if(logProb(currentConf) >= lCutOff)
416 configurations.push_back(currentConf);
420 unsigned int idx = 0;
422 std::unique_ptr<double[]> prob_partials(
new double[
isotopeNo]);
423 std::unique_ptr<double[]> prob_part_acc(
new double[
isotopeNo+1]);
426 while(idx < configurations.size())
428 currentConf = configurations[idx];
432 prob_partials[ii] = minuslogFactorial(currentConf[ii]) + currentConf[ii] *
atom_lProbs[ii];
434 for(
unsigned int ii = 0; ii <
isotopeNo; ii++ )
439 if(currentConf[ii] != 0)
441 double prev_partial_ii = prob_partials[ii];
443 prob_partials[ii] = minuslogFactorial(currentConf[ii]) + currentConf[ii] *
atom_lProbs[ii];
445 for(
unsigned int jj = 0; jj <
isotopeNo; jj++ )
447 prob_part_acc[jj+1] = prob_part_acc[jj] + prob_partials[jj];
454 double logp = prob_part_acc[jj] + minuslogFactorial(1+currentConf[jj]) + (1+currentConf[jj]) *
atom_lProbs[jj];
455 for(
size_t kk = jj+1; kk <
isotopeNo; kk++)
456 logp += prob_partials[kk];
460 auto tmp = allocator.makeCopy(currentConf);
462 configurations.push_back(tmp);
463 lProbs.push_back(logp);
467 prob_part_acc[jj+1] = prob_part_acc[jj] + prob_partials[jj];
473 prob_partials[ii] = prev_partial_ii;
481 no_confs = configurations.size();
482 confs = configurations.data();
484 if(sort && no_confs > 0)
486 std::unique_ptr<size_t[]> order_arr(get_inverse_order(lProbs.data(), no_confs));
487 impose_order(order_arr.get(), no_confs, lProbs.data(), confs);
490 probs =
new double[no_confs];
491 masses =
new double[no_confs];
494 for(
unsigned int ii = 0; ii < no_confs; ii++)
496 probs[ii] = exp(lProbs[ii]);
500 lProbs.push_back(-std::numeric_limits<double>::infinity());
506 if(masses !=
nullptr)
519 :
Marginal(std::move(m)), current_threshold(1.0), allocator(isotopeNo, tabSize),
520 equalizer(isotopeNo), keyHasher(isotopeNo)
523 lProbs.push_back(std::numeric_limits<double>::infinity());
525 lProbs.push_back(-std::numeric_limits<double>::infinity());
526 guarded_lProbs = lProbs.data()+1;
540 while(!fringe.empty())
542 Conf currentConf = fringe.back();
545 double opc = fringe_unn_lprobs.back();
547 fringe_unn_lprobs.pop_back();
548 if(opc < new_threshold)
550 new_fringe.push_back(currentConf);
551 new_fringe_unn_lprobs.push_back(opc);
556 configurations.push_back(currentConf);
558 for(
unsigned int ii = 0; ii <
isotopeNo; ii++ )
563 if(currentConf[ii] > 0)
566 for(
unsigned int jj = 0; jj <
isotopeNo; jj++ )
573 Conf nc = allocator.makeCopy(currentConf);
577 if(lpc >= new_threshold)
579 fringe.push_back(nc);
580 fringe_unn_lprobs.push_back(lpc);
584 new_fringe.push_back(nc);
585 new_fringe_unn_lprobs.push_back(lpc);
601 current_threshold = new_threshold;
602 fringe.swap(new_fringe);
603 fringe_unn_lprobs.swap(new_fringe_unn_lprobs);
607 size_t to_sort_size = configurations.size() - probs.size();
610 std::unique_ptr<size_t[]> order_arr(get_inverse_order(lProbs.data()+1+probs.size(), to_sort_size));
611 double* P = lProbs.data()+1+probs.size();
612 Conf* C = configurations.data()+probs.size();
613 size_t* O = order_arr.get();
614 impose_order(O, to_sort_size, P, C);
618 if(probs.capacity() * 2 < configurations.size() + 2)
621 probs.reserve(configurations.size());
622 masses.reserve(configurations.size());
627 for(
unsigned int ii = probs.size(); ii < configurations.size(); ii++)
629 probs.push_back(exp(lProbs[ii+1]));
633 lProbs.push_back(-std::numeric_limits<double>::infinity());
635 guarded_lProbs = lProbs.data()+1;
643 double ret = std::numeric_limits<double>::infinity();
653 double ret = -std::numeric_limits<double>::infinity();