#include "mltaln.h"
#define DEBUG 0
+#define USEDISTONTREE 1
#if 0
void mdfymtx( char pair[njob][njob], int s1, double **partialmtx, double **mtx )
}
-float score_calc( char **seq, int s ) /* method 3 */
+double score_calc( char **seq, int s ) /* method 3 */
{
int i, j, k, c;
int len = strlen( seq[0] );
- float score;
+ double score;
int tmpscore;
char *mseq1, *mseq2;
score += (double)tmpscore / (double)c;
}
}
- score = (float)score / ( ( (double)s * ((double)s-1.0) ) / 2.0 );
+ score = (double)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 )
+void cpmx_calc( char **seq, double **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( i=0; i<nalphabets; 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;
+ cpmx[(int)amino_n[(unsigned char)seq[k][j]]][j] += (double)eff[k] / totaleff;
}
-void cpmx_calc_new_bk( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
+void cpmx_calc_new_bk( char **seq, double **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
{
int i, j, k;
- float feff;
+ double feff;
- for( i=0; i<26; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
+ for( i=0; i<nalphabets; i++ ) for( j=0; j<lgth; j++ ) cpmx[i][j] = 0.0;
for( k=0; k<clus; k++ )
{
- feff = (float)eff[k];
+ feff = (double)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
+void cpmx_calc_add( char **seq, double **cpmx, double *eff, int lgth, int clus ) // lastmem = newmem; summ eff must be 1.0
+{
+ double neweff, orieff;
+ int newmem, i, j;
+
+ newmem = clus-1;
+ neweff = eff[clus-1];
+ orieff = 1.0-neweff;
+#if 1 // TESTING Feb/1/18:00
+ for( j=0; j<lgth; j++ )
+ {
+ for( i=0;i<nalphabets; i++ ) cpmx[i][j] *= orieff;
+ cpmx[(unsigned char)amino_n[(unsigned char)seq[newmem][j]]][j] += neweff;
+ }
+#else // possibly faster?
+ for( i=0;i<nalphabets; i++ )
+ {
+ for( j=0; j<lgth; j++ ) cpmx[i][j] *= orieff;
+ }
+ for( j=0; j<lgth; j++ ) cpmx[(unsigned char)amino_n[(unsigned char)seq[newmem][j]]][j] += neweff;
+#endif
+}
+
+void cpmx_calc_new( char **seq, double **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
{
int i, j, k;
- float feff;
- float *cpmxpt, **cpmxptpt;
+ double feff;
+ double *cpmxpt, **cpmxptpt;
char *seqpt;
- j = 26;
+ j = nalphabets;
cpmxptpt = cpmx;
while( j-- )
{
}
for( k=0; k<clus; k++ )
{
- feff = (float)eff[k];
+ feff = (double)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;
+ cpmx[(unsigned char)amino_n[(unsigned char)*seqpt++]][j] += feff;
}
}
}
-void MScpmx_calc_new( char **seq, float **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
+void MScpmx_calc_new( char **seq, double **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
{
int i, j, k;
- float feff;
- float **cpmxptpt, *cpmxpt;
+ double feff;
+ double **cpmxptpt, *cpmxpt;
char *seqpt;
j = lgth;
while( j-- )
{
cpmxpt = *cpmxptpt++;
- i = 26;
+ i = nalphabets;
while( i-- )
*cpmxpt++ = 0.0;
}
for( k=0; k<clus; k++ )
{
- feff = (float)eff[k];
+ feff = (double)eff[k];
seqpt = seq[k];
cpmxptpt = cpmx;
j = lgth;
while( j-- )
- (*cpmxptpt++)[(int)amino_n[(int)*seqpt++]] += feff;
+ (*cpmxptpt++)[(int)amino_n[(unsigned char)*seqpt++]] += feff;
}
#if 0
- for( j=0; j<lgth; j++ ) for( i=0; i<26; i++ ) cpmx[j][i] = 0.0;
+ for( j=0; j<lgth; j++ ) for( i=0; i<nalphabets; i++ ) cpmx[j][i] = 0.0;
for( k=0; k<clus; k++ )
{
- feff = (float)eff[k];
+ feff = (double)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
+void cpmx_ribosum( char **seq, char **seqr, char *dir, double **cpmx, double *eff, int lgth, int clus ) // summ eff must be 1.0
{
int i, j, k;
- float feff;
- float **cpmxptpt, *cpmxpt;
+ double feff;
+ double **cpmxptpt, *cpmxpt;
char *seqpt, *seqrpt, *dirpt;
int code, code1, code2;
}
for( k=0; k<clus; k++ )
{
- feff = (float)eff[k];
+ feff = (double)eff[k];
seqpt = seq[k];
seqrpt = seqr[k];
dirpt = dir;
else
code = code1;
#else
- code1 = amino_n[(int)*seqpt];
- code2 = amino_n[(int)*seqrpt];
+ code1 = amino_n[(unsigned char)*seqpt];
+ code2 = amino_n[(unsigned char)*seqrpt];
if( code1 > 3 )
{
code = 36;
}
}
-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 fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d )
+#if 0
+int fastconjuction_target( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d, double mineff, int *targetmap )
{
int m, k, dln;
char b[B];
double total;
+ int *memlistbk = memlist;
#if DEBUG
fprintf( stderr, "s = %d\n", s );
strcat( d, b );
#endif
aseq[k] = seq[m];
- peff[k] = eff[m];
+ if( eff[m] < mineff )
+ peff[k] = mineff;
+ else
+ peff[k] = eff[m];
+
total += peff[k];
}
#if 1
for( m=0; m<k; m++ )
{
-// fprintf( stderr, "peff[%d] = %f\n", m, peff[m] );
+// fprintf( stderr, "Apr17 peff[%d] = %f\n", m, peff[m] );
peff[m] /= total;
}
#endif
+
return( k );
}
+#endif
-int fastconjuction( int *memlist, char **seq, char **aseq, double *peff, double *eff, char name[M][B], char aname[M][B], char *d )
+
+int fastconjuction_noname( int *memlist, char **seq, char **aseq, double *peff, double *eff, char *d, double mineff )
{
int m, k, dln;
char b[B];
strcat( d, b );
#endif
aseq[k] = seq[m];
- peff[k] = eff[m];
+ if( eff[m] < mineff )
+ peff[k] = mineff;
+ else
+ peff[k] = eff[m];
+
total += peff[k];
-#if 0
- strcpy( aname[k], name[m] );
-#endif
}
#if 1
for( m=0; m<k; m++ )
+ {
+// fprintf( stderr, "Apr17 peff[%d] = %f\n", m, peff[m] );
peff[m] /= total;
+ }
#endif
return( k );
}
-
-int conjuctionfortbfast_kozo( double *tmptmptmp, char **pair, int s, char **seq, char **aseq, double *peff, double *eff, double *peff_kozo, double *eff_kozo, char *d )
+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;
+ int m, k, dln;
char b[B];
double total;
- double total_kozo;
#if DEBUG
fprintf( stderr, "s = %d\n", s );
#endif
total = 0.0;
-// total_kozo = 0.0;
- total_kozo = *tmptmptmp; // masaka
d[0] = 0;
- for( m=s, k=0; m<njob; m++ )
+ dln = 0;
+ for( k=0; *memlist!=-1; memlist++, k++ )
{
- if( pair[s][m] != 0 )
- {
- sprintf( b, " %d", m+1 );
+ m = *memlist;
+ dln += sprintf( b, " %d", m+1 );
#if 1
- if( strlen( d ) < 100 ) strcat( d, b );
+ if( dln < 100 ) strcat( d, b );
#else
- strcat( d, b );
+ strcat( d, b );
#endif
- aseq[k] = seq[m];
- peff[k] = eff[m];
- peff_kozo[k] = eff_kozo[m];
- total += peff[k];
- total_kozo += peff_kozo[k];
+ 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;
- if( total_kozo > 0.0 )
- {
- for( m=0; m<k; m++ )
- {
- peff_kozo[m] /= total_kozo;
- if( peff_kozo[m] > 0.0 ) peff_kozo[m] += peff[m];
- }
- }
- else //iranai
- {
- for( m=0; m<k; m++ ) peff_kozo[m] = 0.0;
- }
#endif
-// fprintf( stderr, "\n\ndvtditr_total_kozo = %f\n\n", total_kozo );
- *tmptmptmp = total_kozo;
return( k );
}
-int conjuctionfortbfast( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
+
+
+int conjuctionfortbfast_old( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char *d )
{
int m, k;
char *b;
free( b );
return( k );
}
+
int conjuction( char **pair, int s, char **seq, char **aseq, double *peff, double *eff, char **name, char **aname, char *d )
{
int m, k;
return( k );
}
-void floatdelete( float **cpmx, int d, int len )
+void doubledelete( double **cpmx, int d, int len )
{
int i, j;
for( i=d; i<len-1; i++ )
{
- for( j=0; j<26; j++ )
+ for( j=0; j<nalphabets; j++ )
{
cpmx[j][i] = cpmx[j][i+1];
}
-
-
-#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
+void OneClusterAndTheOther_fast( int locnjob, int *memlist1, int *memlist2, int *s1, int *s2, char *pair, int ***topol, int step, int branch, double **smalldistmtx, double **distmtx, double *distontree )
{
- int i;
+ int i, k, j;
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++ )
+// char *pair;
+
+// pair = calloc( locnjob, sizeof( char ) );
+
+ for( i=0; i<locnjob; i++ ) pair[i] = 0;
+ for( i=0, k=0; (r1=topol[step][branch][i])>-1; i++ )
+ {
+ pair[r1] = 1;
+ memlist1[k++] = r1;
+ }
+ memlist1[k] = -1;
+
+ for( i=0, k=0; i<locnjob; i++ )
{
- if( !pair[*s1][i] )
+ if( !pair[i] )
{
- *s2 = i;
- break;
+ memlist2[k++] = i;
}
}
- for( i=*s2; i<locnjob; i++ )
- {
- if( !pair[*s1][i] )
- pair[*s2][i] = 1;
- }
+ memlist2[k] = -1;
+
+ *s1 = memlist1[0];
+ *s2 = memlist2[0];
+
+ if( smalldistmtx )
+ {
+ int im, jm;
+#if USEDISTONTREE
+ for( i=0; (im=memlist1[i])!=-1; i++ ) for( j=0; (jm=memlist2[j])!=-1; j++ )
+ smalldistmtx[i][j] = distontree[im] + distontree[jm];
+#else
+ for( i=0; (im=memlist1[i])!=-1; i++ ) for( j=0; (jm=memlist2[j])!=-1; j++ )
+ smalldistmtx[i][j] = distmtx[MIN(im,jm)][MAX(im,jm)];
+#endif
+
+#if 0
+ reporterr( "\n" );
+ for( i=0; (im=memlist1[i])!=-1; i++ ) for( j=0; (jm=memlist2[j])!=-1; j++ )
+ reporterr( "smalldistmtx[%d][%d] = %f\n", i, j, smalldistmtx[i][j] );
+
+
+ for( i=0; (im=memlist1[i])!=-1; i++ ) for( j=0; (jm=memlist2[j])!=-1; j++ )
+ smalldistmtx[i][j] = distmtx[MIN(im,jm)][MAX(im,jm)];
+
+ for( i=0; (im=memlist1[i])!=-1; i++ ) for( j=0; (jm=memlist2[j])!=-1; j++ )
+ reporterr( "old smalldistmtx[%d][%d] = %f\n", i, j, smalldistmtx[i][j] );
+if( im > 10 && jm > 10 ) exit( 1 );
+#endif
+ }
}
+
void makeEffMtx( int nseq, double **mtx, double *vec )
{
}
return( 0 );
}
-int msshrinklocalhom( char **pair, int s1, int s2, LocalHom **localhom, LocalHom ***localhomshrink )
+int msshrinklocalhom_fast_target( int *memlist1, int *memlist2, LocalHom **localhom, LocalHom ***localhomshrink, char *swaplist, int *targetmap )
{
- int m1, k1, n1, m2, k2, n2;
+ int m1, k1, m2, k2, t1, i2;
- for( m1=s1, k1=0; m1<njob; m1++ )
+ for( k1=0; (m1=memlist1[k1])!=-1; k1++ )
{
- if( pair[s1][m1] != 0 )
- {
- for( m2=s2, k2=0; m2<njob; m2++ )
+ if( targetmap[m1] == -1 )
+ {
+ swaplist[k1] = 1;
+// swaplist[k1] = 0; // DAME!!!
+ for( k2=0; (m2=memlist2[k2])!=-1; k2++ )
{
- if( pair[s2][m2] != 0 )
+ if( targetmap[m2] == -1 )
{
- 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++;
+ localhomshrink[k1][k2] = NULL;
+ continue;
}
+
+ t1 = targetmap[m2]; // start1 <-> start2, end1 <-> end2
+ i2 = m1;
+
+ if( localhom[t1][i2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[t1]+i2;
+ }
+ }
+ else
+ {
+ swaplist[k1] = 0;
+ for( k2=0; (m2=memlist2[k2])!=-1; k2++ )
+ {
+ t1 = targetmap[m1];
+ i2 = m2;
+
+ if( localhom[t1][i2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[t1]+i2;
}
- k1++;
+ }
+ }
+ return( 0 );
+}
+
+int msshrinklocalhom_fast_half( int *memlist1, int *memlist2, LocalHom **localhom, LocalHom ***localhomshrink )
+{
+ int m1, k1, m2, k2;
+
+ for( k1=0; (m1=memlist1[k1])!=-1; k1++ )
+ {
+ for( k2=0; (m2=memlist2[k2])!=-1; k2++ )
+ {
+ if( m1 < m2 )
+ {
+ if( localhom[m1][m2-m1].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[m1]+m2-m1;
+ }
+ else
+ {
+ if( localhom[m2][m1-m2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[m2]+m1-m2;
+ }
+ }
+ }
+ return( 0 );
+}
+
+int msshrinklocalhom_fast( int *memlist1, int *memlist2, LocalHom **localhom, LocalHom ***localhomshrink )
+{
+ int m1, k1, m2, k2;
+
+ for( k1=0; (m1=memlist1[k1])!=-1; k1++ )
+ {
+ for( k2=0; (m2=memlist2[k2])!=-1; k2++ )
+ {
+ if( localhom[m1][m2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[m1]+m2;
+ }
+ }
+ return( 0 );
+}
+int fastshrinklocalhom_one( int *mem1, int *mem2, int norg, 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( *intpt2 != norg )
+ {
+ fprintf( stderr, "ERROR! *intpt2 = %d\n", *intpt2 );
+ exit( 1 );
+ }
+ if( localhom[*intpt1][0].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[*intpt1];
}
}
return( 0 );
localhomshrink[k1][k2] = NULL;
else
localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2;
+
+// if( localhomshrink[k1][k2] != NULL )
+// printf( "ori localhomshrink[%d][%d].opt = %f\n", k1, k2, localhomshrink[k1][k2]->opt );
+ }
+ }
+ return( 0 );
+}
+
+
+int fastshrinklocalhom_half( 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( *intpt1 < *intpt2 )
+ {
+ if( localhom[*intpt1][*intpt2-*intpt1].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[*intpt1]+*intpt2-*intpt1;
+ }
+ else
+ {
+ if( localhom[*intpt2][*intpt1-*intpt2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[*intpt2]+*intpt1-*intpt2;
+ }
+
+// if( localhomshrink[k1][k2] != NULL )
+// printf( "ori localhomshrink[%d][%d].opt = %f, .importance=%f\n", k1, k2, localhomshrink[k1][k2]->opt, localhomshrink[k1][k2]->importance );
+ }
+ }
+ return( 0 );
+}
+
+
+int fastshrinklocalhom_target( int *mem1, int *mem2, LocalHom **localhom, LocalHom ***localhomshrink, char *swaplist, int *targetmap )
+{
+ int k1, k2;
+ int *intpt1, *intpt2;
+ int t1, i2;
+
+
+ for( intpt1=mem1, k1=0; *intpt1!=-1; intpt1++, k1++ )
+ {
+ if( targetmap[*intpt1] == -1 )
+ {
+ swaplist[k1] = 1;
+// swaplist[k1] = 0; // DAME!!!
+ for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
+ {
+ if( targetmap[*intpt2] == -1 )
+ {
+ localhomshrink[k1][k2] = NULL;
+ continue;
+ }
+
+ t1 = targetmap[*intpt2]; // end1<->end2, start1<->start2
+ i2 = *intpt1;
+
+ if( localhom[t1][i2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[t1]+i2;
+
+// if( localhomshrink[k1][k2] != NULL )
+// printf( "localhomshrink[%d][%d].opt = %f\n", k1, k2, localhomshrink[k1][k2]->opt );
+// else
+// printf( "localhomshrink[%d][%d] = NULL\n", k1, k2 );
+ }
+ }
+ else
+ {
+ swaplist[k1] = 0;
+ for( intpt2=mem2, k2=0; *intpt2!=-1; intpt2++, k2++ )
+ {
+ t1 = targetmap[*intpt1];
+ i2 = *intpt2;
+
+ if( localhom[t1][i2].opt == -1 )
+ localhomshrink[k1][k2] = NULL;
+ else
+ localhomshrink[k1][k2] = localhom[t1]+i2;
+
+// if( localhomshrink[k1][k2] != NULL )
+// printf( "localhomshrink[%d][%d].opt = %f\n", k1, k2, localhomshrink[k1][k2]->opt );
+// else
+// printf( "localhomshrink[%d][%d] = NULL\n", k1, k2 );
+ }
}
}
return( 0 );
git://source.jalview.org / jabaws.git / blobdiff