dynprog.cpp 16.8 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1
/*
2 3 4 5 6 7 8
  This file is part of Vidjil <http://www.vidjil.org>
  Copyright (C) 2011, 2012, 2013, 2014, 2015 by Bonsai bioinformatics 
  at CRIStAL (UMR CNRS 9189, Université Lille) and Inria Lille
  Contributors: 
      Mathieu Giraud <mathieu.giraud@vidjil.org>
      Mikaël Salson <mikael.salson@vidjil.org>
      Marc Duez <marc.duez@vidjil.org>
Mikaël Salson's avatar
Mikaël Salson committed
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

  "Vidjil" is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  "Vidjil" 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.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with "Vidjil". If not, see <http://www.gnu.org/licenses/>
*/

#include "dynprog.h"
#include "tools.h"
26
#include "segment.h"
Mikaël Salson's avatar
Mikaël Salson committed
27 28 29 30
#include <cassert>
#include <list>
#include <cstdlib>
#include <string>
31
#include <cmath>
Mikaël Salson's avatar
Mikaël Salson committed
32

33 34 35 36
#define SUBST '|'
#define MISMATCH '.'
#define INSER 'i'
#define DELET 'd'
37 38
#define INSER_AFFINE 'I'
#define DELET_AFFINE 'D'
39
#define INIT '_'
40 41 42 43
#define HOMO2X 'h'
#define HOMO2Y 'h'
#define FIN ' '
#define BEGIN 'B'
Mikaël Salson's avatar
Mikaël Salson committed
44
#define BACKSIZE 120
Mikaël Salson's avatar
Mikaël Salson committed
45

46 47
#define BAD_BACKTRACK -1

Mikaël Salson's avatar
Mikaël Salson committed
48 49 50 51 52 53 54 55 56 57 58
using namespace std;


Cost::Cost(int match, int mismatch, int indel, int del_end, int homopolymer)
{
  this -> match = match ;
  this -> mismatch = mismatch ;
  this -> insertion = indel ;
  this -> deletion = indel ;
  this -> deletion_end = del_end ;
  this -> homopolymer = (homopolymer == MINUS_INF ? indel: homopolymer);
59

60 61 62
  this -> open_insertion = this -> open_deletion = MINUS_INF ;
  this -> extend_insertion = this -> extend_deletion = MINUS_INF ;
  this -> affine_gap = false ;
63 64

  estimate_K_lambda();
65 66 67 68 69 70 71 72 73 74 75 76 77
}

Cost::Cost(int match, int mismatch, int open_gap, int extend_gap, int del_end, int homopolymer)
{
  this -> match = match ;
  this -> mismatch = mismatch ;
  this -> insertion = MINUS_INF ;
  this -> deletion = MINUS_INF ;
  this -> deletion_end = del_end ;
  this -> homopolymer = homopolymer ;

  this -> open_insertion = this -> open_deletion = open_gap ;
  this -> extend_insertion = this -> extend_deletion = extend_gap ;
78
  this -> affine_gap = true ;
79 80 81 82 83 84 85 86 87

  estimate_K_lambda();
}

void Cost::estimate_K_lambda()
{
  // Rough estimate of K and lambda parameters
  lambda = log(4) / match ;
  K = 0.05 ;
Mikaël Salson's avatar
Mikaël Salson committed
88 89 90 91 92 93 94 95 96
}


ostream& operator<<(ostream& out, const Cost& cost)
{
  out << "(" << cost.match 
      << ", " << cost.mismatch
      << "/" << cost.insertion
      << "/" << cost.deletion
97 98
      << ", " << cost.open_insertion << cost.extend_insertion
      << "/" << cost.open_deletion << cost.extend_deletion
Mikaël Salson's avatar
Mikaël Salson committed
99 100
      << ", " << cost.deletion_end
      << ", " << cost.homopolymer
101 102 103
      << ") "
      << cost.K << "/" << cost.lambda ;

Mikaël Salson's avatar
Mikaël Salson committed
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
  return out;
}

Cost::Cost()
{
}

int Cost::substitution(char a, char b)
{
  return (a == b) ? match : mismatch ;
}

int Cost::homo2(char xa, char xb, char y)
{
  return ((xa == xb) && (xb == y)) ? homopolymer : MINUS_INF ;
}

121
double Cost::toPValue(const int score)
122
{
123 124 125 126 127 128 129
#ifdef DEBUG_EVALUE
  cout << *this
       << " | score " << score
       << ", matches " << (float) score / (float) match
       << " ==> pvalue " << std::scientific << K * exp(-lambda * score) << std::fixed
       << endl;
#endif
130
  return K * exp(-lambda * score);
131 132
}

