3 tree-dependent iteration
4 algorithm A+ when group-to-group, C when group-to-singleSeqence
20 #ifdef enablemultithread
21 typedef struct _threadarg
36 int *generationofinput;
40 int *generationofmastercopypt;
43 LocalHom **localhomtable;
48 float **tscorehistory_detail;
50 pthread_mutex_t *mutex;
51 pthread_cond_t *collection_end;
52 pthread_cond_t *collection_start;
57 static void shuffle( int *arr, int n )
76 static void Writeoption2( FILE *fp, int cycle, double cut )
78 fprintf( fp, "%dth cycle\n", cycle );
79 fprintf( fp, "marginal score to search : current score * (100-%d) / 100\n", (int)cut );
82 static void Writeoptions( FILE *fp )
84 fprintf( fp, "Tree-dependent-iteration\n" );
86 fprintf( fp, "Blosum %d\n", nblosum );
87 else if( scoremtx == -1 )
88 fprintf( fp, "DNA\n" );
89 else if( scoremtx == 2 )
90 fprintf( fp, "Miyata-Yasunaga\n" );
92 fprintf( fp, "JTT %dPAM\n", pamN );
94 if( scoremtx == 0 || scoremtx == 1 )
95 fprintf( fp, "Gap Penalty = %+5.3f, %5.2f, %+5.3f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
97 fprintf( fp, "Gap Penalty = %+5.3f\n", (double)penalty/1000 );
101 fprintf( fp, "score of rnd or sco\n" );
102 else if( scmtd == 4 )
103 fprintf( fp, "score = sigma( score for a pair of homologous amino acids ) / ( number of amino acids pairs )\n" );
104 else if( scmtd == 5 )
105 fprintf( fp, "score : SP\n" );
107 fprintf( fp, "?\n" );
111 fprintf( fp, "weighted rationale-1, geta2 = %f\n", geta2 );
112 else if( weight == 3 )
113 fprintf( fp, "weighted like ClustalW," );
114 else if( weight == 4 )
115 fprintf( fp, "weighted rationale-2, geta2 = %f\n", geta2 );
117 fprintf( fp, "unweighted\n" );
119 if( weight && utree )
120 fprintf( fp, "using tree defined by the file hat2.\n" );
121 if( weight && !utree )
122 fprintf( fp, "using temporary tree.\n" );
124 if( treemethod == 'n' )
125 fprintf( fp, "Tree is calculated with Neighbor-Joining Method.\n" );
126 else if( treemethod == 'q' )
127 fprintf( fp, "Tree is calculated with Minimum linkage.\n" );
128 else if( treemethod == 'X' )
129 fprintf( fp, "Tree is calculated with simplified UPG Method and UPG Method.\n" );
130 else if( treemethod == 'E' )
131 fprintf( fp, "Tree is calculated with UPG Method.\n" );
133 fprintf( fp, "Tree is calculated with unknown Method.\n" );
136 fprintf( fp, "Algorithm A+ / C\n" );
137 else if( alg == 'A' )
138 fprintf( fp, "Algorithm A+ \n" );
139 else if( alg == 'a' )
140 fprintf( fp, "Algorithm A \n" );
142 fprintf( fp, "Algorithm ? \n" );
148 fprintf( fp, "Basis : 4 nucleotides\n" );
153 fprintf( fp, "Basis : Polarity and Volume\n" );
155 fprintf( fp, "Basis : 20 amino acids\n" );
157 fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
158 fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
162 #ifdef enablemultithread
164 static void freelocalarrays(
165 float *tscorehistory,
166 RNApair ***grouprna1, RNApair ***grouprna2,
167 RNApair *rnapairboth,
168 char *indication1, char *indication2,
169 double *effarr, double *effarrforlocalhom, double *effarr1, double *effarr2,
170 char **mseq1, char **mseq2,
172 int *gapmap1, int *gapmap2,
173 double *effarr1_kozo, double *effarr2_kozo, double *effarr_kozo,
175 LocalHom *** localhomshrink
178 // fprintf( stderr, "Skipping freelocalarrays\n" );
181 if( commonIP ) FreeIntMtx( commonIP );
183 Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
184 Falign_localhom( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL,NULL, 0, NULL );
185 if( rnakozo && rnaprediction == 'm' )
187 free( grouprna1 ); // nakamimo?
188 free( grouprna2 ); // nakamimo?
191 free( tscorehistory );
195 free( effarrforlocalhom );
200 FreeCharMtx( localcopy );
204 free( effarr1_kozo );
205 free( effarr2_kozo );
210 if( rnakozo ) free( rnapairboth );
214 for( i=0; i<njob; i++)
216 free( localhomshrink[i] ); // nakamimo??
218 free( localhomshrink );
223 static void *athread( void *arg )
226 threadarg_t *targ = (threadarg_t *)arg;
227 int thread_no = targ->thread_no;
228 int njob = targ->njob;
229 int nbranch = targ->nbranch;
230 int maxiter = targ->maxiter;
231 int *ndonept = targ->ndonept;
232 int *ntrypt = targ->ntrypt;
233 int *collectingpt = targ->collectingpt;
234 int *jobposintpt = targ->jobposintpt;
235 int nkozo = targ->nkozo;
236 float *gainlist = targ->gainlist;
237 float *tscorelist = targ->tscorelist;
238 int *generationofinput = targ->generationofinput;
239 int *subgenerationpt = targ->subgenerationpt;
240 float *basegainpt = targ->basegainpt;
241 char *kozoarivec = targ->kozoarivec;
242 char **mastercopy = targ->mastercopy;
243 char ***candidates = targ->candidates;
244 int *generationofmastercopypt = targ->generationofmastercopypt;
245 int *branchtable = targ->branchtable;
246 RNApair ***singlerna = targ->singlerna;
247 LocalHom **localhomtable = targ->localhomtable;
248 int alloclen = targ->alloclen;
249 Node * stopol = targ->stopol;
250 int ***topol = targ->topol;
251 // double **len = targ->len;
252 float **tscorehistory_detail = targ->tscorehistory_detail;
253 int *finishpt = targ->finishpt;
263 char **mseq1, **mseq2;
264 double *effarr, *effarr_kozo; // re-calc
265 double *effarr1, *effarr2, *effarr1_kozo, *effarr2_kozo;
266 char *indication1, *indication2;
268 RNApair ***grouprna1, ***grouprna2;
269 RNApair *rnapairboth;
270 LocalHom ***localhomshrink;
271 int *gapmap1, *gapmap2;
272 float tscore, mscore, oimpmatch, impmatch;
275 float naivescore0 = 0, naivescore1;
276 double *effarrforlocalhom;
277 float *tscorehistory;
286 int subgenerationatfirst;
287 double unweightedspscore;
297 fprintf( stderr, "Dynamic tree is not supported in the multithread version.\n" );
300 if( score_check == 2 )
302 fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
308 fprintf( stderr, "Weight 2 is not supported in the multithread version.\n" );
311 if( cooling && cut > 0.0 )
313 fprintf( stderr, "Cooling is not supported in the multithread version.\n" );
317 tscorehistory = calloc( maxiter, sizeof( float ) );
319 if( rnakozo && rnaprediction == 'm' )
321 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
322 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
326 grouprna1 = grouprna2 = NULL;
329 indication1 = AllocateCharVec( njob*3+50 );
330 indication2 = AllocateCharVec( njob*3+50 );
331 effarr = AllocateDoubleVec( locnjob );
332 effarrforlocalhom = AllocateDoubleVec( locnjob );
333 effarr1 = AllocateDoubleVec( locnjob );
334 effarr2 = AllocateDoubleVec( locnjob );
335 mseq1 = AllocateCharMtx( locnjob, 0 );
336 mseq2 = AllocateCharMtx( locnjob, 0 );
337 localcopy = AllocateCharMtx( locnjob, alloclen );
338 gapmap1 = AllocateIntVec( alloclen );
339 gapmap2 = AllocateIntVec( alloclen );
341 effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
342 effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
343 effarr_kozo = AllocateDoubleVec( locnjob );
344 for( i=0; i<locnjob; i++ )
345 effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
347 pair = AllocateCharMtx( locnjob, locnjob );
350 if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
354 localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
355 for( i=0; i<njob; i++)
357 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom * ) );
367 for( iterate=0; iterate<maxiter; iterate++ )
369 pthread_mutex_lock( targ->mutex );
371 if( *collectingpt == 1 )
374 *generationofmastercopypt = iterate;
375 *subgenerationpt = 0;
379 for( i=0; i<nwa; i++ ) gainlist[i] = 0;
380 for( i=0; i<nwa; i++ ) tscorelist[i] = 0.0;
381 for( i=0; i<nbranch; i++ ) generationofinput[i] = -1;
382 if( parallelizationstrategy != BESTFIRST && randomseed != 0 ) shuffle( branchtable, nbranch );
383 pthread_cond_broadcast( targ->collection_end );
384 pthread_mutex_unlock( targ->mutex );
388 pthread_cond_broadcast( targ->collection_end );
389 pthread_mutex_unlock( targ->mutex );
393 grouprna1, grouprna2,
395 indication1, indication2,
396 effarr, effarrforlocalhom, effarr1, effarr2,
400 effarr1_kozo, effarr2_kozo, effarr_kozo,
405 pthread_exit( NULL );
408 pthread_mutex_lock( targ->mutex );
409 while( *ndonept < nbranch )
410 pthread_cond_wait( targ->collection_start, targ->mutex );
411 pthread_mutex_unlock( targ->mutex );
412 // fprintf( stderr, "Thread 0 got a signal, *collectionpt = %d\n", *collectingpt );
415 Hoka no thread ga keisan
418 pthread_mutex_lock( targ->mutex );
419 *collectingpt = 1; // chofuku
422 for( i=0; i<nbranch; i++ )
424 if( generationofinput[i] != iterate )
426 fprintf( stderr, "Error! generationofinput[%d] = %d, but iterate=%d\n", i, generationofinput[i], iterate );
433 maxgain = gainlist[1];
435 for( i=2; i<nwa; i++ )
437 if( gainlist[i] > maxgain )
439 maxgain = gainlist[i];
446 // fprintf( stderr, "\nGain = %f\n", maxgain );
447 // fprintf( stderr, "best gain = %f by thread %d\n", gainlist[bestthread], bestthread );
448 // fprintf( stderr, "tscorelist[best] = %f by thread %d\n", tscorelist[bestthread], bestthread );
449 if( parallelizationstrategy == BESTFIRST )
451 for( i=0; i<locnjob; i++ ) strcpy( mastercopy[i], candidates[bestthread][i] );
454 unweightedspscore = plainscore( locnjob, mastercopy );
455 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * nbranch, unweightedspscore );
456 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( locnjob * strlen( mastercopy[0] ) ) );
457 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
458 fprintf( stderr, "\n" );
462 // fprintf( stderr, "gain(%d, by %d) = %f\n", iterate, bestthread, maxgain );
463 for( i=iterate-1; i>0; i-- )
465 // if( iterate-i < 15 ) fprintf( stderr, "hist[%d] = %f\n", i, tscorehistory[i] );
466 if( tscorehistory[i] == tscorelist[bestthread] )
468 fprintf( stderr, "\nOscillating? %f == %f\n", tscorehistory[i], tscorelist[bestthread] );
473 tscorehistory[iterate] = tscorelist[bestthread];
478 fprintf( stderr, "\nConverged.\n" );
480 // pthread_cond_broadcast( targ->collection_end );
481 // pthread_mutex_unlock( targ->mutex );
482 // freelocalarrays();
484 // pthread_exit( NULL );
490 fprintf( stderr, "\nConverged2.\n" );
495 pthread_mutex_unlock( targ->mutex );
497 pthread_mutex_lock( targ->mutex );
498 fprintf( stderr, "\nReached %d\n", maxiter );
500 pthread_cond_broadcast( targ->collection_end );
501 pthread_mutex_unlock( targ->mutex );
505 grouprna1, grouprna2,
507 indication1, indication2,
508 effarr, effarrforlocalhom, effarr1, effarr2,
512 effarr1_kozo, effarr2_kozo, effarr_kozo,
517 pthread_exit( NULL );
524 if( iterate % 2 == 0 )
530 lin = locnjob - 2; ldf = -1;
532 for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
536 pthread_mutex_lock( targ->mutex );
537 while( *collectingpt > 0 )
538 pthread_cond_wait( targ->collection_end, targ->mutex );
539 if( *collectingpt == -1 )
541 pthread_mutex_unlock( targ->mutex );
545 grouprna1, grouprna2,
547 indication1, indication2,
548 effarr, effarrforlocalhom, effarr1, effarr2,
552 effarr1_kozo, effarr2_kozo, effarr_kozo,
557 pthread_exit( NULL );
559 // pthread_mutex_unlock( targ->mutex );
562 // pthread_mutex_lock( targ->mutex );
563 if( *jobposintpt == nbranch )
565 if( *collectingpt != -1 ) *collectingpt = 1; // chofuku
566 pthread_mutex_unlock( targ->mutex );
569 // fprintf( stderr, "JOB jobposintpt=%d\n", *jobposintpt );
570 myjob = branchtable[*jobposintpt];
572 if( l == locnjob-2 ) k = 1;
573 else k = myjob - l * 2;
574 // fprintf( stderr, "JOB l=%d, k=%d\n", l, k );
577 iterate = *generationofmastercopypt;
579 pthread_mutex_unlock( targ->mutex );
583 // fprintf( stderr, "branchpos = %d (thread %d)\n", branchpos, thread_no );
585 // fprintf( stderr, "iterate=%d, l=%d, k=%d (thread %d)\n", iterate, l, k, thread_no );
588 fprintf( stderr, "STEP %03d-%03d-%d (Thread %d) ", iterate+1, l+1, k, thread_no );
589 fprintf( stderr, "STEP %03d-%03d-%d (thread %d) %s ", iterate+1, l+1, k, thread_no, use_fft?"\n":"\n" );
591 for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
593 OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
595 fprintf( stderr, "STEP%d-%d\n", l, k );
596 for( i=0; i<locnjob; i++ )
598 for( j=0; j<locnjob; j++ )
600 fprintf( stderr, "%#3d", pair[i][j] );
602 fprintf( stderr, "\n" );
607 for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
610 for( i=0; i<locnjob; i++ )
613 effarr_kozo[i] = 1.0;
615 effarr_kozo[i] = 0.0;
619 else if( weight == 4 )
621 weightFromABranch( locnjob, effarr, stopol, topol, l, k );
622 if( nkozo ) // hitomadu single weight.
624 for( i=0; i<locnjob; i++ )
626 if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
627 else effarr_kozo[i] = 0.0;
633 fprintf( stderr, "weight error!\n" );
639 pthread_mutex_lock( targ->mutex );
640 for( i=0; i<locnjob; i++ ) strcpy( localcopy[i], mastercopy[i] );
641 subgenerationatfirst = *subgenerationpt;
642 pthread_mutex_unlock( targ->mutex );
643 length = strlen( localcopy[0] );
649 clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
650 for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
652 clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
653 for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
656 fprintf( stderr, "\ngroup1 = %s\n", indication1 );
657 for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[] = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
658 fprintf( stderr, "\ngroup2 = %s\n", indication2 );
659 for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[] = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
664 clus1 = conjuctionfortbfast( pair, s1, localcopy, mseq1, effarr1, effarr, indication1 );
665 clus2 = conjuctionfortbfast( pair, s2, localcopy, mseq2, effarr2, effarr, indication2 );
668 if( rnakozo && rnaprediction == 'm' )
670 makegrouprnait( grouprna1, singlerna, pair, s1 );
671 makegrouprnait( grouprna2, singlerna, pair, s2 );
674 if( score_check == 2 )
676 fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
679 else if( score_check )
681 if( RNAscoremtx == 'r' )
682 intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
684 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
688 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
689 // weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink ); // >>>
695 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
696 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
697 for( i=length-1; i>=0; i-- )
699 oimpmatch += part_imp_match_out_scQ( i, i );
700 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
705 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
706 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
707 for( i=length-1; i>=0; i-- )
709 oimpmatch += part_imp_match_out_sc( i, i );
710 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
713 // fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
719 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
720 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
722 for( i=length-1; i>=0; i-- )
724 oimpmatch += imp_match_out_scQ( i, i );
725 // fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
730 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
732 fprintf( stderr, "not supported\n" );
735 for( i=length-1; i>=0; i-- )
737 oimpmatch += imp_match_out_sc( i, i );
738 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
741 // fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
743 // fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
751 // fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
752 mscore = (double)oimpmatch + tmpdouble;
756 fprintf( stderr, "score_check = %d\n", score_check );
757 fprintf( stderr, "Not supported\n" );
762 // if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
764 // if( !use_fft && !rnakozo )
767 commongappick_record( clus1, mseq1, gapmap1 );
768 commongappick_record( clus2, mseq2, gapmap2 );
772 fprintf( stderr, "##### mscore = %f\n", mscore );
778 fprintf( trap_g, "\n%d-%d-%d\n", iterate+1, l+1, k );
779 fprintf( trap_g, "group1 = %s\n", indication1 );
780 fprintf( trap_g, "group2 = %s\n", indication2 );
786 printf( "STEP %d-%d-%d\n", iterate, l, k );
787 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
790 if( constraint == 2 )
795 // part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
797 // part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
799 Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, subgenerationpt, subgenerationatfirst, &chudanres );
800 // fprintf( stderr, "##### impmatch = %f\n", impmatch );
801 if( chudanres && parallelizationstrategy == BAATARI2 )
803 // fprintf( stderr, "#### yarinaoshi!!! INS-i\n" );
812 // imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap wo tsukuttakara iranai.
813 // if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
815 wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
816 fprintf( stderr, "wm = %f\n", wm );
818 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
819 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
820 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
822 if( naivescore1 > naivescore0 )
823 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
824 else if( naivescore1 < naivescore0 )
825 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
827 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
829 if( abs( wm - naivescore1 ) > 100 )
831 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
833 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
834 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
840 else if( alg == 'R' )
843 imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
844 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
845 // fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
846 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
847 // fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
849 if( naivescore1 > naivescore0 )
850 fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
851 else if( naivescore1 < naivescore0 )
852 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
854 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
856 if( abs( wm - naivescore1 ) > 100 )
858 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
860 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
861 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
866 else if( alg == 'H' )
869 imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
870 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
871 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
872 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
873 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
875 if( naivescore1 > naivescore0 )
876 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
877 else if( naivescore1 < naivescore0 )
878 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
880 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
882 if( abs( wm - naivescore1 ) > 100 )
884 // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
885 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
886 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
893 // imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
894 A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
895 fprintf( stderr, "A__align_gapmap\n" );
896 // fprintf( stderr, "##### impmatch = %f\n", impmatch );
904 totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, subgenerationpt, subgenerationatfirst, &chudanres );
905 if( chudanres && parallelizationstrategy == BAATARI2 )
907 // fprintf( stderr, "#### yarinaoshi!!! FFT-NS-i\n" );
911 // fprintf( stderr, "totalwm = %f\n", totalwm );
915 fprintf( stderr, "totalwm = %f\n", totalwm );
916 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
918 if( naivescore1 > naivescore0 )
919 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
920 else if( naivescore1 < naivescore0 )
921 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
923 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
925 if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
927 // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
928 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
929 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
937 fprintf( stderr, "totalwm = %f\n", totalwm );
938 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
940 if( naivescore1 > naivescore0 )
941 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
942 else if( naivescore1 < naivescore0 )
943 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
945 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
947 if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
949 // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
950 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
951 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
962 MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, subgenerationpt, subgenerationatfirst, &chudanres, outgap, outgap );
963 if( chudanres && parallelizationstrategy == BAATARI2 )
965 // fprintf( stderr, "#### yarinaoshi!!! NW-NS-i\n" );
969 else if( alg == 'A' )
971 A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); //outgap==1
973 else if( alg == 'Q' )
976 wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
977 fprintf( stderr, "wm = %f\n", wm );
978 fprintf( stderr, "impmatch = %f\n", impmatch );
979 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
981 if( naivescore1 > naivescore0 )
982 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
983 else if( naivescore1 < naivescore0 )
984 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
986 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
988 if( abs( wm - naivescore1 ) > 100 )
990 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
993 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
994 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
999 else if( alg == 'R' )
1002 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
1003 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1005 if( naivescore1 > naivescore0 )
1006 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1007 else if( naivescore1 < naivescore0 )
1008 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
1010 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1012 if( abs( wm - naivescore1 ) > 100 )
1014 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
1015 // rewind( stderr );
1016 // rewind( stdout );
1017 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1018 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1023 else if( alg == 'H' )
1026 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
1027 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1029 if( naivescore1 > naivescore0 )
1030 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1031 else if( naivescore1 < naivescore0 )
1033 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
1036 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1039 if( abs( wm - naivescore1 ) > 100 )
1042 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1043 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1048 else if( alg == 'a' )
1050 Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
1052 else ErrorExit( "Sorry!" );
1054 // fprintf( stderr, "## impmatch = %f\n", impmatch );
1057 if( parallelizationstrategy == BAATARI2 && *subgenerationpt != subgenerationatfirst )
1059 // fprintf( stderr, "\nYarinaoshi2!! (Thread %d)\n", thread_no );
1064 identity = !strcmp( localcopy[s1], mastercopy[s1] );
1065 identity *= !strcmp( localcopy[s2], mastercopy[s2] );
1068 /* Bug? : idnetitcal but score change when scoreing mtx != JTT */
1070 length = strlen( mseq1[0] );
1075 // if( !devide ) fprintf( trap_g, "tscore = %f identical.\n", tscore );
1076 // fprintf( stderr, " identical." );
1077 fprintf( stderr, "%03d-%04d-%d (thread %4d) identical \r", iterate+1, *ndonept, k, thread_no );
1083 if( constraint == 2 )
1086 if( RNAscoremtx == 'r' )
1087 intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1089 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1091 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1094 tscore = impmatch + tmpdouble;
1096 // fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
1100 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1103 // fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
1105 for( i=0; i<1; i++ )
1106 fprintf( stderr, "%s\n", mseq1[i] );
1107 fprintf( stderr, "+++++++\n" );
1108 for( i=0; i<1; i++ )
1109 fprintf( stderr, "%s\n", mseq2[i] );
1115 tscore = mscore + 1.0;
1117 // fprintf( stderr, "in line 705, tscore=%f\n", tscore );
1118 // for( i=0; i<length; i++ )
1119 // tscore = tscore + (double)mseq1[0][i];
1120 // mscore = tscore - 1.0;
1123 if( isnan( mscore ) )
1125 fprintf( stderr, "\n\nmscore became NaN\n" );
1128 if( isnan( tscore ) )
1130 fprintf( stderr, "\n\ntscore became NaN\n" );
1136 // fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
1139 if( parallelizationstrategy == BAATARI1 && *subgenerationpt != subgenerationatfirst )
1141 // fprintf( stderr, "\nYarinaoshi1!! (Thread %d)\n", thread_no );
1145 gain = tscore - ( mscore - cut/100.0*mscore );
1148 if( parallelizationstrategy == BESTFIRST )
1150 if( gain > gainlist[thread_no] )
1152 gainlist[thread_no] = gain;
1153 for( i=0; i<locnjob; i++ ) strcpy( candidates[thread_no][i], localcopy[i] );
1154 tscorelist[thread_no] = tscore;
1155 // if( iterate == 0 ) fprintf( stderr, "hist %d-%d-%d, gain=%f (Thread %d)\n", iterate, l, k, gain, thread_no );
1160 pthread_mutex_lock( targ->mutex );
1161 for( i=0; i<locnjob; i++ ) strcpy( mastercopy[i], localcopy[i] );
1162 *subgenerationpt += 1;
1163 gainlist[thread_no] = *basegainpt + gain;
1164 *basegainpt += gain;
1168 unweightedspscore = plainscore( locnjob, localcopy );
1169 fprintf( stderr, "\nSCORE %d = %.0f, ", *ntrypt, unweightedspscore );
1170 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( locnjob * strlen( mastercopy[0] ) ) );
1171 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
1172 fprintf( stderr, "\n" );
1175 pthread_mutex_unlock( targ->mutex );
1176 tscorelist[thread_no] = tscore;
1179 fprintf( stderr, "tscore = %f mscore = %f accepted.\n", tscore, mscore );
1180 fprintf( stderr, "\nbetter! gain = %f (thread %d)\r", gain, thread_no );
1182 fprintf( stderr, "%03d-%04d-%d (thread %4d) better \r", iterate+1, *ndonept, k, thread_no );
1189 fprintf( stderr, "tscore = %f mscore = %f rejected.\r", tscore, mscore );
1190 fprintf( stderr, "worse! gain = %f", gain );
1192 fprintf( stderr, "%03d-%04d-%d (thread %4d) worse \r", iterate+1, *ndonept, k, thread_no );
1198 for( ii=iterate-2; ii>=0; ii-=1 )
1200 // fprintf( stderr, "Checking tscorehistory %f ?= %f\n", tscore, tscorehistory_detail[ii][branchpos] );
1201 if( tscore == tscorehistory_detail[ii][branchpos] )
1207 if( parallelizationstrategy != BESTFIRST && converged2 )
1209 // fprintf( stderr, "\nFINISH!\n" );
1210 pthread_mutex_lock( targ->mutex );
1212 pthread_mutex_unlock( targ->mutex );
1215 tscorehistory_detail[iterate][branchpos] = tscore;
1216 fprintf( stderr, "\r" );
1218 pthread_mutex_lock( targ->mutex );
1220 // fprintf( stderr, "*ndonept = %d, nbranch = %d (thread %d) iterate=%d\n", *ndonept, nbranch, thread_no, iterate );
1221 generationofinput[branchpos] = iterate;
1222 if( *ndonept == nbranch )
1224 if( *collectingpt != -1 ) *collectingpt = 1; // chofuku
1225 // fprintf( stderr, "Thread %d sends a signal, *ndonept = %d\n", thread_no, *ndonept );
1226 pthread_cond_signal( targ->collection_start );
1228 pthread_mutex_unlock( targ->mutex );
1231 } /* for( iterate ) */
1237 int TreeDependentIteration( int locnjob, char **name, int nlen[M],
1238 char **aseq, char **bseq, int ***topol, double **len,
1239 int alloclen, LocalHom **localhomtable,
1240 RNApair ***singlerna,
1241 int nkozo, char *kozoarivec )
1243 int i, j, k, l, iterate, ii, iu, ju;
1244 int lin, ldf, length;
1247 static double **imanoten;
1248 static Node *stopol;
1249 static double *effarrforlocalhom = NULL;
1250 static double *effarr = NULL;
1251 static double *effarr1 = NULL;
1252 static double *effarr2 = NULL;
1253 static double *effarr_kozo = NULL;
1254 static double *effarr1_kozo = NULL;
1255 static double *effarr2_kozo = NULL;
1256 static double **mtx = NULL;
1257 static int **node = NULL;
1258 static int *branchnode = NULL;
1259 static double **branchWeight = NULL;
1260 static char **mseq1, **mseq2;
1261 static float ***history;
1263 double tscore, mscore;
1267 float naivescore0 = 0.0; // by D.Mathog, a guess
1270 char pair[njob][njob];
1275 double score_for_check0, score_for_check1;
1276 static double **effmtx = NULL;
1277 extern double score_calc0();
1279 static char *indication1, *indication2;
1280 static LocalHom ***localhomshrink = NULL;
1281 float impmatch = 0.0, oimpmatch = 0.0;
1282 static int *gapmap1;
1283 static int *gapmap2;
1286 static RNApair *rnapairboth;
1287 RNApair ***grouprna1, ***grouprna2;
1288 double unweightedspscore;
1290 if( rnakozo && rnaprediction == 'm' )
1292 grouprna1 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
1293 grouprna2 = (RNApair ***)calloc( njob, sizeof( RNApair ** ) );
1297 grouprna1 = grouprna2 = NULL;
1300 Writeoptions( trap_g );
1303 if( effarr == NULL ) /* locnjob == njob ni kagiru */
1305 indication1 = AllocateCharVec( njob*3+50 );
1306 indication2 = AllocateCharVec( njob*3+50 );
1307 effarr = AllocateDoubleVec( locnjob );
1308 effarrforlocalhom = AllocateDoubleVec( locnjob );
1309 effarr1 = AllocateDoubleVec( locnjob );
1310 effarr2 = AllocateDoubleVec( locnjob );
1311 mseq1 = AllocateCharMtx( locnjob, 0 );
1312 mseq2 = AllocateCharMtx( locnjob, 0 );
1313 mtx = AllocateDoubleMtx( locnjob, locnjob );
1314 node = AllocateIntMtx( locnjob, locnjob );
1315 branchnode = AllocateIntVec( locnjob );
1316 branchWeight = AllocateDoubleMtx( locnjob, 2 );
1317 history = AllocateFloatCub( niter, locnjob, 2 );
1318 stopol = (Node *)calloc( locnjob * 2, sizeof( Node ) );
1319 gapmap1 = AllocateIntVec( alloclen );
1320 gapmap2 = AllocateIntVec( alloclen );
1321 if( score_check == 2 ) imanoten = AllocateDoubleMtx( njob, njob );
1323 effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
1324 effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
1325 effarr_kozo = AllocateDoubleVec( locnjob );
1326 for( i=0; i<locnjob; i++ )
1327 effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
1331 pair = AllocateCharMtx( locnjob, locnjob );
1332 if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
1336 localhomshrink = (LocalHom ***)calloc( njob, sizeof( LocalHom ** ) );
1337 for( i=0; i<njob; i++)
1339 localhomshrink[i] = (LocalHom **)calloc( njob, sizeof( LocalHom * ) );
1345 if( !effmtx ) effmtx = AllocateDoubleMtx( locnjob, locnjob );
1346 for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) effmtx[i][j] = 1.0;
1349 for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
1351 writePre( locnjob, name, nlen, aseq, 0 );
1357 counteff_simple( locnjob, topol, len, effarrforlocalhom );
1358 calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
1363 countnode_int( locnjob, topol, node );
1366 fprintf( stderr, "Not supported, weight=%d nkozo=%d.\n", weight, nkozo );
1369 else if( weight == 4 )
1371 treeCnv( stopol, locnjob, topol, len, branchWeight );
1372 calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
1376 #ifdef enablemultithread
1381 pthread_mutex_t mutex;
1382 pthread_cond_t collection_end;
1383 pthread_cond_t collection_start;
1385 int generationofmastercopy;
1388 int *generationofinput;
1398 float **tscorehistory_detail;
1402 nbranch = (njob-1) * 2 - 1;
1405 targ = calloc( nwa, sizeof( threadarg_t ) );
1406 handle = calloc( nwa, sizeof( pthread_t ) );
1407 pthread_mutex_init( &mutex, NULL );
1408 pthread_cond_init( &collection_end, NULL );
1409 pthread_cond_init( &collection_start, NULL );
1411 gainlist = calloc( nwa, sizeof( float ) );
1412 tscorelist = calloc( nwa, sizeof( float ) );
1413 branchtable = calloc( nbranch, sizeof( int ) );
1414 generationofinput = calloc( nbranch, sizeof( int ) );
1415 if( parallelizationstrategy == BESTFIRST )
1416 candidates = AllocateCharCub( nwa, locnjob, alloclen );
1417 for( i=0; i<nbranch; i++ ) branchtable[i] = i;
1418 tscorehistory_detail = AllocateFloatMtx( maxiter, nbranch );
1422 for( i=0; i<nwa; i++ )
1424 targ[i].thread_no = i;
1425 targ[i].njob = njob;
1426 targ[i].nbranch = nbranch;
1427 targ[i].maxiter = maxiter;
1428 targ[i].ndonept = &ndone;
1429 targ[i].ntrypt = &ntry;
1430 targ[i].collectingpt = &collecting;
1431 targ[i].jobposintpt = &jobposint;
1432 targ[i].gainlist = gainlist;
1433 targ[i].tscorelist = tscorelist;
1434 targ[i].nkozo = nkozo;
1435 targ[i].kozoarivec = kozoarivec;
1436 targ[i].mastercopy = bseq;
1437 targ[i].candidates = candidates;
1438 targ[i].subgenerationpt = &subgeneration;
1439 targ[i].basegainpt = &basegain;
1440 targ[i].generationofmastercopypt = &generationofmastercopy;
1441 targ[i].generationofinput = generationofinput;
1442 targ[i].branchtable = branchtable;
1443 targ[i].singlerna = singlerna;
1444 targ[i].localhomtable = localhomtable;
1445 targ[i].alloclen = alloclen;
1446 targ[i].stopol = stopol;
1447 targ[i].topol = topol;
1448 // targ[i].len = len;
1449 targ[i].mutex = &mutex;
1450 targ[i].collection_end = &collection_end;
1451 targ[i].collection_start = &collection_start;
1452 targ[i].tscorehistory_detail = tscorehistory_detail;
1453 targ[i].finishpt = &finish;
1455 pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
1458 for( i=0; i<nwa; i++ )
1460 pthread_join( handle[i], NULL );
1463 pthread_mutex_destroy( &mutex );
1464 pthread_cond_destroy( &collection_end );
1465 pthread_cond_destroy( &collection_start );
1471 free( branchtable );
1472 free( generationofinput );
1473 if( parallelizationstrategy == BESTFIRST )
1474 FreeCharCub( candidates );
1475 FreeFloatMtx( tscorehistory_detail );
1486 nbranch = (njob-1) * 2 - 1;
1488 branchtable = calloc( nbranch, sizeof( int ) );
1489 for( i=0; i<nbranch; i++ ) branchtable[i] = i;
1491 srand( randomseed );
1494 if( parallelizationstrategy == BESTFIRST )
1496 fprintf( stderr, "Not implemented. Try --thread 1 --bestfirst\n" );
1500 if( cooling ) cut *= 2.0;
1501 for( iterate = 0; iterate<niter; iterate++ )
1503 if( cooling ) cut *= 0.5; /* ... */
1506 if( randomseed != 0 ) shuffle( branchtable, nbranch );
1509 fprintf( trap_g, "\n" );
1510 Writeoption2( trap_g, iterate, cut );
1511 fprintf( trap_g, "\n" );
1518 fprintf( stderr, "The combination of dynamic tree and kozo is not supported.\n" );
1523 static char *buff1 = NULL;
1524 static char *buff2 = NULL;
1527 buff1 = AllocateCharVec( alloclen );
1528 buff2 = AllocateCharVec( alloclen );
1531 for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
1533 buff1[0] = buff2[0] = 0;
1534 strcat( buff1, res_g[i] ); strcat( buff2, res_g[j] );
1535 strcat( buff1, bseq[i] ); strcat( buff2, bseq[j] );
1536 strcat( buff1, seq_g[i] ); strcat( buff2, seq_g[j] );
1538 mtx[i][j] = (double)substitution_hosei( buff1, buff2 );
1543 for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
1544 mtx[i][j] = (double)substitution_hosei( bseq[i], bseq[j] );
1547 if ( treemethod == 'n' )
1548 nj( locnjob, mtx, topol, len );
1549 else if( treemethod == 's' )
1550 spg( locnjob, mtx, topol, len );
1551 else if( treemethod == 'X' )
1552 supg( locnjob, mtx, topol, len );
1553 else if( treemethod == 'p' )
1554 upg2( locnjob, mtx, topol, len );
1555 /* veryfastsupg
\e$B$O!":#$N$H$3$m;H$($^$;$s!#
\e(B*/
1556 /*
\e$B=gHV$NLdBj$,$"$k$N$G
\e(B */
1559 countnode_int( locnjob, topol, node );
1560 else if( weight == 4 )
1562 treeCnv( stopol, locnjob, topol, len, branchWeight );
1563 calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
1565 trap = fopen( "hat2", "w" );
1566 if( !trap ) ErrorExit( "Cannot open hat2." );
1567 WriteHat2_pointer( trap, locnjob, name, mtx );
1571 counteff_simple( locnjob, topol, len, effarrforlocalhom );
1572 calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
1576 if( iterate % 2 == 0 )
1582 lin = locnjob - 2; ldf = -1;
1585 if( score_check == 2 )
1589 length = strlen( bseq[0] );
1590 for( i=0; i<locnjob-1; i++ )
1591 for( j=i+1; j<locnjob; j++ )
1592 intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
1596 for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
1600 for( k=0; k<2; k++ )
1602 if( l == locnjob-2 ) k = 1;
1605 for( jobpos=0; jobpos<nbranch; jobpos++)
1608 myjob = branchtable[jobpos];
1610 if( l == locnjob-2 ) k = 1;
1611 else k = myjob - l * 2;
1614 fprintf( stderr, "STEP %03d-%03d-%d ", iterate+1, l+1, k );
1617 fprintf( stderr, "STEP %03d-%03d-%d %s", iterate+1, l+1, k, use_fft?"\n":"\n" );
1619 for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
1621 OneClusterAndTheOther( locnjob, pair, &s1, &s2, topol, l, k );
1623 fprintf( stderr, "STEP%d-%d\n", l, k );
1624 for( i=0; i<locnjob; i++ )
1626 for( j=0; j<locnjob; j++ )
1628 fprintf( stderr, "%#3d", pair[i][j] );
1630 fprintf( stderr, "\n" );
1635 for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
1638 for( i=0; i<locnjob; i++ )
1641 effarr_kozo[i] = 1.0;
1643 effarr_kozo[i] = 0.0;
1647 else if( weight == 2 )
1649 nodeFromABranch( locnjob, branchnode, node, topol, len, l, k );
1650 node_eff( locnjob, effarr, branchnode );
1652 else if( weight == 4 )
1654 weightFromABranch( locnjob, effarr, stopol, topol, l, k );
1658 assignstrweight( locnjob, effarr_kozo, stopol, topol, l, k, kozoarivec, effarr );
1662 if( nkozo ) // hitomadu single weight.
1663 for( i=0; i<locnjob; i++ )
1665 if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
1666 else effarr_kozo[i] = 0.0;
1670 fprintf( stderr, "\n" );
1671 fprintf( stderr, "effarr_kozo = \n" );
1672 for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr_kozo[i] );
1673 fprintf( stderr, "\n" );
1674 fprintf( stderr, "effarr = \n" );
1675 for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr[i] );
1676 fprintf( stderr, "\n\n" );
1680 for( i=0; i<locnjob; i++ ) strcpy( aseq[i], bseq[i] );
1681 length = strlen( aseq[0] );
1688 clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1689 for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
1691 clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1692 for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
1695 fprintf( stderr, "\ngroup1 = %s\n", indication1 );
1696 for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[] = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
1697 fprintf( stderr, "\ngroup2 = %s\n", indication2 );
1698 for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[] = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
1706 clus1 = conjuctionfortbfast_kozo_BUG( pair, s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1707 clus2 = conjuctionfortbfast_kozo_BUG( pair, s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1712 clus1 = conjuctionfortbfast( pair, s1, aseq, mseq1, effarr1, effarr, indication1 );
1713 clus2 = conjuctionfortbfast( pair, s2, aseq, mseq2, effarr2, effarr, indication2 );
1718 if( rnakozo && rnaprediction == 'm' )
1720 makegrouprnait( grouprna1, singlerna, pair, s1 );
1721 makegrouprnait( grouprna2, singlerna, pair, s2 );
1724 if( score_check == 2 )
1728 // msshrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1729 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1735 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1736 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1737 for( i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_scQ( i, i );
1741 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1742 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1743 for( i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_sc( i, i );
1750 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1751 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1752 for( i=length-1; i>=0; i-- ) oimpmatch += imp_match_out_scQ( i, i );
1756 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1757 fprintf( stderr, "not supported\n" );
1761 // fprintf( stderr, "### oimpmatch = %f\n", oimpmatch );
1769 for( i=s1; i<locnjob; i++ )
1771 if( !pair[s1][i] ) continue;
1773 for( j=s2; j<locnjob; j++ )
1775 if( !pair[s2][j] ) continue;
1776 // fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
1777 tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(i,j)][MAX(i,j)];
1782 mscore = oimpmatch + tmpdouble;
1784 else if( score_check )
1787 if( RNAscoremtx == 'r' )
1788 intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1790 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1792 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1797 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1798 // weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink ); // >>>
1804 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1805 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1806 for( i=length-1; i>=0; i-- )
1808 oimpmatch += part_imp_match_out_scQ( i, i );
1809 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1814 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1815 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1816 for( i=length-1; i>=0; i-- )
1818 oimpmatch += part_imp_match_out_sc( i, i );
1819 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1822 // fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
1828 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1829 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1831 for( i=length-1; i>=0; i-- )
1833 oimpmatch += imp_match_out_scQ( i, i );
1834 // fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
1839 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1841 fprintf( stderr, "not supported\n" );
1844 for( i=length-1; i>=0; i-- )
1846 oimpmatch += imp_match_out_sc( i, i );
1847 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1850 // fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
1852 // fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
1860 // fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
1861 mscore = (double)oimpmatch + tmpdouble;
1865 // fprintf( stderr, "score_check = %d\n", score_check );
1866 /* atode kousokuka */
1867 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1869 /* atode kousokuka */
1874 shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1879 part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1880 if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1884 part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1885 if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1892 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1893 if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1897 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1898 fprintf( stderr, "Not supported\n" );
1911 imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1912 for( i=length-1; i>=0; i-- )
1914 oimpmatch += imp_match_out_scQ( i, i );
1915 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1920 imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1921 for( i=length-1; i>=0; i-- )
1923 oimpmatch += imp_match_out_sc( i, i );
1924 // fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1931 naivescore0 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1932 else if( alg == 'Q' )
1933 naivescore0 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1934 else if( alg == 'R' )
1935 naivescore0 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1938 // if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
1940 // if( !use_fft && !rnakozo )
1943 commongappick_record( clus1, mseq1, gapmap1 );
1944 commongappick_record( clus2, mseq2, gapmap2 );
1948 fprintf( stderr, "##### mscore = %f\n", mscore );
1954 fprintf( trap_g, "\nSTEP%d-%d-%d\n", iterate+1, l+1, k );
1955 fprintf( trap_g, "group1 = %s\n", indication1 );
1956 fprintf( trap_g, "group2 = %s\n", indication2 );
1962 printf( "STEP %d-%d-%d\n", iterate, l, k );
1963 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
1966 if( constraint == 2 )
1971 // part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1973 // part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1974 Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, NULL, 0, NULL );
1975 // fprintf( stderr, "##### impmatch = %f\n", impmatch );
1982 // imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap wo tsukuttakara iranai.
1983 // if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
1985 wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
1986 fprintf( stderr, "wm = %f\n", wm );
1988 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
1989 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
1990 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
1992 if( naivescore1 > naivescore0 )
1993 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1994 else if( naivescore1 < naivescore0 )
1995 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
1997 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1999 if( abs( wm - naivescore1 ) > 100 )
2001 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2002 // rewind( stdout );
2003 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2004 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2010 else if( alg == 'R' )
2013 imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
2014 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
2015 // fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2016 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2017 // fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2019 if( naivescore1 > naivescore0 )
2020 fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
2021 else if( naivescore1 < naivescore0 )
2022 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2024 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2026 if( abs( wm - naivescore1 ) > 100 )
2028 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2030 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2031 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2036 else if( alg == 'H' )
2039 imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
2040 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
2041 fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2042 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2043 fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2045 if( naivescore1 > naivescore0 )
2046 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2047 else if( naivescore1 < naivescore0 )
2048 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2050 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2052 if( abs( wm - naivescore1 ) > 100 )
2054 // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2055 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2056 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2063 // imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
2064 A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
2065 // fprintf( stderr, "##### impmatch = %f\n", impmatch );
2072 totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
2074 // fprintf( stderr, "totalwm = %f\n", totalwm );
2078 fprintf( stderr, "totalwm = %f\n", totalwm );
2079 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2081 if( naivescore1 > naivescore0 )
2082 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2083 else if( naivescore1 < naivescore0 )
2084 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2086 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2088 if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
2090 // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2091 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2092 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2100 fprintf( stderr, "totalwm = %f\n", totalwm );
2101 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2103 if( naivescore1 > naivescore0 )
2104 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2105 else if( naivescore1 < naivescore0 )
2106 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2108 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2110 if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
2112 // fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2113 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2114 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2124 MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
2126 else if( alg == 'A' )
2128 A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); // outgap==1
2130 else if( alg == 'Q' )
2133 wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2134 fprintf( stderr, "wm = %f\n", wm );
2135 fprintf( stderr, "impmatch = %f\n", impmatch );
2136 naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2138 if( naivescore1 > naivescore0 )
2139 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2140 else if( naivescore1 < naivescore0 )
2141 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
2143 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2145 if( abs( wm - naivescore1 ) > 100 )
2147 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2148 // rewind( stderr );
2149 // rewind( stdout );
2150 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2151 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2156 else if( alg == 'R' )
2159 wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2160 naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2162 if( naivescore1 > naivescore0 )
2163 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2164 else if( naivescore1 < naivescore0 )
2165 fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2167 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2169 if( abs( wm - naivescore1 ) > 100 )
2171 // fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2172 // rewind( stderr );
2173 // rewind( stdout );
2174 // for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2175 // for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2180 else if( alg == 'H' )
2183 wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2184 naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2186 if( naivescore1 > naivescore0 )
2187 fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2188 else if( naivescore1 < naivescore0 )
2190 fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
2193 fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2196 if( abs( wm - naivescore1 ) > 100 )
2199 for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2200 for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2205 else if( alg == 'a' )
2207 Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
2209 else ErrorExit( "Sorry!" );
2211 // fprintf( stderr, "## impmatch = %f\n", impmatch );
2215 extern double DSPscore();
2216 extern double SSPscore();
2221 pre = SSPscore( locnjob, bseq );
2222 cur = SSPscore( locnjob, aseq );
2224 pre = DSPscore( locnjob, bseq );
2225 cur = DSPscore( locnjob, aseq );
2227 fprintf( stderr, "Previous Sscore = %f\n", pre );
2228 fprintf( stderr, "Currnet Sscore = %f\n\n", cur );
2231 // fprintf( stderr, "## impmatch = %f\n", impmatch );
2232 identity = !strcmp( aseq[s1], bseq[s1] );
2233 identity *= !strcmp( aseq[s2], bseq[s2] );
2236 /* Bug? : idnetitcal but score change when scoreing mtx != JTT */
2238 length = strlen( mseq1[0] );
2243 if( !devide ) fprintf( trap_g, "tscore = %f identical.\n", tscore );
2244 fprintf( stderr, " identical." );
2251 if( constraint == 2 )
2254 if( RNAscoremtx == 'r' )
2255 intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2257 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2259 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2262 tscore = impmatch + tmpdouble;
2264 // fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
2268 intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2271 // fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
2273 for( i=0; i<1; i++ )
2274 fprintf( stderr, "%s\n", mseq1[i] );
2275 fprintf( stderr, "+++++++\n" );
2276 for( i=0; i<1; i++ )
2277 fprintf( stderr, "%s\n", mseq2[i] );
2283 tscore = mscore + 1.0;
2285 // fprintf( stderr, "in line 705, tscore=%f\n", tscore );
2286 // for( i=0; i<length; i++ )
2287 // tscore = tscore + (double)mseq1[0][i];
2288 // mscore = tscore - 1.0;
2291 if( isnan( mscore ) )
2293 fprintf( stderr, "\n\nmscore became NaN\n" );
2296 if( isnan( tscore ) )
2298 fprintf( stderr, "\n\ntscore became NaN\n" );
2304 // fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
2306 if( tscore > mscore - cut/100.0*mscore )
2308 writePre( locnjob, name, nlen, aseq, 0 );
2309 for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
2310 if( score_check == 2 )
2314 for( i=0; i<locnjob-1; i++ )
2315 for( j=i+1; j<locnjob; j++ )
2316 intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
2320 fprintf( stderr, "tscore = %f mscore = %f accepted.\n", tscore, mscore );
2322 fprintf( stderr, " accepted." );
2329 fprintf( stderr, "tscore = %f mscore = %f rejected.\n", tscore, mscore );
2331 fprintf( stderr, " rejected." );
2336 fprintf( stderr, "\r" );
2339 history[iterate][l][k] = (float)tscore;
2341 // fprintf( stderr, "tscore = %f\n", tscore );
2343 if( converged >= locnjob * 2 )
2345 fprintf( trap_g, "Converged.\n\n" );
2346 fprintf( stderr, "\nConverged.\n\n" );
2349 unweightedspscore = plainscore( njob, bseq );
2350 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2351 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2352 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
2353 fprintf( stderr, "\n\n" );
2359 /* oscillation check */
2361 for( ii=iterate-2; ii>=0; ii-=2 )
2363 if( (float)tscore == history[ii][l][k] )
2369 if( ( oscillating && !cooling ) || ( oscillating && cut < 0.001 && cooling ) )
2371 fprintf( trap_g, "Oscillating.\n" );
2372 fprintf( stderr, "\nOscillating.\n\n" );
2375 unweightedspscore = plainscore( njob, bseq );
2376 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2377 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2378 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
2379 fprintf( stderr, "\n\n" );
2381 #if 1 /* hujuubun */
2385 } /* if( iterate ) */
2390 unweightedspscore = plainscore( njob, bseq );
2391 fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2392 fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2393 if( weight || constraint ) fprintf( stderr, " (differs from the objective score)" );
2394 fprintf( stderr, "\n\n" );
2396 } /* for( iterate ) */