binaries/src/MSAProbs-0.9.7/MSAProbs/MSAPartProbs.cpp

   1 /***********************************************
   2  * # Copyright 2009-2010. Liu Yongchao
   3  * # Contact: Liu Yongchao, School of Computer Engineering,
   4  * #                     Nanyang Technological University.
   5  * # Emails:     liuy0039@ntu.edu.sg; nkcslyc@hotmail.com
   6  * #
   7  * # GPL version 3.0 applies.
   8  * #
   9  * ************************************************/
  10 #include "SafeVector.h"
  11 #include <iostream>
  12 #include <string.h>
  13 #include <stdlib.h>
  14 #include <stdio.h>
  15 #include <math.h>
  16 #include <time.h>
  17 #include <ctype.h>
  18 #include <assert.h>
  19 #define  TRACE 0                // 0: NOTRACE 1: TRACE
  20 //proba like settings
  21 #define  endgaps 1              // 1: engap penaties enabled 0: disabled
  22 #define  PART_FULL_MEMORY 0     //0: LOW MEM OPTION
  23 #define  REVPART_FULL_MEMORY 0  //0: LOW MEM OPTION
  24 using namespace std;
  25
  26 #ifdef _WIN32
  27 #define OS_HUGE_VALL    HUGE_VAL
  28 #else
  29 #define OS_HUGE_VALL    HUGE_VALL
  30 #endif
  31
  32 typedef struct {
  33         char input[30];
  34         int matrix;
  35         int N;
  36         float T;
  37         float beta;
  38         char opt;                       //can be 'P' or 'M'
  39         float gapopen;
  40         float gapext;
  41 } argument_decl;
  42
  43 typedef struct sequence {
  44         char *title;
  45         char *text;
  46         int length;
  47 } fasta;
  48
  49 typedef struct alignment {
  50         char *title;
  51         char *text;
  52         int length;
  53 } align;
  54
  55 ////////////////////////////////////////////////////////
  56 //externs related to scoring matrix and input arguments
  57 ///////////////////////////////////////////////////////////
  58 extern float g_gap_open1, g_gap_open2, g_gap_ext1, g_gap_ext2;
  59 extern char aminos[26], matrixtype[20], bases[26];
  60
  61 extern double sub_matrix[26][26];
  62 extern int subst_index[26];
  63
  64 extern float TEMPERATURE;
  65 extern int MATRIXTYPE;
  66
  67 extern float GAPOPEN;
  68 extern float GAPEXT;
  69 extern argument_decl argument;
  70
  71 //////////////////////////////////////////////////////////////////////////////
  72 //calculates reverse partition function values based on z matrices
  73 //and also simulaneously calculates the propability of each basepair
  74 //or aminoacid residue pair i,j
  75 //////////////////////////////////////////////////////////////////////////////
  76
  77 VF *revers_partf(fasta sequences[2], const double termgapopen,
  78                 const double termgapextend, long double **Zfm, const double d,
  79                 const double e) {
  80         // printf("revpart\n");
  81         //rest of the declarations
  82         int i, j;
  83         long double **Zm = NULL;
  84         long double **Ze = NULL;
  85         long double **Zf = NULL;
  86         int len0, len1;
  87         float probability;
  88         long double tempvar;
  89         int Si, Tj;
  90         double endgapopen, endgapextend;
  91         FILE *fo;
  92
  93         //Init lengths of sequences
  94         len0 = strlen(sequences[0].text);
  95         len1 = strlen(sequences[1].text);
  96
  97         //Safe vector declared
  98         VF *posteriorPtr = new VF((len0 + 1) * (len1 + 1));
  99         VF & posterior = *posteriorPtr;
 100         VF::iterator ptr = posterior.begin();
 101
 102         if (TRACE)                      //open the trace file
 103                 fo = fopen("revpartdump", "a");
 104
 105         //default:
 106         endgapopen = termgapopen;
 107         endgapextend = termgapextend;
 108
 109         //instantiate the z matrix
 110         if (REVPART_FULL_MEMORY) {
 111
 112                 Ze = new long double *[sequences[1].length + 1];
 113                 Zf = new long double *[sequences[1].length + 1];
 114                 Zm = new long double *[sequences[1].length + 1];
 115
 116                 if (TRACE)
 117                         printf("\n\n %e %e\n", d, e);
 118
 119                 //DYNAMICALLY GROW 2D Zm Zf Ze MARICES (long double)
 120                 for (i = 0; i <= sequences[1].length; i++) {
 121                         Ze[i] = new long double[sequences[0].length + 1];
 122                         Zf[i] = new long double[sequences[0].length + 1];
 123                         Zm[i] = new long double[sequences[0].length + 1];
 124                 }
 125         } else {
 126                 Zm = new long double *[2];
 127                 Ze = new long double *[2];
 128                 Zf = new long double *[2];
 129                 for (i = 0; i <= 1; i++) {
 130                         Zm[i] = new long double[sequences[0].length + 1];
 131                         Ze[i] = new long double[sequences[0].length + 1];
 132                         Zf[i] = new long double[sequences[0].length + 1];
 133                 }
 134
 135         }
 136
 137         if (TRACE) {
 138                 printf("in rev partf---");
 139                 printf("\n\n");
 140         }
 141
 142         if (REVPART_FULL_MEMORY) {
 143                 for (i = 0; i <= len1; i++)
 144                         for (j = 0; j <= len0; j++) {
 145                                 Zm[i][j] = 0.0;
 146                                 Zf[i][j] = 0.0;
 147                                 Ze[i][j] = 0.0;
 148                         }
 149         } else {
 150
 151                 for (j = 0; j <= len0; j++) {
 152                         Zm[0][j] = 0;
 153                         Zf[0][j] = 0;
 154                         Ze[0][j] = 0;
 155                         Zf[1][j] = 0;
 156                         Ze[1][j] = 0;
 157                         Zm[1][j] = 0;
 158                 }
 159         }
 160
 161         //fill the probability matrix with 0s
 162         for (i = 0; i <= len1; i++)
 163                 for (j = 0; j <= len0; j++)
 164                         ptr[j * (len1 + 1) + i] = 0;
 165
 166         if (endgaps == 0) {
 167                 Zm[len1][len0] = 1;
 168                 Ze[len1][len0] = Zf[len1][len0] = 0;
 169                 Zf[len1 - 1][len0] = Zm[len1][len0] * d;
 170                 Ze[len1][len0 - 1] = Zm[len1][len0] * d;
 171
 172                 //>=2ND ROW INIT
 173                 if (REVPART_FULL_MEMORY) {
 174                         for (i = len1 - 2; i >= 0; i--) {
 175                                 Zf[i][len0] = Zf[i + 1][len0] * e;
 176                         }
 177                 }
 178
 179                 //>=2ND COL INIT
 180                 if (REVPART_FULL_MEMORY) {
 181                         for (j = len0 - 2; j >= 0; j--) {
 182                                 Ze[len1][j] = Ze[len1][j + 1] * e;
 183                         }
 184                 } else {
 185                         for (j = len0 - 2; j >= 0; j--) {
 186                                 Ze[0][j] = Ze[0][j + 1] * e;
 187                         }
 188                 }
 189         } else {
 190
 191                 if (REVPART_FULL_MEMORY) {
 192
 193                         Zm[len1][len0] = 1;
 194                         Ze[len1][len0] = Zf[len1][len0] = 0;
 195                         Zf[len1 - 1][len0] = Zm[len1][len0] * endgapopen;
 196                         Ze[len1][len0 - 1] = Zm[len1][len0] * endgapopen;
 197
 198                         //>=2ND ROW INIT
 199                         for (i = len1 - 2; i >= 0; i--) {
 200                                 Zf[i][len0] = Zf[i + 1][len0] * endgapextend;
 201                         }
 202
 203                         //M Iy= d+j*e
 204
 205                         //>=2ND COL INIT
 206                         for (j = len0 - 2; j >= 0; j--) {
 207                                 Ze[len1][j] = Ze[len1][j + 1] * endgapextend;
 208                         }
 209
 210                 } else {
 211                         //in Zm
 212                         //let:
 213                         //  Zm(0) be the current row being filled/computed
 214                         //  Zm(1) be the previous row
 215
 216                         Zm[1][len0] = 1;
 217                         Ze[0][len0] = Zf[0][len0] = 0;
 218                         Zf[1][len0] = Zm[1][len0] * endgapopen;
 219                         Ze[0][len0 - 1] = Zm[1][len0] * endgapopen;
 220
 221                         //>=2ND COL INIT
 222                         for (j = len0 - 2; j >= 0; j--) {
 223                                 Ze[0][j] = Ze[0][j + 1] * endgapextend;
 224                         }
 225
 226                 }                       //END ELSE
 227
 228         }                               //END FULL MEMORY and GAP enablement IF STATEMENT
 229
 230         double scorez, zz = 0;
 231
 232         for (i = len1 - 1; i >= 0; i--) {
 233
 234                 for (j = len0 - 1; j >= 0; j--) {
 235                         Si = subst_index[sequences[1].text[i] - 'A'];
 236                         Tj = subst_index[sequences[0].text[j] - 'A'];
 237                         scorez = sub_matrix[Si][Tj];
 238
 239                         //endgaps modification aug 10
 240                         double open0, extend0, open1, extend1;
 241
 242                         open0 = open1 = d;
 243                         extend0 = extend1 = e;
 244
 245                         if (endgaps == 1) {
 246
 247                                 //check to see if one of the 2 sequences or both reach the end
 248
 249                                 if (i == 0) {
 250                                         open0 = endgapopen;
 251                                         extend0 = endgapextend;
 252
 253                                 }
 254
 255                                 if (j == 0) {
 256                                         open1 = endgapopen;
 257                                         extend1 = endgapextend;
 258                                 }
 259
 260                         }
 261
 262                         if (REVPART_FULL_MEMORY) {
 263                                 //z computation
 264
 265                                 Ze[i][j] = Zm[i][j + 1] * open0 + Ze[i][j + 1] * extend0;
 266                                 Zf[i][j] = Zm[i + 1][j] * open1 + Zf[i + 1][j] * extend1;
 267                                 Zm[i][j] = (Zm[i + 1][j + 1] + Zf[i + 1][j + 1]
 268                                                 + Ze[i + 1][j + 1]) * scorez;
 269                                 zz = Zm[i][j] + Zf[i][j] + Ze[i][j];
 270
 271                         } else {
 272
 273                                 //2 ROW zE zF ALGORITHM GOES...:
 274                                 //Ze[1][j] =Zm[i][j + 1] * exp(beta * open0) + Ze[1][j + 1] *exp(beta * extend0);
 275                                 //Zf[1][j] = Zm[i + 1][j] * exp(beta * open1) + Zf[0][j] * exp(beta * extend1);
 276                                 //Zm[i][j] = (Zm[i + 1][j + 1] + Zf[0][j + 1] + Ze[0][j + 1]) * exp(beta * scorez);
 277                                 //zz = Zm[0][j] + Zf[1][j] + Ze[1][j];
 278
 279                                 //lowmem code for merging probability calculating module
 280                                 //Here we make use of Zm as a 2 row matrix
 281
 282                                 Zf[1][j] = Zm[1][j] * open1 + Zf[0][j] * extend1;
 283                                 Ze[1][j] = Zm[0][j + 1] * open0 + Ze[1][j + 1] * extend0;
 284                                 Zm[0][j] = (Zm[1][j + 1] + Zf[0][j + 1] + Ze[0][j + 1])
 285                                                 * scorez;
 286
 287                                 tempvar = Zfm[i + 1][j + 1] * Zm[0][j];
 288                                 //divide P(i,j) i.e. pairwise probability by denominator
 289                                 tempvar /= (scorez * Zfm[0][0]);
 290                                 probability = (float) tempvar;
 291
 292                                 //store only noticable probabilities
 293                                 if (probability <= 1 && probability >= 0.001) {
 294                                         //algorithm goes...
 295                                         //validprob[i + 1][j + 1] = probability;
 296                                         ptr[(j + 1) * (len1 + 1) + (i + 1)] = probability;
 297                                 }
 298                                 //lowmem code ends here
 299
 300                         }
 301
 302                 }                       //end of for
 303
 304                 if (REVPART_FULL_MEMORY == 0) {
 305                         for (int t = 0; t <= sequences[0].length; t++) {
 306                                 Ze[0][t] = Ze[1][t];
 307                                 Ze[1][t] = 0;
 308
 309                                 Zf[0][t] = Zf[1][t];
 310                                 Zf[1][t] = 0;
 311
 312                                 Zm[1][t] = Zm[0][t];
 313                                 Zm[0][t] = 0;
 314
 315                         }
 316                         Zf[0][len0] = 1;
 317
 318                 }
 319
 320         }                               //end of for
 321
 322         if (TRACE) {
 323                 printf("\n\nrM:....\n\n");
 324                 if (REVPART_FULL_MEMORY) {
 325                         for (i = 0; i <= len1; i++) {
 326                                 for (j = 0; j <= len0; j++)
 327                                         printf("%.2Le ", Zm[i][j]);
 328                                 printf("\n");
 329                         }
 330
 331                         printf("\n\nrE:....\n\n");
 332                         for (i = 0; i <= len1; i++) {
 333                                 for (j = 0; j <= len0; j++)
 334                                         printf("%.2Le ", Ze[i][j]);
 335                                 printf("\n");
 336
 337                         }
 338
 339                         printf("\n\nrF:....\n\n");
 340                         for (i = 0; i <= len1; i++) {
 341                                 for (j = 0; j <= len0; j++)
 342                                         printf("%.2Le ", Zf[i][j]);
 343                                 printf("\n");
 344
 345                         }
 346
 347                 }
 348
 349         }
 350
 351         if (TRACE) {
 352                 fprintf(fo, "\n");
 353                 fclose(fo);
 354         }
 355
 356         //delete unused memory
 357
 358         if (REVPART_FULL_MEMORY) {
 359                 for (i = 0; i <= len1; i++) {
 360                         delete (Zm[i]);
 361                         delete (Zf[i]);
 362                         delete (Ze[i]);
 363                 }
 364         } else {
 365                 delete (Zf[0]);
 366                 delete (Ze[0]);
 367                 delete (Zm[0]);
 368
 369                 delete (Zm[1]);
 370                 delete (Zf[1]);
 371                 delete (Ze[1]);
 372         }
 373
 374         for (i = 0; i <= len1; i++) {
 375                 delete (Zfm[i]);
 376         }
 377
 378         if (Zf != NULL)
 379                 delete (Zf);
 380
 381         if (Ze != NULL)
 382                 delete (Ze);
 383
 384         if (Zm != NULL)
 385                 delete (Zm);
 386
 387         if (Zfm != NULL)
 388                 delete (Zfm);
 389
 390         posterior[0] = 0;
 391         return (posteriorPtr);
 392
 393 }
 394
 395 //////////////////////////////////////////////////////////////
 396 //forward partition function
 397 /////////////////////////////////////////////////////////////
 398
 399 long double **partf(fasta sequences[2], const double termgapopen,
 400                 const double termgapextend, const double d, const double e) {
 401         //printf("partf\n");
 402         int i, j, len1, len0;
 403         long double **Zm = NULL, **Zf = NULL, **Ze = NULL, zz = 0;
 404         double endgapopen, endgapextend;
 405
 406         //default:
 407         endgapopen = termgapopen;
 408         endgapextend = termgapextend;
 409
 410         //the flag endgaps is set at the #define section
 411         if (PART_FULL_MEMORY) {
 412
 413                 Zf = new long double *[sequences[1].length + 1];
 414                 Ze = new long double *[sequences[1].length + 1];
 415                 Zm = new long double *[sequences[1].length + 1];
 416
 417                 //comment
 418                 if (TRACE)
 419                         printf("\nPARTF:====\n");
 420
 421                 //DYNAMICALLY GROW 2D M,IX,IY,PIX,PIY MARICES
 422                 for (i = 0; i <= sequences[1].length; i++) {
 423                         Zf[i] = new long double[sequences[0].length + 1];
 424                         Ze[i] = new long double[sequences[0].length + 1];
 425                         Zm[i] = new long double[sequences[0].length + 1];
 426                 }
 427         } else {
 428                 Zm = new long double *[sequences[1].length + 1];
 429                 Ze = new long double *[2];
 430                 Zf = new long double *[2];
 431                 for (i = 0; i <= sequences[1].length; i++) {
 432                         Zm[i] = new long double[sequences[0].length + 1];
 433                 }
 434                 Ze[0] = new long double[sequences[0].length + 1];
 435                 Zf[0] = new long double[sequences[0].length + 1];
 436                 Ze[1] = new long double[sequences[0].length + 1];
 437                 Zf[1] = new long double[sequences[0].length + 1];
 438         }
 439
 440         len0 = strlen(sequences[0].text);
 441         len1 = strlen(sequences[1].text);
 442
 443         if (PART_FULL_MEMORY) {
 444                 for (i = 0; i <= sequences[1].length; i++)
 445                         for (j = 0; j <= sequences[0].length; j++) {
 446                                 Zm[i][j] = 0.00;
 447                                 Zf[i][j] = 0.00;
 448                                 Ze[i][j] = 0.00;
 449                         }
 450         } else {
 451                 for (i = 0; i <= len1; i++) {
 452                         for (j = 0; j <= len0; j++) {
 453                                 Zm[i][j] = 0;
 454                         }
 455                 }
 456                 for (j = 0; j <= len0; j++) {
 457                         Zf[0][j] = 0;
 458                         Ze[0][j] = 0;
 459                         Zf[1][j] = 0;
 460                         Ze[1][j] = 0;
 461                 }
 462         }
 463
 464         //INTITIALIZE THE DP
 465
 466         if (endgaps == 0) {
 467                 Zm[0][0] = 1.00;
 468
 469                 Zf[0][0] = Ze[0][0] = 0;
 470                 Zf[1][0] = Zm[0][0] * d;
 471                 Ze[0][1] = Zm[0][0] * d;
 472
 473                 //>=2ND ROW INIT
 474                 if (PART_FULL_MEMORY) {
 475                         for (i = 2; i <= sequences[1].length; i++) {
 476                                 Zf[i][0] = Zf[i - 1][0] * e;
 477                         }
 478                 }
 479
 480                 //>=2ND COL INIT
 481                 for (j = 2; j <= sequences[0].length; j++) {
 482                         Ze[0][j] = Ze[0][j - 1] * e;
 483                 }
 484         } else {
 485                 //init z
 486                 Zm[0][0] = 1.00;
 487                 Zf[0][0] = Ze[0][0] = 0;
 488                 Zf[1][0] = Zm[0][0] * endgapopen;
 489                 Ze[0][1] = Zm[0][0] * endgapopen;
 490
 491                 //>=2ND ROW INIT
 492                 if (PART_FULL_MEMORY) {
 493                         for (i = 2; i <= sequences[1].length; i++) {
 494                                 Zf[i][0] = Zf[i - 1][0] * endgapextend;
 495                         }
 496                 }
 497
 498                 //>=2ND COL INIT
 499                 for (j = 2; j <= sequences[0].length; j++) {
 500                         Ze[0][j] = Ze[0][j - 1] * endgapextend;
 501                 }
 502         }
 503
 504         //1ST ROW/COL INIT
 505
 506         int Si, Tj;
 507         double score;
 508
 509         for (i = 1; i <= sequences[1].length; i++) {
 510
 511                 for (j = 1; j <= sequences[0].length; j++) {
 512
 513                         Si = subst_index[sequences[1].text[i - 1] - 'A'];
 514                         Tj = subst_index[sequences[0].text[j - 1] - 'A'];
 515
 516                         score = sub_matrix[Si][Tj];
 517
 518                         double open0, extend0, open1, extend1;
 519
 520                         open0 = open1 = d;
 521                         extend0 = extend1 = e;
 522
 523                         if (endgaps == 1) {
 524                                 //check to see if one of the 2 sequences or both reach the end
 525
 526                                 if (i == sequences[1].length) {
 527                                         open0 = endgapopen;
 528                                         extend0 = endgapextend;
 529
 530                                 }
 531
 532                                 if (j == sequences[0].length) {
 533                                         open1 = endgapopen;
 534                                         extend1 = endgapextend;
 535                                 }
 536                         }
 537
 538                         //
 539                         //z computation using open and extend temp vars
 540                         //open0 is gap open in seq0 and open1 is gap open in seq1
 541                         //entend0 is gap extend in seq0 and extend1 is gap extend in seq1
 542
 543                         if (PART_FULL_MEMORY) {
 544                                 Ze[i][j] = Zm[i][j - 1] * open0 + Ze[i][j - 1] * extend0;
 545
 546                                 if (Ze[i][j] >= OS_HUGE_VALL) {
 547                                         printf("ERROR: huge val error for Ze\n");
 548                                         exit(1);
 549                                 }
 550
 551                                 Zf[i][j] = Zm[i - 1][j] * open1 + Zf[i - 1][j] * extend1;
 552
 553                                 if (Zf[i][j] >= OS_HUGE_VALL) {
 554                                         printf("ERROR: huge val error for Zf\n");
 555                                         exit(1);
 556                                 }
 557
 558                                 Zm[i][j] = (Zm[i - 1][j - 1] + Ze[i - 1][j - 1]
 559                                                 + Zf[i - 1][j - 1]) * score;
 560
 561                                 if (Zm[i][j] >= OS_HUGE_VALL) {
 562                                         printf("ERROR: huge val error for Zm\n");
 563                                         exit(1);
 564                                 }
 565
 566                                 zz = Zm[i][j] + Ze[i][j] + Zf[i][j];
 567                         } else {
 568                                 Ze[1][j] = Zm[i][j - 1] * open0 + Ze[1][j - 1] * extend0;
 569
 570                                 if (Ze[1][j] >= OS_HUGE_VALL) {
 571                                         printf("ERROR: huge val error for zE\n");
 572                                         exit(1);
 573                                 }
 574
 575                                 Zf[1][j] = Zm[i - 1][j] * open1 + Zf[0][j] * extend1;
 576
 577                                 if (Zf[1][j] >= OS_HUGE_VALL) {
 578                                         printf("ERROR: huge val error for zF\n");
 579                                         exit(1);
 580                                 }
 581
 582                                 Zm[i][j] = (Zm[i - 1][j - 1] + Ze[0][j - 1] + Zf[0][j - 1])
 583                                                 * score;
 584
 585                                 if (Zm[i][j] >= OS_HUGE_VALL) {
 586                                         printf("ERROR: huge val error for zM\n");
 587                                         exit(1);
 588                                 }
 589
 590                                 zz = Zm[i][j] + Ze[1][j] + Zf[1][j];
 591                         }
 592
 593                 }                       //end for
 594
 595                 if (!PART_FULL_MEMORY) {
 596                         for (int t = 0; t <= sequences[0].length; t++) {
 597                                 Ze[0][t] = Ze[1][t];
 598                                 Ze[1][t] = 0;
 599
 600                                 Zf[0][t] = Zf[1][t];
 601                                 Zf[1][t] = 0;
 602                         }
 603
 604                         Zf[1][0] = 1;
 605
 606                 }
 607
 608         }                               //end for
 609
 610         //store the sum of zm zf ze (m,n)s in zm's 0,0 th position
 611         Zm[0][0] = zz;
 612
 613         if (TRACE) {
 614                 //debug code aug 3
 615                 //print the 3 Z matrices namely Zm Zf and Ze
 616
 617                 printf("\n\nFINAL Zm:\n");
 618                 for (i = 0; i <= sequences[1].length; i++) {
 619                         for (j = 0; j <= sequences[0].length; j++)
 620                                 printf("%.2Le ", Zm[i][j]);
 621                         printf("\n");
 622                 }
 623
 624                 printf("FINAL Zf \n");
 625                 for (i = 0; i <= sequences[1].length; i++) {
 626                         for (j = 0; j <= sequences[0].length; j++)
 627                                 printf("%.2Le ", Zf[i][j]);
 628                         printf("\n");
 629                 }
 630
 631                 printf("FINAL Ze \n");
 632                 for (i = 0; i <= sequences[1].length; i++) {
 633                         for (j = 0; j <= sequences[0].length; j++)
 634                                 printf("%.2Le ", Ze[i][j]);
 635                         printf("\n");
 636                 }
 637
 638                 //end debug dump code
 639
 640         }
 641
 642         if (PART_FULL_MEMORY) {
 643                 for (i = 0; i <= sequences[1].length; i++) {
 644                         delete (Zf[i]);
 645                         delete (Ze[i]);
 646                 }
 647         } else {
 648                 delete (Zf[0]);
 649                 delete (Ze[0]);
 650                 delete (Zf[1]);
 651                 delete (Ze[1]);
 652         }
 653
 654         delete (Zf);
 655         delete (Ze);
 656
 657         return Zm;
 658
 659 }                               //end of forward partition function
 660
 661 /////////////////////////////////////////////////////////////////////////////////////////
 662 //entry point (was the main function) , returns the posterior probability safe vector
 663 ////////////////////////////////////////////////////////////////////////////////////////
 664 VF *ComputePostProbs(int a, int b, string seq1, string seq2) {
 665         //printf("probamod\n");
 666         double gap_open = -22, gap_ext = -1, beta = 0.2;//T = 5, beta = 1/T = 0.2, by default
 667         int stock_loop = 1;
 668         int le = 160;
 669         double termgapopen = 1.0f;      //exp(0)
 670         double termgapextend = 1.0f;    //exp(0)
 671
 672         //initialize the sequence structure
 673         fasta sequences[2];
 674
 675         sequences[0].length = strlen((char *) seq1.c_str());
 676         sequences[0].text = (char *) seq1.c_str();
 677         sequences[0].title = new char[10];
 678         strcpy(sequences[0].title, "seq0");
 679         sequences[1].length = strlen((char *) seq2.c_str());
 680         sequences[1].text = (char *) seq2.c_str();
 681         sequences[1].title = new char[10];
 682         strcpy(sequences[1].title, "seq1");
 683
 684         if (TRACE)
 685
 686         {
 687                 printf("%d %d %s\n%d %d %s\n--\n", a, sequences[0].length,
 688                                 sequences[0].text, b, sequences[1].length, sequences[1].text);
 689                 printf("after init\n");
 690
 691                 FILE *dump1 = fopen("dump1", "a");
 692                 fprintf(dump1, "%d %d %s\n%d %d %s\n--\n", a, sequences[0].length,
 693                                 sequences[0].text, b, sequences[1].length, sequences[1].text);
 694                 fclose(dump1);
 695         }
 696
 697         gap_open = argument.gapopen;
 698         gap_ext = argument.gapext;
 699         beta = argument.beta;
 700
 701         stock_loop = argument.N;
 702         le = argument.matrix;
 703
 704         //compute the values of exp(beta * ?)
 705         termgapopen = exp(beta * 0.0);
 706         termgapextend = exp(beta * 0.0);
 707         gap_open = exp(beta * gap_open);
 708         gap_ext = exp(beta * gap_ext);
 709
 710         if (TRACE)
 711                 printf("%f %f %f %d\n", gap_open, gap_ext, beta, le);
 712
 713         //call for calculating the posterior probabilities
 714         // 1. call partition function partf
 715         // 2. calculate revpartition using revers_parf
 716         // 3. calculate probabilities
 717         /// MODIFICATION... POPULATE SAFE VECTOR
 718
 719         long double **MAT1;
 720
 721         MAT1 = partf(sequences, termgapopen, termgapextend, gap_open, gap_ext);
 722
 723         return revers_partf(sequences, termgapopen, termgapextend, MAT1, gap_open,
 724                         gap_ext);
 725
 726 }
 727
 728 //end of posterior probability  module