binaries/src/mafft/core/pairlocalalign.c

   1 #include "mltaln.h"
   2
   3 #define DEBUG 0
   4 #define IODEBUG 0
   5 #define SCOREOUT 0
   6 #define SHISHAGONYU 0 // for debug
   7
   8 #define NODIST -9999
   9
  10 static char *whereispairalign;
  11 static char *laraparams;
  12 static char foldalignopt[1000];
  13 static int stdout_align;
  14 static int stdout_dist;
  15 static int store_localhom;
  16 static int store_dist;
  17 static int nadd;
  18 static int laste;
  19 static int lastm;
  20 static int lastsubopt;
  21 static int lastonce;
  22 static int usenaivescoreinsteadofalignmentscore;
  23 static int specifictarget;
  24
  25 typedef struct _lastres
  26 {
  27         int score;
  28         int start1;
  29         int start2;
  30         char *aln1;
  31         char *aln2;
  32 } Lastres;
  33
  34 typedef struct _reg
  35 {
  36         int start;
  37         int end;
  38 } Reg;
  39
  40 typedef struct _aln
  41 {
  42         int nreg;
  43         Reg *reg1;
  44         Reg *reg2;
  45 } Aln;
  46
  47 typedef struct _lastresx
  48 {
  49         int score;
  50         int naln;
  51         Aln *aln;
  52 } Lastresx;
  53
  54 #ifdef enablemultithread
  55 typedef struct _jobtable
  56 {
  57         int i;
  58         int j;
  59 } Jobtable;
  60
  61 typedef struct _thread_arg
  62 {
  63         int thread_no;
  64         int njob;
  65         Jobtable *jobpospt;
  66         char **name;
  67         char **seq;
  68         char **dseq;
  69         int *thereisxineachseq;
  70         LocalHom **localhomtable;
  71         double **distancemtx;
  72         double *selfscore;
  73         char ***bpp;
  74         Lastresx **lastresx;
  75         int alloclen;
  76         int *targetmap;
  77         pthread_mutex_t *mutex_counter;
  78         pthread_mutex_t *mutex_stdout;
  79 } thread_arg_t;
  80 #endif
  81
  82 typedef struct _lastcallthread_arg
  83 {
  84         int nq, nd;
  85         char **dseq;
  86         char **qseq;
  87         Lastresx **lastresx;
  88 #ifdef enablemultithread
  89         int thread_no;
  90         int *kshare;
  91         pthread_mutex_t *mutex;
  92 #endif
  93 } lastcallthread_arg_t;
  94
  95 static void t2u( char *seq )
  96 {
  97         while( *seq )
  98         {
  99                 if     ( *seq == 'A' ) *seq = 'a';
 100                 else if( *seq == 'a' ) *seq = 'a';
 101                 else if( *seq == 'T' ) *seq = 'u';
 102                 else if( *seq == 't' ) *seq = 'u';
 103                 else if( *seq == 'U' ) *seq = 'u';
 104                 else if( *seq == 'u' ) *seq = 'u';
 105                 else if( *seq == 'G' ) *seq = 'g';
 106                 else if( *seq == 'g' ) *seq = 'g';
 107                 else if( *seq == 'C' ) *seq = 'c';
 108                 else if( *seq == 'c' ) *seq = 'c';
 109                 else *seq = 'n';
 110                 seq++;
 111         }
 112 }
 113
 114 static int removex( char *d, char *m )
 115 {
 116         int val = 0;
 117         while( *m != 0 )
 118         {
 119                 if( *m == 'X' || *m == 'x' )
 120                 {
 121                         m++;
 122                         val++;
 123                 }
 124                 else
 125                 {
 126                         *d++ = *m++;
 127                 }
 128         }
 129         *d = 0;
 130         return( val );
 131 }
 132
 133 static void putlocalhom_last( char *s1, char *s2, LocalHom *localhompt, Lastresx *lastresx, char korh )
 134 {
 135         char *pt1, *pt2;
 136         int naln, nreg;
 137         int iscore;
 138         int isumscore;
 139         int sumoverlap;
 140         LocalHom *tmppt = localhompt;
 141         LocalHom *tmppt2;
 142         LocalHom *localhompt0;
 143         Reg *rpt1, *rpt2;
 144         Aln *apt;
 145         int nlocalhom = 0;
 146         int len;
 147
 148 //      fprintf( stderr, "s1=%s\n", s1 );
 149 //      fprintf( stderr, "s2=%s\n", s2 );
 150
 151
 152         naln = lastresx->naln;
 153         apt = lastresx->aln;
 154
 155         if( naln == 0 ) return;
 156         while( naln-- )
 157         {
 158                 rpt1 = apt->reg1;
 159                 rpt2 = apt->reg2;
 160                 nreg = apt->nreg;
 161                 isumscore = 0;
 162                 sumoverlap = 0;
 163                 while( nreg-- )
 164                 {
 165                         if( nlocalhom++ > 0 )
 166                         {
 167 //                              fprintf( stderr, "reallocating ...\n" );
 168                                 tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
 169 //                              fprintf( stderr, "done\n" );
 170                                 tmppt = tmppt->next;
 171                                 tmppt->next = NULL;
 172                         }
 173                         tmppt->start1 = rpt1->start;
 174                         tmppt->start2 = rpt2->start;
 175                         tmppt->end1   = rpt1->end;
 176                         tmppt->end2   = rpt2->end;
 177                         tmppt->korh   = 'h';
 178                         if( rpt1 == apt->reg1 ) localhompt0 = tmppt; // ?
 179
 180 //                      fprintf( stderr, "in putlocalhom, reg1: %d-%d (nreg=%d)\n", rpt1->start, rpt1->end, lastresx->nreg );
 181 //                      fprintf( stderr, "in putlocalhom, reg2: %d-%d (nreg=%d)\n", rpt2->start, rpt2->end, lastresx->nreg );
 182
 183                         len = tmppt->end1 - tmppt->start1 + 1;
 184
 185 //                      fprintf( stderr, "tmppt->start1=%d\n", tmppt->start1 );
 186 //                      fprintf( stderr, "tmppt->start2=%d\n", tmppt->start2 );
 187
 188 //                      fprintf( stderr, "s1+tmppt->start1=%*.*s\n", len, len, s1+tmppt->start1 );
 189 //                      fprintf( stderr, "s2+tmppt->start2=%*.*s\n", len, len, s2+tmppt->start2 );
 190
 191                         pt1 = s1 + tmppt->start1;
 192                         pt2 = s2 + tmppt->start2;
 193                         iscore = 0;
 194                         while( len-- )
 195                         {
 196                                 iscore += n_dis[(int)amino_n[(unsigned char)*pt1++]][(int)amino_n[(unsigned char)*pt2++]]; // - offset \e$B$O$$$i$J$$$+$b\e(B
 197 //                              fprintf( stderr, "len=%d, %c-%c, iscore(0) = %d\n", len, *(pt1-1), *(pt2-1), iscore );
 198                         }
 199
 200                         if( divpairscore )
 201                         {
 202                                 tmppt->overlapaa   = tmppt->end2-tmppt->start2+1;
 203                                 tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
 204                         }
 205                         else
 206                         {
 207                                 isumscore += iscore;
 208                                 sumoverlap += tmppt->end2-tmppt->start2+1;
 209                         }
 210                         rpt1++;
 211                         rpt2++;
 212                 }
 213 #if 0
 214                 fprintf( stderr, "iscore (1)= %d\n", iscore );
 215                 fprintf( stderr, "al1: %d - %d\n", start1, end1 );
 216                 fprintf( stderr, "al2: %d - %d\n", start2, end2 );
 217 #endif
 218
 219                 if( !divpairscore )
 220                 {
 221                         for( tmppt2=localhompt0; tmppt2; tmppt2=tmppt2->next )
 222                         {
 223                                 tmppt2->overlapaa = sumoverlap;
 224                                 tmppt2->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
 225 //                              fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
 226                         }
 227                 }
 228                 apt++;
 229         }
 230 }
 231
 232 static int countcomma( char *s )
 233 {
 234         int v = 0;
 235         while( *s ) if( *s++ == ',' ) v++;
 236         return( v );
 237 }
 238
 239 static double recallpairfoldalign( char **mseq1, char **mseq2, int m1, int m2, int *of1pt, int *of2pt, int alloclen )
 240 {
 241         static FILE *fp = NULL;
 242         double value;
 243         char *aln1;
 244         char *aln2;
 245         int of1tmp, of2tmp;
 246
 247         if( fp == NULL )
 248         {
 249                 fp = fopen( "_foldalignout", "r" );
 250                 if( fp == NULL )
 251                 {
 252                         fprintf( stderr, "Cannot open _foldalignout\n" );
 253                         exit( 1 );
 254                 }
 255         }
 256
 257         aln1 = calloc( alloclen, sizeof( char ) );
 258         aln2 = calloc( alloclen, sizeof( char ) );
 259
 260         readpairfoldalign( fp, *mseq1, *mseq2, aln1, aln2, m1, m2, &of1tmp, &of2tmp, alloclen );
 261
 262         if( strstr( foldalignopt, "-global") )
 263         {
 264                 fprintf( stderr, "Calling G__align11\n" );
 265                 value = G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
 266                 *of1pt = 0;
 267                 *of2pt = 0;
 268         }
 269         else
 270         {
 271                 fprintf( stderr, "Calling L__align11\n" );
 272                 value = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, of1pt, of2pt );
 273         }
 274
 275 //      value = (double)naivepairscore11( *mseq1, *mseq2, penalty ); // nennnotame
 276
 277         if( aln1[0] == 0 )
 278         {
 279                 fprintf( stderr, "FOLDALIGN returned no alignment between %d and %d.  Sequence alignment is used instead.\n", m1+1, m2+1 );
 280         }
 281         else
 282         {
 283                 strcpy( *mseq1, aln1 );
 284                 strcpy( *mseq2, aln2 );
 285                 *of1pt = of1tmp;
 286                 *of2pt = of2tmp;
 287         }
 288
 289 //      value = naivepairscore11( *mseq1, *mseq2, penalty ); // v6.511 ha kore wo tsukau, global nomi dakara.
 290
 291 //      fclose( fp ); // saigo dake yatta houga yoi.
 292
 293 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
 294 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
 295
 296
 297         free( aln1 );
 298         free( aln2 );
 299
 300         return( value );
 301 }
 302
 303 static void block2reg( char *block, Reg *reg1, Reg *reg2, int start1, int start2 )
 304 {
 305         Reg *rpt1, *rpt2;
 306         char *tpt, *npt;
 307         int pos1, pos2;
 308         int len, glen1, glen2;
 309         pos1 = start1;
 310         pos2 = start2;
 311         rpt1 = reg1;
 312         rpt2 = reg2;
 313         while( block )
 314         {
 315                 block++;
 316 //              fprintf( stderr, "block = %s\n", block );
 317                 tpt = strchr( block, ':' );
 318                 npt = strchr( block, ',' );
 319                 if( !tpt || tpt > npt )
 320                 {
 321                         len = atoi( block );
 322                         reg1->start = pos1;
 323                         reg2->start = pos2;
 324                         pos1 += len - 1;
 325                         pos2 += len - 1;
 326                         reg1->end = pos1;
 327                         reg2->end = pos2;
 328 //                      fprintf( stderr, "in loop reg1: %d-%d\n", reg1->start, reg1->end );
 329 //                      fprintf( stderr, "in loop reg2: %d-%d\n", reg2->start, reg2->end );
 330                         reg1++;
 331                         reg2++;
 332                 }
 333                 else
 334                 {
 335                         sscanf( block, "%d:%d", &glen1, &glen2 );
 336                         pos1 += glen1 + 1;
 337                         pos2 += glen2 + 1;
 338                 }
 339                 block = npt;
 340
 341         }
 342         reg1->start = reg1->end = reg2->start = reg2->end = -1;
 343
 344         while( rpt1->start != -1 )
 345         {
 346 //              fprintf( stderr, "reg1: %d-%d\n", rpt1->start, rpt1->end );
 347 //              fprintf( stderr, "reg2: %d-%d\n", rpt2->start, rpt2->end );
 348                 rpt1++;
 349                 rpt2++;
 350         }
 351 //      *apt1 = *apt2 = 0;
 352 //      fprintf( stderr, "aln1 = %s\n", aln1 );
 353 //      fprintf( stderr, "aln2 = %s\n", aln2 );
 354 }
 355
 356
 357 static void readlastresx_singleq( FILE *fp, int n1, int nameq, Lastresx **lastresx )
 358 {
 359         char *gett;
 360         Aln *tmpaln;
 361         int prevnaln, naln, nreg;
 362 #if 0
 363         int i, pstart, pend, end1, end2;
 364 #endif
 365         int score, name1, start1, alnSize1, seqSize1;
 366         int        name2, start2, alnSize2, seqSize2;
 367         char strand1, strand2;
 368         int includeintoscore;
 369         gett = calloc( 10000, sizeof( char ) );
 370
 371 //      fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
 372 //      fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
 373
 374         while( 1 )
 375         {
 376                 fgets( gett, 9999, fp );
 377                 if( feof( fp ) ) break;
 378                 if( gett[0] == '#' ) continue;
 379 //              fprintf( stdout, "gett = %s\n", gett );
 380                 if( gett[strlen(gett)-1] != '\n' )
 381                 {
 382                         fprintf( stderr, "Too long line?\n" );
 383                         exit( 1 );
 384                 }
 385
 386                 sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
 387                                         &score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
 388                                                 &name2, &start2, &alnSize2, &strand2, &seqSize2 );
 389
 390                 if( alg == 'R' && name2 <= name1 ) continue;
 391                 if( name2 != nameq )
 392                 {
 393                         fprintf( stderr, "BUG!!!\n" );
 394                         exit( 1 );
 395                 }
 396
 397 //              if( lastresx[name1][name2].score ) continue; // dame!!!!
 398
 399
 400                 prevnaln = lastresx[name1][name2].naln;
 401 #if 0
 402                 for( i=0; i<prevnaln; i++ )
 403                 {
 404                         nreg = lastresx[name1][name2].aln[i].nreg;
 405
 406                         pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
 407                         pend   = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
 408                         end1 = start1 + alnSize1;
 409 //                      fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
 410                         if( pstart <= start1 && start1 <= pend && pend - start1 > 1 ) break;
 411                         if( pstart <= end1   && end1   <= pend && end1 - pstart > 1 ) break;
 412
 413                         pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
 414                         pend   = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
 415                         end2 = start2 + alnSize2;
 416 //                      fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
 417                         if( pstart <= start2 && start2 <= pend && pend - start2 > 1 ) break;
 418                         if( pstart <= end2   && end2   <= pend && end2 - pstart > 1 ) break;
 419                 }
 420                 includeintoscore = ( i == prevnaln );
 421 #else
 422                 if( prevnaln ) includeintoscore = 0;
 423                 else includeintoscore = 1;
 424 #endif
 425                 if( !includeintoscore && !lastsubopt )
 426                         continue;
 427
 428                 naln = prevnaln + 1;
 429                 lastresx[name1][name2].naln = naln;
 430 //              fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
 431
 432                 if( ( tmpaln = (Aln *)realloc( lastresx[name1][name2].aln, (naln) * sizeof( Aln ) ) ) == NULL ) // yoyu nashi
 433                 {
 434                         fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
 435                         exit( 1 );
 436                 }
 437                 else
 438                         lastresx[name1][name2].aln = tmpaln;
 439
 440                 nreg = countcomma( gett )/2 + 1;
 441                 lastresx[name1][name2].aln[prevnaln].nreg = nreg;
 442 //              lastresx[name1][name2].aln[naln].nreg = -1;
 443 //              lastresx[name1][name2].aln[naln].reg1 = NULL;
 444 //              lastresx[name1][name2].aln[naln].reg2 = NULL;
 445 //              fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
 446
 447                 if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
 448                 {
 449                         fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
 450                         exit( 1 );
 451                 }
 452
 453                 if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
 454                 {
 455                         fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
 456                         exit( 1 );
 457                 }
 458
 459 //              lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
 460 //              lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
 461                 block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
 462
 463                 if( includeintoscore )
 464                 {
 465                         if( lastresx[name1][name2].score ) score += penalty;
 466                         lastresx[name1][name2].score += score;
 467                 }
 468
 469 //              fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
 470         }
 471         free( gett );
 472 }
 473
 474 #ifdef enablemultithread
 475 #if 0
 476 static void readlastresx_group( FILE *fp, Lastresx **lastresx )
 477 {
 478         char *gett;
 479         Aln *tmpaln;
 480         int prevnaln, naln, nreg;
 481 #if 0
 482         int i, pstart, pend, end1, end2;
 483 #endif
 484         int score, name1, start1, alnSize1, seqSize1;
 485         int        name2, start2, alnSize2, seqSize2;
 486         char strand1, strand2;
 487         int includeintoscore;
 488         gett = calloc( 10000, sizeof( char ) );
 489
 490 //      fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
 491 //      fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
 492
 493         while( 1 )
 494         {
 495                 fgets( gett, 9999, fp );
 496                 if( feof( fp ) ) break;
 497                 if( gett[0] == '#' ) continue;
 498 //              fprintf( stdout, "gett = %s\n", gett );
 499                 if( gett[strlen(gett)-1] != '\n' )
 500                 {
 501                         fprintf( stderr, "Too long line?\n" );
 502                         exit( 1 );
 503                 }
 504
 505                 sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
 506                                         &score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
 507                                                 &name2, &start2, &alnSize2, &strand2, &seqSize2 );
 508
 509                 if( alg == 'R' && name2 <= name1 ) continue;
 510
 511 //              if( lastresx[name1][name2].score ) continue; // dame!!!!
 512
 513                 prevnaln = lastresx[name1][name2].naln;
 514 #if 0
 515                 for( i=0; i<prevnaln; i++ )
 516                 {
 517                         nreg = lastresx[name1][name2].aln[i].nreg;
 518
 519                         pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
 520                         pend   = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
 521                         end1 = start1 + alnSize1;
 522 //                      fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
 523                         if( pstart <= start1 && start1 <= pend && pend - start1 > 3 ) break;
 524                         if( pstart <= end1   && end1   <= pend && end1 - pstart > 3 ) break;
 525
 526                         pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
 527                         pend   = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
 528                         end2 = start2 + alnSize2;
 529 //                      fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
 530                         if( pstart <= start2 && start2 <= pend && pend - start2 > 3 ) break;
 531                         if( pstart <= end2   && end2   <= pend && end2 - pstart > 3 ) break;
 532                 }
 533                 includeintoscore = ( i == prevnaln );
 534 #else
 535                 if( prevnaln ) includeintoscore = 0;
 536                 else includeintoscore = 1;
 537 #endif
 538                 if( !includeintoscore && !lastsubopt )
 539                         continue;
 540
 541                 naln = prevnaln + 1;
 542                 lastresx[name1][name2].naln = naln;
 543 //              fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
 544
 545                 if( ( tmpaln = (Aln *)realloc( lastresx[name1][name2].aln, (naln) * sizeof( Aln ) ) ) == NULL ) // yoyu nashi
 546                 {
 547                         fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
 548                         exit( 1 );
 549                 }
 550                 else
 551                         lastresx[name1][name2].aln = tmpaln;
 552
 553
 554
 555                 nreg = countcomma( gett )/2 + 1;
 556                 lastresx[name1][name2].aln[prevnaln].nreg = nreg;
 557 //              lastresx[name1][name2].aln[naln].nreg = -1;
 558 //              lastresx[name1][name2].aln[naln].reg1 = NULL;
 559 //              lastresx[name1][name2].aln[naln].reg2 = NULL;
 560 //              fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
 561
 562                 if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
 563                 {
 564                         fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
 565                         exit( 1 );
 566                 }
 567
 568                 if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
 569                 {
 570                         fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
 571                         exit( 1 );
 572                 }
 573
 574 //              lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
 575 //              lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
 576                 block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
 577
 578                 if( includeintoscore )
 579                 {
 580                         if( lastresx[name1][name2].score ) score += penalty;
 581                         lastresx[name1][name2].score += score;
 582                 }
 583
 584 //              fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
 585         }
 586         free( gett );
 587 }
 588 #endif
 589 #endif
 590
 591 static void readlastresx( FILE *fp, int n1, int n2, Lastresx **lastresx, char **seq1, char **seq2 )
 592 {
 593         char *gett;
 594         Aln *tmpaln;
 595         int prevnaln, naln, nreg;
 596 #if 0
 597         int i, pstart, pend, end1, end2;
 598 #endif
 599         int score, name1, start1, alnSize1, seqSize1;
 600         int        name2, start2, alnSize2, seqSize2;
 601         char strand1, strand2;
 602         int includeintoscore;
 603         gett = calloc( 10000, sizeof( char ) );
 604
 605 //      fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
 606 //      fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
 607
 608         while( 1 )
 609         {
 610                 fgets( gett, 9999, fp );
 611                 if( feof( fp ) ) break;
 612                 if( gett[0] == '#' ) continue;
 613 //              fprintf( stdout, "gett = %s\n", gett );
 614                 if( gett[strlen(gett)-1] != '\n' )
 615                 {
 616                         fprintf( stderr, "Too long line?\n" );
 617                         exit( 1 );
 618                 }
 619
 620                 sscanf( gett, "%d %d %d %d %c %d %d %d %d %c %d",
 621                                         &score, &name1, &start1, &alnSize1, &strand1, &seqSize1,
 622                                                 &name2, &start2, &alnSize2, &strand2, &seqSize2 );
 623
 624                 if( alg == 'R' && name2 <= name1 ) continue;
 625
 626 //              if( lastresx[name1][name2].score ) continue; // dame!!!!
 627
 628                 prevnaln = lastresx[name1][name2].naln;
 629 #if 0
 630                 for( i=0; i<prevnaln; i++ )
 631                 {
 632                         nreg = lastresx[name1][name2].aln[i].nreg;
 633
 634                         pstart = lastresx[name1][name2].aln[i].reg1[0].start + 0;
 635                         pend   = lastresx[name1][name2].aln[i].reg1[nreg-1].end - 0;
 636                         end1 = start1 + alnSize1;
 637 //                      fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
 638                         if( pstart <= start1 && start1 <= pend && pend - start1 > 3 ) break;
 639                         if( pstart <= end1   && end1   <= pend && end1 - pstart > 3 ) break;
 640
 641                         pstart = lastresx[name1][name2].aln[i].reg2[0].start + 0;
 642                         pend   = lastresx[name1][name2].aln[i].reg2[nreg-1].end - 0;
 643                         end2 = start2 + alnSize2;
 644 //                      fprintf( stderr, "pstart = %d, pend = %d\n", pstart, pend );
 645                         if( pstart <= start2 && start2 <= pend && pend - start2 > 3 ) break;
 646                         if( pstart <= end2   && end2   <= pend && end2 - pstart > 3 ) break;
 647                 }
 648                 includeintoscore = ( i == prevnaln );
 649 #else
 650                 if( prevnaln ) includeintoscore = 0;
 651                 else includeintoscore = 1;
 652 #endif
 653                 if( !includeintoscore && !lastsubopt )
 654                         continue;
 655
 656                 naln = prevnaln + 1;
 657                 lastresx[name1][name2].naln = naln;
 658 //              fprintf( stderr, "OK! add this alignment to hat3, %d-%d, naln = %d->%d\n", name1, name2, prevnaln, naln );
 659
 660                 if( ( tmpaln = (Aln *)realloc( lastresx[name1][name2].aln, (naln) * sizeof( Aln ) ) ) == NULL ) // yoyu nashi
 661                 {
 662                         fprintf( stderr, "Cannot reallocate lastresx[][].aln\n" );
 663                         exit( 1 );
 664                 }
 665                 else
 666                         lastresx[name1][name2].aln = tmpaln;
 667
 668
 669
 670                 nreg = countcomma( gett )/2 + 1;
 671                 lastresx[name1][name2].aln[prevnaln].nreg = nreg;
 672 //              lastresx[name1][name2].aln[naln].nreg = -1;
 673 //              lastresx[name1][name2].aln[naln].reg1 = NULL;
 674 //              lastresx[name1][name2].aln[naln].reg2 = NULL;
 675 //              fprintf( stderr, "name1=%d, name2=%d, nreg=%d, prevnaln=%d\n", name1, name2, nreg, prevnaln );
 676
 677                 if( ( lastresx[name1][name2].aln[prevnaln].reg1 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
 678                 {
 679                         fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
 680                         exit( 1 );
 681                 }
 682
 683                 if( ( lastresx[name1][name2].aln[prevnaln].reg2 = (Reg *)calloc( nreg+1, sizeof( Reg ) ) ) == NULL ) // yoyu nashi
 684                 {
 685                         fprintf( stderr, "Cannot reallocate lastresx[][].reg2\n" );
 686                         exit( 1 );
 687                 }
 688
 689 //              lastresx[name1][name2].aln[prevnaln].reg1[0].start = -1; // iranai?
 690 //              lastresx[name1][name2].aln[prevnaln].reg2[0].start = -1; // iranai?
 691                 block2reg( strrchr( gett, '\t' ), lastresx[name1][name2].aln[prevnaln].reg1, lastresx[name1][name2].aln[prevnaln].reg2, start1, start2 );
 692
 693                 if( includeintoscore )
 694                 {
 695                         if( lastresx[name1][name2].score ) score += penalty;
 696                         lastresx[name1][name2].score += score;
 697                 }
 698
 699 //              fprintf( stderr, "score(%d,%d) = %d\n", name1, name2, lastresx[name1][name2].score );
 700         }
 701         free( gett );
 702 }
 703
 704 #ifdef enablemultithread
 705 #if 0
 706 static void *lastcallthread_group( void *arg )
 707 {
 708         lastcallthread_arg_t *targ = (lastcallthread_arg_t *)arg;
 709         int k, i;
 710         int nq = targ->nq;
 711         int nd = targ->nd;
 712 #ifdef enablemultithread
 713         int thread_no = targ->thread_no;
 714         int *kshare = targ->kshare;
 715 #endif
 716         Lastresx **lastresx = targ->lastresx;
 717         char **dseq = targ->dseq;
 718         char **qseq = targ->qseq;
 719         char command[5000];
 720         FILE *lfp;
 721         int msize;
 722         int klim;
 723         int qstart, qend, shou, amari;
 724         char kd[1000];
 725
 726         if( nthread )
 727         {
 728                 shou = nq / nthread;
 729                 amari = nq - shou * nthread;
 730                 fprintf( stderr, "shou: %d, amari: %d\n", shou, amari );
 731
 732                 qstart = thread_no * shou;
 733                 if( thread_no - 1 < amari ) qstart += thread_no;
 734                 else qstart += amari;
 735
 736                 qend = qstart + shou - 1;
 737                 if( thread_no < amari ) qend += 1;
 738                 fprintf( stderr, "%d: %d-%d\n", thread_no, qstart, qend );
 739         }
 740         k = -1;
 741         while( 1 )
 742         {
 743                 if( nthread )
 744                 {
 745                         if( qstart > qend ) break;
 746                         if( k == thread_no ) break;
 747                         fprintf( stderr, "\n%d-%d / %d (thread %d)                    \n", qstart, qend, nq, thread_no );
 748                         k = thread_no;
 749                 }
 750                 else
 751                 {
 752                         k++;
 753                         if( k == nq ) break;
 754                         fprintf( stderr, "\r%d / %d                    \r", k, nq );
 755                 }
 756
 757                 if( alg == 'R' ) // if 'r' -> calllast_fast
 758                 {
 759                         fprintf( stderr, "Not supported\n" );
 760                         exit( 1 );
 761                 }
 762                 else // 'r'
 763                 {
 764                         kd[0] = 0;
 765                 }
 766
 767                 sprintf( command, "_q%d", k );
 768                 lfp = fopen( command, "w" );
 769                 if( !lfp )
 770                 {
 771                         fprintf( stderr, "Cannot open %s", command );
 772                         exit( 1 );
 773                 }
 774                 for( i=qstart; i<=qend; i++ )
 775                         fprintf( lfp, ">%d\n%s\n", i, qseq[i] );
 776                 fclose( lfp );
 777
 778 //              if( alg == 'R' ) msize = MAX(10,k+nq);
 779 //                      else msize = MAX(10,nd+nq);
 780                 if( alg == 'R' ) msize = MAX(10,k*lastm);
 781                         else msize = MAX(10,nd*lastm);
 782
 783 //              fprintf( stderr, "Calling lastal from lastcallthread, msize = %d, k=%d\n", msize, k );
 784 //              sprintf( command, "grep '>' _db%sd", kd );
 785 //              system( command );
 786                 sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q%d > _lastres%d", whereispairalign, msize, laste, -penalty, -penalty_ex, kd, k, k );
 787                 if( system( command ) ) exit( 1 );
 788
 789                 sprintf( command, "_lastres%d", k );
 790                 lfp = fopen( command, "r" );
 791                 if( !lfp )
 792                 {
 793                         fprintf( stderr, "Cannot read _lastres%d", k );
 794                         exit( 1 );
 795                 }
 796 //              readlastres( lfp, nd, nq, lastres, dseq, qseq );
 797 //              fprintf( stderr, "Reading lastres\n" );
 798                 readlastresx_group( lfp, lastresx );
 799                 fclose( lfp );
 800         }
 801         return( NULL );
 802 }
 803 #endif
 804 #endif
 805
 806 static void *lastcallthread( void *arg )
 807 {
 808         lastcallthread_arg_t *targ = (lastcallthread_arg_t *)arg;
 809         int k, i;
 810         int nq = targ->nq;
 811         int nd = targ->nd;
 812 #ifdef enablemultithread
 813         int thread_no = targ->thread_no;
 814         int *kshare = targ->kshare;
 815 #endif
 816         Lastresx **lastresx = targ->lastresx;
 817         char **dseq = targ->dseq;
 818         char **qseq = targ->qseq;
 819         char command[5000];
 820         FILE *lfp;
 821         int msize;
 822         int klim;
 823         char kd[1000];
 824
 825         k = -1;
 826         while( 1 )
 827         {
 828
 829 #ifdef enablemultithread
 830                 if( nthread )
 831                 {
 832                         pthread_mutex_lock( targ->mutex );
 833                         k = *kshare;
 834                         if( k == nq )
 835                         {
 836                                 pthread_mutex_unlock( targ->mutex );
 837                                 break;
 838                         }
 839                         fprintf( stderr, "\r%d / %d (thread %d)                    \r", k, nq, thread_no );
 840                         ++(*kshare);
 841                         pthread_mutex_unlock( targ->mutex );
 842                 }
 843                 else
 844 #endif
 845                 {
 846                         k++;
 847                         if( k == nq ) break;
 848                         fprintf( stderr, "\r%d / %d                    \r", k, nq );
 849                 }
 850
 851                 if( alg == 'R' ) // if 'r' -> calllast_fast
 852                 {
 853                         klim = MIN( k, njob-nadd );
 854 //                      klim = k; // dochira demo yoi
 855                         if( klim == k )
 856                         {
 857                                 sprintf( command, "_db%dd", k );
 858                                 lfp = fopen( command, "w" );
 859                                 if( !lfp )
 860                                 {
 861                                         fprintf( stderr, "Cannot open _db." );
 862                                         exit( 1 );
 863                                 }
 864                                 for( i=0; i<klim; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
 865                                 fclose( lfp );
 866
 867 //                              sprintf( command, "md5sum _db%dd > /dev/tty", k );
 868 //                              system( command );
 869
 870                                 if( dorp == 'd' )
 871                                         sprintf( command, "%s/lastdb _db%dd _db%dd", whereispairalign, k, k );
 872                                 else
 873                                         sprintf( command, "%s/lastdb -p _db%dd _db%dd", whereispairalign, k, k );
 874                                 system( command );
 875                                 sprintf( kd, "%d", k );
 876                         }
 877                         else // calllast_fast de tsukutta nowo riyou
 878                         {
 879                                 kd[0] = 0;
 880 //                              fprintf( stderr, "klim=%d, njob=%d, nadd=%d, skip!\n", klim, njob, nadd );
 881                         }
 882                 }
 883                 else // 'r'
 884                 {
 885                         kd[0] = 0;
 886                 }
 887
 888                 sprintf( command, "_q%d", k );
 889                 lfp = fopen( command, "w" );
 890                 if( !lfp )
 891                 {
 892                         fprintf( stderr, "Cannot open %s", command );
 893                         exit( 1 );
 894                 }
 895                 fprintf( lfp, ">%d\n%s\n", k, qseq[k] );
 896                 fclose( lfp );
 897
 898 //              if( alg == 'R' ) msize = MAX(10,k+nq);
 899 //                      else msize = MAX(10,nd+nq);
 900                 if( alg == 'R' ) msize = MAX(10,k*lastm);
 901                         else msize = MAX(10,nd*lastm);
 902
 903 //              fprintf( stderr, "Calling lastal from lastcallthread, msize = %d, k=%d\n", msize, k );
 904 //              sprintf( command, "grep '>' _db%sd", kd );
 905 //              system( command );
 906                 sprintf( command, "%s/lastal -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db%sd _q%d > _lastres%d", whereispairalign, msize, laste, -penalty, -penalty_ex, kd, k, k );
 907                 if( system( command ) ) exit( 1 );
 908
 909                 sprintf( command, "_lastres%d", k );
 910                 lfp = fopen( command, "r" );
 911                 if( !lfp )
 912                 {
 913                         fprintf( stderr, "Cannot read _lastres%d", k );
 914                         exit( 1 );
 915                 }
 916 //              readlastres( lfp, nd, nq, lastres, dseq, qseq );
 917 //              fprintf( stderr, "Reading lastres\n" );
 918                 readlastresx_singleq( lfp, nd, k, lastresx );
 919                 fclose( lfp );
 920         }
 921         return( NULL );
 922 }
 923
 924
 925 static void calllast_fast( int nd, char **dseq, int nq, char **qseq, Lastresx **lastresx )
 926 {
 927         int i, j;
 928         FILE *lfp;
 929         char command[1000];
 930
 931         lfp = fopen( "_scoringmatrixforlast", "w" );
 932         if( !lfp )
 933         {
 934                 fprintf( stderr, "Cannot open _scoringmatrixforlast" );
 935                 exit( 1 );
 936         }
 937         if( dorp == 'd' )
 938         {
 939                 fprintf( lfp, "      " );
 940                 for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
 941                 fprintf( lfp, "\n" );
 942                 for( i=0; i<4; i++ )
 943                 {
 944                         fprintf( lfp, "%c ", amino[i] );
 945                         for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
 946                         fprintf( lfp, "\n" );
 947                 }
 948         }
 949         else
 950         {
 951                 fprintf( lfp, "      " );
 952                 for( j=0; j<20; j++ ) fprintf( lfp, " %c ", amino[j] );
 953                 fprintf( lfp, "\n" );
 954                 for( i=0; i<20; i++ )
 955                 {
 956                         fprintf( lfp, "%c ", amino[i] );
 957                         for( j=0; j<20; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
 958                         fprintf( lfp, "\n" );
 959                 }
 960         }
 961         fclose( lfp );
 962
 963 //      if( alg == 'r' ) // if 'R' -> lastcallthread, kokonoha nadd>0 no toki nomi shiyou
 964         {
 965                 sprintf( command, "_dbd" );
 966                 lfp = fopen( command, "w" );
 967                 if( !lfp )
 968                 {
 969                         fprintf( stderr, "Cannot open _dbd" );
 970                         exit( 1 );
 971                 }
 972                 if( alg == 'R' )
 973                         j = njob-nadd;
 974                 else
 975                         j = nd;
 976                 for( i=0; i<j; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
 977
 978                 fclose( lfp );
 979                 if( dorp == 'd' )
 980                         sprintf( command, "%s/lastdb _dbd _dbd", whereispairalign );
 981                 else
 982                         sprintf( command, "%s/lastdb -p _dbd _dbd", whereispairalign );
 983                 system( command );
 984         }
 985
 986 #ifdef enablemultithread
 987         if( nthread )
 988         {
 989                 pthread_t *handle;
 990                 pthread_mutex_t mutex;
 991                 lastcallthread_arg_t *targ;
 992                 int *ksharept;
 993                 targ = (lastcallthread_arg_t *)calloc( nthread, sizeof( lastcallthread_arg_t ) );
 994                 handle = calloc( nthread, sizeof( pthread_t ) );
 995                 ksharept = calloc( 1, sizeof(int) );
 996                 *ksharept = 0;
 997                 pthread_mutex_init( &mutex, NULL );
 998                 for( i=0; i<nthread; i++ )
 999                 {
1000                         targ[i].thread_no = i;
1001                         targ[i].kshare = ksharept;
1002                         targ[i].nq = nq;
1003                         targ[i].nd = nd;
1004                         targ[i].dseq = dseq;
1005                         targ[i].qseq = qseq;
1006                         targ[i].lastresx = lastresx;
1007                         targ[i].mutex = &mutex;
1008                         pthread_create( handle+i, NULL, lastcallthread, (void *)(targ+i) );
1009                 }
1010
1011                 for( i=0; i<nthread; i++ )
1012                 {
1013                         pthread_join( handle[i], NULL );
1014                 }
1015                 pthread_mutex_destroy( &mutex );
1016                 free( handle );
1017                 free( targ );
1018                 free( ksharept );
1019         }
1020         else
1021 #endif
1022         {
1023                 lastcallthread_arg_t *targ;
1024                 targ = (lastcallthread_arg_t *)calloc( 1, sizeof( lastcallthread_arg_t ) );
1025                 targ[0].nq = nq;
1026                 targ[0].nd = nd;
1027                 targ[0].dseq = dseq;
1028                 targ[0].qseq = qseq;
1029                 targ[0].lastresx = lastresx;
1030                 lastcallthread( targ );
1031                 free( targ );
1032         }
1033
1034 }
1035
1036 static void calllast_once( int nd, char **dseq, int nq, char **qseq, Lastresx **lastresx )
1037 {
1038         int i, j;
1039         char command[5000];
1040         FILE *lfp;
1041         int msize;
1042         int res;
1043
1044         fprintf( stderr, "nq=%d\n", nq );
1045
1046         lfp = fopen( "_db", "w" );
1047         if( !lfp )
1048         {
1049                 fprintf( stderr, "Cannot open _db" );
1050                 exit( 1 );
1051         }
1052         for( i=0; i<nd; i++ ) fprintf( lfp, ">%d\n%s\n", i, dseq[i] );
1053         fclose( lfp );
1054
1055         if( dorp == 'd' )
1056         {
1057                 sprintf( command, "%s/lastdb _db _db", whereispairalign );
1058                 system( command );
1059                 lfp = fopen( "_scoringmatrixforlast", "w" );
1060                 if( !lfp )
1061                 {
1062                         fprintf( stderr, "Cannot open _scoringmatrixforlast" );
1063                         exit( 1 );
1064                 }
1065                 fprintf( lfp, "      " );
1066                 for( j=0; j<4; j++ ) fprintf( lfp, " %c ", amino[j] );
1067                 fprintf( lfp, "\n" );
1068                 for( i=0; i<4; i++ )
1069                 {
1070                         fprintf( lfp, "%c ", amino[i] );
1071                         for( j=0; j<4; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
1072                         fprintf( lfp, "\n" );
1073                 }
1074                 fclose( lfp );
1075 #if 0
1076                 sprintf( command, "lastex -s 2 -a %d -b %d -p _scoringmatrixforlast -E 10000 _db.prj _db.prj > _lastex", -penalty, -penalty_ex );
1077                 system( command );
1078                 lfp = fopen( "_lastex", "r" );
1079                 fgets( command, 4999, lfp );
1080                 fgets( command, 4999, lfp );
1081                 fgets( command, 4999, lfp );
1082                 fgets( command, 4999, lfp );
1083                 laste = atoi( command );
1084                 fclose( lfp );
1085                 fprintf( stderr, "laste = %d\n", laste );
1086                 sleep( 10 );
1087 #else
1088 //              laste = 5000;
1089 #endif
1090         }
1091         else
1092         {
1093                 sprintf( command, "%s/lastdb -p _db _db", whereispairalign );
1094                 system( command );
1095                 lfp = fopen( "_scoringmatrixforlast", "w" );
1096                 if( !lfp )
1097                 {
1098                         fprintf( stderr, "Cannot open _scoringmatrixforlast" );
1099                         exit( 1 );
1100                 }
1101                 fprintf( lfp, "      " );
1102                 for( j=0; j<20; j++ ) fprintf( lfp, " %c ", amino[j] );
1103                 fprintf( lfp, "\n" );
1104                 for( i=0; i<20; i++ )
1105                 {
1106                         fprintf( lfp, "%c ", amino[i] );
1107                         for( j=0; j<20; j++ ) fprintf( lfp, " %d ", n_dis[i][j] );
1108                         fprintf( lfp, "\n" );
1109                 }
1110                 fclose( lfp );
1111 //              fprintf( stderr, "Not written yet\n" );
1112         }
1113
1114         lfp = fopen( "_q", "w" );
1115         if( !lfp )
1116         {
1117                 fprintf( stderr, "Cannot open _q" );
1118                 exit( 1 );
1119         }
1120         for( i=0; i<nq; i++ )
1121         {
1122                 fprintf( lfp, ">%d\n%s\n", i, qseq[i] );
1123         }
1124         fclose( lfp );
1125
1126         msize = MAX(10,nd*lastm);
1127
1128 //      fprintf( stderr, "Calling lastal from calllast_once, msize=%d\n", msize );
1129         sprintf( command, "%s/lastal -v -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", whereispairalign, msize, laste, -penalty, -penalty_ex );
1130 //      sprintf( command, "lastal -v -m %d -e %d -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", 1, laste, -penalty, -penalty_ex );
1131 //      sprintf( command, "lastal -v -e 40 -f 0 -s 1 -p _scoringmatrixforlast -a %d -b %d _db _q > _lastres", -penalty, -penalty_ex );
1132         res = system( command );
1133         if( res )
1134         {
1135                 fprintf( stderr, "LAST aborted\n" );
1136                 exit( 1 );
1137         }
1138
1139         lfp = fopen( "_lastres", "r" );
1140         if( !lfp )
1141         {
1142                 fprintf( stderr, "Cannot read _lastres" );
1143                 exit( 1 );
1144         }
1145 //      readlastres( lfp, nd, nq, lastres, dseq, qseq );
1146         fprintf( stderr, "Reading lastres\n" );
1147         readlastresx( lfp, nd, nq, lastresx, dseq, qseq );
1148         fclose( lfp );
1149 }
1150
1151 static void callfoldalign( int nseq, char **mseq )
1152 {
1153         FILE *fp;
1154         int i;
1155         int res;
1156         static char com[10000];
1157
1158         for( i=0; i<nseq; i++ )
1159                 t2u( mseq[i] );
1160
1161         fp = fopen( "_foldalignin", "w" );
1162         if( !fp )
1163         {
1164                 fprintf( stderr, "Cannot open _foldalignin\n" );
1165                 exit( 1 );
1166         }
1167         for( i=0; i<nseq; i++ )
1168         {
1169                 fprintf( fp, ">%d\n", i+1 );
1170                 fprintf( fp, "%s\n", mseq[i] );
1171         }
1172         fclose( fp );
1173
1174         sprintf( com, "env PATH=%s  foldalign210 %s _foldalignin > _foldalignout ", whereispairalign, foldalignopt );
1175         res = system( com );
1176         if( res )
1177         {
1178                 fprintf( stderr, "Error in foldalign\n" );
1179                 exit( 1 );
1180         }
1181
1182 }
1183
1184 static void calllara( int nseq, char **mseq, char *laraarg )
1185 {
1186         FILE *fp;
1187         int i;
1188         int res;
1189         static char com[10000];
1190
1191 //      for( i=0; i<nseq; i++ )
1192
1193         fp = fopen( "_larain", "w" );
1194         if( !fp )
1195         {
1196                 fprintf( stderr, "Cannot open _larain\n" );
1197                 exit( 1 );
1198         }
1199         for( i=0; i<nseq; i++ )
1200         {
1201                 fprintf( fp, ">%d\n", i+1 );
1202                 fprintf( fp, "%s\n", mseq[i] );
1203         }
1204         fclose( fp );
1205
1206
1207 //      fprintf( stderr, "calling LaRA\n" );
1208         sprintf( com, "env PATH=%s:/bin:/usr/bin mafft_lara -i _larain -w _laraout -o _lara.params %s", whereispairalign, laraarg );
1209         res = system( com );
1210         if( res )
1211         {
1212                 fprintf( stderr, "Error in lara\n" );
1213                 exit( 1 );
1214         }
1215 }
1216
1217 static double recalllara( char **mseq1, char **mseq2, int alloclen )
1218 {
1219         static FILE *fp = NULL;
1220         static char *ungap1;
1221         static char *ungap2;
1222         static char *ori1;
1223         static char *ori2;
1224 //      int res;
1225         static char com[10000];
1226         double value;
1227
1228
1229         if( fp == NULL )
1230         {
1231                 fp = fopen( "_laraout", "r" );
1232                 if( fp == NULL )
1233                 {
1234                         fprintf( stderr, "Cannot open _laraout\n" );
1235                         exit( 1 );
1236                 }
1237                 ungap1 = AllocateCharVec( alloclen );
1238                 ungap2 = AllocateCharVec( alloclen );
1239                 ori1 = AllocateCharVec( alloclen );
1240                 ori2 = AllocateCharVec( alloclen );
1241         }
1242
1243
1244         strcpy( ori1, *mseq1 );
1245         strcpy( ori2, *mseq2 );
1246
1247         fgets( com, 999, fp );
1248         myfgets( com, 9999, fp );
1249         strcpy( *mseq1, com );
1250         myfgets( com, 9999, fp );
1251         strcpy( *mseq2, com );
1252
1253         gappick0( ungap1, *mseq1 );
1254         gappick0( ungap2, *mseq2 );
1255         t2u( ungap1 );
1256         t2u( ungap2 );
1257         t2u( ori1 );
1258         t2u( ori2 );
1259
1260         if( strcmp( ungap1, ori1 ) || strcmp( ungap2, ori2 ) )
1261         {
1262                 fprintf( stderr, "SEQUENCE CHANGED!!\n" );
1263                 fprintf( stderr, "*mseq1  = %s\n", *mseq1 );
1264                 fprintf( stderr, "ungap1  = %s\n", ungap1 );
1265                 fprintf( stderr, "ori1    = %s\n", ori1 );
1266                 fprintf( stderr, "*mseq2  = %s\n", *mseq2 );
1267                 fprintf( stderr, "ungap2  = %s\n", ungap2 );
1268                 fprintf( stderr, "ori2    = %s\n", ori2 );
1269                 exit( 1 );
1270         }
1271
1272         value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
1273
1274 //      fclose( fp ); // saigo dake yatta houga yoi.
1275
1276         return( value );
1277 }
1278
1279
1280 static double calldafs_giving_bpp( char **mseq1, char **mseq2, char **bpp1, char **bpp2, int alloclen, int i, int j )
1281 {
1282         FILE *fp;
1283         int res;
1284         char *com;
1285         double value;
1286         char *dirname;
1287
1288
1289         dirname = calloc( 100, sizeof( char ) );
1290         com = calloc( 1000, sizeof( char ) );
1291         sprintf( dirname, "_%d-%d", i, j );
1292         sprintf( com, "rm -rf %s", dirname );
1293         system( com );
1294         sprintf( com, "mkdir %s", dirname );
1295         system( com );
1296
1297
1298         sprintf( com, "%s/_bpporg", dirname );
1299         fp = fopen( com, "w" );
1300         if( !fp )
1301         {
1302                 fprintf( stderr, "Cannot write to %s/_bpporg\n", dirname );
1303                 exit( 1 );
1304         }
1305         fprintf( fp, ">a\n" );
1306         while( *bpp1 )
1307                 fprintf( fp, "%s", *bpp1++ );
1308
1309         fprintf( fp, ">b\n" );
1310         while( *bpp2 )
1311                 fprintf( fp, "%s", *bpp2++ );
1312         fclose( fp );
1313
1314         sprintf( com, "tr -d '\\r' < %s/_bpporg > %s/_bpp", dirname, dirname );
1315         system( com ); // for cygwin, wakaran
1316
1317         t2u( *mseq1 );
1318         t2u( *mseq2 );
1319
1320         sprintf( com, "%s/_dafsinorg", dirname );
1321         fp = fopen( com, "w" );
1322         if( !fp )
1323         {
1324                 fprintf( stderr, "Cannot open %s/_dafsinorg\n", dirname );
1325                 exit( 1 );
1326         }
1327         fprintf( fp, ">1\n" );
1328 //      fprintf( fp, "%s\n", *mseq1 );
1329         write1seq( fp, *mseq1 );
1330         fprintf( fp, ">2\n" );
1331 //      fprintf( fp, "%s\n", *mseq2 );
1332         write1seq( fp, *mseq2 );
1333         fclose( fp );
1334
1335         sprintf( com, "tr -d '\\r' < %s/_dafsinorg > %s/_dafsin", dirname, dirname );
1336         system( com ); // for cygwin, wakaran
1337
1338         sprintf( com, "_dafssh%s", dirname );
1339         fp = fopen( com, "w" );
1340         fprintf( fp, "cd %s\n", dirname );
1341         fprintf( fp, "%s/dafs --mafft-in _bpp _dafsin > _dafsout 2>_dum\n", whereispairalign );
1342         fprintf( fp, "exit $tatus\n" );
1343         fclose( fp );
1344
1345         sprintf( com, "tr -d '\\r' < _dafssh%s > _dafssh%s.unix", dirname, dirname );
1346         system( com ); // for cygwin, wakaran
1347
1348         sprintf( com, "sh _dafssh%s.unix 2>_dum%s", dirname, dirname );
1349         res = system( com );
1350         if( res )
1351         {
1352                 fprintf( stderr, "Error in dafs\n" );
1353                 exit( 1 );
1354         }
1355
1356         sprintf( com, "%s/_dafsout", dirname );
1357
1358         fp = fopen( com, "r" );
1359         if( !fp )
1360         {
1361                 fprintf( stderr, "Cannot open %s/_dafsout\n", dirname );
1362                 exit( 1 );
1363         }
1364
1365         myfgets( com, 999, fp ); // nagai kanousei ga arunode
1366         fgets( com, 999, fp );
1367         myfgets( com, 999, fp ); // nagai kanousei ga arunode
1368         fgets( com, 999, fp );
1369         load1SeqWithoutName_new( fp, *mseq1 );
1370         fgets( com, 999, fp );
1371         load1SeqWithoutName_new( fp, *mseq2 );
1372
1373         fclose( fp );
1374
1375 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
1376 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
1377
1378         value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
1379
1380 #if 0
1381         sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
1382         if( system( com ) )
1383         {
1384                 fprintf( stderr, "retrying to rmdir\n" );
1385                 usleep( 2000 );
1386                 system( com );
1387         }
1388 #endif
1389
1390         free( dirname );
1391         free( com );
1392
1393
1394         return( value );
1395 }
1396
1397 static double callmxscarna_giving_bpp( char **mseq1, char **mseq2, char **bpp1, char **bpp2, int alloclen, int i, int j )
1398 {
1399         FILE *fp;
1400         int res;
1401         char *com;
1402         double value;
1403         char *dirname;
1404
1405
1406         dirname = calloc( 100, sizeof( char ) );
1407         com = calloc( 1000, sizeof( char ) );
1408         sprintf( dirname, "_%d-%d", i, j );
1409         sprintf( com, "rm -rf %s", dirname );
1410         system( com );
1411         sprintf( com, "mkdir %s", dirname );
1412         system( com );
1413
1414
1415         sprintf( com, "%s/_bpporg", dirname );
1416         fp = fopen( com, "w" );
1417         if( !fp )
1418         {
1419                 fprintf( stderr, "Cannot write to %s/_bpporg\n", dirname );
1420                 exit( 1 );
1421         }
1422         fprintf( fp, ">a\n" );
1423         while( *bpp1 )
1424                 fprintf( fp, "%s", *bpp1++ );
1425
1426         fprintf( fp, ">b\n" );
1427         while( *bpp2 )
1428                 fprintf( fp, "%s", *bpp2++ );
1429         fclose( fp );
1430
1431         sprintf( com, "tr -d '\\r' < %s/_bpporg > %s/_bpp", dirname, dirname );
1432         system( com ); // for cygwin, wakaran
1433
1434         t2u( *mseq1 );
1435         t2u( *mseq2 );
1436
1437         sprintf( com, "%s/_mxscarnainorg", dirname );
1438         fp = fopen( com, "w" );
1439         if( !fp )
1440         {
1441                 fprintf( stderr, "Cannot open %s/_mxscarnainorg\n", dirname );
1442                 exit( 1 );
1443         }
1444         fprintf( fp, ">1\n" );
1445 //      fprintf( fp, "%s\n", *mseq1 );
1446         write1seq( fp, *mseq1 );
1447         fprintf( fp, ">2\n" );
1448 //      fprintf( fp, "%s\n", *mseq2 );
1449         write1seq( fp, *mseq2 );
1450         fclose( fp );
1451
1452         sprintf( com, "tr -d '\\r' < %s/_mxscarnainorg > %s/_mxscarnain", dirname, dirname );
1453         system( com ); // for cygwin, wakaran
1454
1455 #if 0
1456         sprintf( com, "cd %s; %s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum", dirname, whereispairalign );
1457 #else
1458         sprintf( com, "_mxscarnash%s", dirname );
1459         fp = fopen( com, "w" );
1460         fprintf( fp, "cd %s\n", dirname );
1461         fprintf( fp, "%s/mxscarnamod -readbpp _mxscarnain > _mxscarnaout 2>_dum\n", whereispairalign );
1462         fprintf( fp, "exit $tatus\n" );
1463         fclose( fp );
1464 //sleep( 10000 );
1465
1466         sprintf( com, "tr -d '\\r' < _mxscarnash%s > _mxscarnash%s.unix", dirname, dirname );
1467         system( com ); // for cygwin, wakaran
1468
1469         sprintf( com, "sh _mxscarnash%s.unix 2>_dum%s", dirname, dirname );
1470 #endif
1471         res = system( com );
1472         if( res )
1473         {
1474                 fprintf( stderr, "Error in mxscarna\n" );
1475                 exit( 1 );
1476         }
1477
1478         sprintf( com, "%s/_mxscarnaout", dirname );
1479
1480         fp = fopen( com, "r" );
1481         if( !fp )
1482         {
1483                 fprintf( stderr, "Cannot open %s/_mxscarnaout\n", dirname );
1484                 exit( 1 );
1485         }
1486
1487         fgets( com, 999, fp );
1488         load1SeqWithoutName_new( fp, *mseq1 );
1489         fgets( com, 999, fp );
1490         load1SeqWithoutName_new( fp, *mseq2 );
1491
1492         fclose( fp );
1493
1494 //      fprintf( stderr, "*mseq1 = %s\n", *mseq1 );
1495 //      fprintf( stderr, "*mseq2 = %s\n", *mseq2 );
1496
1497         value = (double)naivepairscore11( *mseq1, *mseq2, penalty );
1498
1499 #if 0
1500         sprintf( com, "rm -rf %s > /dev/null 2>&1", dirname );
1501         if( system( com ) )
1502         {
1503                 fprintf( stderr, "retrying to rmdir\n" );
1504                 usleep( 2000 );
1505                 system( com );
1506         }
1507 #endif
1508
1509         free( dirname );
1510         free( com );
1511
1512
1513         return( value );
1514 }
1515
1516 static void readhat4( FILE *fp, char ***bpp )
1517 {
1518         char oneline[1000];
1519         int bppsize;
1520         int onechar;
1521 //      double prob;
1522 //      int posi, posj;
1523
1524         bppsize = 0;
1525 //      fprintf( stderr, "reading hat4\n" );
1526         onechar = getc(fp);
1527 //      fprintf( stderr, "onechar = %c\n", onechar );
1528         if( onechar != '>' )
1529         {
1530                 fprintf( stderr, "Format error\n" );
1531                 exit( 1 );
1532         }
1533         ungetc( onechar, fp );
1534         fgets( oneline, 999, fp );
1535         while( 1 )
1536         {
1537                 onechar = getc(fp);
1538                 ungetc( onechar, fp );
1539                 if( onechar == '>' || onechar == EOF )
1540                 {
1541 //                      fprintf( stderr, "Next\n" );
1542                         *bpp = realloc( *bpp, (bppsize+2) * sizeof( char * ) );
1543                         (*bpp)[bppsize] = NULL;
1544                         break;
1545                 }
1546                 fgets( oneline, 999, fp );
1547 //              fprintf( stderr, "oneline=%s\n", oneline );
1548 //              sscanf( oneline, "%d %d %lf", &posi, &posj, &prob );
1549 //              fprintf( stderr, "%d %d -> %f\n", posi, posj, prob );
1550                 *bpp = realloc( *bpp, (bppsize+2) * sizeof( char * ) );
1551                 (*bpp)[bppsize] = calloc( 100, sizeof( char ) );
1552                 strcpy( (*bpp)[bppsize], oneline );
1553                 bppsize++;
1554         }
1555 }
1556
1557 static void preparebpp( int nseq, char ***bpp )
1558 {
1559         FILE *fp;
1560         int i;
1561
1562         fp = fopen( "hat4", "r" );
1563         if( !fp )
1564         {
1565                 fprintf( stderr, "Cannot open hat4\n" );
1566                 exit( 1 );
1567         }
1568         for( i=0; i<nseq; i++ )
1569                 readhat4( fp, bpp+i );
1570         fclose( fp );
1571 }
1572
1573 static void arguments( int argc, char *argv[] )
1574 {
1575     int c;
1576
1577         nthread = 1;
1578         laste = 5000;
1579         lastm = 3;
1580         nadd = 0;
1581         lastsubopt = 0;
1582         lastonce = 0;
1583         foldalignopt[0] = 0;
1584         laraparams = NULL;
1585         inputfile = NULL;
1586         fftkeika = 0;
1587         pslocal = -1000.0;
1588         constraint = 0;
1589         nblosum = 62;
1590         fmodel = 0;
1591         calledByXced = 0;
1592         devide = 0;
1593         use_fft = 0;
1594         fftscore = 1;
1595         fftRepeatStop = 0;
1596         fftNoAnchStop = 0;
1597     weight = 3;
1598     utree = 1;
1599         tbutree = 1;
1600     refine = 0;
1601     check = 1;
1602     cut = 0.0;
1603     disp = 0;
1604     outgap = 1;
1605     alg = 'A';
1606     mix = 0;
1607         tbitr = 0;
1608         scmtd = 5;
1609         tbweight = 0;
1610         tbrweight = 3;
1611         checkC = 0;
1612         treemethod = 'x';
1613         contin = 0;
1614         scoremtx = 1;
1615         kobetsubunkatsu = 0;
1616         divpairscore = 0;
1617         stdout_align = 0;
1618         stdout_dist = 0;
1619         store_dist = 1;
1620         store_localhom = 1;
1621 //      dorp = NOTSPECIFIED;
1622         ppenalty = NOTSPECIFIED;
1623         ppenalty_OP = NOTSPECIFIED;
1624         ppenalty_ex = NOTSPECIFIED;
1625         ppenalty_EX = NOTSPECIFIED;
1626         penalty_shift_factor = 1000.0;
1627         poffset = NOTSPECIFIED;
1628         kimuraR = NOTSPECIFIED;
1629         pamN = NOTSPECIFIED;
1630         geta2 = GETA2;
1631         fftWinSize = NOTSPECIFIED;
1632         fftThreshold = NOTSPECIFIED;
1633         RNAppenalty = NOTSPECIFIED;
1634         RNApthr = NOTSPECIFIED;
1635         specificityconsideration = 0.0;
1636         usenaivescoreinsteadofalignmentscore = 0;
1637         specifictarget = 0;
1638         nwildcard = 0;
1639
1640 //      reporterr( "argc=%d\n", argc );
1641 //      reporterr( "*argv=%s\n", *argv );
1642 //      reporterr( "(*argv)[0]=%c\n", (*argv)[0] );
1643     while( --argc > 0 && (*++argv)[0] == '-' )
1644         {
1645 //              reporterr( "(*argv)[0] in while loop = %s\n", (*argv) );
1646         while ( ( c = *++argv[0] ) )
1647                 {
1648             switch( c )
1649             {
1650                                 case 'i':
1651                                         inputfile = *++argv;
1652 //                                      fprintf( stderr, "inputfile = %s\n", inputfile );
1653                                         --argc;
1654                                         goto nextoption;
1655                                 case 'f':
1656                                         ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
1657                                         --argc;
1658                                         goto nextoption;
1659                                 case 'g':
1660                                         ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
1661                                         --argc;
1662                                         goto nextoption;
1663                                 case 'O':
1664                                         ppenalty_OP = (int)( atof( *++argv ) * 1000 - 0.5 );
1665                                         --argc;
1666                                         goto nextoption;
1667                                 case 'E':
1668                                         ppenalty_EX = (int)( atof( *++argv ) * 1000 - 0.5 );
1669                                         --argc;
1670                                         goto nextoption;
1671                                 case 'Q':
1672                                         penalty_shift_factor = atof( *++argv );
1673                                         --argc;
1674                                         goto nextoption;
1675                                 case 'h':
1676                                         poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
1677                                         --argc;
1678                                         goto nextoption;
1679                                 case 'k':
1680                                         kimuraR = myatoi( *++argv );
1681 //                                      fprintf( stderr, "kimuraR = %d\n", kimuraR );
1682                                         --argc;
1683                                         goto nextoption;
1684                                 case 'b':
1685                                         nblosum = myatoi( *++argv );
1686                                         scoremtx = 1;
1687 //                                      fprintf( stderr, "blosum %d\n", nblosum );
1688                                         --argc;
1689                                         goto nextoption;
1690                                 case 'j':
1691                                         pamN = myatoi( *++argv );
1692                                         scoremtx = 0;
1693                                         TMorJTT = JTT;
1694 //                                      fprintf( stderr, "jtt %d\n", pamN );
1695                                         --argc;
1696                                         goto nextoption;
1697                                 case 'm':
1698                                         pamN = myatoi( *++argv );
1699                                         scoremtx = 0;
1700                                         TMorJTT = TM;
1701 //                                      fprintf( stderr, "TM %d\n", pamN );
1702                                         --argc;
1703                                         goto nextoption;
1704 #if 0
1705                                 case 'l':
1706                                         ppslocal = (int)( atof( *++argv ) * 1000 + 0.5 );
1707                                         pslocal = (int)( 600.0 / 1000.0 * ppslocal + 0.5);
1708 //                                      fprintf( stderr, "ppslocal = %d\n", ppslocal );
1709 //                                      fprintf( stderr, "pslocal = %d\n", pslocal );
1710                                         --argc;
1711                                         goto nextoption;
1712 #else
1713                                 case 'l':
1714                                         if( atof( *++argv ) < 0.00001 ) store_localhom = 0;
1715                                         --argc;
1716                                         goto nextoption;
1717 #endif
1718                                 case 'd':
1719                                         whereispairalign = *++argv;
1720                                         fprintf( stderr, "whereispairalign = %s\n", whereispairalign );
1721                                         --argc;
1722                                         goto nextoption;
1723                                 case 'p':
1724                                         laraparams = *++argv;
1725                                         fprintf( stderr, "laraparams = %s\n", laraparams );
1726                                         --argc;
1727                                         goto nextoption;
1728                                 case 'C':
1729                                         nthread = myatoi( *++argv );
1730 //                                      fprintf( stderr, "nthread = %d\n", nthread );
1731                                         --argc;
1732 #ifndef enablemultithread
1733                                         nthread = 0;
1734 #endif
1735                                         goto nextoption;
1736                                 case 'I':
1737                                         nadd = myatoi( *++argv );
1738 //                                      fprintf( stderr, "nadd = %d\n", nadd );
1739                                         --argc;
1740                                         goto nextoption;
1741                                 case 'w':
1742                                         lastm = myatoi( *++argv );
1743                                         fprintf( stderr, "lastm = %d\n", lastm );
1744                                         --argc;
1745                                         goto nextoption;
1746                                 case 'e':
1747                                         laste = myatoi( *++argv );
1748                                         fprintf( stderr, "laste = %d\n", laste );
1749                                         --argc;
1750                                         goto nextoption;
1751                                 case 'u':
1752                                         specificityconsideration = (double)myatof( *++argv );
1753 //                                      fprintf( stderr, "specificityconsideration = %f\n", specificityconsideration );
1754                                         --argc;
1755                                         goto nextoption;
1756                                 case 'K': // Hontou ha iranai. disttbfast.c, tbfast.c to awaserutame.
1757                                         break;
1758                                 case 'c':
1759                                         stdout_dist = 1;
1760                                         break;
1761                                 case 'n':
1762                                         stdout_align = 1;
1763                                         break;
1764                                 case 'x':
1765                                         store_localhom = 0;
1766                                         store_dist = 0;
1767                                         break;
1768 #if 1
1769                                 case 'a':
1770                                         fmodel = 1;
1771                                         break;
1772 #endif
1773 #if 0
1774                                 case 'r':
1775                                         fmodel = -1;
1776                                         break;
1777 #endif
1778                                 case 'D':
1779                                         dorp = 'd';
1780                                         break;
1781                                 case 'P':
1782                                         dorp = 'p';
1783                                         break;
1784 #if 0
1785                                 case 'e':
1786                                         fftscore = 0;
1787                                         break;
1788                                 case 'O':
1789                                         fftNoAnchStop = 1;
1790                                         break;
1791 #endif
1792 #if 0
1793                                 case 'Q':
1794                                         calledByXced = 1;
1795                                         break;
1796                                 case 'x':
1797                                         disp = 1;
1798                                         break;
1799                                 case 'a':
1800                                         alg = 'a';
1801                                         break;
1802                                 case 'S':
1803                                         alg = 'S';
1804                                         break;
1805 #endif
1806                                 case 'U':
1807                                         lastonce = 1;
1808                                         break;
1809                                 case 'S':
1810                                         lastsubopt = 1;
1811                                         break;
1812                                 case 't':
1813                                         alg = 't';
1814                                         store_localhom = 0;
1815                                         break;
1816                                 case 'L':
1817                                         alg = 'L';
1818                                         break;
1819                                 case 'Y':
1820                                         alg = 'Y'; // nadd>0 no toki nomi. moto no hairetsu to atarashii hairetsuno alignmnt -> L;
1821                                         break;
1822                                 case 'Z':
1823                                         usenaivescoreinsteadofalignmentscore = 1;
1824                                         break;
1825                                 case 's':
1826                                         alg = 's';
1827                                         break;
1828                                 case 'G':
1829                                         alg = 'G';
1830                                         break;
1831                                 case 'B':
1832                                         alg = 'B';
1833                                         break;
1834                                 case 'T':
1835                                         alg = 'T';
1836                                         break;
1837                                 case 'H':
1838                                         alg = 'H';
1839                                         break;
1840                                 case 'M':
1841                                         alg = 'M';
1842                                         break;
1843                                 case 'R':
1844                                         alg = 'R';
1845                                         break;
1846                                 case 'r':
1847                                         alg = 'r'; // nadd>0 no toki nomi. moto no hairetsu to atarashii hairetsuno alignmnt -> R, last
1848                                         break;
1849                                 case 'N':
1850                                         alg = 'N';
1851                                         break;
1852                                 case 'A':
1853                                         alg = 'A';
1854                                         break;
1855                                 case 'V':
1856                                         alg = 'V';
1857                                         break;
1858                                 case 'F':
1859                                         use_fft = 1;
1860                                         break;
1861                                 case 'v':
1862                                         tbrweight = 3;
1863                                         break;
1864                                 case 'y':
1865                                         divpairscore = 1;
1866                                         break;
1867                                 case '=':
1868                                         specifictarget = 1;
1869                                         break;
1870                                 case ':':
1871                                         nwildcard = 1;
1872                                         break;
1873 /* Modified 01/08/27, default: user tree */
1874                                 case 'J':
1875                                         tbutree = 0;
1876                                         break;
1877 /* modification end. */
1878                                 case 'o':
1879 //                                      foldalignopt = *++argv;
1880                                         strcat( foldalignopt, " " );
1881                                         strcat( foldalignopt, *++argv );
1882                                         fprintf( stderr, "foldalignopt = %s\n", foldalignopt );
1883                                         --argc;
1884                                         goto nextoption;
1885 #if 0
1886                                 case 'z':
1887                                         fftThreshold = myatoi( *++argv );
1888                                         --argc;
1889                                         goto nextoption;
1890                                 case 'w':
1891                                         fftWinSize = myatoi( *++argv );
1892                                         --argc;
1893                                         goto nextoption;
1894                                 case 'Z':
1895                                         checkC = 1;
1896                                         break;
1897 #endif
1898                 default:
1899                     fprintf( stderr, "illegal option %c\n", c );
1900                     argc = 0;
1901                     break;
1902             }
1903                 }
1904                 nextoption:
1905                         ;
1906         }
1907     if( argc == 1 )
1908     {
1909         cut = atof( (*argv) );
1910         argc--;
1911     }
1912     if( argc != 0 )
1913     {
1914         fprintf( stderr, "pairlocalalign options: Check source file !\n" );
1915         exit( 1 );
1916     }
1917         if( tbitr == 1 && outgap == 0 )
1918         {
1919                 fprintf( stderr, "conflicting options : o, m or u\n" );
1920                 exit( 1 );
1921         }
1922 }
1923
1924 int countamino( char *s, int end )
1925 {
1926         int val = 0;
1927         while( end-- )
1928                 if( *s++ != '-' ) val++;
1929         return( val );
1930 }
1931
1932 static double score2dist( double pscore, double selfscore1, double selfscore2)
1933 {
1934         double val;
1935         double bunbo;
1936 //      fprintf( stderr, "In score2dist\n" );
1937
1938         if( (bunbo=MIN( selfscore1, selfscore2 )) == 0.0 )
1939                 val = 2.0;
1940         else if( bunbo < pscore ) // mondai ari
1941                 val = 0.0;
1942         else
1943                 val = ( 1.0 - pscore / bunbo ) * 2.0;
1944         return( val );
1945 }
1946
1947 #if enablemultithread
1948 static void *athread( void *arg ) // alg='R', alg='r' -> tsukawarenai.
1949 {
1950         thread_arg_t *targ = (thread_arg_t *)arg;
1951         int i, ilim, j, jst;
1952         int off1, off2, dum1, dum2, thereisx;
1953         int intdum;
1954         double pscore = 0.0; // by D.Mathog
1955         double *effarr1;
1956         double *effarr2;
1957         char **mseq1, **mseq2, **distseq1, **distseq2, **dumseq1, **dumseq2;
1958         char **aseq;
1959         double **dynamicmtx = NULL;
1960         double dist;
1961         double scoreoffset;
1962
1963 // thread_arg
1964         int thread_no = targ->thread_no;
1965         int njob = targ->njob;
1966         Jobtable *jobpospt = targ->jobpospt;
1967         char **name = targ->name;
1968         char **seq = targ->seq;
1969         char **dseq = targ->dseq;
1970         int *thereisxineachseq = targ->thereisxineachseq;
1971         LocalHom **localhomtable = targ->localhomtable;
1972         double **distancemtx = targ->distancemtx;
1973         double *selfscore = targ->selfscore;
1974         char ***bpp = targ->bpp;
1975         Lastresx **lastresx = targ->lastresx;
1976         int alloclen = targ->alloclen;
1977         int *targetmap = targ->targetmap;
1978
1979 //      fprintf( stderr, "thread %d start!\n", thread_no );
1980
1981         effarr1 = AllocateDoubleVec( 1 );
1982         effarr2 = AllocateDoubleVec( 1 );
1983         mseq1 = AllocateCharMtx( njob, 0 );
1984         mseq2 = AllocateCharMtx( njob, 0 );
1985         if( alg == 'N' )
1986         {
1987                 dumseq1 = AllocateCharMtx( 1, alloclen+10 );
1988                 dumseq2 = AllocateCharMtx( 1, alloclen+10 );
1989         }
1990         distseq1 = AllocateCharMtx( 1, 0 );
1991         distseq2 = AllocateCharMtx( 1, 0 );
1992         aseq = AllocateCharMtx( 2, alloclen+10 );
1993         if( specificityconsideration > 0.0 ) dynamicmtx = AllocateDoubleMtx( nalphabets, nalphabets );
1994
1995         if( alg == 'Y' || alg == 'r' ) ilim = njob - nadd;
1996         else ilim = njob - 1;
1997
1998
1999         while( 1 )
2000         {
2001                 pthread_mutex_lock( targ->mutex_counter );
2002                 j = jobpospt->j;
2003                 i = jobpospt->i;
2004                 j++;
2005                 if( j == njob )
2006                 {
2007                         i++;
2008
2009                         if( alg == 'Y' || alg == 'r' ) jst = njob - nadd;
2010                         else jst = i + 1;
2011                         j = jst;
2012
2013                         if( i == ilim )
2014                         {
2015 //                              fprintf( stderr, "thread %d end!\n", thread_no );
2016                                 pthread_mutex_unlock( targ->mutex_counter );
2017
2018                                 if( commonIP ) FreeIntMtx( commonIP );
2019                                 commonIP = NULL;
2020                                 if( commonJP ) FreeIntMtx( commonJP );
2021                                 commonJP = NULL;
2022                                 Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
2023                                 G__align11( NULL, NULL, NULL, 0, 0, 0 ); // 20130603
2024                                 G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
2025                                 L__align11( NULL, 0.0, NULL, NULL, 0, NULL, NULL );
2026                                 L__align11_noalign( NULL, NULL, NULL );
2027                                 genL__align11( NULL, NULL, NULL, 0, NULL, NULL );
2028                                 free( effarr1 );
2029                                 free( effarr2 );
2030                                 free( mseq1 );
2031                                 free( mseq2 );
2032                                 if( alg == 'N' )
2033                                 {
2034                                         FreeCharMtx( dumseq1 );
2035                                         FreeCharMtx( dumseq2 );
2036                                 }
2037                                 free( distseq1 );
2038                                 free( distseq2 );
2039                                 FreeCharMtx( aseq  );
2040                                 if( dynamicmtx ) FreeDoubleMtx( dynamicmtx  );
2041                                 return( NULL );
2042                         }
2043                 }
2044                 jobpospt->j = j;
2045                 jobpospt->i = i;
2046                 pthread_mutex_unlock( targ->mutex_counter );
2047
2048
2049 //              if( j == i+1 || j % 100 == 0 )
2050                 if( j == i+1 && i % 10 == 0 )
2051                 {
2052                         fprintf( stderr, "% 5d / %d (by thread %3d) \r", i, njob-nadd, thread_no );
2053 //                      fprintf( stderr, "% 5d - %5d / %d (thread %d)\n", i, j, njob, thread_no );
2054                 }
2055
2056
2057                 if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
2058                 {
2059                         if( store_dist )
2060                         {
2061                                 if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = 3.0;
2062                                 else distancemtx[i][j-i] = 3.0;
2063                         }
2064                         if( stdout_dist)
2065                         {
2066                                 pthread_mutex_lock( targ->mutex_stdout );
2067                                 fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 3.0 );
2068                                 pthread_mutex_unlock( targ->mutex_stdout );
2069                         }
2070                         continue;
2071                 }
2072
2073                 strcpy( aseq[0], seq[i] );
2074                 strcpy( aseq[1], seq[j] );
2075 //              clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
2076 //              clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
2077 //              fprintf( stderr, "Skipping conjuction..\n" );
2078
2079                 effarr1[0] = 1.0;
2080                 effarr2[0] = 1.0;
2081                 mseq1[0] = aseq[0];
2082                 mseq2[0] = aseq[1];
2083
2084                 thereisx = thereisxineachseq[i] + thereisxineachseq[j];
2085 //              strcpy( distseq1[0], dseq[i] ); // nen no tame
2086 //              strcpy( distseq2[0], dseq[j] ); // nen no tame
2087                 distseq1[0] = dseq[i];
2088                 distseq2[0] = dseq[j];
2089
2090 //              fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
2091 //              fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
2092
2093 #if 0
2094                 fprintf( stderr, "group1 = %.66s", indication1 );
2095                 fprintf( stderr, "\n" );
2096                 fprintf( stderr, "group2 = %.66s", indication2 );
2097                 fprintf( stderr, "\n" );
2098 #endif
2099 //              for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
2100
2101                 if( use_fft )
2102                 {
2103                         pscore = Falign( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, 1, 1, alloclen, &intdum, NULL, 0, NULL );
2104 //                      fprintf( stderr, "pscore (fft) = %f\n", pscore );
2105                         off1 = off2 = 0;
2106                 }
2107                 else
2108                 {
2109                         switch( alg )
2110                         {
2111                                 case( 'R' ):
2112                                         if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2113                                                 pscore = 0.0;
2114                                         else
2115                                                 pscore = (double)lastresx[i][j].score; // all pair
2116                                         break;
2117                                 case( 'r' ):
2118                                         if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) )
2119                                                 pscore = (double)lastresx[i][j-(njob-nadd)].score;
2120                                         else
2121                                                 pscore = 0.0;
2122                                         break;
2123                                 case( 'L' ):
2124                                         if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2125                                                 pscore = 0.0;
2126                                         else
2127                                         {
2128                                                 if( usenaivescoreinsteadofalignmentscore )
2129                                                 {
2130                                                         L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2131                                                         pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2132                                                 }
2133                                                 else
2134                                                 {
2135 //                                                      if( store_localhom )
2136                                                         if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
2137                                                         {
2138                                                                 pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2139                                                                 if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // uwagaki
2140 #if 1
2141                                                                 if( specificityconsideration > 0.0 )
2142                                                                 {
2143                                                                         dist = score2dist( pscore, selfscore[i], selfscore[j] );
2144                                                                         if( ( scoreoffset = dist2offset( dist ) ) < 0.0 )
2145                                                                         {
2146                                                                                 makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
2147                                                                                 strcpy( mseq1[0], seq[i] );
2148                                                                                 strcpy( mseq2[0], seq[j] );
2149                                                                                 L__align11( dynamicmtx, scoreoffset, mseq1, mseq2, alloclen, &off1, &off2 );
2150                                                                         }
2151                                                                 }
2152 #endif
2153                                                         }
2154                                                         else
2155                                                                 pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 );
2156                                                 }
2157                                         }
2158 //                                      pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // CHUUI!!!!!!
2159                                         break;
2160                                 case( 'Y' ):
2161                                         if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) ) // new sequence vs exiting sequence nomi keisan
2162                                         {
2163                                                 if( usenaivescoreinsteadofalignmentscore )
2164                                                 {
2165                                                         L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2166                                                         pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2167                                                 }
2168                                                 else
2169                                                 {
2170                                                         if( store_localhom )
2171                                                         {
2172                                                                 pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2173                                                                 if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // uwagaki
2174                                                         }
2175                                                         else
2176                                                                 pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 );
2177                                                 }
2178                                         }
2179                                         else
2180                                                 pscore = 0.0;
2181                                         break;
2182                                 case( 'A' ):
2183                                         if( usenaivescoreinsteadofalignmentscore )
2184                                         {
2185                                                 G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
2186                                                 pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2187                                         }
2188                                         else
2189                                         {
2190 //                                              if( store_localhom )
2191                                                 if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
2192                                                 {
2193                                                         pscore = G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
2194                                                         if( thereisx ) pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2195 #if 1
2196                                                         if( specificityconsideration > 0.0 )
2197                                                         {
2198                                                                 dist = score2dist( pscore, selfscore[i], selfscore[j] );
2199 //                                                              dist = score2dist( L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ), selfscore[i], selfscore[j] ); // 2014/Feb/20
2200                                                                 if( dist2offset( dist ) < 0.0 )
2201                                                                 {
2202                                                                         makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
2203                                                                         strcpy( mseq1[0], seq[i] );
2204                                                                         strcpy( mseq2[0], seq[j] );
2205                                                                         G__align11( dynamicmtx, mseq1, mseq2, alloclen, outgap, outgap );
2206
2207                                                                 }
2208 //                                                              pscore = (double)naivepairscore11( *mseq1, *mseq2, 0.0 );
2209                                                         }
2210 #endif
2211                                                 }
2212                                                 else
2213                                                         pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2214                                         }
2215                                         off1 = off2 = 0;
2216                                         break;
2217                                 case( 'N' ):
2218                                         if( usenaivescoreinsteadofalignmentscore )
2219                                         {
2220                                                 genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
2221                                                 pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2222                                         }
2223                                         else
2224                                         {
2225 //                                              pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
2226                                                 pscore = genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
2227                                                 if( thereisx )
2228                                                 {
2229                                                         strcpy( dumseq1[0], distseq1[0] );
2230                                                         strcpy( dumseq2[0], distseq2[0] );
2231                                                         pscore = genL__align11( n_dis_consweight_multi, dumseq1, dumseq2, alloclen, &dum1, &dum2 ); // uwagaki
2232                                                 }
2233 #if 1
2234                                                 if( specificityconsideration > 0.0 )
2235                                                 {
2236                                                         dist = score2dist( pscore, selfscore[i], selfscore[j] );
2237                                                         if( dist2offset( dist ) < 0.0 )
2238                                                         {
2239                                                                 makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
2240                                                                 strcpy( mseq1[0], seq[i] );
2241                                                                 strcpy( mseq2[0], seq[j] );
2242                                                                 genL__align11( dynamicmtx, mseq1, mseq2, alloclen, &off1, &off2 );
2243                                                         }
2244                                                 }
2245 #endif
2246                                         }
2247                                         break;
2248                                 case( 't' ):
2249                                         off1 = off2 = 0;
2250 //                                      pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
2251                                         pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // tsuneni distseq shiyou
2252                                         break;
2253                                 case( 's' ):
2254                                         pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2255                                         off1 = off2 = 0;
2256                                         break;
2257                                 case( 'G' ):
2258                                         pscore = calldafs_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2259                                         off1 = off2 = 0;
2260                                         break;
2261 #if 0
2262                                 case( 'a' ):
2263                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
2264                                         off1 = off2 = 0;
2265                                         break;
2266                                 case( 'K' ):
2267                                         pscore = genG__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen );
2268                                         off1 = off2 = 0;
2269                                         break;
2270                                 case( 'H' ):
2271                                         pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
2272                                         break;
2273                                 case( 'B' ):
2274                                 case( 'T' ):
2275                                         pscore = recalllara( mseq1, mseq2, alloclen );
2276                                         off1 = off2 = 0;
2277                                         break;
2278                                 case( 'M' ):
2279                                         pscore = MSalign11( mseq1, mseq2, alloclen );
2280                                         break;
2281 #endif
2282                                 default:
2283                                         ErrorExit( "\n\nERROR IN SOURCE FILE\n\n" );
2284                         }
2285                 }
2286
2287                 if( alg == 't' || ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) ) // 't' no jouken ha iranai to omou. if( ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) )
2288                 {
2289 #if SCOREOUT
2290                         fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
2291 #endif
2292 //                      if( pscore > 0.0 && ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) ) x-ins-i de seido teika
2293                         if( ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) )
2294                         {
2295                                 if( !store_localhom )
2296                                         ;
2297                                 else if( specifictarget && targetmap[i] == -1 && targetmap[j] == -1)
2298                                         ;
2299                                 else if( alg == 'R' )
2300                                         putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j, lastresx[i]+j, 'h' );
2301                                 else if( alg == 'r' )
2302                                         putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), lastresx[i]+j-(njob-nadd), 'h' );// ?????
2303                                 else if( alg == 'H' )
2304                                         putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2305                                 else if( alg == 'Y' )
2306                                         putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2307                                 else if( !specifictarget && alg != 'S' && alg != 'V' )
2308                                         putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-i, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2309                                 else
2310 //                                      putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ) );
2311                                 {
2312                                         if( targetmap[i] != -1 && targetmap[j] != -1 )
2313                                         {
2314                                                 putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
2315                                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' ); // sukoshi muda.
2316                                         }
2317                                         else if( targetmap[j] != -1 )
2318                                                 putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
2319                                         else if( targetmap[i] != -1 )
2320                                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2321 #if 0
2322                                         if( targetmap[i] != -1 )
2323                                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2324
2325                                         else if( targetmap[j] != -1 )
2326                                                 putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
2327 #endif
2328                                         else
2329                                         {
2330                                                 reporterr( "okashii\n" );
2331                                                 exit( 1 );
2332                                         }
2333                                 }
2334                         }
2335                         pscore = score2dist( pscore, selfscore[i], selfscore[j] );
2336
2337 //                      pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 );
2338 //                      pscore = score2dist( pscore, selfscore[i], selfscore[j] );
2339 //                      reporterr( "->pscore = %f\n", pscore );
2340
2341                 }
2342                 else
2343                 {
2344                         pscore = 2.0;
2345                 }
2346
2347 #if 1 // mutex
2348                 if( stdout_align )
2349                 {
2350                         pthread_mutex_lock( targ->mutex_stdout );
2351                         if( alg != 't' )
2352                         {
2353                                 fprintf( stdout, "sequence %d - sequence %d, pairwise distance = %10.5f\n", i+1, j+1, pscore );
2354                                 fprintf( stdout, ">%s\n", name[i] );
2355                                 write1seq( stdout, mseq1[0] );
2356                                 fprintf( stdout, ">%s\n", name[j] );
2357                                 write1seq( stdout, mseq2[0] );
2358                                 fprintf( stdout, "\n" );
2359                         }
2360                         pthread_mutex_unlock( targ->mutex_stdout );
2361                 }
2362                 if( stdout_dist )
2363                 {
2364                         pthread_mutex_lock( targ->mutex_stdout );
2365                         if( j == i+1 ) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
2366                         fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
2367                         pthread_mutex_unlock( targ->mutex_stdout );
2368                 }
2369 #endif // mutex
2370                 if( store_dist )
2371                 {
2372                         if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = pscore;
2373                         else distancemtx[i][j-i] = pscore;
2374                 }
2375         }
2376 }
2377 #endif
2378
2379 static void pairalign( char **name, int *nlen, char **seq, char **aseq, char **dseq, int *thereisxineachseq, char **mseq1, char **mseq2, int alloclen, Lastresx **lastresx, double **distancemtx, LocalHom **localhomtable, int ngui )
2380 {
2381         int i, j, ilim, jst, jj;
2382         int off1, off2, dum1, dum2, thereisx;
2383         double pscore = 0.0; // by D.Mathog
2384         FILE *hat2p, *hat3p;
2385 //      double **distancemtx;
2386         double *selfscore;
2387         double *effarr1;
2388         double *effarr2;
2389         char *pt;
2390         char *hat2file = "hat2";
2391 //      LocalHom **localhomtable = NULL,
2392         LocalHom *tmpptr;
2393         int intdum;
2394         char ***bpp = NULL; // mxscarna no toki dake
2395         char **distseq1, **distseq2;
2396         char **dumseq1, **dumseq2;
2397         double dist;
2398         double scoreoffset;
2399         int ntarget;
2400         int *targetmap, *targetmapr;
2401
2402
2403         if( specifictarget )
2404         {
2405                 targetmap = calloc( njob, sizeof( int ) );
2406                 ntarget = 0;
2407                 for( i=0; i<njob; i++ )
2408                 {
2409                         targetmap[i] = -1;
2410                         if( !strncmp( name[i]+1, "_focus_", 7 ) )
2411                                 targetmap[i] = ntarget++;
2412                 }
2413                 targetmapr = calloc( ntarget, sizeof( int ) );
2414                 for( i=0; i<njob; i++ )
2415                         if( targetmap[i] != -1 ) targetmapr[targetmap[i]] = i;
2416
2417                 if( ntarget == 0 )
2418                 {
2419                         reporterr( "\n\nAdd '>_focus_' to the title lines of the sequences to be focused on.\n\n" );
2420                         exit( 1 );
2421                 }
2422                 else
2423                 {
2424                         reporterr( "nfocus = %d \n", ntarget );
2425                 }
2426         }
2427         else
2428         {
2429                 ntarget = njob;
2430                 targetmap = calloc( njob, sizeof( int ) );
2431                 targetmapr = calloc( njob, sizeof( int ) );
2432                 for( i=0; i<njob; i++ )
2433                         targetmap[i] = targetmapr[i] = i;
2434         }
2435
2436 #if 0
2437         for( i=0; i<njob; i++ )
2438                 reporterr( "targetmap[%d] = %d\n", i, targetmap[i] );
2439         for( i=0; i<ntarget; i++ )
2440                 reporterr( "targetmapr[%d] = %d\n", i, targetmapr[i] );
2441 #endif
2442
2443         if( store_localhom && localhomtable == NULL )
2444         {
2445                 if( alg == 'Y' || alg == 'r' )
2446                 {
2447                         ilim = njob - nadd;
2448                         jst = nadd;
2449                 }
2450                 else
2451                 {
2452                         ilim = ntarget;
2453                         jst = njob;
2454                 }
2455                 localhomtable = (LocalHom **)calloc( ilim, sizeof( LocalHom *) );
2456                 for( i=0; i<ilim; i++)
2457                 {
2458                         localhomtable[i] = (LocalHom *)calloc( jst, sizeof( LocalHom ) );
2459                         for( j=0; j<jst; j++)
2460                         {
2461                                 localhomtable[i][j].start1 = -1;
2462                                 localhomtable[i][j].end1 = -1;
2463                                 localhomtable[i][j].start2 = -1;
2464                                 localhomtable[i][j].end2 = -1;
2465                                 localhomtable[i][j].opt = -1.0;
2466                                 localhomtable[i][j].next = NULL;
2467                                 localhomtable[i][j].nokori = 0;
2468                                 localhomtable[i][j].extended = -1;
2469                                 localhomtable[i][j].last = localhomtable[i]+j;
2470                                 localhomtable[i][j].korh = 'h';
2471                         }
2472                         if( !specifictarget && alg != 'Y' && alg != 'r' ) jst--;
2473                 }
2474         }
2475
2476         if( store_dist )
2477         {
2478                 if( ngui == 0 )
2479                 {
2480                         if( alg == 'Y' || alg == 'r' )
2481                                 distancemtx = AllocateDoubleMtx( njob, nadd );
2482                         else
2483                                 distancemtx = AllocateDoubleHalfMtx( njob );
2484 //                              distancemtx = AllocateDoubleMtx( njob, njob );
2485                 }
2486         }
2487         else distancemtx = NULL;
2488
2489         if( alg == 'N' )
2490         {
2491                 dumseq1 = AllocateCharMtx( 1, alloclen+10 );
2492                 dumseq2 = AllocateCharMtx( 1, alloclen+10 );
2493         }
2494         distseq1 = AllocateCharMtx( 1, 0 ); // muda
2495         distseq2 = AllocateCharMtx( 1, 0 ); // muda
2496
2497         selfscore = AllocateDoubleVec( njob );
2498         effarr1 = AllocateDoubleVec( njob );
2499         effarr2 = AllocateDoubleVec( njob );
2500
2501 #if 0
2502         fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
2503 #endif
2504
2505 #if 0
2506         for( i=0; i<njob; i++ )
2507                 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
2508 #endif
2509
2510
2511 //      writePre( njob, name, nlen, aseq, 0 );
2512
2513         reporterr( "All-to-all alignment.\n" );
2514         if( alg == 'R' )
2515         {
2516                 fprintf( stderr, "Calling last (http://last.cbrc.jp/)\n" );
2517                 if( lastonce )
2518                         calllast_once( njob, seq, njob, seq, lastresx );
2519                 else
2520                         calllast_fast( njob, seq, njob, seq, lastresx );
2521                 fprintf( stderr, "done.\n" );
2522 //              nthread = 0; // igo multithread nashi
2523         }
2524         if( alg == 'r' )
2525         {
2526                 fprintf( stderr, "Calling last (http://last.cbrc.jp/)\n" );
2527                 fprintf( stderr, "nadd=%d\n", nadd );
2528 #if 1 // last_fast ha, lastdb ga muda
2529                 if( lastonce )
2530                         calllast_once( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
2531                 else
2532                         calllast_fast( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
2533 #else
2534                 calllast_once( njob-nadd, seq, nadd, seq+njob-nadd, lastresx );
2535 #endif
2536
2537                 fprintf( stderr, "nadd=%d\n", nadd );
2538                 fprintf( stderr, "done.\n" );
2539 //              nthread = 0; // igo multithread nashi
2540         }
2541
2542         if( alg == 'H' )
2543         {
2544                 fprintf( stderr, "Calling FOLDALIGN with option '%s'\n", foldalignopt );
2545                 callfoldalign( njob, seq );
2546                 fprintf( stderr, "done.\n" );
2547         }
2548         if( alg == 'B' )
2549         {
2550                 fprintf( stderr, "Running LARA (Bauer et al. http://www.planet-lisa.net/)\n" );
2551                 calllara( njob, seq, "" );
2552                 fprintf( stderr, "done.\n" );
2553         }
2554         if( alg == 'T' )
2555         {
2556                 fprintf( stderr, "Running SLARA (Bauer et al. http://www.planet-lisa.net/)\n" );
2557                 calllara( njob, seq, "-s" );
2558                 fprintf( stderr, "done.\n" );
2559         }
2560         if( alg == 's' )
2561         {
2562                 fprintf( stderr, "Preparing bpp\n" );
2563 //              bpp = AllocateCharCub( njob, nlenmax, 0 );
2564                 bpp = calloc( njob, sizeof( char ** ) );
2565                 preparebpp( njob, bpp );
2566                 fprintf( stderr, "done.\n" );
2567                 fprintf( stderr, "Running MXSCARNA (Tabei et al. http://www.ncrna.org/software/mxscarna)\n" );
2568         }
2569         if( alg == 'G' )
2570         {
2571                 fprintf( stderr, "Preparing bpp\n" );
2572 //              bpp = AllocateCharCub( njob, nlenmax, 0 );
2573                 bpp = calloc( njob, sizeof( char ** ) );
2574                 preparebpp( njob, bpp );
2575                 fprintf( stderr, "done.\n" );
2576                 fprintf( stderr, "Running DAFS (Sato et al. http://www.ncrna.org/)\n" );
2577         }
2578
2579         for( i=0; i<njob; i++ )
2580         {
2581                 pscore = 0.0;
2582                 for( pt=seq[i]; *pt; pt++ )
2583                         pscore += amino_dis[(unsigned char)*pt][(unsigned char)*pt];
2584                 selfscore[i] = pscore;
2585 //              fprintf( stderr, "selfscore[%d] = %f\n", i, selfscore[i] );
2586         }
2587
2588 #if enablemultithread
2589         if( nthread > 0 ) // alg=='r' || alg=='R' -> nthread:=0 (sukoshi ue)
2590         {
2591                 Jobtable jobpos;
2592                 pthread_t *handle;
2593                 pthread_mutex_t mutex_counter;
2594                 pthread_mutex_t mutex_stdout;
2595                 thread_arg_t *targ;
2596
2597                 if( alg == 'Y' || alg == 'r' ) jobpos.j = njob - nadd - 1;
2598                 else jobpos.j = 0;
2599                 jobpos.i = 0;
2600
2601                 targ = calloc( nthread, sizeof( thread_arg_t ) );
2602                 handle = calloc( nthread, sizeof( pthread_t ) );
2603                 pthread_mutex_init( &mutex_counter, NULL );
2604                 pthread_mutex_init( &mutex_stdout, NULL );
2605
2606                 for( i=0; i<nthread; i++ )
2607                 {
2608                         targ[i].thread_no = i;
2609                         targ[i].njob = njob;
2610                         targ[i].jobpospt = &jobpos;
2611                         targ[i].name = name;
2612                         targ[i].seq = seq;
2613                         targ[i].dseq = dseq;
2614                         targ[i].thereisxineachseq = thereisxineachseq;
2615                         targ[i].localhomtable = localhomtable;
2616                         targ[i].distancemtx = distancemtx;
2617                         targ[i].selfscore = selfscore;
2618                         targ[i].bpp = bpp;
2619                         targ[i].lastresx = lastresx;
2620                         targ[i].alloclen = alloclen;
2621                         targ[i].targetmap = targetmap;
2622                         targ[i].mutex_counter = &mutex_counter;
2623                         targ[i].mutex_stdout = &mutex_stdout;
2624
2625 //                      athread( (void *)targ );
2626                         pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
2627 //                      pthread_create( handle+i, NULL, bthread, (void *)(targ+i) );
2628                 }
2629
2630
2631                 for( i=0; i<nthread; i++ )
2632                 {
2633                         pthread_join( handle[i], NULL );
2634                 }
2635                 pthread_mutex_destroy( &mutex_counter );
2636                 pthread_mutex_destroy( &mutex_stdout );
2637                 free( handle );
2638                 free( targ );
2639         }
2640         else
2641 #endif
2642         {
2643                 double **dynamicmtx = NULL;
2644                 if( specificityconsideration > 0.0 ) dynamicmtx = AllocateDoubleMtx( nalphabets, nalphabets );
2645
2646                 if( alg == 'Y' || alg == 'r' ) ilim = njob - nadd;
2647                 else ilim = njob - 1;
2648                 for( i=0; i<ilim; i++ )
2649                 {
2650                         if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, i+1, 0.0 );
2651                         fprintf( stderr, "% 5d / %d\r", i, njob-nadd );
2652                         fflush( stderr );
2653
2654                         if( alg == 'Y' || alg == 'r' ) jst = njob - nadd;
2655                         else jst = i + 1;
2656                         for( j=jst; j<njob; j++ )
2657                         {
2658
2659                                 if( strlen( seq[i] ) == 0 || strlen( seq[j] ) == 0 )
2660                                 {
2661                                         if( store_dist )
2662                                         {
2663                                                 if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = 3.0;
2664                                                 else distancemtx[i][j-i] = 3.0;
2665                                         }
2666                                         if( stdout_dist) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, 3.0 );
2667                                         continue;
2668                                 }
2669
2670                                 strcpy( aseq[0], seq[i] );
2671                                 strcpy( aseq[1], seq[j] );
2672 //                              clus1 = conjuctionfortbfast( pair, i, aseq, mseq1, effarr1, effarr, indication1 );
2673 //                              clus2 = conjuctionfortbfast( pair, j, aseq, mseq2, effarr2, effarr, indication2 );
2674 //                              fprintf( stderr, "Skipping conjuction..\n" );
2675
2676                                 effarr1[0] = 1.0;
2677                                 effarr2[0] = 1.0;
2678                                 mseq1[0] = aseq[0];
2679                                 mseq2[0] = aseq[1];
2680
2681                                 thereisx = thereisxineachseq[i] + thereisxineachseq[j];
2682 //                              strcpy( distseq1[0], dseq[i] ); // nen no tame
2683 //                              strcpy( distseq2[0], dseq[j] ); // nen no tame
2684                                 distseq1[0] = dseq[i];
2685                                 distseq2[0] = dseq[j];
2686
2687         //                      fprintf( stderr, "mseq1 = %s\n", mseq1[0] );
2688         //                      fprintf( stderr, "mseq2 = %s\n", mseq2[0] );
2689
2690 #if 0
2691                                 fprintf( stderr, "group1 = %.66s", indication1 );
2692                                 fprintf( stderr, "\n" );
2693                                 fprintf( stderr, "group2 = %.66s", indication2 );
2694                                 fprintf( stderr, "\n" );
2695 #endif
2696         //                      for( l=0; l<clus1; l++ ) fprintf( stderr, "## STEP-eff for mseq1-%d %f\n", l, effarr1[l] );
2697
2698                                 if( use_fft )
2699                                 {
2700                                         pscore = Falign( NULL, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, 1, 1, alloclen, &intdum, NULL, 0, NULL );
2701 //                                      fprintf( stderr, "pscore (fft) = %f\n", pscore );
2702                                         off1 = off2 = 0;
2703                                 }
2704                                 else
2705                                 {
2706                                         switch( alg )
2707                                         {
2708                                                 case( 't' ):
2709 //                                                      pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
2710                                                         pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // tsuneni distseq shiyou
2711                                                         off1 = off2 = 0;
2712                                                         break;
2713                                                 case( 'A' ):
2714                                                         if( usenaivescoreinsteadofalignmentscore )
2715                                                         {
2716                                                                 G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
2717                                                                 pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2718                                                         }
2719                                                         else
2720                                                         {
2721 //                                                              if( store_localhom )
2722                                                                 if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
2723                                                                 {
2724                                                                         pscore = G__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, outgap, outgap );
2725                                                                         if( thereisx ) pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2726 #if 1
2727                                                                         if( specificityconsideration > 0.0 )
2728                                                                         {
2729                                                                                 dist = score2dist( pscore, selfscore[i], selfscore[j] );
2730 //                                                                              dist = score2dist( L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ), selfscore[i], selfscore[j] ); // 2014/Feb/20
2731                                                                                 if( dist2offset( dist ) < 0.0 )
2732                                                                                 {
2733                                                                                         makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
2734                                                                                         strcpy( mseq1[0], seq[i] );
2735                                                                                         strcpy( mseq2[0], seq[j] );
2736                                                                                         G__align11( dynamicmtx, mseq1, mseq2, alloclen, outgap, outgap );
2737                                                                                 }
2738 //                                                                              pscore = (double)naivepairscore11( *mseq1, *mseq2, 0.0 );
2739                                                                         }
2740 #endif
2741                                                                 }
2742                                                                 else
2743                                                                         pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // uwagaki
2744                                                         }
2745                                                         off1 = off2 = 0;
2746                                                         break;
2747                                                 case( 'N' ):
2748 //                                                      pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, alloclen );
2749                                                         if( usenaivescoreinsteadofalignmentscore )
2750                                                         {
2751                                                                 genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
2752                                                                 pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2753                                                         }
2754                                                         else
2755                                                         {
2756                                                                 pscore = genL__align11( n_dis_consweight_multi, mseq1, mseq2, alloclen, &off1, &off2 );
2757                                                                 if( thereisx )
2758                                                                 {
2759                                                                         strcpy( dumseq1[0], distseq1[0] );
2760                                                                         strcpy( dumseq2[0], distseq2[0] );
2761                                                                         pscore = genL__align11( n_dis_consweight_multi, dumseq1, dumseq2, alloclen, &dum1, &dum2 ); // uwagaki
2762                                                                 }
2763 #if 1
2764                                                                 if( specificityconsideration > 0.0 )
2765                                                                 {
2766 //                                                                      fprintf( stderr, "dist = %f\n", score2dist( pscore, selfscore[i], selfscore[j] ) );
2767                                                                         dist = score2dist( pscore, selfscore[i], selfscore[j] );
2768                                                                         if( dist2offset( dist ) < 0.0 )
2769                                                                         {
2770                                                                                 makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
2771                                                                                 strcpy( mseq1[0], seq[i] );
2772                                                                                 strcpy( mseq2[0], seq[j] );
2773                                                                                 genL__align11( dynamicmtx, mseq1, mseq2, alloclen, &off1, &off2 );
2774                                                                         }
2775                                                                 }
2776 #endif
2777                                                         }
2778                                                         break;
2779                                                 case( 'R' ):
2780                                                         if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2781                                                                 pscore = 0.0;
2782                                                         else
2783                                                                 pscore = (double)lastresx[i][j].score; // all pair
2784                                                         break;
2785                                                 case( 'r' ):
2786                                                         if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) )
2787                                                                 pscore = (double)lastresx[i][j-(njob-nadd)].score;
2788                                                         else
2789                                                                 pscore = 0.0;
2790                                                         break;
2791                                                 case( 'L' ):
2792                                                         if( nadd && njob-nadd <= j && njob-nadd <= i ) // new sequence doushi ha mushi
2793                                                                 pscore = 0.0;
2794                                                         else
2795                                                         {
2796                                                                 if( usenaivescoreinsteadofalignmentscore )
2797                                                                 {
2798                                                                         L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2799                                                                         pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2800                                                                 }
2801                                                                 else
2802                                                                 {
2803 //                                                                      if( store_localhom )
2804                                                                         if( store_localhom && ( targetmap[i] != -1 || targetmap[j] != -1 ) )
2805                                                                         {
2806                                                                                 pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 ); // all pair
2807                                                                                 if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // all pair
2808 #if 1
2809                                                                                 if( specificityconsideration > 0.0 )
2810                                                                                 {
2811                                                                                         dist = score2dist( pscore, selfscore[i], selfscore[j] );
2812                                                                                         if( ( scoreoffset = dist2offset( dist ) ) < 0.0 )
2813                                                                                         {
2814                                                                                                 makedynamicmtx( dynamicmtx, n_dis_consweight_multi, 0.5 * dist ); // upgma ni awaseru.
2815                                                                                                 strcpy( mseq1[0], seq[i] );
2816                                                                                                 strcpy( mseq2[0], seq[j] );
2817                                                                                                 L__align11( dynamicmtx, scoreoffset, mseq1, mseq2, alloclen, &off1, &off2 );
2818                                                                                         }
2819                                                                                 }
2820 #endif
2821                                                                         }
2822                                                                         else
2823                                                                                 pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // all pair
2824                                                                 }
2825                                                         }
2826 //                                                      pscore = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, distseq1, distseq2, alloclen ); // CHUUI!!!!!!
2827                                                         break;
2828                                                 case( 'Y' ):
2829                                                         if( nadd == 0 || ( i < njob-nadd && njob-nadd <= j ) ) // new sequence vs exiting sequence nomi keisan
2830                                                         {
2831                                                                 if( usenaivescoreinsteadofalignmentscore )
2832                                                                 {
2833                                                                         L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2834                                                                         pscore = (double)naivepairscore11( mseq1[0], mseq2[0], 0.0 ); // uwagaki
2835                                                                 }
2836                                                                 else
2837                                                                 {
2838                                                                         if( store_localhom )
2839                                                                         {
2840                                                                                 pscore = L__align11( n_dis_consweight_multi, 0.0, mseq1, mseq2, alloclen, &off1, &off2 );
2841                                                                                 if( thereisx ) pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 ); // uwagaki
2842                                                                         }
2843                                                                         else
2844                                                                                 pscore = L__align11_noalign( n_dis_consweight_multi, distseq1, distseq2 );
2845                                                                 }
2846                                                         }
2847                                                         else
2848                                                                 pscore = 0.0;
2849                                                         break;
2850                                                 case( 'a' ):
2851                                                         pscore = Aalign( mseq1, mseq2, effarr1, effarr2, 1, 1, alloclen );
2852                                                         off1 = off2 = 0;
2853                                                         break;
2854 #if 0
2855                                                 case( 'K' ):
2856                                                         pscore = genG__align11( mseq1, mseq2, alloclen );
2857                                                         off1 = off2 = 0;
2858                                                         break;
2859 #endif
2860                                                 case( 'H' ):
2861                                                         pscore = recallpairfoldalign( mseq1, mseq2, i, j, &off1, &off2, alloclen );
2862                                                         break;
2863                                                 case( 'B' ):
2864                                                 case( 'T' ):
2865                                                         pscore = recalllara( mseq1, mseq2, alloclen );
2866                                                         off1 = off2 = 0;
2867                                                         break;
2868                                                 case( 's' ):
2869                                                         pscore = callmxscarna_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2870                                                         off1 = off2 = 0;
2871                                                         break;
2872                                                 case( 'G' ):
2873                                                         pscore = calldafs_giving_bpp( mseq1, mseq2, bpp[i], bpp[j], alloclen, i, j );
2874                                                         off1 = off2 = 0;
2875                                                         break;
2876                                                 case( 'M' ):
2877                                                         pscore = MSalign11( mseq1, mseq2, alloclen );
2878                                                         break;
2879                                                 default:
2880                                                         ErrorExit( "ERROR IN SOURCE FILE" );
2881                                         }
2882                                 }
2883
2884                                 if( alg == 't' || ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) ) // 't' no jouken ha iranai to omou. if( ( mseq1[0][0] != 0 && mseq2[0][0] != 0  ) )
2885                                 {
2886 #if SCOREOUT
2887                                         fprintf( stderr, "score = %10.2f (%d,%d)\n", pscore, i, j );
2888 #endif
2889 //                                      if( pscore > 0.0 && ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) ) // x-ins-i de seido teika
2890                                         if( ( nadd == 0 || ( alg != 'Y' && alg != 'r' ) || ( i < njob-nadd && njob-nadd <= j ) ) )
2891                                         {
2892                                                 if( !store_localhom )
2893                                                         ;
2894                                                 else if( specifictarget && targetmap[i] == -1 && targetmap[j] == -1)
2895                                                         ;
2896                                                 else if( alg == 'R' )
2897                                                         putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j, lastresx[i]+j, 'h' );
2898                                                 else if( alg == 'r' )
2899                                                         putlocalhom_last( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), lastresx[i]+j-(njob-nadd), 'h' );// ?????
2900                                                 else if( alg == 'H' )
2901                                                         putlocalhom_ext( mseq1[0], mseq2[0], localhomtable[i]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2902                                                 else if( alg == 'Y' )
2903                                                         putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-(njob-nadd), off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2904                                                 else if( !specifictarget && alg != 'S' && alg != 'V' )
2905                                                         putlocalhom2( mseq1[0], mseq2[0], localhomtable[i]+j-i, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2906                                                 else
2907                                                 {
2908                                                         if( targetmap[i] != -1 && targetmap[j] != -1 )
2909                                                         {
2910                                                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2911                                                                 putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' ); // sukoshi muda.
2912                                                         }
2913                                                         else if( targetmap[j] != -1 )
2914                                                                 putlocalhom2( mseq2[0], mseq1[0], localhomtable[targetmap[j]]+i, off2, off1, (int)pscore, strlen( mseq2[0] ), 'h' );
2915                                                         else if( targetmap[i] != -1 )
2916                                                                 putlocalhom2( mseq1[0], mseq2[0], localhomtable[targetmap[i]]+j, off1, off2, (int)pscore, strlen( mseq1[0] ), 'h' );
2917                                                         else
2918                                                         {
2919                                                                 reporterr( "okashii\n" );
2920                                                                 exit( 1 );
2921                                                         }
2922                                                 }
2923                                         }
2924
2925                                         pscore = score2dist( pscore, selfscore[i], selfscore[j] );
2926                                 }
2927                                 else
2928                                 {
2929                                         pscore = 2.0;
2930                                 }
2931
2932                                 if( stdout_align )
2933                                 {
2934                                         if( alg != 't' )
2935                                         {
2936                                                 fprintf( stdout, "sequence %d - sequence %d, pairwise distance = %10.5f\n", i+1, j+1, pscore );
2937                                                 fprintf( stdout, ">%s\n", name[i] );
2938                                                 write1seq( stdout, mseq1[0] );
2939                                                 fprintf( stdout, ">%s\n", name[j] );
2940                                                 write1seq( stdout, mseq2[0] );
2941                                                 fprintf( stdout, "\n" );
2942                                         }
2943                                 }
2944                                 if( stdout_dist ) fprintf( stdout, "%d %d d=%.3f\n", i+1, j+1, pscore );
2945                                 if( store_dist)
2946                                 {
2947                                         if( alg == 'Y' || alg == 'r' ) distancemtx[i][j-(njob-nadd)] = pscore;
2948                                         else distancemtx[i][j-i] = pscore;
2949                                 }
2950                         }
2951                 }
2952                 if( dynamicmtx ) FreeDoubleMtx( dynamicmtx );
2953         }
2954
2955
2956         if( store_dist && ngui == 0 )
2957         {
2958                 hat2p = fopen( hat2file, "w" );
2959                 if( !hat2p ) ErrorExit( "Cannot open hat2." );
2960                 if( alg == 'Y' || alg == 'r' )
2961                         WriteHat2_part_pointer( hat2p, njob, nadd, name, distancemtx );
2962                 else
2963 //                      WriteHat2_pointer( hat2p, njob, name, distancemtx );
2964                         WriteFloatHat2_pointer_halfmtx( hat2p, njob, name, distancemtx ); // jissiha double
2965                 fclose( hat2p );
2966         }
2967
2968         hat3p = fopen( "hat3", "w" );
2969         if( !hat3p ) ErrorExit( "Cannot open hat3." );
2970         if( store_localhom && ngui == 0 )
2971         {
2972
2973                 fprintf( stderr, "\n\n##### writing hat3\n" );
2974                 if( alg == 'Y' || alg == 'r' )
2975                         ilim = njob-nadd;
2976                 else if( specifictarget )
2977                         ilim = ntarget;
2978                 else
2979                         ilim = njob-1;
2980                 for( i=0; i<ilim; i++ )
2981                 {
2982                         if( alg == 'Y' || alg == 'r' )
2983                         {
2984                                 jst = njob-nadd;
2985                                 jj = 0;
2986                         }
2987                         else if( specifictarget )
2988                         {
2989                                 jst = 0;
2990                                 jj = 0;
2991                         }
2992                         else
2993                         {
2994                                 jst = i;
2995                                 jj = 0;
2996                         }
2997                         for( j=jst; j<njob; j++, jj++ )
2998                         {
2999                                 for( tmpptr=localhomtable[i]+jj; tmpptr; tmpptr=tmpptr->next )
3000                                 {
3001 //                                      fprintf( stderr, "j=%d, jj=%d\n", j, jj );
3002                                         if( tmpptr->opt == -1.0 ) continue;
3003 // tmptmptmptmptmp
3004 //                                      if( alg == 'B' || alg == 'T' )
3005 //                                              fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d %p\n", i, j, tmpptr->overlapaa, 1.0, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, (void *)tmpptr->next );
3006 //                                      else
3007                                         if( targetmap[j] == -1 || targetmap[i] < targetmap[j] )
3008                                                 fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", targetmapr[i], j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2 );
3009 //                                              fprintf( hat3p, "%d %d %d %7.5f %d %d %d %d h\n", i, j, tmpptr->overlapaa, tmpptr->opt, tmpptr->start1, tmpptr->end1, tmpptr->start2+1, tmpptr->end2+1 ); // zettai dame!!!!
3010                                 }
3011                         }
3012                 }
3013 //              if( ngui == 0 )
3014 //              {
3015 #if DEBUG
3016                         fprintf( stderr, "calling FreeLocalHomTable\n" );
3017 #endif
3018                         if( alg == 'Y' || alg == 'r' )
3019                                 FreeLocalHomTable_part( localhomtable, (njob-nadd), nadd );
3020                         else if( specifictarget )
3021                                 FreeLocalHomTable_part( localhomtable, ntarget, njob );
3022                         else
3023                                 FreeLocalHomTable_half( localhomtable, njob );
3024 #if DEBUG
3025                         fprintf( stderr, "done. FreeLocalHomTable\n" );
3026 #endif
3027 //              }
3028         }
3029         fclose( hat3p );
3030
3031         if( alg == 's' )
3032         {
3033                 char **ptpt;
3034                 for( i=0; i<njob; i++ )
3035                 {
3036                         ptpt = bpp[i];
3037                         while( 1 )
3038                         {
3039                                 if( *ptpt ) free( *ptpt );
3040                                 else break;
3041                                 ptpt++;
3042                         }
3043                         free( bpp[i] );
3044                 }
3045                 free( bpp );
3046         }
3047         free( selfscore );
3048         free( effarr1 );
3049         free( effarr2 );
3050         if( alg == 'N' )
3051         {
3052                 FreeCharMtx( dumseq1 );
3053                 FreeCharMtx( dumseq2 );
3054         }
3055         free( distseq1 );
3056         free( distseq2 );
3057         if( store_dist && ngui == 0 ) FreeDoubleHalfMtx( distancemtx, njob );
3058
3059         free( targetmap );
3060         free( targetmapr );
3061 }
3062
3063
3064 int pairlocalalign( int ngui, int lgui, char **namegui, char **seqgui, double **distancemtx, LocalHom **localhomtable, int argc, char **argv )
3065 {
3066         int  *nlen, *thereisxineachseq;
3067         char **name, **seq;
3068         char **mseq1, **mseq2;
3069         char **aseq;
3070         char **bseq;
3071         char **dseq;
3072         int i, j, k;
3073         FILE *infp;
3074         char c;
3075         int alloclen;
3076         Lastresx **lastresx;
3077
3078 //      reporterr( "argc=%d, argv[0]=%s\n", argc, argv[0] );
3079
3080         arguments( argc, argv );
3081
3082
3083         if( !ngui )
3084         {
3085                 if( inputfile )
3086                 {
3087                         infp = fopen( inputfile, "r" );
3088                         if( !infp )
3089                         {
3090                                 fprintf( stderr, "Cannot open %s\n", inputfile );
3091                                 exit( 1 );
3092                         }
3093                 }
3094                 else
3095                         infp = stdin;
3096
3097                 getnumlen( infp );
3098                 rewind( infp );
3099
3100                 if( njob < 2 )
3101                 {
3102                         fprintf( stderr, "At least 2 sequences should be input!\n"
3103                                                          "Only %d sequence found.\n", njob );
3104                         exit( 1 );
3105                 }
3106                 if( njob > M )
3107                 {
3108                         fprintf( stderr, "The number of sequences must be < %d\n", M );
3109                         fprintf( stderr, "Please try the splittbfast program for such large data.\n" );
3110                         exit( 1 );
3111                 }
3112         }
3113
3114         if( ( alg == 'r' || alg == 'R' ) && dorp == 'p' )
3115         {
3116                 fprintf( stderr, "Not yet supported\n" );
3117                 exit( 1 );
3118         }
3119
3120         alloclen = nlenmax*2;
3121         if( ngui )
3122         {
3123                 seq = seqgui;
3124                 name = namegui;
3125         }
3126         else
3127         {
3128                 seq = AllocateCharMtx( njob, alloclen+10 );
3129                 name = AllocateCharMtx( njob, B );
3130         }
3131
3132         aseq = AllocateCharMtx( 2, alloclen+10 );
3133         bseq = AllocateCharMtx( njob, alloclen+10 );
3134         dseq = AllocateCharMtx( njob, alloclen+10 );
3135         mseq1 = AllocateCharMtx( njob, 0 );
3136         mseq2 = AllocateCharMtx( njob, 0 );
3137         nlen = AllocateIntVec( njob );
3138         thereisxineachseq = AllocateIntVec( njob );
3139
3140
3141         if( alg == 'R' )
3142         {
3143                 lastresx = calloc( njob+1, sizeof( Lastresx * ) );
3144                 for( i=0; i<njob; i++ )
3145                 {
3146                         lastresx[i] = calloc( njob+1, sizeof( Lastresx ) ); // muda
3147                         for( j=0; j<njob; j++ )
3148                         {
3149                                 lastresx[i][j].score = 0;
3150                                 lastresx[i][j].naln = 0;
3151                                 lastresx[i][j].aln = NULL;
3152                         }
3153                         lastresx[i][njob].naln = -1;
3154                 }
3155                 lastresx[njob] = NULL;
3156         }
3157         else if( alg == 'r' )
3158         {
3159 //              fprintf( stderr, "Allocating lastresx (%d), njob=%d, nadd=%d\n", njob-nadd+1, njob, nadd );
3160                 lastresx = calloc( njob-nadd+1, sizeof( Lastresx * ) );
3161                 for( i=0; i<njob-nadd; i++ )
3162                 {
3163 //                      fprintf( stderr, "Allocating lastresx[%d]\n", i );
3164                         lastresx[i] = calloc( nadd+1, sizeof( Lastresx ) );
3165                         for( j=0; j<nadd; j++ )
3166                         {
3167 //                              fprintf( stderr, "Initializing lastresx[%d][%d]\n", i, j );
3168                                 lastresx[i][j].score = 0;
3169                                 lastresx[i][j].naln = 0;
3170                                 lastresx[i][j].aln = NULL;
3171                         }
3172                         lastresx[i][nadd].naln = -1;
3173                 }
3174                 lastresx[njob-nadd] = NULL;
3175         }
3176         else
3177                 lastresx = NULL;
3178
3179 #if 0
3180         Read( name, nlen, seq );
3181 #else
3182         if( !ngui )
3183         {
3184                 readData_pointer( infp, name, nlen, seq );
3185                 fclose( infp );
3186         }
3187 #endif
3188
3189         constants( njob, seq );
3190
3191 #if 0
3192         fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
3193 #endif
3194
3195         initSignalSM();
3196
3197         initFiles();
3198
3199 //      WriteOptions( trap_g );
3200
3201         c = seqcheck( seq );
3202         if( c )
3203         {
3204                 fprintf( stderr, "Illegal character %c\n", c );
3205                 exit( 1 );
3206         }
3207
3208 //      writePre( njob, name, nlen, seq, 0 );
3209
3210
3211
3212
3213         for( i=0; i<njob; i++ )
3214         {
3215                 gappick0( bseq[i], seq[i] );
3216                 thereisxineachseq[i] = removex( dseq[i], bseq[i] );
3217         }
3218
3219         pairalign( name, nlen, bseq, aseq, dseq, thereisxineachseq, mseq1, mseq2, alloclen, lastresx, distancemtx, localhomtable, ngui );
3220
3221         fprintf( trap_g, "done.\n" );
3222 #if DEBUG
3223         fprintf( stderr, "closing trap_g\n" );
3224 #endif
3225         fclose( trap_g );
3226         fclose( prep_g );
3227
3228 //      writePre( njob, name, nlen, aseq, !contin );
3229 #if 0
3230         writeData( stdout, njob, name, nlen, aseq );
3231 #endif
3232 #if IODEBUG
3233         fprintf( stderr, "OSHIMAI\n" );
3234 #endif
3235         SHOWVERSION;
3236
3237         if( stdout_dist && nthread > 1 )
3238         {
3239                 fprintf( stderr, "\nThe order of distances is not identical to that in the input file, because of the parallel calculation.  Reorder them by yourself, using sort -n -k 2 | sort -n -k 1 -s\n" );
3240         }
3241         if( stdout_align && nthread > 1 )
3242         {
3243                 fprintf( stderr, "\nThe order of pairwise alignments is not identical to that in the input file, because of the parallel calculation.  Reorder them by yourself.\n" );
3244         }
3245
3246 #if 1
3247         if( lastresx )
3248         {
3249                 for( i=0; lastresx[i]; i++ )
3250                 {
3251                         for( j=0; lastresx[i][j].naln!=-1; j++ )
3252                         {
3253                                 for( k=0; k<lastresx[i][j].naln; k++ )
3254                                 {
3255                                         free( lastresx[i][j].aln[k].reg1 );
3256                                         free( lastresx[i][j].aln[k].reg2 );
3257                                 }
3258                                 free( lastresx[i][j].aln );
3259                         }
3260                         free( lastresx[i] );
3261                 }
3262                 free( lastresx );
3263         }
3264 #endif
3265         if( ngui == 0 )
3266         {
3267                 FreeCharMtx( seq );
3268                 FreeCharMtx( name );
3269         }
3270         FreeCharMtx( aseq );
3271         FreeCharMtx( bseq );
3272         FreeCharMtx( dseq );
3273         free( mseq1 );
3274         free( mseq2 );
3275         free( nlen );
3276         free( thereisxineachseq );
3277         freeconstants();
3278
3279         if( !ngui )
3280         {
3281                 FreeCommonIP();
3282         }
3283         Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
3284         G__align11( NULL, NULL, NULL, 0, 0, 0 ); // 20130603
3285         G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
3286         L__align11( NULL, 0.0, NULL, NULL, 0, NULL, NULL );
3287         L__align11_noalign( NULL, NULL, NULL );
3288         genL__align11( NULL, NULL, NULL, 0, NULL, NULL );
3289
3290 #if SHISHAGONYU
3291         if( ngui )
3292         {
3293                 char buf[100];
3294                 for( i=0; i<njob-1; i++ ) for( j=i+1; j<njob; j++ )
3295                 {
3296                         sprintf( buf, "%5.3f", distancemtx[i][j-i] );
3297                         distancemtx[i][j-i] = 0.0;
3298                         sscanf( buf, "%lf", distancemtx[i]+j-i );
3299 //                      distancemtx[i][j-i] = 0.001 * (int)(distancemtx[i][j-i] * 1000 + 0.5);
3300                 }
3301
3302         }
3303 #endif
3304
3305
3306         return( 0 );
3307 }
3308