dynprog.cpp 16.2 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
#define INIT '_'
38 39 40 41
#define HOMO2X 'h'
#define HOMO2Y 'h'
#define FIN ' '
#define BEGIN 'B'
Mikaël Salson's avatar
Mikaël Salson committed
42
#define BACKSIZE 120
Mikaël Salson's avatar
Mikaël Salson committed
43

44 45
#define BAD_BACKTRACK -1

Mikaël Salson's avatar
Mikaël Salson committed
46 47 48 49 50 51 52 53 54 55 56
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);
57

58 59 60
  this -> open_insertion = this -> open_deletion = MINUS_INF ;
  this -> extend_insertion = this -> extend_deletion = MINUS_INF ;
  this -> affine_gap = false ;
61 62

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

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 ;
76
  this -> affine_gap = true ;
77 78 79 80 81 82 83 84 85

  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
86 87 88 89 90 91 92 93 94
}


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

Mikaël Salson's avatar
Mikaël Salson committed
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
  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 ;
}

119
double Cost::toPValue(const int score)
120
{
121
  return K * exp(-lambda * score);
122 123
}

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

125 126 127
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
128 129 130 131 132 133 134
{
  this -> x = reverse_x ? reverse(x) : x ;
  this -> y = reverse_y ? reverse(y) : y ;

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

135 136 137
  this -> marked_pos_j = marked_pos_j;
  this -> marked_pos_i = 0;

Mikaël Salson's avatar
Mikaël Salson committed
138 139 140 141
  m = x.size();
  n = y.size();
  // cout << "Dynprog " << x << "(" << m << ") / " << y << "(" << n << ")" << endl ;

142
  B = new operation*[m+1];
Mikaël Salson's avatar
Mikaël Salson committed
143
  for (int i = 0; i <= m; i++) {
144
    B[i] = new operation[n+1];
Mikaël Salson's avatar
Mikaël Salson committed
145
  }
146 147 148 149 150 151 152 153
  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
154 155 156 157 158 159 160
  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
161 162 163 164
  
  gap1=NULL;
  gap2=NULL;
  linkgap=NULL;
Mikaël Salson's avatar
Mikaël Salson committed
165 166 167 168 169

  init();
}

DynProg::~DynProg() {
170

Mikaël Salson's avatar
Mikaël Salson committed
171 172
  for (int i = 0; i <= m; i++) {
    delete [] B[i];
173 174
    delete [] Bins[i];
    delete [] Bdel[i];
Mikaël Salson's avatar
Mikaël Salson committed
175
  }
176
  delete [] B;  
177 178
  delete [] Bins;
  delete [] Bdel;  
Mikaël Salson's avatar
Mikaël Salson committed
179 180 181
  delete [] gap1;
  delete [] gap2;
  delete [] linkgap;
Mikaël Salson's avatar
Mikaël Salson committed
182 183 184
}

void DynProg::init()
185
{
186 187 188 189 190 191 192 193 194 195 196 197 198 199
  // 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;
200 201
          Bins[i][j].type = INIT;
          Bdel[i][j].type = INIT;
202 203
        }
      }
204

205
  // Initialize backtrack labels (B[0][0] labels are never used)
206 207 208 209 210 211 212 213 214 215 216 217 218 219
  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 ;
    }

220
  if (!(mode == Local || mode == LocalEndWithSomeDeletions)) // Global, SemiGlobal
221
    for (int i=0; i<=m; i++) {
222
      B[i][0].score = i * cost.insertion ;
223 224 225
      if (cost.affine_gap) {
        Bins[i][0].score = cost.open_insertion + cost.extend_insertion * i;
        Bdel[i][0].score = MINUS_INF;
226
        B[i][0].score = Bins[i][0].score ;
227 228 229 230 231 232
      }

      if (mode == GlobalButMostlyLocal)
        B[i][0].score /= 2;
    }

233
  if (!(mode == SemiGlobal || mode == SemiGlobalTrans || mode == Local || mode == LocalEndWithSomeDeletions)) // Global
234
    for (int j=0; j<=n; j++) {
235
      B[0][j].score = j * cost.deletion ;
236
       if (cost.affine_gap) {
237
         if (j) Bins[0][j].score = MINUS_INF;
238
         Bdel[0][j].score = cost.open_deletion + cost.extend_deletion * j;
239
         B[0][j].score = Bdel[0][j].score;
240 241 242 243 244
       }

       if (mode == GlobalButMostlyLocal)
         B[0][j].score /= 2;
    }
Mikaël Salson's avatar
Mikaël Salson committed
245 246
}

247

248 249 250 251 252 253 254 255 256 257 258
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 ;
    }
}

