BackP.cc 11.8 KB
Newer Older
Marie-Josee Cros's avatar
Marie-Josee Cros committed
1
2
// ------------------------------------------------------------------
// Copyright (C) 2004 INRA <eugene@ossau.toulouse.inra.fr>
Thomas Schiex's avatar
Thomas Schiex committed
3
//
Marie-Josee Cros's avatar
Marie-Josee Cros committed
4
5
// This program is open source; you can redistribute it and/or modify
// it under the terms of the Artistic License (see LICENSE file).
Thomas Schiex's avatar
Thomas Schiex committed
6
//
Marie-Josee Cros's avatar
Marie-Josee Cros committed
7
8
9
// This program 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. 
Thomas Schiex's avatar
Thomas Schiex committed
10
//
Marie-Josee Cros's avatar
Marie-Josee Cros committed
11
12
// You should have received a copy of Artistic License along with
// this program; if not, please see http://www.opensource.org
Thomas Schiex's avatar
Thomas Schiex committed
13
//
Marie-Josee Cros's avatar
Marie-Josee Cros committed
14
15
16
17
18
// $Id$
// ------------------------------------------------------------------
// File:     BackP.cc
// Contents: Definitions for a linear time shortest path with constraint alg.
// ------------------------------------------------------------------
Thomas Schiex's avatar
Thomas Schiex committed
19
20

#include <stdio.h>
21
#include <math.h>
22
#include "assert.h"
Marie-Josee Cros's avatar
Marie-Josee Cros committed
23
24

#include "BackP.h"
25
#include "System.h"
Thomas Schiex's avatar
Thomas Schiex committed
26

27
28
29
30
31
32
33
34
35
36
37
38
// Bit masks in the Status field
const char OptimalMask = 0x01;
const char MarkedMask  = 0x02;
// ----------------------------------------------------------------
// Status handling
// ----------------------------------------------------------------
inline char BackPoint :: IsOptimal() { return (Status & OptimalMask); }
inline char BackPoint :: IsMarked() { return (Status & MarkedMask); }
inline void BackPoint :: SetOptimal() { Status |= OptimalMask; }
inline void BackPoint :: SetMark() { Status |= MarkedMask; }
inline void BackPoint :: ClearOptimal(){ Status &= (0xff^OptimalMask); }
inline void BackPoint :: ClearMark(){ Status &= (0xff^MarkedMask); }
Thomas Schiex's avatar
Thomas Schiex committed
39
40
41
42
// ----------------------------------------------------------------
// Default constructor.
// ----------------------------------------------------------------
BackPoint :: BackPoint  ()
43
44
{
  State = -1;
45
  StartPos = 0;
46
47
48
49
  Cost = 0.0;
  Additional = 0.0;
  Next = Prev = this;
  Origin = NULL;
50
  Status = 0;
51
}
Thomas Schiex's avatar
Thomas Schiex committed
52
53
54
55
56
57
58
59
60
61
// ----------------------------------------------------------------
//  Default constructor.
// ----------------------------------------------------------------
BackPoint :: BackPoint  (char state, int pos, double cost)
{
  State = state;
  StartPos = pos;
  Cost = cost;
  Additional = 0.0;
  Next = Prev = Origin = NULL;
62
  Status = 0;
Thomas Schiex's avatar
Thomas Schiex committed
63
}
Philippe Bardou's avatar
Philippe Bardou committed
64

Thomas Schiex's avatar
Thomas Schiex committed
65
66
67
68
69
70
71
72
// ----------------------------------------------------------------
//  Default destructor.
// ----------------------------------------------------------------
BackPoint :: ~BackPoint  ()
{
  Prev->Next = Next;
  Next->Prev = Prev;
}
Philippe Bardou's avatar
Philippe Bardou committed
73

sallet's avatar
sallet committed
74
75
// ----------------------------------------------------------------
//  Reinit Backpoint
76
//   NOT USED
sallet's avatar
sallet committed
77
78
79
80
81
82
83
// ----------------------------------------------------------------
void BackPoint :: Clean  ()
{
    State = -1;
    StartPos = 0;
    Cost = 0.0;
    Additional = 0.0;
sallet's avatar
sallet committed
84
85
86
    //if (Next != NULL) delete Next;
    //if (Prev != NULL) delete Prev;
    //if (Prev != NULL) delete Origin;
sallet's avatar
sallet committed
87
88
89
90
91
    Next = Prev = this;
    Origin = NULL;
    Status = 0;
}

