static int treebase( int *nlen, char **aseq, int nadd, char *mergeoralign, char **mseq1, char **mseq2, int ***topol, Treedep *dep, int **memhist, double ***cpmxhist, double *effarr, double **newdistmtx, int *selfscore, int *alloclen, int (*callback)(int, int, char*) ) { int l, len1, len2, i, m, immin, immax; int len1nocommongap, len2nocommongap; int clus1, clus2; double pscore, tscore; char *indication1 = NULL, *indication2 = NULL; double *effarr1 = NULL; double *effarr2 = NULL; int *fftlog = NULL; // fixed at 2006/07/26 // double dumfl = 0.0; double dumdb = 0.0; int ffttry; int m1, m2; int *gaplen = NULL; int *gapmap = NULL; int *alreadyaligned = NULL; double **dynamicmtx = NULL; double ssi, ssm, bunbo; int tm, ti; int gapmaplen; int **localmem = NULL; double **cpmxchild0, **cpmxchild1; double orieff1, orieff2; #if SKIP int **skiptable1 = NULL, **skiptable2 = NULL; #endif #if 0 int i, j; #endif if( effarr1 == NULL ) { effarr1 = AllocateDoubleVec( njob ); effarr2 = AllocateDoubleVec( njob ); indication1 = AllocateCharVec( 150 ); indication2 = AllocateCharVec( 150 ); fftlog = AllocateIntVec( njob ); gaplen = AllocateIntVec( *alloclen+10 ); gapmap = AllocateIntVec( *alloclen+10 ); alreadyaligned = AllocateIntVec( njob ); if( specificityconsideration ) dynamicmtx = AllocateDoubleMtx( nalphabets, nalphabets ); localmem = calloc( sizeof( int * ), 2 ); } for( i=0; i 0 && dep[l].child0 == l-1 && dep[l].child1 == -1 && dep[dep[l].child0].child1 == -1 ) { localmem[0][clus1] = topol[l-1][1][0]; localmem[0][clus1+1] = -1; localmem[1][0] = topol[l][1][0]; localmem[1][1] = -1; } else { localmem[0][0] = -1; posinmem = topolorderz( localmem[0], topol, dep, l, 0 ) - localmem[0]; localmem[1][0] = -1; posinmem = topolorderz( localmem[1], topol, dep, l, 1 ) - localmem[1]; } #else if( dep[l].child0 == -1 ) { localmem[0] = calloc( sizeof( int ), 2 ); localmem[0][0] = m1; localmem[0][1] = -1; clus1 = 1; } else { localmem[0] = memhist[dep[l].child0]; clus1 = intlen( localmem[0] ); } if( dep[l].child1 == -1 ) { localmem[1] = calloc( sizeof( int ), 2 ); localmem[1][0] = m2; localmem[1][1] = -1; clus2 = 1; } else { localmem[1] = memhist[dep[l].child1]; clus2 = intlen( localmem[1] ); } if( l != njob-2 ) { memhist[l] = calloc( sizeof( int ), clus1+clus2+1 ); intcpy( memhist[l], localmem[0] ); intcpy( memhist[l]+clus1, localmem[1] ); memhist[l][clus1+clus2] = -1; } #endif if( mergeoralign[l] == 'n' ) { // reporterr( "SKIP!\n" ); // free( topol[l][0] ); topol[l][0] = NULL; // free( topol[l][1] ); topol[l][1] = NULL; // free( topol[l] ); topol[l] = NULL; continue; } // reporterr( "\ndistfromtip = %f\n", dep[l].distfromtip ); if( specificityconsideration ) makedynamicmtx( dynamicmtx, n_dis_consweight_multi, dep[l].distfromtip ); else dynamicmtx = n_dis_consweight_multi; // makedynamicmtx( dynamicmtx, n_dis_consweight_multi, ( dep[l].distfromtip - 0.2 ) * 3 ); len1 = strlen( aseq[m1] ); len2 = strlen( aseq[m2] ); if( *alloclen < len1 + len2 ) { reporterr( "\nReallocating.." ); *alloclen = ( len1 + len2 ) + 1000; ReallocateCharMtx( aseq, njob, *alloclen + 10 ); gaplen = realloc( gaplen, ( *alloclen + 10 ) * sizeof( int ) ); if( gaplen == NULL ) { reporterr( "Cannot realloc gaplen\n" ); exit( 1 ); } gapmap = realloc( gapmap, ( *alloclen + 10 ) * sizeof( int ) ); if( gapmap == NULL ) { reporterr( "Cannot realloc gapmap\n" ); exit( 1 ); } reporterr( "done. *alloclen = %d\n", *alloclen ); } #if 1 // CHUUI@@@@ clus1 = fastconjuction_noname( localmem[0], aseq, mseq1, effarr1, effarr, indication1, 0.0, &orieff1 ); clus2 = fastconjuction_noname( localmem[1], aseq, mseq2, effarr2, effarr, indication2, 0.0, &orieff2 ); #else clus1 = fastconjuction_noname( topol[l][0], aseq, mseq1, effarr1, effarr, indication1, 0.0 ); clus2 = fastconjuction_noname( topol[l][1], aseq, mseq2, effarr2, effarr, indication2, 0.0 ); // clus1 = fastconjuction_noweight( topol[l][0], aseq, mseq1, effarr1, indication1 ); // clus2 = fastconjuction_noweight( topol[l][1], aseq, mseq2, effarr2, indication2 ); #endif if( mergeoralign[l] == '1' || mergeoralign[l] == '2' ) { newgapstr = "="; } else newgapstr = "-"; len1nocommongap = len1; len2nocommongap = len2; if( mergeoralign[l] == '1' ) // nai { findcommongaps( clus2, mseq2, gapmap ); commongappick( clus2, mseq2 ); len2nocommongap = strlen( mseq2[0] ); } else if( mergeoralign[l] == '2' ) { findcommongaps( clus1, mseq1, gapmap ); commongappick( clus1, mseq1 ); len1nocommongap = strlen( mseq1[0] ); } #if 0 for( i=0; i 66 ) reporterr( "..." ); reporterr( "\n" ); reporterr( "group2 = %.66s", indication2 ); if( strlen( indication2 ) > 66 ) reporterr( "..." ); reporterr( "\n" ); #endif /* reporterr( "before align all\n" ); display( aseq, njob ); reporterr( "\n" ); reporterr( "before align 1 %s \n", indication1 ); display( mseq1, clus1 ); reporterr( "\n" ); reporterr( "before align 2 %s \n", indication2 ); display( mseq2, clus2 ); reporterr( "\n" ); */ if( !nevermemsave && ( alg != 'M' && ( len1 > 30000 || len2 > 30000 ) ) ) { reporterr( "\nlen1=%d, len2=%d, Switching to the memsave mode\n", len1, len2 ); alg = 'M'; if( commonIP ) FreeIntMtx( commonIP ); commonIP = NULL; commonAlloc1 = 0; commonAlloc2 = 0; } // if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 ); if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000); else ffttry = 0; // ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000); // v6.708 // reporterr( "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (double)len1/fftlog[m1], clus1, (double)len2/fftlog[m2], clus2 ); if( force_fft || ( use_fft && ffttry ) ) { if( l < 500 || l % 100 == 0 ) reporterr( " f\b\b" ); if( alg == 'M' ) { if( l < 500 || l % 100 == 0 ) reporterr( "m" ); pscore = Falign_udpari_long( NULL, NULL, dynamicmtx, mseq1, mseq2, effarr1, effarr2, NULL, NULL, clus1, clus2, *alloclen, fftlog+m1 ); } else { pscore = Falign( NULL, NULL, dynamicmtx, mseq1, mseq2, effarr1, effarr2, NULL, NULL, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL ); // reporterr( "######### mseq1[0] = %s\n", mseq1[0] ); } } else { if( l < 500 || l % 100 == 0 ) reporterr( " d\b\b" ); fftlog[m1] = 0; switch( alg ) { case( 'a' ): pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen ); break; case( 'M' ): if( l < 500 || l % 100 == 0 ) reporterr( "m" ); if( l < 500 || l % 100 == 0 ) if( cpmxchild1 || cpmxchild0 ) reporterr( " h" ); // reporterr( "%d-%d", clus1, clus2 ); pscore = MSalignmm( dynamicmtx, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap, cpmxchild0, cpmxchild1, cpmxhist+l, orieff1, orieff2 ); break; case( 'd' ): if( 1 && clus1 == 1 && clus2 == 1 ) { // reporterr( "%d-%d", clus1, clus2 ); pscore = G__align11( dynamicmtx, mseq1, mseq2, *alloclen, outgap, outgap ); } else { // reporterr( "%d-%d", clus1, clus2 ); pscore = D__align_ls( dynamicmtx, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, 0, &dumdb, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap ); } break; case( 'A' ): if( clus1 == 1 && clus2 == 1 ) { // reporterr( "%d-%d", clus1, clus2 ); pscore = G__align11( dynamicmtx, mseq1, mseq2, *alloclen, outgap, outgap ); } else { if( l < 500 || l % 100 == 0 ) if( cpmxchild1 || cpmxchild0 ) reporterr( " h" ); // reporterr( "\n\n %d - %d (%d x %d) : \n", topol[l][0][0], topol[l][1][0], clus1, clus2 ); pscore = A__align( dynamicmtx, penalty, penalty_ex, mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, 0, &dumdb, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap, localmem[0][0], 1, cpmxchild0, cpmxchild1, cpmxhist+l, orieff1, orieff2 ); } break; default: ErrorExit( "ERROR IN SOURCE FILE" ); } } #if SCOREOUT reporterr( "score = %10.2f\n", pscore ); #endif tscore += pscore; nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] ); // writePre( njob, name, nlen, aseq, 0 ); if( disp ) display( aseq, njob ); // reporterr( "\n" ); if( mergeoralign[l] == '1' ) // jissainiha nai. atarashii hairetsu ha saigo dakara. { reporterr( "Check source!!!\n" ); exit( 1 ); } if( mergeoralign[l] == '2' ) { // if( localkeeplength ) ndeleted += deletenewinsertions( clus1, clus2, mseq1, mseq2, NULL ); // for( i=0; iSTEP0 mseq1[%d] = \n%s\n", i, mseq1[i] ); // for( i=0; iSTEP0 mseq2[%d] = \n%s\n", i, mseq2[i] ); gapmaplen = strlen( mseq1[0] )-len1nocommongap+len1; adjustgapmap( gapmaplen, gapmap, mseq1[0] ); #if 0 reporterr( "\n" ); for( i=0; iSTEP1 mseq1[%d] = \n%s\n", i, mseq1[i] ); for( i=0; iSTEP1 mseq2[%d] = \n%s\n", i, mseq2[i] ); #endif // if( clus1 + clus2 < njob ) restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen, '-' ); if( smoothing ) { restorecommongapssmoothly( njob, njob-(clus1+clus2), aseq, localmem[0], localmem[1], gapmap, *alloclen, '-' ); findnewgaps( clus1, 0, mseq1, gaplen ); insertnewgaps_bothorders( njob, alreadyaligned, aseq, localmem[0], localmem[1], gaplen, gapmap, gapmaplen, *alloclen, alg, '-' ); } else { restorecommongaps( njob, njob-(clus1+clus2), aseq, localmem[0], localmem[1], gapmap, *alloclen, '-' ); findnewgaps( clus1, 0, mseq1, gaplen ); insertnewgaps( njob, alreadyaligned, aseq, localmem[0], localmem[1], gaplen, gapmap, *alloclen, alg, '-' ); } #if 0 reporterr( "\n" ); for( i=0; iSTEP3 mseq1[%d] = \n%s\n", i, mseq1[i] ); for( i=0; iSTEP3 mseq2[%d] = \n%s\n", i, mseq2[i] ); #endif #if 0 for( i=0; i-1; i++ ) alreadyaligned[m] = 1; } if( newdistmtx ) // tsukawanai { #if 0 reporterr( "group1 = " ); for( i=0; i514\n%s\n", aseq[514] ); free( localmem[0] ); free( localmem[1] ); } #if SCOREOUT reporterr( "totalscore = %10.2f\n\n", tscore ); #endif Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL ); Falign_udpari_long( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL ); A__align( NULL, 0, 0, NULL, NULL, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0, -1, -1, NULL, NULL, NULL, 0.0, 0.0 ); D__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); G__align11( NULL, NULL, NULL, 0, 0, 0 ); // iru? free( effarr1 ); free( effarr2 ); free( indication1 ); free( indication2 ); free( fftlog ); free( gaplen ); free( gapmap ); if( specificityconsideration ) FreeDoubleMtx( dynamicmtx ); free( alreadyaligned ); free( localmem ); effarr1 = NULL; return( 0 ); chudan_tbfast: Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL ); Falign_udpari_long( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL ); A__align( NULL, 0, 0, NULL, NULL, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0, -1, -1, NULL, NULL, NULL, 0.0, 0.0 ); D__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 ); G__align11( NULL, NULL, NULL, 0, 0, 0 ); // iru? if( effarr1 ) free( effarr1 ); effarr1 = NULL; if( effarr2 ) free( effarr2 ); effarr2 = NULL; if( indication1 ) free( indication1 ); indication1 = NULL; if( indication2 ) free( indication2 ); indication2 = NULL; if( fftlog ) free( fftlog ); fftlog = NULL; if( gaplen ) free( gaplen ); gaplen = NULL; if( gapmap ) free( gapmap ); gapmap = NULL; if( alreadyaligned ) free( alreadyaligned ); alreadyaligned = NULL; if( specificityconsideration ) { if( dynamicmtx ) FreeDoubleMtx( dynamicmtx ); dynamicmtx = NULL; } if( localmem ) free( localmem ); localmem = NULL; #if SKIP if( skiptable1 ) FreeIntMtx( skiptable1 ); skiptable1 = NULL; if( skiptable2 ) FreeIntMtx( skiptable2 ); skiptable2 = NULL; #endif return( 1 ); }