259
int DynProg::compute(bool onlyBottomTriangle, int onlyBottomTriangleShift)
Mikaël Salson's avatar
Mikaël Salson committed
260 261
{
  best_score = MINUS_INF ;
262 263 264 265
  best_i = 0 ;
  best_j = 0 ;
  
  operation best ;
Mikaël Salson's avatar
Mikaël Salson committed
266 267 268 269

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

272 273 274 275 276 277 278 279 280 281 282 283 284
      // 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 ;
        }

285
      // The edit operations, with their backtracking information and their score
286 287

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

290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
      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);
          try_operation(best, INSER, Bins[i][j].i, j ,Bins[i][j].score);

          // 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);
307
          try_operation(Bdel[i][j], 'x', i, Bdel[i][j-1].j, Bdel[i][j-1].score + cost.extend_deletion);
308 309 310 311
          try_operation(best, DELET, i, Bdel[i][j].j,  Bdel[i][j].score);
        }

      // Homopolymers
312 313
      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);
314
      
315
      // Local alignment
316 317
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	try_operation(best, BEGIN, 0, 0, 0);
Mikaël Salson's avatar
Mikaël Salson committed
318

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

322 323
      // Fill the score and the backtrack information
      B[i][j] = best ;
Mikaël Salson's avatar
Mikaël Salson committed
324
      
325 326
      // 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
327 328
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
329
	  int tbest = best.score ;
Mikaël Salson's avatar
Mikaël Salson committed
330 331 332 333 334 335 336 337 338 339 340

	  if (mode == LocalEndWithSomeDeletions)
	    tbest += cost.deletion_end*(n-j);   
     
	  if (tbest > best_score)
	    {
	      best_score = tbest ;
	      best_i = i ;
	      best_j = j ;
	    }
	    
341
	  if (best.score == 0){
342
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
343 344
	  }
	}
Mikaël Salson's avatar
Mikaël Salson committed
345

Mikaël Salson's avatar
Mikaël Salson committed
346 347
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
348
      	  if (best.score == 0){
349
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
350 351 352 353 354
	  }
	}
      
    }

355 356
  // End. Find best_i and best_j, put FIN keywords where the backtrack should stop
  
