|
1
|
1 /*
|
|
|
2 Copyright 2009-2010 Stéphane De Mita, Mathieu Siol
|
|
|
3
|
|
|
4 This file is part of the EggLib library.
|
|
|
5
|
|
|
6 EggLib is free software: you can redistribute it and/or modify
|
|
|
7 it under the terms of the GNU General Public License as published by
|
|
|
8 the Free Software Foundation, either version 3 of the License, or
|
|
|
9 (at your option) any later version.
|
|
|
10
|
|
|
11 EggLib is distributed in the hope that it will be useful,
|
|
|
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
14 GNU General Public License for more details.
|
|
|
15
|
|
|
16 You should have received a copy of the GNU General Public License
|
|
|
17 along with EggLib. If not, see <http://www.gnu.org/licenses/>.
|
|
|
18 */
|
|
|
19
|
|
|
20 #ifndef EGGLIB_EDGE_HPP
|
|
|
21 #define EGGLIB_EDGE_HPP
|
|
|
22
|
|
|
23 #include <vector>
|
|
|
24 #include <climits>
|
|
|
25 #include "EggException.hpp"
|
|
|
26
|
|
|
27 namespace egglib {
|
|
|
28
|
|
|
29 class Random;
|
|
|
30
|
|
|
31 /** \brief Edge of the ancestral recombination graph
|
|
|
32 *
|
|
|
33 * \ingroup coalesce
|
|
|
34 *
|
|
|
35 * Each Edge instance provides access to its 0, 1 or 2 descendants
|
|
|
36 * (the former holds for a terminal node, the middle for the parent
|
|
|
37 * of a recombined node and the latter for the parent of a coalesced
|
|
|
38 * node (most classical node in the coalescent).The Edge also
|
|
|
39 * provides to the edge length. Note that the Edge instance must be
|
|
|
40 * understood as an ARG node and the branch above it (latter in the
|
|
|
41 * coalescence process). Edge instances also keep track of the list
|
|
|
42 * of descendants descending from this node (which may differ along
|
|
|
43 * recombining segment). Edge instances *must* be created through one
|
|
|
44 * of the "default" and "coalescence" constructors or through the
|
|
|
45 * recombination method. Edge instances should never be copied but
|
|
|
46 * manipulated by references.
|
|
|
47 *
|
|
|
48 */
|
|
|
49 class Edge {
|
|
|
50
|
|
|
51 public:
|
|
|
52
|
|
|
53 /// Destructor
|
|
|
54 virtual ~Edge();
|
|
|
55
|
|
|
56 /** \brief Constructor
|
|
|
57 *
|
|
|
58 * \param numberOfSegments the number of recombining segments
|
|
|
59 * (one for a non-recombining region).
|
|
|
60 *
|
|
|
61 * Use the Pool, instead. Objects are delivered with a
|
|
|
62 * complete coverage.
|
|
|
63 *
|
|
|
64 */
|
|
|
65 Edge(unsigned int numberOfSegments);
|
|
|
66
|
|
|
67
|
|
|
68 /// Restore object to `factory` state
|
|
|
69 void reset();
|
|
|
70
|
|
|
71
|
|
|
72 /** \brief Builds for internal node
|
|
|
73 *
|
|
|
74 * \param date the date of creation of the edge.
|
|
|
75 * \param son1 first edge descending from this edge.
|
|
|
76 * \param son2 second edge descending from this edge.
|
|
|
77 * \param edgesPerSegments counts the current number of
|
|
|
78 * (non-coalesced lineages for each lineages); must have the
|
|
|
79 * appropriate size and will be updated.
|
|
|
80 * \param MRCA the list where to place the address of segment
|
|
|
81 * MRCA, if it occurs.
|
|
|
82 * \param totalLength the total length of the tree.
|
|
|
83 * \param segmentLengths the table of tree lengths per
|
|
|
84 * segment.
|
|
|
85 *
|
|
|
86 * Assumes the current object has the correct number of
|
|
|
87 * segments.
|
|
|
88 *
|
|
|
89 */
|
|
|
90 void coalescence(double date, Edge* son1, Edge* son2,
|
|
|
91 unsigned int* edgesPerSegments, Edge** MRCA,
|
|
|
92 double& totalLength, double* segmentLengths);
|
|
|
93
|
|
|
94
|
|
|
95 /** \brief Generates a recombination event
|
|
|
96 *
|
|
|
97 * \param date the date of the event.
|
|
|
98 * \param dest1 destination for the first resulting edge.
|
|
|
99 * \param dest2 destination for the second resulting edge.
|
|
|
100 * \param random pointer to the Random instance used by the
|
|
|
101 * simulator.
|
|
|
102 * \param totalLength the total length of the tree.
|
|
|
103 * \param segmentLengths the table of tree lengths per
|
|
|
104 * segment.
|
|
|
105 *
|
|
|
106 * dest1 and dest2 must be Edge address initialized with the
|
|
|
107 * appropriate number of segments.
|
|
|
108 *
|
|
|
109 */
|
|
|
110 void recombination(double date,
|
|
|
111 Edge* dest1, Edge* dest2, Random* random,
|
|
|
112 double& totalLength, double* segmentLengths);
|
|
|
113
|
|
|
114
|
|
|
115 /** \brief Define this edge to be terminal
|
|
|
116 *
|
|
|
117 * The edge have only non-covered segments
|
|
|
118 *
|
|
|
119 */
|
|
|
120 void set_terminal(unsigned int leaf_index);
|
|
|
121
|
|
|
122 /// Branch's raw length (doesn't take account segment coverage)
|
|
|
123 double length;
|
|
|
124
|
|
|
125 /// Number of covered segments
|
|
|
126 unsigned int coverage;
|
|
|
127
|
|
|
128 /// Time position of the branch's bottom
|
|
|
129 double bottom;
|
|
|
130
|
|
|
131 /// Time position of the branch's top
|
|
|
132 double top;
|
|
|
133
|
|
|
134 /// Address of the first son
|
|
|
135 Edge* son1;
|
|
|
136
|
|
|
137 /// Address of the second son
|
|
|
138 Edge* son2;
|
|
|
139
|
|
|
140 /// Number of sons (0, 1 or 2)
|
|
|
141 unsigned int numberOfSons;
|
|
|
142
|
|
|
143 /// Checks if a given segment is covered
|
|
|
144 inline bool segment(unsigned int segmentIndex) const {
|
|
|
145 //while(segments[segmentIndex]==0) segmentIndex--;
|
|
|
146 if (segments[segmentIndex]==0) return false;
|
|
|
147 if (segments[segmentIndex]==UINT_MAX) return true;
|
|
|
148 return segbools[segmentIndex];
|
|
|
149 }
|
|
|
150
|
|
|
151 /// leaf index (0 for internal nodes)
|
|
|
152 inline unsigned int label() const {
|
|
|
153 return _label;
|
|
|
154 }
|
|
|
155
|
|
|
156 /// Number of mutations per segment
|
|
|
157 unsigned int* numberOfMutationsPerActualSite;
|
|
|
158
|
|
|
159 /// Sets the actual number of sites
|
|
|
160 void set_actualNumberOfSites(unsigned int actualNumberOfSites);
|
|
|
161
|
|
|
162
|
|
|
163 private:
|
|
|
164
|
|
|
165 unsigned int _label;
|
|
|
166
|
|
|
167 /// Default constructor is not available
|
|
|
168 Edge();
|
|
|
169
|
|
|
170 /// Copy constructor is not available
|
|
|
171 Edge(const Edge& edge) {}
|
|
|
172
|
|
|
173 /// Assignment operator is not available
|
|
|
174 Edge& operator=(const Edge& edge) { return *this; }
|
|
|
175
|
|
|
176 void init(unsigned int numberOfSegments);
|
|
|
177
|
|
|
178 unsigned int numberOfSegments;
|
|
|
179 unsigned int* segments;
|
|
|
180 bool* segbools;
|
|
|
181
|
|
|
182 /// complete the covered ranges with UINT_MAX's and the non-covered by 0's
|
|
|
183 inline void fill() {
|
|
|
184 coverage=0.;
|
|
|
185 unsigned int i=0,j;
|
|
|
186 while (i<numberOfSegments) {
|
|
|
187 if (segbools[i]==true) {
|
|
|
188 coverage += segments[i];
|
|
|
189 for (j=1; j<segments[i]; j++) {
|
|
|
190 segments[i+j] = UINT_MAX;
|
|
|
191 }
|
|
|
192 }
|
|
|
193 i+=segments[i];
|
|
|
194 }
|
|
|
195 }
|
|
|
196
|
|
|
197
|
|
|
198
|
|
|
199 /// update containing Arg's branch lengths
|
|
|
200 inline void update_lengths(double& totalLength, double* segmentLengths) {
|
|
|
201 unsigned int i=0,j;
|
|
|
202 while (i<numberOfSegments) {
|
|
|
203 if (segbools[i]==true) {
|
|
|
204 totalLength += segments[i]*length;
|
|
|
205 for (j=0; j<segments[i]; j++) {
|
|
|
206 segmentLengths[i+j] += length;
|
|
|
207 }
|
|
|
208 }
|
|
|
209 i+=segments[i];
|
|
|
210 }
|
|
|
211 }
|
|
|
212
|
|
|
213
|
|
|
214 };
|
|
|
215
|
|
|
216
|
|
|
217
|
|
|
218
|
|
|
219
|
|
|
220
|
|
|
221
|
|
|
222
|
|
|
223
|
|
|
224
|
|
|
225
|
|
|
226
|
|
|
227 /** \brief Pool of Edge objects
|
|
|
228 *
|
|
|
229 * \ingroup coalesce
|
|
|
230 *
|
|
|
231 * Holds a pool of Edge objects that can be recycled to spare the
|
|
|
232 * building burden. A construction time, a number of Edge objects
|
|
|
233 * equals to the predicted number of needed instances should be
|
|
|
234 * requested. The Edge's will be prebuilt immediately and delivered
|
|
|
235 * upon request. After use, the Edge's should be released. It is only
|
|
|
236 * possible to release the last issued Edge instance or all of them
|
|
|
237 * at once.
|
|
|
238 *
|
|
|
239 */
|
|
|
240 class EdgePool {
|
|
|
241
|
|
|
242 public:
|
|
|
243
|
|
|
244 /// Default constructor (nothing allocated)
|
|
|
245 EdgePool();
|
|
|
246
|
|
|
247
|
|
|
248 /// Destructor
|
|
|
249 virtual ~EdgePool();
|
|
|
250
|
|
|
251
|
|
|
252 /** \brief Configure pool
|
|
|
253 *
|
|
|
254 * Pre-allocates a given number of Edge objects. The objects
|
|
|
255 * will be immediately available.
|
|
|
256 *
|
|
|
257 * Data previously allocated (by a previous call of this
|
|
|
258 * function or by the deliver() method) will be lost so it
|
|
|
259 * can be required to use clear() before.
|
|
|
260 *
|
|
|
261 * \param numberOfSegments the number of segments of the
|
|
|
262 * simulation; all Edge instances will use this value.
|
|
|
263 *
|
|
|
264 * \param numberOfPreAllocated the number of Edge that should
|
|
|
265 * be kept ready for immediate use.
|
|
|
266 *
|
|
|
267 */
|
|
|
268 void set(unsigned int numberOfSegments, unsigned numberOfPreAllocated);
|
|
|
269
|
|
|
270
|
|
|
271 /** \brief Frees internally stored memory
|
|
|
272 *
|
|
|
273 * This invalidate all points that have been delivered
|
|
|
274 * previously. However, any previously set number of segments
|
|
|
275 * (0, by default) is retained.
|
|
|
276 *
|
|
|
277 */
|
|
|
278 void clear();
|
|
|
279
|
|
|
280
|
|
|
281 /** \brief Deliver an Edge
|
|
|
282 *
|
|
|
283 * The object must not be freed by the client! This object is
|
|
|
284 * allocated on the heap if the cache is not large enough,
|
|
|
285 * only reset if it was previously released, or just delivered
|
|
|
286 * if it is one of the initially allocated instances.
|
|
|
287 *
|
|
|
288 */
|
|
|
289 Edge* deliver();
|
|
|
290
|
|
|
291
|
|
|
292 /** \brief Release an Edge
|
|
|
293 *
|
|
|
294 * Release the last delivered Edge. The instance is only
|
|
|
295 * cached for a potential future use; it is not freed nor
|
|
|
296 * reset immediately. If no Edge's are in use, nothing is
|
|
|
297 * done.
|
|
|
298 *
|
|
|
299 */
|
|
|
300 void releaseLast();
|
|
|
301
|
|
|
302
|
|
|
303 /** \brief Release all Edge's
|
|
|
304 *
|
|
|
305 * Release all delivered Edges. The instances are only
|
|
|
306 * cached for a potential future use; they are not freed nor
|
|
|
307 * reset immediately. If no Edge's are in use, nothing is
|
|
|
308 * done.
|
|
|
309 *
|
|
|
310 */
|
|
|
311 void releaseAll();
|
|
|
312
|
|
|
313 private:
|
|
|
314
|
|
|
315 /// Not available
|
|
|
316 EdgePool(const EdgePool& ep) {}
|
|
|
317
|
|
|
318 /// Not available
|
|
|
319 EdgePool& operator=(const EdgePool& ep) { return *this; }
|
|
|
320
|
|
|
321 unsigned int numberOfSegments;
|
|
|
322 unsigned int used;
|
|
|
323 unsigned int released;
|
|
|
324 unsigned int ready;
|
|
|
325 Edge** cache;
|
|
|
326
|
|
|
327 };
|
|
|
328
|
|
|
329 }
|
|
|
330
|
|
|
331 #endif
|
