23 #include <unordered_map> 24 #include <unordered_set> 34 #include "marginalTrek++.h" 36 #include "allocator.h" 37 #include "operators.h" 39 #include "element_tables.h" 55 Conf
initialConfigure(
const int atomCnt,
const int isotopeNo,
const double* probs,
const double* lprobs)
61 Conf res =
new int[isotopeNo];
64 for(
int i = 0; i < isotopeNo; ++i )
65 res[i] =
int( atomCnt * probs[i] ) + 1;
70 for(
int i = 0; i < isotopeNo; ++i) s += res[i];
72 int diff = atomCnt - s;
81 int i = 0, coordDiff = 0;
84 coordDiff = res[i] - diff;
92 diff = abs(coordDiff);
100 double LP = unnormalized_logProb(res, lprobs, isotopeNo);
106 for(
int ii = 0; ii<isotopeNo; ii++)
107 for(
int jj = 0; jj<isotopeNo; jj++)
108 if(ii != jj && res[ii] > 0)
112 NLP = unnormalized_logProb(res, lprobs, isotopeNo);
113 if(NLP>LP || (NLP==LP && ii>jj))
132 #if !ISOSPEC_BUILDING_R 133 void printMarginal(
const std::tuple<double*,double*,int*,int>& results,
int dim)
135 for(
int i=0; i<std::get<3>(results); i++){
137 std::cout <<
"Mass = " << std::get<0>(results)[i] <<
138 " log-prob =\t" << std::get<1>(results)[i] <<
139 " prob =\t" << exp(std::get<1>(results)[i]) <<
140 "\tand configuration =\t";
142 for(
int j=0; j<dim; j++) std::cout << std::get<2>(results)[i*dim + j] <<
" ";
144 std::cout << std::endl;
157 int curr_method = fegetround();
158 fesetround(FE_UPWARD);
159 double* ret =
new double[isoNo];
162 for(
int i = 0; i < isoNo; i++)
164 ret[i] = log(probs[i]);
165 for(
int j=0; j<ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES; j++)
166 if(elem_table_probability[j] == probs[i])
168 ret[i] = elem_table_log_probability[j];
172 fesetround(curr_method);
176 double get_loggamma_nominator(
int x)
179 int curr_method = fegetround();
180 fesetround(FE_UPWARD);
181 double ret = lgamma(x+1);
182 fesetround(curr_method);
188 const double* _masses,
189 const double* _probs,
194 isotopeNo(_isotopeNo),
196 atom_masses(array_copy<double>(_masses, _isotopeNo)),
198 loggamma_nominator(get_loggamma_nominator(_atomCnt)),
200 mode_lprob(loggamma_nominator+unnormalized_logProb(mode_conf, atom_lProbs, isotopeNo)),
201 mode_mass(mass(mode_conf, atom_masses, isotopeNo)),
202 mode_prob(exp(mode_lprob)),
203 smallest_lprob(atomCnt * *
std::min_element(atom_lProbs, atom_lProbs+isotopeNo))
207 if(ISOSPEC_G_FACT_TABLE_SIZE-1 <=
atomCnt)
208 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));
210 if(_probs[ii] <= 0.0 || _probs[ii] > 1.0)
211 throw std::invalid_argument(
"All isotope probabilities p must fulfill: 0.0 < p <= 1.0");
224 disowned(other.disowned),
236 other.disowned =
true;
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++)
279 visited(hashSize,keyHasher,equalizer),
285 int* initialConf = allocator.makeCopy(
mode_conf);
287 pq.push(initialConf);
288 visited[initialConf] = 0;
298 bool MarginalTrek::add_next_conf()
304 if(pq.size() < 1)
return false;
306 Conf topConf = pq.top();
309 visited[topConf] = current_count;
311 _confs.push_back(topConf);
313 double logprob =
logProb(topConf);
314 _conf_lprobs.push_back(logprob);
317 totalProb.add( exp( logprob ) );
319 for(
unsigned int i = 0; i <
isotopeNo; ++i )
321 for(
unsigned int j = 0; j <
isotopeNo; ++j )
325 if( i != j && topConf[j] > 0 ){
326 copyConf(topConf, candidate, isotopeNo);
332 if( visited.count( candidate ) == 0 )
334 Conf acceptedCandidate = allocator.makeCopy(candidate);
335 pq.push(acceptedCandidate);
337 visited[acceptedCandidate] = 0;
350 for(
unsigned int i=0; i<_conf_lprobs.size(); i++)
352 s.add(_conf_lprobs[i]);
353 if(s.get() >= cutoff)
362 while(totalProb.get() < cutoff && add_next_conf()) {}
363 return _conf_lprobs.size();
367 MarginalTrek::~MarginalTrek()
387 std::unordered_set<Conf,KeyHasher,ConfEqual> visited(hashSize,keyHasher,equalizer);
389 Conf currentConf = allocator.makeCopy(
mode_conf);
390 if(
logProb(currentConf) >= lCutOff)
394 auto tmp = allocator.makeCopy(currentConf);
395 configurations.push_back(tmp);
399 unsigned int idx = 0;
401 while(idx < configurations.size())
403 memcpy(currentConf, configurations[idx],
sizeof(
int)*
isotopeNo);
405 for(
unsigned int ii = 0; ii <
isotopeNo; ii++ )
406 for(
unsigned int jj = 0; jj <
isotopeNo; jj++ )
407 if( ii != jj && currentConf[jj] > 0)
412 if (visited.count(currentConf) == 0 &&
logProb(currentConf) >= lCutOff)
416 auto tmp = allocator.makeCopy(currentConf);
418 configurations.push_back(tmp);
431 std::sort(configurations.begin(), configurations.end(), orderMarginal);
434 confs = &configurations[0];
435 no_confs = configurations.size();
436 lProbs =
new double[no_confs+1];
437 probs =
new double[no_confs];
438 masses =
new double[no_confs];
441 for(
unsigned int ii=0; ii < no_confs; ii++)
443 lProbs[ii] =
logProb(confs[ii]);
444 probs[ii] = exp(lProbs[ii]);
447 lProbs[no_confs] = -std::numeric_limits<double>::infinity();
453 if(lProbs !=
nullptr)
455 if(masses !=
nullptr)
double getLightestConfMass() const
Get the mass of the lightest subisotopologue.
virtual ~PrecalculatedMarginal()
Destructor.
The marginal distribution class (a subisotopologue).
int processUntilCutoff(double cutoff)
Calculate subisotopologues with probability above or equal to the cut-off.
MarginalTrek(Marginal &&m, int tabSize=1000, int hashSize=1000)
Move constructor: specializes the Marginal class.
const unsigned int isotopeNo
const unsigned int atomCnt
virtual ~Marginal()
Destructor.
const double *const atom_lProbs
double getHeaviestConfMass() const
Get the mass of the heaviest subisotopologue.
const double loggamma_nominator
Conf initialConfigure(const int atomCnt, const int isotopeNo, const double *probs, const double *lprobs)
Find one of the most probable subisotopologues.
const double smallest_lprob
Marginal(const double *_masses, const double *_probs, int _isotopeNo, int _atomCnt)
Class constructor.
PrecalculatedMarginal(Marginal &&m, double lCutOff, bool sort=true, int tabSize=1000, int hashSize=1000)
The move constructor (disowns the Marginal).
const double *const atom_masses
double * getMLogProbs(const double *probs, int isoNo)
double logProb(Conf conf) const
Calculate the log-probability of a given subisotopologue.