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),
225 isotopeNo(other.isotopeNo),
226 atomCnt(other.atomCnt),
227 atom_masses(other.atom_masses),
228 atom_lProbs(other.atom_lProbs),
229 loggamma_nominator(other.loggamma_nominator),
230 mode_conf(other.mode_conf),
231 mode_lprob(other.mode_lprob),
232 mode_mass(other.mode_mass),
233 mode_prob(other.mode_prob),
234 smallest_lprob(other.smallest_lprob)
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++)
276 keyHasher(isotopeNo),
277 equalizer(isotopeNo),
278 orderMarginal(atom_lProbs, isotopeNo),
279 visited(hashSize,keyHasher,equalizer),
282 candidate(new int[isotopeNo]),
283 allocator(isotopeNo, tabSize)
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 ){
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()
381 allocator(isotopeNo, tabSize)
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)