#include "mltaln.h"
+#define FULLSCORE 0
+
#define DEBUG 0
#define RECORD 0
static int nwa;
+
#ifdef enablemultithread
typedef struct _threadarg
{
int maxiter;
int nkozo;
int *subgenerationpt;
- float *basegainpt;
- float *gainlist;
- float *tscorelist;
+ double *basegainpt;
+ double *gainlist;
+ double *tscorelist;
int *generationofinput;
char *kozoarivec;
char **mastercopy;
int alloclen;
Node *stopol;
int ***topol;
-// double **len;
- float **tscorehistory_detail;
+ double **len;
+ double **tscorehistory_detail;
int *finishpt;
+ int **skipthisbranch;
+ double **distmtx;
+ int ntarget;
+ int *targetmap;
pthread_mutex_t *mutex;
pthread_cond_t *collection_end;
pthread_cond_t *collection_start;
}
#endif
+
+static void makescoringmatrices( double ***matrices, double **originalmtx )
+{
+ int c;
+ double rep;
+ for( c=0; c<maxdistclass; c++ )
+ {
+ rep = (double) 2 * c / ndistclass; // rep:0..2
+// fprintf( stderr, "rep = %f\n", rep );
+ makedynamicmtx( matrices[c], originalmtx, rep * 0.5 ); // upgma ni awaseru node, 0..1
+// fprintf( stderr, "c=%d, score for %c-%c = %f\n", c, 'W', 'W', matrices[c][amino_n['W']][amino_n['W']] );
+ }
+}
+
+static void classifypairs( int n1, double **eff1s, double *eff1, int n2, double **eff2s, double *eff2, double **smalldistmtx, int **matnum, int maxdistclass )
+{
+ int i, j, c;
+ for( c=0; c<maxdistclass; c++ )
+ {
+ for( i=0; i<n1; i++ ) eff1s[c][i] = 0.0;
+ for( j=0; j<n2; j++ ) eff2s[c][j] = 0.0;
+ }
+
+#if 0
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
+ {
+ c = (int)( smalldistmtx[i][j] / 2.0 * ndistclass ); // dist:0..2
+// if( c >= ndistclass ) c = ndistclass-1;
+ if( c >= maxdistclass ) c = maxdistclass-1;
+ fprintf( stderr, "pair %d-%d (%f), dist=%f -> c=%d\n", i, j, eff1[i] * eff2[j], smalldistmtx[i][j], c );
+ eff1s[c][i] += 1.0;
+ eff2s[c][j] += 1.0;
+ matnum[i][j] = c;
+ }
+ for( c=0; c<maxdistclass; c++ ) for( i=0; i<n1; i++ ) if(eff1s[c][i]) eff1s[c][i] = eff1[i]/eff1s[c][i];
+ for( c=0; c<maxdistclass; c++ ) for( i=0; i<n2; i++ ) if(eff2s[c][i]) eff2s[c][i] = eff2[i]/eff2s[c][i];
+#else
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
+ {
+ c = (int)( smalldistmtx[i][j] / 2.0 * ndistclass ); // dist:0..2
+// if( c >= ndistclass ) c = ndistclass-1;
+ if( c >= maxdistclass ) c = maxdistclass-1;
+// fprintf( stderr, "pair %d-%d (%f), dist=%f -> c=%d\n", i, j, eff1[i] * eff2[j], smalldistmtx[i][j], c );
+ eff1s[c][i] = eff1[i];
+ eff2s[c][j] = eff2[j];
+ matnum[i][j] = c;
+ }
+#endif
+#if 0
+ double totaleff;
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ ) reporterr( "whichmtx[%d][%d] = %d\n", i, j, matnum[i][j] );
+ for( c=0; c<maxdistclass; c++ )
+ {
+ totaleff = 0.0;
+ for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ ) totaleff += eff1s[c][i] * eff2s[c][j];
+ fprintf( stderr, "c=%d, sum totaleff1s-2s = %f\n", c, totaleff );
+ }
+ totaleff = 0.0; for( c=0; c<maxdistclass; c++ ) for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ ) totaleff += eff1s[c][i] * eff2s[c][j];
+ fprintf( stderr, "totaleff1s-2s = %f\n", totaleff );
+ totaleff = 0.0; for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ ) totaleff += eff1[i] * eff2[j];
+ fprintf( stderr, "totaleff1-2 = %f\n", totaleff );
+
+ totaleff = 0.0; for( c=0; c<maxdistclass; c++ ) for( i=0; i<n1; i++ ) totaleff += eff1s[c][i];
+ fprintf( stderr, "totaleff1s = %f\n", totaleff );
+ totaleff = 0.0; for( c=0; c<maxdistclass; c++ ) for( i=0; i<n2; i++ ) totaleff += eff2s[c][i];
+ fprintf( stderr, "totaleff2s = %f\n", totaleff );
+ totaleff = 0.0; for( i=0; i<n1; i++ ) totaleff += eff1[i];
+ fprintf( stderr, "totaleff1 = %f\n", totaleff );
+ totaleff = 0.0; for( i=0; i<n2; i++ ) totaleff += eff2[i];
+ fprintf( stderr, "totaleff2 = %f\n", totaleff );
+ for( c=0; c<maxdistclass; c++ )
+ {
+ for( i=0; i<n1; i++ ) fprintf( stderr, "eff1s[%d][%d] = %f\n", c, i, eff1s[c][i] );
+ for( i=0; i<n2; i++ ) fprintf( stderr, "eff2s[%d][%d] = %f\n", c, i, eff2s[c][i] );
+ fprintf( stderr, "\n" );
+ }
+exit( 1 );
+#endif
+}
+
static void Writeoption2( FILE *fp, int cycle, double cut )
{
fprintf( fp, "%dth cycle\n", cycle );
#ifdef enablemultithread
static void freelocalarrays(
- float *tscorehistory,
+ double *tscorehistory,
RNApair ***grouprna1, RNApair ***grouprna2,
RNApair *rnapairboth,
char *indication1, char *indication2,
+ double *distarr,
double *effarr, double *effarrforlocalhom, double *effarr1, double *effarr2,
char **mseq1, char **mseq2,
char **localcopy,
int *gapmap1, int *gapmap2,
double *effarr1_kozo, double *effarr2_kozo, double *effarr_kozo,
- char **pair,
- LocalHom *** localhomshrink
+ int **memlist,
+ char *pairbuf,
+ LocalHom *** localhomshrink,
+ char *swaplist,
+ double **smalldistmtx,
+ double ***scoringmatrices,
+ double **eff1s, double **eff2s,
+ int **whichmtx
)
{
// fprintf( stderr, "Skipping freelocalarrays\n" );
int i;
if( commonIP ) FreeIntMtx( commonIP );
commonIP = NULL;
- Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
- Falign_localhom( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL,NULL, 0, NULL );
+ Falign( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
+ Falign_localhom( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL,NULL, 0, NULL );
+ D__align( NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
+ partA__align_variousdist( NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL );
+
+
if( rnakozo && rnaprediction == 'm' )
{
free( grouprna1 ); // nakamimo?
free( grouprna2 ); // nakamimo?
}
- free( tscorehistory );
+ if( tscorehistory ) free( tscorehistory );
free( indication1 );
free( indication2 );
+ free( distarr );
free( effarr );
free( effarrforlocalhom );
free( effarr1 );
free( effarr2 );
free( mseq1 );
free( mseq2 );
- FreeCharMtx( localcopy );
+ if( localcopy ) FreeCharMtx( localcopy );
free( gapmap1 );
free( gapmap2 );
free( effarr2_kozo );
free( effarr_kozo );
- FreeCharMtx( pair );
+ FreeIntMtx( memlist );
+ free( pairbuf );
+
+ if( smalldistmtx ) FreeDoubleMtx( smalldistmtx );
+ if( scoringmatrices ) FreeDoubleCub( scoringmatrices );
+ if( eff1s ) FreeDoubleMtx( eff1s );
+ if( eff2s ) FreeDoubleMtx( eff2s );
+ if( whichmtx ) FreeIntMtx( whichmtx );
if( rnakozo ) free( rnapairboth );
free( localhomshrink[i] ); // nakamimo??
}
free( localhomshrink );
+ free( swaplist );
}
}
int *collectingpt = targ->collectingpt;
int *jobposintpt = targ->jobposintpt;
int nkozo = targ->nkozo;
- float *gainlist = targ->gainlist;
- float *tscorelist = targ->tscorelist;
+ double *gainlist = targ->gainlist;
+ double *tscorelist = targ->tscorelist;
int *generationofinput = targ->generationofinput;
int *subgenerationpt = targ->subgenerationpt;
- float *basegainpt = targ->basegainpt;
+ double *basegainpt = targ->basegainpt;
char *kozoarivec = targ->kozoarivec;
char **mastercopy = targ->mastercopy;
char ***candidates = targ->candidates;
int alloclen = targ->alloclen;
Node * stopol = targ->stopol;
int ***topol = targ->topol;
-// double **len = targ->len;
- float **tscorehistory_detail = targ->tscorehistory_detail;
+ double **len = targ->len;
+ double **tscorehistory_detail = targ->tscorehistory_detail;
int *finishpt = targ->finishpt;
+ int **skipthisbranch = targ->skipthisbranch;
+ double **distmtx = targ->distmtx;
+ int ntarget = targ->ntarget;
+ int *targetmap = targ->targetmap;
- int i, j, k, l, ii;
- float gain;
+ int i, k, l, ii;
+ double gain;
int iterate;
- char **pair;
+ int **memlist;
+ char *pairbuf;
int locnjob;
int s1, s2;
int clus1, clus2;
char **localcopy;
char **mseq1, **mseq2;
+ double *distarr; // re-calc
double *effarr, *effarr_kozo; // re-calc
double *effarr1, *effarr2, *effarr1_kozo, *effarr2_kozo;
char *indication1, *indication2;
RNApair ***grouprna1, ***grouprna2;
RNApair *rnapairboth;
LocalHom ***localhomshrink;
+ char *swaplist;
int *gapmap1, *gapmap2;
- float tscore, mscore, oimpmatch, impmatch;
+ double tscore, mscore;
+ double oimpmatchdouble;
+ double impmatchdouble;
int identity;
double tmpdouble;
- float naivescore0 = 0, naivescore1;
+// double naivescore0 = 0, naivescore1;
double *effarrforlocalhom;
- float *tscorehistory;
+ double *tscorehistory;
int intdum;
#if 0
int oscillating;
int lin, ldf;
#endif
- float maxgain;
+ double maxgain;
int bestthread;
int branchpos;
int subgenerationatfirst;
int myjob;
int converged2 = 0;
int chudanres;
+ double **smalldistmtx;
+ double ***scoringmatrices;
+ double **eff1s, **eff2s;
+ int **whichmtx;
locnjob = njob;
exit( 1 );
}
- tscorehistory = calloc( maxiter, sizeof( float ) );
+ tscorehistory = calloc( maxiter, sizeof( double ) );
if( rnakozo && rnaprediction == 'm' )
{
grouprna1 = grouprna2 = NULL;
}
- indication1 = AllocateCharVec( njob*3+50 );
- indication2 = AllocateCharVec( njob*3+50 );
+ indication1 = AllocateCharVec( 150 );
+ indication2 = AllocateCharVec( 150 );
+ distarr = AllocateDoubleVec( locnjob );
effarr = AllocateDoubleVec( locnjob );
effarrforlocalhom = AllocateDoubleVec( locnjob );
effarr1 = AllocateDoubleVec( locnjob );
localcopy = AllocateCharMtx( locnjob, alloclen );
gapmap1 = AllocateIntVec( alloclen );
gapmap2 = AllocateIntVec( alloclen );
+ if( specificityconsideration != 0 )
+ {
+ smalldistmtx = AllocateDoubleMtx( locnjob, locnjob ); // ookii?
+ scoringmatrices = AllocateDoubleCub( maxdistclass, nalphabets, nalphabets );
+ makescoringmatrices( scoringmatrices, n_dis_consweight_multi );
+ eff1s = AllocateDoubleMtx( maxdistclass, locnjob );
+ eff2s = AllocateDoubleMtx( maxdistclass, locnjob );
+ whichmtx = AllocateIntMtx( locnjob, locnjob );
+ }
+ else
+ {
+ smalldistmtx = NULL;
+ scoringmatrices = NULL;
+ eff1s = eff2s = NULL;
+ whichmtx = NULL;
+ }
effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
for( i=0; i<locnjob; i++ )
effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
- pair = AllocateCharMtx( locnjob, locnjob );
+ memlist = AllocateIntMtx( 2, locnjob );
+ pairbuf = AllocateCharVec( locnjob );
if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
+ swaplist = NULL;
if( constraint )
{
+ if( ntarget < locnjob ) swaplist = calloc( njob, sizeof( char ) );
localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
for( i=0; i<njob; i++)
{
}
}
-
if( thread_no == 0 )
{
*ntrypt = 0;
grouprna1, grouprna2,
rnapairboth,
indication1, indication2,
+ distarr,
effarr, effarrforlocalhom, effarr1, effarr2,
mseq1, mseq2,
localcopy,
gapmap1, gapmap2,
effarr1_kozo, effarr2_kozo, effarr_kozo,
- pair,
- localhomshrink
+ memlist, pairbuf,
+ localhomshrink,
+ swaplist,
+ smalldistmtx,
+ scoringmatrices,
+ eff1s, eff2s,
+ whichmtx
);
// return( NULL );
pthread_exit( NULL );
grouprna1, grouprna2,
rnapairboth,
indication1, indication2,
+ distarr,
effarr, effarrforlocalhom, effarr1, effarr2,
mseq1, mseq2,
localcopy,
gapmap1, gapmap2,
effarr1_kozo, effarr2_kozo, effarr_kozo,
- pair,
- localhomshrink
+ memlist, pairbuf,
+ localhomshrink,
+ swaplist,
+ smalldistmtx,
+ scoringmatrices,
+ eff1s, eff2s,
+ whichmtx
);
return( NULL );
pthread_exit( NULL );
grouprna1, grouprna2,
rnapairboth,
indication1, indication2,
+ distarr,
effarr, effarrforlocalhom, effarr1, effarr2,
mseq1, mseq2,
localcopy,
gapmap1, gapmap2,
effarr1_kozo, effarr2_kozo, effarr_kozo,
- pair,
- localhomshrink
+ memlist, pairbuf,
+ localhomshrink,
+ swaplist,
+ smalldistmtx,
+ scoringmatrices,
+ eff1s, eff2s,
+ whichmtx
);
return( NULL );
pthread_exit( NULL );
if( l == locnjob-2 ) k = 1;
else k = myjob - l * 2;
// fprintf( stderr, "JOB l=%d, k=%d\n", l, k );
+
+
branchpos = myjob;
(*jobposintpt)++;
iterate = *generationofmastercopypt;
(*ntrypt)++;
pthread_mutex_unlock( targ->mutex );
-
+// fprintf( stderr, "\n IRANAI IRANAI *jobposintpt=%d, nbranch = %d\n", *jobposintpt, nbranch );
// fprintf( stderr, "branchpos = %d (thread %d)\n", branchpos, thread_no );
fprintf( stderr, "STEP %03d-%03d-%d (Thread %d) ", iterate+1, l+1, k, thread_no );
fprintf( stderr, "STEP %03d-%03d-%d (thread %d) %s ", iterate+1, l+1, k, thread_no, use_fft?"\n":"\n" );
#endif
- for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
- OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
+// for( i=0; i<2; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
+ distFromABranch( locnjob, distarr, stopol, topol, len, l, k ); // ato de idou
+ OneClusterAndTheOther_fast( locnjob, memlist[0], memlist[1], &s1, &s2, pairbuf, topol, l, k, smalldistmtx, distmtx, distarr );
+
+
+// reporterr( "\n\n\n\n##### memlist[0][0], memlist[1][0] = %d, %d\n", memlist[0][0]+1, memlist[1][0]+1 );
+
+
+
#if 0
fprintf( stderr, "STEP%d-%d\n", l, k );
for( i=0; i<locnjob; i++ )
}
else if( weight == 4 )
{
+
weightFromABranch( locnjob, effarr, stopol, topol, l, k );
if( nkozo ) // hitomadu single weight.
{
if( nkozo )
{
- double tmptmptmp;
- tmptmptmp = 0.0;
- clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
+// double tmptmptmp;
+// tmptmptmp = 0.0;
+// clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair[0], s1, localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
+ clus1 = fastconjuction_noname_kozo( memlist[0], localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
- tmptmptmp = 0.0;
- clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
+// tmptmptmp = 0.0;
+// clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair[1], s2, localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
+ clus2 = fastconjuction_noname_kozo( memlist[1], localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
#if 0
}
else
{
- clus1 = conjuctionfortbfast( pair, s1, localcopy, mseq1, effarr1, effarr, indication1 );
- clus2 = conjuctionfortbfast( pair, s2, localcopy, mseq2, effarr2, effarr, indication2 );
+// clus1 = conjuctionfortbfast( pair[0], s1, localcopy, mseq1, effarr1, effarr, indication1 );
+// clus2 = conjuctionfortbfast( pair[1], s2, localcopy, mseq2, effarr2, effarr, indication2 );
+ clus1 = fastconjuction_noname( memlist[0], localcopy, mseq1, effarr1, effarr, indication1, minimumweight ); // 2015/Apr/18
+ clus2 = fastconjuction_noname( memlist[1], localcopy, mseq2, effarr2, effarr, indication2, minimumweight ); // 2015/Apr/18
}
if( rnakozo && rnaprediction == 'm' )
{
- makegrouprnait( grouprna1, singlerna, pair, s1 );
- makegrouprnait( grouprna2, singlerna, pair, s2 );
+// makegrouprnait( grouprna1, singlerna, pair[0], s1 );
+// makegrouprnait( grouprna2, singlerna, pair[1], s2 );
+ makegrouprna( grouprna1, singlerna, memlist[0] );
+ makegrouprna( grouprna2, singlerna, memlist[1] );
}
+ if( smalldistmtx )
+ {
+ classifypairs( clus1, eff1s, effarr1, clus2, eff2s, effarr2, smalldistmtx, whichmtx, maxdistclass );
+ }
+
+ if( alg == 'd' ) // Ichijiteki, koredeha scorecheck ga dekinai. gapmap ni taiou surumade
+ {
+ commongappick( clus1, mseq1 );
+ commongappick( clus2, mseq2 );
+ }
+
if( score_check == 2 )
{
fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
}
else if( score_check )
{
- if( RNAscoremtx == 'r' )
- intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
- else
- intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
-
if( constraint )
{
- shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
-// weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink ); // >>>
- oimpmatch = 0.0;
+ if( RNAscoremtx == 'r' )
+ intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+ else
+ intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+
+// shrinklocalhom( pair, int s1, int s2, localhomtable, localhomshrink );
+// msshrinklocalhom( pair[0], pair[1], s1, s2, localhomtable, localhomshrink );
+ if( ntarget < njob )
+ msshrinklocalhom_fast_target( memlist[0], memlist[1], localhomtable, localhomshrink, swaplist, targetmap ); // swaplist hitsuyou!!
+ else
+ msshrinklocalhom_fast_half( memlist[0], memlist[1], localhomtable, localhomshrink );
+ oimpmatchdouble = 0.0;
if( use_fft )
{
if( alg == 'Q' )
{
- part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
+ }
+ else if( alg == 'd' )
+ {
+ imp_match_init_strictD( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
+ if( rnakozo ) imp_rnaD( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
for( i=length-1; i>=0; i-- )
{
- oimpmatch += part_imp_match_out_scQ( i, i );
+ oimpmatchdouble += (double)imp_match_out_scD( i, i );
// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
}
}
else
{
- part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
for( i=length-1; i>=0; i-- )
{
- oimpmatch += part_imp_match_out_sc( i, i );
+ oimpmatchdouble += (double)part_imp_match_out_sc( i, i );
// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
}
}
{
if( alg == 'Q' )
{
- imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
-
- for( i=length-1; i>=0; i-- )
- {
- oimpmatch += imp_match_out_scQ( i, i );
-// fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
- }
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
}
else
{
- imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
fprintf( stderr, "not supported\n" );
exit( 1 );
for( i=length-1; i>=0; i-- )
{
- oimpmatch += imp_match_out_sc( i, i );
+ oimpmatchdouble += (double)imp_match_out_sc( i, i );
// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
}
}
}
else
{
- oimpmatch = 0.0;
+ if( RNAscoremtx == 'r' )
+ intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+ else
+ {
+ if( smalldistmtx )
+#if 1
+ intergroup_score_multimtx( whichmtx, scoringmatrices, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+#else
+ intergroup_score_dynmtx( smalldistmtx, amino_dis, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+#endif
+ else
+ intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+ }
+ oimpmatchdouble = 0.0;
}
// fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
- mscore = (double)oimpmatch + tmpdouble;
+ mscore = oimpmatchdouble + tmpdouble;
}
else
{
fprintf( stderr, "score_check = %d\n", score_check );
- fprintf( stderr, "Not supported\n" );
+ fprintf( stderr, "Not supported. Please add --threadit 0 to disable the multithreading in the iterative refinement calculation.\n" );
exit( 1 );
}
// if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
// if( !use_fft && !rnakozo )
+// if( !use_fft )
if( !use_fft )
{
commongappick_record( clus1, mseq1, gapmap1 );
for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
printf( "\n" );
#endif
- if( constraint == 2 )
+ if( !skipthisbranch[l][k] )
{
- if( use_fft )
- {
-// if( alg == 'Q' )
-// part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
-// else
-// part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
- chudanres = 0;
- Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, subgenerationpt, subgenerationatfirst, &chudanres );
-// fprintf( stderr, "##### impmatch = %f\n", impmatch );
- if( chudanres && parallelizationstrategy == BAATARI2 )
- {
-// fprintf( stderr, "#### yarinaoshi!!! INS-i\n" );
- goto yarinaoshi;
- }
- }
- else
- {
- if( alg == 'Q' )
- {
- float wm;
-// imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap wo tsukuttakara iranai.
-// if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
- wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
- fprintf( stderr, "wm = %f\n", wm );
-#if 0
- fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
- naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
- fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
-// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
-// rewind( stdout );
-// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
-// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
-// exit( 1 );
- }
-#endif
-#endif
- }
- else if( alg == 'R' )
- {
- float wm;
- imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
- wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
-// fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
- naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
-// fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
-// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
- rewind( stdout );
- for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- exit( 1 );
- }
-#endif
- }
- else if( alg == 'H' )
+ if( constraint == 2 )
+ {
+ if( use_fft )
{
- float wm;
- imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
- wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
- fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
- naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
- fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
+// if( alg == 'Q' )
+// part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
+// else
+// part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ chudanres = 0;
+ if( alg == 'd' )
+ D__align( n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatchdouble, NULL, NULL, NULL, NULL, subgenerationpt, subgenerationatfirst, &chudanres, 1, 1 );
+
else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
+ Falign_localhom( whichmtx, scoringmatrices, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, localhomshrink, &impmatchdouble, gapmap1, gapmap2, subgenerationpt, subgenerationatfirst, &chudanres );
+// fprintf( stderr, "##### impmatch = %f\n", impmatch );
+ if( chudanres && parallelizationstrategy == BAATARI2 )
{
-// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
-// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
-// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
-// exit( 1 );
+// fprintf( stderr, "#### yarinaoshi!!! INS-i\n" );
+ goto yarinaoshi;
}
-#endif
}
else
{
-// imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
- A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
- fprintf( stderr, "A__align_gapmap\n" );
-// fprintf( stderr, "##### impmatch = %f\n", impmatch );
+ fprintf( stderr, "Not supported\n" );
+ exit( 1 );
}
}
- }
- else if( use_fft )
- {
- float totalwm;
- chudanres = 0;
- totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, subgenerationpt, subgenerationatfirst, &chudanres );
- if( chudanres && parallelizationstrategy == BAATARI2 )
+ else if( use_fft )
{
-// fprintf( stderr, "#### yarinaoshi!!! FFT-NS-i\n" );
- goto yarinaoshi;
- }
-
-// fprintf( stderr, "totalwm = %f\n", totalwm );
#if 0
- if( alg == 'Q' )
- {
- fprintf( stderr, "totalwm = %f\n", totalwm );
- naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 1 // chuui
- if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
+ if( alg == 'd' )
{
-// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
-// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
-// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
-// exit( 1 );
+ D__align_variousdist( whichmtx, scoringmatrices, NULL, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, NULL, &dumdb, NULL, NULL, NULL, NULL, subgenerationpt, subgenerationatfirst, &chudanres, 1, 1 );
}
-#endif
- }
-#endif
- if( alg == 'R' )
- {
- fprintf( stderr, "totalwm = %f\n", totalwm );
- naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 1 // chuui
- if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
+#endif
{
-// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
-// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
-// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
-// exit( 1 );
+ double totalwm;
+ chudanres = 0;
+// totalwm = Falign( smalldistmtx, NULL, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, NULL, NULL, clus1, clus2, alloclen, &intdum, subgenerationpt, subgenerationatfirst, &chudanres );
+ totalwm = Falign( whichmtx, scoringmatrices, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, &intdum, subgenerationpt, subgenerationatfirst, &chudanres );
}
- }
-#endif
- }
- else
- {
- if( alg == 'M' )
- {
- chudanres = 0;
- MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, subgenerationpt, subgenerationatfirst, &chudanres, outgap, outgap );
if( chudanres && parallelizationstrategy == BAATARI2 )
{
-// fprintf( stderr, "#### yarinaoshi!!! NW-NS-i\n" );
+// fprintf( stderr, "#### yarinaoshi!!! FFT-NS-i\n" );
goto yarinaoshi;
}
+
}
- else if( alg == 'A' )
+ else
{
- A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); //outgap==1
+ fprintf( stderr, "\n\nUnexpected error. Please contact kazutaka.katoh@aist.go.jp\n\n\n" );
+ exit( 1 );
}
- else if( alg == 'Q' )
+// fprintf( stderr, "## impmatch = %f\n", impmatch );
+
+#if 1
+ if( parallelizationstrategy == BAATARI2 && *subgenerationpt != subgenerationatfirst )
{
- float wm;
- wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
- fprintf( stderr, "wm = %f\n", wm );
- fprintf( stderr, "impmatch = %f\n", impmatch );
- naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 1 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
-// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
-// rewind( stderr );
-// rewind( stdout );
-// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
-// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
-// exit( 1 );
- }
-#endif
+// fprintf( stderr, "\nYarinaoshi2!! (Thread %d)\n", thread_no );
+ goto yarinaoshi;
}
- else if( alg == 'R' )
- {
- float wm;
- wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
- naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-#if 1 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
-// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
-// rewind( stderr );
-// rewind( stdout );
-// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
-// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
-// exit( 1 );
- }
#endif
- }
- else if( alg == 'H' )
- {
- float wm;
- wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
- naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- {
- fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
- }
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
+
+ identity = !strcmp( localcopy[s1], mastercopy[s1] );
+ identity *= !strcmp( localcopy[s2], mastercopy[s2] );
+ fprintf( stderr, "%03d-%04d-%d (thread %4d) identical \r", iterate+1, *ndonept, k, thread_no );
-#if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
- rewind( stdout );
- for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- exit( 1 );
- }
-#endif
- }
- else if( alg == 'a' )
- {
- Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
- }
- else ErrorExit( "Sorry!" );
}
-// fprintf( stderr, "## impmatch = %f\n", impmatch );
-
-#if 1
- if( parallelizationstrategy == BAATARI2 && *subgenerationpt != subgenerationatfirst )
+ else
{
-// fprintf( stderr, "\nYarinaoshi2!! (Thread %d)\n", thread_no );
- goto yarinaoshi;
+ identity = 1;
+ fprintf( stderr, "%03d-%04d-%d (thread %4d) skip \r", iterate+1, *ndonept, k, thread_no );
}
-#endif
-
- identity = !strcmp( localcopy[s1], mastercopy[s1] );
- identity *= !strcmp( localcopy[s2], mastercopy[s2] );
/* Bug? : idnetitcal but score change when scoreing mtx != JTT */
if( identity )
{
tscore = mscore;
-// if( !devide ) fprintf( trap_g, "tscore = %f identical.\n", tscore );
-// fprintf( stderr, " identical." );
- fprintf( stderr, "%03d-%04d-%d (thread %4d) identical \r", iterate+1, *ndonept, k, thread_no );
}
else
{
intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
#endif
- tscore = impmatch + tmpdouble;
+ tscore = impmatchdouble + tmpdouble;
// fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
}
else
{
- intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+ if( smalldistmtx )
+#if 1
+ intergroup_score_multimtx( whichmtx, scoringmatrices, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+#else
+ intergroup_score_dynmtx( smalldistmtx, amino_dis, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+#endif
+ else
+ intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
tscore = tmpdouble;
}
// fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
goto yarinaoshi;
}
#endif
+// reporterr( "tscore = %f, mscore = %f\n", tscore, mscore );
gain = tscore - ( mscore - cut/100.0*mscore );
if( gain > 0 )
{
tscore = mscore;
}
}
+#if FULLSCORE
+ {
+ int j;
+ double fullscore = 0.0;
+ for( i=1; i<locnjob; i++ ) for( j=0; j<i; j++ )
+ fullscore += (double)naivepairscore11( localcopy[i], localcopy[j], penalty );
+ reporterr( "\n######## fullscore = %f\n", fullscore / (locnjob*(locnjob-1)/2) );
+ }
+#endif
+
+
converged2 = 0;
for( ii=iterate-2; ii>=0; ii-=1 )
{
int TreeDependentIteration( int locnjob, char **name, int nlen[M],
char **aseq, char **bseq, int ***topol, double **len,
+ double **distmtx,
+ int **skipthisbranch,
int alloclen, LocalHom **localhomtable,
RNApair ***singlerna,
- int nkozo, char *kozoarivec )
+ int nkozo, char *kozoarivec,
+ int ntarget, int *targetmap, int *targetmapr )
{
int i, j, k, l, iterate, ii, iu, ju;
int lin, ldf, length;
int s1, s2;
static double **imanoten;
static Node *stopol;
+ static double *distarr = NULL;
static double *effarrforlocalhom = NULL;
static double *effarr = NULL;
static double *effarr1 = NULL;
static int *branchnode = NULL;
static double **branchWeight = NULL;
static char **mseq1, **mseq2;
- static float ***history;
+ static double ***history;
FILE *trap;
double tscore, mscore;
int identity;
int converged;
int oscillating;
- float naivescore0 = 0.0; // by D.Mathog, a guess
- float naivescore1;
+// double naivescore0 = 0.0; // by D.Mathog, a guess
+// double naivescore1;
#if 0
char pair[njob][njob];
#else
- static char **pair;
+ static int **memlist;
+ static char *pairbuf;
#endif
#if DEBUG + RECORD
double score_for_check0, score_for_check1;
#endif
static char *indication1, *indication2;
static LocalHom ***localhomshrink = NULL;
- float impmatch = 0.0, oimpmatch = 0.0;
+ static char *swaplist = NULL;
+ double impmatchdouble = 0.0;
+ double oimpmatchdouble = 0.0;
static int *gapmap1;
static int *gapmap2;
double tmpdouble;
static RNApair *rnapairboth;
RNApair ***grouprna1, ***grouprna2;
double unweightedspscore;
+ static double **smalldistmtx;
+ static double ***scoringmatrices;
+ static double **eff1s, **eff2s;
+ static int **whichmtx;
+ int value;
if( rnakozo && rnaprediction == 'm' )
{
if( effarr == NULL ) /* locnjob == njob ni kagiru */
{
- indication1 = AllocateCharVec( njob*3+50 );
- indication2 = AllocateCharVec( njob*3+50 );
+ indication1 = AllocateCharVec( 150 );
+ indication2 = AllocateCharVec( 150 );
effarr = AllocateDoubleVec( locnjob );
+ distarr = AllocateDoubleVec( locnjob );
effarrforlocalhom = AllocateDoubleVec( locnjob );
effarr1 = AllocateDoubleVec( locnjob );
effarr2 = AllocateDoubleVec( locnjob );
gapmap1 = AllocateIntVec( alloclen );
gapmap2 = AllocateIntVec( alloclen );
if( score_check == 2 ) imanoten = AllocateDoubleMtx( njob, njob );
+ if( specificityconsideration != 0 )
+ {
+ smalldistmtx = AllocateDoubleMtx( locnjob, locnjob ); // ookii?
+ scoringmatrices = AllocateDoubleCub( maxdistclass, nalphabets, nalphabets );
+ makescoringmatrices( scoringmatrices, n_dis_consweight_multi );
+ eff1s = AllocateDoubleMtx( maxdistclass, locnjob );
+ eff2s = AllocateDoubleMtx( maxdistclass, locnjob );
+ whichmtx = AllocateIntMtx( locnjob, locnjob );
+ }
+ else
+ {
+ smalldistmtx = NULL;
+ scoringmatrices = NULL;
+ eff1s = eff2s = NULL;
+ whichmtx = NULL;
+ }
effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
#if 0
#else
- pair = AllocateCharMtx( locnjob, locnjob );
+ pairbuf = AllocateCharVec( locnjob );
+ memlist = AllocateIntMtx( 2, locnjob );
if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
+ swaplist = NULL;
if( constraint )
{
+ if( ntarget < locnjob ) swaplist = calloc( njob, sizeof( char ) );
localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
for( i=0; i<njob; i++)
{
if( constraint )
{
counteff_simple( locnjob, topol, len, effarrforlocalhom );
- calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
+ if( ntarget < locnjob )
+ calcimportance_target( locnjob, ntarget, effarrforlocalhom, aseq, localhomtable, targetmap, targetmapr );
+ else
+ calcimportance_half( locnjob, effarrforlocalhom, aseq, localhomtable );
}
if( weight == 2 )
fprintf( stderr, "Not supported, weight=%d nkozo=%d.\n", weight, nkozo );
}
}
- else if( weight == 4 )
+// else if( weight == 4 )
+// else if( weight == 4 || weight == 0 )
+ else if( locnjob > 2 && ( weight == 4 || weight == 0 ) )
{
treeCnv( stopol, locnjob, topol, len, branchWeight );
- calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
+ calcBranchWeight( branchWeight, locnjob, stopol, topol, len ); // IRU!!!
}
}
int jobposint;
int generationofmastercopy;
int subgeneration;
- float basegain;
+ double basegain;
int *generationofinput;
- float *gainlist;
- float *tscorelist;
+ double *gainlist;
+ double *tscorelist;
int ndone;
int ntry;
int collecting;
int maxiter;
char ***candidates;
int *branchtable;
- float **tscorehistory_detail;
+ double **tscorehistory_detail;
int finish;
nwa = nthread + 1;
pthread_cond_init( &collection_end, NULL );
pthread_cond_init( &collection_start, NULL );
- gainlist = calloc( nwa, sizeof( float ) );
- tscorelist = calloc( nwa, sizeof( float ) );
+ gainlist = calloc( nwa, sizeof( double ) );
+ tscorelist = calloc( nwa, sizeof( double ) );
branchtable = calloc( nbranch, sizeof( int ) );
generationofinput = calloc( nbranch, sizeof( int ) );
if( parallelizationstrategy == BESTFIRST )
targ[i].alloclen = alloclen;
targ[i].stopol = stopol;
targ[i].topol = topol;
-// targ[i].len = len;
+ targ[i].skipthisbranch = skipthisbranch;
+ targ[i].distmtx = distmtx;
+ targ[i].len = len;
targ[i].mutex = &mutex;
targ[i].collection_end = &collection_end;
targ[i].collection_start = &collection_start;
targ[i].tscorehistory_detail = tscorehistory_detail;
+ targ[i].ntarget = ntarget;
+ targ[i].targetmap = targetmap;
targ[i].finishpt = &finish;
pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
if( weight == 2 )
countnode_int( locnjob, topol, node );
- else if( weight == 4 )
+// else if( weight == 4 )
+// else if( weight == 4 || weight == 0 )
+ else if( locnjob > 2 && ( weight == 4 || weight == 0 ) )
{
treeCnv( stopol, locnjob, topol, len, branchWeight );
- calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
+ calcBranchWeight( branchWeight, locnjob, stopol, topol, len ); // IRU!!!
}
trap = fopen( "hat2", "w" );
if( !trap ) ErrorExit( "Cannot open hat2." );
if( constraint )
{
counteff_simple( locnjob, topol, len, effarrforlocalhom );
- calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
+ if( ntarget < locnjob )
+ calcimportance_target( locnjob, ntarget, effarrforlocalhom, aseq, localhomtable, targetmap, targetmapr );
+ else
+ calcimportance_half( locnjob, effarrforlocalhom, aseq, localhomtable );
}
}
if( l == locnjob-2 ) k = 1;
else k = myjob - l * 2;
#endif
+ #if 0 // IRANAI!!!!
+ fprintf( stderr, "\nSTEP %d-%d\n", l, k );
+ for( i=0; topol[l][k][i]!=-1; i++ ) fprintf( stderr, " %d ", topol[l][k][i]+1 );
+ fprintf( stderr, "\n" );
+ fprintf( stderr, "SKIP %d\n", skipthisbranch[l][k] );
+ #endif
#if 1
fprintf( stderr, "STEP %03d-%03d-%d ", iterate+1, l+1, k );
fflush( stderr );
#else
fprintf( stderr, "STEP %03d-%03d-%d %s", iterate+1, l+1, k, use_fft?"\n":"\n" );
#endif
- for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
-
- OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
- #if 0
+ if( skipthisbranch[l][k] )
+ {
+ fprintf( stderr, " skip. \r" );
+ continue;
+ }
+
+ distFromABranch( locnjob, distarr, stopol, topol, len, l, k ); // ato de dou. weight 4 igai demo tsukaeruyou ni.
+// for( i=0; i<locnjob; i++ ) // BUG??
+// for( i=0; i<2; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
+ OneClusterAndTheOther_fast( locnjob, memlist[0], memlist[1], &s1, &s2, pairbuf, topol, l, k, smalldistmtx, distmtx, distarr );
+
+
+
+
+// reporterr( "\n\n##### memlist[0][0] = %d\n", memlist[0][0]+1 );
+// reporterr( "##### memlist[1][0] = %d\n\n", memlist[1][0]+1 );
+
+
+
+ #if 0 // IRANAI!!!!
fprintf( stderr, "STEP%d-%d\n", l, k );
- for( i=0; i<locnjob; i++ )
+ for( i=0; i<2; i++ )
{
for( j=0; j<locnjob; j++ )
{
}
else if( weight == 4 )
{
+
weightFromABranch( locnjob, effarr, stopol, topol, l, k );
+
#if 0
if( nkozo )
{
if( nkozo )
{
- #if 1
- double tmptmptmp;
- tmptmptmp = 0.0;
- clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
+// double tmptmptmp;
+// tmptmptmp = 0.0;
+// clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair[0], s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
+ clus1 = fastconjuction_noname_kozo( memlist[0], aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
- tmptmptmp = 0.0;
- clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
+// tmptmptmp = 0.0;
+// clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair[1], s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
+ clus2 = fastconjuction_noname_kozo( memlist[1], aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
#if 0
for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[] = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
#endif
-
-
-
-
- #else
- clus1 = conjuctionfortbfast_kozo_BUG( pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
- clus2 = conjuctionfortbfast_kozo_BUG( pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
- #endif
}
else
{
- clus1 = conjuctionfortbfast( pair, s1, aseq, mseq1, effarr1, effarr, indication1 );
- clus2 = conjuctionfortbfast( pair, s2, aseq, mseq2, effarr2, effarr, indication2 );
+// clus1 = conjuctionfortbfast( pair[0], s1, aseq, mseq1, effarr1, effarr, indication1 );
+// clus2 = conjuctionfortbfast( pair[1], s2, aseq, mseq2, effarr2, effarr, indication2 );
+ clus1 = fastconjuction_noname( memlist[0], aseq, mseq1, effarr1, effarr, indication1, minimumweight ); // 2015/Apr/18
+ clus2 = fastconjuction_noname( memlist[1], aseq, mseq2, effarr2, effarr, indication2, minimumweight ); // 2015/Apr/18
}
-
if( rnakozo && rnaprediction == 'm' )
{
- makegrouprnait( grouprna1, singlerna, pair, s1 );
- makegrouprnait( grouprna2, singlerna, pair, s2 );
+// makegrouprnait( grouprna1, singlerna, pair[0], s1 );
+// makegrouprnait( grouprna2, singlerna, pair[1], s2 );
+ makegrouprna( grouprna1, singlerna, memlist[0] );
+ makegrouprna( grouprna2, singlerna, memlist[1] );
}
+ if( smalldistmtx )
+ {
+ classifypairs( clus1, eff1s, effarr1, clus2, eff2s, effarr2, smalldistmtx, whichmtx, maxdistclass );
+ }
+
+ if( alg == 'd' ) // Ichijiteki, koredeha scorecheck ga dekinai. gapmap ni taiou surumade
+ {
+ commongappick( clus1, mseq1 );
+ commongappick( clus2, mseq2 );
+ }
+
+
+
if( score_check == 2 )
{
+ fprintf( stderr, "Not supported\n" );
+ exit( 1 );
if( constraint )
{
- // msshrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
- shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
- oimpmatch = 0.0;
+// msshrinklocalhom( pair[0], pair[1], s1, s2, localhomtable, localhomshrink );
+ if( ntarget < locnjob )
+ msshrinklocalhom_fast_target( memlist[0], memlist[1], localhomtable, localhomshrink, swaplist, targetmap );
+ else
+ msshrinklocalhom_fast_half( memlist[0], memlist[1], localhomtable, localhomshrink );
+ oimpmatchdouble = 0.0;
if( use_fft )
{
if( alg == 'Q' )
{
- part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
- for( i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_scQ( i, i );
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
}
else
{
- part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
- for( i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_sc( i, i );
+ for( i=length-1; i>=0; i-- ) oimpmatchdouble += (double)part_imp_match_out_sc( i, i );
}
}
else
{
if( alg == 'Q' )
{
- imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
- for( i=length-1; i>=0; i-- ) oimpmatch += imp_match_out_scQ( i, i );
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
}
else
{
- imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
fprintf( stderr, "not supported\n" );
exit( 1 );
}
}
else
{
- oimpmatch = 0.0;
+ oimpmatchdouble = 0.0;
}
+#if 0
tmpdouble = 0.0;
iu=0;
for( i=s1; i<locnjob; i++ )
{
- if( !pair[s1][i] ) continue;
+ if( !pair[0][i] ) continue;
ju=0;
for( j=s2; j<locnjob; j++ )
{
- if( !pair[s2][j] ) continue;
+ if( !pair[1][j] ) continue;
// fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(i,j)][MAX(i,j)];
ju++;
}
iu++;
}
- mscore = oimpmatch + tmpdouble;
+#else // not yet checked
+ fprintf( stderr, "##### NOT YET CHECKED!!!!\n" );
+ exit( 1 );
+ tmpdouble = 0.0;
+ iu=0;
+ for( i=0; (s1=memlist[0][i])!=-1; i++ )
+ {
+ ju=0;
+ for( j=0; (s2=memlist[1][j])!=-1; j++ )
+ {
+ // fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
+ tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(s1,s2)][MAX(s1,s2)];
+ ju++;
+ }
+ iu++;
+ }
+#endif
+ mscore = oimpmatchdouble + tmpdouble;
}
else if( score_check )
{
- #if 1
- if( RNAscoremtx == 'r' )
- intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
- else
- intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
- #else
- intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
- #endif
-
if( constraint )
{
- shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
+ if( RNAscoremtx == 'r' )
+ intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+ else
+ intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+
+// shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
+ if( ntarget < locnjob )
+ msshrinklocalhom_fast_target( memlist[0], memlist[1], localhomtable, localhomshrink, swaplist, targetmap );
+ else
+ msshrinklocalhom_fast_half( memlist[0], memlist[1], localhomtable, localhomshrink );
// weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink ); // >>>
- oimpmatch = 0.0;
+ oimpmatchdouble = 0.0;
if( use_fft )
{
if( alg == 'Q' )
{
- part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
+ }
+ else if( alg == 'd' )
+ {
+ imp_match_init_strictD( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
+
for( i=length-1; i>=0; i-- )
{
- oimpmatch += part_imp_match_out_scQ( i, i );
- // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
+ oimpmatchdouble += (double)imp_match_out_scD( i, i );
+ // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
}
}
else
{
- part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
for( i=length-1; i>=0; i-- )
{
- oimpmatch += part_imp_match_out_sc( i, i );
+ oimpmatchdouble += (double)part_imp_match_out_sc( i, i );
// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
}
}
{
if( alg == 'Q' )
{
- imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
-
- for( i=length-1; i>=0; i-- )
- {
- oimpmatch += imp_match_out_scQ( i, i );
- // fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
- }
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
}
else
{
- imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
fprintf( stderr, "not supported\n" );
exit( 1 );
for( i=length-1; i>=0; i-- )
{
- oimpmatch += imp_match_out_sc( i, i );
+ oimpmatchdouble += (double)imp_match_out_sc( i, i );
// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
}
}
}
else
{
- oimpmatch = 0.0;
+ if( RNAscoremtx == 'r' )
+ intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+ else
+ {
+ if( smalldistmtx )
+#if 1
+ intergroup_score_multimtx( whichmtx, scoringmatrices, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+#else
+ intergroup_score_dynmtx( offsetmtx, amino_dis, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // n_dis ha machigai
+#endif
+ else
+ intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
+ }
+ oimpmatchdouble = 0.0;
}
// fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
- mscore = (double)oimpmatch + tmpdouble;
+ mscore = oimpmatchdouble + tmpdouble;
}
else
{
- // fprintf( stderr, "score_check = %d\n", score_check );
- /* atode kousokuka */
+// fprintf( stderr, "score_check = %d\n" );
+#if 1
+ /* Oscilation check no tame hitsuyou! atode kousokuka */
intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
mscore = tmpdouble;
/* atode kousokuka */
+#else
+ mscore = 0.0;
+#endif
if( constraint )
{
- oimpmatch = 0.0;
- shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
+ oimpmatchdouble = 0.0;
+// shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
+ if( ntarget < locnjob )
+ msshrinklocalhom_fast_target( memlist[0], memlist[1], localhomtable, localhomshrink, swaplist, targetmap );
+ else
+ msshrinklocalhom_fast_half( memlist[0], memlist[1], localhomtable, localhomshrink );
if( use_fft )
{
if( alg == 'Q' )
{
- part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
+ }
+ else if( alg == 'd' )
+ {
+ imp_match_init_strictD( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
+ if( rnakozo ) imp_rnaD( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
}
else
{
- part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
}
}
{
if( alg == 'Q' )
{
- imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
- if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
+ fprintf( stderr, "'Q' is no longer supported\n" );
+ exit( 1 );
}
else
{
- imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
+ imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, swaplist, 1, memlist[0], memlist[1] );
fprintf( stderr, "Not supported\n" );
exit( 1 );
}
// if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
// if( !use_fft && !rnakozo )
+// if( !use_fft )
if( !use_fft )
{
commongappick_record( clus1, mseq1, gapmap1 );
// part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
// else
// part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
- Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, NULL, 0, NULL );
- // fprintf( stderr, "##### impmatch = %f\n", impmatch );
- }
- else
- {
- if( alg == 'Q' )
+ if( alg == 'd' )
{
- float wm;
- // imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap wo tsukuttakara iranai.
- // if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
-
- wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
- fprintf( stderr, "wm = %f\n", wm );
- #if 0
- fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
- naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
- fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
+ if( scoringmatrices ) // called by tditeration.c
{
- // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
- // rewind( stdout );
- // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- // exit( 1 );
+ D__align_variousdist( whichmtx, scoringmatrices, NULL, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, localhomshrink, &impmatchdouble, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 );
}
- #endif
- #endif
- }
- else if( alg == 'R' )
- {
- float wm;
- imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
- wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
- // fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
- naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
- // fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
{
- // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
- rewind( stdout );
- for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- exit( 1 );
+ D__align( n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatchdouble, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 );
}
- #endif
- }
- else if( alg == 'H' )
- {
- float wm;
- imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
- wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
- fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
- naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
- fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
- // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
- // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- // exit( 1 );
- }
- #endif
}
else
- {
- // imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
- A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
- // fprintf( stderr, "##### impmatch = %f\n", impmatch );
- }
+ Falign_localhom( whichmtx, scoringmatrices, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, localhomshrink, &impmatchdouble, gapmap1, gapmap2, NULL, 0, NULL );
+ // fprintf( stderr, "##### impmatch = %f\n", impmatch );
+ }
+ else
+ {
+ fprintf( stderr, "\n\nUnexpected error. Please contact kazutaka.katoh@aist.go.jp\n\n\n" );
+ exit( 1 );
}
}
else if( use_fft )
{
- float totalwm;
- totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
-
- // fprintf( stderr, "totalwm = %f\n", totalwm );
- #if 0
- if( alg == 'Q' )
+ double totalwm;
+#if 0
+ double dumdb;
+ // D ha Falign wo toshite yobareru.
+ if( alg == 'd' )
{
- fprintf( stderr, "totalwm = %f\n", totalwm );
- naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 1 // chuui
- if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
- {
- // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
- // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- // exit( 1 );
- }
- #endif
+ D__align_variousdist( whichmtx, scoringmatrices, NULL, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, NULL, &dumdb, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 );
}
- #endif
- if( alg == 'R' )
+ else
+#endif
{
- fprintf( stderr, "totalwm = %f\n", totalwm );
- naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 1 // chuui
- if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
- {
- // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
- // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- // exit( 1 );
- }
+ totalwm = Falign( whichmtx, scoringmatrices, n_dis_consweight_multi, mseq1, mseq2, effarr1, effarr2, eff1s, eff2s, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
+
+ // fprintf( stderr, "totalwm = %f\n", totalwm );
}
- #endif
}
else
{
- if( alg == 'M' )
- {
- MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
- }
- else if( alg == 'A' )
- {
- A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); // outgap==1
- }
- else if( alg == 'Q' )
- {
- float wm;
- wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
- fprintf( stderr, "wm = %f\n", wm );
- fprintf( stderr, "impmatch = %f\n", impmatch );
- naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 1 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
- // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
- // rewind( stderr );
- // rewind( stdout );
- // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- // exit( 1 );
- }
- #endif
- }
- else if( alg == 'R' )
- {
- float wm;
- wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
- naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
- #if 1 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
- // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
- // rewind( stderr );
- // rewind( stdout );
- // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- // exit( 1 );
- }
- #endif
- }
- else if( alg == 'H' )
- {
- float wm;
- wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
- naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
-
- if( naivescore1 > naivescore0 )
- fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
- else if( naivescore1 < naivescore0 )
- {
- fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
- }
- else
- fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
-
- #if 0 // chuui
- if( abs( wm - naivescore1 ) > 100 )
- {
- rewind( stdout );
- for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
- for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
- exit( 1 );
- }
- #endif
- }
- else if( alg == 'a' )
- {
- Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
- }
- else ErrorExit( "Sorry!" );
+ fprintf( stderr, "\n\nUnexpected error. Please contact kazutaka.katoh@aist.go.jp\n\n\n" );
+ exit( 1 );
}
// fprintf( stderr, "## impmatch = %f\n", impmatch );
/* Bug? : idnetitcal but score change when scoreing mtx != JTT */
+#if FULLSCORE
+ {
+ double fullscore = 0.0, fullscoreori = 0.0;
+ for( i=1; i<locnjob; i++ ) for( j=0; j<i; j++ )
+ fullscore += (double)naivepairscore11( aseq[i], aseq[j], penalty );
+ reporterr( "\nfullscore = %f\n", fullscore / (locnjob*(locnjob-1)/2) );
+ for( i=1; i<locnjob; i++ ) for( j=0; j<i; j++ )
+ fullscoreori += (double)naivepairscore11( bseq[i], bseq[j], penalty );
+ reporterr( "\nfullscoreori = %f\n", fullscoreori / (locnjob*(locnjob-1)/2) );
+ if( 0 && fullscoreori > fullscore )
+ {
+ for( i=0; i<clus1; i++ )
+ fprintf( stdout, ">group1\n%s\n", mseq1[i] );
+ for( i=0; i<clus2; i++ )
+ fprintf( stdout, ">group2\n%s\n", mseq2[i] );
+
+
+ for( i=0; i<locnjob; i++ )
+ fprintf( stdout, ">better alignment\n%s\n", bseq[i] );
+ exit( 1 );
+ }
+ }
+#endif
length = strlen( mseq1[0] );
{
tscore = mscore;
if( !devide ) fprintf( trap_g, "tscore = %f identical.\n", tscore );
- fprintf( stderr, " identical." );
+ fprintf( stderr, " identical. " );
converged++;
}
else
intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
#endif
- tscore = impmatch + tmpdouble;
+ tscore = impmatchdouble + tmpdouble;
// fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
}
else
{
- intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+ if( smalldistmtx )
+#if 1
+ intergroup_score_multimtx( whichmtx, scoringmatrices, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
+#else
+ intergroup_score_dynmtx( offsetmtx, amino_dis, mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // n_dis ha machigai
+#endif
+ else
+ intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
tscore = tmpdouble;
}
// fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
converged++;
}
}
- fprintf( stderr, "\r" );
+ if( alg == 'd' )
+ fprintf( stderr, "\n" );
+ else
+ fprintf( stderr, "\r" );
- history[iterate][l][k] = (float)tscore;
+ history[iterate][l][k] = (double)tscore;
// fprintf( stderr, "tscore = %f\n", tscore );
if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
fprintf( stderr, "\n\n" );
}
- return( 0 );
+ value = 0;
+ goto end;
}
if( iterate >= 1 )
{
oscillating = 0;
for( ii=iterate-2; ii>=0; ii-=2 )
{
- if( (float)tscore == history[ii][l][k] )
+ if( (double)tscore == history[ii][l][k] )
{
oscillating = 1;
break;
fprintf( stderr, "\n\n" );
}
#if 1 /* hujuubun */
- return( -1 );
+ value = -1;
+ goto end;
#endif
}
} /* if( iterate ) */
}
} /* for( iterate ) */
}
- return( 2 );
+ value = 2;
+
+ end:
+ if( grouprna1 ) free( grouprna1 );
+ if( grouprna2 ) free( grouprna2 );
+// freelocalarrays
+// (
+// NULL,
+// grouprna1, grouprna2,
+// rnapairboth,
+// indication1, indication2,
+// distarr,
+// effarr, effarrforlocalhom, effarr1, effarr2,
+// mseq1, mseq2,
+// NULL,
+// gapmap1, gapmap2,
+// effarr1_kozo, effarr2_kozo, effarr_kozo,
+// memlist, pairbuf,
+// localhomshrink,
+// smalldistmtx,
+// scoringmatrices,
+// eff1s, eff2s,
+// whichmtx
+// );
+// free( branchnode );
+// free( stopol );
+ return( value );
} /* int Tree... */
git://source.jalview.org / jabaws.git / blobdiff