Thomas Schiex's avatar
Thomas Schiex committed
92
// ----------------------------------------------------------------
93
94
95
96
97
98
// Prints the BackPoint contents
// ----------------------------------------------------------------
void BackPoint :: Print()
{
  printf("pos = %d, state = %d\n", StartPos,State);
}
99
100
101
102
103
104
105
106


// ----------------------------------------------------------------
// Class Track
// STATIC INIT
// ----------------------------------------------------------------
unsigned int Track::NumBPAlloc = 0; 
unsigned int Track::NumBPCollect = 0;
Thomas Schiex's avatar
Thomas Schiex committed
107
// ----------------------------------------------------------------
108
109
110
// Track creator
// ----------------------------------------------------------------
Track :: Track()
Thomas Schiex's avatar
Thomas Schiex committed
111
{
112
113
  Optimal = 0.0;
  OptPos =0;
114
115
116
117
118
119
120
121
122
123
  PenD = NULL;
}
// ----------------------------------------------------------------
// Track creator: from another Track
// ----------------------------------------------------------------
Track :: Track(Track *other)
{
  Optimal 		= other->Optimal;
  OptPos 		= other->OptPos;
  PenD 		= other->PenD;
124
}
125

126
127
128
129
130
131
// ----------------------------------------------------------------
// Track destructor
// ----------------------------------------------------------------
Track :: ~Track()
{
  Zap();
Thomas Schiex's avatar
Thomas Schiex committed
132
  delete PenD;
133
134
135
136
137
138
139
}
// ----------------------------------------------------------------
// Zap  the path data of a  whole track
// ----------------------------------------------------------------
void Track :: Zap()
{
  BackPoint *Dead = Path.Next,*Skip;
Thomas Schiex's avatar
Thomas Schiex committed
140

141
142
143
144
145
146
  while (Dead != &Path) {
    Skip = Dead->Next;
    delete Dead;
    Dead = Skip;
  }
}
sallet's avatar
sallet committed
147
148
149

// ----------------------------------------------------------------
// Reinit the track
150
// NOT USED
sallet's avatar
sallet committed
151
152
153
// ----------------------------------------------------------------
void Track :: Clean()
{
sallet's avatar
sallet committed
154
    Zap();
sallet's avatar
sallet committed
155
156
157
158
159
    Path.Clean();
    Optimal = 0.0;
    OptPos  = 0;
}

160
// ----------------------------------------------------------------
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
// Clear Mark of all BP in a whole track
// ----------------------------------------------------------------
void Track :: ClearMark(int pos)
{
  BackPoint *It = Path.Next;

  while ((It != &Path) && (It->StartPos >= pos)) {
    It->ClearMark();
    It = It->Next;
  }
}
// ----------------------------------------------------------------
// Mark all backpoint from a useful one following Origin
// ----------------------------------------------------------------
void MarkTrace(int pos, BackPoint *Source) {

  while ((Source) &&  (Source->StartPos >= pos) && 
	 (!Source->IsMarked())) {
    Source->SetMark();
    Source = Source->Origin;
  }
}
// ----------------------------------------------------------------
// Mark all useful BP in a whole track
// ----------------------------------------------------------------
void Track :: Mark(int pos)
{
  BackPoint *It = Path.Next;

  do {
    MarkTrace(pos,It);
    if (isinf(It->Additional)) break;

    if (It->IsOptimal() &&  
195
	(abs(pos-It->StartPos) > PenD->MaxLen)) break;
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213

    It = It -> Next;
  } while (It != Path.Next);
}
// ----------------------------------------------------------------
// Sweep for Mark and Sweep
// ---------------------------------------------------------------- We
// have marked any BP accessible from Origin of "young" (not in the
// linear tail). So remaining BP are old BP unaccessible by any Origin
// backtracing from young BP. They can be deleted once the information
// inside is taken into account in existing BP.
// Among Optimal/Mark/State/StartPos/Cost/additional/Origin/Prev/Next
// only Additional will be "saved" by including it in the Additional
// field of the next BackP structure