Mikaël Salson's avatar
Mikaël Salson committed
133

134 135 136
DynProg::DynProg(const string &x, const string &y, DynProgMode mode, const Cost& cost,
                 const bool reverse_x, const bool reverse_y,
                 const int marked_pos_j)
Mikaël Salson's avatar
Mikaël Salson committed
137 138 139 140 141 142 143
{
  this -> x = reverse_x ? reverse(x) : x ;
  this -> y = reverse_y ? reverse(y) : y ;

  this -> reverse_x = reverse_x ;
  this -> reverse_y = reverse_y ;

144 145 146
  this -> marked_pos_j = marked_pos_j;
  this -> marked_pos_i = 0;

Mikaël Salson's avatar
Mikaël Salson committed
147 148 149 150
  m = x.size();
  n = y.size();
  // cout << "Dynprog " << x << "(" << m << ") / " << y << "(" << n << ")" << endl ;

151
  B = new operation*[m+1];
Mikaël Salson's avatar
Mikaël Salson committed
152
  for (int i = 0; i <= m; i++) {
153
    B[i] = new operation[n+1];
Mikaël Salson's avatar
Mikaël Salson committed
154
  }
155 156 157 158 159 160 161 162
  Bins = new operation*[m+1];
  for (int i = 0; i <= m; i++) {
    Bins[i] = new operation[n+1];
  }
  Bdel = new operation*[m+1];
  for (int i = 0; i <= m; i++) {
    Bdel[i] = new operation[n+1];
  }  
Mikaël Salson's avatar
Mikaël Salson committed
163 164 165 166 167 168 169
  this -> mode = mode;
  this -> cost = cost;

  this -> best_i = -1 ;
  this -> best_j = -1 ;
  this -> first_i = -1 ;
  this -> first_j = -1 ;
Mikaël Salson's avatar
Mikaël Salson committed
170 171 172 173
  
  gap1=NULL;
  gap2=NULL;
  linkgap=NULL;
Mikaël Salson's avatar
Mikaël Salson committed
174 175 176 177 178

  init();
}

DynProg::~DynProg() {
179

Mikaël Salson's avatar
Mikaël Salson committed
180 181
  for (int i = 0; i <= m; i++) {
    delete [] B[i];
182 183
    delete [] Bins[i];
    delete [] Bdel[i];
Mikaël Salson's avatar
Mikaël Salson committed
184
  }
185
  delete [] B;  
186 187
  delete [] Bins;
  delete [] Bdel;  
Mikaël Salson's avatar
Mikaël Salson committed
188 189 190
  delete [] gap1;
  delete [] gap2;
  delete [] linkgap;
Mikaël Salson's avatar
Mikaël Salson committed
191 192 193
}

void DynProg::init()
194
{
195 196 197 198 199 200 201 202 203 204 205 206 207 208
  // Initialize scores
  for (int i=0; i<=m; i++)
    for (int j=0; j<=n; j++)
      {
        B[i][j].score = 0 ;
        B[i][j].i = BAD_BACKTRACK ;
        B[i][j].j = BAD_BACKTRACK ;
        if (cost.affine_gap) {
          Bins[i][j].score = 0;
          Bdel[i][j].score = 0;
          Bins[i][j].i = BAD_BACKTRACK;
          Bdel[i][j].i = BAD_BACKTRACK;
          Bins[i][j].j = BAD_BACKTRACK;
          Bdel[i][j].j = BAD_BACKTRACK;
209 210
          Bins[i][j].type = INIT;
          Bdel[i][j].type = INIT;
211 212
        }
      }
213

214
  // Initialize backtrack labels (B[0][0] labels are never used)
215 216 217 218 219 220 221 222 223 224 225 226 227 228
  for (int i=1; i<=m; i++)
    {
      B[i][0].type = INSER ;
      B[i][0].i = i-1 ;
      B[i][0].j = 0 ;
    }

  for (int j=1; j<=n; j++)
    {
      B[0][j].type = DELET ;    
      B[0][j].i = 0 ;
      B[0][j].j = j-1 ;
    }

229
  if (!(mode == Local || mode == LocalEndWithSomeDeletions)) // Global, SemiGlobal
230
    for (int i=0; i<=m; i++) {
231
      B[i][0].score = i * cost.insertion ;
232 233 234
      if (cost.affine_gap) {
        Bins[i][0].score = cost.open_insertion + cost.extend_insertion * i;
        Bdel[i][0].score = MINUS_INF;
235
        B[i][0].score = Bins[i][0].score ;
236
        B[i][0].type = INSER_AFFINE;
237 238
      }

239
      if (mode == GlobalButMostlyLocal) {
240
        B[i][0].score /= 2;
241 242 243 244 245
        if (cost.affine_gap) {
          Bins[i][0].score -= cost.open_insertion / 2 ;
          B[i][0].score = Bins[i][0].score ;
        }
      }
246 247
    }

248
  if (!(mode == SemiGlobal || mode == SemiGlobalTrans || mode == Local || mode == LocalEndWithSomeDeletions)) // Global
249
    for (int j=0; j<=n; j++) {
250
      B[0][j].score = j * cost.deletion ;
251
       if (cost.affine_gap) {
252
         if (j) Bins[0][j].score = MINUS_INF;
253
         Bdel[0][j].score = cost.open_deletion + cost.extend_deletion * j;
254
         B[0][j].score = Bdel[0][j].score;
255
         B[0][j].type = DELET_AFFINE;
256 257
       }

258
       if (mode == GlobalButMostlyLocal) {
259
         B[0][j].score /= 2;
260 261 262 263 264
         if (cost.affine_gap) {
           Bdel[0][j].score -= cost.open_deletion / 2 ;
           B[0][j].score = Bdel[0][j].score;
         }
       }
265
    }
Mikaël Salson's avatar
Mikaël Salson committed
266 267
}

