dynprog.cpp 15.7 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 45 46 47 48 49 50 51 52 53

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);
54

55 56 57
  this -> open_insertion = this -> open_deletion = MINUS_INF ;
  this -> extend_insertion = this -> extend_deletion = MINUS_INF ;
  this -> affine_gap = false ;
58 59

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

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

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


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

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

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

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

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

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

132 133 134
  this -> marked_pos_j = marked_pos_j;
  this -> marked_pos_i = 0;

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

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

  init();
}

DynProg::~DynProg() {
167

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

void DynProg::init()
182
{
183 184
  // Initialize backtrack labels (B[0][0] labels are never used)

185 186 187 188 189 190 191 192 193 194 195 196 197 198
  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 ;
    }

199 200
  // Initialize scores

201 202
  for (int i=0; i<=m; i++)
    for (int j=0; j<=n; j++)
203
      {
204
        B[i][j].score = 0 ;
205
        if (cost.affine_gap) {
206 207
          Bins[i][j].score = 0;
          Bdel[i][j].score = 0;
208 209
        }
      }
210 211

  if (!(mode == Local || mode == LocalEndWithSomeDeletions)) // Global, SemiGlobal
212
    for (int i=0; i<=m; i++) {
213
      B[i][0].score = i * cost.insertion ;
214 215 216 217 218 219 220 221 222
      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;
    }

223
  if (!(mode == SemiGlobal || mode == SemiGlobalTrans || mode == Local || mode == LocalEndWithSomeDeletions)) // Global
224
    for (int j=0; j<=n; j++) {
225
      B[0][j].score = j * cost.deletion ;
226 227 228 229 230 231 232 233
       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
234 235
}

236

237 238 239 240 241 242 243 244 245 246 247
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 ;
    }
}

248
int DynProg::compute(bool onlyBottomTriangle, int onlyBottomTriangleShift)
Mikaël Salson's avatar
Mikaël Salson committed
249 250
{
  best_score = MINUS_INF ;
251 252 253 254
  best_i = 0 ;
  best_j = 0 ;
  
  operation best ;
Mikaël Salson's avatar
Mikaël Salson committed
255 256 257 258

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

261 262 263 264 265 266 267 268 269 270 271 272 273
      // 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 ;
        }

274
      // The edit operations, with their backtracking information and their score
275 276

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

279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
      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
301 302
      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);
303
      
304
      // Local alignment
305 306
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	try_operation(best, BEGIN, 0, 0, 0);
Mikaël Salson's avatar
Mikaël Salson committed
307

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

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

	  if (mode == LocalEndWithSomeDeletions)
	    tbest += cost.deletion_end*(n-j);   
     
	  if (tbest > best_score)
	    {
	      best_score = tbest ;
	      best_i = i ;
	      best_j = j ;
	    }
	    
330
	  if (best.score == 0){
331
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
332 333
	  }
	}
Mikaël Salson's avatar
Mikaël Salson committed
334

Mikaël Salson's avatar
Mikaël Salson committed
335 336
      if (mode == Local || mode == LocalEndWithSomeDeletions)
	{
337
      	  if (best.score == 0){
338
	    B[i][j].type = FIN;
Mikaël Salson's avatar
Mikaël Salson committed
339 340 341 342 343
	  }
	}
      
    }

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

  if (mode == SemiGlobalTrans)
    {
      best_j = n ;
      
      for (int i=1; i<=m; i++)
376
	if (B[i][n].score > best_score)
Mikaël Salson's avatar
Mikaël Salson committed
377
	  {
378
	    best_score = B[i][n].score ;
Mikaël Salson's avatar
Mikaël Salson committed
379 380 381 382
	    best_i = i ;
	  }
	  
      for (int i=0; i<=n; i++){
383
	B[0][i].type = FIN;
384
      }     
Mikaël Salson's avatar
Mikaël Salson committed
385 386
    }

Mathieu Giraud's avatar
Mathieu Giraud committed
387
  if ((mode == Global) || (mode == GlobalButMostlyLocal))
Mikaël Salson's avatar
Mikaël Salson committed
388 389 390
    {
      best_i = m ;
      best_j = n ;
391
      best_score = B[m][n].score;
Mikaël Salson's avatar
Mikaël Salson committed
392 393 394 395 396 397 398 399
    }

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

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

400 401
  B[0][0].type = FIN;
  
Mikaël Salson's avatar
Mikaël Salson committed
402 403 404 405 406 407 408
  // In the matrix positions start at 1, but start positions may start at 0
  best_i -= 1;   
  best_j -= 1;

  return best_score ;
}

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

444 445 446 447 448 449 450
  // Retake good i/j when there were reversed strings
  if (reverse_x)
    i = m - i + 1 ;

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