void Track :: Sweep(int pos) {

  BackPoint *It = Path.Next;
214
  BackPoint *NextIt;
215
  while ( (It != &Path)){// && (It->StartPos >= pos)) {
216
    NextIt = It->Next;
217
218
219
220
221
222
223
    if (!It->IsMarked()) {
      It->Next->Additional += It->Additional;
      It->Next->Prev = It->Prev;
      It->Prev->Next = It->Next;
      delete It;
      NumBPCollect++;
    }
224
    It = NextIt;
225
226
227
  }
}
// ----------------------------------------------------------------
228
229
230
231
232
233
234
// Insert  a new backpoint
// ----------------------------------------------------------------
void Track :: InsertNew(char state, int pos, double cost, BackPoint *Or)
{
  BackPoint *It = Path.Next;
  
  //  if (cost > NINFINITY) {
235
236
237
238
239
  //  if (cost > Optimal-PenD->MaxDelta) {
  //  if (cost > Optimal-PenD->GetDelta(pos-OptPos)) {
  //  if (cost > Path.Next->Cost+Path.Next->Additional-PenD->GetDelta(pos-Path.Next->StartPos)) {
  if (cost > Optimal-PenD->GetDelta(abs(pos-OptPos)) && 
      cost > Path.Next->Cost+Path.Next->Additional-PenD->GetDelta(abs(pos-Path.Next->StartPos))) {
240
    NumBPAlloc++;
Thomas Schiex's avatar
Thomas Schiex committed
241
    It =  new BackPoint(state,pos,cost);
242
    if (cost > Optimal) { 
243
244
      Optimal = cost; 
      It->SetOptimal();
245
246
247
248
      OptPos = pos;
    }
    It->Next = Path.Next;
    It->Prev = Path.Next->Prev;
Thomas Schiex's avatar
Thomas Schiex committed
249
    It->Origin = Or;
250
251
    Path.Next->Prev = It;
    Path.Next = It;
Thomas Schiex's avatar
Thomas Schiex committed
252
253
254
  }
}
// ----------------------------------------------------------------
255
// Insert  a new backpoint
Thomas Schiex's avatar
Thomas Schiex committed
256
// ----------------------------------------------------------------
257
void Track :: ForceNew(char state, int pos, double cost, BackPoint *Or)
Thomas Schiex's avatar
Thomas Schiex committed
258
{
259
260
  BackPoint *It = Path.Next;
  
261
  NumBPAlloc++;
262
  It =  new BackPoint(state,pos,cost);
263
  if (cost > Optimal) { Optimal =cost; It->SetOptimal();}
264
265
266
267
268
  It->Next = Path.Next;
  It->Prev = Path.Next->Prev;
  It->Origin = Or;
  Path.Next->Prev = It;
  Path.Next = It;
Thomas Schiex's avatar
Thomas Schiex committed
269
270
}
// ----------------------------------------------------------------
271
// Returns the best BackPoint taking into account length penalties
Erika Sallet's avatar
Erika Sallet committed
272
// Update value of *cost
Thomas Schiex's avatar
Thomas Schiex committed
273
// ----------------------------------------------------------------
274
BackPoint *Track :: BestUsable(int pos, double *cost, int pen)
Thomas Schiex's avatar
Thomas Schiex committed
275
{
276
277
278
279
280
281
  BackPoint *BestBP = NULL;
  double BestCost = NINFINITY;
  int Len;
  BackPoint *It = Path.Next;
  double LenPen,Add = 0.0;

Thomas Schiex's avatar
Thomas Schiex committed
282
  do {
283
284
    if (isinf(It->Additional)) break;
    Add += It->Additional;
285
    Len = abs(pos-It->StartPos);
286

287
    // when pen == 0 or Origin is NULL, this means we reach the
288
289
    // extremities of the sequence. Therefore we account for an
    // optimistic penality (MinPen) given that the length can only be
290
291
292
    // longer than the actual length. In all cases, we discount the
    // Slope penalty on the length. We further discount one on extremities 
    // (because -1 and Data_Len+1 are used).
293

294
    if (pen && It->Origin) 
295
      LenPen = (*PenD)[Len] - (PenD->FinalSlope)*Len;
296
    else 
297
      LenPen = PenD->MinPen(Len) - PenD->FinalSlope*(Len-1);
298
299
300
301
302

    if ((Add + It->Cost - LenPen) > BestCost) {
      BestCost = Add+It->Cost - LenPen;
      BestBP = It;
    }
303
    if (It->IsOptimal() && Len >= PenD->MaxLen) break;
Thomas Schiex's avatar
Thomas Schiex committed
304

305
306
307
308
309
310
311
    It = It -> Next;
  } while (It != Path.Next);
  

  *cost = BestCost;
  return BestBP;
}
Thomas Schiex's avatar
Thomas Schiex committed
312
// ----------------------------------------------------------------
313
// BackTrace and build a prediction object
Thomas Schiex's avatar
Thomas Schiex committed
314
// ----------------------------------------------------------------
315
Prediction* Track :: BackTrace (int From, int To, int Forward)
Thomas Schiex's avatar
Thomas Schiex committed
316
{
317
318
319
  std::vector <int>         vPos;
  std::vector <signed char> vState;
  Prediction *pred;
320
  BackPoint  *It;
Philippe Bardou's avatar
Philippe Bardou committed
321
  int  pos;
Thomas Schiex's avatar
Thomas Schiex committed
322
  char etat;
323
  int  prevpos, prevstate;
324

325
326
327
  // put state back on correct transitions for backward predictions
  if (!Forward) {
    It = Path.Next;
328
    etat = It->State;
329
330
331
332
333
334
335
336
337
338
339
    It = It->Origin;
    
    while (It != NULL) {
      prevstate = etat;
      etat = It->State;
      It->State = prevstate;
      It = It->Origin;
    }
  }

  It = Path.Next;
340

341
  // initialisation by the terminal state
342
  pos  = It->StartPos;
Philippe Bardou's avatar
Philippe Bardou committed
343
344
  etat = It->State;
  It   = It->Origin;
345
  if (pos > From) {
346
347
348
    vPos.push_back  ( pos  );
    vState.push_back( etat );
  }
349
350

  prevpos = pos;
351
  prevstate = etat;
352
353

  //  printf("pos %d etat %d CDSlen %d prevpos %d\n", pos,etat,CDSlen,prevpos);
354

355
  while (It != NULL) {
356

357
    pos  = It->StartPos;
358
    etat = It->State;
Philippe Bardou's avatar
Philippe Bardou committed
359
    It   = It->Origin;
360

361
    //    printf("pos %d etat %d CDSlen %d prevpos %d\n", pos,etat,CDSlen,prevpos);
362
363
364
365
    
    // Ignore case where pos=From to not display an initial feature with a length = 0
    // NOTE: if we want to detect incomplete genes, we will need this information!
    if (pos > From) {   
366
367
368
      vPos.push_back  ( pos  );
      vState.push_back( etat );
    }
369
    
370
    prevpos   = pos;
371
    prevstate = etat;
372
  }
373
  
374
  if (Forward) {
375
    vPos[0] -=  1;
376
377
378
    reverse(vState.begin(), vState.end());
    reverse(vPos.begin(),   vPos.end());
  }
379
  pred = new Prediction(From, To, vPos, vState);
380
381
  vPos.clear();
  vState.clear();
Thomas Schiex's avatar
Thomas Schiex committed
382

Philippe Bardou's avatar
Philippe Bardou committed
383
384
  return pred;
}
Thomas Schiex's avatar
Thomas Schiex committed
385
// ----------------------------------------------------------------
386
// Dumps the contents of all the Backpoints in the path
Thomas Schiex's avatar
Thomas Schiex committed
387
// ----------------------------------------------------------------
388
void Track :: Dump ()
Thomas Schiex's avatar
Thomas Schiex committed
389
{
390
  BackPoint* It = Path.Next;
Thomas Schiex's avatar
Thomas Schiex committed
391

392
  printf("Number of BP allocated: %d\n",NumBPAlloc);
393
  printf("Number of BP garbage collected: %d\n",NumBPCollect);
394
  printf("Number of active BP: %d\n\n",NumBPAlloc-NumBPCollect);
Thomas Schiex's avatar
Thomas Schiex committed
395
396

  do {
397
398
    printf("pos = %d, state = %d, cost = %f%s, additional = %f",
	   It->StartPos,It->State,It->Cost,
399
	   (It->IsOptimal() ? "*" : ""), It->Additional);
400
401
402
403
404
    if (It->Origin)
      printf(" Or.state = %d Or.pos = %d",It->Origin->State,It->Origin->StartPos);
    printf("\n");
    It = It->Next;
  }  while (It != Path.Next);
Thomas Schiex's avatar
Thomas Schiex committed
405
}