[llvm-commits] [test-suite] r65043 [2/2] - in /test-suite/trunk/MultiSource/Benchmarks: ./ mafft/
Torok Edwin
edwintorok at gmail.com
Thu Feb 19 04:38:03 PST 2009
Added: test-suite/trunk/MultiSource/Benchmarks/mafft/rna.c
URL: http://llvm.org/viewvc/llvm-project/test-suite/trunk/MultiSource/Benchmarks/mafft/rna.c?rev=65043&view=auto
==============================================================================
--- test-suite/trunk/MultiSource/Benchmarks/mafft/rna.c (added)
+++ test-suite/trunk/MultiSource/Benchmarks/mafft/rna.c Thu Feb 19 06:37:59 2009
@@ -0,0 +1,560 @@
+#include "mltaln.h"
+#include "dp.h"
+
+#define MEMSAVE 1
+
+#define DEBUG 1
+#define USE_PENALTY_EX 1
+#define STOREWM 1
+
+
+
+static float singleribosumscore( int n1, int n2, char **s1, char **s2, double *eff1, double *eff2, int p1, int p2 )
+{
+ float val;
+ int i, j;
+ int code1, code2;
+
+ val = 0.0;
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
+ {
+ code1 = amino_n[(int)s1[i][p1]];
+ if( code1 > 3 ) code1 = 36;
+ code2 = amino_n[(int)s2[j][p2]];
+ if( code2 > 3 ) code2 = 36;
+
+// fprintf( stderr, "'l'%c-%c: %f\n", s1[i][p1], s2[j][p2], (float)ribosumdis[code1][code2] );
+
+ val += (float)ribosumdis[code1][code2] * eff1[i] * eff2[j];
+ }
+ return( val );
+}
+static float pairedribosumscore53( int n1, int n2, char **s1, char **s2, double *eff1, double *eff2, int p1, int p2, int c1, int c2 )
+{
+ float val;
+ int i, j;
+ int code1o, code1u, code2o, code2u, code1, code2;
+
+ val = 0.0;
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
+ {
+ code1o = amino_n[(int)s1[i][p1]];
+ code1u = amino_n[(int)s1[i][c1]];
+ if( code1o > 3 ) code1 = code1o = 36;
+ else if( code1u > 3 ) code1 = 36;
+ else code1 = 4 + code1o * 4 + code1u;
+
+ code2o = amino_n[(int)s2[j][p2]];
+ code2u = amino_n[(int)s2[j][c2]];
+ if( code2o > 3 ) code2 = code1o = 36;
+ else if( code2u > 3 ) code2 = 36;
+ else code2 = 4 + code2o * 4 + code2u;
+
+
+// fprintf( stderr, "%c%c-%c%c: %f\n", s1[i][p1], s1[i][c1], s2[j][p2], s2[j][c2], (float)ribosumdis[code1][code2] );
+
+ if( code1 == 36 || code2 == 36 )
+ val += (float)n_dis[code1o][code2o] * eff1[i] * eff2[j];
+ else
+ val += (float)ribosumdis[code1][code2] * eff1[i] * eff2[j];
+ }
+ return( val );
+}
+
+static float pairedribosumscore35( int n1, int n2, char **s1, char **s2, double *eff1, double *eff2, int p1, int p2, int c1, int c2 )
+{
+ float val;
+ int i, j;
+ int code1o, code1u, code2o, code2u, code1, code2;
+
+ val = 0.0;
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
+ {
+ code1o = amino_n[(int)s1[i][p1]];
+ code1u = amino_n[(int)s1[i][c1]];
+ if( code1o > 3 ) code1 = code1o = 36;
+ else if( code1u > 3 ) code1 = 36;
+ else code1 = 4 + code1u * 4 + code1o;
+
+ code2o = amino_n[(int)s2[j][p2]];
+ code2u = amino_n[(int)s2[j][c2]];
+ if( code2o > 3 ) code2 = code1o = 36;
+ else if( code2u > 3 ) code2 = 36;
+ else code2 = 4 + code2u * 4 + code2o;
+
+
+// fprintf( stderr, "%c%c-%c%c: %f\n", s1[i][p1], s1[i][c1], s2[j][p2], s2[j][c2], (float)ribosumdis[code1][code2] );
+
+ if( code1 == 36 || code2 == 36 )
+ val += (float)n_dis[code1o][code2o] * eff1[i] * eff2[j];
+ else
+ val += (float)ribosumdis[code1][code2] * eff1[i] * eff2[j];
+ }
+ return( val );
+}
+
+
+static void mccaskillextract( char **seq, char **nogap, int nseq, RNApair **pairprob, RNApair ***single, int **sgapmap, double *eff )
+{
+ int lgth;
+ int nogaplgth;
+ int i, j;
+ int left, right, adpos;
+ float prob;
+ static int *pairnum;
+ RNApair *pt, *pt2;
+
+ lgth = strlen( seq[0] );
+ pairnum = calloc( lgth, sizeof( int ) );
+ for( i=0; i<lgth; i++ ) pairnum[i] = 0;
+
+ for( i=0; i<nseq; i++ )
+ {
+ nogaplgth = strlen( nogap[i] );
+ for( j=0; j<nogaplgth; j++ ) for( pt=single[i][j]; pt->bestpos!=-1; pt++ )
+ {
+ left = sgapmap[i][j];
+ right = sgapmap[i][pt->bestpos];
+ prob = pt->bestscore;
+
+
+ for( pt2=pairprob[left]; pt2->bestpos!=-1; pt2++ )
+ if( pt2->bestpos == right ) break;
+
+// fprintf( stderr, "i,j=%d,%d, left=%d, right=%d, pt=%d, pt2->bestpos = %d\n", i, j, left, right, pt-single[i][j], pt2->bestpos );
+ if( pt2->bestpos == -1 )
+ {
+ pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
+ adpos = pairnum[left];
+ pairnum[left]++;
+ pairprob[left][adpos].bestscore = 0.0;
+ pairprob[left][adpos].bestpos = right;
+ pairprob[left][adpos+1].bestscore = -1.0;
+ pairprob[left][adpos+1].bestpos = -1;
+ pt2 = pairprob[left]+adpos;
+ }
+ else
+ adpos = pt2-pairprob[left];
+
+ pt2->bestscore += prob * eff[i];
+
+ if( pt2->bestpos != right )
+ {
+ fprintf( stderr, "okashii!\n" );
+ exit( 1 );
+ }
+// fprintf( stderr, "adding %d-%d, %f\n", left, right, prob );
+// fprintf( stderr, "pairprob[0][0].bestpos=%d\n", pairprob[0][0].bestpos );
+// fprintf( stderr, "pairprob[0][0].bestscore=%f\n", pairprob[0][0].bestscore );
+ }
+ }
+
+// fprintf( stderr, "before taikakuka\n" );
+ for( i=0; i<lgth; i++ ) for( j=0; j<pairnum[i]; j++ )
+ {
+ if( pairprob[i][j].bestpos > -1 )
+ {
+// pairprob[i][j].bestscore /= (float)nseq;
+// fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, pairprob[i][j].bestpos, pairprob[i][j].bestscore, seq[0][i], seq[0][pairprob[i][j].bestpos] );
+ }
+ }
+
+#if 0
+ for( i=0; i<lgth; i++ ) for( j=0; j<pairnum[i]; j++ )
+ {
+ right=pairprob[i][j].bestpos;
+ if( right < i ) continue;
+ fprintf( stderr, "no%d-%d, adding %d -> %d\n", i, j, right, i );
+ pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
+ pairprob[right][pairnum[right]].bestscore = pairprob[i][j].bestscore;
+ pairprob[right][pairnum[right]].bestpos = i;
+ pairnum[right]++;
+ pairprob[right][pairnum[right]].bestscore = -1.0;
+ pairprob[right][pairnum[right]].bestpos = -1;
+ }
+#endif
+
+ free( pairnum );
+
+}
+
+
+void rnaalifoldcall( char **seq, int nseq, RNApair **pairprob )
+{
+ int lgth;
+ int i;
+ static int *order = NULL;
+ static char **name = NULL;
+ char gett[1000];
+ FILE *fp;
+ int left, right, dumm;
+ float prob;
+ static int pid;
+ static char fnamein[100];
+ static char cmd[1000];
+ static int *pairnum;
+
+ lgth = strlen( seq[0] );
+ if( order == NULL )
+ {
+ pid = (int)getpid();
+ sprintf( fnamein, "/tmp/_rnaalifoldin.%d", pid );
+ order = AllocateIntVec( njob );
+ name = AllocateCharMtx( njob, 10 );
+ for( i=0; i<njob; i++ )
+ {
+ order[i] = i;
+ sprintf( name[i], "seq%d", i );
+ }
+ }
+ pairnum = calloc( lgth, sizeof( int ) );
+ for( i=0; i<lgth; i++ ) pairnum[i] = 0;
+
+ fp = fopen( fnamein, "w" );
+ if( !fp )
+ {
+ fprintf( stderr, "Cannot open /tmp/_rnaalifoldin\n" );
+ exit( 1 );
+ }
+ clustalout_pointer( fp, nseq, lgth, seq, name, NULL, NULL, order );
+ fclose( fp );
+
+ sprintf( cmd, "RNAalifold -p %s", fnamein );
+ system( cmd );
+
+ fp = fopen( "alifold.out", "r" );
+ if( !fp )
+ {
+ fprintf( stderr, "Cannot open /tmp/_rnaalifoldin\n" );
+ exit( 1 );
+ }
+
+#if 0
+ for( i=0; i<lgth; i++ ) // atode kesu
+ {
+ pairprob[i] = (RNApair *)realloc( pairprob[i], (2) * sizeof( RNApair ) ); // atode kesu
+ pairprob[i][1].bestscore = -1.0;
+ pairprob[i][1].bestpos = -1;
+ }
+#endif
+
+ while( 1 )
+ {
+ fgets( gett, 999, fp );
+ if( gett[0] == '(' ) break;
+ if( gett[0] == '{' ) break;
+ if( gett[0] == '.' ) break;
+ if( gett[0] == ',' ) break;
+ if( gett[0] != ' ' ) continue;
+
+ sscanf( gett, "%d %d %d %f", &left, &right, &dumm, &prob );
+ left--;
+ right--;
+
+
+#if 0
+ if( prob > 50.0 && prob > pairprob[left][0].bestscore )
+ {
+ pairprob[left][0].bestscore = prob;
+ pairprob[left][0].bestpos = right;
+#else
+ if( prob > 0.0 )
+ {
+ pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
+ pairprob[left][pairnum[left]].bestscore = prob / 100.0;
+ pairprob[left][pairnum[left]].bestpos = right;
+ pairnum[left]++;
+ pairprob[left][pairnum[left]].bestscore = -1.0;
+ pairprob[left][pairnum[left]].bestpos = -1;
+ fprintf( stderr, "%d-%d, %f\n", left, right, prob );
+
+ pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
+ pairprob[right][pairnum[right]].bestscore = prob / 100.0;
+ pairprob[right][pairnum[right]].bestpos = left;
+ pairnum[right]++;
+ pairprob[right][pairnum[right]].bestscore = -1.0;
+ pairprob[right][pairnum[right]].bestpos = -1;
+ fprintf( stderr, "%d-%d, %f\n", left, right, prob );
+#endif
+ }
+ }
+ fclose( fp );
+ sprintf( cmd, "rm -f %s", fnamein );
+ system( cmd );
+
+ for( i=0; i<lgth; i++ )
+ {
+ if( (right=pairprob[i][0].bestpos) > -1 )
+ {
+ pairprob[right][0].bestpos = i;
+ pairprob[right][0].bestscore = pairprob[i][0].bestscore;
+ }
+ }
+
+#if 0
+ for( i=0; i<lgth; i++ ) // atode kesu
+ if( pairprob[i][0].bestscore > -1 ) pairprob[i][0].bestscore = 1.0; // atode kesu
+#endif
+
+// fprintf( stderr, "after taikakuka in rnaalifoldcall\n" );
+// for( i=0; i<lgth; i++ )
+// {
+// fprintf( stderr, "pair of %d = %d (%f) %c:%c\n", i, pairprob[i][0].bestpos, pairprob[i][0].bestscore, seq[0][i], seq[0][pairprob[i][0].bestpos] );
+// }
+
+ free( pairnum );
+}
+
+
+static void utot( int n, int l, char **s )
+{
+ int i, j;
+ for( i=0; i<l; i++ )
+ {
+ for( j=0; j<n; j++ )
+ {
+ if ( s[j][i] == 'a' ) s[j][i] = 'a';
+ else if( s[j][i] == 't' ) s[j][i] = 't';
+ else if( s[j][i] == 'u' ) s[j][i] = 't';
+ else if( s[j][i] == 'g' ) s[j][i] = 'g';
+ else if( s[j][i] == 'c' ) s[j][i] = 'c';
+ else if( s[j][i] == '-' ) s[j][i] = '-';
+ else s[j][i] = 'n';
+ }
+ }
+}
+
+
+void foldrna( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, float **impmtx, int *gapmap1, int *gapmap2, RNApair *additionalpair )
+{
+ int i, j;
+// int ui, uj;
+// int uiup, ujup;
+ int uido, ujdo;
+ static char **useq1, **useq2;
+ static char **oseq1, **oseq2, **oseq1r, **oseq2r, *odir1, *odir2;
+ static RNApair **pairprob1, **pairprob2;
+ static RNApair *pairpt1, *pairpt2;
+ int lgth1 = strlen( seq1[0] );
+ int lgth2 = strlen( seq2[0] );
+ static float **impmtx2;
+ static float **map;
+// double lenfac;
+ float prob;
+ int **sgapmap1, **sgapmap2;
+ char *nogapdum;
+ float **tbppmtx;
+
+
+// fprintf( stderr, "nseq1=%d, lgth1=%d\n", nseq1, lgth1 );
+ useq1 = AllocateCharMtx( nseq1, lgth1+10 );
+ useq2 = AllocateCharMtx( nseq2, lgth2+10 );
+ oseq1 = AllocateCharMtx( nseq1, lgth1+10 );
+ oseq2 = AllocateCharMtx( nseq2, lgth2+10 );
+ oseq1r = AllocateCharMtx( nseq1, lgth1+10 );
+ oseq2r = AllocateCharMtx( nseq2, lgth2+10 );
+ odir1 = AllocateCharVec( lgth1+10 );
+ odir2 = AllocateCharVec( lgth2+10 );
+ sgapmap1 = AllocateIntMtx( nseq1, lgth1+1 );
+ sgapmap2 = AllocateIntMtx( nseq2, lgth2+1 );
+ nogapdum = AllocateCharVec( MAX( lgth1, lgth2 ) );
+ pairprob1 = (RNApair **)calloc( lgth1, sizeof( RNApair *) );
+ pairprob2 = (RNApair **)calloc( lgth2, sizeof( RNApair *) );
+ map = AllocateFloatMtx( lgth1, lgth2 );
+ impmtx2 = AllocateFloatMtx( lgth1, lgth2 );
+ tbppmtx = AllocateFloatMtx( lgth1, lgth2 );
+
+ for( i=0; i<nseq1; i++ ) strcpy( useq1[i], seq1[i] );
+ for( i=0; i<nseq2; i++ ) strcpy( useq2[i], seq2[i] );
+ for( i=0; i<nseq1; i++ ) strcpy( oseq1[i], seq1[i] );
+ for( i=0; i<nseq2; i++ ) strcpy( oseq2[i], seq2[i] );
+
+ for( i=0; i<nseq1; i++ ) commongappick_record( 1, useq1+i, sgapmap1[i] );
+ for( i=0; i<nseq2; i++ ) commongappick_record( 1, useq2+i, sgapmap2[i] );
+
+ for( i=0; i<lgth1; i++ )
+ {
+ pairprob1[i] = (RNApair *)calloc( 1, sizeof( RNApair ) );
+ pairprob1[i][0].bestpos = -1;
+ pairprob1[i][0].bestscore = -1;
+ }
+ for( i=0; i<lgth2; i++ )
+ {
+ pairprob2[i] = (RNApair *)calloc( 1, sizeof( RNApair ) );
+ pairprob2[i][0].bestpos = -1;
+ pairprob2[i][0].bestscore = -1;
+ }
+
+ utot( nseq1, lgth1, oseq1 );
+ utot( nseq2, lgth2, oseq2 );
+
+// fprintf( stderr, "folding group1\n" );
+// rnalocal( oseq1, useq1, eff1, eff1, nseq1, nseq1, lgth1+10, pair1 );
+
+/* base-pairing probability of group 1 */
+ if( rnaprediction == 'r' )
+ rnaalifoldcall( oseq1, nseq1, pairprob1 );
+ else
+ mccaskillextract( oseq1, useq1, nseq1, pairprob1, grouprna1, sgapmap1, eff1 );
+
+
+// fprintf( stderr, "folding group2\n" );
+// rnalocal( oseq2, useq2, eff2, eff2, nseq2, nseq2, lgth2+10, pair2 );
+
+/* base-pairing probability of group 2 */
+ if( rnaprediction == 'r' )
+ rnaalifoldcall( oseq2, nseq2, pairprob2 );
+ else
+ mccaskillextract( oseq2, useq2, nseq2, pairprob2, grouprna2, sgapmap2, eff2 );
+
+
+
+#if 0
+ makerseq( oseq1, oseq1r, odir1, pairprob1, nseq1, lgth1 );
+ makerseq( oseq2, oseq2r, odir2, pairprob2, nseq2, lgth2 );
+
+ fprintf( stderr, "%s\n", odir2 );
+
+ for( i=0; i<nseq1; i++ )
+ {
+ fprintf( stdout, ">ori\n%s\n", oseq1[0] );
+ fprintf( stdout, ">rev\n%s\n", oseq1r[0] );
+ }
+#endif
+
+/* similarity score */
+ Lalignmm_hmout( oseq1, oseq2, eff1, eff2, nseq1, nseq2, 10000, NULL, NULL, NULL, NULL, map );
+
+ if( 1 )
+ {
+ if( RNAscoremtx == 'n' )
+ {
+ for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
+ {
+// impmtx2[i][j] = osoiaveragescore( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j ) * consweight_multi;
+ impmtx2[i][j] = 0.0;
+ }
+ }
+ else if( RNAscoremtx == 'r' )
+ {
+ for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
+ {
+ tbppmtx[i][j] = 1.0;
+ impmtx2[i][j] = 0.0;
+ }
+ for( i=0; i<lgth1; i++ ) for( pairpt1=pairprob1[i]; pairpt1->bestpos!=-1; pairpt1++ )
+ {
+ for( j=0; j<lgth2; j++ ) for( pairpt2=pairprob2[j]; pairpt2->bestpos!=-1; pairpt2++ )
+ {
+ uido = pairpt1->bestpos;
+ ujdo = pairpt2->bestpos;
+ prob = pairpt1->bestscore * pairpt2->bestscore;
+ if( uido > -1 && ujdo > -1 )
+ {
+ if( uido > i && j > ujdo )
+ {
+ impmtx2[i][j] += prob * pairedribosumscore53( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j, uido, ujdo ) * consweight_multi;
+ tbppmtx[i][j] -= prob;
+ }
+ else if( i < uido && j < ujdo )
+ {
+ impmtx2[i][j] += prob * pairedribosumscore35( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j, uido, ujdo ) * consweight_multi;
+ tbppmtx[i][j] -= prob;
+ }
+ }
+ }
+ }
+
+
+ for( i=0; i<lgth1; i++ )
+ {
+ for( j=0; j<lgth2; j++ )
+ {
+ impmtx2[i][j] += tbppmtx[i][j] * singleribosumscore( nseq1, nseq2, oseq1, oseq2, eff1, eff2, i, j ) * consweight_multi;
+ }
+ }
+ }
+
+
+/* four-way consistency */
+
+ for( i=0; i<lgth1; i++ ) for( pairpt1=pairprob1[i]; pairpt1->bestpos!=-1; pairpt1++ )
+ {
+
+// if( pairprob1[i] == NULL ) continue;
+
+ for( j=0; j<lgth2; j++ ) for( pairpt2=pairprob2[j]; pairpt2->bestpos!=-1; pairpt2++ )
+ {
+// fprintf( stderr, "i=%d, j=%d, pn1=%d, pn2=%d\n", i, j, pairpt1-pairprob1[i], pairpt2-pairprob2[j] );
+// if( pairprob2[j] == NULL ) continue;
+
+ uido = pairpt1->bestpos;
+ ujdo = pairpt2->bestpos;
+ prob = pairpt1->bestscore * pairpt2->bestscore;
+// prob = 1.0;
+// fprintf( stderr, "i=%d->uido=%d, j=%d->ujdo=%d\n", i, uido, j, ujdo );
+
+// fprintf( stderr, "impmtx2[%d][%d] = %f\n", i, j, impmtx2[i][j] );
+
+// if( i < uido && j > ujdo ) continue;
+// if( i > uido && j < ujdo ) continue;
+
+
+// posdistj = abs( ujdo-j );
+
+// if( uido > -1 && ujdo > -1 )
+ if( uido > -1 && ujdo > -1 && ( ( i > uido && j > ujdo ) || ( i < uido && j < ujdo ) ) )
+ {
+ {
+ impmtx2[i][j] += MAX( 0, map[uido][ujdo] ) * consweight_rna * 600 * prob; // osoi
+ }
+ }
+
+ }
+ }
+ for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
+ {
+ impmtx[i][j] += impmtx2[i][j];
+// fprintf( stderr, "fastathreshold=%f, consweight_multi=%f, consweight_rna=%f\n", fastathreshold, consweight_multi, consweight_rna );
+// impmtx[i][j] *= 0.5;
+ }
+
+// impmtx[0][0] += 10000.0;
+// impmtx[lgth1-1][lgth2-1] += 10000.0;
+
+
+
+#if 0
+ fprintf( stdout, "#impmtx2 = \n" );
+ for( i=0; i<lgth1; i++ )
+ {
+ for( j=0; j<lgth2; j++ )
+ {
+ fprintf( stdout, "%d %d %f\n", i, j, impmtx2[i][j] );
+ }
+ fprintf( stdout, "\n" );
+ }
+ exit( 1 );
+#endif
+ }
+
+ FreeCharMtx( useq1 );
+ FreeCharMtx( useq2 );
+ FreeCharMtx( oseq1 );
+ FreeCharMtx( oseq2 );
+ FreeCharMtx( oseq1r );
+ FreeCharMtx( oseq2r );
+ free( odir1 );
+ free( odir2 );
+ FreeFloatMtx( impmtx2 );
+ FreeFloatMtx( map );
+ FreeIntMtx( sgapmap1 );
+ FreeIntMtx( sgapmap2 );
+ FreeFloatMtx( tbppmtx );
+
+ for( i=0; i<lgth1; i++ ) free( pairprob1[i] );
+ for( i=0; i<lgth2; i++ ) free( pairprob2[i] );
+ free( pairprob1 );
+ free( pairprob2 );
+}
+
Added: test-suite/trunk/MultiSource/Benchmarks/mafft/share.h
URL: http://llvm.org/viewvc/llvm-project/test-suite/trunk/MultiSource/Benchmarks/mafft/share.h?rev=65043&view=auto
==============================================================================
--- test-suite/trunk/MultiSource/Benchmarks/mafft/share.h (added)
+++ test-suite/trunk/MultiSource/Benchmarks/mafft/share.h Thu Feb 19 06:37:59 2009
@@ -0,0 +1,16 @@
+#if 0
+#include <sys/types.h>
+#include <sys/ipc.h>
+#include <sys/shm.h>
+
+#endif
+#define IMA_YONDERU 'x' /* iranai */
+#define IMA_KAITERU 0 /* iranai */
+#define KAKIOWATTA 'w'
+#define YOMIOWATTA 'r'
+#define OSHIMAI 'd'
+#define ISRUNNING 0
+#define SEMAPHORE 1
+#define STATUS 2
+
+#define IPC_ALLOC 0100000
Added: test-suite/trunk/MultiSource/Benchmarks/mafft/suboptalign11.c
URL: http://llvm.org/viewvc/llvm-project/test-suite/trunk/MultiSource/Benchmarks/mafft/suboptalign11.c?rev=65043&view=auto
==============================================================================
--- test-suite/trunk/MultiSource/Benchmarks/mafft/suboptalign11.c (added)
+++ test-suite/trunk/MultiSource/Benchmarks/mafft/suboptalign11.c Thu Feb 19 06:37:59 2009
@@ -0,0 +1,676 @@
+#include "mltaln.h"
+#include "dp.h"
+
+#define DEBUG 0
+#define DEBUG2 0
+#define XXXXXXX 0
+#define USE_PENALTY_EX 1
+
+typedef struct _shuryoten
+{
+ int i;
+ int j;
+ float wm;
+ struct _shuryoten *next;
+ struct _shuryoten *prev;
+} Shuryoten;
+
+
+static int localstop;
+
+static int compshuryo( Shuryoten *s1_arg, Shuryoten *s2_arg )
+{
+ Shuryoten *s1 = (Shuryoten *)s1_arg;
+ Shuryoten *s2 = (Shuryoten *)s2_arg;
+ if ( s1->wm > s2->wm ) return( -1 );
+ else if ( s1->wm < s2->wm ) return( 1 );
+ else return( 0 );
+}
+
+static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
+{
+ int j;
+
+ for( j=0; j<lgth2; j++ )
+ match[j] = amino_dis[(int)(*s1)[i1]][(int)(*s2)[j]];
+}
+
+static float gentracking( int **used,
+ char **seq1, char **seq2,
+ char **mseq1, char **mseq2,
+ float **cpmx1, float **cpmx2,
+ int **ijpi, int **ijpj, int *off1pt, int *off2pt, int endi, int endj )
+{
+ int l, iin, jin, lgth1, lgth2, k, limk;
+ int ifi=0, jfi=0;
+ char gap[] = "-";
+ static char *res1 = NULL, *res2 = NULL;
+ char *mspt1, *mspt2;
+ if( res1 == NULL )
+ {
+ res1 = (char *)calloc( N, sizeof( char ) );
+ res2 = (char *)calloc( N, sizeof( char ) );
+ }
+
+ lgth1 = strlen( seq1[0] );
+ lgth2 = strlen( seq2[0] );
+
+ mspt1 = res1 + lgth1+lgth2;
+ *mspt1 = 0;
+ mspt2 = res2 + lgth1+lgth2;
+ *mspt2 = 0;
+ iin = endi; jin = endj;
+
+ limk = lgth1+lgth2;
+ if( used[iin][jin] ) return( -1.0 );
+ for( k=0; k<=limk; k++ )
+ {
+ ifi = ( ijpi[iin][jin] );
+ jfi = ( ijpj[iin][jin] );
+
+ if( used[ifi][jfi] ) return( -1.0 );
+
+ l = iin - ifi;
+ while( --l )
+ {
+ *--mspt1 = seq1[0][ifi+l];
+ *--mspt2 = *gap;
+ k++;
+ }
+ l= jin - jfi;
+ while( --l )
+ {
+ *--mspt1 = *gap;
+ *--mspt2 = seq2[0][jfi+l];
+ k++;
+ }
+
+ if( iin <= 0 || jin <= 0 ) break;
+ *--mspt1 = seq1[0][ifi];
+ *--mspt2 = seq2[0][jfi];
+ if( ijpi[ifi][jfi] == localstop ) break;
+ if( ijpj[ifi][jfi] == localstop ) break;
+ k++;
+ iin = ifi; jin = jfi;
+ }
+ if( ifi == -1 ) *off1pt = 0; else *off1pt = ifi;
+ if( jfi == -1 ) *off2pt = 0; else *off2pt = jfi;
+
+// fprintf( stderr, "ifn = %d, jfn = %d\n", ifi, jfi );
+
+ iin = endi; jin = endj;
+ limk = lgth1+lgth2;
+ for( k=0; k<=limk; k++ )
+ {
+ ifi = ( ijpi[iin][jin] );
+ jfi = ( ijpj[iin][jin] );
+
+ used[ifi][jfi] = 1;
+ if( iin <= 0 || jin <= 0 ) break;
+ if( ijpi[ifi][jfi] == localstop ) break;
+ if( ijpj[ifi][jfi] == localstop ) break;
+
+ k++;
+ iin = ifi; jin = jfi;
+ }
+
+
+ strcpy( mseq1[0], mspt1 );
+ strcpy( mseq2[0], mspt2 );
+
+ fprintf( stderr, "mseq1=%s\nmseq2=%s\n", mspt1, mspt2 );
+
+ return( 0.0 );
+}
+
+
+float suboptalign11( char **seq1, char **seq2, int alloclen, int *off1pt, int *off2pt, LocalHom *lhmpt )
+/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
+{
+ int k;
+ static int **used;
+ register int i, j;
+ int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
+ int lgth1, lgth2;
+ int resultlen;
+ float wm = 0.0; // by D.Mathog,
+ float g;
+ float *currentw, *previousw;
+#if 1
+ float *wtmp;
+ int *ijpipt;
+ int *ijpjpt;
+ float *mjpt, *Mjpt, *prept, *curpt;
+ int *mpjpt, *Mpjpt;
+#endif
+ static float mi, *m;
+ static float Mi, *largeM;
+ static int **ijpi;
+ static int **ijpj;
+ static int mpi, *mp;
+ 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 float **cpmx1;
+ static float **cpmx2;
+ static int **intwork;
+ static float **floatwork;
+ static int orlgth1 = 0, orlgth2 = 0;
+ float maxwm;
+ float tbk;
+ int tbki, tbkj;
+ int endali, endalj;
+// float localthr = 0.0;
+// float localthr2 = 0.0;
+ float fpenalty = (float)penalty;
+ float fpenalty_OP = (float)penalty_OP;
+ float fpenalty_ex = (float)penalty_ex;
+// float fpenalty_EX = (float)penalty_EX;
+ float foffset = (float)offset;
+ float localthr = -foffset;
+ float localthr2 = -foffset;
+ static Shuryoten *shuryo = NULL;
+ int numshuryo;
+ float minshuryowm = 0.0; // by D.Mathog
+ int minshuryopos = 0; // by D.Mathog
+ float resf;
+
+
+// fprintf( stderr, "@@@@@@@@@@@@@ penalty_OP = %f, penalty_EX = %f, pelanty = %f\n", fpenalty_OP, fpenalty_EX, fpenalty );
+
+ fprintf( stderr, "in suboptalign11\n" );
+ if( !shuryo )
+ {
+ shuryo = (Shuryoten *)calloc( 100, sizeof( Shuryoten ) );
+ }
+ for( i=0; i<100; i++ )
+ {
+ shuryo[i].i = -1;
+ shuryo[i].j = -1;
+ shuryo[i].wm = 0.0;
+ }
+ numshuryo = 0;
+
+ if( orlgth1 == 0 )
+ {
+ }
+
+
+ lgth1 = strlen( seq1[0] );
+ lgth2 = strlen( seq2[0] );
+
+ fprintf( stderr, "in suboptalign11 step 1\n" );
+
+ if( lgth1 > orlgth1 || lgth2 > orlgth2 )
+ {
+ int ll1, ll2;
+
+ fprintf( stderr, "in suboptalign11 step 1.3\n" );
+ if( orlgth1 > 0 && orlgth2 > 0 )
+ {
+ fprintf( stderr, "in suboptalign11 step 1.4\n" );
+ FreeFloatVec( w1 );
+ FreeFloatVec( w2 );
+// FreeFloatVec( match );
+ FreeFloatVec( initverticalw );
+ FreeFloatVec( lastverticalw );
+ fprintf( stderr, "in suboptalign11 step 1.5\n" );
+
+ FreeFloatVec( m );
+ FreeIntVec( mp );
+ FreeFloatVec( largeM );
+ FreeIntVec( Mp );
+ fprintf( stderr, "in suboptalign11 step 1.6\n" );
+
+
+ FreeFloatMtx( cpmx1 );
+ FreeFloatMtx( cpmx2 );
+
+ fprintf( stderr, "in suboptalign11 step 1.7\n" );
+ 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 );
+ largeM = AllocateFloatVec( ll2+2 );
+ Mp = AllocateIntVec( ll2+2 );
+
+ cpmx1 = AllocateFloatMtx( 26, ll1+2 );
+ cpmx2 = AllocateFloatMtx( 26, ll2+2 );
+
+ floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
+ intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
+
+ mseq1 = AllocateCharMtx( njob, ll1+ll2 );
+ mseq2 = AllocateCharMtx( njob, ll1+ll2 );
+
+#if DEBUG
+ fprintf( stderr, "succeeded\n" );
+#endif
+
+ orlgth1 = ll1 - 100;
+ orlgth2 = ll2 - 100;
+ }
+ fprintf( stderr, "in suboptalign11 step 1.6\n" );
+
+
+
+ fprintf( stderr, "in suboptalign11 step 2\n" );
+
+ if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
+ {
+ int ll1, ll2;
+
+ if( commonAlloc1 && commonAlloc2 )
+ {
+ FreeIntMtx( commonIP );
+ FreeIntMtx( commonJP );
+ FreeIntMtx( used );
+ }
+
+ ll1 = MAX( orlgth1, commonAlloc1 );
+ ll2 = MAX( orlgth2, commonAlloc2 );
+
+#if DEBUG
+ fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
+#endif
+
+ used = AllocateIntMtx( ll1+10, ll2+10 );
+ commonIP = AllocateIntMtx( ll1+10, ll2+10 );
+ commonJP = AllocateIntMtx( ll1+10, ll2+10 );
+
+#if DEBUG
+ fprintf( stderr, "succeeded\n\n" );
+#endif
+
+ commonAlloc1 = ll1;
+ commonAlloc2 = ll2;
+ }
+ ijpi = commonIP;
+ ijpj = commonJP;
+
+
+#if 0
+ for( i=0; i<lgth1; i++ )
+ fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
+#endif
+
+ fprintf( stderr, "in suboptalign11 step 3\n" );
+ currentw = w1;
+ previousw = w2;
+
+ match_calc( initverticalw, seq2, seq1, 0, lgth1 );
+
+ match_calc( currentw, seq1, seq2, 0, lgth2 );
+
+
+ lasti = lgth2+1;
+ for( j=1; j<lasti; ++j )
+ {
+ m[j] = currentw[j-1]; mp[j] = 0;
+ largeM[j] = currentw[j-1]; Mp[j] = 0;
+ }
+
+ lastverticalw[0] = currentw[lgth2-1];
+
+
+#if 0
+fprintf( stderr, "currentw = \n" );
+for( i=0; i<lgth1+1; i++ )
+{
+ fprintf( stderr, "%5.2f ", currentw[i] );
+}
+fprintf( stderr, "\n" );
+fprintf( stderr, "initverticalw = \n" );
+for( i=0; i<lgth2+1; i++ )
+{
+ fprintf( stderr, "%5.2f ", initverticalw[i] );
+}
+fprintf( stderr, "\n" );
+#endif
+#if DEBUG2
+ fprintf( stderr, "\n" );
+ fprintf( stderr, " " );
+ for( j=0; j<lgth2+1; j++ )
+ fprintf( stderr, "%c ", seq2[0][j] );
+ fprintf( stderr, "\n" );
+#endif
+
+ localstop = lgth1+lgth2+1;
+ maxwm = -999.9;
+ endali = endalj = 0;
+#if DEBUG2
+ fprintf( stderr, "\n" );
+ fprintf( stderr, "%c ", seq1[0][0] );
+
+ for( j=0; j<lgth2+1; j++ )
+ fprintf( stderr, "%5.0f ", currentw[j] );
+ fprintf( stderr, "\n" );
+#endif
+
+ lasti = lgth1+1;
+ for( i=1; i<lasti; i++ )
+ {
+ wtmp = previousw;
+ previousw = currentw;
+ currentw = wtmp;
+
+ previousw[0] = initverticalw[i-1];
+
+ match_calc( currentw, seq1, seq2, i, lgth2 );
+#if DEBUG2
+ fprintf( stderr, "%c ", seq1[0][i] );
+ fprintf( stderr, "%5.0f ", currentw[0] );
+#endif
+
+#if XXXXXXX
+fprintf( stderr, "\n" );
+fprintf( stderr, "i=%d\n", i );
+fprintf( stderr, "currentw = \n" );
+for( j=0; j<lgth2; j++ )
+{
+ fprintf( stderr, "%5.2f ", currentw[j] );
+}
+fprintf( stderr, "\n" );
+#endif
+#if XXXXXXX
+fprintf( stderr, "\n" );
+fprintf( stderr, "i=%d\n", i );
+fprintf( stderr, "currentw = \n" );
+for( j=0; j<lgth2; j++ )
+{
+ fprintf( stderr, "%5.2f ", currentw[j] );
+}
+fprintf( stderr, "\n" );
+#endif
+ currentw[0] = initverticalw[i];
+
+ mi = previousw[0]; mpi = 0;
+ Mi = previousw[0]; Mpi = 0;
+
+#if 0
+ if( mi < localthr ) mi = localthr2;
+#endif
+
+ ijpipt = ijpi[i] + 1;
+ ijpjpt = ijpj[i] + 1;
+ mjpt = m + 1;
+ Mjpt = largeM + 1;
+ prept = previousw;
+ curpt = currentw + 1;
+ mpjpt = mp + 1;
+ Mpjpt = Mp + 1;
+ tbk = -999999.9;
+ tbki = 0;
+ tbkj = 0;
+ lastj = lgth2+1;
+ for( j=1; j<lastj; j++ )
+ {
+ wm = *prept;
+ *ijpipt = i-1;
+ *ijpjpt = j-1;
+
+
+// fprintf( stderr, "i,j=%d,%d %c-%c\n", i, j, seq1[0][i], seq2[0][j] );
+// fprintf( stderr, "wm=%f\n", wm );
+#if 0
+ fprintf( stderr, "%5.0f->", wm );
+#endif
+ g = mi + fpenalty;
+#if 0
+ fprintf( stderr, "%5.0f?", g );
+#endif
+ if( g > wm )
+ {
+ wm = g;
+// *ijpipt = i - 1;
+ *ijpjpt = mpi;
+ }
+ g = *prept;
+ if( g > mi )
+ {
+ mi = g;
+ mpi = j-1;
+ }
+
+#if USE_PENALTY_EX
+ mi += fpenalty_ex;
+#endif
+
+#if 0
+ fprintf( stderr, "%5.0f->", wm );
+#endif
+ g = *mjpt + fpenalty;
+#if 0
+ fprintf( stderr, "m%5.0f?", g );
+#endif
+ if( g > wm )
+ {
+ wm = g;
+ *ijpipt = *mpjpt;
+// *ijpjpt = j - 1;
+ }
+ g = *prept;
+ if( g > *mjpt )
+ {
+ *mjpt = g;
+ *mpjpt = i-1;
+ }
+#if USE_PENALTY_EX
+ *mjpt += fpenalty_ex;
+#endif
+
+
+ g = tbk + fpenalty_OP;
+// g = tbk;
+ if( g > wm )
+ {
+ wm = g;
+ *ijpipt = tbki;
+ *ijpjpt = tbkj;
+// fprintf( stderr, "hit! i%d, j%d, ijpi = %d, ijpj = %d\n", i, j, *ijpipt, *ijpjpt );
+ }
+ g = Mi;
+ if( g > tbk )
+ {
+ tbk = g;
+ tbki = i-1;
+ tbkj = Mpi;
+ }
+ g = *Mjpt;
+ if( g > tbk )
+ {
+ tbk = g;
+ tbki = *Mpjpt;
+ tbkj = j-1;
+ }
+// tbk += fpenalty_EX;// + foffset;
+
+ g = *prept;
+ if( g > *Mjpt )
+ {
+ *Mjpt = g;
+ *Mpjpt = i-1;
+ }
+// *Mjpt += fpenalty_EX;// + foffset;
+
+ g = *prept;
+ if( g > Mi )
+ {
+ Mi = g;
+ Mpi = j-1;
+ }
+// Mi += fpenalty_EX;// + foffset;
+
+
+// fprintf( stderr, "wm=%f, tbk=%f(%c-%c), mi=%f, *mjpt=%f\n", wm, tbk, seq1[0][tbki], seq2[0][tbkj], mi, *mjpt );
+// fprintf( stderr, "ijp = %c,%c\n", seq1[0][abs(*ijpipt)], seq2[0][abs(*ijpjpt)] );
+
+
+ if( maxwm < wm )
+ {
+ maxwm = wm;
+ endali = i;
+ endalj = j;
+ }
+
+#if 1
+ if( numshuryo < 100 )
+ {
+ shuryo[numshuryo].i = i;
+ shuryo[numshuryo].j = j;
+ shuryo[numshuryo].wm = wm;
+
+ if( minshuryowm > wm )
+ {
+ minshuryowm = wm;
+ minshuryopos = numshuryo;
+ }
+ numshuryo++;
+ }
+ else
+ {
+ if( wm > minshuryowm )
+ {
+ shuryo[minshuryopos].i = i;
+ shuryo[minshuryopos].j = j;
+ shuryo[minshuryopos].wm = wm;
+ minshuryowm = wm;
+ for( k=0; k<100; k++ ) // muda
+ {
+ if( shuryo[k].wm < minshuryowm )
+ {
+ minshuryowm = shuryo[k].wm;
+ minshuryopos = k;
+ break;
+ }
+ }
+ }
+ }
+#endif
+#if 1
+ if( wm < localthr )
+ {
+// fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
+ *ijpipt = localstop;
+// *ijpjpt = localstop;
+ wm = localthr2;
+ }
+#endif
+#if 0
+ fprintf( stderr, "%5.0f ", *curpt );
+#endif
+#if DEBUG2
+ fprintf( stderr, "%5.0f ", wm );
+// fprintf( stderr, "%c-%c *ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], *ijppt, localstop );
+#endif
+
+ *curpt += wm;
+ ijpipt++;
+ ijpjpt++;
+ mjpt++;
+ Mjpt++;
+ prept++;
+ mpjpt++;
+ Mpjpt++;
+ curpt++;
+ }
+#if DEBUG2
+ fprintf( stderr, "\n" );
+#endif
+
+ lastverticalw[i] = currentw[lgth2-1];
+ }
+
+ for( k=0; k<100; k++ )
+ {
+ fprintf( stderr, "shuryo[%d].i,j,wm = %d,%d,%f\n", k, shuryo[k].i, shuryo[k].j, shuryo[k].wm );
+ }
+
+
+#if 1
+ fprintf( stderr, "maxwm = %f\n", maxwm );
+ fprintf( stderr, "endali = %d\n", endali );
+ fprintf( stderr, "endalj = %d\n", endalj );
+#endif
+
+ qsort( shuryo, 100, sizeof( Shuryoten ), (int (*)())compshuryo );
+ for( k=0; k<100; k++ )
+ {
+ fprintf( stderr, "shuryo[%d].i,j,wm = %d,%d,%f\n", k, shuryo[k].i, shuryo[k].j, shuryo[k].wm );
+ }
+
+
+ lasti = lgth1+1;
+ for( i=0; i<lasti; i++ )
+ {
+ ijpi[i][0] = localstop;
+ ijpj[i][0] = localstop;
+ }
+ lastj = lgth2+1;
+ for( j=0; j<lastj; j++ )
+ {
+ ijpi[0][j] = localstop;
+ ijpj[0][j] = localstop;
+ }
+
+ for( i=0; i<lasti; i++ ) for( j=0; j<lastj; j++ ) used[i][j] = 0;
+
+ for( k=0; k<numshuryo; k++ )
+ {
+ if( shuryo[k].wm < shuryo[0].wm * 0.3 ) break;
+ fprintf( stderr, "k=%d, shuryo[k].i,j,wm=%d,%d,%f go\n", k, shuryo[k].i, shuryo[k].j, shuryo[k].wm );
+ resf = gentracking( used, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijpi, ijpj, off1pt, off2pt, shuryo[k].i, shuryo[k].j );
+ if( resf == -1.0 ) continue;
+ putlocalhom3( mseq1[0], mseq2[0], lhmpt, *off1pt, *off2pt, (int)shuryo[k].wm, strlen( mseq1[0] ) );
+#if 0
+ fprintf( stderr, "\n" );
+ fprintf( stderr, ">\n%s\n", mseq1[0] );
+ fprintf( stderr, ">\n%s\n", mseq2[0] );
+#endif
+ }
+ for( i=0; i<20; i++ )
+ {
+ for( j=0; j<20; j++ )
+ {
+ fprintf( stderr, "%2d ", used[i][j] );
+ }
+ fprintf( stderr, "\n" );
+ }
+
+
+// fprintf( stderr, "### impmatch = %f\n", *impmatch );
+
+ resultlen = strlen( mseq1[0] );
+ if( alloclen < resultlen || resultlen > N )
+ {
+ fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
+ ErrorExit( "LENGTH OVER!\n" );
+ }
+
+
+
+
+
+ return( wm );
+}
+
Added: test-suite/trunk/MultiSource/Benchmarks/mafft/tddis.c
URL: http://llvm.org/viewvc/llvm-project/test-suite/trunk/MultiSource/Benchmarks/mafft/tddis.c?rev=65043&view=auto
==============================================================================
--- test-suite/trunk/MultiSource/Benchmarks/mafft/tddis.c (added)
+++ test-suite/trunk/MultiSource/Benchmarks/mafft/tddis.c Thu Feb 19 06:37:59 2009
@@ -0,0 +1,827 @@
+#include "mltaln.h"
+
+#define DEBUG 0
+
+#if 0
+void mdfymtx( char pair[njob][njob], int s1, double **partialmtx, double **mtx )
+#else
+void mdfymtx( char **pair, int s1, double **partialmtx, double **mtx )
+#endif
+{
+ int i, j;
+ int icount, jcount;
+#if DEBUG
+ FILE *fp;
+ static char name[M][B];
+
+ for( i=0; i<M; i++ ) name[i][0] = 0;
+ fprintf( stdout, "s1 = %d\n", s1 );
+ for( i=0; i<njob; i++ )
+ {
+ for( j=0; j<njob; j++ )
+ {
+ printf( "%#2d", pair[i][j] );
+ }
+ printf( "\n" );
+ }
+#endif
+
+ for( i=0, icount=0; i<njob-1; i++ )
+ {
+ if( !pair[s1][i] ) continue;
+ for( j=i+1, jcount=icount+1; j<njob; j++ )
+ {
+ if( !pair[s1][j] ) continue;
+ partialmtx[icount][jcount] = mtx[i][j];
+ jcount++;
+ }
+ icount++;
+ }
+#if DEBUG
+ fp = fopen( "hat2.org", "w" );
+ WriteHat2( fp, njob, name, mtx );
+ fclose( fp );
+ fp = fopen( "hat2.mdf", "w" );
+ WriteHat2( fp, icount, name, partialmtx );
+ fclose( fp );
+#endif
+
+}
+
+
+float score_calc( char **seq, int s ) /* method 3 */
+{
+ int i, j, k, c;
+ int len = strlen( seq[0] );
+ float score;
+ int tmpscore;
+ char *mseq1, *mseq2;
+
+ score = 0.0;
+ for( i=0; i<s-1; i++ )
+ {
+ for( j=i+1; j<s; j++ )
+ {
+ mseq1 = seq[i];
+ mseq2 = seq[j];
+ tmpscore = 0;
+ c = 0;
+ for( k=0; k<len; k++ )
+ {
+ if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
+ c++;
+ tmpscore += amino_dis[(int)mseq1[k]][(int)mseq2[k]];
+ if( mseq1[k] == '-' )
+ {
+ tmpscore += penalty;
+ while( mseq1[++k] == '-' )
+ ;
+ k--;
+ if( k > len-2 ) break;
+ continue;
+ }
+ if( mseq2[k] == '-' )
+ {
+ tmpscore += penalty;
+ while( mseq2[++k] == '-' )
+ ;
+ k--;
+ if( k > len-2 ) break;
+ continue;
+ }
+ }
+ /*
+ if( mseq1[0] == '-' || mseq2[0] == '-' )
+ {
+ for( k=0; k<len; k++ )
+ {
+ if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
+ if( !( mseq1[k] != '-' && mseq2[k] != '-' ) )
+ {
+ c--;
+ tmpscore -= penalty;
+ break;
+ }
+ else break;
+ }
+ }
+ if( mseq1[len-1] == '-' || mseq2[len-1] == '-' )
+ {
+ for( k=0; k<len; k++ )
+ {
+ if( mseq1[k] == '-' && mseq2[k] == '-' ) continue;
+ if( !( mseq1[k] != '-' && mseq2[k] != '-' ) )
+ {
+ c--;
+ tmpscore -= penalty;
+ break;
+ }
+ else break;
+ }
+ }
+ */
+ score += (double)tmpscore / (double)c;
+ }
+ }
+ score = (float)score / ( ( (double)s * ((double)s-1.0) ) / 2.0 );
+ fprintf( stderr, "score in score_calc = %f\n", score );
+ return( score );
+}
+
+void cpmx_calc( char **seq, float **cpmx, double *eff, int lgth, int clus )
+{
+ int i, j, k;
+ double totaleff = 0.0;
+
+ for( i=0; i<clus; i++ ) totaleff += eff[i];
+ for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
+ for( j=0; j<lgth; j++ ) for( k=0; k<clus; k++ )
+ cpmx[(int)amino_n[(int)seq[k][j]]][j] += (float)eff[k] / totaleff;
+}
+
+
+void cpmx_calc_new_bk( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
+{
+ int i, j, k;
+ float feff;
+
+ for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
+ for( k=0; k<clus; k++ )
+ {
+ feff = (float)eff[k];
+ for( j=0; j<lgth; j++ )
+ {
+ cpmx[(int)amino_n[(int)seq[k][j]]][j] += feff;
+ }
+ }
+}
+
+void cpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
+{
+ int i, j, k;
+ float feff;
+ float *cpmxpt, **cpmxptpt;
+ char *seqpt;
+
+ j = 26;
+ cpmxptpt = cpmx;
+ while( j-- )
+ {
+ cpmxpt = *cpmxptpt++;
+ i = lgth;
+ while( i-- )
+ *cpmxpt++ = 0.0;
+ }
+ for( k=0; k<clus; k++ )
+ {
+ feff = (float)eff[k];
+ seqpt = seq[k];
+// fprintf( stderr, "seqpt = %s, lgth=%d\n", seqpt, lgth );
+ for( j=0; j<lgth; j++ )
+ {
+ cpmx[(int)amino_n[(int)*seqpt++]][j] += feff;
+ }
+ }
+}
+void MScpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
+{
+ int i, j, k;
+ float feff;
+ float **cpmxptpt, *cpmxpt;
+ char *seqpt;
+
+ j = lgth;
+ cpmxptpt = cpmx;
+ while( j-- )
+ {
+ cpmxpt = *cpmxptpt++;
+ i = 26;
+ while( i-- )
+ *cpmxpt++ = 0.0;
+ }
+ for( k=0; k<clus; k++ )
+ {
+ feff = (float)eff[k];
+ seqpt = seq[k];
+ cpmxptpt = cpmx;
+ j = lgth;
+ while( j-- )
+ (*cpmxptpt++)[(int)amino_n[(int)*seqpt++]] += feff;
+ }
+#if 0
+ for( j=0; j<lgth; j++ ) for( i=0; i<26; i++ ) cpmx[j][i] = 0.0;
+ for( k=0; k<clus; k++ )
+ {
+ feff = (float)eff[k];
+ for( j=0; j<lgth; j++ )
+ cpmx[j][(int)amino_n[(int)seq[k][j]]] += feff;
+ }
+#endif
+}
+
+void cpmx_ribosum( char **seq, char **seqr, char *dir, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
+{
+ int i, j, k;
+ float feff;
+ float **cpmxptpt, *cpmxpt;
+ char *seqpt, *seqrpt, *dirpt;
+ int code, code1, code2;
+
+ j = lgth;
+ cpmxptpt = cpmx;
+ while( j-- )
+ {
+ cpmxpt = *cpmxptpt++;
+ i = 37;
+ while( i-- )
+ *cpmxpt++ = 0.0;
+ }
+ for( k=0; k<clus; k++ )
+ {
+ feff = (float)eff[k];
+ seqpt = seq[k];
+ seqrpt = seqr[k];
+ dirpt = dir;
+ cpmxptpt = cpmx;
+ j = lgth;
+ while( j-- )
+ {
+#if 0
+ code1 = amino_n[(int)*seqpt];
+ if( code1 > 3 ) code = 36;
+ else
+ code = code1;
+#else
+ code1 = amino_n[(int)*seqpt];
+ code2 = amino_n[(int)*seqrpt];
+ if( code1 > 3 )
+ {
+ code = 36;
+ }
+ else if( code2 > 3 )
+ {
+ code = code1;
+ }
+ else if( *dirpt == '5' )
+ {
+ code = 4 + code2 * 4 + code1;
+ }
+ else if( *dirpt == '3' )
+ {
+ code = 20 + code2 * 4 + code1;
+ }
+ else // if( *dirpt == 'o' ) // nai
+ {
+ code = code1;
+ }
+#endif
+
+// fprintf( stderr, "%c -> code=%d toa=%d, tog=%d, toc=%d, tot=%d, ton=%d, efee=%f\n", *seqpt, code%4, ribosumdis[code][4+0], ribosumdis[code][4+1], ribosumdis[code][4+2], ribosumdis[code][20+3], ribosumdis[code][36], feff );
+
+ seqpt++;
+ seqrpt++;
+ dirpt++;
+
+ (*cpmxptpt++)[code] += feff;
+ }
+ }
+}
+
+void mseqcat( char **seq1, char **seq2, double **eff, double *effarr1, double *effarr2, char name1[M][B], char name2[M][B], int clus1, int clus2 )
+{
+ int i, j;
+ for( i=0; i<clus2; i++ )
+ seq1[clus1+i] = seq2[i];
+ for( i=0; i<clus2; i++ ) strcpy( name1[clus1+i], name2[i] );
+
+ for( i=0; i<clus1; i++ ) for( j=0; j<clus1; j++ ) eff[i][j] = effarr1[i]* effarr1[j];
+ for( i=0; i<clus1; i++ ) for( j=clus1; j<clus1+clus2; j++ ) eff[i][j] = effarr1[i]* effarr2[j-clus1];
+ for( i=clus1; i<clus1+clus2; i++ ) for( j=0; j<clus1; j++ ) eff[i][j] = effarr2[i-clus1] * effarr1[j];
+ for( i=clus1; i<clus1+clus2; i++ ) for( j=clus1; j<clus1+clus2; j++ ) eff[i][j] = effarr2[i-clus1] * effarr2[j-clus1];
+}
+
+
+
+void strnbcat( char *s1, char *s2, int m ) /* s1 + s2 ---> s2 */
+{
+ static char b[N];
+ strncpy( b, s1, m ); b[m] = 0;
+ strcat( b, s2 );
+ strcpy( s2, b );
+}
+
+#if 0
+int conjuction( char pair[njob][njob], int s, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
+#else
+int conjuctionforgaln( int s0, int s1, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
+#endif
+{
+ int m, k;
+ char b[B];
+ double total;
+
+#if 0
+ fprintf( stderr, "s0 = %d, s1 = %d\n", s0, s1 );
+#endif
+
+ total = 0.0;
+ d[0] = 0;
+ for( m=s0, k=0; m<s1; m++ )
+ {
+ {
+ sprintf( b, " %d", m+1 );
+#if 1
+ if( strlen( d ) < 100 ) strcat( d, b );
+#else
+ strcat( d, b );
+#endif
+ aseq[k] = seq[m];
+ peff[k] = eff[m];
+ total += peff[k];
+#if 0
+ strcpy( aname[k], name[m] );
+#endif
+ k++;
+ }
+ }
+#if 1
+ for( m=0; m<k; m++ )
+ {
+ peff[m] /= total;
+// fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
+ }
+#endif
+ return( k );
+}
+
+void makegrouprna( RNApair ***group, RNApair ***all, int *memlist )
+{
+ int k, m;
+ for( k=0; (m=*memlist)!=-1; memlist++, k++ )
+ {
+ group[k] = all[m];
+ }
+}
+
+void makegrouprnait( RNApair ***group, RNApair ***all, char **pair, int s )
+{
+ int k, m;
+ for( m=s, k=0; m<njob; m++ )
+ {
+ if( pair[s][m] != 0 )
+ {
+ group[k++] = all[m];
+ }
+ }
+}
+
+int fastconjuction_noweight( int *memlist, char **seq, char **aseq, double *peff, char *d )
+{
+ int m, k, dln;
+ char b[B];
+ double total;
+
+#if DEBUG
+ fprintf( stderr, "s = %d\n", s );
+#endif
+
+ total = 0.0;
+ d[0] = 0;
+ dln = 0;
+ for( k=0; *memlist!=-1; memlist++, k++ )
+ {
+ m = *memlist;
+ dln += sprintf( b, " %d", m+1 );
+#if 1
+ if( dln < 100 ) strcat( d, b );
+#else
+ strcat( d, b );
+#endif
+ aseq[k] = seq[m];
+ peff[k] = 1.0;
+ total += peff[k];
+ }
+#if 1
+ for( m=0; m<k; m++ )
+ peff[m] /= total;
+#endif
+ return( k );
+}
+int fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d )
+{
+ int m, k, dln;
+ char b[B];
+ double total;
+
+#if DEBUG
+ fprintf( stderr, "s = %d\n", s );
+#endif
+
+ total = 0.0;
+ d[0] = 0;
+ dln = 0;
+ for( k=0; *memlist!=-1; memlist++, k++ )
+ {
+ m = *memlist;
+ dln += sprintf( b, " %d", m+1 );
+#if 1
+ if( dln < 100 ) strcat( d, b );
+#else
+ strcat( d, b );
+#endif
+ aseq[k] = seq[m];
+ peff[k] = eff[m];
+ total += peff[k];
+ }
+#if 1
+ for( m=0; m<k; m++ )
+ {
+// fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
+ peff[m] /= total;
+ }
+#endif
+ return( k );
+}
+
+int fastconjuction( int *memlist, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
+{
+ int m, k, dln;
+ char b[B];
+ double total;
+
+#if DEBUG
+ fprintf( stderr, "s = %d\n", s );
+#endif
+
+ total = 0.0;
+ d[0] = 0;
+ dln = 0;
+ for( k=0; *memlist!=-1; memlist++, k++ )
+ {
+ m = *memlist;
+ dln += sprintf( b, " %d", m+1 );
+#if 1
+ if( dln < 100 ) strcat( d, b );
+#else
+ strcat( d, b );
+#endif
+ aseq[k] = seq[m];
+ peff[k] = eff[m];
+ total += peff[k];
+#if 0
+ strcpy( aname[k], name[m] );
+#endif
+ }
+#if 1
+ for( m=0; m<k; m++ )
+ peff[m] /= total;
+#endif
+ return( k );
+}
+int conjuctionfortbfast( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
+{
+ int m, k;
+ char b[B];
+ double total;
+
+#if DEBUG
+ fprintf( stderr, "s = %d\n", s );
+#endif
+
+ total = 0.0;
+ d[0] = 0;
+ for( m=s, k=0; m<njob; m++ )
+ {
+ if( pair[s][m] != 0 )
+ {
+ sprintf( b, " %d", m+1 );
+#if 1
+ if( strlen( d ) < 100 ) strcat( d, b );
+#else
+ strcat( d, b );
+#endif
+ aseq[k] = seq[m];
+ peff[k] = eff[m];
+ total += peff[k];
+#if 0
+ strcpy( aname[k], name[m] );
+#endif
+ k++;
+ }
+ }
+#if 1
+ for( m=0; m<k; m++ )
+ peff[m] /= total;
+#endif
+ return( k );
+}
+int conjuction( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
+{
+ int m, k;
+ char b[B];
+ double total;
+
+#if DEBUG
+ fprintf( stderr, "s = %d\n", s );
+#endif
+
+ total = 0.0;
+ d[0] = 0;
+ for( m=s, k=0; m<njob; m++ )
+ {
+ if( pair[s][m] != 0 )
+ {
+ sprintf( b, " %d", m+1 );
+#if 1
+ if( strlen( d ) < 100 ) strcat( d, b );
+#else
+ strcat( d, b );
+#endif
+ aseq[k] = seq[m];
+ peff[k] = eff[m];
+ total += peff[k];
+#if 0
+ strcpy( aname[k], name[m] );
+#endif
+ k++;
+ }
+ }
+#if 0
+ for( m=0; m<k; m++ )
+ peff[m] /= total;
+#endif
+ return( k );
+}
+
+void floatdelete( float **cpmx, int d, int len )
+{
+ int i, j;
+
+ for( i=d; i<len-1; i++ )
+ {
+ for( j=0; j<26; j++ )
+ {
+ cpmx[j][i] = cpmx[j][i+1];
+ }
+ }
+}
+
+void chardelete( char *seq, int d )
+{
+ char b[N];
+
+ strcpy( b, seq+d+1 );
+ strcpy( seq+d, b );
+}
+
+int RootBranchNode( int nseq, int ***topol, int step, int branch )
+{
+ int i, j, s1, s2, value;
+
+ value = 1;
+ for( i=step+1; i<nseq-2; i++ )
+ {
+ for( j=0; (s1=topol[i][0][j])>-1; j++ )
+ if( s1 == topol[step][branch][0] ) value++;
+ for( j=0; (s2=topol[i][1][j])>-1; j++ )
+ if( s2 == topol[step][branch][0] ) value++;
+ }
+ return( value );
+}
+
+void BranchLeafNode( int nseq, int ***topol, int *node, int step, int branch )
+{
+ int i, j, k, s;
+
+ for( i=0; i<nseq; i++ ) node[i] = 0;
+ for( i=0; i<step-1; i++ )
+ for( k=0; k<2; k++ )
+ for( j=0; (s=topol[i][k][j])>-1; j++ )
+ node[s]++;
+ for( k=0; k<branch+1; k++ )
+ for( j=0; (s=topol[step][k][j])>-1; j++ )
+ node[s]++;
+}
+
+void RootLeafNode( int nseq, int ***topol, int *node )
+{
+ int i, j, k, s;
+ for( i=0; i<nseq; i++ ) node[i] = 0;
+ for( i=0; i<nseq-2; i++ )
+ for( k=0; k<2; k++ )
+ for( j=0; (s=topol[i][k][j])>-1; j++ )
+ node[s]++;
+}
+
+void nodeFromABranch( int nseq, int *result, int **pairwisenode, int ***topol, double **len, int step, int num )
+{
+ int i, s, count;
+ int *innergroup;
+ int *outergroup1;
+#if 0
+ int outergroup2[nseq];
+ int table[nseq];
+#else
+ static int *outergroup2 = NULL;
+ static int *table = NULL;
+ if( outergroup2 == NULL )
+ {
+ outergroup2 = AllocateIntVec( nseq );
+ table = AllocateIntVec( nseq );
+ }
+#endif
+ innergroup = topol[step][num];
+ outergroup1 = topol[step][!num];
+ for( i=0; i<nseq; i++ ) table[i] = 1;
+ for( i=0; (s=innergroup[i])>-1; i++ ) table[s] = 0;
+ for( i=0; (s=outergroup1[i])>-1; i++ ) table[s] = 0;
+ for( i=0, count=0; i<nseq; i++ )
+ {
+ if( table[i] )
+ {
+ outergroup2[count] = i;
+ count++;
+ }
+ }
+ outergroup2[count] = -1;
+
+ for( i=0; (s=innergroup[i])>-1; i++ )
+ {
+ result[s] = pairwisenode[s][outergroup1[0]]
+ + pairwisenode[s][outergroup2[0]]
+ - pairwisenode[outergroup1[0]][outergroup2[0]] - 1;
+ result[s] /= 2;
+ }
+ for( i=0; (s=outergroup1[i])>-1; i++ )
+ {
+ result[s] = pairwisenode[s][outergroup2[0]]
+ + pairwisenode[s][innergroup[0]]
+ - pairwisenode[innergroup[0]][outergroup2[0]] + 1;
+ result[s] /= 2;
+ }
+ for( i=0; (s=outergroup2[i])>-1; i++ )
+ {
+ result[s] = pairwisenode[s][outergroup1[0]]
+ + pairwisenode[s][innergroup[0]]
+ - pairwisenode[innergroup[0]][outergroup1[0]] + 1;
+ result[s] /= 2;
+ }
+
+#if 0
+ for( i=0; i<nseq; i++ )
+ fprintf( stderr, "result[%d] = %d\n", i+1, result[i] );
+#endif
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#if 0
+void OneClusterAndTheOther( int locnjob, char pair[njob][njob], int *s1, int *s2, int ***topol, int step, int branch )
+#else
+void OneClusterAndTheOther( int locnjob, char **pair, int *s1, int *s2, int ***topol, int step, int branch )
+#endif
+{
+ int i;
+ int r1;
+
+ *s1 = topol[step][branch][0];
+ for( i=0; (r1=topol[step][branch][i])>-1; i++ )
+ pair[*s1][r1] = 1;
+ for( i=0; i<locnjob; i++ )
+ {
+ if( !pair[*s1][i] )
+ {
+ *s2 = i;
+ break;
+ }
+ }
+ for( i=*s2; i<locnjob; i++ )
+ {
+ if( !pair[*s1][i] )
+ pair[*s2][i] = 1;
+ }
+}
+
+void makeEffMtx( int nseq, double **mtx, double *vec )
+{
+ int i, j;
+ for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
+ mtx[i][j] = vec[i] * vec[j];
+}
+
+void node_eff( int nseq, double *eff, int *node )
+{
+ int i;
+ extern double ipower( double, int );
+ for( i=0; i<nseq; i++ )
+ eff[i] = ipower( 0.5, node[i] ) + geta2;
+ /*
+ eff[i] = ipower( 0.5, node[i] ) + 0;
+ */
+#if DEBUG
+ for( i=0; i<nseq; i++ )
+#endif
+}
+
+
+int shrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
+{
+ int m1, k1, m2, k2;
+
+ for( m1=s1, k1=0; m1<njob; m1++ )
+ {
+ if( pair[s1][m1] != 0 )
+ {
+ for( m2=s2, k2=0; m2<njob; m2++ )
+ {
+ if( pair[s2][m2] != 0 )
+ {
+ if( localhom[m1][m2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[m1]+m2;
+ k2++;
+ }
+ }
+ k1++;
+ }
+ }
+ return( 0 );
+}
+int msshrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
+{
+ int m1, k1, n1, m2, k2, n2;
+
+ for( m1=s1, k1=0; m1<njob; m1++ )
+ {
+ if( pair[s1][m1] != 0 )
+ {
+ for( m2=s2, k2=0; m2<njob; m2++ )
+ {
+ if( pair[s2][m2] != 0 )
+ {
+ n1 = MIN(m1,m2); n2=MAX(m1,m2);
+ if( localhom[m1][m2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[n1]+n2;
+ k2++;
+ }
+ }
+ k1++;
+ }
+ }
+ return( 0 );
+}
+
+int fastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
+{
+ int k1, k2;
+ int *intpt1, *intpt2;
+
+
+ for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
+ {
+ for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
+ {
+ if( localhom[*intpt1][*intpt2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2;
+ }
+ }
+ return( 0 );
+}
+
+int msfastshrinklocalhom( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink )
+{
+ int k1, k2;
+ int *intpt1, *intpt2;
+ int m1, m2;
+
+ for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
+ {
+ for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
+ {
+ m1 = MIN(*intpt1,*intpt2); m2 = MAX(*intpt1,*intpt2);
+ if( localhom[m1][m2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[m1]+m2;
+ }
+ }
+ return( 0 );
+}
+
More information about the llvm-commits
mailing list