451 452 453 454 455
  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
456
  
457
  // Compute backtrack strings
458 459
  ostringstream back_x;
  ostringstream back_y;
Mikaël Salson's avatar
Mikaël Salson committed
460 461
  ostringstream back_tr;
  
462 463
  while (1) {

464 465 466

    if ((!reverse_y && (j == marked_pos_j))
        || (reverse_y && (n-j+1 == marked_pos_j)))
467
      {
468
        marked_pos_i = reverse_x ? m-i+1 : i ;
469 470
      }

471 472
    if  (B[i][j].type == FIN)
      break ;
Mikaël Salson's avatar
Mikaël Salson committed
473
      
474
    // cout << "bt " << i << "/" << j << " " << B[i][j].type << " " << B[i][j].i << "," << B[i][j].j << endl ;
475 476 477
    
    int next_i = B[i][j].i;
    int next_j = B[i][j].j;
Mikaël Salson's avatar
Mikaël Salson committed
478

479 480 481 482 483 484 485 486 487
    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);
      }

488 489
    // 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
490
      
491 492 493
    for (int k=0; k<max_k; k++)
      {
        linkgap[g1]=g2;
494 495 496

        // Some character(s) from x, then fill with spaces
        if (i-k > next_i)
497
          {
498
            back_x << x[i-1-k];
499 500 501 502
            g1--;
          }
        else
          {
503
            back_x << " " ;
504 505
            gap1[g1]++ ;
          }
506 507 508

        // Some character(s) from y, then fill with spaces
        if (j-k > next_j)
509
          {
510
            back_y << y[j-1-k];
511 512 513 514
            g2--;
          }
        else
          {
515
            back_y << " " ;
516 517
            gap2[g2]++ ;
          }
518 519 520 521 522 523

        // The operation character, then fill with spaces
        if (k == max_k-1)          
          back_tr << B[i][j].type ;
        else
          back_tr << " " ;
524
      }
Mikaël Salson's avatar
Mikaël Salson committed
525
    
526 527
    i = next_i;
    j = next_j;
Mikaël Salson's avatar
Mikaël Salson committed
528 529
  }

530 531 532
  // Format backtrack strings
  
  string str1, str2, str3;
533
  str1 = back_x.str();
Mikaël Salson's avatar
Mikaël Salson committed
534 535 536 537
  str1 =string (str1.rbegin(), str1.rend());
  str1 = str1;
  str2=back_tr.str();
  str2 = string (str2.rbegin(), str2.rend());
538
  str3 = back_y.str();
Mikaël Salson's avatar
Mikaël Salson committed
539
  str3 = string (str3.rbegin(), str3.rend());
540 541

  ostringstream back;    
542 543 544
  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
545 546 547 548 549
  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;
  }
550
  back << "score: " << best_score << endl;
Mikaël Salson's avatar
Mikaël Salson committed
551
  
552 553
  first_i=i;
  first_j=j;
Mikaël Salson's avatar
Mikaël Salson committed
554 555
  
  str_back=back.str();
Mathieu Giraud's avatar
Mathieu Giraud committed
556 557

  // cout << str_back << endl ;
Mikaël Salson's avatar
Mikaël Salson committed
558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574
}


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


575
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
576
{
577
  out << label << "  " ;
Mikaël Salson's avatar
Mikaël Salson committed
578

579 580
  for (int j=0; j<n; j++)
    out << setw(4) << y[j] << " ";
Mikaël Salson's avatar
Mikaël Salson committed
581 582 583

  out << endl ;

584
  for (int i=0; i<=m; i++)
Mikaël Salson's avatar
Mikaël Salson committed
585 586
    {
      if (i)
587
	out << x[i-1] << " " ;
Mikaël Salson's avatar
Mikaël Salson committed
588 589 590
      else
	out << "  " ;

591
      for (int j=0; j<=n; j++)
592
        {
593
          if (B[i][j].score)
594
            {
595
              if (B[i][j].score <= MINUS_INF)
596 597
                out << "-INF" ;
              else
598
                out << setw(4) << B[i][j].score ;
599
            }
600
          else
601 602
            out << "    ";

603
          out << B[i][j].type ;
604
        }
Mikaël Salson's avatar
Mikaël Salson committed
605 606
      out << endl ;      
    }
607 608

  out << endl;
Mikaël Salson's avatar
Mikaël Salson committed
609 610 611 612
  
  return out ;
}

613 614
ostream& operator<<(ostream& out, const DynProg& dp)
{
615 616 617 618 619 620
  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);
621 622 623 624
  out << "best: " << dp.best_i << "," << dp.best_j << endl;

  return out;
}
Mikaël Salson's avatar
Mikaël Salson committed
625 626 627 628 629 630 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

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