1 /*
2 Copyright 2007, 2008 Daniel Zerbino (zerbino@ebi.ac.uk)
3
4 This file is part of Velvet.
5
6 Velvet 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 2 of the License, or
9 (at your option) any later version.
10
11 Velvet 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 Velvet; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19
20 */
21 #include <stdlib.h>
22 #include <stdio.h>
23
24 #include "globals.h"
25 #include "preGraph.h"
26 #include "utility.h"
27
28 // Replaces two consecutive preNodes into a single equivalent preNode
29 // The extra memory is freed
concatenatePreNodes(IDnum preNodeAID,PreArc * oldPreArc,PreGraph * preGraph)30 static void concatenatePreNodes(IDnum preNodeAID, PreArc * oldPreArc,
31 PreGraph * preGraph)
32 {
33 IDnum preNodeBID = preNodeAID;
34 IDnum currentPreNodeID, nextPreNodeID;
35 PreArc *preArc = oldPreArc;
36 Coordinate totalLength = 0;
37 Coordinate arrayLength;
38 Descriptor * descr, * ptr;
39 int writeOffset = 0;
40 int wordLength = getWordLength_pg(preGraph);
41
42 //printf("Concatenating nodes %li and %li\n", preNodeAID, preNodeBID);
43
44 while(hasSinglePreArc_pg(preNodeBID, preGraph)
45 &&
46 hasSinglePreArc_pg(getOtherEnd_pg
47 (preArc, preNodeBID),
48 preGraph)
49 && !isLoop_pg(preArc)
50 && getDestination_pg(preArc, preNodeBID) != preNodeAID) {
51
52 totalLength += getPreNodeLength_pg(preNodeBID, preGraph);
53 preNodeBID = getDestination_pg(preArc, preNodeBID);
54 preArc = getPreArc_pg(preNodeBID, preGraph);
55 }
56 totalLength += getPreNodeLength_pg(preNodeBID, preGraph);
57 totalLength += wordLength - 1;
58
59 // Descriptor management (preNode)
60 arrayLength = totalLength / 4;
61 if (totalLength % 4)
62 arrayLength++;
63 descr = callocOrExit(arrayLength, Descriptor);
64 ptr = descr;
65 if (preNodeAID > 0) {
66 currentPreNodeID = preNodeAID;
67 appendDescriptors_pg(&ptr, &writeOffset, currentPreNodeID, preGraph, true);
68 preArc = getPreArc_pg(currentPreNodeID, preGraph);
69 currentPreNodeID = getDestination_pg(preArc, currentPreNodeID);
70 while (currentPreNodeID != preNodeBID) {
71 appendDescriptors_pg(&ptr, &writeOffset, currentPreNodeID, preGraph, false);
72 preArc = getPreArc_pg(currentPreNodeID, preGraph);
73 currentPreNodeID = getDestination_pg(preArc, currentPreNodeID);
74 }
75 appendDescriptors_pg(&ptr, &writeOffset, currentPreNodeID, preGraph, false);
76 } else {
77 currentPreNodeID = -preNodeBID;
78 appendDescriptors_pg(&ptr, &writeOffset ,currentPreNodeID, preGraph, true);
79 preArc = getPreArc_pg(currentPreNodeID, preGraph);
80 currentPreNodeID = getDestination_pg(preArc, currentPreNodeID);
81 while (currentPreNodeID != -preNodeAID) {
82 appendDescriptors_pg(&ptr, &writeOffset ,currentPreNodeID, preGraph, false);
83 preArc = getPreArc_pg(currentPreNodeID, preGraph);
84 currentPreNodeID = getDestination_pg(preArc, currentPreNodeID);
85 }
86 appendDescriptors_pg(&ptr, &writeOffset ,currentPreNodeID, preGraph, false);
87 }
88
89 if (writeOffset != 0)
90 while (writeOffset++ != 4)
91 (*ptr) >>= 2;
92
93 setPreNodeDescriptor_pg(descr, totalLength - wordLength + 1, preNodeAID, preGraph);
94
95 // Correct preArcs
96 for (preArc = getPreArc_pg(preNodeBID, preGraph); preArc != NULL;
97 preArc = getNextPreArc_pg(preArc, preNodeBID)) {
98 if (getDestination_pg(preArc, preNodeBID) != -preNodeBID)
99 createAnalogousPreArc_pg(preNodeAID,
100 getDestination_pg(preArc,
101 preNodeBID),
102 preArc, preGraph);
103 else
104 createAnalogousPreArc_pg(preNodeAID, -preNodeAID,
105 preArc, preGraph);
106 }
107
108 // Freeing gobbled preNode
109 currentPreNodeID = -preNodeBID;
110 while (currentPreNodeID != -preNodeAID) {
111 preArc = getPreArc_pg(currentPreNodeID, preGraph);
112 nextPreNodeID = getDestination_pg(preArc, currentPreNodeID);
113 destroyPreNode_pg(currentPreNodeID, preGraph);
114 currentPreNodeID = nextPreNodeID;
115 }
116 }
117
118 // Detects sequences that could be simplified through concatentation
119 // Iterates till preGraph cannot be more simplified
120 // Useless preNodes are freed from memory and remaining ones are renumbered
concatenatePreGraph_pg(PreGraph * preGraph)121 void concatenatePreGraph_pg(PreGraph * preGraph)
122 {
123 IDnum preNodeIndex;
124 PreArc *preArc;
125 PreNode *preNode;
126
127 puts("Concatenation...");
128
129 for (preNodeIndex = 1; preNodeIndex < preNodeCount_pg(preGraph);
130 preNodeIndex++) {
131 preNode = getPreNodeInPreGraph_pg(preGraph, preNodeIndex);
132
133 if (preNode == NULL)
134 continue;
135
136 preArc = getPreArc_pg(preNodeIndex, preGraph);
137
138 while (hasSinglePreArc_pg(preNodeIndex, preGraph)
139 &&
140 hasSinglePreArc_pg(getOtherEnd_pg
141 (preArc, preNodeIndex),
142 preGraph)) {
143 if (isLoop_pg(preArc))
144 break;
145 concatenatePreNodes(preNodeIndex, preArc,
146 preGraph);
147 preArc = getPreArc_pg(preNodeIndex, preGraph);
148 }
149
150 preArc = getPreArc_pg(-preNodeIndex, preGraph);
151
152 while (hasSinglePreArc_pg(-preNodeIndex, preGraph)
153 &&
154 hasSinglePreArc_pg(getOtherEnd_pg
155 (preArc, -preNodeIndex),
156 preGraph)) {
157 if (isLoop_pg(preArc))
158 break;
159 concatenatePreNodes(-preNodeIndex, preArc,
160 preGraph);
161 preArc = getPreArc_pg(-preNodeIndex, preGraph);
162 }
163 }
164
165 renumberPreNodes_pg(preGraph);
166 puts("Concatenation over!");
167 }
168
isEligibleTip(IDnum index,PreGraph * preGraph,Coordinate cutoffLength)169 static boolean isEligibleTip(IDnum index, PreGraph * preGraph, Coordinate
170 cutoffLength)
171 {
172 IDnum currentIndex = -index;
173 Coordinate totalLength = 0;
174 PreArc *activeArc = NULL;
175 PreArc *arc;
176 IDnum mult = 0;
177
178 if (getPreArc_pg(index, preGraph) != NULL)
179 return false;
180
181 // Finding first tangle
182 while (currentIndex != 0
183 && simplePreArcCount_pg(-currentIndex, preGraph) < 2
184 && simplePreArcCount_pg(currentIndex, preGraph) < 2) {
185 totalLength += getPreNodeLength_pg(currentIndex, preGraph);
186 activeArc = getPreArc_pg(currentIndex, preGraph);
187 currentIndex = getDestination_pg(activeArc, currentIndex);
188 }
189
190 // If too long
191 if (totalLength >= cutoffLength)
192 return false;
193
194 // If isolated snippet:
195 if (currentIndex == 0)
196 return true;
197
198 // Joined tips
199 if (simplePreArcCount_pg(-currentIndex, preGraph) < 2)
200 return false;
201
202 // If unique event
203 if (getMultiplicity_pg(activeArc) == 1)
204 return true;
205
206 // Computing max arc
207 for (arc = getPreArc_pg(-currentIndex, preGraph); arc != NULL;
208 arc = getNextPreArc_pg(arc, -currentIndex))
209 if (getMultiplicity_pg(arc) > mult)
210 mult = getMultiplicity_pg(arc);
211
212 // Testing for minority
213 return mult > getMultiplicity_pg(activeArc);
214 }
215
clipTips_pg(PreGraph * preGraph)216 void clipTips_pg(PreGraph * preGraph)
217 {
218 IDnum index;
219 PreNode *current;
220 boolean modified = true;
221 Coordinate cutoffLength = getWordLength_pg(preGraph) * 2;
222 IDnum counter = 0;
223
224 puts("Clipping short tips off preGraph");
225
226 while (modified) {
227 modified = false;
228 for (index = 1; index <= preNodeCount_pg(preGraph);
229 index++) {
230 current = getPreNodeInPreGraph_pg(preGraph, index);
231
232 if (current == NULL)
233 continue;
234
235 if (isEligibleTip(index, preGraph, cutoffLength)
236 || isEligibleTip(-index, preGraph,
237 cutoffLength)) {
238 counter++;
239 destroyPreNode_pg(index, preGraph);
240 modified = true;
241 }
242 }
243 }
244
245 concatenatePreGraph_pg(preGraph);
246 printf("%d tips cut off\n", counter);
247 printf("%d nodes left\n", preNodeCount_pg(preGraph));
248 }
249