dynprog.cpp 17.1 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
  else
    out << "\"" << cost.match
        << ", " << cost.mismatch
        << ", " << cost.insertion
        << ", " << cost.deletion_end
109
        << ", " << cost.homopolymer
110
        << "\"" ;
111

Mikaël Salson's avatar
Mikaël Salson committed
112 113 114
  return out;
}

115 116 117 118 119 120 121 122 123
string string_of_cost(const Cost cost)
{
   stringstream ss;
   ss << cost ;
   return ss.str();
}



Mikaël Salson's avatar
Mikaël Salson committed
124 125 126 127 128 129 130 131 132 133 134 135 136 137
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 ;
}

138
double Cost::toPValue(const int score)
139
{
140 141 142 143 144 145 146
#ifdef DEBUG_EVALUE
  cout << *this
       << " | score " << score
       << ", matches " << (float) score / (float) match
       << " ==> pvalue " << std::scientific << K * exp(-lambda * score) << std::fixed
       << endl;
#endif
147
  return K * exp(-lambda * score);
148 149
}

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

151 152 153
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
154 155 156 157 158 159 160
{
  this -> x = reverse_x ? reverse(x) : x ;
  this -> y = reverse_y ? reverse(y) : y ;

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

161 162 163
  this -> marked_pos_j = marked_pos_j;
  this -> marked_pos_i = 0;

Mikaël Salson's avatar
Mikaël Salson committed
164 165 166 167
  m = x.size();
  n = y.size();
  // cout << "Dynprog " << x << "(" << m << ") / " << y << "(" << n << ")" << endl ;

168
  B = new operation*[m+1];
Mikaël Salson's avatar
Mikaël Salson committed
169
  for (int i = 0; i <= m; i++) {
170
    B[i] = new operation[n+1];
Mikaël Salson's avatar
Mikaël Salson committed
171
  }
172 173 174 175 176 177 178 179
  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
180 181 182 183 184 185 186
  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
187 188 189 190
  
  gap1=NULL;
  gap2=NULL;
  linkgap=NULL;
Mikaël Salson's avatar
Mikaël Salson committed
191 192 193 194 195

  init();
}

DynProg::~DynProg() {
196

Mikaël Salson's avatar
Mikaël Salson committed
197 198
  for (int i = 0; i <= m; i++) {
    delete [] B[i];
199 200
    delete [] Bins[i];
    delete [] Bdel[i];
Mikaël Salson's avatar
Mikaël Salson committed
201
  }
202
  delete [] B;  
203 204
  delete [] Bins;
  delete [] Bdel;  
Mikaël Salson's avatar
Mikaël Salson committed
205 206 207
  delete [] gap1;
  delete [] gap2;
  delete [] linkgap;
Mikaël Salson's avatar
Mikaël Salson committed
208 209 210
}

void DynProg::init()
211
{
212 213 214 215 216 217 218 219 220 221 222 223 224 225
  // 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;
226 227
          Bins[i][j].type = INIT;
          Bdel[i][j].type = INIT;
228 229
        }
      }
230

231
  // Initialize backtrack labels (B[0][0] labels are never used)
232 233 234 235 236 237 238 239 240 241 242 243 244 245
  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 ;
    }

246
  if (!(mode == Local || mode == LocalEndWithSomeDeletions)) // Global, SemiGlobal
247
    for (int i=0; i<=m; i++) {
248
      B[i][0].score = i * cost.insertion ;
249 250 251
      if (cost.affine_gap) {
        Bins[i][0].score = cost.open_insertion + cost.extend_insertion * i;
        Bdel[i][0].score = MINUS_INF;
252
        B[i][0].score = Bins[i][0].score ;
253
        B[i][0].type = INSER_AFFINE;
254 255
      }

256
      if (mode == GlobalButMostlyLocal) {
257
        B[i][0].score /= 2;
258 259 260 261 262
        if (cost.affine_gap) {
          Bins[i][0].score -= cost.open_insertion / 2 ;
          B[i][0].score = Bins[i][0].score ;
        }
      }
263 264
    }

265
  if (!(mode == SemiGlobal || mode == SemiGlobalTrans || mode == Local || mode == LocalEndWithSomeDeletions)) // Global
266
    for (int j=0; j<=n; j++) {
267
      B[0][j].score = j * cost.deletion ;
268
       if (cost.affine_gap) {
269
         if (j) Bins[0][j].score = MINUS_INF;
270
         Bdel[0][j].score = cost.open_deletion + cost.extend_deletion * j;
271
         B[0][j].score = Bdel[0][j].score;
272
         B[0][j].type = DELET_AFFINE;
273 274
       }

275
       if (mode == GlobalButMostlyLocal) {
276
         B[0][j].score /= 2;
277 278 279 280 281
         if (cost.affine_gap) {
           Bdel[0][j].score -= cost.open_deletion / 2 ;
           B[0][j].score = Bdel[0][j].score;
         }
       }
282
    }
Mikaël Salson's avatar
Mikaël Salson committed
283 284
}

285

