dynprog.cpp 16.1 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 226 227 228 229 230 231
      if (cost.affine_gap) {
        Bins[i][0].score = cost.open_insertion + cost.extend_insertion * i;
        Bdel[i][0].score = MINUS_INF;
      }

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

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

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

245

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

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

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

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

283
      // The edit operations, with their backtracking information and their score
284 285

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

288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
      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);
305
          try_operation(Bdel[i][j], 'x', i, Bdel[i][j-1].j, Bdel[i][j-1].score + cost.extend_deletion);
306 307 308 309
          try_operation(best, DELET, i, Bdel[i][j].j,  Bdel[i][j].score);
        }

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

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

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

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

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

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

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

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

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

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

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

  return best_score ;
}

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

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

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

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

473 474 475

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

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

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

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

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

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

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

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

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

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


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


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

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

  out << endl ;

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

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

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

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

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

  return out;
}
Mikaël Salson's avatar
Mikaël Salson committed
634 635 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

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