268

269 270 271 272 273 274 275 276 277 278 279
inline void try_operation(operation &best, int type, int i, int j, int score)
{
  if (score > best.score)
    {
      best.type = type ;
      best.i = i ;
      best.j = j ;
      best.score = score ;
    }
}

280
int DynProg::compute(bool onlyBottomTriangle, int onlyBottomTriangleShift)
Mikaël Salson's avatar
Mikaël Salson committed
281 282
{
  best_score = MINUS_INF ;
283 284 285 286
  best_i = 0 ;
  best_j = 0 ;
  
  operation best ;
Mikaël Salson's avatar
Mikaël Salson committed
287 288 289 290

  for (int i=1; i<=m; i++)
    for (int j=1; j<=n; j++)
    {
291
      best.score = MINUS_INF ;
Mikaël Salson's avatar
Mikaël Salson committed
292

293 294 295 296 297 298 299 300 301 302 303 304 305
      // If onlyBottomTriangle, stops when we are not in the bottom triangle
      if (onlyBottomTriangle && ((n-j) >= (m-i) + onlyBottomTriangleShift))
        {
          best.type = 'k';
          B[i][j] = best;
          if (cost.affine_gap)
            {
              Bins[i][j] = best;
              Bdel[i][j] = best;
            }
          continue ;
        }

306
      // The edit operations, with their backtracking information and their score
307 308

      // Match, mismatch
309
      try_operation(best, SUBST, i-1, j-1, B[i-1][j-1].score + cost.substitution(x[i-1], y[j-1]));
310

311 312 313 314 315 316 317 318 319 320 321 322
      if (!cost.affine_gap)
        {
          // Regular indels
          try_operation(best, INSER, i-1, j  , B[i-1][j  ].score + cost.insertion);
          try_operation(best, DELET, i  , j-1, B[i  ][j-1].score + cost.deletion);
        }
      else
        {
          // Gotoh affine gaps - insertion
          Bins[i][j].score = MINUS_INF ;
          try_operation(Bins[i][j], 'o', i-1, j, B[i-1][j].score + cost.open_insertion);
          try_operation(Bins[i][j], 'x', Bins[i-1][j].i, j, Bins[i-1][j].score + cost.extend_insertion);
323
          try_operation(best, INSER_AFFINE, Bins[i][j].i, j ,Bins[i][j].score);
324 325 326 327

          // Gotoh affine gaps - deletions
          Bdel[i][j].score = MINUS_INF ;
          try_operation(Bdel[i][j], 'o', i, j-1, B[i][j-1].score + cost.open_deletion);
328
          try_operation(Bdel[i][j], 'x', i, Bdel[i][j-1].j, Bdel[i][j-1].score + cost.extend_deletion);
329
          try_operation(best, DELET_AFFINE, i, Bdel[i][j].j,  Bdel[i][j].score);
330 331 332
        }

      // Homopolymers
333 334
      try_operation(best, HOMO2X, i-2, j-1, i > 1 ? B[i-2][j-1].score + cost.homo2(x[i-2], x[i-1], y[j-1]) : MINUS_INF);
      try_operation(best, HOMO2Y, i-1, j-2, j > 1 ? B[i-1][j-2].score + cost.homo2(y[j-2], y[j-1], x[i-1]) : MINUS_INF);
335
      
336
      // Local alignment
337 338
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	try_operation(best, BEGIN, 0, 0, 0);
Mikaël Salson's avatar
Mikaël Salson committed
339

340 341
      if ((best.type == SUBST) && (x[i-1] != y[j-1]))
        best.type = MISMATCH ;
Mikaël Salson's avatar
Mikaël Salson committed
342

343 344
      // Fill the score and the backtrack information
      B[i][j] = best ;
Mikaël Salson's avatar
Mikaël Salson committed
345
      
346 347
      // cout << " " << i << "," << j << " " << best_op.type << " " <<  best_op.i << "," << best_op.j << " "  << best_op.score << " " << best << endl ;

Mikaël Salson's avatar
Mikaël Salson committed
348 349
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
350
	  int tbest = best.score ;
Mikaël Salson's avatar
Mikaël Salson committed
351 352 353 354 355 356 357 358 359 360 361

	  if (mode == LocalEndWithSomeDeletions)
	    tbest += cost.deletion_end*(n-j);   
     
	  if (tbest > best_score)
	    {
	      best_score = tbest ;
	      best_i = i ;
	      best_j = j ;
	    }
	    
362
	  if (best.score == 0){
363
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
364 365
	  }
	}
Mikaël Salson's avatar
Mikaël Salson committed
366

Mikaël Salson's avatar
Mikaël Salson committed
367 368
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
369
      	  if (best.score == 0){
370
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
371 372 373 374 375
	  }
	}
      
    }