286 287 288 289 290 291 292 293 294 295 296
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 ;
    }
}

297
int DynProg::compute(bool onlyBottomTriangle, int onlyBottomTriangleShift)
Mikaël Salson's avatar
Mikaël Salson committed
298 299
{
  best_score = MINUS_INF ;
300 301 302 303
  best_i = 0 ;
  best_j = 0 ;
  
  operation best ;
Mikaël Salson's avatar
Mikaël Salson committed
304 305 306 307

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

310 311 312 313 314 315 316 317 318 319 320 321 322
      // 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 ;
        }

323
      // The edit operations, with their backtracking information and their score
324 325

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

328 329 330 331 332 333 334 335 336 337 338 339
      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);
340
          try_operation(best, INSER_AFFINE, Bins[i][j].i, j ,Bins[i][j].score);
341 342 343 344

          // 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);
345
          try_operation(Bdel[i][j], 'x', i, Bdel[i][j-1].j, Bdel[i][j-1].score + cost.extend_deletion);
346
          try_operation(best, DELET_AFFINE, i, Bdel[i][j].j,  Bdel[i][j].score);
347 348 349
        }

      // Homopolymers
350 351
      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);
352
      
353
      // Local alignment
354 355
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	try_operation(best, BEGIN, 0, 0, 0);
Mikaël Salson's avatar
Mikaël Salson committed
356

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

360 361
      // Fill the score and the backtrack information
      B[i][j] = best ;
Mikaël Salson's avatar
Mikaël Salson committed
362
      
363 364
      // 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
365 366
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
367
	  int tbest = best.score ;
Mikaël Salson's avatar
Mikaël Salson committed
368 369 370 371 372 373 374 375 376 377 378

	  if (mode == LocalEndWithSomeDeletions)
	    tbest += cost.deletion_end*(n-j);   
     
	  if (tbest > best_score)
	    {
	      best_score = tbest ;
	      best_i = i ;
	      best_j = j ;
	    }
	    
379
	  if (best.score == 0){
380
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
381 382
	  }
	}
Mikaël Salson's avatar
Mikaël Salson committed
383

Mikaël Salson's avatar
Mikaël Salson committed
384 385
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
386
      	  if (best.score == 0){
387
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
388 389 390 391 392
	  }
	}
      
    }

393

394 395
  // End. Find best_i and best_j, put FIN keywords where the backtrack should stop
  
