Loading...
Searching...
No Matches
fixedEnvelopes.cpp
/*
* Copyright (C) 2015-2020 Mateusz Łącki and Michał Startek.
*
* This file is part of IsoSpec.
*
* IsoSpec is free software: you can redistribute it and/or modify
* it under the terms of the Simplified ("2-clause") BSD licence.
*
* IsoSpec is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*
* You should have received a copy of the Simplified BSD Licence
* along with IsoSpec. If not, see <https://opensource.org/licenses/BSD-2-Clause>.
*/
#include "fixedEnvelopes.h"
#include <limits>
#include <cassert>
#include "isoMath.h"
namespace IsoSpec
{
FixedEnvelope::FixedEnvelope(const FixedEnvelope& other) :
_masses(array_copy<double>(other._masses, other._confs_no)),
_probs(array_copy<double>(other._probs, other._confs_no)),
_confs(array_copy_nptr<int>(other._confs, other._confs_no*other.allDim)),
_confs_no(other._confs_no),
allDim(other.allDim),
sorted_by_mass(other.sorted_by_mass),
sorted_by_prob(other.sorted_by_prob),
total_prob(other.total_prob)
{}
FixedEnvelope::FixedEnvelope(FixedEnvelope&& other) :
_masses(other._masses),
_probs(other._probs),
_confs(other._confs),
_confs_no(other._confs_no),
allDim(other.allDim),
sorted_by_mass(other.sorted_by_mass),
sorted_by_prob(other.sorted_by_prob),
total_prob(other.total_prob)
{
other._masses = nullptr;
other._probs = nullptr;
other._confs = nullptr;
other._confs_no = 0;
other.total_prob = 0.0;
}
FixedEnvelope::FixedEnvelope(double* in_masses, double* in_probs, size_t in_confs_no, bool masses_sorted, bool probs_sorted, double _total_prob) :
_masses(in_masses),
_probs(in_probs),
_confs(nullptr),
_confs_no(in_confs_no),
allDim(0),
sorted_by_mass(masses_sorted),
sorted_by_prob(probs_sorted),
total_prob(_total_prob)
{}
FixedEnvelope FixedEnvelope::operator+(const FixedEnvelope& other) const
{
double* nprobs = reinterpret_cast<double*>(malloc(sizeof(double) * (_confs_no+other._confs_no)));
if(nprobs == nullptr)
throw std::bad_alloc();
double* nmasses = reinterpret_cast<double*>(malloc(sizeof(double) * (_confs_no+other._confs_no)));
if(nmasses == nullptr)
{
free(nprobs);
throw std::bad_alloc();
}
memcpy(nprobs, _probs, sizeof(double) * _confs_no);
memcpy(nmasses, _masses, sizeof(double) * _confs_no);
memcpy(nprobs+_confs_no, other._probs, sizeof(double) * other._confs_no);
memcpy(nmasses+_confs_no, other._masses, sizeof(double) * other._confs_no);
return FixedEnvelope(nmasses, nprobs, _confs_no + other._confs_no);
}
FixedEnvelope FixedEnvelope::operator*(const FixedEnvelope& other) const
{
double* nprobs = reinterpret_cast<double*>(malloc(sizeof(double) * _confs_no * other._confs_no));
if(nprobs == nullptr)
throw std::bad_alloc();
// deepcode ignore CMemoryLeak: It's not a memleak: the memory is passed to FixedEnvelope which
// deepcode ignore CMemoryLeak: takes ownership of it, and will properly free() it in destructor.
double* nmasses = reinterpret_cast<double*>(malloc(sizeof(double) * _confs_no * other._confs_no));
if(nmasses == nullptr)
{
free(nprobs);
throw std::bad_alloc();
}
size_t tgt_idx = 0;
for(size_t ii = 0; ii < _confs_no; ii++)
for(size_t jj = 0; jj < other._confs_no; jj++)
{
nprobs[tgt_idx] = _probs[ii] * other._probs[jj];
nmasses[tgt_idx] = _masses[ii] + other._masses[jj];
tgt_idx++;
}
return FixedEnvelope(nmasses, nprobs, tgt_idx);
}
void FixedEnvelope::sort_by_mass()
{
if(sorted_by_mass)
return;
sort_by(_masses);
sorted_by_mass = true;
sorted_by_prob = false;
}
void FixedEnvelope::sort_by_prob()
{
if(sorted_by_prob)
return;
sort_by(_probs);
sorted_by_prob = true;
sorted_by_mass = false;
}
template<typename T> void reorder_array(T* arr, size_t* order, size_t size, bool can_destroy = false)
{
if(!can_destroy)
{
size_t* order_c = new size_t[size];
memcpy(order_c, order, sizeof(size_t)*size);
order = order_c;
}
for(size_t ii = 0; ii < size; ii++)
while(order[ii] != ii)
{
std::swap(arr[ii], arr[order[ii]]);
std::swap(order[order[ii]], order[ii]);
}
if(!can_destroy)
delete[] order;
}
void FixedEnvelope::sort_by(double* order)
{
if(_confs_no <= 1)
return;
size_t* indices = new size_t[_confs_no];
for(size_t ii = 0; ii < _confs_no; ii++)
indices[ii] = ii;
std::sort<size_t*>(indices, indices + _confs_no, TableOrder<double>(order));
size_t* inverse = new size_t[_confs_no];
for(size_t ii = 0; ii < _confs_no; ii++)
inverse[indices[ii]] = ii;
delete[] indices;
reorder_array(_masses, inverse, _confs_no);
reorder_array(_probs, inverse, _confs_no, _confs == nullptr);
if(_confs != nullptr)
{
int* swapspace = new int[allDim];
for(size_t ii = 0; ii < _confs_no; ii++)
while(inverse[ii] != ii)
{
memcpy(swapspace, &_confs[ii*allDim], allDimSizeofInt);
memcpy(&_confs[ii*allDim], &_confs[inverse[ii]*allDim], allDimSizeofInt);
memcpy(&_confs[inverse[ii]*allDim], swapspace, allDimSizeofInt);
std::swap(inverse[inverse[ii]], inverse[ii]);
}
delete[] swapspace;
}
delete[] inverse;
}
double FixedEnvelope::get_total_prob()
{
if(std::isnan(total_prob))
{
total_prob = 0.0;
for(size_t ii = 0; ii < _confs_no; ii++)
total_prob += _probs[ii];
}
return total_prob;
}
void FixedEnvelope::scale(double factor)
{
for(size_t ii = 0; ii < _confs_no; ii++)
_probs[ii] *= factor;
total_prob *= factor;
}
void FixedEnvelope::normalize()
{
double tp = get_total_prob();
if(tp != 1.0)
{
scale(1.0/tp);
total_prob = 1.0;
}
}
void FixedEnvelope::shift_mass(double value)
{
for(size_t ii = 0; ii < _confs_no; ii++)
_masses[ii] += value;
}
void FixedEnvelope::resample(size_t samples, double beta_bias)
{
if(_confs_no == 0)
throw std::logic_error("Resample called on an empty spectrum");
double pprob = 0.0;
double cprob = 0.0;
size_t pidx = -1; // Overflows - but it doesn't matter.
_probs[_confs_no-1] = (std::numeric_limits<double>::max)();
while(samples > 0)
{
pprob += _probs[++pidx];
_probs[pidx] = 0.0;
double covered_part = (pprob - cprob) / (1.0 - cprob);
while(samples * covered_part < beta_bias && samples > 0)
{
cprob += rdvariate_beta_1_b(samples) * (1.0 - cprob);
while(pprob < cprob)
{
pprob += _probs[++pidx];
_probs[pidx] = 0.0;
}
_probs[pidx] += 1.0;
samples--;
covered_part = (pprob - cprob) / (1.0 - cprob);
}
if(samples <= 0)
break;
size_t nrtaken = rdvariate_binom(samples, covered_part);
_probs[pidx] += static_cast<double>(nrtaken);
samples -= nrtaken;
cprob = pprob;
}
pidx++;
memset(_probs + pidx, 0, sizeof(double)*(_confs_no - pidx));
}
FixedEnvelope FixedEnvelope::LinearCombination(const std::vector<const FixedEnvelope*>& spectra, const std::vector<double>& intensities)
{
return LinearCombination(spectra.data(), intensities.data(), spectra.size());
}
FixedEnvelope FixedEnvelope::LinearCombination(const FixedEnvelope* const * spectra, const double* intensities, size_t size)
{
size_t ret_size = 0;
for(size_t ii = 0; ii < size; ii++)
ret_size += spectra[ii]->_confs_no;
double* newprobs = reinterpret_cast<double*>(malloc(sizeof(double)*ret_size));
if(newprobs == nullptr)
throw std::bad_alloc();
double* newmasses = reinterpret_cast<double*>(malloc(sizeof(double)*ret_size));
if(newmasses == nullptr)
{
free(newprobs);
throw std::bad_alloc();
}
size_t cntr = 0;
for(size_t ii = 0; ii < size; ii++)
{
double mul = intensities[ii];
for(size_t jj = 0; jj < spectra[ii]->_confs_no; jj++)
newprobs[jj+cntr] = spectra[ii]->_probs[jj] * mul;
memcpy(newmasses + cntr, spectra[ii]->_masses, sizeof(double) * spectra[ii]->_confs_no);
cntr += spectra[ii]->_confs_no;
}
return FixedEnvelope(newmasses, newprobs, cntr);
}
double FixedEnvelope::WassersteinDistance(FixedEnvelope& other)
{
double ret = 0.0;
if((get_total_prob()*0.999 > other.get_total_prob()) || (other.get_total_prob() > get_total_prob()*1.001))
throw std::logic_error("Spectra must be normalized before computing Wasserstein Distance");
if(_confs_no == 0 || other._confs_no == 0)
return 0.0;
sort_by_mass();
other.sort_by_mass();
size_t idx_this = 0;
size_t idx_other = 0;
double acc_prob = 0.0;
double last_point = 0.0;
while(idx_this < _confs_no && idx_other < other._confs_no)
{
if(_masses[idx_this] < other._masses[idx_other])
{
ret += (_masses[idx_this] - last_point) * std::abs(acc_prob);
acc_prob += _probs[idx_this];
last_point = _masses[idx_this];
idx_this++;
}
else
{
ret += (other._masses[idx_other] - last_point) * std::abs(acc_prob);
acc_prob -= other._probs[idx_other];
last_point = other._masses[idx_other];
idx_other++;
}
}
acc_prob = std::abs(acc_prob);
while(idx_this < _confs_no)
{
ret += (_masses[idx_this] - last_point) * acc_prob;
acc_prob -= _probs[idx_this];
last_point = _masses[idx_this];
idx_this++;
}
while(idx_other < other._confs_no)
{
ret += (other._masses[idx_other] - last_point) * acc_prob;
acc_prob -= other._probs[idx_other];
last_point = other._masses[idx_other];
idx_other++;
}
return ret;
}
double FixedEnvelope::OrientedWassersteinDistance(FixedEnvelope& other)
{
double ret = 0.0;
if((get_total_prob()*0.999 > other.get_total_prob()) || (other.get_total_prob() > get_total_prob()*1.001))
throw std::logic_error("Spectra must be normalized before computing Wasserstein Distance");
if(_confs_no == 0 || other._confs_no == 0)
return 0.0;
sort_by_mass();
other.sort_by_mass();
size_t idx_this = 0;
size_t idx_other = 0;
double acc_prob = 0.0;
double last_point = 0.0;
while(idx_this < _confs_no && idx_other < other._confs_no)
{
if(_masses[idx_this] < other._masses[idx_other])
{
ret += (_masses[idx_this] - last_point) * acc_prob;
acc_prob += _probs[idx_this];
last_point = _masses[idx_this];
idx_this++;
}
else
{
ret += (other._masses[idx_other] - last_point) * acc_prob;
acc_prob -= other._probs[idx_other];
last_point = other._masses[idx_other];
idx_other++;
}
}
while(idx_this < _confs_no)
{
ret += (_masses[idx_this] - last_point) * acc_prob;
acc_prob -= _probs[idx_this];
last_point = _masses[idx_this];
idx_this++;
}
while(idx_other < other._confs_no)
{
ret += (other._masses[idx_other] - last_point) * acc_prob;
acc_prob -= other._probs[idx_other];
last_point = other._masses[idx_other];
idx_other++;
}
return ret;
}
double FixedEnvelope::AbyssalWassersteinDistance(FixedEnvelope& other, double abyss_depth, double other_scale)
{
sort_by_mass();
other.sort_by_mass();
std::vector<std::pair<double, double>> carried;
size_t idx_this = 0;
size_t idx_other = 0;
//std::cout << "AAA" << std::endl;
auto finished = [&]() -> bool { return idx_this >= _confs_no && idx_other >= other._confs_no; };
auto next = [&]() -> std::pair<double, double> {
if(idx_this >= _confs_no || (idx_other < other._confs_no && _masses[idx_this] > other._masses[idx_other]))
{
std::pair<double, double> res = std::pair<double, double>(other._masses[idx_other], other._probs[idx_other]*other_scale);
idx_other++;
return res;
}
else
{
std::pair<double, double> res = std::pair<double, double>(_masses[idx_this], -_probs[idx_this]);
idx_this++;
return res;
}
};
double accd = 0.0;
double condemned = 0.0;
while(!finished())
{
auto pair = next();
double m = pair.first;
double p = pair.second;
if(!carried.empty() && carried[0].second * p > 0.0)
{
carried.emplace_back(m, p);
continue;
}
while(!carried.empty())
{
double cm = carried.back().first;
double cp = carried.back().second;
if(m - cm >= abyss_depth)
{
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(it->second);
carried.clear();
break;
}
if((cp+p)*p > 0.0)
{
accd += fabs((m-cm)*cp);
p += cp;
carried.pop_back();
}
else
{
accd += fabs((m-cm)*p);
cp += p;
p = 0.0;
break;
}
}
if(p != 0.0)
carried.emplace_back(m, p);
//std::cout << m << " " << p << std::endl;
}
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(it->second);
return accd + condemned * abyss_depth * 0.5;
}
#if 0
double FixedEnvelope::ScaledAbyssalWassersteinDistance(FixedEnvelope * const * others, double abyss_depth, const double* other_scales, const size_t N)
{
sort_by_mass();
std::priority_queue<std::pair<double, size_t>> PQ;
std::unique_ptr<size_t[]> indices = std::make_unique<size_t[]>(N);
memset(indices.get(), 0, sizeof(size_t)*N);
for(size_t ii=0; ii<N; ii++)
{
others[ii]->sort_by_mass();
if(others[ii]->_confs_no > 0)
PQ.emplace_back({others._probs[0], ii});
}
std::vector<std::pair<double, double>> carried;
size_t idx_this = 0;
size_t idx_other = 0;
//std::cout << "AAA" << std::endl;
auto finished = [&]() -> bool { return idx_this >= _confs_no && PQ.empty(); };
auto next = [&]() -> std::tuple<double, double, size_t> {
if(idx_this >= _confs_no || (idx_other < other._confs_no && _masses[idx_this] > other._masses[idx_other]))
{
std::pair<double, double> res = std::pair<double, double>(other._masses[idx_other], other._probs[idx_other]*other_scale);
idx_other++;
return res;
}
else
{
std::pair<double, double> res = std::pair<double, double>(_masses[idx_this], -_probs[idx_this]);
idx_this++;
return res;
}
};
double accd = 0.0;
double condemned = 0.0;
while(!finished())
{
auto [m, p] = next();
if(!carried.empty() && carried[0].second * p > 0.0)
{
carried.emplace_back(m, p);
continue;
}
while(!carried.empty())
{
auto& [cm, cp] = carried.back();
if(m - cm >= abyss_depth)
{
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(it->second);
carried.clear();
break;
}
if((cp+p)*p > 0.0)
{
accd += fabs((m-cm)*cp);
p += cp;
carried.pop_back();
}
else
{
accd += fabs((m-cm)*p);
cp += p;
p = 0.0;
break;
}
}
if(p != 0.0)
carried.emplace_back(m, p);
//std::cout << m << " " << p << std::endl;
}
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(it->second);
return accd + condemned * abyss_depth * 0.5;
}
#endif
#if 0
double AbyssalWassersteinDistanceGrad(FixedEnvelope* const* envelopes, const double* scales, double* ret_gradient, size_t N, double abyss_depth_exp, double abyss_depth_the)
{
return 0.0;
std::unique_ptr<size_t[]> env_idx = std::make_unique<size_t[]>(N+1);
memset(env_idx.get(), 0, (N+1)*sizeof(size_t));
memset(ret_gradient, 0, (N+1)*sizeof(double));
auto pq_cmp = [](std::pair<double, size_t>& p1, std::pair<double, size_t>& p2) { return p1.first > p2.first; };
std::priority_queue<std::pair<double, size_t>, std::vector<std::pair<double, size_t>>, decltype(pq_cmp)> PQ(pq_cmp);
for(size_t ii=0; ii<=N; ii++)
{
envelopes[ii]->sort_by_mass();
if(envelopes[ii]->_confs_no > 0)
{
std::cout << "Initial push: " << envelopes[ii]->_masses[0] << " " << ii << '\n';
PQ.push({envelopes[ii]->_masses[0], ii});
}
}
std::vector<std::tuple<double, double, size_t>> carried;
auto next = [&]() -> std::tuple<double, double, size_t> {
auto [mass, eidx] = PQ.top();
double prob = envelopes[eidx]->_probs[env_idx[eidx]];
PQ.pop();
if(eidx == 0)
prob = -prob;
else
prob = prob * scales[eidx];
std::cout << "Next: " << mass << " " << prob << " " << eidx << '\n';
env_idx[eidx]++;
if(env_idx[eidx] < envelopes[eidx]->_confs_no)
PQ.push({envelopes[eidx]->_masses[env_idx[eidx]], eidx});
return {mass, prob, eidx};
};
double accd = 0.0;
double condemned = 0.0;
//double flow;
const double max_flow_dist = abyss_depth_exp + abyss_depth_the;
max_flow_dist *= 2.0;
while(!PQ.empty())
{
auto [m, p, eidx] = next();
if(!carried.empty() && std::get<1>(carried[0]) * p > 0.0)
{
carried.emplace_back(m, p, eidx);
continue;
}
while(!carried.empty())
{
auto& [cm, cp, ceidx] = carried.back();
if(m - cm >= max_flow_dist)
{
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(std::get<1>(*it));
carried.clear();
break;
}
std::cout << "accing\n";
if((cp+p)*p > 0.0)
{
accd += fabs((m-cm)*cp);
p += cp;
carried.pop_back();
std::cout << "cprob:" << cp << '\n';
}
else
{
accd += fabs((m-cm)*p);
cp += p;
std::cout << "prob:" << p << '\n';
p = 0.0;
break;
}
}
if(p != 0.0)
carried.emplace_back(m, p, eidx);
//std::cout << m << " " << p << std::endl;
}
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(std::get<1>(*it));
std::cout << "ISO:" << accd << " " << condemned << '\n';
return accd + condemned * max_flow_dist * 0.5;
while(!PQ.empty())
{
auto [m, p, eidx] = next();
if(!carried.empty() && (std::get<1>(carried[0]) * p > 0.0))
{
carried.emplace_back(m, p, eidx);
continue;
}
while(!carried.empty())
{
auto& [cm, cp, ceidx] = carried.back();
if(m - cm >= max_flow_dist)
{
for(auto it = carried.cbegin(); it != carried.cend(); it++)
{
flow = fabs(std::get<1>(*it));
const size_t target_idx = std::get<2>(*it);
flow *= target_idx == 0 ? abyss_depth_exp : abyss_depth_the;
ret_gradient[target_idx] += flow;
condemned += flow;
std::cout << "condemnin': " << m << " " << p << " " << eidx << " | " << cm << " " << cp << " " << ceidx << '\n';
}
carried.clear();
break;
}
if((cp+p)*p > 0.0)
{
flow = fabs((m-cm)*cp);
accd += flow;
p += cp;
ret_gradient[ceidx] += flow;
carried.pop_back();
}
else
{
flow = fabs((m-cm)*p);
accd += flow;
cp += p;
ret_gradient[eidx] += flow;
p = 0.0;
break;
}
}
if(p != 0.0)
carried.emplace_back(m, p, eidx);
//std::cout << m << " " << p << std::endl;
}
for(auto it = carried.cbegin(); it != carried.cend(); it++)
condemned += fabs(std::get<1>(*it));
return accd + condemned * (abyss_depth_exp + abyss_depth_the) * 0.5;
}
#endif
std::tuple<double, double, double> FixedEnvelope::WassersteinMatch(FixedEnvelope& other, double flow_distance, double other_scale)
{
if(_confs_no == 0)
return {0.0, other.get_total_prob() * other_scale, 0.0};
double unmatched1 = 0.0;
double unmatched2 = 0.0;
double massflow = 0.0;
sort_by_mass();
other.sort_by_mass();
size_t idx_this = 0;
size_t idx_other = 0;
double used_prob_this = 0.0;
double used_prob_other = 0.0;
while(idx_this < _confs_no && idx_other < other._confs_no)
{
bool moved = true;
while(moved && idx_this < _confs_no && idx_other < other._confs_no)
{
moved = false;
if(_masses[idx_this] < other._masses[idx_other] - flow_distance)
{
unmatched1 += _probs[idx_this] - used_prob_this;
used_prob_this = 0.0;
idx_this++;
moved = true;
}
if(other._masses[idx_other] < _masses[idx_this] - flow_distance)
{
unmatched2 += other._probs[idx_other]*other_scale - used_prob_other;
used_prob_other = 0.0;
idx_other++;
moved = true;
}
}
if(idx_this < _confs_no && idx_other < other._confs_no)
{
assert(_probs[idx_this] - used_prob_this >= 0.0);
assert(other._probs[idx_other]*other_scale - used_prob_other >= 0.0);
if(_probs[idx_this] - used_prob_this < other._probs[idx_other]*other_scale - used_prob_other)
{
massflow += _probs[idx_this] - used_prob_this;
used_prob_other += _probs[idx_this] - used_prob_this;
assert(used_prob_other >= 0.0);
used_prob_this = 0.0;
idx_this++;
}
else
{
massflow += other._probs[idx_other]*other_scale - used_prob_other;
used_prob_this += other._probs[idx_other]*other_scale - used_prob_other;
assert(used_prob_this >= 0.0);
used_prob_other = 0.0;
idx_other++;
}
}
}
unmatched1 -= used_prob_this;
unmatched2 -= used_prob_other;
for(; idx_this < _confs_no; idx_this++)
unmatched1 += _probs[idx_this];
for(; idx_other < other._confs_no; idx_other++)
unmatched2 += other._probs[idx_other]*other_scale;
return {unmatched1, unmatched2, massflow};
}
FixedEnvelope FixedEnvelope::bin(double bin_width, double middle)
{
sort_by_mass();
FixedEnvelope ret;
if(_confs_no == 0)
return ret;
ret.reallocate_memory<false>(ISOSPEC_INIT_TABLE_SIZE);
if(bin_width == 0)
{
double curr_mass = _masses[0];
double accd_prob = _probs[0];
for(size_t ii = 1; ii<_confs_no; ii++)
{
if(curr_mass != _masses[ii])
{
ret.store_conf(curr_mass, accd_prob);
curr_mass = _masses[ii];
accd_prob = _probs[ii];
}
else
accd_prob += _probs[ii];
}
ret.store_conf(curr_mass, accd_prob);
return ret;
}
size_t ii = 0;
double half_width = 0.5*bin_width;
double hwmm = half_width-middle;
while(ii < _confs_no)
{
double current_bin_middle = floor((_masses[ii]+hwmm)/bin_width)*bin_width + middle;
double current_bin_end = current_bin_middle + half_width;
double bin_prob = 0.0;
while(ii < _confs_no && _masses[ii] <= current_bin_end)
{
bin_prob += _probs[ii];
ii++;
}
ret.store_conf(current_bin_middle, bin_prob);
}
return ret;
}
template<bool tgetConfs> void FixedEnvelope::reallocate_memory(size_t new_size)
{
current_size = new_size;
// FIXME: Handle overflow gracefully here. It definitely could happen for people still stuck on 32 bits...
_masses = reinterpret_cast<double*>(realloc(_masses, new_size * sizeof(double)));
if(_masses == nullptr)
throw std::bad_alloc();
tmasses = _masses + _confs_no;
_probs = reinterpret_cast<double*>(realloc(_probs, new_size * sizeof(double)));
if(_probs == nullptr)
throw std::bad_alloc();
tprobs = _probs + _confs_no;
constexpr_if(tgetConfs)
{
_confs = reinterpret_cast<int*>(realloc(_confs, new_size * allDimSizeofInt));
if(_confs == nullptr)
throw std::bad_alloc();
tconfs = _confs + (allDim * _confs_no);
}
}
void FixedEnvelope::slow_reallocate_memory(size_t new_size)
{
current_size = new_size;
// FIXME: Handle overflow gracefully here. It definitely could happen for people still stuck on 32 bits...
_masses = reinterpret_cast<double*>(realloc(_masses, new_size * sizeof(double)));
if(_masses == nullptr)
throw std::bad_alloc();
tmasses = _masses + _confs_no;
_probs = reinterpret_cast<double*>(realloc(_probs, new_size * sizeof(double)));
if(_probs == nullptr)
throw std::bad_alloc();
tprobs = _probs + _confs_no;
if(_confs != nullptr)
{
_confs = reinterpret_cast<int*>(realloc(_confs, new_size * allDimSizeofInt));
if(_confs == nullptr)
throw std::bad_alloc();
tconfs = _confs + (allDim * _confs_no);
}
}
template<bool tgetConfs> void FixedEnvelope::threshold_init(Iso&& iso, double threshold, bool absolute)
{
IsoThresholdGenerator generator(std::move(iso), threshold, absolute);
size_t tab_size = generator.count_confs();
this->allDim = generator.getAllDim();
this->allDimSizeofInt = this->allDim * sizeof(int);
this->reallocate_memory<tgetConfs>(tab_size);
double* ttmasses = this->_masses;
double* ttprobs = this->_probs;
ISOSPEC_MAYBE_UNUSED int* ttconfs;
constexpr_if(tgetConfs)
ttconfs = _confs;
while(generator.advanceToNextConfiguration())
{
*ttmasses = generator.mass(); ttmasses++;
*ttprobs = generator.prob(); ttprobs++;
constexpr_if(tgetConfs) { generator.get_conf_signature(ttconfs); ttconfs += allDim; }
}
this->_confs_no = tab_size;
}
template void FixedEnvelope::threshold_init<true>(Iso&& iso, double threshold, bool absolute);
template void FixedEnvelope::threshold_init<false>(Iso&& iso, double threshold, bool absolute);
template<bool tgetConfs> void FixedEnvelope::total_prob_init(Iso&& iso, double target_total_prob, bool optimize)
{
if(target_total_prob <= 0.0)
return;
if(target_total_prob >= 1.0)
{
threshold_init<tgetConfs>(std::move(iso), 0.0, true);
return;
}
IsoLayeredGenerator generator(std::move(iso), 1000, 1000, true, std::min<double>(target_total_prob, 0.9999));
this->allDim = generator.getAllDim();
this->allDimSizeofInt = this->allDim*sizeof(int);
this->reallocate_memory<tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);
size_t last_switch = 0;
double prob_at_last_switch = 0.0;
double prob_so_far = 0.0;
double layer_delta;
const double sum_above = log1p(-target_total_prob) - 2.3025850929940455; // log(0.1);
do
{ // Store confs until we accumulate more prob than needed - and, if optimizing,
// store also the rest of the last layer
while(generator.advanceToNextConfigurationWithinLayer())
{
this->template addConfILG<tgetConfs>(generator);
prob_so_far += *(tprobs-1); // The just-stored probability
if(prob_so_far >= target_total_prob)
{
if (optimize)
{
while(generator.advanceToNextConfigurationWithinLayer())
this->template addConfILG<tgetConfs>(generator);
break;
}
else
return;
}
}
if(prob_so_far >= target_total_prob)
break;
last_switch = this->_confs_no;
prob_at_last_switch = prob_so_far;
layer_delta = sum_above - log1p(-prob_so_far);
layer_delta = (std::max)((std::min)(layer_delta, -0.1), -5.0);
} while(generator.nextLayer(layer_delta));
if(!optimize || prob_so_far <= target_total_prob)
return;
// Right. We have extra configurations and we have been asked to produce an optimal p-set, so
// now we shall trim unneeded configurations, using an algorithm dubbed "quicktrim"
// - similar to the quickselect algorithm, except that we use the cumulative sum of elements
// left of pivot to decide whether to go left or right, instead of the positional index.
// We'll be sorting by the prob array, permuting the other ones in parallel.
int* conf_swapspace = nullptr;
constexpr_if(tgetConfs)
conf_swapspace = reinterpret_cast<int*>(malloc(this->allDimSizeofInt));
size_t start = last_switch;
size_t end = this->_confs_no;
double sum_to_start = prob_at_last_switch;
while(start < end)
{
// Partition part
size_t len = end - start;
#if ISOSPEC_BUILDING_R
size_t pivot = len/2 + start;
#else
size_t pivot = random_gen() % len + start; // Using Mersenne twister directly - we don't
// need a very uniform distribution just for pivot
// selection
#endif
double pprob = this->_probs[pivot];
swap<tgetConfs>(pivot, end-1, conf_swapspace);
double new_csum = sum_to_start;
size_t loweridx = start;
for(size_t ii = start; ii < end-1; ii++)
if(this->_probs[ii] > pprob)
{
swap<tgetConfs>(ii, loweridx, conf_swapspace);
new_csum += this->_probs[loweridx];
loweridx++;
}
swap<tgetConfs>(end-1, loweridx, conf_swapspace);
// Selection part
if(new_csum < target_total_prob)
{
start = loweridx + 1;
sum_to_start = new_csum + this->_probs[loweridx];
}
else
end = loweridx;
}
constexpr_if(tgetConfs)
free(conf_swapspace);
if(end <= current_size/2)
// Overhead in memory of 2x or more, shrink to fit
this->template reallocate_memory<tgetConfs>(end);
this->_confs_no = end;
}
template void FixedEnvelope::total_prob_init<true>(Iso&& iso, double target_total_prob, bool optimize);
template void FixedEnvelope::total_prob_init<false>(Iso&& iso, double target_total_prob, bool optimize);
template<bool tgetConfs> void FixedEnvelope::stochastic_init(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias)
{
IsoStochasticGenerator generator(std::move(iso), _no_molecules, _precision, _beta_bias);
this->allDim = generator.getAllDim();
this->allDimSizeofInt = this->allDim * sizeof(int);
this->reallocate_memory<tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);
while(generator.advanceToNextConfiguration())
addConfILG<tgetConfs, IsoStochasticGenerator>(generator);
}
template void FixedEnvelope::stochastic_init<true>(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias);
template void FixedEnvelope::stochastic_init<false>(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias);
double FixedEnvelope::empiric_average_mass()
{
double ret = 0.0;
for(size_t ii = 0; ii < _confs_no; ii++)
{
ret += _masses[ii] * _probs[ii];
}
return ret / get_total_prob();
}
double FixedEnvelope::empiric_variance()
{
double ret = 0.0;
double avg = empiric_average_mass();
for(size_t ii = 0; ii < _confs_no; ii++)
{
double msq = _masses[ii] - avg;
ret += msq * msq * _probs[ii];
}
return ret / get_total_prob();
}
FixedEnvelope FixedEnvelope::Binned(Iso&& iso, double target_total_prob, double bin_width, double bin_middle)
{
FixedEnvelope ret;
double min_mass = iso.getLightestPeakMass();
double range_len = iso.getHeaviestPeakMass() - min_mass;
size_t no_bins = static_cast<size_t>(range_len / bin_width) + 2;
double half_width = 0.5*bin_width;
double hwmm = half_width-bin_middle;
size_t idx_min = floor((min_mass + hwmm) / bin_width);
size_t idx_max = idx_min + no_bins;
double* acc;
# if ISOSPEC_GOT_MMAN
acc = reinterpret_cast<double*>(mmap(nullptr, sizeof(double)*no_bins, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
#else
// This will probably crash for large molecules and high resolutions...
acc = reinterpret_cast<double*>(calloc(no_bins, sizeof(double)));
#endif
if(acc == NULL)
throw std::bad_alloc();
acc -= idx_min;
IsoLayeredGenerator ITG(std::move(iso));
bool non_empty;
while((non_empty = ITG.advanceToNextConfiguration()) && ITG.prob() == 0.0)
{}
if(non_empty)
{
double accum_prob = ITG.prob();
size_t nonzero_idx = floor((ITG.mass() + hwmm)/bin_width);
acc[nonzero_idx] = accum_prob;
while(ITG.advanceToNextConfiguration() && accum_prob < target_total_prob)
{
double prob = ITG.prob();
accum_prob += prob;
size_t bin_idx = floor((ITG.mass() + hwmm)/bin_width);
acc[bin_idx] += prob;
}
// Making the assumption that there won't be gaps of more than 10 Da in the spectrum. This is true for all
// molecules made of natural elements.
size_t distance_10da = static_cast<size_t>(10.0/bin_width) + 1;
size_t empty_steps = 0;
ret.reallocate_memory<false>(ISOSPEC_INIT_TABLE_SIZE);
for(size_t ii = nonzero_idx; ii >= idx_min && empty_steps < distance_10da; ii--)
{
if(acc[ii] > 0.0)
{
empty_steps = 0;
ret.store_conf(static_cast<double>(ii)*bin_width + bin_middle, acc[ii]);
}
else
empty_steps++;
}
empty_steps = 0;
for(size_t ii = nonzero_idx+1; ii < idx_max && empty_steps < distance_10da; ii++)
{
if(acc[ii] > 0.0)
{
empty_steps = 0;
ret.store_conf(static_cast<double>(ii)*bin_width + bin_middle, acc[ii]);
}
else
empty_steps++;
}
}
acc += idx_min;
# if ISOSPEC_GOT_MMAN
munmap(acc, sizeof(double)*no_bins);
#else
free(acc);
#endif
return ret;
}
} // namespace IsoSpec
IsoSpec::FixedEnvelope
Definition fixedEnvelopes.h:40
IsoSpec
Definition allocator.cpp:21