IsoSpec  1.95
fixedEnvelopes.cpp
1 /*
2  * Copyright (C) 2015-2019 Mateusz Łącki and Michał Startek.
3  *
4  * This file is part of IsoSpec.
5  *
6  * IsoSpec is free software: you can redistribute it and/or modify
7  * it under the terms of the Simplified ("2-clause") BSD licence.
8  *
9  * IsoSpec is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12  *
13  * You should have received a copy of the Simplified BSD Licence
14  * along with IsoSpec. If not, see <https://opensource.org/licenses/BSD-2-Clause>.
15  */
16 
17 #include "fixedEnvelopes.h"
18 
19 namespace IsoSpec
20 {
21 
22 template<bool tgetlProbs, bool tgetMasses, bool tgetProbs, bool tgetConfs> void FixedEnvelope::reallocate_memory(size_t new_size)
23 {
24  // FIXME: Handle overflow gracefully here. It definitely could happen for people still stuck on 32 bits...
25  constexpr_if(tgetlProbs) { _lprobs = (double*) realloc(_lprobs, new_size * sizeof(double)); tlprobs = _lprobs + _confs_no; }
26  constexpr_if(tgetMasses) { _masses = (double*) realloc(_masses, new_size * sizeof(double)); tmasses = _masses + _confs_no; }
27  constexpr_if(tgetProbs) { _probs = (double*) realloc(_probs, new_size * sizeof(double)); tprobs = _probs + _confs_no; }
28  constexpr_if(tgetConfs) { _confs = (int*) realloc(_confs, new_size * allDimSizeofInt); tconfs = _confs + (allDim * _confs_no); }
29 }
30 
31 
32 template<bool tgetlProbs, bool tgetMasses, bool tgetProbs, bool tgetConfs> void ThresholdFixedEnvelope::init(Iso&& iso)
33 {
34  IsoThresholdGenerator generator(std::move(iso), threshold, absolute);
35 
36  size_t tab_size = generator.count_confs();
37  this->allDim = generator.getAllDim();
38  this->allDimSizeofInt = this->allDim * sizeof(int);
39 
40  this->reallocate_memory<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(tab_size);
41 
42  while(generator.advanceToNextConfiguration())
43  store_conf<IsoThresholdGenerator, tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(generator);
44 
45  this->_confs_no = tab_size;
46 }
47 
48 
49 template<bool tgetlProbs, bool tgetMasses, bool tgetProbs, bool tgetConfs> void TotalProbFixedEnvelope::init(Iso&& iso)
50 {
51  if(optimize && !tgetProbs)
52  // If we want to optimize, we need the probs
53  throw std::logic_error("Cannot perform quicktrim if we're not computing probabilities");
54 
55  IsoLayeredGenerator generator(std::move(iso), 1000, 1000, true, std::min<double>(target_total_prob, 0.9999));
56 
57  this->allDim = generator.getAllDim();
58  this->allDimSizeofInt = this->allDim*sizeof(int);
59 
60 
61  this->reallocate_memory<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);
62 
63  size_t last_switch = 0;
64  double prob_at_last_switch = 0.0;
65  double prob_so_far = 0.0;
66  do
67  { // Store confs until we accumulate more prob than needed - and, if optimizing,
68  // store also the rest of the last layer
69  while(generator.advanceToNextConfigurationWithinLayer())
70  {
71  this->template addConf<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(generator);
72  prob_so_far += generator.prob();
73  if(prob_so_far >= target_total_prob)
74  {
75  if (optimize)
76  {
77  while(generator.advanceToNextConfigurationWithinLayer())
78  this->template addConf<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(generator);
79  break;
80  }
81  else
82  return;
83  }
84  }
85  if(prob_so_far >= target_total_prob)
86  break;
87 
88  last_switch = this->_confs_no;
89  prob_at_last_switch = prob_so_far;
90  } while(generator.nextLayer(-3.0));
91 
92  if(!optimize || prob_so_far <= target_total_prob)
93  return;
94 
95  // Right. We have extra configurations and we have been asked to produce an optimal p-set, so
96  // now we shall trim unneeded configurations, using an algorithm dubbed "quicktrim"
97  // - similar to the quickselect algorithm, except that we use the cumulative sum of elements
98  // left of pivot to decide whether to go left or right, instead of the positional index.
99  // We'll be sorting by the prob array, permuting the other ones in parallel.
100 
101  int* conf_swapspace = nullptr;
102  constexpr_if(tgetConfs)
103  conf_swapspace = (int*) malloc(this->allDimSizeofInt);
104 
105  size_t start = last_switch;
106  size_t end = this->_confs_no;
107  double sum_to_start = prob_at_last_switch;
108 
109  while(start < end)
110  {
111  // Partition part
112  size_t len = end - start;
113 #if ISOSPEC_BUILDING_R
114  size_t pivot = len/2 + start;
115 #else
116  size_t pivot = rand() % len + start;
117 #endif
118  double pprob = this->_probs[pivot];
119  swap<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(pivot, end-1, conf_swapspace);
120 
121  double new_csum = sum_to_start;
122 
123  size_t loweridx = start;
124  for(size_t ii=start; ii<end-1; ii++)
125  if(this->_probs[ii] > pprob)
126  {
127  swap<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(ii, loweridx, conf_swapspace);
128  new_csum += this->_probs[loweridx];
129  loweridx++;
130  }
131 
132  swap<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(end-1, loweridx, conf_swapspace);
133 
134  // Selection part
135  if(new_csum < target_total_prob)
136  {
137  start = loweridx + 1;
138  sum_to_start = new_csum + this->_probs[loweridx];
139  }
140  else
141  end = loweridx;
142  }
143 
144  constexpr_if(tgetConfs)
145  free(conf_swapspace);
146 
147  if(end <= current_size/2)
148  // Overhead in memory of 2x or more, shrink to fit
149  this->template reallocate_memory<tgetlProbs, tgetMasses, tgetProbs, tgetConfs>(end);
150 
151  this->_confs_no = end;
152 }
153 
154 template<typename T> void call_init(T* tabulator, Iso&& iso, bool tgetlProbs, bool tgetMasses, bool tgetProbs, bool tgetConfs)
155 {
156  if(tgetlProbs)
157  {
158  if(tgetMasses)
159  {
160  if(tgetProbs)
161  {
162  if(tgetConfs)
163  tabulator->template init<true, true, true, true>(std::move(iso));
164  else
165  tabulator->template init<true, true, true, false>(std::move(iso));
166  }
167  else
168  {
169  if(tgetConfs)
170  tabulator->template init<true, true, false, true>(std::move(iso));
171  else
172  tabulator->template init<true, true, false, false>(std::move(iso));
173  }
174  }
175  else
176  {
177  if(tgetProbs)
178  {
179  if(tgetConfs)
180  tabulator->template init<true, false, true, true>(std::move(iso));
181  else
182  tabulator->template init<true, false, true, false>(std::move(iso));
183  }
184  else
185  {
186  if(tgetConfs)
187  tabulator->template init<true, false, false, true>(std::move(iso));
188  else
189  tabulator->template init<true, false, false, false>(std::move(iso));
190  }
191  }
192  }
193  else
194  {
195  if(tgetMasses)
196  {
197  if(tgetProbs)
198  {
199  if(tgetConfs)
200  tabulator->template init<false, true, true, true>(std::move(iso));
201  else
202  tabulator->template init<false, true, true, false>(std::move(iso));
203  }
204  else
205  {
206  if(tgetConfs)
207  tabulator->template init<false, true, false, true>(std::move(iso));
208  else
209  tabulator->template init<false, true, false, false>(std::move(iso));
210  }
211  }
212  else
213  {
214  if(tgetProbs)
215  {
216  if(tgetConfs)
217  tabulator->template init<false, false, true, true>(std::move(iso));
218  else
219  tabulator->template init<false, false, true, false>(std::move(iso));
220  }
221  else
222  {
223  if(tgetConfs)
224  tabulator->template init<false, false, false, true>(std::move(iso));
225  else
226  tabulator->template init<false, false, false, false>(std::move(iso));
227  }
228  }
229  }
230 }
231 
232 template void call_init<TotalProbFixedEnvelope>(TotalProbFixedEnvelope* tabulator, Iso&& iso, bool tgetlProbs, bool tgetMasses, bool tgetProbs, bool tgetConfs);
233 template void call_init<ThresholdFixedEnvelope>(ThresholdFixedEnvelope* tabulator, Iso&& iso, bool tgetlProbs, bool tgetMasses, bool tgetProbs, bool tgetConfs);
234 
235 
236 } // namespace IsoSpec
IsoSpec
Definition: allocator.cpp:21