14 typedef struct _jobtable
20 #ifdef enablemultithread
21 typedef struct _distancematrixthread_arg
29 pthread_mutex_t *mutex;
30 } distancematrixthread_arg_t;
32 typedef struct _treebasethread_arg
44 LocalHom **localhomtable;
48 pthread_mutex_t *mutex;
49 pthread_cond_t *treecond;
50 } treebasethread_arg_t;
53 void arguments( int argc, char *argv[] )
99 ppenalty = NOTSPECIFIED;
100 ppenalty_ex = NOTSPECIFIED;
101 poffset = NOTSPECIFIED;
102 kimuraR = NOTSPECIFIED;
105 fftWinSize = NOTSPECIFIED;
106 fftThreshold = NOTSPECIFIED;
107 RNAppenalty = NOTSPECIFIED;
108 RNAppenalty_ex = NOTSPECIFIED;
109 RNApthr = NOTSPECIFIED;
111 consweight_multi = 1.0;
112 consweight_rna = 0.0;
114 while( --argc > 0 && (*++argv)[0] == '-' )
116 while ( ( c = *++argv[0] ) )
122 fprintf( stderr, "inputfile = %s\n", inputfile );
126 nadd = atoi( *++argv );
127 fprintf( stderr, "nadd = %d\n", nadd );
131 RNApthr = (int)( atof( *++argv ) * 1000 - 0.5 );
135 RNAppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
139 ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
140 // fprintf( stderr, "ppenalty = %d\n", ppenalty );
144 ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
145 fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
149 poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
150 // fprintf( stderr, "poffset = %d\n", poffset );
154 kimuraR = atoi( *++argv );
155 fprintf( stderr, "kappa = %d\n", kimuraR );
159 nblosum = atoi( *++argv );
161 fprintf( stderr, "blosum %d / kimura 200\n", nblosum );
165 pamN = atoi( *++argv );
168 fprintf( stderr, "jtt/kimura %d\n", pamN );
172 pamN = atoi( *++argv );
175 fprintf( stderr, "tm %d\n", pamN );
179 fastathreshold = atof( *++argv );
184 consweight_rna = atof( *++argv );
189 consweight_multi = atof( *++argv );
193 nthread = atoi( *++argv );
194 fprintf( stderr, "nthread = %d\n", nthread );
309 /* Modified 01/08/27, default: user tree */
313 /* modification end. */
315 fftThreshold = atoi( *++argv );
319 fftWinSize = atoi( *++argv );
326 fprintf( stderr, "illegal option %c\n", c );
336 cut = atof( (*argv) );
341 fprintf( stderr, "options: Check source file !\n" );
344 if( tbitr == 1 && outgap == 0 )
346 fprintf( stderr, "conflicting options : o, m or u\n" );
349 if( alg == 'C' && outgap == 0 )
351 fprintf( stderr, "conflicting options : C, o\n" );
357 static void *distancematrixthread2( void *arg )
359 distancematrixthread_arg_t *targ = (distancematrixthread_arg_t *)arg;
360 int njob = targ->njob;
361 int thread_no = targ->thread_no;
362 float *selfscore = targ->selfscore;
363 float **iscore = targ->iscore;
364 char **seq = targ->seq;
365 Jobtable *jobpospt = targ->jobpospt;
367 float ssi, ssj, bunbo;
372 pthread_mutex_lock( targ->mutex );
377 pthread_mutex_unlock( targ->mutex );
381 pthread_mutex_unlock( targ->mutex );
384 if( i % 10 == 0 ) fprintf( stderr, "\r% 5d / %d (thread %4d)", i, njob, thread_no );
385 for( j=i+1; j<njob; j++)
388 bunbo = MIN( ssi, ssj );
390 iscore[i][j-i] = 1.0;
392 iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
398 #ifdef enablemultithread
399 static void *distancematrixthread( void *arg )
401 distancematrixthread_arg_t *targ = (distancematrixthread_arg_t *)arg;
402 int njob = targ->njob;
403 int thread_no = targ->thread_no;
404 float *selfscore = targ->selfscore;
405 float **iscore = targ->iscore;
406 char **seq = targ->seq;
407 Jobtable *jobpospt = targ->jobpospt;
409 float ssi, ssj, bunbo;
414 pthread_mutex_lock( targ->mutex );
424 pthread_mutex_unlock( targ->mutex );
430 pthread_mutex_unlock( targ->mutex );
433 if( j==i+1 && i % 10 == 0 ) fprintf( stderr, "\r% 5d / %d (thread %4d)", i, njob, thread_no );
436 bunbo = MIN( ssi, ssj );
438 iscore[i][j-i] = 1.0;
440 iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
444 static void *treebasethread( void *arg )
446 treebasethread_arg_t *targ = (treebasethread_arg_t *)arg;
447 int *nrunpt = targ->nrunpt;
448 int thread_no = targ->thread_no;
449 int njob = targ->njob;
450 int *nlen = targ->nlen;
451 int *jobpospt = targ->jobpospt;
452 int ***topol = targ->topol;
453 Treedep *dep = targ->dep;
454 char **aseq = targ->aseq;
455 double *effarr = targ->effarr;
456 int *alloclen = targ->alloclenpt;
457 LocalHom **localhomtable = targ->localhomtable;
458 RNApair ***singlerna = targ->singlerna;
459 double *effarr_kozo = targ->effarr_kozo;
460 int *fftlog = targ->fftlog;
462 char **mseq1, **mseq2;
468 char *indication1, *indication2;
469 double *effarr1 = NULL;
470 double *effarr2 = NULL;
471 double *effarr1_kozo = NULL;
472 double *effarr2_kozo = NULL;
473 LocalHom ***localhomshrink = NULL;
477 RNApair ***grouprna1, ***grouprna2;
480 mseq1 = AllocateCharMtx( njob, 0 );
481 mseq2 = AllocateCharMtx( njob, 0 );
482 localcopy = calloc( njob, sizeof( char * ) );
484 if( rnakozo && rnaprediction == 'm' )
486 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
487 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
491 grouprna1 = grouprna2 = NULL;
494 effarr1 = AllocateDoubleVec( njob );
495 effarr2 = AllocateDoubleVec( njob );
496 indication1 = AllocateCharVec( 150 );
497 indication2 = AllocateCharVec( 150 );
502 localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
503 for( i=0; i<njob; i++)
504 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
507 effarr1_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
508 effarr2_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
509 for( i=0; i<njob; i++ ) effarr1_kozo[i] = 0.0;
510 for( i=0; i<njob; i++ ) effarr2_kozo[i] = 0.0;
517 for( i=0; i<njob; i++ )
518 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
521 #if 0 //-> main thread
523 calcimportance( njob, effarr, aseq, localhomtable );
527 // writePre( njob, name, nlen, aseq, 0 );
530 // for( l=0; l<njob-1; l++ )
534 pthread_mutex_lock( targ->mutex );
538 pthread_mutex_unlock( targ->mutex );
539 if( commonIP ) FreeIntMtx( commonIP );
541 Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
542 A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
548 free( effarr1_kozo );
549 free( effarr2_kozo );
554 for( i=0; i<njob; i++)
555 free( localhomshrink[i] );
556 free( localhomshrink );
562 if( dep[l].child0 != -1 )
564 while( dep[dep[l].child0].done == 0 )
565 pthread_cond_wait( targ->treecond, targ->mutex );
567 if( dep[l].child1 != -1 )
569 while( dep[dep[l].child1].done == 0 )
570 pthread_cond_wait( targ->treecond, targ->mutex );
572 // while( *nrunpt >= nthread )
573 // pthread_cond_wait( targ->treecond, targ->mutex );
576 // pthread_mutex_unlock( targ->mutex );
582 // pthread_mutex_lock( targ->mutex );
584 len1 = strlen( aseq[m1] );
585 len2 = strlen( aseq[m2] );
586 if( *alloclen <= len1 + len2 )
588 fprintf( stderr, "\nReallocating (by thread %d) ..", thread_no );
589 *alloclen = ( len1 + len2 ) + 1000;
590 ReallocateCharMtx( aseq, njob, *alloclen + 10 );
591 fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
594 for( i=0; (j=topol[l][0][i])!=-1; i++ )
596 localcopy[j] = calloc( *alloclen, sizeof( char ) );
597 strcpy( localcopy[j], aseq[j] );
599 for( i=0; (j=topol[l][1][i])!=-1; i++ )
601 localcopy[j] = calloc( *alloclen, sizeof( char ) );
602 strcpy( localcopy[j], aseq[j] );
605 pthread_mutex_unlock( targ->mutex );
609 clus1 = fastconjuction_noname_kozo( topol[l][0], localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
610 clus2 = fastconjuction_noname_kozo( topol[l][1], localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
614 clus1 = fastconjuction_noname( topol[l][0], localcopy, mseq1, effarr1, effarr, indication1 );
615 clus2 = fastconjuction_noname( topol[l][1], localcopy, mseq2, effarr2, effarr, indication2 );
621 fprintf( stderr, "\rSTEP % 5d /%d (thread %4d) ", l+1, njob-1, thread_no );
623 fprintf( stderr, "STEP %d /%d (thread %d) \n", l+1, njob-1, thread_no );
624 fprintf( stderr, "group1 = %.66s", indication1 );
625 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
626 fprintf( stderr, ", child1 = %d\n", dep[l].child0 );
627 fprintf( stderr, "group2 = %.66s", indication2 );
628 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
629 fprintf( stderr, ", child2 = %d\n", dep[l].child1 );
631 fprintf( stderr, "Group1's lengths = " );
632 for( i=0; i<clus1; i++ ) fprintf( stderr, "%d ", strlen( mseq1[i] ) );
633 fprintf( stderr, "\n" );
634 fprintf( stderr, "Group2's lengths = " );
635 for( i=0; i<clus2; i++ ) fprintf( stderr, "%d ", strlen( mseq2[i] ) );
636 fprintf( stderr, "\n" );
642 // for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
646 fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
647 // msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
648 // fprintf( stderr, "localhomshrink =\n" );
649 // outlocalhompt( localhomshrink, clus1, clus2 );
650 // weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
651 // fprintf( stderr, "after weight =\n" );
652 // outlocalhompt( localhomshrink, clus1, clus2 );
654 if( rnakozo && rnaprediction == 'm' )
656 makegrouprna( grouprna1, singlerna, topol[l][0] );
657 makegrouprna( grouprna2, singlerna, topol[l][1] );
662 fprintf( stderr, "before align all\n" );
663 display( localcopy, njob );
664 fprintf( stderr, "\n" );
665 fprintf( stderr, "before align 1 %s \n", indication1 );
666 display( mseq1, clus1 );
667 fprintf( stderr, "\n" );
668 fprintf( stderr, "before align 2 %s \n", indication2 );
669 display( mseq2, clus2 );
670 fprintf( stderr, "\n" );
673 if( !nevermemsave && ( constraint != 2 && alg != 'M' && ( len1 > 30000 || len2 > 30000 ) ) )
675 fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
677 if( commonIP ) FreeIntMtx( commonIP );
683 // if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
684 if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
686 // ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
687 // fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
688 // fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
689 if( constraint == 2 )
693 fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
696 fprintf( stderr, "c" );
699 imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
700 if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
701 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
703 else if( alg == 'H' )
705 imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
706 pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
708 else if( alg == 'Q' )
710 imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
711 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
712 pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
714 else if( alg == 'R' )
716 imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
717 pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
720 else if( force_fft || ( use_fft && ffttry ) )
722 fprintf( stderr, "f" );
725 fprintf( stderr, "m" );
726 pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
729 pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
733 fprintf( stderr, "d" );
738 pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
741 fprintf( stderr, "m" );
742 pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
745 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
748 pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
751 pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
754 pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
757 ErrorExit( "ERROR IN SOURCE FILE" );
761 nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
764 fprintf( stderr, "score = %10.2f\n", pscore );
768 fprintf( stderr, "after align 1 %s \n", indication1 );
769 display( mseq1, clus1 );
770 fprintf( stderr, "\n" );
771 fprintf( stderr, "after align 2 %s \n", indication2 );
772 display( mseq2, clus2 );
773 fprintf( stderr, "\n" );
776 // writePre( njob, name, nlen, localcopy, 0 );
778 if( disp ) display( localcopy, njob );
780 pthread_mutex_lock( targ->mutex );
783 pthread_cond_broadcast( targ->treecond );
785 // pthread_mutex_unlock( targ->mutex );
786 // pthread_mutex_lock( targ->mutex );
788 for( i=0; (j=topol[l][0][i])!=-1; i++ )
789 strcpy( aseq[j], localcopy[j] );
790 for( i=0; (j=topol[l][1][i])!=-1; i++ )
791 strcpy( aseq[j], localcopy[j] );
792 pthread_mutex_unlock( targ->mutex );
794 for( i=0; (j=topol[l][0][i])!=-1; i++ )
795 free( localcopy[j] );
796 for( i=0; (j=topol[l][1][i])!=-1; i++ )
797 free( localcopy[j] );
806 void treebase( int *nlen, char **aseq, int nadd, char *mergeoralign, char **mseq1, char **mseq2, int ***topol, double *effarr, int *alloclen, LocalHom **localhomtable, RNApair ***singlerna, double *effarr_kozo )
809 int len1nocommongap, len2nocommongap;
812 float pscore, tscore;
813 static char *indication1, *indication2;
814 static double *effarr1 = NULL;
815 static double *effarr2 = NULL;
816 static double *effarr1_kozo = NULL;
817 static double *effarr2_kozo = NULL;
818 static LocalHom ***localhomshrink = NULL;
823 static int *alreadyaligned;
826 RNApair ***grouprna1, ***grouprna2;
828 if( rnakozo && rnaprediction == 'm' )
830 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
831 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
835 grouprna1 = grouprna2 = NULL;
838 if( effarr1 == NULL )
840 fftlog = AllocateIntVec( njob );
841 effarr1 = AllocateDoubleVec( njob );
842 effarr2 = AllocateDoubleVec( njob );
843 indication1 = AllocateCharVec( 150 );
844 indication2 = AllocateCharVec( 150 );
845 gaplen = AllocateIntVec( *alloclen+10 );
846 gapmap = AllocateIntVec( *alloclen+10 );
847 alreadyaligned = AllocateIntVec( njob );
852 localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
853 for( i=0; i<njob; i++)
854 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
857 effarr1_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
858 effarr2_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
859 for( i=0; i<njob; i++ ) effarr1_kozo[i] = 0.0;
860 for( i=0; i<njob; i++ ) effarr2_kozo[i] = 0.0;
863 for( i=0; i<njob-nadd; i++ ) alreadyaligned[i] = 1;
864 for( i=njob-nadd; i<njob; i++ ) alreadyaligned[i] = 0;
866 for( l=0; l<njob; l++ ) fftlog[l] = 1;
869 fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
873 for( i=0; i<njob; i++ )
874 fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
879 calcimportance( njob, effarr, aseq, localhomtable );
882 // writePre( njob, name, nlen, aseq, 0 );
885 for( l=0; l<njob-1; l++ )
887 if( mergeoralign[l] == 'n' )
889 // fprintf( stderr, "SKIP!\n" );
898 len1 = strlen( aseq[m1] );
899 len2 = strlen( aseq[m2] );
900 if( *alloclen < len1 + len2 )
902 fprintf( stderr, "\nReallocating.." );
903 *alloclen = ( len1 + len2 ) + 1000;
904 ReallocateCharMtx( aseq, njob, *alloclen + 10 );
905 gaplen = realloc( gaplen, ( *alloclen + 10 ) * sizeof( int ) );
908 fprintf( stderr, "Cannot realloc gaplen\n" );
911 gapmap = realloc( gapmap, ( *alloclen + 10 ) * sizeof( int ) );
914 fprintf( stderr, "Cannot realloc gapmap\n" );
917 fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
922 clus1 = fastconjuction_noname_kozo( topol[l][0], aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
923 clus2 = fastconjuction_noname_kozo( topol[l][1], aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
927 clus1 = fastconjuction_noname( topol[l][0], aseq, mseq1, effarr1, effarr, indication1 );
928 clus2 = fastconjuction_noname( topol[l][1], aseq, mseq2, effarr2, effarr, indication2 );
931 if( mergeoralign[l] == '1' || mergeoralign[l] == '2' )
939 len1nocommongap = len1;
940 len2nocommongap = len2;
941 if( mergeoralign[l] == '1' ) // nai
943 findcommongaps( clus2, mseq2, gapmap );
944 commongappick( clus2, mseq2 );
945 len2nocommongap = strlen( mseq2[0] );
947 else if( mergeoralign[l] == '2' )
949 findcommongaps( clus1, mseq1, gapmap );
950 commongappick( clus1, mseq1 );
951 len1nocommongap = strlen( mseq1[0] );
955 fprintf( trap_g, "\nSTEP-%d\n", l );
956 fprintf( trap_g, "group1 = %s\n", indication1 );
957 fprintf( trap_g, "group2 = %s\n", indication2 );
960 fprintf( stderr, "\rSTEP % 5d /%d ", l+1, njob-1 );
963 fprintf( stdout, "STEP %d /%d\n", l+1, njob-1 );
964 fprintf( stderr, "STEP %d /%d\n", l+1, njob-1 );
965 fprintf( stderr, "group1 = %.66s", indication1 );
966 if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
967 fprintf( stderr, "\n" );
968 fprintf( stderr, "group2 = %.66s", indication2 );
969 if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
970 fprintf( stderr, "\n" );
975 // for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
979 fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
980 // msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
981 // fprintf( stderr, "localhomshrink =\n" );
982 // outlocalhompt( localhomshrink, clus1, clus2 );
983 // weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
984 // fprintf( stderr, "after weight =\n" );
985 // outlocalhompt( localhomshrink, clus1, clus2 );
987 if( rnakozo && rnaprediction == 'm' )
989 makegrouprna( grouprna1, singlerna, topol[l][0] );
990 makegrouprna( grouprna2, singlerna, topol[l][1] );
995 fprintf( stderr, "before align all\n" );
996 display( aseq, njob );
997 fprintf( stderr, "\n" );
998 fprintf( stderr, "before align 1 %s \n", indication1 );
999 display( mseq1, clus1 );
1000 fprintf( stderr, "\n" );
1001 fprintf( stderr, "before align 2 %s \n", indication2 );
1002 display( mseq2, clus2 );
1003 fprintf( stderr, "\n" );
1006 if( !nevermemsave && ( constraint != 2 && alg != 'M' && ( len1 > 30000 || len2 > 30000 ) ) )
1008 fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
1010 if( commonIP ) FreeIntMtx( commonIP );
1016 // if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
1017 if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
1019 // ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
1020 // fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
1021 // fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
1022 if( constraint == 2 )
1026 fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
1029 fprintf( stderr, "c" );
1032 imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1033 if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1034 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1036 else if( alg == 'H' )
1038 imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1039 pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1041 else if( alg == 'Q' )
1043 imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1044 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1045 pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1047 else if( alg == 'R' )
1049 imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1050 pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1053 else if( force_fft || ( use_fft && ffttry ) )
1055 fprintf( stderr, "f" );
1058 fprintf( stderr, "m" );
1059 pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
1062 pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
1066 fprintf( stderr, "d" );
1071 pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
1074 fprintf( stderr, "m" );
1075 pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1078 pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1081 pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1084 pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1087 pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1090 ErrorExit( "ERROR IN SOURCE FILE" );
1094 nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
1097 fprintf( stderr, "score = %10.2f\n", pscore );
1101 fprintf( stderr, "after align 1 %s \n", indication1 );
1102 display( mseq1, clus1 );
1103 fprintf( stderr, "\n" );
1104 fprintf( stderr, "after align 2 %s \n", indication2 );
1105 display( mseq2, clus2 );
1106 fprintf( stderr, "\n" );
1109 // writePre( njob, name, nlen, aseq, 0 );
1111 if( disp ) display( aseq, njob );
1113 if( mergeoralign[l] == '1' ) // jissainiha nai. atarashii hairetsu ha saigo dakara.
1115 adjustgapmap( strlen( mseq2[0] )-len2nocommongap+len2, gapmap, mseq2[0] );
1116 restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen );
1117 findnewgaps( clus2, mseq2, gaplen );
1118 insertnewgaps( njob, alreadyaligned, aseq, topol[l][1], topol[l][0], gaplen, gapmap, *alloclen, alg );
1119 for( i=0; i<njob; i++ ) eq2dash( aseq[i] );
1120 for( i=0; (m=topol[l][0][i])>-1; i++ ) alreadyaligned[m] = 1;
1122 if( mergeoralign[l] == '2' )
1124 // fprintf( stderr, ">mseq1[0] = \n%s\n", mseq1[0] );
1125 // fprintf( stderr, ">mseq2[0] = \n%s\n", mseq2[0] );
1126 adjustgapmap( strlen( mseq1[0] )-len1nocommongap+len1, gapmap, mseq1[0] );
1127 restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen );
1128 findnewgaps( clus1, mseq1, gaplen );
1129 insertnewgaps( njob, alreadyaligned, aseq, topol[l][0], topol[l][1], gaplen, gapmap, *alloclen, alg );
1130 for( i=0; i<njob; i++ ) eq2dash( aseq[i] );
1131 for( i=0; (m=topol[l][1][i])>-1; i++ ) alreadyaligned[m] = 1;
1134 free( topol[l][0] );
1135 free( topol[l][1] );
1139 fprintf( stderr, "totalscore = %10.2f\n\n", tscore );
1143 static void WriteOptions( FILE *fp )
1146 if( dorp == 'd' ) fprintf( fp, "DNA\n" );
1149 if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
1150 else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
1151 else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
1153 fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1154 if( use_fft ) fprintf( fp, "FFT on\n" );
1156 fprintf( fp, "tree-base method\n" );
1157 if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
1158 else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
1159 if( tbitr || tbweight )
1161 fprintf( fp, "iterate at each step\n" );
1162 if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
1163 if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
1164 if( tbweight ) fprintf( fp, " weighted\n" );
1165 fprintf( fp, "\n" );
1168 fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1171 fprintf( fp, "Algorithm A\n" );
1172 else if( alg == 'A' )
1173 fprintf( fp, "Algorithm A+\n" );
1174 else if( alg == 'C' )
1175 fprintf( fp, "Apgorithm A+/C\n" );
1177 fprintf( fp, "Unknown algorithm\n" );
1179 if( treemethod == 'X' )
1180 fprintf( fp, "Tree = UPGMA (mix).\n" );
1181 else if( treemethod == 'E' )
1182 fprintf( fp, "Tree = UPGMA (average).\n" );
1183 else if( treemethod == 'q' )
1184 fprintf( fp, "Tree = Minimum linkage.\n" );
1186 fprintf( fp, "Unknown tree.\n" );
1190 fprintf( fp, "FFT on\n" );
1192 fprintf( fp, "Basis : 4 nucleotides\n" );
1196 fprintf( fp, "Basis : Polarity and Volume\n" );
1198 fprintf( fp, "Basis : 20 amino acids\n" );
1200 fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
1201 fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
1204 fprintf( fp, "FFT off\n" );
1209 int main( int argc, char *argv[] )
1212 static float *selfscore;
1214 static char **name, **seq;
1215 static char **mseq1, **mseq2;
1217 static float **iscore, **iscore_kozo;
1218 static double *eff, *eff_kozo, *eff_kozo_mapped = NULL;
1219 int i, j, ien, ik, jk;
1220 static int ***topol, ***topol_kozo;
1222 static Treedep *dep;
1223 static float **len, **len_kozo;
1228 double unweightedspscore;
1229 int alignmentlength;
1235 LocalHom **localhomtable = NULL;
1236 RNApair ***singlerna;
1237 float ssi, ssj, bunbo;
1238 static char *kozoarivec;
1241 arguments( argc, argv );
1242 #ifndef enablemultithread
1248 infp = fopen( inputfile, "r" );
1251 fprintf( stderr, "Cannot open %s\n", inputfile );
1266 fprintf( stderr, "At least 2 sequences should be input!\n"
1267 "Only %d sequence found.\n", njob );
1271 seq = AllocateCharMtx( njob, nlenmax+1 );
1272 mseq1 = AllocateCharMtx( njob, 0 );
1273 mseq2 = AllocateCharMtx( njob, 0 );
1275 name = AllocateCharMtx( njob, B+1 );
1276 nlen = AllocateIntVec( njob );
1277 selfscore = AllocateFloatVec( njob );
1279 topol = AllocateIntCub( njob, 2, 0 );
1280 len = AllocateFloatMtx( njob, 2 );
1281 iscore = AllocateFloatHalfMtx( njob );
1282 eff = AllocateDoubleVec( njob );
1283 kozoarivec = AllocateCharVec( njob );
1285 mergeoralign = AllocateCharVec( njob );
1287 dep = (Treedep *)calloc( njob, sizeof( Treedep ) );
1288 if( nadd ) addmem = AllocateIntVec( nadd+1 );
1292 localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
1293 for( i=0; i<njob; i++)
1295 localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1296 for( j=0; j<njob; j++)
1298 localhomtable[i][j].start1 = -1;
1299 localhomtable[i][j].end1 = -1;
1300 localhomtable[i][j].start2 = -1;
1301 localhomtable[i][j].end2 = -1;
1302 localhomtable[i][j].overlapaa = -1.0;
1303 localhomtable[i][j].opt = -1.0;
1304 localhomtable[i][j].importance = -1.0;
1305 localhomtable[i][j].next = NULL;
1306 localhomtable[i][j].korh = 'h';
1310 fprintf( stderr, "Loading 'hat3' ... " );
1311 prep = fopen( "hat3", "r" );
1312 if( prep == NULL ) ErrorExit( "Make hat3." );
1313 readlocalhomtable( prep, njob, localhomtable, kozoarivec );
1315 fprintf( stderr, "\ndone.\n" );
1319 for( i=0; i<njob; i++ )
1321 // fprintf( stderr, "kozoarivec[%d] = %d\n", i, kozoarivec[i] );
1322 if( kozoarivec[i] ) nkozo++;
1326 topol_kozo = AllocateIntCub( nkozo, 2, 0 );
1327 len_kozo = AllocateFloatMtx( nkozo, 2 );
1328 iscore_kozo = AllocateFloatHalfMtx( nkozo );
1329 eff_kozo = AllocateDoubleVec( nkozo );
1330 eff_kozo_mapped = AllocateDoubleVec( njob );
1334 // outlocalhom( localhomtable, njob );
1337 fprintf( stderr, "Extending localhom ... " );
1338 extendlocalhom2( njob, localhomtable );
1339 fprintf( stderr, "done.\n" );
1344 readData( infp, name, nlen, seq );
1346 readData_pointer( infp, name, nlen, seq );
1350 constants( njob, seq );
1353 fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
1360 WriteOptions( trap_g );
1362 c = seqcheck( seq );
1365 fprintf( stderr, "Illegal character %c\n", c );
1369 // writePre( njob, name, nlen, seq, 0 );
1376 fprintf( stderr, "Both structure and user tree have been given. Not yet supported!\n" );
1380 fprintf( stderr, "Loading a tree ... " );
1381 loadtree( njob, topol, len, name, nlen, dep );
1382 fprintf( stderr, "\ndone.\n\n" );
1388 for( i=1; i<njob; i++ )
1390 if( nlen[i] != nlen[0] )
1392 fprintf( stderr, "Input pre-aligned seqences or make hat2.\n" );
1397 fprintf( stderr, "Making a distance matrix .. \n" );
1400 for( i=0; i<njob; i++ )
1402 selfscore[i] = naivepairscore11( seq[i], seq[i], penalty );
1404 #ifdef enablemultithread
1407 distancematrixthread_arg_t *targ;
1410 pthread_mutex_t mutex;
1415 targ = calloc( nthread, sizeof( distancematrixthread_arg_t ) );
1416 handle = calloc( nthread, sizeof( pthread_t ) );
1417 pthread_mutex_init( &mutex, NULL );
1419 for( i=0; i<nthread; i++ )
1421 targ[i].thread_no = i;
1422 targ[i].njob = njob;
1423 targ[i].selfscore = selfscore;
1424 targ[i].iscore = iscore;
1426 targ[i].jobpospt = &jobpos;
1427 targ[i].mutex = &mutex;
1429 pthread_create( handle+i, NULL, distancematrixthread, (void *)(targ+i) );
1432 for( i=0; i<nthread; i++ )
1434 pthread_join( handle[i], NULL );
1436 pthread_mutex_destroy( &mutex );
1443 for( i=0; i<ien; i++ )
1447 fprintf( stderr, "\r% 5d / %d", i, ien );
1451 for( j=i+1; j<njob; j++ )
1454 bunbo = MIN( ssi, ssj );
1456 iscore[i][j-i] = 1.0;
1458 // iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / MIN( selfscore[i], selfscore[j] );
1459 iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
1462 fprintf( stderr, "### ssj = %f\n", ssj );
1463 fprintf( stderr, "### selfscore[i] = %f\n", selfscore[i] );
1464 fprintf( stderr, "### selfscore[j] = %f\n", selfscore[j] );
1465 fprintf( stderr, "### rawscore = %f\n", naivepairscore11( seq[i], seq[j], penalty ) );
1470 fprintf( stderr, "\ndone.\n\n" );
1475 fprintf( stderr, "Loading 'hat2' ... " );
1476 prep = fopen( "hat2", "r" );
1477 if( prep == NULL ) ErrorExit( "Make hat2." );
1478 readhat2_floathalf_pointer( prep, njob, name, iscore );
1480 fprintf( stderr, "done.\n" );
1485 hat2p = fopen( "hat2", "w" );
1486 WriteFloatHat2_pointer_halfmtx( hat2p, njob, name, iscore );
1494 for( i=0; i<ien; i++ )
1497 for( j=i+1; j<njob; j++ )
1499 if( kozoarivec[i] && kozoarivec[j] )
1501 iscore_kozo[ik][jk-ik] = iscore[i][j-i];
1503 if( kozoarivec[j] ) jk++;
1505 if( kozoarivec[i] ) ik++;
1509 fprintf( stderr, "Constructing a UPGMA tree ... " );
1513 fprintf( stderr, "Loading a topology ... " );
1514 loadtop( njob, iscore, topol, len );
1515 fprintf( stderr, "\ndone.\n\n" );
1519 fixed_musclesupg_float_realloc_nobk_halfmtx_treeout( njob, iscore, topol, len, name, nlen, dep );
1523 fixed_musclesupg_float_realloc_nobk_halfmtx( njob, iscore, topol, len, dep );
1526 // ErrorExit( "Incorrect tree\n" );
1530 // for( i=0; i<nkozo-1; i++ )
1531 // for( j=i+1; j<nkozo; j++ )
1532 // fprintf( stderr, "iscore_kozo[%d][%d] =~ %f\n", i, j, iscore_kozo[i][j-i] );
1533 fixed_musclesupg_float_realloc_nobk_halfmtx( nkozo, iscore_kozo, topol_kozo, len_kozo, NULL );
1535 fprintf( stderr, "\ndone.\n\n" );
1540 orderfp = fopen( "order", "w" );
1543 fprintf( stderr, "Cannot open 'order'\n" );
1546 for( i=0; (j=topol[njob-2][0][i])!=-1; i++ )
1548 fprintf( orderfp, "%d\n", j );
1550 for( i=0; (j=topol[njob-2][1][i])!=-1; i++ )
1552 fprintf( orderfp, "%d\n", j );
1556 if( treeout && noalign )
1558 writeData_pointer( prep_g, njob, name, nlen, seq );
1559 fprintf( stderr, "\n" );
1564 // countnode( njob, topol, node0 );
1569 utree = 0; counteff( njob, topol, len, eff ); utree = 1;
1571 counteff_simple_float( njob, topol, len, eff );
1575 // counteff_simple_float( nkozo, topol_kozo, len_kozo, eff_kozo ); // single weight nanode iranai
1576 for( i=0,j=0; i<njob; i++ )
1580 // eff_kozo_mapped[i] = eff_kozo[j]; //
1581 eff_kozo_mapped[i] = eff[i]; // single weight
1585 eff_kozo_mapped[i] = 0.0;
1586 // fprintf( stderr, "eff_kozo_mapped[%d] = %f\n", i, eff_kozo_mapped[i] );
1587 // fprintf( stderr, " eff[%d] = %f\n", i, eff[i] );
1596 for( i=0; i<njob; i++ ) eff[i] = 1.0;
1599 for( i=0; i<njob; i++ )
1602 eff_kozo_mapped[i] = 1.0;
1604 eff_kozo_mapped[i] = 0.0;
1609 FreeFloatHalfMtx( iscore, njob );
1610 FreeFloatMtx( len );
1612 alloclen = nlenmax*2+1; //chuui!
1613 bseq = AllocateCharMtx( njob, alloclen );
1617 alignmentlength = strlen( seq[0] );
1618 for( i=0; i<njob-nadd; i++ )
1620 if( alignmentlength != strlen( seq[i] ) )
1622 fprintf( stderr, "#################################################################################\n" );
1623 fprintf( stderr, "# ERROR! #\n" );
1624 fprintf( stderr, "# The original%4d sequences must be aligned #\n", njob-nadd );
1625 fprintf( stderr, "#################################################################################\n" );
1631 alignmentlength = strlen( seq[njob-nadd] );
1632 for( i=njob-nadd; i<njob; i++ )
1634 if( alignmentlength != strlen( seq[i] ) )
1636 fprintf( stderr, "###############################################################################\n" );
1637 fprintf( stderr, "# ERROR! #\n" );
1638 fprintf( stderr, "# The%4d additional sequences must be aligned #\n", nadd );
1639 fprintf( stderr, "# Otherwise, try the '--add' option, instead of '--addprofile' option. #\n" );
1640 fprintf( stderr, "###############################################################################\n" );
1644 for( i=0; i<nadd; i++ ) addmem[i] = njob-nadd+i;
1647 for( i=0; i<njob-1; i++ )
1649 if( samemember( topol[i][0], addmem ) ) // jissainiha nai
1651 mergeoralign[i] = '1';
1654 else if( samemember( topol[i][1], addmem ) )
1656 mergeoralign[i] = '2';
1661 mergeoralign[i] = 'n';
1664 if( !foundthebranch )
1666 fprintf( stderr, "###############################################################################\n" );
1667 fprintf( stderr, "# ERROR! #\n" );
1668 fprintf( stderr, "# There is no appropriate position to add the%4d sequences in the guide tree.#\n", nadd );
1669 fprintf( stderr, "# Check whether the%4d sequences form a monophyletic cluster. #\n", nadd );
1670 fprintf( stderr, "# If not, try the '--add' option, instead of the '--addprofile' option. #\n" );
1671 fprintf( stderr, "############################################################################### \n" );
1674 commongappick( nadd, seq+njob-nadd );
1675 for( i=njob-nadd; i<njob; i++ ) strcpy( bseq[i], seq[i] );
1679 for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'n';
1680 for( j=njob-nadd; j<njob; j++ )
1684 for( i=0; i<njob-1; i++ )
1686 if( samemember( topol[i][0], addmem ) ) // arieru
1688 // fprintf( stderr, "HIT!\n" );
1689 if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1690 else mergeoralign[i] = '1';
1692 else if( samemember( topol[i][1], addmem ) )
1694 // fprintf( stderr, "HIT!\n" );
1695 if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1696 else mergeoralign[i] = '2';
1701 for( i=0; i<nadd; i++ ) addmem[i] = njob-nadd+i;
1703 for( i=0; i<njob-1; i++ )
1705 if( includemember( topol[i][0], addmem ) && includemember( topol[i][1], addmem ) )
1707 mergeoralign[i] = 'w';
1709 else if( includemember( topol[i][0], addmem ) )
1711 mergeoralign[i] = '1';
1713 else if( includemember( topol[i][1], addmem ) )
1715 mergeoralign[i] = '2';
1719 for( i=0; i<njob-1; i++ )
1721 fprintf( stderr, "mem0 = " );
1722 for( j=0; topol[i][0][j]>-1; j++ ) fprintf( stderr, "%d ", topol[i][0][j] );
1723 fprintf( stderr, "\n" );
1724 fprintf( stderr, "mem1 = " );
1725 for( j=0; topol[i][1][j]>-1; j++ ) fprintf( stderr, "%d ", topol[i][1][j] );
1726 fprintf( stderr, "\n" );
1727 fprintf( stderr, "i=%d, mergeoralign[] = %c\n", i, mergeoralign[i] );
1730 for( i=njob-nadd; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1733 commongappick( njob-nadd, seq );
1734 for( i=0; i<njob-nadd; i++ ) strcpy( bseq[i], seq[i] );
1738 for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1739 for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'a';
1742 if( rnakozo && rnaprediction == 'm' )
1744 singlerna = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
1745 prep = fopen( "hat4", "r" );
1746 if( prep == NULL ) ErrorExit( "Make hat4 using mccaskill." );
1747 fprintf( stderr, "Loading 'hat4' ... " );
1748 for( i=0; i<njob; i++ )
1750 nogaplen = strlen( bseq[i] );
1751 singlerna[i] = (RNApair **)calloc( nogaplen, sizeof( RNApair * ) );
1752 for( j=0; j<nogaplen; j++ )
1754 singlerna[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
1755 singlerna[i][j][0].bestpos = -1;
1756 singlerna[i][j][0].bestscore = -1.0;
1758 readmccaskill( prep, singlerna[i], nogaplen );
1761 fprintf( stderr, "\ndone.\n" );
1767 fprintf( stderr, "Progressive alignment ... \n" );
1769 #ifdef enablemultithread
1770 if( nthread > 0 && nadd == 0 )
1772 treebasethread_arg_t *targ;
1775 pthread_mutex_t mutex;
1776 pthread_cond_t treecond;
1781 nthread_yoyu = nthread * 1;
1784 targ = calloc( nthread_yoyu, sizeof( treebasethread_arg_t ) );
1785 fftlog = AllocateIntVec( njob );
1786 handle = calloc( nthread_yoyu, sizeof( pthread_t ) );
1787 pthread_mutex_init( &mutex, NULL );
1788 pthread_cond_init( &treecond, NULL );
1790 for( i=0; i<njob; i++ ) dep[i].done = 0;
1791 for( i=0; i<njob; i++ ) fftlog[i] = 1;
1794 calcimportance( njob, eff, bseq, localhomtable );
1796 for( i=0; i<nthread_yoyu; i++ )
1798 targ[i].thread_no = i;
1799 targ[i].nrunpt = &nrun;
1800 targ[i].njob = njob;
1801 targ[i].nlen = nlen;
1802 targ[i].jobpospt = &jobpos;
1803 targ[i].topol = topol;
1805 targ[i].aseq = bseq;
1806 targ[i].effarr = eff;
1807 targ[i].alloclenpt = &alloclen;
1808 targ[i].localhomtable = localhomtable;
1809 targ[i].singlerna = singlerna;
1810 targ[i].effarr_kozo = eff_kozo_mapped;
1811 targ[i].fftlog = fftlog;
1812 targ[i].mutex = &mutex;
1813 targ[i].treecond = &treecond;
1815 pthread_create( handle+i, NULL, treebasethread, (void *)(targ+i) );
1818 for( i=0; i<nthread_yoyu; i++ )
1820 pthread_join( handle[i], NULL );
1822 pthread_mutex_destroy( &mutex );
1823 pthread_cond_destroy( &treecond );
1830 treebase( nlen, bseq, nadd, mergeoralign, mseq1, mseq2, topol, eff, &alloclen, localhomtable, singlerna, eff_kozo_mapped );
1831 fprintf( stderr, "\ndone.\n" );
1834 unweightedspscore = plainscore( njob, bseq );
1835 fprintf( stderr, "\nSCORE %s = %.0f, ", "(treebase)", unweightedspscore );
1836 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
1837 fprintf( stderr, "\n\n" );
1843 LocalHom *tmppt1, *tmppt2;
1844 for( i=0; i<njob; i++)
1846 for( j=0; j<njob; j++)
1848 tmppt1 = localhomtable[i]+j;
1849 while( tmppt2 = tmppt1->next )
1851 free( (void *)tmppt1 );
1854 free( (void *)tmppt1 );
1856 free( (void *)(localhomtable[i]+j) );
1858 free( (void *)localhomtable );
1862 fprintf( trap_g, "done.\n" );
1864 free( mergeoralign );
1866 writeData_pointer( prep_g, njob, name, nlen, bseq );
1868 writeData( stdout, njob, name, nlen, bseq );
1869 writePre( njob, name, nlen, bseq, !contin );
1870 writeData_pointer( prep_g, njob, name, nlen, aseq );
1873 fprintf( stderr, "OSHIMAI\n" );
1876 if( constraint ) FreeLocalHomTable( localhomtable, njob );