17 #include "fixedEnvelopes.h"
24 FixedEnvelope::FixedEnvelope(
const FixedEnvelope& other) :
25 _masses(array_copy<double>(other._masses, other._confs_no)),
26 _probs(array_copy<double>(other._probs, other._confs_no)),
27 _confs(array_copy_nptr<int>(other._confs, other._confs_no*other.allDim)),
28 _confs_no(other._confs_no),
30 sorted_by_mass(other.sorted_by_mass),
31 sorted_by_prob(other.sorted_by_prob),
32 total_prob(other.total_prob)
35 FixedEnvelope::FixedEnvelope(FixedEnvelope&& other) :
36 _masses(other._masses),
39 _confs_no(other._confs_no),
41 sorted_by_mass(other.sorted_by_mass),
42 sorted_by_prob(other.sorted_by_prob),
43 total_prob(other.total_prob)
45 other._masses =
nullptr;
46 other._probs =
nullptr;
47 other._confs =
nullptr;
49 other.total_prob = 0.0;
52 FixedEnvelope::FixedEnvelope(
double* in_masses,
double* in_probs,
size_t in_confs_no,
bool masses_sorted,
bool probs_sorted,
double _total_prob) :
56 _confs_no(in_confs_no),
58 sorted_by_mass(masses_sorted),
59 sorted_by_prob(probs_sorted),
60 total_prob(_total_prob)
63 FixedEnvelope FixedEnvelope::operator+(
const FixedEnvelope& other)
const
65 double* nprobs =
reinterpret_cast<double*
>(malloc(
sizeof(
double) * (_confs_no+other._confs_no)));
67 throw std::bad_alloc();
68 double* nmasses =
reinterpret_cast<double*
>(malloc(
sizeof(
double) * (_confs_no+other._confs_no)));
69 if(nmasses ==
nullptr)
72 throw std::bad_alloc();
75 memcpy(nprobs, _probs,
sizeof(
double) * _confs_no);
76 memcpy(nmasses, _masses,
sizeof(
double) * _confs_no);
78 memcpy(nprobs+_confs_no, other._probs,
sizeof(
double) * other._confs_no);
79 memcpy(nmasses+_confs_no, other._masses,
sizeof(
double) * other._confs_no);
81 return FixedEnvelope(nmasses, nprobs, _confs_no + other._confs_no);
84 FixedEnvelope FixedEnvelope::operator*(
const FixedEnvelope& other)
const
86 double* nprobs =
reinterpret_cast<double*
>(malloc(
sizeof(
double) * _confs_no * other._confs_no));
88 throw std::bad_alloc();
91 double* nmasses =
reinterpret_cast<double*
>(malloc(
sizeof(
double) * _confs_no * other._confs_no));
92 if(nmasses ==
nullptr)
95 throw std::bad_alloc();
100 for(
size_t ii = 0; ii < _confs_no; ii++)
101 for(
size_t jj = 0; jj < other._confs_no; jj++)
103 nprobs[tgt_idx] = _probs[ii] * other._probs[jj];
104 nmasses[tgt_idx] = _masses[ii] + other._masses[jj];
108 return FixedEnvelope(nmasses, nprobs, tgt_idx);
111 void FixedEnvelope::sort_by_mass()
118 sorted_by_mass =
true;
119 sorted_by_prob =
false;
123 void FixedEnvelope::sort_by_prob()
130 sorted_by_prob =
true;
131 sorted_by_mass =
false;
134 template<
typename T>
void reorder_array(T* arr,
size_t* order,
size_t size,
bool can_destroy =
false)
138 size_t* order_c =
new size_t[size];
139 memcpy(order_c, order,
sizeof(
size_t)*size);
143 for(
size_t ii = 0; ii < size; ii++)
144 while(order[ii] != ii)
146 std::swap(arr[ii], arr[order[ii]]);
147 std::swap(order[order[ii]], order[ii]);
154 void FixedEnvelope::sort_by(
double* order)
156 size_t* indices =
new size_t[_confs_no];
161 for(
size_t ii = 0; ii < _confs_no; ii++)
164 std::sort<size_t*>(indices, indices + _confs_no, TableOrder<double>(order));
166 size_t* inverse =
new size_t[_confs_no];
168 for(
size_t ii = 0; ii < _confs_no; ii++)
169 inverse[indices[ii]] = ii;
173 reorder_array(_masses, inverse, _confs_no);
174 reorder_array(_probs, inverse, _confs_no, _confs ==
nullptr);
175 if(_confs !=
nullptr)
177 int* swapspace =
new int[allDim];
178 for(
size_t ii = 0; ii < _confs_no; ii++)
179 while(inverse[ii] != ii)
181 memcpy(swapspace, &_confs[ii*allDim], allDimSizeofInt);
182 memcpy(&_confs[ii*allDim], &_confs[inverse[ii]*allDim], allDimSizeofInt);
183 memcpy(&_confs[inverse[ii]*allDim], swapspace, allDimSizeofInt);
184 std::swap(inverse[inverse[ii]], inverse[ii]);
192 double FixedEnvelope::get_total_prob()
194 if(std::isnan(total_prob))
197 for(
size_t ii = 0; ii < _confs_no; ii++)
198 total_prob += _probs[ii];
203 void FixedEnvelope::scale(
double factor)
205 for(
size_t ii = 0; ii < _confs_no; ii++)
206 _probs[ii] *= factor;
207 total_prob *= factor;
210 void FixedEnvelope::normalize()
212 double tp = get_total_prob();
220 void FixedEnvelope::shift_mass(
double value)
222 for(
size_t ii = 0; ii < _confs_no; ii++)
223 _masses[ii] += value;
226 void FixedEnvelope::resample(
size_t samples,
double beta_bias)
229 throw std::logic_error(
"Resample called on an empty spectrum");
235 _probs[_confs_no-1] = (std::numeric_limits<double>::max)();
239 pprob += _probs[++pidx];
241 double covered_part = (pprob - cprob) / (1.0 - cprob);
242 while(samples * covered_part < beta_bias && samples > 0)
244 cprob += rdvariate_beta_1_b(samples) * (1.0 - cprob);
247 pprob += _probs[++pidx];
252 covered_part = (pprob - cprob) / (1.0 - cprob);
256 size_t nrtaken = rdvariate_binom(samples, covered_part);
257 _probs[pidx] +=
static_cast<double>(nrtaken);
263 memset(_probs + pidx, 0,
sizeof(
double)*(_confs_no - pidx));
266 FixedEnvelope FixedEnvelope::LinearCombination(
const std::vector<const FixedEnvelope*>& spectra,
const std::vector<double>& intensities)
268 return LinearCombination(spectra.data(), intensities.data(), spectra.size());
271 FixedEnvelope FixedEnvelope::LinearCombination(
const FixedEnvelope*
const * spectra,
const double* intensities,
size_t size)
274 for(
size_t ii = 0; ii < size; ii++)
275 ret_size += spectra[ii]->_confs_no;
277 double* newprobs =
reinterpret_cast<double*
>(malloc(
sizeof(
double)*ret_size));
278 if(newprobs ==
nullptr)
279 throw std::bad_alloc();
280 double* newmasses =
reinterpret_cast<double*
>(malloc(
sizeof(
double)*ret_size));
281 if(newmasses ==
nullptr)
284 throw std::bad_alloc();
288 for(
size_t ii = 0; ii < size; ii++)
290 double mul = intensities[ii];
291 for(
size_t jj = 0; jj < spectra[ii]->_confs_no; jj++)
292 newprobs[jj+cntr] = spectra[ii]->_probs[jj] * mul;
293 memcpy(newmasses + cntr, spectra[ii]->_masses,
sizeof(
double) * spectra[ii]->_confs_no);
294 cntr += spectra[ii]->_confs_no;
296 return FixedEnvelope(newmasses, newprobs, cntr);
299 double FixedEnvelope::WassersteinDistance(FixedEnvelope& other)
302 if((get_total_prob()*0.999 > other.get_total_prob()) || (other.get_total_prob() > get_total_prob()*1.001))
303 throw std::logic_error(
"Spectra must be normalized before computing Wasserstein Distance");
305 if(_confs_no == 0 || other._confs_no == 0)
309 other.sort_by_mass();
312 size_t idx_other = 0;
314 double acc_prob = 0.0;
315 double last_point = 0.0;
318 while(idx_this < _confs_no && idx_other < other._confs_no)
320 if(_masses[idx_this] < other._masses[idx_other])
322 ret += (_masses[idx_this] - last_point) * std::abs(acc_prob);
323 acc_prob += _probs[idx_this];
324 last_point = _masses[idx_this];
329 ret += (other._masses[idx_other] - last_point) * std::abs(acc_prob);
330 acc_prob -= other._probs[idx_other];
331 last_point = other._masses[idx_other];
336 acc_prob = std::abs(acc_prob);
338 while(idx_this < _confs_no)
340 ret += (_masses[idx_this] - last_point) * acc_prob;
341 acc_prob -= _probs[idx_this];
342 last_point = _masses[idx_this];
346 while(idx_other < other._confs_no)
348 ret += (other._masses[idx_other] - last_point) * acc_prob;
349 acc_prob -= other._probs[idx_other];
350 last_point = other._masses[idx_other];
358 double FixedEnvelope::OrientedWassersteinDistance(FixedEnvelope& other)
361 if((get_total_prob()*0.999 > other.get_total_prob()) || (other.get_total_prob() > get_total_prob()*1.001))
362 throw std::logic_error(
"Spectra must be normalized before computing Wasserstein Distance");
364 if(_confs_no == 0 || other._confs_no == 0)
368 other.sort_by_mass();
371 size_t idx_other = 0;
373 double acc_prob = 0.0;
374 double last_point = 0.0;
377 while(idx_this < _confs_no && idx_other < other._confs_no)
379 if(_masses[idx_this] < other._masses[idx_other])
381 ret += (_masses[idx_this] - last_point) * acc_prob;
382 acc_prob += _probs[idx_this];
383 last_point = _masses[idx_this];
388 ret += (other._masses[idx_other] - last_point) * acc_prob;
389 acc_prob -= other._probs[idx_other];
390 last_point = other._masses[idx_other];
395 while(idx_this < _confs_no)
397 ret += (_masses[idx_this] - last_point) * acc_prob;
398 acc_prob -= _probs[idx_this];
399 last_point = _masses[idx_this];
403 while(idx_other < other._confs_no)
405 ret += (other._masses[idx_other] - last_point) * acc_prob;
406 acc_prob -= other._probs[idx_other];
407 last_point = other._masses[idx_other];
414 double FixedEnvelope::AbyssalWassersteinDistance(FixedEnvelope& other,
double abyss_depth,
double other_scale)
417 other.sort_by_mass();
419 std::vector<std::pair<double, double>> carried;
422 size_t idx_other = 0;
426 auto finished = [&]() ->
bool {
return idx_this >= _confs_no && idx_other >= other._confs_no; };
427 auto next = [&]() -> std::pair<double, double> {
428 if(idx_this >= _confs_no || (idx_other < other._confs_no && _masses[idx_this] > other._masses[idx_other]))
430 std::pair<double, double> res = std::pair<double, double>(other._masses[idx_other], other._probs[idx_other]*other_scale);
436 std::pair<double, double> res = std::pair<double, double>(_masses[idx_this], -_probs[idx_this]);
442 double condemned = 0.0;
446 auto [m, p] = next();
447 if(!carried.empty() && carried[0].second * p > 0.0)
449 carried.emplace_back(m, p);
453 while(!carried.empty())
455 auto& [cm, cp] = carried.back();
456 if(m - cm >= abyss_depth)
458 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
459 condemned += fabs(it->second);
465 accd += fabs((m-cm)*cp);
471 accd += fabs((m-cm)*p);
478 carried.emplace_back(m, p);
482 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
483 condemned += fabs(it->second);
485 return accd + condemned * abyss_depth * 0.5;
489 double FixedEnvelope::ScaledAbyssalWassersteinDistance(FixedEnvelope *
const * others,
double abyss_depth,
const double* other_scales,
const size_t N)
493 std::priority_queue<std::pair<double, size_t>> PQ;
494 std::unique_ptr<size_t[]> indices = std::make_unique<size_t[]>(N);
495 memset(indices.get(), 0,
sizeof(
size_t)*N);
497 for(
size_t ii=0; ii<N; ii++)
499 others[ii]->sort_by_mass();
500 if(others[ii]->_confs_no > 0)
501 PQ.emplace_back({others._probs[0], ii});
507 std::vector<std::pair<double, double>> carried;
510 size_t idx_other = 0;
514 auto finished = [&]() ->
bool {
return idx_this >= _confs_no && PQ.empty(); };
515 auto next = [&]() -> std::tuple<double, double, size_t> {
516 if(idx_this >= _confs_no || (idx_other < other._confs_no && _masses[idx_this] > other._masses[idx_other]))
518 std::pair<double, double> res = std::pair<double, double>(other._masses[idx_other], other._probs[idx_other]*other_scale);
524 std::pair<double, double> res = std::pair<double, double>(_masses[idx_this], -_probs[idx_this]);
530 double condemned = 0.0;
534 auto [m, p] = next();
535 if(!carried.empty() && carried[0].second * p > 0.0)
537 carried.emplace_back(m, p);
541 while(!carried.empty())
543 auto& [cm, cp] = carried.back();
544 if(m - cm >= abyss_depth)
546 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
547 condemned += fabs(it->second);
553 accd += fabs((m-cm)*cp);
559 accd += fabs((m-cm)*p);
566 carried.emplace_back(m, p);
570 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
571 condemned += fabs(it->second);
573 return accd + condemned * abyss_depth * 0.5;
579 double AbyssalWassersteinDistanceGrad(FixedEnvelope*
const* envelopes,
const double* scales,
double* ret_gradient,
size_t N,
double abyss_depth_exp,
double abyss_depth_the)
582 std::unique_ptr<size_t[]> env_idx = std::make_unique<size_t[]>(N+1);
583 memset(env_idx.get(), 0, (N+1)*
sizeof(
size_t));
584 memset(ret_gradient, 0, (N+1)*
sizeof(
double));
586 auto pq_cmp = [](std::pair<double, size_t>& p1, std::pair<double, size_t>& p2) {
return p1.first > p2.first; };
587 std::priority_queue<std::pair<double, size_t>, std::vector<std::pair<double, size_t>>, decltype(pq_cmp)> PQ(pq_cmp);
589 for(
size_t ii=0; ii<=N; ii++)
591 envelopes[ii]->sort_by_mass();
592 if(envelopes[ii]->_confs_no > 0)
594 std::cout <<
"Initial push: " << envelopes[ii]->_masses[0] <<
" " << ii <<
'\n';
595 PQ.push({envelopes[ii]->_masses[0], ii});
599 std::vector<std::tuple<double, double, size_t>> carried;
601 auto next = [&]() -> std::tuple<double, double, size_t> {
602 auto [mass, eidx] = PQ.top();
603 double prob = envelopes[eidx]->_probs[env_idx[eidx]];
608 prob = prob * scales[eidx];
609 std::cout <<
"Next: " << mass <<
" " << prob <<
" " << eidx <<
'\n';
611 if(env_idx[eidx] < envelopes[eidx]->_confs_no)
612 PQ.push({envelopes[eidx]->_masses[env_idx[eidx]], eidx});
614 return {mass, prob, eidx};
617 double condemned = 0.0;
619 const double max_flow_dist = abyss_depth_exp + abyss_depth_the;
620 max_flow_dist *= 2.0;
624 auto [m, p, eidx] = next();
625 if(!carried.empty() && std::get<1>(carried[0]) * p > 0.0)
627 carried.emplace_back(m, p, eidx);
631 while(!carried.empty())
633 auto& [cm, cp, ceidx] = carried.back();
634 if(m - cm >= max_flow_dist)
636 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
637 condemned += fabs(std::get<1>(*it));
641 std::cout <<
"accing\n";
644 accd += fabs((m-cm)*cp);
647 std::cout <<
"cprob:" << cp <<
'\n';
651 accd += fabs((m-cm)*p);
653 std::cout <<
"prob:" << p <<
'\n';
659 carried.emplace_back(m, p, eidx);
663 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
664 condemned += fabs(std::get<1>(*it));
666 std::cout <<
"ISO:" << accd <<
" " << condemned <<
'\n';
667 return accd + condemned * max_flow_dist * 0.5;
670 auto [m, p, eidx] = next();
671 if(!carried.empty() && (std::get<1>(carried[0]) * p > 0.0))
673 carried.emplace_back(m, p, eidx);
677 while(!carried.empty())
679 auto& [cm, cp, ceidx] = carried.back();
680 if(m - cm >= max_flow_dist)
682 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
684 flow = fabs(std::get<1>(*it));
685 const size_t target_idx = std::get<2>(*it);
686 flow *= target_idx == 0 ? abyss_depth_exp : abyss_depth_the;
687 ret_gradient[target_idx] += flow;
689 std::cout <<
"condemnin': " << m <<
" " << p <<
" " << eidx <<
" | " << cm <<
" " << cp <<
" " << ceidx <<
'\n';
696 flow = fabs((m-cm)*cp);
699 ret_gradient[ceidx] += flow;
704 flow = fabs((m-cm)*p);
707 ret_gradient[eidx] += flow;
713 carried.emplace_back(m, p, eidx);
717 for(
auto it = carried.cbegin(); it != carried.cend(); it++)
718 condemned += fabs(std::get<1>(*it));
721 return accd + condemned * (abyss_depth_exp + abyss_depth_the) * 0.5;
726 std::tuple<double, double, double> FixedEnvelope::WassersteinMatch(FixedEnvelope& other,
double flow_distance,
double other_scale)
729 return {0.0, other.get_total_prob() * other_scale, 0.0};
731 double unmatched1 = 0.0;
732 double unmatched2 = 0.0;
733 double massflow = 0.0;
736 other.sort_by_mass();
739 size_t idx_other = 0;
740 double used_prob_this = 0.0;
741 double used_prob_other = 0.0;
743 while(idx_this < _confs_no && idx_other < other._confs_no)
746 while(moved && idx_this < _confs_no && idx_other < other._confs_no)
749 if(_masses[idx_this] < other._masses[idx_other] - flow_distance)
751 unmatched1 += _probs[idx_this] - used_prob_this;
752 used_prob_this = 0.0;
756 if(other._masses[idx_other] < _masses[idx_this] - flow_distance)
758 unmatched2 += other._probs[idx_other]*other_scale - used_prob_other;
759 used_prob_other = 0.0;
764 if(idx_this < _confs_no && idx_other < other._confs_no)
766 assert(_probs[idx_this] - used_prob_this >= 0.0);
767 assert(other._probs[idx_other]*other_scale - used_prob_other >= 0.0);
769 if(_probs[idx_this] - used_prob_this < other._probs[idx_other]*other_scale - used_prob_other)
771 massflow += _probs[idx_this] - used_prob_this;
772 used_prob_other += _probs[idx_this] - used_prob_this;
773 assert(used_prob_other >= 0.0);
774 used_prob_this = 0.0;
779 massflow += other._probs[idx_other]*other_scale - used_prob_other;
780 used_prob_this += other._probs[idx_other]*other_scale - used_prob_other;
781 assert(used_prob_this >= 0.0);
782 used_prob_other = 0.0;
788 unmatched1 -= used_prob_this;
789 unmatched2 -= used_prob_other;
791 for(; idx_this < _confs_no; idx_this++)
792 unmatched1 += _probs[idx_this];
793 for(; idx_other < other._confs_no; idx_other++)
794 unmatched2 += other._probs[idx_other]*other_scale;
796 return {unmatched1, unmatched2, massflow};
799 FixedEnvelope FixedEnvelope::bin(
double bin_width,
double middle)
808 ret.reallocate_memory<
false>(ISOSPEC_INIT_TABLE_SIZE);
812 double curr_mass = _masses[0];
813 double accd_prob = _probs[0];
814 for(
size_t ii = 1; ii<_confs_no; ii++)
816 if(curr_mass != _masses[ii])
818 ret.store_conf(curr_mass, accd_prob);
819 curr_mass = _masses[ii];
820 accd_prob = _probs[ii];
823 accd_prob += _probs[ii];
825 ret.store_conf(curr_mass, accd_prob);
831 double half_width = 0.5*bin_width;
832 double hwmm = half_width-middle;
834 while(ii < _confs_no)
836 double current_bin_middle = floor((_masses[ii]+hwmm)/bin_width)*bin_width + middle;
837 double current_bin_end = current_bin_middle + half_width;
838 double bin_prob = 0.0;
840 while(ii < _confs_no && _masses[ii] <= current_bin_end)
842 bin_prob += _probs[ii];
845 ret.store_conf(current_bin_middle, bin_prob);
851 template<
bool tgetConfs>
void FixedEnvelope::reallocate_memory(
size_t new_size)
853 current_size = new_size;
855 _masses =
reinterpret_cast<double*
>(realloc(_masses, new_size *
sizeof(
double)));
856 if(_masses ==
nullptr)
857 throw std::bad_alloc();
858 tmasses = _masses + _confs_no;
860 _probs =
reinterpret_cast<double*
>(realloc(_probs, new_size *
sizeof(
double)));
861 if(_probs ==
nullptr)
862 throw std::bad_alloc();
863 tprobs = _probs + _confs_no;
865 constexpr_if(tgetConfs)
867 _confs =
reinterpret_cast<int*
>(realloc(_confs, new_size * allDimSizeofInt));
868 if(_confs ==
nullptr)
869 throw std::bad_alloc();
870 tconfs = _confs + (allDim * _confs_no);
874 void FixedEnvelope::slow_reallocate_memory(
size_t new_size)
876 current_size = new_size;
878 _masses =
reinterpret_cast<double*
>(realloc(_masses, new_size *
sizeof(
double)));
879 if(_masses ==
nullptr)
880 throw std::bad_alloc();
881 tmasses = _masses + _confs_no;
883 _probs =
reinterpret_cast<double*
>(realloc(_probs, new_size *
sizeof(
double)));
884 if(_probs ==
nullptr)
885 throw std::bad_alloc();
886 tprobs = _probs + _confs_no;
888 if(_confs !=
nullptr)
890 _confs =
reinterpret_cast<int*
>(realloc(_confs, new_size * allDimSizeofInt));
891 if(_confs ==
nullptr)
892 throw std::bad_alloc();
893 tconfs = _confs + (allDim * _confs_no);
897 template<
bool tgetConfs>
void FixedEnvelope::threshold_init(Iso&& iso,
double threshold,
bool absolute)
899 IsoThresholdGenerator generator(std::move(iso), threshold, absolute);
901 size_t tab_size = generator.count_confs();
902 this->allDim = generator.getAllDim();
903 this->allDimSizeofInt = this->allDim *
sizeof(int);
905 this->reallocate_memory<tgetConfs>(tab_size);
907 double* ttmasses = this->_masses;
908 double* ttprobs = this->_probs;
909 ISOSPEC_MAYBE_UNUSED
int* ttconfs;
910 constexpr_if(tgetConfs)
913 while(generator.advanceToNextConfiguration())
915 *ttmasses = generator.mass(); ttmasses++;
916 *ttprobs = generator.prob(); ttprobs++;
917 constexpr_if(tgetConfs) { generator.get_conf_signature(ttconfs); ttconfs += allDim; }
920 this->_confs_no = tab_size;
923 template void FixedEnvelope::threshold_init<true>(Iso&& iso,
double threshold,
bool absolute);
924 template void FixedEnvelope::threshold_init<false>(Iso&& iso,
double threshold,
bool absolute);
927 template<
bool tgetConfs>
void FixedEnvelope::total_prob_init(Iso&& iso,
double target_total_prob,
bool optimize)
929 if(target_total_prob <= 0.0)
932 if(target_total_prob >= 1.0)
934 threshold_init<tgetConfs>(std::move(iso), 0.0,
true);
938 IsoLayeredGenerator generator(std::move(iso), 1000, 1000,
true, std::min<double>(target_total_prob, 0.9999));
940 this->allDim = generator.getAllDim();
941 this->allDimSizeofInt = this->allDim*
sizeof(int);
944 this->reallocate_memory<tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);
946 size_t last_switch = 0;
947 double prob_at_last_switch = 0.0;
948 double prob_so_far = 0.0;
951 const double sum_above = log1p(-target_total_prob) - 2.3025850929940455;
956 while(generator.advanceToNextConfigurationWithinLayer())
958 this->
template addConfILG<tgetConfs>(generator);
959 prob_so_far += *(tprobs-1);
960 if(prob_so_far >= target_total_prob)
964 while(generator.advanceToNextConfigurationWithinLayer())
965 this->
template addConfILG<tgetConfs>(generator);
972 if(prob_so_far >= target_total_prob)
975 last_switch = this->_confs_no;
976 prob_at_last_switch = prob_so_far;
978 layer_delta = sum_above - log1p(-prob_so_far);
979 layer_delta = (std::max)((std::min)(layer_delta, -0.1), -5.0);
980 }
while(generator.nextLayer(layer_delta));
982 if(!optimize || prob_so_far <= target_total_prob)
991 int* conf_swapspace =
nullptr;
992 constexpr_if(tgetConfs)
993 conf_swapspace =
reinterpret_cast<int*
>(malloc(this->allDimSizeofInt));
995 size_t start = last_switch;
996 size_t end = this->_confs_no;
997 double sum_to_start = prob_at_last_switch;
1002 size_t len = end - start;
1003 #if ISOSPEC_BUILDING_R
1004 size_t pivot = len/2 + start;
1006 size_t pivot = random_gen() % len + start;
1010 double pprob = this->_probs[pivot];
1011 swap<tgetConfs>(pivot, end-1, conf_swapspace);
1013 double new_csum = sum_to_start;
1015 size_t loweridx = start;
1016 for(
size_t ii = start; ii < end-1; ii++)
1017 if(this->_probs[ii] > pprob)
1019 swap<tgetConfs>(ii, loweridx, conf_swapspace);
1020 new_csum += this->_probs[loweridx];
1024 swap<tgetConfs>(end-1, loweridx, conf_swapspace);
1027 if(new_csum < target_total_prob)
1029 start = loweridx + 1;
1030 sum_to_start = new_csum + this->_probs[loweridx];
1036 constexpr_if(tgetConfs)
1037 free(conf_swapspace);
1039 if(end <= current_size/2)
1041 this->
template reallocate_memory<tgetConfs>(end);
1043 this->_confs_no = end;
1046 template void FixedEnvelope::total_prob_init<true>(Iso&& iso,
double target_total_prob,
bool optimize);
1047 template void FixedEnvelope::total_prob_init<false>(Iso&& iso,
double target_total_prob,
bool optimize);
1049 template<
bool tgetConfs>
void FixedEnvelope::stochastic_init(Iso&& iso,
size_t _no_molecules,
double _precision,
double _beta_bias)
1051 IsoStochasticGenerator generator(std::move(iso), _no_molecules, _precision, _beta_bias);
1053 this->allDim = generator.getAllDim();
1054 this->allDimSizeofInt = this->allDim *
sizeof(int);
1056 this->reallocate_memory<tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);
1058 while(generator.advanceToNextConfiguration())
1059 addConfILG<tgetConfs, IsoStochasticGenerator>(generator);
1062 template void FixedEnvelope::stochastic_init<true>(Iso&& iso,
size_t _no_molecules,
double _precision,
double _beta_bias);
1063 template void FixedEnvelope::stochastic_init<false>(Iso&& iso,
size_t _no_molecules,
double _precision,
double _beta_bias);
1065 double FixedEnvelope::empiric_average_mass()
1068 for(
size_t ii = 0; ii < _confs_no; ii++)
1070 ret += _masses[ii] * _probs[ii];
1072 return ret / get_total_prob();
1075 double FixedEnvelope::empiric_variance()
1078 double avg = empiric_average_mass();
1079 for(
size_t ii = 0; ii < _confs_no; ii++)
1081 double msq = _masses[ii] - avg;
1082 ret += msq * msq * _probs[ii];
1085 return ret / get_total_prob();
1088 FixedEnvelope FixedEnvelope::Binned(Iso&& iso,
double target_total_prob,
double bin_width,
double bin_middle)
1092 double min_mass = iso.getLightestPeakMass();
1093 double range_len = iso.getHeaviestPeakMass() - min_mass;
1094 size_t no_bins =
static_cast<size_t>(range_len / bin_width) + 2;
1095 double half_width = 0.5*bin_width;
1096 double hwmm = half_width-bin_middle;
1097 size_t idx_min = floor((min_mass + hwmm) / bin_width);
1098 size_t idx_max = idx_min + no_bins;
1101 # if ISOSPEC_GOT_MMAN
1102 acc =
reinterpret_cast<double*
>(mmap(
nullptr,
sizeof(
double)*no_bins, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
1105 acc =
reinterpret_cast<double*
>(calloc(no_bins,
sizeof(
double)));
1108 throw std::bad_alloc();
1112 IsoLayeredGenerator ITG(std::move(iso));
1116 while((non_empty = ITG.advanceToNextConfiguration()) && ITG.prob() == 0.0)
1121 double accum_prob = ITG.prob();
1122 size_t nonzero_idx = floor((ITG.mass() + hwmm)/bin_width);
1123 acc[nonzero_idx] = accum_prob;
1125 while(ITG.advanceToNextConfiguration() && accum_prob < target_total_prob)
1127 double prob = ITG.prob();
1129 size_t bin_idx = floor((ITG.mass() + hwmm)/bin_width);
1130 acc[bin_idx] += prob;
1135 size_t distance_10da =
static_cast<size_t>(10.0/bin_width) + 1;
1137 size_t empty_steps = 0;
1139 ret.reallocate_memory<
false>(ISOSPEC_INIT_TABLE_SIZE);
1141 for(
size_t ii = nonzero_idx; ii >= idx_min && empty_steps < distance_10da; ii--)
1146 ret.store_conf(
static_cast<double>(ii)*bin_width + bin_middle, acc[ii]);
1153 for(
size_t ii = nonzero_idx+1; ii < idx_max && empty_steps < distance_10da; ii++)
1158 ret.store_conf(
static_cast<double>(ii)*bin_width + bin_middle, acc[ii]);
1167 # if ISOSPEC_GOT_MMAN
1168 munmap(acc,
sizeof(
double)*no_bins);