+++ /dev/null
-#include "mltaln.h"
-
-#define SEGMENTSIZE 150
-#define TMPTMPTMP 0
-
-#define DEBUG 0
-
-void keika( char *str, int current, int all )
-{
- if( current == 0 )
- fprintf( stderr, "%s : ", str );
-
- fprintf( stderr, "\b\b\b\b\b\b\b\b" );
- fprintf( stderr, "%3d /%3d", current+1, all+1 );
-
- if( current+1 == all )
- fprintf( stderr, "\b\b\b\b\b\b\b\bdone. \n" );
-}
-
-double maxItch( double *soukan, int size )
-{
- int i;
- double value = 0.0;
- double cand;
- for( i=0; i<size; i++ )
- if( ( cand = soukan[i] ) > value ) value = cand;
- return( value );
-}
-
-void calcNaiseki( Fukusosuu *value, Fukusosuu *x, Fukusosuu *y )
-{
- value->R = x->R * y->R + x->I * y->I;
- value->I = -x->R * y->I + x->I * y->R;
-}
-
-Fukusosuu *AllocateFukusosuuVec( int l1 )
-{
- Fukusosuu *value;
- value = (Fukusosuu *)calloc( l1, sizeof( Fukusosuu ) );
- if( !value )
- {
- fprintf( stderr, "Cannot allocate %d FukusosuuVec\n", l1 );
- return( NULL );
- }
- return( value );
-}
-
-Fukusosuu **AllocateFukusosuuMtx( int l1, int l2 )
-{
- Fukusosuu **value;
- int j;
-// fprintf( stderr, "allocating %d x %d FukusosuuMtx\n", l1, l2 );
- value = (Fukusosuu **)calloc( l1+1, sizeof( Fukusosuu * ) );
- if( !value )
- {
- fprintf( stderr, "Cannot allocate %d x %d FukusosuuVecMtx\n", l1, l2 );
- exit( 1 );
- }
- for( j=0; j<l1; j++ )
- {
- value[j] = AllocateFukusosuuVec( l2 );
- if( !value[j] )
- {
- fprintf( stderr, "Cannot allocate %d x %d FukusosuuVecMtx\n", l1, l2 );
- exit( 1 );
- }
- }
- value[l1] = NULL;
- return( value );
-}
-
-Fukusosuu ***AllocateFukusosuuCub( int l1, int l2, int l3 )
-{
- Fukusosuu ***value;
- int i;
- value = calloc( l1+1, sizeof( Fukusosuu ** ) );
- if( !value ) ErrorExit( "Cannot allocate Fukusosuu" );
- for( i=0; i<l1; i++ ) value[i] = AllocateFukusosuuMtx( l2, l3 );
- value[l1] = NULL;
- return( value );
-}
-
-void FreeFukusosuuVec( Fukusosuu *vec )
-{
- free( (void *)vec );
-}
-
-void FreeFukusosuuMtx( Fukusosuu **mtx )
-{
- int i;
-
- for( i=0; mtx[i]; i++ )
- free( (void *)mtx[i] );
- free( (void *)mtx );
-}
-
-int getKouho( int *kouho, int nkouho, double *soukan, int nlen2 )
-{
- int i, j;
- int nlen4 = nlen2 / 2;
- double max;
- double tmp;
- int ikouho = 0; // by D.Mathog, iinoka?
- for( j=0; j<nkouho; j++ )
- {
- max = -9999.9;
- for( i=0; i<nlen2; i++ )
- {
- if( ( tmp = soukan[i] ) > max )
- {
- ikouho = i;
- max = tmp;
- }
- }
-#if 0
- if( max < 0.15 )
- {
- break;
- }
-#endif
-#if 0
- fprintf( stderr, "Kouho No.%d, pos=%d, score=%f, lag=%d\n", j, ikouho, soukan[ikouho], ikouho-nlen4 );
-#endif
- soukan[ikouho] = -9999.9;
- kouho[j] = ( ikouho - nlen4 );
- }
- return( j );
-}
-
-void zurasu2( int lag, int clus1, int clus2,
- char **seq1, char **seq2,
- char **aseq1, char **aseq2 )
-{
- int i;
-#if 0
- fprintf( stderr, "### lag = %d\n", lag );
-#endif
- if( lag > 0 )
- {
- for( i=0; i<clus1; i++ ) aseq1[i] = seq1[i];
- for( i=0; i<clus2; i++ ) aseq2[i] = seq2[i]+lag;
- }
- else
- {
- for( i=0; i<clus1; i++ ) aseq1[i] = seq1[i]-lag;
- for( i=0; i<clus2; i++ ) aseq2[i] = seq2[i];
- }
-}
-
-void zurasu( int lag, int clus1, int clus2,
- char **seq1, char **seq2,
- char **aseq1, char **aseq2 )
-{
- int i;
-#if DEBUG
- fprintf( stderr, "lag = %d\n", lag );
-#endif
- if( lag > 0 )
- {
- for( i=0; i<clus1; i++ ) strcpy( aseq1[i], seq1[i] );
- for( i=0; i<clus2; i++ ) strcpy( aseq2[i], seq2[i]+lag );
- }
- else
- {
- for( i=0; i<clus1; i++ ) strcpy( aseq1[i], seq1[i]-lag );
- for( i=0; i<clus2; i++ ) strcpy( aseq2[i], seq2[i] );
- }
-}
-
-
-int alignableReagion( int clus1, int clus2,
- char **seq1, char **seq2,
- double *eff1, double *eff2,
- Segment *seg )
-{
- int i, j, k;
- int status, starttmp = 0; // by D.Mathog, a gess
- double score;
- int value = 0;
- int len, maxlen;
- int length = 0; // by D.Mathog, a gess
- static TLS double *stra = NULL;
- static TLS int alloclen = 0;
- double totaleff;
- double cumscore;
- static TLS double threshold;
- static TLS double *prf1 = NULL;
- static TLS double *prf2 = NULL;
- static TLS int *hat1 = NULL;
- static TLS int *hat2 = NULL;
- int pre1, pre2;
-#if 0
- char **seq1pt;
- char **seq2pt;
- double *eff1pt;
- double *eff2pt;
-#endif
-
-#if 0
- fprintf( stderr, "### In alignableRegion, clus1=%d, clus2=%d \n", clus1, clus2 );
- fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
- fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
- fprintf( stderr, "eff1[0] = %f\n", eff1[0] );
- fprintf( stderr, "eff2[0] = %f\n", eff2[0] );
-#endif
-
- if( clus1 == 0 )
- {
- FreeDoubleVec( stra ); stra = NULL;
- FreeDoubleVec( prf1 ); prf1 = NULL;
- FreeDoubleVec( prf2 ); prf2 = NULL;
- FreeIntVec( hat1 ); hat1 = NULL;
- FreeIntVec( hat2 ); hat2 = NULL;
- return( 0 );
- }
-
- if( prf1 == NULL )
- {
- prf1 = AllocateDoubleVec( 26 );
- prf2 = AllocateDoubleVec( 26 );
- hat1 = AllocateIntVec( 27 );
- hat2 = AllocateIntVec( 27 );
- }
-
- len = MIN( strlen( seq1[0] ), strlen( seq2[0] ) );
- maxlen = MAX( strlen( seq1[0] ), strlen( seq2[0] ) ) + fftWinSize;
- if( alloclen < maxlen )
- {
- if( alloclen )
- {
- FreeDoubleVec( stra );
- }
- else
- {
- threshold = (int)fftThreshold / 100.0 * 600.0 * fftWinSize;
- }
- stra = AllocateDoubleVec( maxlen );
- alloclen = maxlen;
- }
-
-
- totaleff = 0.0;
- for( i=0; i<clus1; i++ ) for( j=0; j<clus2; j++ ) totaleff += eff1[i] * eff2[j];
- for( i=0; i<len; i++ )
- {
- /* make prfs */
- for( j=0; j<26; j++ )
- {
- prf1[j] = 0.0;
- prf2[j] = 0.0;
- }
-#if 0
- seq1pt = seq1;
- eff1pt = eff1;
- j = clus1;
- while( j-- ) prf1[amino_n[(*seq1pt++)[i]]] += *eff1pt++;
-#else
- for( j=0; j<clus1; j++ ) prf1[amino_n[(int)seq1[j][i]]] += eff1[j];
-#endif
- for( j=0; j<clus2; j++ ) prf2[amino_n[(int)seq2[j][i]]] += eff2[j];
-
- /* make hats */
- pre1 = pre2 = 26;
- for( j=25; j>=0; j-- )
- {
- if( prf1[j] )
- {
- hat1[pre1] = j;
- pre1 = j;
- }
- if( prf2[j] )
- {
- hat2[pre2] = j;
- pre2 = j;
- }
- }
- hat1[pre1] = -1;
- hat2[pre2] = -1;
-
- /* make site score */
- stra[i] = 0.0;
- for( k=hat1[26]; k!=-1; k=hat1[k] )
- for( j=hat2[26]; j!=-1; j=hat2[j] )
-// stra[i] += n_dis[k][j] * prf1[k] * prf2[j];
- stra[i] += n_disFFT[k][j] * prf1[k] * prf2[j];
- stra[i] /= totaleff;
- }
-
- (seg+0)->skipForeward = 0;
- (seg+1)->skipBackward = 0;
- status = 0;
- cumscore = 0.0;
- score = 0.0;
- for( j=0; j<fftWinSize; j++ ) score += stra[j];
-
- for( i=1; i<len-fftWinSize; i++ )
- {
- score = score - stra[i-1] + stra[i+fftWinSize-1];
-#if TMPTMPTMP
- fprintf( stderr, "%d %10.0f ? %10.0f\n", i, score, threshold );
-#endif
-
- if( score > threshold )
- {
-#if 0
- seg->start = i;
- seg->end = i;
- seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
- seg->score = score;
- status = 0;
- value++;
-#else
- if( !status )
- {
- status = 1;
- starttmp = i;
- length = 0;
- cumscore = 0.0;
- }
- length++;
- cumscore += score;
-#endif
- }
- if( score <= threshold || length > SEGMENTSIZE )
- {
- if( status )
- {
- if( length > fftWinSize )
- {
- seg->start = starttmp;
- seg->end = i;
- seg->center = ( seg->start + seg->end + fftWinSize ) / 2 ;
- seg->score = cumscore;
-#if 0
- fprintf( stderr, "%d-%d length = %d, score = %f, value = %d\n", seg->start, seg->end, length, cumscore, value );
-#endif
- if( length > SEGMENTSIZE )
- {
- (seg+0)->skipForeward = 1;
- (seg+1)->skipBackward = 1;
- }
- else
- {
- (seg+0)->skipForeward = 0;
- (seg+1)->skipBackward = 0;
- }
- value++;
- seg++;
- }
- length = 0;
- cumscore = 0.0;
- status = 0;
- starttmp = i;
- if( value > MAXSEG - 3 ) ErrorExit( "TOO MANY SEGMENTS!");
- }
- }
- }
- if( status && length > fftWinSize )
- {
- seg->end = i;
- seg->start = starttmp;
- seg->center = ( starttmp + i + fftWinSize ) / 2 ;
- seg->score = cumscore;
-#if 0
-fprintf( stderr, "%d-%d length = %d\n", seg->start, seg->end, length );
-#endif
- value++;
- }
-#if TMPTMPTMP
- exit( 0 );
-#endif
-// fprintf( stderr, "returning %d\n", value );
- return( value );
-}
-
-
-static int permit( Segment *seg1, Segment *seg2 )
-{
- return( 0 );
- if( seg1->end >= seg2->start ) return( 0 );
- if( seg1->pair->end >= seg2->pair->start ) return( 0 );
- else return( 1 );
-}
-
-void blockAlign2( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
-{
- int i, j, k, shift, cur1, cur2, count, klim;
- static TLS int crossscoresize = 0;
- static TLS int *result1 = NULL;
- static TLS int *result2 = NULL;
- static TLS int *ocut1 = NULL;
- static TLS int *ocut2 = NULL;
- double maximum;
- static TLS double **crossscore = NULL;
- static TLS int **track = NULL;
- static TLS double maxj, maxi;
- static TLS int pointj, pointi;
-
- if( cut1 == NULL)
- {
- if( result1 )
- {
- free( result1 );
- free( result2 );
- free( ocut1 );
- free( ocut2 );
- FreeIntMtx( track );
- FreeDoubleMtx( crossscore );
- }
- return;
- }
-
- if( result1 == NULL )
- {
- result1 = AllocateIntVec( MAXSEG );
- result2 = AllocateIntVec( MAXSEG );
- ocut1 = AllocateIntVec( MAXSEG );
- ocut2 = AllocateIntVec( MAXSEG );
- }
-
- if( crossscoresize < *ncut+2 )
- {
- crossscoresize = *ncut+2;
- if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
- if( track ) FreeIntMtx( track );
- if( crossscore ) FreeDoubleMtx( crossscore );
- track = AllocateIntMtx( crossscoresize, crossscoresize );
- crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
- }
-
-#if 0
- for( i=0; i<*ncut-2; i++ )
- fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
-
- for( i=0; i<*ncut; i++ )
- fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
- for( i=0; i<*ncut; i++ )
- {
- for( j=0; j<*ncut; j++ )
- fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
- fprintf( stderr, "\n" );
- }
-#endif
-
- for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ ) /* mudadanaa */
- crossscore[i][j] = ocrossscore[i][j];
- for( i=0; i<*ncut; i++ )
- {
- ocut1[i] = cut1[i];
- ocut2[i] = cut2[i];
- }
-
- for( i=1; i<*ncut; i++ )
- {
-#if 0
- fprintf( stderr, "### i=%d/%d\n", i,*ncut );
-#endif
- for( j=1; j<*ncut; j++ )
- {
- pointi = 0; maxi = 0.0;
- klim = j-2;
- for( k=0; k<klim; k++ )
- {
-/*
- fprintf( stderr, "k=%d, i=%d\n", k, i );
-*/
- if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
- if( crossscore[i-1][k] > maxj )
- {
- pointi = k;
- maxi = crossscore[i-1][k];
- }
- }
-
- pointj = 0; maxj = 0.0;
- klim = i-2;
- for( k=0; k<klim; k++ )
- {
- if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
- if( crossscore[k][j-1] > maxj )
- {
- pointj = k;
- maxj = crossscore[k][j-1];
- }
- }
-
- maxi += penalty;
- maxj += penalty;
-
- maximum = crossscore[i-1][j-1];
- track[i][j] = 0;
-
- if( maximum < maxi )
- {
- maximum = maxi ;
- track[i][j] = j - pointi;
- }
-
- if( maximum < maxj )
- {
- maximum = maxj ;
- track[i][j] = pointj - i;
- }
-
- crossscore[i][j] += maximum;
- }
- }
-#if 0
- for( i=0; i<*ncut; i++ )
- {
- for( j=0; j<*ncut; j++ )
- fprintf( stderr, "%3d ", track[i][j] );
- fprintf( stderr, "\n" );
- }
-#endif
-
-
- result1[MAXSEG-1] = *ncut-1;
- result2[MAXSEG-1] = *ncut-1;
-
- for( i=MAXSEG-1; i>=1; i-- )
- {
- cur1 = result1[i];
- cur2 = result2[i];
- if( cur1 == 0 || cur2 == 0 ) break;
- shift = track[cur1][cur2];
- if( shift == 0 )
- {
- result1[i-1] = cur1 - 1;
- result2[i-1] = cur2 - 1;
- continue;
- }
- else if( shift > 0 )
- {
- result1[i-1] = cur1 - 1;
- result2[i-1] = cur2 - shift;
- }
- else if( shift < 0 )
- {
- result1[i-1] = cur1 + shift;
- result2[i-1] = cur2 - 1;
- }
- }
-
- count = 0;
- for( j=i; j<MAXSEG; j++ )
- {
- if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
-
- if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
- if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
- count--;
-
- cut1[count] = ocut1[result1[j]];
- cut2[count] = ocut2[result2[j]];
-
- count++;
- }
-
- *ncut = count;
-#if 0
- for( i=0; i<*ncut; i++ )
- fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
-#endif
-}
-
-void blockAlign3( int *cut1, int *cut2, Segment **seg1, Segment **seg2, double **ocrossscore, int *ncut )
-// memory complexity = O(n^3), time complexity = O(n^2)
-{
- int i, j, shift, cur1, cur2, count;
- static TLS int crossscoresize = 0;
- static TLS int jumpposi, *jumppos;
- static TLS double jumpscorei, *jumpscore;
- static TLS int *result1 = NULL;
- static TLS int *result2 = NULL;
- static TLS int *ocut1 = NULL;
- static TLS int *ocut2 = NULL;
- double maximum;
- static TLS double **crossscore = NULL;
- static TLS int **track = NULL;
-
- if( result1 == NULL )
- {
- result1 = AllocateIntVec( MAXSEG );
- result2 = AllocateIntVec( MAXSEG );
- ocut1 = AllocateIntVec( MAXSEG );
- ocut2 = AllocateIntVec( MAXSEG );
- }
- if( crossscoresize < *ncut+2 )
- {
- crossscoresize = *ncut+2;
- if( fftkeika ) fprintf( stderr, "allocating crossscore and track, size = %d\n", crossscoresize );
- if( track ) FreeIntMtx( track );
- if( crossscore ) FreeDoubleMtx( crossscore );
- if( jumppos ) FreeIntVec( jumppos );
- if( jumpscore ) FreeDoubleVec( jumpscore );
- track = AllocateIntMtx( crossscoresize, crossscoresize );
- crossscore = AllocateDoubleMtx( crossscoresize, crossscoresize );
- jumppos = AllocateIntVec( crossscoresize );
- jumpscore = AllocateDoubleVec( crossscoresize );
- }
-
-#if 0
- for( i=0; i<*ncut-2; i++ )
- fprintf( stderr, "%d.start = %d, score = %f\n", i, seg1[i]->start, seg1[i]->score );
-
- for( i=0; i<*ncut; i++ )
- fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
- for( i=0; i<*ncut; i++ )
- {
- for( j=0; j<*ncut; j++ )
- fprintf( stderr, "%#4.0f ", ocrossscore[i][j] );
- fprintf( stderr, "\n" );
- }
-#endif
-
- for( i=0; i<*ncut; i++ ) for( j=0; j<*ncut; j++ ) /* mudadanaa */
- crossscore[i][j] = ocrossscore[i][j];
- for( i=0; i<*ncut; i++ )
- {
- ocut1[i] = cut1[i];
- ocut2[i] = cut2[i];
- }
- for( j=0; j<*ncut; j++ )
- {
- jumpscore[j] = -999.999;
- jumppos[j] = -1;
- }
-
- for( i=1; i<*ncut; i++ )
- {
-
- jumpscorei = -999.999;
- jumpposi = -1;
-
- for( j=1; j<*ncut; j++ )
- {
-#if 1
- fprintf( stderr, "in blockalign3, ### i=%d, j=%d\n", i, j );
-#endif
-
-
-#if 0
- for( k=0; k<j-2; k++ )
- {
-/*
- fprintf( stderr, "k=%d, i=%d\n", k, i );
-*/
- if( k && k<*ncut-1 && j<*ncut-1 && !permit( seg1[k-1], seg1[j-1] ) ) continue;
- if( crossscore[i-1][k] > maxj )
- {
- pointi = k;
- maxi = crossscore[i-1][k];
- }
- }
-
- pointj = 0; maxj = 0.0;
- for( k=0; k<i-2; k++ )
- {
- if( k && k<*ncut-1 && i<*ncut-1 && !permit( seg2[k-1], seg2[i-1] ) ) continue;
- if( crossscore[k][j-1] > maxj )
- {
- pointj = k;
- maxj = crossscore[k][j-1];
- }
- }
-
-
- maxi += penalty;
- maxj += penalty;
-#endif
- maximum = crossscore[i-1][j-1];
- track[i][j] = 0;
-
- if( maximum < jumpscorei && permit( seg1[jumpposi], seg1[i] ) )
- {
- maximum = jumpscorei;
- track[i][j] = j - jumpposi;
- }
-
- if( maximum < jumpscore[j] && permit( seg2[jumppos[j]], seg2[j] ) )
- {
- maximum = jumpscore[j];
- track[i][j] = jumpscore[j] - i;
- }
-
- crossscore[i][j] += maximum;
-
- if( jumpscorei < crossscore[i-1][j] )
- {
- jumpscorei = crossscore[i-1][j];
- jumpposi = j;
- }
-
- if( jumpscore[j] < crossscore[i][j-1] )
- {
- jumpscore[j] = crossscore[i][j-1];
- jumppos[j] = i;
- }
- }
- }
-#if 0
- for( i=0; i<*ncut; i++ )
- {
- for( j=0; j<*ncut; j++ )
- fprintf( stderr, "%3d ", track[i][j] );
- fprintf( stderr, "\n" );
- }
-#endif
-
-
- result1[MAXSEG-1] = *ncut-1;
- result2[MAXSEG-1] = *ncut-1;
-
- for( i=MAXSEG-1; i>=1; i-- )
- {
- cur1 = result1[i];
- cur2 = result2[i];
- if( cur1 == 0 || cur2 == 0 ) break;
- shift = track[cur1][cur2];
- if( shift == 0 )
- {
- result1[i-1] = cur1 - 1;
- result2[i-1] = cur2 - 1;
- continue;
- }
- else if( shift > 0 )
- {
- result1[i-1] = cur1 - 1;
- result2[i-1] = cur2 - shift;
- }
- else if( shift < 0 )
- {
- result1[i-1] = cur1 + shift;
- result2[i-1] = cur2 - 1;
- }
- }
-
- count = 0;
- for( j=i; j<MAXSEG; j++ )
- {
- if( ocrossscore[result1[j]][result2[j]] == 0.0 ) continue;
-
- if( result1[j] == result1[j-1] || result2[j] == result2[j-1] )
- if( ocrossscore[result1[j]][result2[j]] > ocrossscore[result1[j-1]][result2[j-1]] )
- count--;
-
- cut1[count] = ocut1[result1[j]];
- cut2[count] = ocut2[result2[j]];
-
- count++;
- }
-
- *ncut = count;
-#if 0
- for( i=0; i<*ncut; i++ )
- fprintf( stderr, "i=%d, cut1 = %d, cut2 = %d\n", i, cut1[i], cut2[i] );
-#endif
-}
-
git://source.jalview.org / jabaws.git / blobdiff