Mikaël Salson's avatar
Mikaël Salson committed
357 358
  if (mode == Local || mode == LocalEndWithSomeDeletions)
    {
359 360
      for (int j=0; j<=n; j++){
	B[0][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
361 362
      }
      for (int i=0; i<=m; i++){
363
	B[i][0].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
364 365 366 367 368 369 370 371
      }
      
    }
  if (mode == SemiGlobal)
    {
      best_i = m ;
      
      for (int j=1; j<=n; j++)
372
	if (B[m][j].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
373
	  {
374
	    best_score = B[m][j].score ;
Mikaël Salson's avatar
Mikaël Salson committed
375 376 377
	    best_j = j ;
	  }
      for (int i=0; i<=m; i++){
378
	B[i][0].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
379 380 381 382 383 384 385 386
      }
    }

  if (mode == SemiGlobalTrans)
    {
      best_j = n ;
      
      for (int i=1; i<=m; i++)
387
	if (B[i][n].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
388
	  {
389
	    best_score = B[i][n].score ;
Mikaël Salson's avatar
Mikaël Salson committed
390 391 392 393
	    best_i = i ;
	  }
	  
      for (int i=0; i<=n; i++){
394
	B[0][i].type = FIN;
395
      }     
Mikaël Salson's avatar
Mikaël Salson committed
396 397
    }

Mathieu Giraud's avatar
Mathieu Giraud committed
398
  if ((mode == Global) || (mode == GlobalButMostlyLocal))
Mikaël Salson's avatar
Mikaël Salson committed
399 400 401
    {
      best_i = m ;
      best_j = n ;
402
      best_score = B[m][n].score;
Mikaël Salson's avatar
Mikaël Salson committed
403 404 405 406 407 408 409 410
    }

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

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

411 412
  B[0][0].type = FIN;
  
Mikaël Salson's avatar
Mikaël Salson committed
413 414 415 416 417 418 419
  // In the matrix positions start at 1, but start positions may start at 0
  best_i -= 1;   
  best_j -= 1;

  return best_score ;
}

420 421 422 423 424 425 426 427 428 429 430 431 432
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
433 434
void DynProg::backtrack()
{
435
  // Tables for displaying the alignment
Mikaël Salson's avatar
Mikaël Salson committed
436 437 438 439
  gap1 = new int[x.size()+1];
  linkgap = new int[x.size()+1];
  gap2 = new int[y.size()+1];
      
440
  for (unsigned int i = 0; i <=x.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
441 442 443
    gap1[i] = 0;
    linkgap[i] = 0;
  }
444
  for (unsigned int i = 0; i <= y.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
445 446 447 448 449 450 451
    gap2[i] = 0;
  }
  
  
  int g1=x.size();
  int g2=y.size();
  
452 453
  int i=best_i+1;
  int j=best_j+1;
454

455 456 457 458 459 460 461
  // Retake good i/j when there were reversed strings
  if (reverse_x)
    i = m - i + 1 ;

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

462 463 464 465 466
  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
467
  
468
  // Compute backtrack strings
469 470
  ostringstream back_x;
  ostringstream back_y;
Mikaël Salson's avatar
Mikaël Salson committed
471 472
  ostringstream back_tr;
  
473 474
  while (1) {

475 476 477

    if ((!reverse_y && (j == marked_pos_j))
        || (reverse_y && (n-j+1 == marked_pos_j)))
478
      {
479
        marked_pos_i = reverse_x ? m-i+1 : i ;
480 481
      }

482 483
    if  (B[i][j].type == FIN)
      break ;
Mikaël Salson's avatar
Mikaël Salson committed
484
      
485
    // cout << "bt " << i << "/" << j << " " << B[i][j].type << " " << B[i][j].i << "," << B[i][j].j << endl ;
486 487 488
    
    int next_i = B[i][j].i;
    int next_j = B[i][j].j;
Mikaël Salson's avatar
Mikaël Salson committed
489

490 491 492 493 494 495 496 497 498
    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);
      }

499 500
    // 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
501
      
502 503 504
    for (int k=0; k<max_k; k++)
      {
        linkgap[g1]=g2;
505 506 507

        // Some character(s) from x, then fill with spaces
        if (i-k > next_i)
508
          {
509
            back_x << x[i-1-k];
510 511 512 513
            g1--;
          }
        else
          {
514
            back_x << " " ;
515 516
            gap1[g1]++ ;
          }
517 518 519

        // Some character(s) from y, then fill with spaces
        if (j-k > next_j)
520
          {
521
            back_y << y[j-1-k];
522 523 524 525
            g2--;
          }
        else
          {
526
            back_y << " " ;
527 528
            gap2[g2]++ ;
          }
529 530 531 532 533 534

        // The operation character, then fill with spaces
        if (k == max_k-1)          
          back_tr << B[i][j].type ;
        else
          back_tr << " " ;
535
      }
Mikaël Salson's avatar
Mikaël Salson committed
536
    
537 538
    i = next_i;
    j = next_j;
Mikaël Salson's avatar
Mikaël Salson committed
539 540
  }

541 542 543
  // Format backtrack strings
  
  string str1, str2, str3;
544
  str1 = back_x.str();
Mikaël Salson's avatar
Mikaël Salson committed
545 546 547 548
  str1 =string (str1.rbegin(), str1.rend());
  str1 = str1;
  str2=back_tr.str();
  str2 = string (str2.rbegin(), str2.rend());
549
  str3 = back_y.str();
Mikaël Salson's avatar
Mikaël Salson committed
550
  str3 = string (str3.rbegin(), str3.rend());
551 552

  ostringstream back;    
553 554 555
  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
556 557 558 559 560
  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;
  }
561
  back << "score: " << best_score << endl;
Mikaël Salson's avatar
Mikaël Salson committed
562
  
563 564
  first_i=i;
  first_j=j;
Mikaël Salson's avatar
Mikaël Salson committed
565 566
  
  str_back=back.str();
Mathieu Giraud's avatar
Mathieu Giraud committed
567 568

  // cout << str_back << endl ;
Mikaël Salson's avatar
Mikaël Salson committed
569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585
}


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 ;
}


586
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
587
{
588
  out << label << "  " ;
Mikaël Salson's avatar
Mikaël Salson committed
589

590 591
  for (int j=0; j<n; j++)
    out << setw(4) << y[j] << " ";
Mikaël Salson's avatar
Mikaël Salson committed
592 593 594

  out << endl ;

595
  for (int i=0; i<=m; i++)
Mikaël Salson's avatar
Mikaël Salson committed
596 597
    {
      if (i)
598
	out << x[i-1] << " " ;
Mikaël Salson's avatar
Mikaël Salson committed
599 600 601
      else
	out << "  " ;

602
      for (int j=0; j<=n; j++)
603
        {
604
          if (B[i][j].score)
605
            {
606
              if (B[i][j].score <= MINUS_INF)
607 608
                out << "-INF" ;
              else
609
                out << setw(4) << B[i][j].score ;
610
            }
611
          else
612 613
            out << "    ";

614
          out << B[i][j].type ;
615
        }
Mikaël Salson's avatar
Mikaël Salson committed
616 617
      out << endl ;      
    }
618 619

  out << endl;
Mikaël Salson's avatar
Mikaël Salson committed
620 621 622 623
  
  return out ;
}

624 625
ostream& operator<<(ostream& out, const DynProg& dp)
{
626 627 628 629 630 631
  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);
632 633 634 635
  out << "best: " << dp.best_i << "," << dp.best_j << endl;

  return out;
}
Mikaël Salson's avatar
Mikaël Salson committed
636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675

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;
}