Mikaël Salson's avatar
Mikaël Salson committed
396 397
  if (mode == Local || mode == LocalEndWithSomeDeletions)
    {
398 399
      for (int j=0; j<=n; j++){
	B[0][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
400 401
      }
      for (int i=0; i<=m; i++){
402
	B[i][0].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
403 404 405 406 407 408 409 410
      }
      
    }
  if (mode == SemiGlobal)
    {
      best_i = m ;
      
      for (int j=1; j<=n; j++)
411
	if (B[m][j].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
412
	  {
413
	    best_score = B[m][j].score ;
Mikaël Salson's avatar
Mikaël Salson committed
414 415 416
	    best_j = j ;
	  }
      for (int i=0; i<=m; i++){
417
	B[i][0].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
418 419 420 421 422 423 424 425
      }
    }

  if (mode == SemiGlobalTrans)
    {
      best_j = n ;
      
      for (int i=1; i<=m; i++)
426
	if (B[i][n].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
427
	  {
428
	    best_score = B[i][n].score ;
Mikaël Salson's avatar
Mikaël Salson committed
429 430 431 432
	    best_i = i ;
	  }
	  
      for (int i=0; i<=n; i++){
433
	B[0][i].type = FIN;
434
      }     
Mikaël Salson's avatar
Mikaël Salson committed
435 436
    }

Mathieu Giraud's avatar
Mathieu Giraud committed
437
  if ((mode == Global) || (mode == GlobalButMostlyLocal))
Mikaël Salson's avatar
Mikaël Salson committed
438 439 440
    {
      best_i = m ;
      best_j = n ;
441
      best_score = B[m][n].score;
Mikaël Salson's avatar
Mikaël Salson committed
442 443 444 445 446 447 448 449
    }

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

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

450 451
  B[0][0].type = FIN;
  
Mikaël Salson's avatar
Mikaël Salson committed
452 453 454 455 456 457 458
  // In the matrix positions start at 1, but start positions may start at 0
  best_i -= 1;   
  best_j -= 1;

  return best_score ;
}

459 460 461 462 463 464 465 466 467 468 469 470 471
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
472 473
void DynProg::backtrack()
{
474
  // Tables for displaying the alignment
Mikaël Salson's avatar
Mikaël Salson committed
475 476 477 478
  gap1 = new int[x.size()+1];
  linkgap = new int[x.size()+1];
  gap2 = new int[y.size()+1];
      
479
  for (unsigned int i = 0; i <=x.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
480 481 482
    gap1[i] = 0;
    linkgap[i] = 0;
  }
483
  for (unsigned int i = 0; i <= y.size(); i++) {
Mikaël Salson's avatar
Mikaël Salson committed
484 485 486 487 488 489 490
    gap2[i] = 0;
  }
  
  
  int g1=x.size();
  int g2=y.size();
  
491 492
  int i=best_i+1;
  int j=best_j+1;
493

494 495 496 497 498 499 500
  // Retake good i/j when there were reversed strings
  if (reverse_x)
    i = m - i + 1 ;

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

501 502 503 504 505
  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
506
  
507
  // Compute backtrack strings
508 509
  ostringstream back_x;
  ostringstream back_y;
Mikaël Salson's avatar
Mikaël Salson committed
510 511
  ostringstream back_tr;
  
512 513
  while (1) {

514 515 516

    if ((!reverse_y && (j == marked_pos_j))
        || (reverse_y && (n-j+1 == marked_pos_j)))
517
      {
518
        marked_pos_i = reverse_x ? m-i+1 : i ;
519 520
      }

521 522
    if  (B[i][j].type == FIN)
      break ;
Mikaël Salson's avatar
Mikaël Salson committed
523
      
524
    // cout << "bt " << i << "/" << j << " " << B[i][j].type << " " << B[i][j].i << "," << B[i][j].j << endl ;
525 526 527
    
    int next_i = B[i][j].i;
    int next_j = B[i][j].j;
Mikaël Salson's avatar
Mikaël Salson committed
528

529 530 531 532 533 534 535 536 537
    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);
      }

538 539
    // 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
540
      
541 542 543
    for (int k=0; k<max_k; k++)
      {
        linkgap[g1]=g2;
544 545 546

        // Some character(s) from x, then fill with spaces
        if (i-k > next_i)
547
          {
548
            back_x << x[i-1-k];
549 550 551 552
            g1--;
          }
        else
          {
553
            back_x << " " ;
554 555
            gap1[g1]++ ;
          }
556 557 558

        // Some character(s) from y, then fill with spaces
        if (j-k > next_j)
559
          {
560
            back_y << y[j-1-k];
561 562 563 564
            g2--;
          }
        else
          {
565
            back_y << " " ;
566 567
            gap2[g2]++ ;
          }
568 569 570 571 572 573

        // The operation character, then fill with spaces
        if (k == max_k-1)          
          back_tr << B[i][j].type ;
        else
          back_tr << " " ;
574
      }
Mikaël Salson's avatar
Mikaël Salson committed
575
    
576 577
    i = next_i;
    j = next_j;
Mikaël Salson's avatar
Mikaël Salson committed
578 579
  }

580 581 582
  // Format backtrack strings
  
  string str1, str2, str3;
583
  str1 = back_x.str();
Mikaël Salson's avatar
Mikaël Salson committed
584 585 586 587
  str1 =string (str1.rbegin(), str1.rend());
  str1 = str1;
  str2=back_tr.str();
  str2 = string (str2.rbegin(), str2.rend());
588
  str3 = back_y.str();
Mikaël Salson's avatar
Mikaël Salson committed
589
  str3 = string (str3.rbegin(), str3.rend());
590 591

  ostringstream back;    
592 593 594
  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
595 596 597 598 599
  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;
  }
600
  back << "score: " << best_score << endl;
Mikaël Salson's avatar
Mikaël Salson committed
601
  
602 603
  first_i=i;
  first_j=j;
Mikaël Salson's avatar
Mikaël Salson committed
604 605
  
  str_back=back.str();
Mathieu Giraud's avatar
Mathieu Giraud committed
606 607

  // cout << str_back << endl ;
Mikaël Salson's avatar
Mikaël Salson committed
608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
}


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


625
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
626
{
627
  out << label << "  " ;
Mikaël Salson's avatar
Mikaël Salson committed
628

629 630
  for (int j=0; j<n; j++)
    out << setw(4) << y[j] << " ";
Mikaël Salson's avatar
Mikaël Salson committed
631 632 633

  out << endl ;

634
  for (int i=0; i<=m; i++)
Mikaël Salson's avatar
Mikaël Salson committed
635 636
    {
      if (i)
637
	out << x[i-1] << " " ;
Mikaël Salson's avatar
Mikaël Salson committed
638 639 640
      else
	out << "  " ;

641
      for (int j=0; j<=n; j++)
642
        {
643
          if (B[i][j].score)
644
            {
645
              if (B[i][j].score <= MINUS_INF)
646 647
                out << "-INF" ;
              else
648
                out << setw(4) << B[i][j].score ;
649
            }
650
          else
651 652
            out << "    ";

653
          out << B[i][j].type ;
654
        }
Mikaël Salson's avatar
Mikaël Salson committed
655 656
      out << endl ;      
    }
657 658

  out << endl;
Mikaël Salson's avatar
Mikaël Salson committed
659 660 661 662
  
  return out ;
}

663 664
ostream& operator<<(ostream& out, const DynProg& dp)
{
665 666 667 668 669 670
  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);
671 672 673 674
  out << "best: " << dp.best_i << "," << dp.best_j << endl;

  return out;
}
Mikaël Salson's avatar
Mikaël Salson committed
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 705 706 707 708 709 710 711 712 713 714

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