X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=binaries%2Fsrc%2Fmafft%2Fcore%2FpartQalignmm.c;fp=binaries%2Fsrc%2Fmafft%2Fcore%2FpartQalignmm.c;h=0000000000000000000000000000000000000000;hb=6e0ce943f09b5ac30f3eb8dc0f20bc75114669ce;hp=162df7a030db2812783cdf2546dbaec7dffff884;hpb=48f0451143441afce738d887b0cde1fb5cdbc02e;p=jabaws.git diff --git a/binaries/src/mafft/core/partQalignmm.c b/binaries/src/mafft/core/partQalignmm.c deleted file mode 100644 index 162df7a..0000000 --- a/binaries/src/mafft/core/partQalignmm.c +++ /dev/null @@ -1,1258 +0,0 @@ -#include "mltaln.h" -#include "dp.h" - -#define MACHIGAI 0 -#define OUTGAP0TRY 1 -#define DEBUG 0 -#define XXXXXXX 0 -#define USE_PENALTY_EX 0 -#define FASTMATCHCALC 1 - - -static int impalloclen = 0; -static float **impmtx = NULL; -float part_imp_match_out_scQ( int i1, int j1 ) -{ -// fprintf( stderr, "impalloclen = %d\n", impalloclen ); -// fprintf( stderr, "i1,j1=%d,%d -> impmtx=%f\n", i1, j1, impmtx[i1][j1] ); - return( impmtx[i1][j1] ); -#if 0 - if( i1 == l1 || j1 == l2 ) return( 0.0 ); - return( impmtx[i1+start1][j1+start2] ); -#endif -} -static void part_imp_match_out_vead_gapmapQ( float *imp, int i1, int lgth2, int start2, int *gapmap2 ) -{ -#if FASTMACHCALC - float *pt = imp; - int *gapmappt = gapmap2; - while( lgth2-- ) - *pt++ += impmtx[i1][start2+*gapmappt++]; -#else - int j; - for( j=0; jstart1 ); -// fprintf( stderr, "end1 = %d\n", tmpptr->end1 ); -// fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] ); -// fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] ); - pt = seq1[i]; - tmpint = -1; - while( *pt != 0 ) - { - if( *pt++ != '-' ) tmpint++; - if( tmpint == tmpptr->start1 ) break; - } - start1 = (int)( pt - seq1[i] ) - 1; - - if( tmpptr->start1 == tmpptr->end1 ) end1 = start1; - else - { -#if MACHIGAI - while( *pt != 0 ) - { - if( tmpint == tmpptr->end1 ) break; - if( *pt++ != '-' ) tmpint++; - } - end1 = (int)( pt - seq1[i] ) - 1; -#else - while( *pt != 0 ) - { -// fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] ); - if( *pt++ != '-' ) tmpint++; - if( tmpint == tmpptr->end1 ) break; - } - end1 = (int)( pt - seq1[i] ) - 1; -#endif - } - - pt = seq2[j]; - tmpint = -1; - while( *pt != 0 ) - { - if( *pt++ != '-' ) tmpint++; - if( tmpint == tmpptr->start2 ) break; - } - start2 = (int)( pt - seq2[j] ) - 1; - if( tmpptr->start2 == tmpptr->end2 ) end2 = start2; - else - { -#if MACHIGAI - while( *pt != 0 ) - { - if( tmpint == tmpptr->end2 ) break; - if( *pt++ != '-' ) tmpint++; - } - end2 = (int)( pt - seq2[j] ) - 1; -#else - while( *pt != 0 ) - { - if( *pt++ != '-' ) tmpint++; - if( tmpint == tmpptr->end2 ) break; - } - end2 = (int)( pt - seq2[j] ) - 1; -#endif - } -// fprintf( stderr, "start1 = %d (%c), end1 = %d (%c), start2 = %d (%c), end2 = %d (%c)\n", start1, seq1[i][start1], end1, seq1[i][end1], start2, seq2[j][start2], end2, seq2[j][end2] ); -// fprintf( stderr, "step 0\n" ); - if( end1 - start1 != end2 - start2 ) - { -// fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 ); - } - - k1 = start1; k2 = start2; - pt1 = seq1[i] + k1; - pt2 = seq2[j] + k2; - while( *pt1 && *pt2 ) - { - if( *pt1 != '-' && *pt2 != '-' ) - { -// ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃí°Õ¤·¤Æ²¼¤µ¤¤¡£ -// impmtx[k1][k2] += tmpptr->wimportance * fastathreshold; -// impmtx[k1][k2] += tmpptr->importance * effij; - impmtx[k1][k2] += tmpptr->fimportance * effij; -// fprintf( stderr, "k1=%d, k2=%d, impalloclen=%d\n", k1, k2, impalloclen ); -// fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); - k1++; k2++; - pt1++; pt2++; - } - else if( *pt1 != '-' && *pt2 == '-' ) - { -// fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); - k2++; pt2++; - } - else if( *pt1 == '-' && *pt2 != '-' ) - { -// fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); - k1++; pt1++; - } - else if( *pt1 == '-' && *pt2 == '-' ) - { -// fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); - k1++; pt1++; - k2++; pt2++; - } - if( k1 > end1 || k2 > end2 ) break; - } - tmpptr = tmpptr->next; - } - } - } -#if 0 - fprintf( stderr, "impmtx = \n" ); - for( k2=0; k2start1 ); - fprintf( stderr, "end1 = %d\n", localhom[i][j]->end1 ); - fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] ); - pt = seq1[i]; - tmpint = -1; - while( *pt != 0 ) - { - if( *pt++ != '-' ) tmpint++; - if( tmpint == localhom[i][j]->start1 ) break; - } - start1 = pt - seq1[i] - 1; - - while( *pt != 0 ) - { -// fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] ); - if( *pt++ != '-' ) tmpint++; - if( tmpint == localhom[i][j]->end1 ) break; - } - end1 = pt - seq1[i] - 1; - - pt = seq2[j]; - tmpint = -1; - while( *pt != 0 ) - { - if( *pt++ != '-' ) tmpint++; - if( tmpint == localhom[i][j]->start2 ) break; - } - start2 = pt - seq2[j] - 1; - while( *pt != 0 ) - { - if( *pt++ != '-' ) tmpint++; - if( tmpint == localhom[i][j]->end2 ) break; - } - end2 = pt - seq2[j] - 1; -// fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 ); - k1 = start1; - k2 = start2; - fprintf( stderr, "step 0\n" ); - while( k1 <= end1 && k2 <= end2 ) - { -#if 0 - if( !nocount1[k1] && !nocount2[k2] ) - impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j]; - k1++; k2++; -#else - if( !nocount1[k1] && !nocount2[k2] ) - impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j]; - k1++; k2++; -#endif - } - - dif = ( end1 - start1 ) - ( end2 - start2 ); - fprintf( stderr, "dif = %d\n", dif ); - if( dif > 0 ) - { - do - { - fprintf( stderr, "dif = %d\n", dif ); - k1 = start1; - k2 = start2 - dif; - while( k1 <= end1 && k2 <= end2 ) - { - if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] ) - impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j]; - k1++; k2++; - } - } - while( dif-- ); - } - else - { - do - { - k1 = start1 + dif; - k2 = start2; - while( k1 <= end1 ) - { - if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] ) - impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j]; - k1++; k2++; - } - } - while( dif++ ); - } - } - } -#if 0 - fprintf( stderr, "impmtx = \n" ); - for( k2=0; k2-1 ) - *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++; - matchpt++; - } - } -#else - int j, k, l; - float scarr[26]; - float **cpmxpd = floatwork; - int **cpmxpdn = intwork; - // simple - if( initialize ) - { - int count = 0; - for( j=0; j-1; k++ ) - match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j]; - } -#endif -} - -static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw, - char **seq1, char **seq2, - char **mseq1, char **mseq2, - float **cpmx1, float **cpmx2, - int **ijp, int icyc, int jcyc, - int start1, int end1, int start2, int end2, - int *gapmap1, int *gapmap2 ) -{ - int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k; - char gap[] = "-"; - float wm; - lgth1 = strlen( seq1[0] ); - lgth2 = strlen( seq2[0] ); - -#if 0 - for( i=0; i= wm ) - { - wm = lastverticalw[i]; - iin = i; jin = lgth2-1; - ijp[lgth1][lgth2] = +( lgth1 - i ); - } - } - for( j=0; j= wm ) - { - wm = lasthorizontalw[j]; - iin = lgth1-1; jin = j; - ijp[lgth1][lgth2] = -( lgth2 - j ); - } - } - } - - for( i=0; i 0 ) - { - ifi = iin-ijp[iin][jin]; jfi = jin-1; - } - else - { - ifi = iin-1; jfi = jin-1; - } - l = iin - ifi; - while( --l ) - { - for( i=0; i= wm ) - { - wm = lastverticalw[i]; - iin = i; jin = lgth2-1; - ijp[lgth1][lgth2] = +( lgth1 - i ); - } - } - for( j=0; j= wm ) - { - wm = lasthorizontalw[j]; - iin = lgth1-1; jin = j; - ijp[lgth1][lgth2] = -( lgth2 - j ); - } - } - } - - for( i=0; i 0 ) - { - ifi = iin-ijp[iin][jin]; jfi = jin-1; - } - else - { - ifi = iin-1; jfi = jin-1; - } - l = iin - ifi; - while( --l ) - { - for( i=0; i lgth1, outgap == 1 -> lgth1+1 */ - int lgth1, lgth2; - int resultlen; - float wm = 0.0; /* int ?????? */ - float g; - float *currentw, *previousw; -#if 1 - float *wtmp; - int *ijppt; - float *mjpt, *prept, *curpt; - int *mpjpt; -#endif - static float mi, *m; - static int **ijp; - static int mpi, *mp; - static float *w1, *w2; - static float *match; - static float *initverticalw; /* kufuu sureba iranai */ - static float *lastverticalw; /* kufuu sureba iranai */ - static char **mseq1; - static char **mseq2; - static char **mseq; - static float *digf1; - static float *digf2; - static float *diaf1; - static float *diaf2; - static float *gapz1; - static float *gapz2; - static float *gapf1; - static float *gapf2; - static float *ogcp1g; - static float *ogcp2g; - static float *fgcp1g; - static float *fgcp2g; - static float *og_h_dg_n1_p; - static float *og_h_dg_n2_p; - static float *fg_h_dg_n1_p; - static float *fg_h_dg_n2_p; - static float *og_t_fg_h_dg_n1_p; - static float *og_t_fg_h_dg_n2_p; - static float *fg_t_og_h_dg_n1_p; - static float *fg_t_og_h_dg_n2_p; - static float *gapz_n1; - static float *gapz_n2; - static float **cpmx1; - static float **cpmx2; - static int **intwork; - static float **floatwork; - static int orlgth1 = 0, orlgth2 = 0; - float fpenalty = (float)penalty; -#if USE_PENALTY_EX - float fpenalty_ex = (float)penalty_ex; -#endif - float tmppenal; - float *fg_t_og_h_dg_n2_p_pt; - float *og_t_fg_h_dg_n2_p_pt; - float *og_h_dg_n2_p_pt; - float *fg_h_dg_n2_p_pt; - float *gapz_n2_pt0; - float *gapz_n2_pt1; - float *fgcp2pt; - float *ogcp2pt; - float fg_t_og_h_dg_n1_p_va; - float og_t_fg_h_dg_n1_p_va; - float og_h_dg_n1_p_va; - float fg_h_dg_n1_p_va; - float gapz_n1_va0; - float gapz_n1_va1; - float fgcp1va; - float ogcp1va; - - - -#if 0 - fprintf( stderr, "eff in SA+++align\n" ); - for( i=0; i orlgth1 || lgth2 > orlgth2 ) - { - int ll1, ll2; - - if( orlgth1 > 0 && orlgth2 > 0 ) - { - FreeFloatVec( w1 ); - FreeFloatVec( w2 ); - FreeFloatVec( match ); - FreeFloatVec( initverticalw ); - FreeFloatVec( lastverticalw ); - - FreeFloatVec( m ); - FreeIntVec( mp ); - - FreeCharMtx( mseq ); - - FreeFloatVec( digf1 ); - FreeFloatVec( digf2 ); - FreeFloatVec( diaf1 ); - FreeFloatVec( diaf2 ); - FreeFloatVec( gapz1 ); - FreeFloatVec( gapz2 ); - FreeFloatVec( gapf1 ); - FreeFloatVec( gapf2 ); - FreeFloatVec( ogcp1g ); - FreeFloatVec( ogcp2g ); - FreeFloatVec( fgcp1g ); - FreeFloatVec( fgcp2g ); - FreeFloatVec( og_h_dg_n1_p ); - FreeFloatVec( og_h_dg_n2_p ); - FreeFloatVec( fg_h_dg_n1_p ); - FreeFloatVec( fg_h_dg_n2_p ); - FreeFloatVec( og_t_fg_h_dg_n1_p ); - FreeFloatVec( og_t_fg_h_dg_n2_p ); - FreeFloatVec( fg_t_og_h_dg_n1_p ); - FreeFloatVec( fg_t_og_h_dg_n2_p ); - FreeFloatVec( gapz_n1 ); - FreeFloatVec( gapz_n2 ); - - FreeFloatMtx( cpmx1 ); - FreeFloatMtx( cpmx2 ); - - FreeFloatMtx( floatwork ); - FreeIntMtx( intwork ); - } - - ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100; - ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100; - -#if DEBUG - fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 ); -#endif - - w1 = AllocateFloatVec( ll2+2 ); - w2 = AllocateFloatVec( ll2+2 ); - match = AllocateFloatVec( ll2+2 ); - - initverticalw = AllocateFloatVec( ll1+2 ); - lastverticalw = AllocateFloatVec( ll1+2 ); - - m = AllocateFloatVec( ll2+2 ); - mp = AllocateIntVec( ll2+2 ); - - mseq = AllocateCharMtx( njob, ll1+ll2 ); - - digf1 = AllocateFloatVec( ll1+2 ); - digf2 = AllocateFloatVec( ll2+2 ); - diaf1 = AllocateFloatVec( ll1+2 ); - diaf2 = AllocateFloatVec( ll2+2 ); - gapz1 = AllocateFloatVec( ll1+2 ); - gapz2 = AllocateFloatVec( ll2+2 ); - gapf1 = AllocateFloatVec( ll1+2 ); - gapf2 = AllocateFloatVec( ll2+2 ); - ogcp1g = AllocateFloatVec( ll1+2 ); - ogcp2g = AllocateFloatVec( ll2+2 ); - fgcp1g = AllocateFloatVec( ll1+2 ); - fgcp2g = AllocateFloatVec( ll2+2 ); - og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 ); - og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 ); - fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 ); - fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 ); - og_t_fg_h_dg_n1_p = AllocateFloatVec( ll1 + 2 ); - og_t_fg_h_dg_n2_p = AllocateFloatVec( ll2 + 2 ); - fg_t_og_h_dg_n1_p = AllocateFloatVec( ll1 + 2 ); - fg_t_og_h_dg_n2_p = AllocateFloatVec( ll2 + 2 ); - gapz_n1 = AllocateFloatVec( ll1+2 ); - gapz_n2 = AllocateFloatVec( ll2+2 ); - - cpmx1 = AllocateFloatMtx( 26, ll1+2 ); - cpmx2 = AllocateFloatMtx( 26, ll2+2 ); - -#if FASTMATCHCALC - floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); - intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); -#else - floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 ); - intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 ); -#endif - -#if DEBUG - fprintf( stderr, "succeeded\n" ); -#endif - - orlgth1 = ll1 - 100; - orlgth2 = ll2 - 100; - } - - - for( i=0; i commonAlloc1 || orlgth2 > commonAlloc2 ) - { - int ll1, ll2; - - if( commonAlloc1 && commonAlloc2 ) - { - FreeIntMtx( commonIP ); - } - - ll1 = MAX( orlgth1, commonAlloc1 ); - ll2 = MAX( orlgth2, commonAlloc2 ); - -#if DEBUG - fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 ); -#endif - - commonIP = AllocateIntMtx( ll1+10, ll2+10 ); - -#if DEBUG - fprintf( stderr, "succeeded\n\n" ); -#endif - - commonAlloc1 = ll1; - commonAlloc2 = ll2; - } - ijp = commonIP; - - cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc ); - cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc ); - - if( sgap1 ) - { - new_OpeningGapCount_zure( ogcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 ); - new_OpeningGapCount_zure( ogcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 ); - new_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1, sgap1, egap1 ); - new_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2, sgap2, egap2 ); - getdigapfreq_part( digf1, icyc, seq1, eff1, lgth1, sgap1, egap1 ); - getdigapfreq_part( digf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 ); - getdiaminofreq_part( diaf1, icyc, seq1, eff1, lgth1, sgap1, egap1 ); - getdiaminofreq_part( diaf2, jcyc, seq2, eff2, lgth2, sgap2, egap2 ); - getgapfreq( gapf1, icyc, seq1, eff1, lgth1 ); - getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 ); - getgapfreq_zure_part( gapz1, icyc, seq1, eff1, lgth1, sgap1 ); - getgapfreq_zure_part( gapz2, jcyc, seq2, eff2, lgth2, sgap1 ); - } - else - { - st_OpeningGapCount( ogcp1g, icyc, seq1, eff1, lgth1 ); - st_OpeningGapCount( ogcp2g, jcyc, seq2, eff2, lgth2 ); - st_FinalGapCount_zure( fgcp1g, icyc, seq1, eff1, lgth1 ); - st_FinalGapCount_zure( fgcp2g, jcyc, seq2, eff2, lgth2 ); - getdigapfreq_st( digf1, icyc, seq1, eff1, lgth1 ); - getdigapfreq_st( digf2, jcyc, seq2, eff2, lgth2 ); - getdiaminofreq_x( diaf1, icyc, seq1, eff1, lgth1 ); - getdiaminofreq_x( diaf2, jcyc, seq2, eff2, lgth2 ); - getgapfreq( gapf1, icyc, seq1, eff1, lgth1 ); - getgapfreq( gapf2, jcyc, seq2, eff2, lgth2 ); - getgapfreq_zure( gapz1, icyc, seq1, eff1, lgth1 ); - getgapfreq_zure( gapz2, jcyc, seq2, eff2, lgth2 ); - } - -#if 1 - lastj = lgth2+2; - for( i=0; i tbfast.c - if( localhom ) - imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 ); - -#endif - - if( outgap == 1 ) - { - g = 0.0; - - g += ogcp1g[0] * og_h_dg_n2_p[0]; - - g += ogcp2g[0] * og_h_dg_n1_p[0]; - - g += fgcp1g[0] * fg_h_dg_n2_p[0]; - - g += fgcp2g[0] * fg_h_dg_n1_p[0]; - - initverticalw[0] += g; - currentw[0] += g; - - for( i=1; i", wm ); -#endif - tmppenal = gapz_n1_va1 * *fg_t_og_h_dg_n2_p_pt; -#if 0 - fprintf( stderr, "%5.0f?", g ); -#endif - if( (g=mi+tmppenal) > wm ) - { - wm = g; - *ijppt = -( j - mpi ); - } - - tmppenal = gapz_n1_va0 * *og_t_fg_h_dg_n2_p_pt; - if( (g=*prept+tmppenal) >= mi ) - { - mi = g; - mpi = j-1; - } -#if USE_PENALTY_EX - mi += fpenalty_ex; -#endif - - tmppenal = *gapz_n2_pt1 * fg_t_og_h_dg_n1_p_va; -#if 0 - fprintf( stderr, "%5.0f?", g ); -#endif - if( (g=*mjpt+tmppenal) > wm ) - { - wm = g; - *ijppt = +( i - *mpjpt ); - } - - tmppenal = *gapz_n2_pt0 * og_t_fg_h_dg_n1_p_va; - if( (g=*prept+tmppenal) >= *mjpt ) - { - *mjpt = g; - *mpjpt = i-1; - } -#if USE_PENALTY_EX - m[j] += fpenalty_ex; -#endif - -#if 0 - fprintf( stderr, "%5.0f ", wm ); -#endif - *curpt++ += wm; - ijppt++; - mjpt++; - prept++; - mpjpt++; - fg_t_og_h_dg_n2_p_pt++; - og_t_fg_h_dg_n2_p_pt++; - og_h_dg_n2_p_pt++; - fg_h_dg_n2_p_pt++; - gapz_n2_pt0++; - gapz_n2_pt1++; - fgcp2pt++; - ogcp2pt++; - } - lastverticalw[i] = currentw[lgth2-1]; - } - -#if OUTGAP0TRY - if( !outgap ) - { - for( j=1; j" ); - for( i=0; i N ) - { - fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N ); - ErrorExit( "LENGTH OVER!\n" ); - } - - - for( i=0; i