dynprog.cpp 16 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 37 38 39 40
#define SUBST '|'
#define MISMATCH '.'
#define INSER 'i'
#define DELET 'd'
#define HOMO2X 'h'
#define HOMO2Y 'h'
#define FIN ' '
#define BEGIN 'B'
Mikaël Salson's avatar
Mikaël Salson committed
41
#define BACKSIZE 120
Mikaël Salson's avatar
Mikaël Salson committed
42

43 44
#define BAD_BACKTRACK -1

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

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

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

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

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


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

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

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

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

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

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

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

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

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

  init();
}

DynProg::~DynProg() {
169

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

void DynProg::init()
184
{
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
  // 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;
        }
      }
201

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

217
  if (!(mode == Local || mode == LocalEndWithSomeDeletions)) // Global, SemiGlobal
218
    for (int i=0; i<=m; i++) {
219
      B[i][0].score = i * cost.insertion ;
220 221 222 223 224 225 226 227 228
      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;
    }

229
  if (!(mode == SemiGlobal || mode == SemiGlobalTrans || mode == Local || mode == LocalEndWithSomeDeletions)) // Global
230
    for (int j=0; j<=n; j++) {
231
      B[0][j].score = j * cost.deletion ;
232 233 234 235 236 237 238 239
       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
240 241
}

242

243 244 245 246 247 248 249 250 251 252 253
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 ;
    }
}

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

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

267 268 269 270 271 272 273 274 275 276 277 278 279
      // 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 ;
        }

280
      // The edit operations, with their backtracking information and their score
281 282

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

285 286 287 288 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);
          try_operation(Bdel[i][j], 'x', i, Bdel[i][j-1].j, Bdel[i-1][j].score + cost.extend_deletion);
          try_operation(best, DELET, i, Bdel[i][j].j,  Bdel[i][j].score);
        }

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

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

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

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

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

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

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

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

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

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

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

  return best_score ;
}

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

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

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

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

470 471 472

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

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

485 486 487 488 489 490 491 492 493
    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);
      }

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

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

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

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

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

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

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


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


581
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
582
{
583
  out << label << "  " ;
Mikaël Salson's avatar
Mikaël Salson committed
584

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

  out << endl ;

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

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

609
          out << B[i][j].type ;
610
        }
Mikaël Salson's avatar
Mikaël Salson committed
611 612
      out << endl ;      
    }
613 614

  out << endl;
Mikaël Salson's avatar
Mikaël Salson committed
615 616 617 618
  
  return out ;
}

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

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

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