dynprog.cpp 17 KB
Newer Older
Mikaël Salson's avatar
Mikaël Salson committed
1
/*
2
  This file is part of Vidjil <http://www.vidjil.org>
3
  Copyright (C) 2011-2017 by Bonsai bioinformatics
4 5 6 7 8
  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
}


ostream& operator<<(ostream& out, const Cost& cost)
{
93
  if (cost.debug)
Mikaël Salson's avatar
Mikaël Salson committed
94 95 96 97
  out << "(" << cost.match 
      << ", " << cost.mismatch
      << "/" << cost.insertion
      << "/" << cost.deletion
98 99
      << ", " << cost.open_insertion << cost.extend_insertion
      << "/" << cost.open_deletion << cost.extend_deletion
Mikaël Salson's avatar
Mikaël Salson committed
100 101
      << ", " << cost.deletion_end
      << ", " << cost.homopolymer
102 103
      << ") "
      << cost.K << "/" << cost.lambda ;
104 105 106 107 108 109 110
  else
    out << "\"" << cost.match
        << ", " << cost.mismatch
        << ", " << cost.insertion
        << ", " << cost.homopolymer
        << ", " << cost.deletion_end
        << "\"" ;
111

Mikaël Salson's avatar
Mikaël Salson committed
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
  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 ;
}

129
double Cost::toPValue(const int score)
130
{
131 132 133 134 135 136 137
#ifdef DEBUG_EVALUE
  cout << *this
       << " | score " << score
       << ", matches " << (float) score / (float) match
       << " ==> pvalue " << std::scientific << K * exp(-lambda * score) << std::fixed
       << endl;
#endif
138
  return K * exp(-lambda * score);
139 140
}

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

142 143 144
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
145 146 147 148 149 150 151
{
  this -> x = reverse_x ? reverse(x) : x ;
  this -> y = reverse_y ? reverse(y) : y ;

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

152 153 154
  this -> marked_pos_j = marked_pos_j;
  this -> marked_pos_i = 0;

Mikaël Salson's avatar
Mikaël Salson committed
155 156 157 158
  m = x.size();
  n = y.size();
  // cout << "Dynprog " << x << "(" << m << ") / " << y << "(" << n << ")" << endl ;

159
  B = new operation*[m+1];
Mikaël Salson's avatar
Mikaël Salson committed
160
  for (int i = 0; i <= m; i++) {
161
    B[i] = new operation[n+1];
Mikaël Salson's avatar
Mikaël Salson committed
162
  }
163 164 165 166 167 168 169 170
  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
171 172 173 174 175 176 177
  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
178 179 180 181
  
  gap1=NULL;
  gap2=NULL;
  linkgap=NULL;
Mikaël Salson's avatar
Mikaël Salson committed
182 183 184 185 186

  init();
}

DynProg::~DynProg() {
187

Mikaël Salson's avatar
Mikaël Salson committed
188 189
  for (int i = 0; i <= m; i++) {
    delete [] B[i];
190 191
    delete [] Bins[i];
    delete [] Bdel[i];
Mikaël Salson's avatar
Mikaël Salson committed
192
  }
193
  delete [] B;  
194 195
  delete [] Bins;
  delete [] Bdel;  
Mikaël Salson's avatar
Mikaël Salson committed
196 197 198
  delete [] gap1;
  delete [] gap2;
  delete [] linkgap;
Mikaël Salson's avatar
Mikaël Salson committed
199 200 201
}

void DynProg::init()
202
{
203 204 205 206 207 208 209 210 211 212 213 214 215 216
  // 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;
217 218
          Bins[i][j].type = INIT;
          Bdel[i][j].type = INIT;
219 220
        }
      }
221

222
  // Initialize backtrack labels (B[0][0] labels are never used)
223 224 225 226 227 228 229 230 231 232 233 234 235 236
  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 ;
    }

237
  if (!(mode == Local || mode == LocalEndWithSomeDeletions)) // Global, SemiGlobal
238
    for (int i=0; i<=m; i++) {
239
      B[i][0].score = i * cost.insertion ;
240 241 242
      if (cost.affine_gap) {
        Bins[i][0].score = cost.open_insertion + cost.extend_insertion * i;
        Bdel[i][0].score = MINUS_INF;
243
        B[i][0].score = Bins[i][0].score ;
244
        B[i][0].type = INSER_AFFINE;
245 246
      }

247
      if (mode == GlobalButMostlyLocal) {
248
        B[i][0].score /= 2;
249 250 251 252 253
        if (cost.affine_gap) {
          Bins[i][0].score -= cost.open_insertion / 2 ;
          B[i][0].score = Bins[i][0].score ;
        }
      }
254 255
    }

256
  if (!(mode == SemiGlobal || mode == SemiGlobalTrans || mode == Local || mode == LocalEndWithSomeDeletions)) // Global
257
    for (int j=0; j<=n; j++) {
258
      B[0][j].score = j * cost.deletion ;
259
       if (cost.affine_gap) {
260
         if (j) Bins[0][j].score = MINUS_INF;
261
         Bdel[0][j].score = cost.open_deletion + cost.extend_deletion * j;
262
         B[0][j].score = Bdel[0][j].score;
263
         B[0][j].type = DELET_AFFINE;
264 265
       }

266
       if (mode == GlobalButMostlyLocal) {
267
         B[0][j].score /= 2;
268 269 270 271 272
         if (cost.affine_gap) {
           Bdel[0][j].score -= cost.open_deletion / 2 ;
           B[0][j].score = Bdel[0][j].score;
         }
       }
273
    }
Mikaël Salson's avatar
Mikaël Salson committed
274 275
}

276

277 278 279 280 281 282 283 284 285 286 287
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 ;
    }
}

288
int DynProg::compute(bool onlyBottomTriangle, int onlyBottomTriangleShift)
Mikaël Salson's avatar
Mikaël Salson committed
289 290
{
  best_score = MINUS_INF ;
291 292 293 294
  best_i = 0 ;
  best_j = 0 ;
  
  operation best ;
Mikaël Salson's avatar
Mikaël Salson committed
295 296 297 298

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

301 302 303 304 305 306 307 308 309 310 311 312 313
      // 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 ;
        }

314
      // The edit operations, with their backtracking information and their score
315 316

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

319 320 321 322 323 324 325 326 327 328 329 330
      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);
331
          try_operation(best, INSER_AFFINE, Bins[i][j].i, j ,Bins[i][j].score);
332 333 334 335

          // 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);
336
          try_operation(Bdel[i][j], 'x', i, Bdel[i][j-1].j, Bdel[i][j-1].score + cost.extend_deletion);
337
          try_operation(best, DELET_AFFINE, i, Bdel[i][j].j,  Bdel[i][j].score);
338 339 340
        }

      // Homopolymers
341 342
      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);
343
      
344
      // Local alignment
345 346
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	try_operation(best, BEGIN, 0, 0, 0);
Mikaël Salson's avatar
Mikaël Salson committed
347

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

351 352
      // Fill the score and the backtrack information
      B[i][j] = best ;
Mikaël Salson's avatar
Mikaël Salson committed
353
      
354 355
      // 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
356 357
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
358
	  int tbest = best.score ;
Mikaël Salson's avatar
Mikaël Salson committed
359 360 361 362 363 364 365 366 367 368 369

	  if (mode == LocalEndWithSomeDeletions)
	    tbest += cost.deletion_end*(n-j);   
     
	  if (tbest > best_score)
	    {
	      best_score = tbest ;
	      best_i = i ;
	      best_j = j ;
	    }
	    
370
	  if (best.score == 0){
371
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
372 373
	  }
	}
Mikaël Salson's avatar
Mikaël Salson committed
374

Mikaël Salson's avatar
Mikaël Salson committed
375 376
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
377
      	  if (best.score == 0){
378
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
379 380 381 382 383
	  }
	}
      
    }

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

  if (mode == SemiGlobalTrans)
    {
      best_j = n ;
      
      for (int i=1; i<=m; i++)
416
	if (B[i][n].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
417
	  {
418
	    best_score = B[i][n].score ;
Mikaël Salson's avatar
Mikaël Salson committed
419 420 421 422
	    best_i = i ;
	  }
	  
      for (int i=0; i<=n; i++){
423
	B[0][i].type = FIN;
424
      }     
Mikaël Salson's avatar
Mikaël Salson committed
425 426
    }

Mathieu Giraud's avatar
Mathieu Giraud committed
427
  if ((mode == Global) || (mode == GlobalButMostlyLocal))
Mikaël Salson's avatar
Mikaël Salson committed
428 429 430
    {
      best_i = m ;
      best_j = n ;
431
      best_score = B[m][n].score;
Mikaël Salson's avatar
Mikaël Salson committed
432 433 434 435 436 437 438 439
    }

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

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

440 441
  B[0][0].type = FIN;
  
Mikaël Salson's avatar
Mikaël Salson committed
442 443 444 445 446 447 448
  // In the matrix positions start at 1, but start positions may start at 0
  best_i -= 1;   
  best_j -= 1;

  return best_score ;
}

449 450 451 452 453 454 455 456 457 458 459 460 461
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
462 463
void DynProg::backtrack()
{
464
  // Tables for displaying the alignment
Mikaël Salson's avatar
Mikaël Salson committed
465 466 467 468
  gap1 = new int[x.size()+1];
  linkgap = new int[x.size()+1];
  gap2 = new int[y.size()+1];
      
469
  for (unsigned int i = 0; i <=x.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
470 471 472
    gap1[i] = 0;
    linkgap[i] = 0;
  }
473
  for (unsigned int i = 0; i <= y.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
474 475 476 477 478 479 480
    gap2[i] = 0;
  }
  
  
  int g1=x.size();
  int g2=y.size();
  
481 482
  int i=best_i+1;
  int j=best_j+1;
483

484 485 486 487 488 489 490
  // Retake good i/j when there were reversed strings
  if (reverse_x)
    i = m - i + 1 ;

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

491 492 493 494 495
  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
496
  
497
  // Compute backtrack strings
498 499
  ostringstream back_x;
  ostringstream back_y;
Mikaël Salson's avatar
Mikaël Salson committed
500 501
  ostringstream back_tr;
  
502 503
  while (1) {

504 505 506

    if ((!reverse_y && (j == marked_pos_j))
        || (reverse_y && (n-j+1 == marked_pos_j)))
507
      {
508
        marked_pos_i = reverse_x ? m-i+1 : i ;
509 510
      }

511 512
    if  (B[i][j].type == FIN)
      break ;
Mikaël Salson's avatar
Mikaël Salson committed
513
      
514
    // cout << "bt " << i << "/" << j << " " << B[i][j].type << " " << B[i][j].i << "," << B[i][j].j << endl ;
515 516 517
    
    int next_i = B[i][j].i;
    int next_j = B[i][j].j;
Mikaël Salson's avatar
Mikaël Salson committed
518

519 520 521 522 523 524 525 526 527
    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);
      }

528 529
    // 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
530
      
531 532 533
    for (int k=0; k<max_k; k++)
      {
        linkgap[g1]=g2;
534 535 536

        // Some character(s) from x, then fill with spaces
        if (i-k > next_i)
537
          {
538
            back_x << x[i-1-k];
539 540 541 542
            g1--;
          }
        else
          {
543
            back_x << " " ;
544 545
            gap1[g1]++ ;
          }
546 547 548

        // Some character(s) from y, then fill with spaces
        if (j-k > next_j)
549
          {
550
            back_y << y[j-1-k];
551 552 553 554
            g2--;
          }
        else
          {
555
            back_y << " " ;
556 557
            gap2[g2]++ ;
          }
558 559 560 561 562 563

        // The operation character, then fill with spaces
        if (k == max_k-1)          
          back_tr << B[i][j].type ;
        else
          back_tr << " " ;
564
      }
Mikaël Salson's avatar
Mikaël Salson committed
565
    
566 567
    i = next_i;
    j = next_j;
Mikaël Salson's avatar
Mikaël Salson committed
568 569
  }

570 571 572
  // Format backtrack strings
  
  string str1, str2, str3;
573
  str1 = back_x.str();
Mikaël Salson's avatar
Mikaël Salson committed
574 575 576 577
  str1 =string (str1.rbegin(), str1.rend());
  str1 = str1;
  str2=back_tr.str();
  str2 = string (str2.rbegin(), str2.rend());
578
  str3 = back_y.str();
Mikaël Salson's avatar
Mikaël Salson committed
579
  str3 = string (str3.rbegin(), str3.rend());
580 581

  ostringstream back;    
582 583 584
  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
585 586 587 588 589
  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;
  }
590
  back << "score: " << best_score << endl;
Mikaël Salson's avatar
Mikaël Salson committed
591
  
592 593
  first_i=i;
  first_j=j;
Mikaël Salson's avatar
Mikaël Salson committed
594 595
  
  str_back=back.str();
Mathieu Giraud's avatar
Mathieu Giraud committed
596 597

  // cout << str_back << endl ;
Mikaël Salson's avatar
Mikaël Salson committed
598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614
}


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


615
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
616
{
617
  out << label << "  " ;
Mikaël Salson's avatar
Mikaël Salson committed
618

619 620
  for (int j=0; j<n; j++)
    out << setw(4) << y[j] << " ";
Mikaël Salson's avatar
Mikaël Salson committed
621 622 623

  out << endl ;

624
  for (int i=0; i<=m; i++)
Mikaël Salson's avatar
Mikaël Salson committed
625 626
    {
      if (i)
627
	out << x[i-1] << " " ;
Mikaël Salson's avatar
Mikaël Salson committed
628 629 630
      else
	out << "  " ;

631
      for (int j=0; j<=n; j++)
632
        {
633
          if (B[i][j].score)
634
            {
635
              if (B[i][j].score <= MINUS_INF)
636 637
                out << "-INF" ;
              else
638
                out << setw(4) << B[i][j].score ;
639
            }
640
          else
641 642
            out << "    ";

643
          out << B[i][j].type ;
644
        }
Mikaël Salson's avatar
Mikaël Salson committed
645 646
      out << endl ;      
    }
647 648

  out << endl;
Mikaël Salson's avatar
Mikaël Salson committed
649 650 651 652
  
  return out ;
}

653 654
ostream& operator<<(ostream& out, const DynProg& dp)
{
655 656 657 658 659 660
  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);
661 662 663 664
  out << "best: " << dp.best_i << "," << dp.best_j << endl;

  return out;
}
Mikaël Salson's avatar
Mikaël Salson committed
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 697 698 699 700 701 702 703 704

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