376 377
  // End. Find best_i and best_j, put FIN keywords where the backtrack should stop
  
Mikaël Salson's avatar
Mikaël Salson committed
378 379
  if (mode == Local || mode == LocalEndWithSomeDeletions)
    {
380 381
      for (int j=0; j<=n; j++){
	B[0][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
382 383
      }
      for (int i=0; i<=m; i++){
384
	B[i][0].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
385 386 387 388 389 390 391 392
      }
      
    }
  if (mode == SemiGlobal)
    {
      best_i = m ;
      
      for (int j=1; j<=n; j++)
393
	if (B[m][j].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
394
	  {
395
	    best_score = B[m][j].score ;
Mikaël Salson's avatar
Mikaël Salson committed
396 397 398
	    best_j = j ;
	  }
      for (int i=0; i<=m; i++){
399
	B[i][0].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
400 401 402 403 404 405 406 407
      }
    }

  if (mode == SemiGlobalTrans)
    {
      best_j = n ;
      
      for (int i=1; i<=m; i++)
408
	if (B[i][n].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
409
	  {
410
	    best_score = B[i][n].score ;
Mikaël Salson's avatar
Mikaël Salson committed
411 412 413 414
	    best_i = i ;
	  }
	  
      for (int i=0; i<=n; i++){
415
	B[0][i].type = FIN;
416
      }     
Mikaël Salson's avatar
Mikaël Salson committed
417 418
    }

Mathieu Giraud's avatar
Mathieu Giraud committed
419
  if ((mode == Global) || (mode == GlobalButMostlyLocal))
Mikaël Salson's avatar
Mikaël Salson committed
420 421 422
    {
      best_i = m ;
      best_j = n ;
423
      best_score = B[m][n].score;
Mikaël Salson's avatar
Mikaël Salson committed
424 425 426 427 428 429 430 431
    }

  if (reverse_x)
    best_i = m - best_i + 1 ;

  if (reverse_y)
    best_j = n - best_j + 1;

432 433
  B[0][0].type = FIN;
  
Mikaël Salson's avatar
Mikaël Salson committed
434 435 436 437 438 439 440
  // In the matrix positions start at 1, but start positions may start at 0
  best_i -= 1;   
  best_j -= 1;

  return best_score ;
}

441 442 443 444 445 446 447 448 449 450 451 452 453
int DynProg::best_score_on_i(int i, int *best_j)
{
  int best_score = MINUS_INF ;

  for (int j=0; j<=n; j++)
    if (B[i][j].score > best_score) {
      *best_j = j ;
      best_score = B[i][j].score;
    }

  return best_score ;
}

Mikaël Salson's avatar
Mikaël Salson committed
454 455
void DynProg::backtrack()
{
456
  // Tables for displaying the alignment
Mikaël Salson's avatar
Mikaël Salson committed
457 458 459 460
  gap1 = new int[x.size()+1];
  linkgap = new int[x.size()+1];
  gap2 = new int[y.size()+1];
      
461
  for (unsigned int i = 0; i <=x.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
462 463 464
    gap1[i] = 0;
    linkgap[i] = 0;
  }
465
  for (unsigned int i = 0; i <= y.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
466 467 468 469 470 471 472
    gap2[i] = 0;
  }
  
  
  int g1=x.size();
  int g2=y.size();
  
473 474
  int i=best_i+1;
  int j=best_j+1;
475

476 477 478 479 480 481 482
  // Retake good i/j when there were reversed strings
  if (reverse_x)
    i = m - i + 1 ;

  if (reverse_y)
    j = n - j + 1;

483 484 485 486 487
  if ((i<0) || (i>m) || (j<0) || (j>n))
    {
      cout << "Invalid backtrack starting point: " << i << "," << j << endl ;
      exit(1);
    }
Mikaël Salson's avatar
Mikaël Salson committed
488
  
489
  // Compute backtrack strings
490 491
  ostringstream back_x;
  ostringstream back_y;
Mikaël Salson's avatar
Mikaël Salson committed
492 493
  ostringstream back_tr;
  
494 495
  while (1) {

496 497 498

    if ((!reverse_y && (j == marked_pos_j))
        || (reverse_y && (n-j+1 == marked_pos_j)))
499
      {
500
        marked_pos_i = reverse_x ? m-i+1 : i ;
501 502
      }

503 504
    if  (B[i][j].type == FIN)
      break ;
Mikaël Salson's avatar
Mikaël Salson committed
505
      
506
    // cout << "bt " << i << "/" << j << " " << B[i][j].type << " " << B[i][j].i << "," << B[i][j].j << endl ;
507 508 509
    
    int next_i = B[i][j].i;
    int next_j = B[i][j].j;
Mikaël Salson's avatar
Mikaël Salson committed
510

511 512 513 514 515 516 517 518 519
    if ((next_i<0) || (next_i>m) || (next_j<0) || (next_j>n))
      {
        cout << "Invalid backtrack: " << i << "," << j << " -> " << next_i << "," << next_j << endl ;
        cout << back_x.str() << endl ;
        cout << back_tr.str() << endl ;
        cout << back_y.str() << endl ;
        exit(1);
      }

520 521
    // The maximum number of characters implied by this operation
    int max_k = max(i - next_i, j - next_j);
Mikaël Salson's avatar
Mikaël Salson committed
522
      
523 524 525
    for (int k=0; k<max_k; k++)
      {
        linkgap[g1]=g2;
526 527 528

        // Some character(s) from x, then fill with spaces
        if (i-k > next_i)
529
          {
530
            back_x << x[i-1-k];
531 532 533 534
            g1--;
          }
        else
          {
535
            back_x << " " ;
536 537
            gap1[g1]++ ;
          }
538 539 540

        // Some character(s) from y, then fill with spaces
        if (j-k > next_j)
541
          {
542
            back_y << y[j-1-k];
543 544 545 546
            g2--;
          }
        else
          {
547
            back_y << " " ;
548 549
            gap2[g2]++ ;
          }
550 551 552 553 554 555

        // The operation character, then fill with spaces
        if (k == max_k-1)          
          back_tr << B[i][j].type ;
        else
          back_tr << " " ;
556
      }
Mikaël Salson's avatar
Mikaël Salson committed
557
    
558 559
    i = next_i;
    j = next_j;
Mikaël Salson's avatar
Mikaël Salson committed
560 561
  }

562 563 564
  // Format backtrack strings
  
  string str1, str2, str3;
565
  str1 = back_x.str();
Mikaël Salson's avatar
Mikaël Salson committed
566 567 568 569
  str1 =string (str1.rbegin(), str1.rend());
  str1 = str1;
  str2=back_tr.str();
  str2 = string (str2.rbegin(), str2.rend());
570
  str3 = back_y.str();
Mikaël Salson's avatar
Mikaël Salson committed
571
  str3 = string (str3.rbegin(), str3.rend());
572 573

  ostringstream back;    
574 575 576
  back << setw(3) << i   << " " << str1.substr(0, BACKSIZE-8) << endl;
  back << setw(3) << " " << " " << str2.substr(0, BACKSIZE-8) << endl;
  back << setw(3) << j   << " " << str3.substr(0, BACKSIZE-8) << endl << endl;
Mikaël Salson's avatar
Mikaël Salson committed
577 578 579 580 581
  for (size_t k=0 ; (BACKSIZE-8+k*BACKSIZE)< str1.length() ; k++){
    back << str1.substr(BACKSIZE-8+k*BACKSIZE, BACKSIZE) << endl;
    back << str2.substr(BACKSIZE-8+k*BACKSIZE, BACKSIZE) << endl;
    back << str3.substr(BACKSIZE-8+k*BACKSIZE, BACKSIZE) << endl << endl;
  }
582
  back << "score: " << best_score << endl;
Mikaël Salson's avatar
Mikaël Salson committed
583
  
584 585
  first_i=i;
  first_j=j;
Mikaël Salson's avatar
Mikaël Salson committed
586 587
  
  str_back=back.str();
Mathieu Giraud's avatar
Mathieu Giraud committed
588 589

  // cout << str_back << endl ;
Mikaël Salson's avatar
Mikaël Salson committed
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606
}


float identity_percent(int score)
{
  int match = (score + IdentityDirty.match/2) / IdentityDirty.match ;

  int mismatch = (score - match * IdentityDirty.match) / IdentityDirty.mismatch ;

  float pcent = (100. * (float) match) / (match + mismatch) ;

  // cout << score << "  " << match << "  " << mismatch << "  " << pcent << endl ; 

  return pcent ;
}


607
ostream& to_ostream(ostream& out, string label, operation** B, const string x, const string y, int m, int n)
Mikaël Salson's avatar
Mikaël Salson committed
608
{
609
  out << label << "  " ;
Mikaël Salson's avatar
Mikaël Salson committed
610

611 612
  for (int j=0; j<n; j++)
    out << setw(4) << y[j] << " ";
Mikaël Salson's avatar
Mikaël Salson committed
613 614 615

  out << endl ;

616
  for (int i=0; i<=m; i++)
Mikaël Salson's avatar
Mikaël Salson committed
617 618
    {
      if (i)
619
	out << x[i-1] << " " ;
Mikaël Salson's avatar
Mikaël Salson committed
620 621 622
      else
	out << "  " ;

623
      for (int j=0; j<=n; j++)
624
        {
625
          if (B[i][j].score)
626
            {
627
              if (B[i][j].score <= MINUS_INF)
628 629
                out << "-INF" ;
              else
630
                out << setw(4) << B[i][j].score ;
631
            }
632
          else
633 634
            out << "    ";

635
          out << B[i][j].type ;
636
        }
Mikaël Salson's avatar
Mikaël Salson committed
637 638
      out << endl ;      
    }
639 640

  out << endl;
Mikaël Salson's avatar
Mikaël Salson committed
641 642 643 644
  
  return out ;
}

645 646
ostream& operator<<(ostream& out, const DynProg& dp)
{
647 648 649 650 651 652
  if (dp.cost.affine_gap)
    {
      to_ostream(out, "[Bins]", dp.Bins, dp.x, dp.y, dp.m, dp.n);
      to_ostream(out, "[Bdel]", dp.Bdel, dp.x, dp.y, dp.m, dp.n);
    }
  to_ostream(out, "[B]   ",dp.B, dp.x, dp.y, dp.m, dp.n);
653 654 655 656
  out << "best: " << dp.best_i << "," << dp.best_j << endl;

  return out;
}
Mikaël Salson's avatar
Mikaël Salson committed
657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696

Cost strToCost(string str, Cost default_cost){

  if (str.length()!=0){
    string::size_type stTemp = str.find(',');
	  
    std::list<int> val;
    
    while(stTemp != string::npos)
    {
	  val.push_back(atoi(str.substr(0, stTemp).c_str() ) );
	  str = str.substr(stTemp + 1);
	  stTemp = str.find(',');
    }
    val.push_back(atoi(str.c_str() ) );

    if (val.size()==5){ 
      int val1=val.front();
      val.pop_front();
      int val2=val.front();
      val.pop_front();
      int val3=val.front();
      val.pop_front();
      int val4=val.front();
      val.pop_front();
      int val5=val.front();
      val.pop_front();
      Cost result = Cost(val1, val2, val3, val4, val5);
      cout << "use custom Cost "<< result << endl;
      return result;
    }
    
    cout << "incorrect Cost format, use default "<< default_cost<<endl;
  }
  
  return default_cost;
}