source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/io.c

Last change on this file was 19480, checked in by westram, 15 months ago
File size: 116.2 KB
Line 
1#include "mltaln.h"
2
3static int upperCase = 0;
4
5#define DEBUG   0
6#define IODEBUG 0
7
8char creverse( char f )
9{
10        static char *table = NULL;
11        if( table == NULL )
12        {
13                int i;
14                table = AllocateCharVec(0x80);
15                for( i=0; i<0x80; i++ ) table[i] = i;
16                table['A'] = 'T';
17                table['C'] = 'G';
18                table['G'] = 'C';
19                table['T'] = 'A';
20                table['U'] = 'A';
21                table['M'] = 'K';
22                table['R'] = 'Y';
23                table['W'] = 'W';
24                table['S'] = 'S';
25                table['Y'] = 'R';
26                table['K'] = 'M';
27                table['V'] = 'B';
28                table['H'] = 'D';
29                table['D'] = 'H';
30                table['B'] = 'V';
31                table['N'] = 'N';
32                table['a'] = 't';
33                table['c'] = 'g';
34                table['g'] = 'c';
35                table['t'] = 'a';
36                table['u'] = 'a';
37                table['m'] = 'k';
38                table['r'] = 'y';
39                table['w'] = 'w';
40                table['s'] = 's';
41                table['y'] = 'r';
42                table['k'] = 'm';
43                table['v'] = 'b';
44                table['h'] = 'd';
45                table['d'] = 'h';
46                table['b'] = 'v';
47                table['n'] = 'n';
48//              table['-'] = '-';
49//              table['.'] = '.';
50//              table['*'] = '*';
51        }
52        return( table[(int)f] );
53}
54
55void sreverse( char *r, char *s )
56{
57        r += strlen( s );
58        *r-- = 0;
59        while( *s )
60                *r-- = creverse( *s++ );
61//              *r-- = ( *s++ );
62}
63
64void gappick_samestring( char *seq )
65{
66        char *aseq = seq;
67
68        for( ; *seq != 0; seq++ )
69        {
70                if( *seq != '-' )
71                        *aseq++ = *seq;
72        }
73        *aseq = 0;
74}
75
76#if 0
77
78static int addlocalhom2( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa, int skip )
79{
80        int pos1, pos2, start1, start2, end1, end2;
81        char *pt1, *pt2;
82        int iscore;
83        int isumscore;
84        int sumoverlap;
85        LocalHom *tmppt;
86        int st;
87        int nlocalhom = 0;
88        pt1 = al1; pt2 = al2;
89        pos1 = off1; pos2 = off2;
90
91        isumscore = 0;
92        sumoverlap = 0;
93
94#if 0
95        fprintf( stderr, "nlocalhom = %d in addlocalhom\n", nlocalhom );
96        fprintf( stderr, "al1 = %s, al2 = %s\n", al1, al2 );
97        fprintf( stderr, "off1 = %d, off2 = %d\n", off1, off2 );
98        fprintf( stderr, "localhopt = %p, skip = %d\n", localhompt, skip );
99        fprintf( stderr, "pt1 = \n%s\n, pt2 = \n%s\n", pt1, pt2 );
100#endif
101
102        if( skip )
103        {
104                while( --skip > 0 ) localhompt = localhompt->next;
105                localhompt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
106                localhompt = localhompt->next;
107//              fprintf( stderr, "tmppt = %p, localhompt = %p\n", tmppt, localhompt );
108        }
109        tmppt = localhompt;
110
111        st = 0;
112        iscore = 0;
113        while( *pt1 != 0 )
114        {
115//              fprintf( stderr, "In in while loop\n" );
116//              fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
117                if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
118                {
119                        end1 = pos1 - 1;
120                        end2 = pos2 - 1;
121
122                        if( nlocalhom++ > 0 )
123                        {
124//                              fprintf( stderr, "reallocating ...\n" );
125                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
126//                              fprintf( stderr, "done\n" );
127                                tmppt = tmppt->next;
128                                tmppt->next = NULL;
129                        }
130                        tmppt->start1 = start1;
131                        tmppt->start2 = start2;
132                        tmppt->end1   = end1  ;
133                        tmppt->end2   = end2  ;
134
135#if 1
136                        isumscore += iscore;
137                        sumoverlap += end2-start2+1;
138#else
139                        tmppt->overlapaa   = end2-start2+1;
140                        tmppt->opt = iscore * 5.8 / 600;
141                        tmppt->overlapaa   = overlapaa;
142                        tmppt->opt = (double)opt;
143#endif
144
145#if 0
146                        fprintf( stderr, "iscore (1)= %d\n", iscore );
147                        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
148                        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
149#endif
150                        iscore = 0;
151                        st = 0;
152                }
153                else if( *pt1 != '-' && *pt2 != '-' )
154                {
155                        if( st == 0 )
156                        {
157                                start1 = pos1; start2 = pos2;
158                                st = 1;
159                        }
160                        iscore += n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]];
161//                      fprintf( stderr, "%c-%c, score(0) = %d\n", *pt1, *pt2, iscore );
162                }
163                if( *pt1++ != '-' ) pos1++;
164                if( *pt2++ != '-' ) pos2++;
165        }
166
167        if( st )
168        {
169                if( nlocalhom++ > 0 )
170                {
171//                      fprintf( stderr, "reallocating ...\n" );
172                        tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
173//                      fprintf( stderr, "done\n" );
174                        tmppt = tmppt->next;
175                        tmppt->next = NULL;
176                }
177                end1 = pos1 - 1;
178                end2 = pos2 - 1;
179                tmppt->start1 = start1;
180                tmppt->start2 = start2;
181                tmppt->end1   = end1  ;
182                tmppt->end2   = end2  ;
183
184#if 1
185                isumscore += iscore;
186                sumoverlap += end2-start2+1;
187#else
188                tmppt->overlapaa   = end2-start2+1;
189                tmppt->opt = (double)iscore * 5.8 / 600;
190                tmppt->overlapaa   = overlapaa;
191                tmppt->opt = (double)opt;
192#endif
193#if 0
194                fprintf( stderr, "score (2)= %d\n", iscore );
195                fprintf( stderr, "al1: %d - %d\n", start1, end1 );
196                fprintf( stderr, "al2: %d - %d\n", start2, end2 );
197#endif
198        }
199
200        for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
201        {
202                tmppt->overlapaa = sumoverlap;
203                tmppt->opt = (double)sumscore * 5.8 / 600 / sumoverlap;
204        }
205        return( nlocalhom );
206}
207
208#endif
209
210
211
212static int addlocalhom_r( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa, int skip )
213{
214        int pos1, pos2, start1, start2, end1, end2;
215        char *pt1, *pt2;
216        double score;
217        double sumscore;
218        int sumoverlap;
219        LocalHom *tmppt = NULL; // by D.Mathog, a guess
220        int st;
221        int nlocalhom = 0;
222        pt1 = al1; pt2 = al2;
223        pos1 = off1; pos2 = off2;
224
225        sumscore = 0.0;
226        sumoverlap = 0;
227        start1 = 0; // by D.Mathog, a guess
228        start2 = 0; // by D.Mathog, a guess
229
230#if 0
231        fprintf( stderr, "nlocalhom = %d in addlocalhom\n", nlocalhom );
232        fprintf( stderr, "al1 = %s, al2 = %s\n", al1, al2 );
233        fprintf( stderr, "off1 = %d, off2 = %d\n", off1, off2 );
234        fprintf( stderr, "localhopt = %p, skip = %d\n", localhompt, skip );
235#endif
236        fprintf( stderr, "pt1 = \n%s\n, pt2 = \n%s\n", pt1, pt2 );
237
238        if( skip )
239        {
240                while( --skip > 0 ) localhompt = localhompt->next;
241                localhompt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
242                localhompt = localhompt->next;
243                fprintf( stderr, "tmppt = %p, localhompt = %p\n", (void *)tmppt, (void *)localhompt );
244        }
245        tmppt = localhompt;
246
247        st = 0;
248        score = 0.0;
249        while( *pt1 != 0 )
250        {
251                fprintf( stderr, "In in while loop\n" );
252                fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
253                if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
254                {
255                        end1 = pos1 - 1;
256                        end2 = pos2 - 1;
257
258                        if( nlocalhom++ > 0 )
259                        {
260//                              fprintf( stderr, "reallocating ...\n" );
261                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
262//                              fprintf( stderr, "done\n" );
263                                tmppt = tmppt->next;
264                                tmppt->next = NULL;
265                        }
266                        tmppt->start1 = start1;
267                        tmppt->start2 = start2;
268                        tmppt->end1   = end1  ;
269                        tmppt->end2   = end2  ;
270
271#if 1
272                        sumscore += score;
273                        sumoverlap += end2-start2+1;
274#else
275                        tmppt->overlapaa   = end2-start2+1;
276                        tmppt->opt = score * 5.8 / 600;
277                        tmppt->overlapaa   = overlapaa;
278                        tmppt->opt = (double)opt;
279#endif
280
281                        fprintf( stderr, "score (1)= %f\n", score );
282                        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
283                        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
284                        score = 0.0;
285                        st = 0;
286                }
287                else if( *pt1 != '-' && *pt2 != '-' )
288                {
289                        if( st == 0 )
290                        {
291                                start1 = pos1; start2 = pos2;
292                                st = 1;
293                        }
294                        score += (double)n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]];
295//                      fprintf( stderr, "%c-%c, score(0) = %f\n", *pt1, *pt2, score );
296                }
297                if( *pt1++ != '-' ) pos1++;
298                if( *pt2++ != '-' ) pos2++;
299        }
300        if( nlocalhom++ > 0 )
301        {
302//              fprintf( stderr, "reallocating ...\n" );
303                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
304//              fprintf( stderr, "done\n" );
305                tmppt = tmppt->next;
306                tmppt->next = NULL;
307        }
308        end1 = pos1 - 1;
309        end2 = pos2 - 1;
310        tmppt->start1 = start1;
311        tmppt->start2 = start2;
312        tmppt->end1   = end1  ;
313        tmppt->end2   = end2  ;
314
315#if 1
316        sumscore += score;
317        sumoverlap += end2-start2+1;
318#else
319        tmppt->overlapaa   = end2-start2+1;
320        tmppt->opt = score * 5.8 / 600;
321        tmppt->overlapaa   = overlapaa;
322        tmppt->opt = (double)opt;
323#endif
324
325        fprintf( stderr, "score (2)= %f\n", score );
326        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
327        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
328
329        for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
330        {
331                tmppt->overlapaa = sumoverlap;
332                tmppt->opt = sumscore * 5.8 / 600 / sumoverlap;
333        }
334        return( nlocalhom );
335}
336void putlocalhom3( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
337{
338        int pos1, pos2, start1, start2, end1, end2;
339        char *pt1, *pt2;
340        double score;
341        double sumscore;
342        int sumoverlap;
343        LocalHom *tmppt;
344        LocalHom *subnosento;
345        int st;
346        int saisho;
347
348        pt1 = al1; pt2 = al2;
349        pos1 = off1; pos2 = off2;
350
351        sumscore = 0.0;
352        sumoverlap = 0;
353        start1 = 0; // by Mathog, a guess
354        start2 = 0; // by Mathog, a guess
355
356        subnosento = localhompt;
357        while( subnosento->next ) subnosento = subnosento->next;
358        tmppt = subnosento;
359
360        saisho = ( localhompt->nokori == 0 );
361
362        fprintf( stderr, "localhompt = %p\n", (void *)localhompt );
363        fprintf( stderr, "tmppt = %p\n", (void *)tmppt );
364        fprintf( stderr, "subnosento = %p\n", (void *)subnosento );
365
366        st = 0;
367        score = 0.0;
368        while( *pt1 != 0 )
369        {
370//              fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
371                if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
372                {
373                        end1 = pos1 - 1;
374                        end2 = pos2 - 1;
375
376                        if( localhompt->nokori++ > 0 )
377                        {
378//                              fprintf( stderr, "reallocating ...\n" );
379                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
380//                              fprintf( stderr, "done\n" );
381                                tmppt = tmppt->next;
382                                tmppt->next = NULL;
383                        }
384                        tmppt->start1 = start1;
385                        tmppt->start2 = start2;
386                        tmppt->end1   = end1  ;
387                        tmppt->end2   = end2  ;
388
389#if 1
390                        if( divpairscore )
391                        {
392                                tmppt->overlapaa   = end2-start2+1;
393                                tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
394                        }
395                        else
396                        {
397                                sumscore += score;
398                                sumoverlap += end2-start2+1;
399                        }
400#else
401                        tmppt->overlapaa   = overlapaa;
402                        tmppt->opt = (double)opt;
403#endif
404
405#if 0
406                        fprintf( stderr, "score (1)= %f\n", score );
407                        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
408                        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
409#endif
410                        score = 0.0;
411                        st = 0;
412                }
413                else if( *pt1 != '-' && *pt2 != '-' )
414                {
415                        if( st == 0 )
416                        {
417                                start1 = pos1; start2 = pos2;
418                                st = 1;
419                        }
420                        score += (double)n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset €Ï€€€é€Ê€€€«€â
421//                      fprintf( stderr, "%c-%c, score(0) = %f\n", *pt1, *pt2, score );
422                }
423                if( *pt1++ != '-' ) pos1++;
424                if( *pt2++ != '-' ) pos2++;
425        }
426        if( *(pt1-1) != '-' && *(pt2-1) != '-'  )
427        {
428                if( localhompt->nokori++ > 0 )
429                {
430//                      fprintf( stderr, "reallocating ...\n" );
431                        tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
432//                      fprintf( stderr, "done\n" );
433                        tmppt = tmppt->next;
434                        tmppt->next = NULL;
435                }
436
437                end1 = pos1 - 1;
438                end2 = pos2 - 1;
439                tmppt->start1 = start1;
440                tmppt->start2 = start2;
441                tmppt->end1   = end1  ;
442                tmppt->end2   = end2  ;
443
444
445#if 1
446                if( divpairscore )
447                {
448                        tmppt->overlapaa   = end2-start2+1;
449                        tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
450                }
451                else
452                {
453                        sumscore += score;
454                        sumoverlap += end2-start2+1;
455                }
456#else
457                tmppt->overlapaa   = overlapaa;
458                tmppt->opt = (double)opt;
459#endif
460
461#if 0
462                fprintf( stderr, "score (2)= %f\n", score );
463                fprintf( stderr, "al1: %d - %d\n", start1, end1 );
464                fprintf( stderr, "al2: %d - %d\n", start2, end2 );
465#endif
466        }
467
468        fprintf( stderr, "sumscore = %f\n", sumscore );
469        if( !divpairscore )
470        {
471
472                if( !saisho ) subnosento = subnosento->next;
473                for( tmppt=subnosento; tmppt; tmppt=tmppt->next )
474                {
475                        tmppt->overlapaa = sumoverlap;
476                        tmppt->opt = sumscore * 5.8 / 600 / sumoverlap;
477                        fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
478                }
479        }
480}
481void putlocalhom_ext( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
482{
483        int pos1, pos2, start1, start2, end1, end2;
484        char *pt1, *pt2;
485        int iscore;
486        int isumscore;
487        int sumoverlap;
488        LocalHom *tmppt = localhompt;
489        int nlocalhom = 0;
490        int st;
491        pt1 = al1; pt2 = al2;
492        pos1 = off1; pos2 = off2;
493
494
495        isumscore = 0;
496        sumoverlap = 0;
497        start1 = 0; // by D.Mathog, a guess
498        start2 = 0; // by D.Mathog, a guess
499
500        st = 0;
501        iscore = 0;
502        while( *pt1 != 0 )
503        {
504//              fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
505                if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
506                {
507                        end1 = pos1 - 1;
508                        end2 = pos2 - 1;
509
510                        if( nlocalhom++ > 0 )
511                        {
512//                              fprintf( stderr, "reallocating ...\n" );
513                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
514//                              fprintf( stderr, "done\n" );
515                                tmppt = tmppt->next;
516                                tmppt->next = NULL;
517                        }
518                        tmppt->start1 = start1;
519                        tmppt->start2 = start2;
520                        tmppt->end1   = end1  ;
521                        tmppt->end2   = end2  ;
522
523#if 1
524                        if( divpairscore )
525                        {
526                                tmppt->overlapaa   = end2-start2+1;
527                                tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
528                        }
529                        else
530                        {
531                                isumscore += iscore;
532                                sumoverlap += end2-start2+1;
533                        }
534#else
535                        tmppt->overlapaa   = overlapaa;
536                        tmppt->opt = (double)opt;
537#endif
538
539#if 0
540                        fprintf( stderr, "iscore (1)= %d\n", iscore );
541                        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
542                        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
543#endif
544                        iscore = 0;
545                        st = 0;
546                }
547                else if( *pt1 != '-' && *pt2 != '-' )
548                {
549                        if( st == 0 )
550                        {
551                                start1 = pos1; start2 = pos2;
552                                st = 1;
553                        }
554                        iscore += n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset €Ï€€€é€Ê€€€«€â
555//                      fprintf( stderr, "%c-%c, iscore(0) = %d\n", *pt1, *pt2, iscore );
556                }
557                if( *pt1++ != '-' ) pos1++;
558                if( *pt2++ != '-' ) pos2++;
559        }
560        if( *(pt1-1) != '-' && *(pt2-1) != '-'  )
561        {
562                if( nlocalhom++ > 0 )
563                {
564//                      fprintf( stderr, "reallocating ...\n" );
565                        tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
566//                      fprintf( stderr, "done\n" );
567                        tmppt = tmppt->next;
568                        tmppt->next = NULL;
569                }
570                end1 = pos1 - 1;
571                end2 = pos2 - 1;
572                tmppt->start1 = start1;
573                tmppt->start2 = start2;
574                tmppt->end1   = end1  ;
575                tmppt->end2   = end2  ;
576       
577#if 1
578                if( divpairscore )
579                {
580                        tmppt->overlapaa   = end2-start2+1;
581                        tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
582                }
583                else
584                {
585                        isumscore += iscore;
586                        sumoverlap += end2-start2+1;
587                }
588#else
589                tmppt->overlapaa   = overlapaa;
590                tmppt->opt = (double)opt;
591#endif
592
593#if 0
594                fprintf( stderr, "iscore (2)= %d\n", iscore );
595                fprintf( stderr, "al1: %d - %d\n", start1, end1 );
596                fprintf( stderr, "al2: %d - %d\n", start2, end2 );
597#endif
598        }
599
600        if( !divpairscore )
601        {
602                for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
603                {
604                        tmppt->overlapaa = sumoverlap;
605//                      tmppt->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
606                        tmppt->opt = (double)600 * 5.8 / 600;
607//                      fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
608                }
609        }
610}
611
612void putlocalhom_str( char *al1, char *al2, double *equiv, double scale, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
613{
614        int posinaln, pos1, pos2, start1, start2, end1, end2;
615        char *pt1, *pt2;
616        int isumscore;
617        int sumoverlap;
618        LocalHom *tmppt = localhompt;
619        int nlocalhom = 0;
620//      int st;
621        pt1 = al1; pt2 = al2;
622        pos1 = off1; pos2 = off2;
623
624        isumscore = 0;
625        sumoverlap = 0;
626        start1 = 0; // by D.Mathog, a guess
627        start2 = 0; // by D.Mathog, a guess
628
629        posinaln = 0;
630        while( *pt1 != 0 )
631        {
632                if( *pt1 != '-' && *pt2 != '-' && equiv[posinaln] > 0.0 )
633                {
634                        start1 = end1 = pos1; start2 = end2 = pos2;
635                        if( nlocalhom++ > 0 )
636                        {
637//                              fprintf( stderr, "reallocating ... (posinaln=%d)\n", posinaln );
638                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
639//                              fprintf( stderr, "done\n" );
640                                tmppt = tmppt->next;
641                                tmppt->next = NULL;
642                        }
643                        tmppt->start1 = start1;
644                        tmppt->start2 = start2;
645                        tmppt->end1   = end1  ;
646                        tmppt->end2   = end2  ;
647
648                        tmppt->overlapaa   = 1;
649//                      tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
650                        tmppt->opt = equiv[posinaln] * scale;
651//                      fprintf( stdout, "*pt1=%c, *pt2=%c, equiv=%f\n", *pt1, *pt2, equiv[posinaln] );
652
653                }
654                if( *pt1++ != '-' ) pos1++;
655                if( *pt2++ != '-' ) pos2++;
656                posinaln++;
657        }
658}
659
660
661void putlocalhom2( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
662{
663        int pos1, pos2, start1, start2, end1, end2;
664        char *pt1, *pt2;
665        int iscore;
666        int isumscore;
667        int sumoverlap;
668        LocalHom *tmppt = localhompt;
669        int nlocalhom = 0;
670        int st;
671        pt1 = al1; pt2 = al2;
672        pos1 = off1; pos2 = off2;
673
674
675        isumscore = 0;
676        sumoverlap = 0;
677        start1 = 0; // by D.Mathog, a guess
678        start2 = 0; // by D.Mathog, a guess
679
680        st = 0;
681        iscore = 0;
682        while( *pt1 != 0 )
683        {
684//              fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
685                if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
686                {
687                        end1 = pos1 - 1;
688                        end2 = pos2 - 1;
689
690                        if( nlocalhom++ > 0 )
691                        {
692//                              fprintf( stderr, "reallocating ...\n" );
693                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
694//                              fprintf( stderr, "done\n" );
695                                tmppt = tmppt->next;
696                                tmppt->next = NULL;
697                        }
698                        tmppt->start1 = start1;
699                        tmppt->start2 = start2;
700                        tmppt->end1   = end1  ;
701                        tmppt->end2   = end2  ;
702
703#if 1
704                        if( divpairscore )
705                        {
706                                tmppt->overlapaa   = end2-start2+1;
707                                tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
708                        }
709                        else
710                        {
711                                isumscore += iscore;
712                                sumoverlap += end2-start2+1;
713                        }
714#else
715                        tmppt->overlapaa   = overlapaa;
716                        tmppt->opt = (double)opt;
717#endif
718
719#if 0
720                        fprintf( stderr, "iscore (1)= %d\n", iscore );
721                        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
722                        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
723#endif
724                        iscore = 0;
725                        st = 0;
726                }
727                else if( *pt1 != '-' && *pt2 != '-' )
728                {
729                        if( st == 0 )
730                        {
731                                start1 = pos1; start2 = pos2;
732                                st = 1;
733                        }
734                        iscore += n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset €Ï€€€é€Ê€€€«€â
735//                      fprintf( stderr, "%c-%c, iscore(0) = %d\n", *pt1, *pt2, iscore );
736                }
737                if( *pt1++ != '-' ) pos1++;
738                if( *pt2++ != '-' ) pos2++;
739        }
740        if( *(pt1-1) != '-' && *(pt2-1) != '-'  )
741        {
742                if( nlocalhom++ > 0 )
743                {
744//                      fprintf( stderr, "reallocating ...\n" );
745                        tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
746//                      fprintf( stderr, "done\n" );
747                        tmppt = tmppt->next;
748                        tmppt->next = NULL;
749                }
750                end1 = pos1 - 1;
751                end2 = pos2 - 1;
752                tmppt->start1 = start1;
753                tmppt->start2 = start2;
754                tmppt->end1   = end1  ;
755                tmppt->end2   = end2  ;
756       
757#if 1
758                if( divpairscore )
759                {
760                        tmppt->overlapaa   = end2-start2+1;
761                        tmppt->opt = (double)iscore / tmppt->overlapaa * 5.8 / 600;
762                }
763                else
764                {
765                        isumscore += iscore;
766                        sumoverlap += end2-start2+1;
767                }
768#else
769                tmppt->overlapaa   = overlapaa;
770                tmppt->opt = (double)opt;
771#endif
772
773#if 0
774                fprintf( stderr, "iscore (2)= %d\n", iscore );
775                fprintf( stderr, "al1: %d - %d\n", start1, end1 );
776                fprintf( stderr, "al2: %d - %d\n", start2, end2 );
777#endif
778        }
779
780        if( !divpairscore )
781        {
782                for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
783                {
784                        tmppt->overlapaa = sumoverlap;
785                        tmppt->opt = (double)isumscore * 5.8 / ( 600 * sumoverlap );
786//                      fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
787                }
788        }
789}
790void putlocalhom( char *al1, char *al2, LocalHom *localhompt, int off1, int off2, int opt, int overlapaa )
791{
792        int pos1, pos2, start1, start2, end1, end2;
793        char *pt1, *pt2;
794        double score;
795        double sumscore;
796        int sumoverlap;
797        LocalHom *tmppt = localhompt;
798        int nlocalhom = 0;
799        int st;
800        pt1 = al1; pt2 = al2;
801        pos1 = off1; pos2 = off2;
802
803
804        sumscore = 0.0;
805        sumoverlap = 0;
806        start1 = 0; // by D.Mathog, a guess
807        start2 = 0; // by D.Mathog, a guess
808
809        st = 0;
810        score = 0.0;
811        while( *pt1 != 0 )
812        {
813//              fprintf( stderr, "pt = %c, %c, st=%d\n", *pt1, *pt2, st );
814                if( st == 1 && ( *pt1 == '-' || *pt2 == '-' ) )
815                {
816                        end1 = pos1 - 1;
817                        end2 = pos2 - 1;
818
819                        if( nlocalhom++ > 0 )
820                        {
821//                              fprintf( stderr, "reallocating ...\n" );
822                                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
823//                              fprintf( stderr, "done\n" );
824                                tmppt = tmppt->next;
825                                tmppt->next = NULL;
826                        }
827                        tmppt->start1 = start1;
828                        tmppt->start2 = start2;
829                        tmppt->end1   = end1  ;
830                        tmppt->end2   = end2  ;
831
832#if 1
833                        if( divpairscore )
834                        {
835                                tmppt->overlapaa   = end2-start2+1;
836                                tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
837                        }
838                        else
839                        {
840                                sumscore += score;
841                                sumoverlap += end2-start2+1;
842                        }
843#else
844                        tmppt->overlapaa   = overlapaa;
845                        tmppt->opt = (double)opt;
846#endif
847
848#if 0
849                        fprintf( stderr, "score (1)= %f\n", score );
850                        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
851                        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
852#endif
853                        score = 0.0;
854                        st = 0;
855                }
856                else if( *pt1 != '-' && *pt2 != '-' )
857                {
858                        if( st == 0 )
859                        {
860                                start1 = pos1; start2 = pos2;
861                                st = 1;
862                        }
863                        score += (double)n_dis[(int)amino_n[(int)*pt1]][(int)amino_n[(int)*pt2]]; // - offset €Ï€€€é€Ê€€€«€â
864//                      fprintf( stderr, "%c-%c, score(0) = %f\n", *pt1, *pt2, score );
865                }
866                if( *pt1++ != '-' ) pos1++;
867                if( *pt2++ != '-' ) pos2++;
868        }
869        if( nlocalhom++ > 0 )
870        {
871//              fprintf( stderr, "reallocating ...\n" );
872                tmppt->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
873//              fprintf( stderr, "done\n" );
874                tmppt = tmppt->next;
875                tmppt->next = NULL;
876        }
877        end1 = pos1 - 1;
878        end2 = pos2 - 1;
879        tmppt->start1 = start1;
880        tmppt->start2 = start2;
881        tmppt->end1   = end1  ;
882        tmppt->end2   = end2  ;
883
884#if 1
885        if( divpairscore )
886        {
887                tmppt->overlapaa   = end2-start2+1;
888                tmppt->opt = score / tmppt->overlapaa * 5.8 / 600;
889        }
890        else
891        {
892                sumscore += score;
893                sumoverlap += end2-start2+1;
894        }
895#else
896        tmppt->overlapaa   = overlapaa;
897        tmppt->opt = (double)opt;
898#endif
899
900#if 0
901        fprintf( stderr, "score (2)= %f\n", score );
902        fprintf( stderr, "al1: %d - %d\n", start1, end1 );
903        fprintf( stderr, "al2: %d - %d\n", start2, end2 );
904#endif
905
906        if( !divpairscore )
907        {
908                for( tmppt=localhompt; tmppt; tmppt=tmppt->next )
909                {
910                        tmppt->overlapaa = sumoverlap;
911                        tmppt->opt = sumscore * 5.8 / 600 / sumoverlap;
912//                      fprintf( stderr, "tmpptr->opt = %f\n", tmppt->opt );
913                }
914        }
915}
916
917char *cutal( char *al, int al_display_start, int start, int end )
918{
919        int pos;
920        char *pt = al;
921        char *val = NULL;
922
923        pos = al_display_start;
924        do
925        {
926                if( start == pos ) val = pt;
927                if( end == pos ) break;
928//              fprintf( stderr, "pos=%d, *pt=%c, val=%p\n", pos, *pt, val );
929                if( *pt != '-' ) pos++;
930        } while( *pt++ != 0 );
931        *(pt+1) = 0;
932        return( val );
933}
934
935void ErrorExit( char *message )
936{
937        fprintf( stderr, "%s\n", message );
938        exit( 1 );
939}
940
941void strncpy_caseC( char *str1, char *str2, int len )
942{
943        if( dorp == 'd' && upperCase > 0 ) 
944        {
945                while( len-- )
946                        *str1++ = toupper( *str2++ );
947        }
948        else strncpy( str1, str2, len );
949}
950       
951void seqUpper( int nseq, char **seq ) /* not used */
952{
953        int i, j, len;
954        for( i=0; i<nseq; i++ ) 
955        {
956                len = strlen( seq[i] );
957                for( j=0; j<len; j++ ) 
958                        seq[i][j] = toupper( seq[i][j] );
959        }
960}
961
962void seqLower( int nseq, char **seq )
963{
964        int i, j, len;
965        for( i=0; i<nseq; i++ ) 
966        {
967                len = strlen( seq[i] );
968                for( j=0; j<len; j++ ) 
969                        seq[i][j] = tolower( seq[i][j] );
970        }
971}
972
973int getaline_fp_eof( char *s, int l, FILE *fp )  /* end of file -> return 1 */
974{
975    int c, i = 0 ;
976    int noteofflag = 0;
977    for( i=0; i<l && ( noteofflag = ( (c=getc(fp)) != EOF ) ) && c != '\n'; i++ ) 
978        *s++ = c;
979    *s = '\0' ;
980     return( !noteofflag );
981}
982
983int getaline_fp_eof_new(s, l, fp)  /* end of file -> return 1 */
984char    s[] ; int l ; FILE *fp ;
985{
986        int c = 0, i = 0 ;
987                int noteofflag = 0;
988
989                if( feof( fp ) ) return( 1 );
990
991                for( i=0; i<l && ( noteofflag = ( (c=getc(fp)) != EOF ) ) && c != '\n'; i++ ) 
992        { *s++ = c ; }
993        *s = '\0' ;
994                if( c != '\n' && c != EOF ) while( getc(fp) != '\n' )
995                        ;
996                return( !noteofflag );
997}
998
999int myfgets(s, l, fp)  /* l°ÊŸå€Ï¡¢¹ÔËö€Þ€ÇÆÉ€ßÈô€Ð€¹ */
1000char    s[] ; int l ; FILE *fp ;
1001{
1002        int     c = 0, i = 0 ;
1003
1004                if( feof( fp ) ) return( 1 );
1005
1006                for( i=0; i<l && ( c=getc( fp ) ) != '\n'; i++ ) 
1007                *s++ = c;
1008        *s = '\0' ;
1009                if( c != '\n' ) 
1010                        while( getc(fp) != '\n' )
1011                                ;
1012                return( 0 );
1013}
1014
1015float input_new( FILE *fp, int d )
1016{
1017        char mojiretsu[10];
1018        int i, c;
1019
1020        c = getc( fp );
1021        if( c != '\n' )
1022                ungetc( c, fp );
1023
1024        for( i=0; i<d; i++ )
1025                mojiretsu[i] = getc( fp );
1026        mojiretsu[i] = 0;
1027
1028        return( atof( mojiretsu ) );
1029}
1030
1031void fgets_IGNORE_RESULT(char *buf, int size, FILE *fp) {
1032    char *result = fgets(buf, size, fp);
1033    if (!result) {
1034        fprintf(stderr, "fgets() fails, but we ignore that fact[1].\n");
1035    }
1036}
1037
1038void PreRead( FILE *fp, int *locnjob, int *locnlenmax )
1039{
1040    int i, nleni;
1041        char b[B];
1042
1043        fgets_IGNORE_RESULT( b, B-1, fp ); *locnjob = atoi( b );
1044        *locnlenmax = 0;
1045        i=0; 
1046        while( i<*locnjob )
1047        {
1048                fgets_IGNORE_RESULT( b, B-1, fp );
1049                if( !strncmp( b, "=", 1 ) )
1050                {
1051                        i++;
1052                        fgets_IGNORE_RESULT( b, B-1, fp ); nleni = atoi( b );
1053                        if( nleni > *locnlenmax ) *locnlenmax = nleni;
1054                }
1055        }
1056        if( *locnlenmax > N )
1057        {
1058                fprintf( stderr, "TOO LONG SEQUENCE!\n" );
1059                exit( 1 );
1060        }
1061        if( njob > M  ) 
1062        {
1063                fprintf( stderr, "TOO MANY SEQUENCE!\n" );
1064                fprintf( stderr, "%d > %d\n", njob, M );
1065                exit( 1 );
1066        }
1067}       
1068
1069int allSpace( char *str )
1070{
1071        int value = 1;
1072        while( *str ) value *= ( !isdigit( *str++ ) );
1073        return( value );
1074}
1075       
1076void Read( char name[M][B], int nlen[M], char **seq )
1077{
1078        extern void FRead( FILE *x, char y[M][B], int z[M], char **w );
1079        FRead( stdin, name, nlen, seq );
1080}
1081
1082
1083void FRead( FILE *fp, char name[][B], int nlen[], char **seq )
1084{
1085        int i, j; 
1086        char b[B];
1087
1088        fgets_IGNORE_RESULT( b, B-1, fp );
1089#if DEBUG
1090        fprintf( stderr, "b = %s\n", b );
1091#endif
1092
1093    if( strstr( b, "onnet" ) ) scoremtx = 1;
1094    else if( strstr( b, "DnA" ) ) 
1095        {
1096                scoremtx = -1; 
1097                upperCase = -1;
1098        }
1099    else if( strstr( b, "dna" ) ) 
1100        {
1101                scoremtx = -1; 
1102                upperCase = 0;
1103        }
1104        else if( strstr( b, "DNA" ) )
1105        {
1106                scoremtx = -1; 
1107                upperCase = 1;
1108        }
1109    else if( strstr( b, "M-Y" ) || strstr( b, "iyata" ) ) scoremtx = 2; 
1110    else scoremtx = 0;
1111#if DEBUG
1112        fprintf( stderr, " %s->scoremtx = %d\n", b, scoremtx );
1113#endif
1114
1115        geta2 = GETA2;
1116
1117#if 0
1118        if( strlen( b ) >=25 )
1119        {
1120                b[25] = 0;
1121        #if DEBUG
1122                fprintf( stderr, "kimuraR = %s\n", b+20 );
1123        #endif
1124                kimuraR = atoi( b+20 );
1125
1126                if( kimuraR < 0 || 20 < kimuraR ) ErrorExit( "Illeagal kimuraR value.\n" );
1127                if( allSpace( b+20 ) ) kimuraR = NOTSPECIFIED;
1128        }
1129        else kimuraR = NOTSPECIFIED;
1130        #if DEBUG
1131        fprintf( stderr, "kimuraR = %d\n", kimuraR );
1132        #endif
1133
1134        if( strlen( b ) >=20 )
1135        {
1136                b[20] = 0;
1137        #if DEBUG
1138                fprintf( stderr, "pamN = %s\n", b+15 );
1139        #endif
1140                pamN = atoi( b+15 );
1141                if( pamN < 0 || 400 < pamN ) ErrorExit( "Illeagal pam value.\n" );
1142                if( allSpace( b+15 ) ) pamN = NOTSPECIFIED;
1143        }
1144        else pamN = NOTSPECIFIED;
1145
1146        if( strlen( b ) >= 15 )
1147        {
1148                b[15] = 0;
1149        #if DEBUG
1150                fprintf( stderr, "poffset = %s\n", b+10 );
1151        #endif
1152                poffset = atoi( b+10 );
1153                if( poffset > 500 ) ErrorExit( "Illegal extending gap ppenalty\n" );
1154                if( allSpace( b+10 ) ) poffset = NOTSPECIFIED;
1155        }
1156        else poffset = NOTSPECIFIED;
1157
1158        if( strlen( b ) >= 10 )
1159        {
1160                b[10] = 0;
1161        #if DEBUG
1162                fprintf( stderr, "ppenalty = %s\n", b+5 );
1163        #endif
1164                ppenalty = atoi( b+5 );
1165                if( ppenalty > 0 ) ErrorExit( "Illegal opening gap ppenalty\n" );
1166                if( allSpace( b+5 ) ) ppenalty = NOTSPECIFIED;
1167        }
1168        else ppenalty = NOTSPECIFIED;
1169#endif
1170
1171        for( i=0; i<njob; i++ )
1172        {
1173                getaline_fp_eof_new( b, B-1, fp );
1174                strcpy( name[i], b );
1175#if DEBUG
1176                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1177#endif
1178                fgets_IGNORE_RESULT( b, B-1, fp ); nlen[i] = atoi( b ); /* seq i no nagasa */
1179                seq[i][0] = 0;
1180                if( nlen[i] ) for( j=0; j <= (nlen[i]-1)/C; j++ )
1181                {
1182                        getaline_fp_eof_new( b, B-1, fp );
1183                        /*      b[C] = 0;  */
1184                        strcat( seq[i], b );
1185                } 
1186                seq[i][nlen[i]] = 0;
1187        }
1188        if( scoremtx == -1 && upperCase != -1 ) seqLower( njob, seq );
1189}
1190
1191
1192static int countKUorWA( FILE *fp )
1193{
1194        int value;
1195        int c, b;
1196
1197        value= 0;
1198        b = '\n';
1199        while( ( c = getc( fp ) ) != EOF )
1200        {
1201                if( b == '\n' && ( c == '>' ) )
1202                        value++;
1203                b = c;
1204        }
1205        rewind( fp );
1206        return( value );
1207}
1208
1209void searchKUorWA( FILE *fp )
1210{
1211        int c, b;
1212        b = '\n';
1213        while( !( ( ( c = getc( fp ) ) == '>' || c == EOF ) && b == '\n' ) )
1214                b = c;
1215        ungetc( c, fp );
1216}
1217
1218static int onlyGraph( char *str )
1219{
1220        char tmp;
1221        char *res = str;
1222        char *bk = str;
1223
1224//      while( (tmp=*str++) ) if( isgraph( tmp ) ) *res++ = tmp;
1225        while( (tmp=*str++) ) 
1226        {
1227                if( 0x20 < tmp && tmp < 0x7f ) *res++ = tmp;
1228                if( tmp == '>' )
1229                {
1230                        fprintf( stderr, "========================================================\n" );
1231                        fprintf( stderr, "========================================================\n" );
1232                        fprintf( stderr, "=== \n" );
1233                        fprintf( stderr, "=== ERROR!! \n" );
1234                        fprintf( stderr, "=== In the '--anysymbol' and '--preservecase' modes, \n" );
1235                        fprintf( stderr, "=== '>' in sequence is unacceptable.\n" );
1236                        fprintf( stderr, "=== \n" );
1237                        fprintf( stderr, "========================================================\n" );
1238                        fprintf( stderr, "========================================================\n" );
1239                        exit( 1 );
1240                }
1241        }
1242        *res = 0;
1243        return( res - bk );
1244}
1245
1246static int onlyAlpha_lower( char *str )
1247{
1248        char tmp;
1249        char *res = str;
1250        char *bk = str;
1251
1252        while( (tmp=*str++) )
1253                if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1254                        *res++ = tolower( tmp );
1255        *res = 0;
1256        return( res - bk );
1257}
1258static int onlyAlpha_upper( char *str )
1259{
1260        char tmp;
1261        char *res = str;
1262        char *bk = str;
1263
1264        while( (tmp=*str++) )
1265                if( isalpha( tmp ) || tmp == '-' || tmp == '*' || tmp == '.' )
1266                        *res++ = toupper( tmp );
1267        *res = 0;
1268        return( res - bk );
1269}
1270
1271void kake2hiku( char *str )
1272{
1273        do
1274                if( *str == '*' ) *str = '-';
1275        while( *str++ );
1276}
1277
1278char *load1SeqWithoutName_realloc_casepreserve( FILE *fpp )
1279{
1280        int c, b;
1281        char *cbuf;
1282        int size = N;
1283        char *val;
1284
1285        val = malloc( (size+1) * sizeof( char ) );
1286        cbuf = val;
1287
1288        b = '\n';
1289        while( ( c = getc( fpp ) ) != EOF &&           
1290          !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1291        {
1292                *cbuf++ = (char)c;  /* Ĺ€¹€®€Æ€â€·€é€Ê€€ */
1293                if( cbuf - val == size )
1294                {
1295                        size += N;
1296                        fprintf( stderr, "reallocating...\n" );
1297                        val = (char *)realloc( val, (size+1) * sizeof( char ) );
1298                        if( !val )
1299                        {
1300                                fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1301                                exit( 1 );
1302                        }
1303                        fprintf( stderr, "done.\n" );
1304                        cbuf = val + size-N;
1305                }
1306                b = c;
1307        }
1308        ungetc( c, fpp );
1309        *cbuf = 0;
1310        onlyGraph( val );
1311//      kake2hiku( val );
1312        return( val );
1313}
1314
1315char *load1SeqWithoutName_realloc( FILE *fpp )
1316{
1317        int c, b;
1318        char *cbuf;
1319        int size = N;
1320        char *val;
1321
1322        val = malloc( (size+1) * sizeof( char ) );
1323        cbuf = val;
1324
1325        b = '\n';
1326        while( ( c = getc( fpp ) ) != EOF &&           
1327          !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1328        {
1329                *cbuf++ = (char)c;  /* Ĺ€¹€®€Æ€â€·€é€Ê€€ */
1330                if( cbuf - val == size )
1331                {
1332                        size += N;
1333                        fprintf( stderr, "reallocating...\n" );
1334                        val = (char *)realloc( val, (size+1) * sizeof( char ) );
1335                        if( !val )
1336                        {
1337                                fprintf( stderr, "Allocation error in load1SeqWithoutName_realloc \n" );
1338                                exit( 1 );
1339                        }
1340                        fprintf( stderr, "done.\n" );
1341                        cbuf = val + size-N;
1342                }
1343                b = c;
1344        }
1345        ungetc( c, fpp );
1346        *cbuf = 0;
1347        if( dorp == 'd' )
1348                onlyAlpha_lower( val );
1349        else
1350                onlyAlpha_upper( val );
1351        kake2hiku( val );
1352        return( val );
1353}
1354
1355int load1SeqWithoutName_new( FILE *fpp, char *cbuf )
1356{
1357        int c, b;
1358        char *bk = cbuf;
1359
1360        b = '\n';
1361        while( ( c = getc( fpp ) ) != EOF &&                    /* by T. Nishiyama */
1362          !( ( c == '>' || c == '(' || c == EOF ) && b == '\n' ) )
1363        {
1364                *cbuf++ = (char)c;  /* Ĺ€¹€®€Æ€â€·€é€Ê€€ */
1365                b = c;
1366        }
1367        ungetc( c, fpp );
1368        *cbuf = 0;
1369        if( dorp == 'd' )
1370                onlyAlpha_lower( bk );
1371        else
1372                onlyAlpha_upper( bk );
1373        kake2hiku( bk );
1374        return( 0 );
1375}
1376
1377
1378void readDataforgaln( FILE *fp, char **name, int *nlen, char **seq )
1379{
1380        int i; 
1381        static char *tmpseq = NULL;
1382
1383#if 0
1384        if( !tmpseq )
1385        {
1386                tmpseq = AllocateCharVec( N );
1387        }
1388#endif
1389
1390        rewind( fp );
1391        searchKUorWA( fp );
1392
1393        for( i=0; i<njob; i++ )
1394        {
1395                name[i][0] = '='; getc( fp ); 
1396#if 0
1397                fgets( name[i]+1, B-2, fp );
1398                j = strlen( name[i] );
1399                if( name[i][j-1] != '\n' )
1400                        ErrorExit( "Too long name\n" );
1401                name[i][j-1] = 0;
1402#else
1403                myfgets( name[i]+1, B-2, fp ); 
1404#endif
1405#if 0
1406                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1407#endif
1408                tmpseq = load1SeqWithoutName_realloc( fp );
1409                strcpy( seq[i], tmpseq );
1410                nlen[i] = strlen( seq[i] );
1411                free( tmpseq );
1412        }
1413        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1414#if 0
1415        free( tmpseq );
1416#endif
1417}
1418
1419void readData_varlen( FILE *fp, char **name, int *nlen, char **seq )
1420{
1421        int i; 
1422        static char *tmpseq = NULL;
1423
1424        rewind( fp );
1425        searchKUorWA( fp );
1426
1427        for( i=0; i<njob; i++ )
1428        {
1429                name[i][0] = '='; getc( fp ); 
1430#if 0
1431                fgets( name[i]+1, B-2, fp );
1432                j = strlen( name[i] );
1433                if( name[i][j-1] != '\n' )
1434                        ErrorExit( "Too long name\n" );
1435                name[i][j-1] = 0;
1436#else
1437                myfgets( name[i]+1, B-2, fp ); 
1438#endif
1439#if 0
1440                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1441#endif
1442                tmpseq = load1SeqWithoutName_realloc( fp );
1443                nlen[i] = strlen( tmpseq );
1444//              fprintf( stderr, "nlen[%d] = %d\n", i+1, nlen[i] );
1445                seq[i] = calloc( nlen[i]+1, sizeof( char ) );
1446                strcpy( seq[i], tmpseq );
1447                free( tmpseq );
1448        }
1449        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1450#if 0
1451        free( tmpseq );
1452#endif
1453}
1454
1455void readData_pointer2( FILE *fp, int nseq, char **name, int *nlen, char **seq )
1456{
1457        int i; 
1458        static char *tmpseq = NULL;
1459
1460#if 0
1461        if( !tmpseq )
1462        {
1463                tmpseq = AllocateCharVec( N );
1464        }
1465#endif
1466
1467        rewind( fp );
1468        searchKUorWA( fp );
1469
1470        for( i=0; i<nseq; i++ )
1471        {
1472                name[i][0] = '='; getc( fp ); 
1473#if 0
1474                fgets( name[i]+1, B-2, fp );
1475                j = strlen( name[i] );
1476                if( name[i][j-1] != '\n' )
1477                        ErrorExit( "Too long name\n" );
1478                name[i][j-1] = 0;
1479#else
1480                myfgets( name[i]+1, B-2, fp ); 
1481#endif
1482#if 0
1483                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1484#endif
1485                tmpseq = load1SeqWithoutName_realloc( fp );
1486                strcpy( seq[i], tmpseq );
1487                free( tmpseq );
1488                nlen[i] = strlen( seq[i] );
1489        }
1490        if( dorp == 'd' && upperCase != -1 ) seqLower( nseq, seq );
1491#if 0
1492        free( tmpseq );
1493#endif
1494        if( outnumber )
1495        {
1496                char *namebuf;
1497                char *cptr;
1498                namebuf = calloc( B+100, sizeof( char ) );
1499                for( i=0; i<nseq; i++ )
1500                {
1501                        namebuf[0] = '=';
1502                        cptr = strstr( name[i], "_numo_e_" );
1503                        if( cptr )
1504                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1505                        else
1506                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1507                        strncpy( name[i], namebuf, B );
1508                        name[i][B-1] = 0;
1509                }
1510                free( namebuf );
1511//              exit( 1 );
1512        }
1513}
1514
1515
1516void readData_pointer_casepreserve( FILE *fp, char **name, int *nlen, char **seq )
1517{
1518        int i; 
1519        static char *tmpseq = NULL;
1520
1521#if 0
1522        if( !tmpseq )
1523        {
1524                tmpseq = AllocateCharVec( N );
1525        }
1526#endif
1527
1528        rewind( fp );
1529        searchKUorWA( fp );
1530
1531        for( i=0; i<njob; i++ )
1532        {
1533                name[i][0] = '='; getc( fp ); 
1534#if 0
1535                fgets( name[i]+1, B-2, fp );
1536                j = strlen( name[i] );
1537                if( name[i][j-1] != '\n' )
1538                        ErrorExit( "Too long name\n" );
1539                name[i][j-1] = 0;
1540#else
1541                myfgets( name[i]+1, B-2, fp ); 
1542#endif
1543#if 0
1544                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1545#endif
1546                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1547                strcpy( seq[i], tmpseq );
1548                free( tmpseq );
1549                nlen[i] = strlen( seq[i] );
1550        }
1551}
1552
1553
1554void readData_pointer( FILE *fp, char **name, int *nlen, char **seq )
1555{
1556        int i; 
1557        static char *tmpseq = NULL;
1558
1559#if 0
1560        if( !tmpseq )
1561        {
1562                tmpseq = AllocateCharVec( N );
1563        }
1564#endif
1565
1566        rewind( fp );
1567        searchKUorWA( fp );
1568
1569        for( i=0; i<njob; i++ )
1570        {
1571                name[i][0] = '='; getc( fp ); 
1572#if 0
1573                fgets( name[i]+1, B-2, fp );
1574                j = strlen( name[i] );
1575                if( name[i][j-1] != '\n' )
1576                        ErrorExit( "Too long name\n" );
1577                name[i][j-1] = 0;
1578#else
1579                myfgets( name[i]+1, B-2, fp ); 
1580#endif
1581#if 0
1582                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1583#endif
1584                tmpseq = load1SeqWithoutName_realloc( fp );
1585                strcpy( seq[i], tmpseq );
1586                free( tmpseq );
1587                nlen[i] = strlen( seq[i] );
1588        }
1589        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1590#if 0
1591        free( tmpseq );
1592#endif
1593        if( outnumber )
1594        {
1595                char *namebuf;
1596                char *cptr;
1597                namebuf = calloc( B+100, sizeof( char ) );
1598                for( i=0; i<njob; i++ )
1599                {
1600                        namebuf[0] = '=';
1601                        cptr = strstr( name[i], "_numo_e_" );
1602                        if( cptr )
1603                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, cptr+8 );
1604                        else
1605                                sprintf( namebuf+1, "_numo_s_%08d_numo_e_%s", i+1, name[i]+1 );
1606                        strncpy( name[i], namebuf, B );
1607                        name[i][B-1] = 0;
1608                }
1609                free( namebuf );
1610//              exit( 1 );
1611        }
1612}
1613
1614void readData( FILE *fp, char name[][B], int nlen[], char **seq )
1615{
1616        int i; 
1617        static char *tmpseq = NULL;
1618
1619#if 0
1620        if( !tmpseq )
1621        {
1622                tmpseq = AllocateCharVec( N );
1623        }
1624#endif
1625
1626        rewind( fp );
1627        searchKUorWA( fp );
1628
1629        for( i=0; i<njob; i++ )
1630        {
1631                name[i][0] = '='; getc( fp ); 
1632#if 0
1633                fgets( name[i]+1, B-2, fp );
1634                j = strlen( name[i] );
1635                if( name[i][j-1] != '\n' )
1636                        ErrorExit( "Too long name\n" );
1637                name[i][j-1] = 0;
1638#else
1639                myfgets( name[i]+1, B-2, fp ); 
1640#endif
1641#if 0
1642                fprintf( stderr, "name[%d] = %s\n", i, name[i] );
1643#endif
1644                tmpseq = load1SeqWithoutName_realloc( fp );
1645                strcpy( seq[i], tmpseq );
1646                nlen[i] = strlen( seq[i] );
1647                free( tmpseq );
1648        }
1649        if( dorp == 'd' && upperCase != -1 ) seqLower( njob, seq );
1650#if 0
1651        free( tmpseq );
1652#endif
1653}
1654
1655void cutAlignment( FILE *fp, int **regtable, char **revtable, int *outtable, char **name, char **outseq )
1656{
1657        int i, j; 
1658        int outlen;
1659        static char *tmpseq = NULL;
1660        static char *dumname = NULL;
1661        char *fs, *rs;
1662        int npos, lpos;
1663        int startpos, endpos, seqlen;
1664
1665        if( dumname == NULL )
1666        {
1667                dumname = AllocateCharVec( N );
1668        }
1669
1670        rewind( fp );
1671        searchKUorWA( fp );
1672
1673
1674        npos = 0;
1675        for( i=0; i<njob; i++ )
1676        {
1677                dumname[0] = '>'; getc( fp ); 
1678                myfgets( dumname+1, B-1, fp ); 
1679                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1680
1681                if( outtable[i] )
1682                {
1683//                      putc( '>', stdout );
1684//                      puts( dumname+1 );
1685
1686       
1687                        strncat( name[npos], dumname, B-1 );
1688                        name[npos][B-1] = 0;
1689       
1690                        if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1691                        seqlen = strlen( tmpseq );
1692                        lpos = 0;
1693                        for( j=0; j<5; j++ )
1694                        {
1695                                if( regtable[0][j*2] == -1 && regtable[0][j*2+1] == -1 ) continue;
1696
1697                                startpos = regtable[0][j*2];
1698                                endpos   = regtable[0][j*2+1];
1699                                if( startpos > endpos )
1700                                {
1701                                        endpos   = regtable[0][j*2];
1702                                        startpos = regtable[0][j*2+1];
1703                                }
1704
1705                                if( startpos < 0 ) startpos = 0;
1706                                if( endpos   < 0 ) endpos   = 0;
1707                                if( endpos   >= seqlen ) endpos   = seqlen-1;
1708                                if( startpos >= seqlen ) startpos = seqlen-1;
1709
1710//                              fprintf( stderr, "startpos = %d, endpos = %d\n", startpos, endpos );
1711
1712                                outlen = endpos - startpos+1;
1713                                if( revtable[0][j] == 'f' )
1714                                {
1715//                                      fprintf( stderr, "regtable[%d][st] = %d\n", i, regtable[0][j*2+0] );
1716//                                      fprintf( stderr, "regtable[%d][en] = %d\n", i, regtable[0][j*2+1] );
1717//                                      fprintf( stderr, "outlen = %d\n", outlen );
1718//                                      fprintf( stdout, "%.*s\n", outlen, tmpseq+regtable[0][j*2] );
1719                                        strncpy( outseq[npos] + lpos, tmpseq+startpos, outlen );
1720                                        lpos += outlen;
1721                                }
1722                                else
1723                                {
1724                                        fs = AllocateCharVec( outlen+1 );
1725                                        rs = AllocateCharVec( outlen+1 );
1726
1727                                        fs[outlen] = 0;
1728                                        strncpy( fs, tmpseq+startpos, outlen );
1729                                        sreverse( rs, fs );
1730//                                      fprintf( stdout, "%s\n", rs );
1731                                        strncpy( outseq[npos] + lpos, rs, outlen );
1732                                        lpos += outlen;
1733                                        free( fs );
1734                                        free( rs );
1735                                }
1736                                outseq[npos][lpos] = 0;
1737                        }
1738                        npos++;
1739                }
1740                free( tmpseq );
1741        }
1742}
1743
1744void cutData( FILE *fp, int **regtable, char **revtable, int *outtable )
1745{
1746        int i, j; 
1747        int outlen, seqlen, startpos, endpos;
1748        static char *tmpseq = NULL;
1749        static char *dumname = NULL;
1750        char *fs, *rs;
1751
1752        if( dumname == NULL )
1753        {
1754                dumname = AllocateCharVec( N );
1755        }
1756
1757        rewind( fp );
1758        searchKUorWA( fp );
1759
1760        for( i=0; i<njob; i++ )
1761        {
1762                dumname[0] = '='; getc( fp ); 
1763                myfgets( dumname+1, B-2, fp ); 
1764                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1765
1766                if( outtable[i] )
1767                {
1768                        gappick_samestring( tmpseq );
1769                        putc( '>', stdout );
1770                        puts( dumname+1 );
1771
1772                        seqlen = strlen( tmpseq );
1773
1774                        if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1775                        if( outtable[i] == 2 )
1776                        {
1777                                startpos = 0;
1778                                endpos   = seqlen-1;
1779                                outlen = endpos - startpos + 1;
1780                                fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1781                        }
1782                        else
1783                        {
1784                                for( j=0; j<5; j++ )
1785                                {
1786                                        if( regtable[i][j*2] == -1 && regtable[i][j*2+1] == -1 ) continue;
1787       
1788                                        startpos = regtable[i][j*2];
1789                                        endpos   = regtable[i][j*2+1];
1790       
1791                                        if( startpos > endpos )
1792                                        {
1793                                                endpos   = regtable[i][j*2];
1794                                                startpos = regtable[i][j*2+1];
1795                                        }
1796       
1797                                        if( startpos < 0 ) startpos = 0;
1798                                        if( endpos   < 0 ) endpos   = 0;
1799                                        if( endpos   >= seqlen ) endpos   = seqlen-1;
1800                                        if( startpos >= seqlen ) startpos = seqlen-1;
1801       
1802                                        outlen = endpos - startpos + 1;
1803                                        if( revtable[i][j] == 'f' )
1804                                        {
1805                                                fprintf( stderr, "startpos = %d\n", startpos );
1806                                                fprintf( stderr, "endpos   = %d\n", endpos );
1807                                                fprintf( stderr, "outlen = %d\n", outlen );
1808                                                fprintf( stdout, "%.*s\n", outlen, tmpseq+startpos );
1809                                        }
1810                                        else
1811                                        {
1812                                                fs = AllocateCharVec( outlen+1 );
1813                                                rs = AllocateCharVec( outlen+1 );
1814       
1815                                                fs[outlen] = 0;
1816                                                strncpy( fs, tmpseq+startpos, outlen );
1817                                                sreverse( rs, fs );
1818                                                fprintf( stdout, "%s\n", rs );
1819                                                free( fs );
1820                                                free( rs );
1821                                        }
1822                                }
1823                        }
1824                }
1825                free( tmpseq );
1826        }
1827}
1828
1829void catData( FILE *fp )
1830{
1831        int i; 
1832        static char *tmpseq = NULL;
1833        static char *dumname = NULL;
1834//      char *cptr;
1835
1836        if( dumname == NULL )
1837        {
1838                dumname = AllocateCharVec( N );
1839        }
1840
1841        rewind( fp );
1842        searchKUorWA( fp );
1843
1844        for( i=0; i<njob; i++ )
1845        {
1846                dumname[0] = '='; getc( fp ); 
1847                myfgets( dumname+1, B-2, fp ); 
1848                if( outnumber )
1849                {
1850                        fprintf( stdout, ">_numo_s_%08d_numo_e_", i+1 );
1851                }
1852                else
1853                {
1854                        putc( '>', stdout );
1855                }
1856                puts( dumname+1 );
1857                tmpseq = load1SeqWithoutName_realloc( fp );
1858                if( dorp == 'd' && upperCase != -1 ) seqLower( 1, &tmpseq );
1859                puts( tmpseq );
1860                free( tmpseq );
1861        }
1862}
1863
1864int countATGC( char *s, int *total )
1865{
1866        int nATGC;
1867        int nChar;
1868        char c;
1869        nATGC = nChar = 0;
1870
1871        if( *s == 0 ) 
1872        {
1873                *total = 0;
1874                return( 0 );
1875        }
1876
1877        do
1878        {
1879                c = tolower( *s );
1880                if( isalpha( c ) )
1881                {
1882                        nChar++;
1883                        if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1884                                nATGC++;
1885                }
1886        }
1887        while( *++s );
1888
1889        *total = nChar;
1890        return( nATGC );
1891}
1892
1893double countATGCbk( char *s )
1894{
1895        int nATGC;
1896        int nChar;
1897        char c;
1898        nATGC = nChar = 0;
1899
1900        do
1901        {
1902                c = tolower( *s );
1903                if( isalpha( c ) )
1904                {
1905                        nChar++;
1906                        if( c == 'a' || c == 't' || c == 'g' || c == 'c' || c == 'u' || c == 'n' )
1907                                nATGC++;
1908                }
1909        }
1910        while( *++s );
1911        return( (double)nATGC / nChar );
1912}
1913
1914
1915int countnogaplen( char *seq )
1916{
1917        int val = 0;
1918        while( *seq )
1919                if( *seq++ != '-' ) val++;
1920        return( val );
1921}
1922
1923int countnormalletters( char *seq, char *ref )
1924{
1925        int val = 0;
1926        while( *seq )
1927                if( strchr( ref, *seq++ ) ) val++;
1928        return( val );
1929}
1930
1931void getnumlen_casepreserve( FILE *fp, int *nlenminpt )
1932{
1933        int total;
1934        int nsite = 0;
1935        int atgcnum;
1936        int i, tmp;
1937        char *tmpseq, *tmpname;
1938        double atgcfreq;
1939        tmpname = AllocateCharVec( N );
1940        njob = countKUorWA( fp );
1941        searchKUorWA( fp );
1942        nlenmax = 0;
1943        *nlenminpt = 99999999;
1944        atgcnum = 0;
1945        total = 0;
1946        for( i=0; i<njob; i++ )
1947        {
1948                myfgets( tmpname, N-1, fp );
1949                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
1950                tmp = strlen( tmpseq );
1951                if( tmp > nlenmax ) nlenmax  = tmp;
1952                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1953                atgcnum += countATGC( tmpseq, &nsite );
1954                total += nsite;
1955                free( tmpseq );
1956        }
1957        free( tmpname );
1958        atgcfreq = (double)atgcnum / total;
1959//      fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
1960        if( dorp == NOTSPECIFIED )
1961        {
1962                if( atgcfreq > 0.75 )   
1963                {
1964                        dorp = 'd';
1965                        upperCase = -1;
1966                }
1967                else                 
1968                {
1969                        dorp = 'p';
1970                        upperCase = 0;
1971                }
1972        }
1973}
1974
1975void getnumlen_nogap( FILE *fp, int *nlenminpt )
1976{
1977        int total;
1978        int nsite = 0;
1979        int atgcnum;
1980        int i, tmp;
1981        char *tmpseq, *tmpname;
1982        double atgcfreq;
1983        tmpname = AllocateCharVec( N );
1984        njob = countKUorWA( fp );
1985        searchKUorWA( fp );
1986        nlenmax = 0;
1987        *nlenminpt = 99999999;
1988        atgcnum = 0;
1989        total = 0;
1990        for( i=0; i<njob; i++ )
1991        {
1992                myfgets( tmpname, N-1, fp );
1993                tmpseq = load1SeqWithoutName_realloc( fp );
1994                tmp = countnogaplen( tmpseq );
1995                if( tmp > nlenmax ) nlenmax  = tmp;
1996                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
1997                atgcnum += countATGC( tmpseq, &nsite );
1998                total += nsite;
1999                free( tmpseq );
2000        }
2001        free( tmpname );
2002        atgcfreq = (double)atgcnum / total;
2003        fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2004        if( dorp == NOTSPECIFIED )
2005        {
2006                if( atgcfreq > 0.75 )   
2007                {
2008                        dorp = 'd';
2009                        upperCase = -1;
2010                }
2011                else                 
2012                {
2013                        dorp = 'p';
2014                        upperCase = 0;
2015                }
2016        }
2017}
2018
2019
2020void getnumlen_nogap_outallreg( FILE *fp, int *nlenminpt )
2021{
2022        int total;
2023        int nsite = 0;
2024        int atgcnum;
2025        int i, tmp;
2026        char *tmpseq, *tmpname;
2027        double atgcfreq;
2028        tmpname = AllocateCharVec( N );
2029        njob = countKUorWA( fp );
2030        searchKUorWA( fp );
2031        nlenmax = 0;
2032        *nlenminpt = 99999999;
2033        atgcnum = 0;
2034        total = 0;
2035        for( i=0; i<njob; i++ )
2036        {
2037                myfgets( tmpname, N-1, fp );
2038                fprintf( stdout, "%s\n", tmpname );
2039                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2040                tmp = countnogaplen( tmpseq );
2041                fprintf( stdout, "%d\n", tmp );
2042                if( tmp > nlenmax ) nlenmax  = tmp;
2043                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2044                atgcnum += countATGC( tmpseq, &nsite );
2045                total += nsite;
2046                free( tmpseq );
2047        }
2048        free( tmpname );
2049        atgcfreq = (double)atgcnum / total;
2050        fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2051        if( dorp == NOTSPECIFIED )
2052        {
2053                if( atgcfreq > 0.75 )   
2054                {
2055                        dorp = 'd';
2056                        upperCase = -1;
2057                }
2058                else                 
2059                {
2060                        dorp = 'p';
2061                        upperCase = 0;
2062                }
2063        }
2064}
2065
2066static void escapehtml( char *res, char *ori, int maxlen )
2067{
2068        char *res0 = res;
2069        while( *ori )
2070        {
2071                if( *ori == '<' ) 
2072                {
2073                        strcpy( res, "&lt;" );
2074                        res += 3;
2075                }
2076                else if( *ori == '>' ) 
2077                {
2078                        strcpy( res, "&gt;" );
2079                        res += 3;
2080                }
2081                else if( *ori == '&' ) 
2082                {
2083                        strcpy( res, "&amp;" );
2084                        res += 4;
2085                }
2086                else if( *ori == '"' ) 
2087                {
2088                        strcpy( res, "&quot;" );
2089                        res += 5;
2090                }
2091                else if( *ori == ' ' ) 
2092                {
2093                        strcpy( res, "&nbsp;" );
2094                        res += 5;
2095                }
2096                else
2097                {
2098                        *res = *ori;
2099                }
2100                res++;
2101                ori++;
2102
2103                if( res - res0 -10 > N ) break;
2104        }
2105        *res = 0;
2106}
2107
2108void getnumlen_nogap_outallreg_web( FILE *fp, FILE *ofp, int *nlenminpt, int *isalignedpt )
2109{
2110        int total;
2111        int nsite = 0;
2112        int atgcnum;
2113        int alnlen = 0, alnlen_prev;
2114        int i, tmp, lennormalchar;
2115        char *tmpseq, *tmpname, *tmpname2;
2116        double atgcfreq;
2117        tmpname = AllocateCharVec( N );
2118        tmpname2 = AllocateCharVec( N );
2119        njob = countKUorWA( fp );
2120        searchKUorWA( fp );
2121        nlenmax = 0;
2122        *nlenminpt = 99999999;
2123        atgcnum = 0;
2124        total = 0;
2125        alnlen_prev = -1;
2126        *isalignedpt = 1;
2127        for( i=0; i<njob; i++ )
2128        {
2129                myfgets( tmpname, N-1, fp );
2130                tmpname2[0] = tmpname[0];
2131                escapehtml( tmpname2+1, tmpname+1, N );
2132//              fprintf( stdout, "%s\n", tmpname );
2133//              fprintf( stdout, "%s\n", tmpname2 );
2134//              exit(1);
2135                tmpseq = load1SeqWithoutName_realloc_casepreserve( fp );
2136                tmp = countnogaplen( tmpseq );
2137//              fprintf( stdout, "%d\n", tmp );
2138                if( tmp > nlenmax ) nlenmax  = tmp;
2139                if( tmp < *nlenminpt ) *nlenminpt  = tmp;
2140                atgcnum += countATGC( tmpseq, &nsite );
2141                total += nsite;
2142
2143                alnlen = strlen( tmpseq );
2144//              fprintf( stdout, "##### alnlen, alnlen_prev = %d, %d\n", alnlen, alnlen_prev );
2145                if( i>0 && alnlen_prev != alnlen ) *isalignedpt = 0;
2146                alnlen_prev = alnlen;
2147
2148                atgcfreq = (double)atgcnum / total;
2149//              fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2150//              if( dorp == NOTSPECIFIED ) // you kentou
2151                {
2152                        if( atgcfreq > 0.75 )   
2153                        {
2154                                dorp = 'd';
2155                                upperCase = -1;
2156                        }
2157                        else                 
2158                        {
2159                                dorp = 'p';
2160                                upperCase = 0;
2161                        }
2162                }
2163       
2164                if( dorp == 'd' ) lennormalchar = countnormalletters( tmpseq, "atgcuATGCU" );
2165                else              lennormalchar = countnormalletters( tmpseq, "ARNDCQEGHILKMFPSTWYVarndcqeghilkmfpstwyv" );
2166                free( tmpseq );
2167
2168                fprintf( ofp, " <label for='s%d'><span id='ss%d'><input type='checkbox' id='s%d' name='s%d' checked></span> <input type='text' class='ll' id='ll%d' style='display:none' size='6' value='%d' readonly='readonly'>%s</label>\n", i, i, i, i, i, lennormalchar, tmpname2 );
2169                fprintf( ofp, "<span id='t%d-0' style='display:none'>", i );
2170                fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"t%d\")'>+reg</a>", i );
2171                fprintf( ofp, " Begin:<input type='text' name='b%d-0' size='8' value='1' class='ie'> End:<input type='text' name='e%d-0' size='8' value='%d' class='ie'>", i, i, tmp );
2172                if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-0'><input type='checkbox' name='r%d-0' id='r%d-0'>Reverse</label>", i, i, i );
2173//              fprintf( ofp, "  Sequence Length:<input type='text' name='l%d' size='8' value='%d' readonly='readonly'>", i, tmp );
2174                fprintf( ofp, "\n</span>" );
2175                fprintf( ofp, "<span id='t%d-1' style='display:none'>", i );
2176                fprintf( ofp, "      Begin:<input type='text' name='b%d-1' size='8' value='' class='ie'> End:<input type='text' name='e%d-1' size='8' value='' class='ie'>", i, i );
2177                if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-1'><input type='checkbox' name='r%d-1' id='r%d-1'>Reverse</label>", i, i, i );
2178                fprintf( ofp, "\n</span>" );
2179                fprintf( ofp, "<span id='t%d-2' style='display:none'>", i );
2180                fprintf( ofp, "      Begin:<input type='text' name='b%d-2' size='8' value='' class='ie'> End:<input type='text' name='e%d-2' size='8' value='' class='ie'>", i, i );
2181                if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-2'><input type='checkbox' name='r%d-2' id='r%d-2'>Reverse</label>", i, i, i );
2182                fprintf( ofp, "\n</span>" );
2183                fprintf( ofp, "<span id='t%d-3' style='display:none'>", i );
2184                fprintf( ofp, "      Begin:<input type='text' name='b%d-3' size='8' value='' class='ie'> End:<input type='text' name='e%d-3' size='8' value='' class='ie'>", i, i );
2185                if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-3'><input type='checkbox' name='r%d-3' id='r%d-3'>Reverse</label>", i, i, i );
2186                fprintf( ofp, "\n</span>" );
2187                fprintf( ofp, "<span id='t%d-4' style='display:none'>", i );
2188                fprintf( ofp, "      Begin:<input type='text' name='b%d-4' size='8' value='' class='ie'> End:<input type='text' name='e%d-4' size='8' value='' class='ie'>", i, i );
2189                if( dorp == 'd' ) fprintf( ofp, " <label for='r%d-4'><input type='checkbox' name='r%d-4' id='r%d-4'>Reverse</label>", i, i, i );
2190                fprintf( ofp, "\n</span>" );
2191        }
2192        free( tmpname );
2193        free( tmpname2 );
2194        atgcfreq = (double)atgcnum / total;
2195        fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2196//      if( dorp == NOTSPECIFIED ) // you kentou
2197        {
2198                if( atgcfreq > 0.75 )   
2199                {
2200                        dorp = 'd';
2201                        upperCase = -1;
2202                }
2203                else                 
2204                {
2205                        dorp = 'p';
2206                        upperCase = 0;
2207                }
2208        }
2209        fprintf( ofp, "\n" );
2210        if( *isalignedpt )
2211        {
2212                fprintf( ofp, "<span id='tall-0' style='display:none'>" );
2213                fprintf( ofp, "Cut the alignment\n" );
2214                fprintf( ofp, " <a href='javascript:void(0)' onclick='ddcycle(this.form,\"tall\")'>+reg</a>" );
2215                fprintf( ofp, " Begin:<input type='text' name='ball-0' size='8' value='1'> End:<input type='text' name='eall-0' size='8' value='%d'>", alnlen );
2216                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-0'><input type='checkbox' name='rall-0' id='rall-0'>Reverse</label>" );
2217                fprintf( ofp, "  Alignment length:<input type='text' name='lall' size='8' value='%d' readonly='readonly'>", alnlen );
2218                fprintf( ofp, "\n</span>" );
2219                fprintf( ofp, "<span id='tall-1' style='display:none'>" );
2220                fprintf( ofp, "      Begin:<input type='text' name='ball-1' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2221                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-1'><input type='checkbox' name='rall-1' id='rall-1'>Reverse</label>" );
2222                fprintf( ofp, "\n</span>" );
2223                fprintf( ofp, "<span id='tall-2' style='display:none'>" );
2224                fprintf( ofp, "      Begin:<input type='text' name='ball-2' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2225                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-2'><input type='checkbox' name='rall-2' id='rall-2'>Reverse</label>" );
2226                fprintf( ofp, "\n</span>" );
2227                fprintf( ofp, "<span id='tall-3' style='display:none'>" );
2228                fprintf( ofp, "      Begin:<input type='text' name='ball-3' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2229                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-3'><input type='checkbox' name='rall-3' id='rall-3'>Reverse</label>" );
2230                fprintf( ofp, "\n</span>" );
2231                fprintf( ofp, "<span id='tall-4' style='display:none'>" );
2232                fprintf( ofp, "      Begin:<input type='text' name='ball-4' size='8' value=''> End:<input type='text' name='eall-1' size='8' value=''>" );
2233                if( dorp == 'd' ) fprintf( ofp, " <label for='rall-4'><input type='checkbox' name='rall-4' id='rall-4'>Reverse</label>" );
2234                fprintf( ofp, "\n</span>" );
2235        }
2236
2237}
2238
2239void getnumlen( FILE *fp )
2240{
2241        int total;
2242        int nsite = 0;
2243        int atgcnum;
2244        int i, tmp;
2245        char *tmpseq;
2246        char *tmpname;
2247        double atgcfreq;
2248        tmpname = AllocateCharVec( N );
2249        njob = countKUorWA( fp );
2250        searchKUorWA( fp );
2251        nlenmax = 0;
2252        atgcnum = 0;
2253        total = 0;
2254        for( i=0; i<njob; i++ )
2255        {
2256                myfgets( tmpname, N-1, fp );
2257                tmpseq = load1SeqWithoutName_realloc( fp );
2258                tmp = strlen( tmpseq );
2259                if( tmp > nlenmax ) nlenmax  = tmp;
2260                atgcnum += countATGC( tmpseq, &nsite );
2261                total += nsite;
2262//              fprintf( stderr, "##### total = %d\n", total );
2263                free( tmpseq );
2264        }
2265
2266
2267        atgcfreq = (double)atgcnum / total;
2268//      fprintf( stderr, "##### atgcfreq = %f\n", atgcfreq );
2269        if( dorp == NOTSPECIFIED )
2270        {
2271                if( atgcfreq > 0.75 )   
2272                {
2273                        dorp = 'd';
2274                        upperCase = -1;
2275                }
2276                else                 
2277                {
2278                        dorp = 'p';
2279                        upperCase = 0;
2280                }
2281        }
2282        free( tmpname );
2283}
2284       
2285
2286
2287void WriteGapFill( FILE *fp, int locnjob, char name[][B], int nlen[M], char **aseq )
2288{
2289        static char b[N];
2290        int i, j;
2291        int nalen[M];
2292        static char gap[N];
2293        static char buff[N];
2294
2295#if IODEBUG
2296        fprintf( stderr, "IMAKARA KAKU\n" );
2297#endif
2298        nlenmax = 0;
2299        for( i=0; i<locnjob; i++ )
2300        {
2301                int len = strlen( aseq[i] );
2302                if( nlenmax < len ) nlenmax = len;
2303        }
2304
2305        for( i=0; i<nlenmax; i++ ) gap[i] = '-';
2306        gap[nlenmax] = 0;
2307
2308        fprintf( fp, "%5d", locnjob );
2309        fprintf( fp, "\n" );
2310
2311        for( i=0; i<locnjob; i++ )
2312        {
2313                strcpy( buff, aseq[i] );
2314                strncat( buff, gap, nlenmax-strlen( aseq[i] ) );
2315                buff[nlenmax] = 0;
2316                nalen[i] = strlen( buff );
2317                fprintf( fp, "%s\n", name[i] );
2318                fprintf( fp, "%5d\n", nalen[i] );
2319                for( j=0; j<nalen[i]; j=j+C )
2320                {
2321                        strncpy_caseC( b, buff+j, C ); b[C] = 0;
2322                        fprintf( fp, "%s\n",b );
2323                }
2324        }
2325#if DEBUG
2326        fprintf( stderr, "nalen[0] = %d\n", nalen[0] );
2327#endif
2328#if IODEBUG
2329        fprintf( stderr, "KAKIOWATTA\n" );
2330#endif
2331}
2332
2333void writeDataforgaln( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2334{
2335        int i, j;
2336        int nalen;
2337
2338        for( i=0; i<locnjob; i++ )
2339        {
2340                nalen = strlen( aseq[i] );
2341                fprintf( fp, ">%s\n", name[i]+1 );
2342                for( j=0; j<nalen; j=j+C )
2343                {
2344#if 0
2345                        strncpy( b, aseq[i]+j, C ); b[C] = 0;
2346                        fprintf( fp, "%s\n",b );
2347#else
2348                        fprintf( fp, "%.*s\n", C, aseq[i]+j );
2349#endif
2350                }
2351        }
2352}
2353
2354void writeData_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq )
2355{
2356        int i, j;
2357        int nalen;
2358
2359        for( i=0; i<locnjob; i++ )
2360        {
2361#if DEBUG
2362                fprintf( stderr, "i = %d in writeData\n", i );
2363#endif
2364                nalen = strlen( aseq[i] );
2365                fprintf( fp, ">%s\n", name[i]+1 );
2366                for( j=0; j<nalen; j=j+C )
2367                {
2368#if 0
2369                        strncpy( b, aseq[i]+j, C ); b[C] = 0;
2370                        fprintf( fp, "%s\n",b );
2371#else
2372                        fprintf( fp, "%.*s\n", C, aseq[i]+j );
2373#endif
2374                }
2375        }
2376}
2377
2378void writeData( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq )
2379{
2380        int i, j;
2381        int nalen;
2382
2383        for( i=0; i<locnjob; i++ )
2384        {
2385#if DEBUG
2386                fprintf( stderr, "i = %d in writeData\n", i );
2387#endif
2388                nalen = strlen( aseq[i] );
2389                fprintf( fp, ">%s\n", name[i]+1 );
2390                for( j=0; j<nalen; j=j+C )
2391                {
2392#if 0
2393                        strncpy( b, aseq[i]+j, C ); b[C] = 0;
2394                        fprintf( fp, "%s\n",b );
2395#else
2396                        fprintf( fp, "%.*s\n", C, aseq[i]+j );
2397#endif
2398                }
2399        }
2400}
2401
2402
2403void write1seq( FILE *fp, char *aseq )
2404{
2405        int j;
2406        int nalen;
2407
2408        nalen = strlen( aseq );
2409        for( j=0; j<nalen; j=j+C )
2410                fprintf( fp, "%.*s\n", C, aseq+j );
2411}
2412
2413void readhat2_floathalf_part_pointer( FILE *fp, int nseq, int nadd, char **name, float **mtx )
2414{
2415    int i, j, nseq0, norg;
2416    char b[B];
2417
2418    fgets_IGNORE_RESULT( b, B, fp );
2419    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) 
2420        {
2421                fprintf( stderr, "%d != %d\n", nseq, nseq0 );
2422                ErrorExit( "hat2 is wrong." );
2423        }
2424    fgets_IGNORE_RESULT( b, B, fp );
2425    for( i=0; i<nseq; i++ )
2426    {
2427#if 0
2428        getaline_fp_eof( b, B, fp );
2429#else
2430                myfgets( b, B-2, fp );
2431#endif
2432#if 0
2433                j = MIN( strlen( b+6 ), 10 );
2434        if( strncmp( name[i], b+6 , j ) )
2435                {
2436                        fprintf( stderr, "Error in hat2\n" );
2437                        fprintf( stderr, "%s != %s\n", b, name[i] );
2438                        exit( 1 );
2439                }
2440#endif
2441    }
2442        norg = nseq-nadd;
2443    for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ )
2444    {
2445        mtx[i][j] = ( input_new( fp, D ) );
2446    }
2447}
2448
2449void readhat2_floathalf_pointer( FILE *fp, int nseq, char **name, float **mtx )
2450{
2451    int i, j, nseq0;
2452    char b[B];
2453
2454    fgets_IGNORE_RESULT( b, B, fp );
2455    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) 
2456        {
2457                fprintf( stderr, "%d != %d\n", nseq, nseq0 );
2458                ErrorExit( "hat2 is wrong." );
2459        }
2460    fgets_IGNORE_RESULT( b, B, fp );
2461    for( i=0; i<nseq; i++ )
2462    {
2463#if 0
2464        getaline_fp_eof( b, B, fp );
2465#else
2466                myfgets( b, B-2, fp );
2467#endif
2468#if 0
2469                j = MIN( strlen( b+6 ), 10 );
2470        if( strncmp( name[i], b+6 , j ) )
2471                {
2472                        fprintf( stderr, "Error in hat2\n" );
2473                        fprintf( stderr, "%s != %s\n", b, name[i] );
2474                        exit( 1 );
2475                }
2476#endif
2477    }
2478    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2479    {
2480        mtx[i][j-i] = ( input_new( fp, D ) );
2481    }
2482}
2483void readhat2_floathalf( FILE *fp, int nseq, char name[M][B], float **mtx )
2484{
2485    int i, j, nseq0;
2486    char b[B];
2487
2488    fgets_IGNORE_RESULT( b, B, fp );
2489    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2490    fgets_IGNORE_RESULT( b, B, fp );
2491    for( i=0; i<nseq; i++ )
2492    {
2493#if 0
2494        getaline_fp_eof( b, B, fp );
2495#else
2496                myfgets( b, B-2, fp );
2497#endif
2498#if 0
2499                j = MIN( strlen( b+6 ), 10 );
2500        if( strncmp( name[i], b+6 , j ) )
2501                {
2502                        fprintf( stderr, "Error in hat2\n" );
2503                        fprintf( stderr, "%s != %s\n", b, name[i] );
2504                        exit( 1 );
2505                }
2506#endif
2507    }
2508    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2509    {
2510        mtx[i][j-i] = ( input_new( fp, D ) );
2511    }
2512}
2513void readhat2_float( FILE *fp, int nseq, char name[M][B], float **mtx )
2514{
2515    int i, j, nseq0;
2516    char b[B];
2517
2518    fgets_IGNORE_RESULT( b, B, fp );
2519    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2520    fgets_IGNORE_RESULT( b, B, fp );
2521    for( i=0; i<nseq; i++ )
2522    {
2523#if 0
2524        getaline_fp_eof( b, B, fp );
2525#else
2526                myfgets( b, B-2, fp );
2527#endif
2528#if 0
2529                j = MIN( strlen( b+6 ), 10 );
2530        if( strncmp( name[i], b+6 , j ) )
2531                {
2532                        fprintf( stderr, "Error in hat2\n" );
2533                        fprintf( stderr, "%s != %s\n", b, name[i] );
2534                        exit( 1 );
2535                }
2536#endif
2537    }
2538    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2539    {
2540        mtx[i][j] = ( input_new( fp, D ) );
2541    }
2542}
2543void readhat2_int( FILE *fp, int nseq, char name[M][B], int **mtx )
2544{
2545    int i, j, nseq0;
2546    char b[B];
2547
2548    fgets_IGNORE_RESULT( b, B, fp );
2549    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2550    fgets_IGNORE_RESULT( b, B, fp );
2551    for( i=0; i<nseq; i++ )
2552    {
2553#if 0
2554        getaline_fp_eof( b, B, fp );
2555#else
2556                myfgets( b, B-2, fp );
2557#endif
2558#if 0
2559                j = MIN( strlen( b+6 ), 10 );
2560        if( strncmp( name[i], b+6 , j ) )
2561                {
2562                        fprintf( stderr, "Error in hat2\n" );
2563                        fprintf( stderr, "%s != %s\n", b, name[i] );
2564                        exit( 1 );
2565                }
2566#endif
2567    }
2568    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2569    {
2570        mtx[i][j] = (int)( input_new( fp, D ) * INTMTXSCALE + 0.5 );
2571    }
2572}
2573
2574void readhat2_pointer( FILE *fp, int nseq, char **name, double **mtx )
2575{
2576    int i, j, nseq0;
2577    char b[B];
2578
2579    fgets_IGNORE_RESULT( b, B, fp );
2580    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2581    fgets_IGNORE_RESULT( b, B, fp );
2582    for( i=0; i<nseq; i++ )
2583    {
2584#if 0
2585        getaline_fp_eof( b, B, fp );
2586#else
2587                myfgets( b, B-2, fp );
2588#endif
2589#if 0
2590                j = MIN( strlen( b+6 ), 10 );
2591        if( strncmp( name[i], b+6 , j ) )
2592                {
2593                        fprintf( stderr, "Error in hat2\n" );
2594                        fprintf( stderr, "%s != %s\n", b, name[i] );
2595                        exit( 1 );
2596                }
2597#endif
2598    }
2599    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2600    {
2601        mtx[i][j] = (double)input_new( fp, D);
2602    }
2603}
2604void readhat2( FILE *fp, int nseq, char name[M][B], double **mtx )
2605{
2606    int i, j, nseq0;
2607    char b[B];
2608
2609    fgets_IGNORE_RESULT( b, B, fp );
2610    fgets_IGNORE_RESULT( b, B, fp ); b[5] = 0; nseq0 = atoi( b ); if( nseq != nseq0 ) ErrorExit( "hat2 is wrong." );
2611    fgets_IGNORE_RESULT( b, B, fp );
2612    for( i=0; i<nseq; i++ )
2613    {
2614#if 0
2615        getaline_fp_eof( b, B, fp );
2616#else
2617                myfgets( b, B-2, fp );
2618#endif
2619#if 0
2620                j = MIN( strlen( b+6 ), 10 );
2621        if( strncmp( name[i], b+6 , j ) )
2622                {
2623                        fprintf( stderr, "Error in hat2\n" );
2624                        fprintf( stderr, "%s != %s\n", b, name[i] );
2625                        exit( 1 );
2626                }
2627#endif
2628    }
2629    for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ )
2630    {
2631        mtx[i][j] = (double)input_new( fp, D);
2632    }
2633}
2634
2635void WriteFloatHat2_pointer_halfmtx( FILE *hat2p, int locnjob, char **name, float **mtx )
2636{
2637        int i, j, ijsa;
2638        double max = 0.0;
2639        for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2640
2641        fprintf( hat2p, "%5d\n", 1 );
2642        fprintf( hat2p, "%5d\n", locnjob );
2643        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2644
2645        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2646        for( i=0; i<locnjob; i++ )
2647        {
2648                for( j=i+1; j<njob; j++ )
2649                {
2650                        fprintf( hat2p, "%#6.3f", mtx[i][j-i] );
2651                        ijsa = j-i;
2652                        if( ijsa % 12 == 0 || ijsa == locnjob-i-1 ) fprintf( hat2p, "\n" );
2653                }
2654        }
2655}
2656
2657void WriteFloatHat2_pointer( FILE *hat2p, int locnjob, char **name, float **mtx )
2658{
2659        int i, j;
2660        double max = 0.0;
2661        for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2662
2663        fprintf( hat2p, "%5d\n", 1 );
2664        fprintf( hat2p, "%5d\n", locnjob );
2665        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2666
2667        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2668        for( i=0; i<locnjob; i++ )
2669        {
2670                for( j=1; j<locnjob-i; j++ ) 
2671                {
2672                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2673                        if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2674                }
2675        }
2676}
2677
2678void WriteFloatHat2( FILE *hat2p, int locnjob, char name[M][B], float **mtx )
2679{
2680        int i, j;
2681        double max = 0.0;
2682        for( i=0; i<locnjob-1; i++ ) for( j=1; j<locnjob-i; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2683
2684        fprintf( hat2p, "%5d\n", 1 );
2685        fprintf( hat2p, "%5d\n", locnjob );
2686        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2687
2688        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2689        for( i=0; i<locnjob; i++ )
2690        {
2691                for( j=1; j<locnjob-i; j++ ) 
2692                {
2693                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2694                        if( j % 12 == 0 || j == locnjob-i-1 ) fprintf( hat2p, "\n" );
2695                }
2696        }
2697}
2698
2699void WriteHat2_int( FILE *hat2p, int locnjob, char name[M][B], int **mtx )
2700{
2701        int i, j;
2702        double max = 0.0;
2703        for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2704        max /= INTMTXSCALE;
2705
2706        fprintf( hat2p, "%5d\n", 1 );
2707        fprintf( hat2p, "%5d\n", locnjob );
2708        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2709
2710        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2711        for( i=0; i<locnjob-1; i++ )
2712        {
2713                for( j=i+1; j<locnjob; j++ ) 
2714                {
2715                        fprintf( hat2p, "%#6.3f", (float)mtx[i][j] / INTMTXSCALE );
2716                        if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2717                }
2718        }
2719}
2720
2721void WriteHat2_part_pointer( FILE *hat2p, int locnjob, int nadd, char **name, double **mtx )
2722{
2723        int i, j;
2724        int norg = locnjob-nadd;
2725        double max = 0.0;
2726//      for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2727
2728        fprintf( hat2p, "%5d\n", 1 );
2729        fprintf( hat2p, "%5d\n", locnjob );
2730        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2731
2732        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2733        for( i=0; i<norg; i++ )
2734        {
2735                for( j=0; j<nadd; j++ ) 
2736                {
2737                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2738                        if( (j+1) % 12 == 0 || j == nadd-1 ) fprintf( hat2p, "\n" );
2739                }
2740        }
2741}
2742
2743void WriteHat2_pointer( FILE *hat2p, int locnjob, char **name, double **mtx )
2744{
2745        int i, j;
2746        double max = 0.0;
2747        for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2748
2749        fprintf( hat2p, "%5d\n", 1 );
2750        fprintf( hat2p, "%5d\n", locnjob );
2751        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2752
2753        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2754        for( i=0; i<locnjob-1; i++ )
2755        {
2756                for( j=i+1; j<locnjob; j++ ) 
2757                {
2758                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2759                        if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2760                }
2761        }
2762}
2763
2764void WriteHat2( FILE *hat2p, int locnjob, char name[M][B], double **mtx )
2765{
2766        int i, j;
2767        double max = 0.0;
2768        for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ ) if( mtx[i][j] > max ) max = mtx[i][j];
2769
2770        fprintf( hat2p, "%5d\n", 1 );
2771        fprintf( hat2p, "%5d\n", locnjob );
2772        fprintf( hat2p, " %#6.3f\n", max * 2.5 );
2773
2774        for( i=0; i<locnjob; i++ ) fprintf( hat2p, "%4d. %s\n", i+1, name[i] );
2775        for( i=0; i<locnjob-1; i++ )
2776        {
2777                for( j=i+1; j<locnjob; j++ ) 
2778                {
2779                        fprintf( hat2p, "%#6.3f", mtx[i][j] );
2780                        if( (j-i) % 12 == 0 || j == locnjob-1 ) fprintf( hat2p, "\n" );
2781                }
2782        }
2783}
2784
2785#if 0
2786void WriteHat2plain( FILE *hat2p, int locnjob, double **mtx )
2787{
2788        int i, j, ilim;
2789
2790        ilim = locnjob-1;
2791        for( i=0; i<ilim; i++ )
2792        {
2793                fprintf( hat2p, "%d-%d d=%.3f\n", i+1, i+1, 0.0 );
2794                for( j=i+1; j<locnjob; j++ )
2795                {
2796                        fprintf( hat2p, "%d-%d d=%.3f\n", i+1, j+1, mtx[i][j] );
2797                }
2798        }
2799}
2800#endif
2801
2802int ReadFasta_sub( FILE *fp, double *dis, int nseq, char name[M][B] )
2803{
2804    int i, count=0;
2805    char b[B];
2806    int junban[M];
2807
2808    count = 0;
2809    for( i=0; i<10000000 && count<nseq; i++ )
2810    {
2811        fgets_IGNORE_RESULT( b, B-1, fp );
2812        if( !strncmp( "+==========+", b, 12 ) )
2813        {
2814            junban[count] = atoi( b+12 );
2815            count++;
2816        }
2817    }
2818
2819        for( i=0; i<nseq; i++ ) dis[i] = 0.0;
2820    count = 0;
2821    for( i=0; i<100000 && count<nseq; i++ )
2822    {
2823                if( fgets( b, B-1, fp ) ) break;
2824        if( !strncmp( name[junban[count]], b, 20  ) )
2825        {
2826            fgets_IGNORE_RESULT( b, B-1, fp );
2827            dis[junban[count]] = atof( b );
2828            count++;
2829        }
2830    }
2831    return 0;
2832}
2833
2834
2835int ReadSsearch( FILE *fp, double *dis, int nseq, char name[M][B] )
2836{
2837    int i, count=0;
2838    char b[B];
2839    int junban[M];
2840        int opt;
2841
2842    count = 0;
2843    for( i=0; i<10000000 && count<nseq; i++ )
2844    {
2845        fgets_IGNORE_RESULT( b, B-1, fp );
2846        if( !strncmp( "+==========+", b, 12 ) )
2847        {
2848            junban[count] = atoi( b+12 );
2849                        sscanf( b+75, "%d", &opt ); 
2850            dis[junban[count]] = (double)opt;
2851            count++;
2852        }
2853    }
2854
2855/*
2856    count = 0;
2857    for( i=0; i<100000 && count<nseq; i++ )
2858    {
2859        fgets( b, B-1, fp );
2860        if( !strncmp( name[junban[count]], b, 20  ) )
2861        {
2862            dis[junban[count]] = atof( b+65 );
2863            count++;
2864        }
2865    }
2866*/
2867    return 0;
2868}
2869
2870int ReadBlastm7_avscore( FILE *fp, double *dis, int nin )
2871{
2872    int count=0;
2873    char b[B];
2874        char *pt;
2875    int *junban;
2876        double score, sumscore;
2877        double len, sumlen;
2878        int qstart, qend, tstart, tend;
2879        double scorepersite;
2880        static char qal[N], tal[N], al[N];
2881        int nlocalhom;
2882
2883        junban = calloc( nin, sizeof( int ) );
2884
2885        count = 0;
2886        sumscore = 0.0;
2887        sumlen = 0.0;
2888        score = 0.0;
2889        len = 0.0;
2890        scorepersite = 0.0; // by D.Mathog, a guess
2891    while( 1 )
2892        {
2893
2894                if( feof( fp ) ) break;
2895
2896                while( fgets( b, B-1, fp ) )
2897                {
2898                        if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
2899                }
2900
2901                if( !strncmp( "          <Hit_def>", b, 19 ) )
2902                {
2903                        junban[count] = atoi( b+31 );
2904                        nlocalhom = 0;
2905                }
2906
2907
2908                while( fgets( b, B-1, fp ) )
2909                        if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
2910                pt = b + 25;
2911                score = atof( pt );
2912                sumscore += score;
2913
2914
2915                while( fgets( b, B-1, fp ) )
2916                        if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
2917                pt = b + 30;
2918                qstart = atoi( pt ) - 1;
2919
2920
2921                while( fgets( b, B-1, fp ) )
2922                        if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
2923                pt = b + 28;
2924                qend = atoi( pt ) - 1;
2925
2926
2927                while( fgets( b, B-1, fp ) )
2928                        if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
2929                pt = b + 28;
2930                tstart = atoi( pt ) - 1;
2931
2932
2933                while( fgets( b, B-1, fp ) )
2934                        if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
2935                pt = b + 26;
2936                tend = atoi( pt ) - 1;
2937
2938
2939                while( fgets( b, B-1, fp ) )
2940                        if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
2941                pt = b + 29;
2942                len = atoi( pt );
2943                sumlen += len;
2944
2945
2946                while( fgets( al, N-100, fp ) )
2947                        if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
2948
2949                strcpy( qal, al+24 );
2950                pt = qal;
2951                while( *++pt != '<' )
2952                        ;
2953                *pt = 0;
2954
2955
2956                while( fgets( al, N-100, fp ) )
2957                        if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
2958
2959                strcpy( tal, al+24 );
2960                pt = tal;
2961                while( *++pt != '<' )
2962                        ;
2963                *pt = 0;
2964
2965
2966//              fprintf( stderr, "t=%d, score = %f, qstart=%d, qend=%d, tstart=%d, tend=%d, overlapaa=%d\n", junban[count], score, qstart, qend, tstart, tend, overlapaa );
2967
2968
2969                while( fgets( b, B-1, fp ) )
2970                        if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
2971
2972
2973                fgets_IGNORE_RESULT( b, B-1, fp );
2974
2975
2976                if( !strncmp( "          </Hit_hsps>", b, 21 ) )
2977                {
2978                        dis[junban[count++]] = sumscore;
2979                        sumscore = 0.0;
2980                        fgets_IGNORE_RESULT( b, B-1, fp );
2981                        fgets_IGNORE_RESULT( b, B-1, fp );
2982                        scorepersite = sumscore / sumlen;
2983                        if( scorepersite != (int)scorepersite )
2984                        {
2985                                fprintf( stderr, "ERROR! sumscore=%f, sumlen=%f, and scorepersite=%f\n", sumscore, sumlen, scorepersite );
2986                                exit( 1 );
2987                        }
2988
2989                        if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
2990                }
2991        }
2992
2993        free( junban );
2994
2995    return (int)scorepersite;
2996}
2997int ReadBlastm7_scoreonly( FILE *fp, double *dis, int nin )
2998{
2999    int count=0;
3000    char b[B];
3001        char *pt;
3002    int *junban;
3003        int overlapaa;
3004        double score, sumscore;
3005        int qstart, qend, tstart, tend;
3006        static char qal[N], tal[N], al[N];
3007        int nlocalhom;
3008
3009        junban = calloc( nin, sizeof( int ) );
3010
3011        count = 0;
3012        sumscore = 0.0;
3013        score = 0.0;
3014    while( 1 )
3015        {
3016
3017                if( feof( fp ) ) break;
3018
3019                while( fgets( b, B-1, fp ) )
3020                {
3021                        if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
3022                }
3023
3024                if( !strncmp( "          <Hit_def>", b, 19 ) )
3025                {
3026                        junban[count] = atoi( b+31 );
3027                        nlocalhom = 0;
3028                }
3029
3030
3031                while( fgets( b, B-1, fp ) )
3032                        if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
3033                pt = b + 25;
3034                score = atof( pt );
3035                sumscore += score;
3036
3037
3038                while( fgets( b, B-1, fp ) )
3039                        if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
3040                pt = b + 30;
3041                qstart = atoi( pt ) - 1;
3042
3043
3044                while( fgets( b, B-1, fp ) )
3045                        if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
3046                pt = b + 28;
3047                qend = atoi( pt ) - 1;
3048
3049
3050                while( fgets( b, B-1, fp ) )
3051                        if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3052                pt = b + 28;
3053                tstart = atoi( pt ) - 1;
3054
3055
3056                while( fgets( b, B-1, fp ) )
3057                        if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3058                pt = b + 26;
3059                tend = atoi( pt ) - 1;
3060
3061
3062                while( fgets( b, B-1, fp ) )
3063                        if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3064                pt = b + 29;
3065                overlapaa = atoi( pt );
3066
3067
3068                while( fgets( al, N-100, fp ) )
3069                        if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3070
3071                strcpy( qal, al+24 );
3072                pt = qal;
3073                while( *++pt != '<' )
3074                        ;
3075                *pt = 0;
3076
3077
3078                while( fgets( al, N-100, fp ) )
3079                        if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3080
3081                strcpy( tal, al+24 );
3082                pt = tal;
3083                while( *++pt != '<' )
3084                        ;
3085                *pt = 0;
3086
3087
3088//              fprintf( stderr, "t=%d, score = %f, qstart=%d, qend=%d, tstart=%d, tend=%d, overlapaa=%d\n", junban[count], score, qstart, qend, tstart, tend, overlapaa );
3089
3090//              nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
3091
3092                while( fgets( b, B-1, fp ) )
3093                        if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3094
3095
3096                fgets_IGNORE_RESULT( b, B-1, fp );
3097
3098
3099                if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3100                {
3101                        dis[junban[count++]] = sumscore;
3102                        sumscore = 0.0;
3103                        fgets_IGNORE_RESULT( b, B-1, fp );
3104                        fgets_IGNORE_RESULT( b, B-1, fp );
3105                        if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3106                }
3107        }
3108
3109        free( junban );
3110
3111    return count;
3112}
3113
3114int ReadBlastm7( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3115{
3116    int count=0;
3117    char b[B];
3118        char *pt;
3119    static int junban[M];
3120        int overlapaa;
3121        double score, sumscore;
3122        int qstart, qend, tstart, tend;
3123        static char qal[N], tal[N], al[N];
3124        int nlocalhom;
3125
3126
3127
3128        count = 0;
3129        sumscore = 0.0;
3130        score = 0.0;
3131        nlocalhom = 0;
3132    while( 1 )
3133        {
3134
3135                if( feof( fp ) ) break;
3136
3137                while( fgets( b, B-1, fp ) )
3138                {
3139                        if( !strncmp( "          <Hit_def>", b, 19 ) || !strncmp( "              <Hsp_num>", b, 23 ) ) break;
3140                }
3141
3142                if( !strncmp( "          <Hit_def>", b, 19 ) )
3143                {
3144                        junban[count] = atoi( b+31 );
3145                        nlocalhom = 0;
3146                }
3147
3148
3149                while( fgets( b, B-1, fp ) )
3150                        if( !strncmp( "              <Hsp_score>", b, 25 ) ) break;
3151                pt = b + 25;
3152                score = atof( pt );
3153                sumscore += score;
3154
3155
3156                while( fgets( b, B-1, fp ) )
3157                        if( !strncmp( "              <Hsp_query-from>", b, 30 ) ) break;
3158                pt = b + 30;
3159                qstart = atoi( pt ) - 1;
3160
3161
3162                while( fgets( b, B-1, fp ) )
3163                        if( !strncmp( "              <Hsp_query-to>", b, 28 ) ) break;
3164                pt = b + 28;
3165                qend = atoi( pt ) - 1;
3166
3167
3168                while( fgets( b, B-1, fp ) )
3169                        if( !strncmp( "              <Hsp_hit-from>", b, 28 ) ) break;
3170                pt = b + 28;
3171                tstart = atoi( pt ) - 1;
3172
3173
3174                while( fgets( b, B-1, fp ) )
3175                        if( !strncmp( "              <Hsp_hit-to>", b, 26 ) ) break;
3176                pt = b + 26;
3177                tend = atoi( pt ) - 1;
3178
3179
3180                while( fgets( b, B-1, fp ) )
3181                        if( !strncmp( "              <Hsp_align-len>", b, 29 ) ) break;
3182                pt = b + 29;
3183                overlapaa = atoi( pt );
3184
3185
3186                while( fgets( al, N-100, fp ) )
3187                        if( !strncmp( "              <Hsp_qseq>", al, 24 ) ) break;
3188
3189                strcpy( qal, al+24 );
3190                pt = qal;
3191                while( *++pt != '<' )
3192                        ;
3193                *pt = 0;
3194
3195
3196                while( fgets( al, N-100, fp ) )
3197                        if( !strncmp( "              <Hsp_hseq>", al, 24 ) ) break;
3198
3199                strcpy( tal, al+24 );
3200                pt = tal;
3201                while( *++pt != '<' )
3202                        ;
3203                *pt = 0;
3204
3205
3206//              fprintf( stderr, "t=%d, score = %f, qstart=%d, qend=%d, tstart=%d, tend=%d, overlapaa=%d\n", junban[count], score, qstart, qend, tstart, tend, overlapaa );
3207
3208                nlocalhom += addlocalhom_r( qal, tal, localhomlist+junban[count], qstart, tstart, score, overlapaa, nlocalhom );
3209
3210                while( fgets( b, B-1, fp ) )
3211                        if( !strncmp( "            </Hsp>:", b, 18 ) ) break;
3212
3213
3214                fgets_IGNORE_RESULT( b, B-1, fp );
3215
3216
3217                if( !strncmp( "          </Hit_hsps>", b, 21 ) )
3218                {
3219                        dis[junban[count++]] = sumscore;
3220                        sumscore = 0.0;
3221                        fgets_IGNORE_RESULT( b, B-1, fp );
3222                        fgets_IGNORE_RESULT( b, B-1, fp );
3223                        if( !strncmp( "      </Iteration_hits>", b, 23 ) ) break;
3224                }
3225        }
3226    return count;
3227}
3228
3229int ReadFasta34noalign( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3230{
3231    int count=0;
3232    char b[B];
3233        char *pt;
3234    static int junban[M];
3235        int opt;
3236        double z, bits;
3237
3238
3239    count = 0;
3240#if 0
3241    for( i=0; i<10000000 && count<nseq; i++ )
3242#else
3243    while( !feof( fp ) )
3244#endif
3245    {
3246        fgets_IGNORE_RESULT( b, B-1, fp );
3247        if( !strncmp( "+==========+", b, 12 ) )
3248        {
3249            junban[count] = atoi( b+12 );
3250
3251                        pt = strchr( b, ')' ) + 1;
3252                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3253            dis[junban[count]] = (double)opt;
3254            count++;
3255
3256        }
3257    }
3258
3259    return count;
3260}
3261int ReadFasta34m10_nuc( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3262{
3263    int count=0;
3264    char b[B];
3265        char *pt;
3266    static int junban[M];
3267        int overlapaa;
3268        int opt, qstart, qend, tstart, tend;
3269        double z, bits;
3270        int qal_display_start, tal_display_start;
3271        static char qal[N], tal[N];
3272        char *qal2, *tal2;
3273        int c;
3274
3275
3276    count = 0;
3277#if 0
3278    for( i=0; i<10000000 && count<nseq; i++ )
3279#else
3280    while( !feof( fp ) )
3281#endif
3282    {
3283        fgets_IGNORE_RESULT( b, B-1, fp );
3284        if( !strncmp( "+==========+", b, 12 ) )
3285        {
3286            junban[count] = atoi( b+12 );
3287
3288                        if( strchr( b, 'r' ) ) continue;
3289
3290                        pt = strchr( b, ']' ) + 1;
3291                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3292            dis[junban[count]] = (double)opt;
3293            count++;
3294
3295        }
3296                else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3297                {
3298                        break;
3299                }
3300
3301    }
3302        if( !count ) return -1;
3303
3304        count = 0;
3305    while( 1 )
3306        {
3307                if( strncmp( ">>+==========+", b, 14 ) )
3308                {
3309                        fgets_IGNORE_RESULT( b, B-1, fp );
3310                        if( feof( fp ) ) break;
3311                        continue;
3312                }
3313                junban[count++] = atoi( b+14 );
3314//              fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3315                while( fgets( b, B-1, fp ) )
3316                        if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3317                pt = strstr( b, ":" ) +1;
3318                opt = atoi( pt );
3319
3320
3321                while( fgets( b, B-1, fp ) )
3322                        if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3323                pt = strstr( b, ":" ) +1;
3324                overlapaa = atoi( pt );
3325
3326                while( fgets( b, B-1, fp ) )
3327                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3328                pt = strstr( b, ":" ) +1;
3329                qstart = atoi( pt ) - 1;
3330
3331                while( fgets( b, B-1, fp ) )
3332                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3333                pt = strstr( b, ":" ) +1;
3334                qend = atoi( pt ) - 1;
3335
3336                while( fgets( b, B-1, fp ) )
3337                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3338                pt = strstr( b, ":" ) +1;
3339                qal_display_start = atoi( pt ) - 1;
3340
3341                pt = qal;
3342                while( (c = fgetc( fp )) )
3343                {
3344                        if( c == '>' ) 
3345                        {
3346                                ungetc( c, fp );
3347                                break;
3348                        }
3349                        if( isalpha( c ) || c == '-' ) 
3350                        *pt++ = c;
3351                }
3352                *pt = 0;
3353
3354                while( fgets( b, B-1, fp ) )
3355                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3356                pt = strstr( b, ":" ) + 1;
3357                tstart = atoi( pt ) - 1;
3358
3359                while( fgets( b, B-1, fp ) )
3360                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3361                pt = strstr( b, ":" ) + 1;
3362                tend = atoi( pt ) - 1;
3363
3364                while( fgets( b, B-1, fp ) )
3365                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3366                pt = strstr( b, ":" ) + 1;
3367                tal_display_start = atoi( pt ) - 1;
3368
3369                pt = tal;
3370                while( ( c = fgetc( fp ) ) )
3371                {
3372                        if( c == '>' ) 
3373                        {
3374                                ungetc( c, fp );
3375                                break;
3376                        }
3377                        if( isalpha( c ) || c == '-' ) 
3378                        *pt++ = c;
3379                }
3380                *pt = 0;
3381
3382//              fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3383//              fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3384
3385//              fprintf( stderr, "qal = %s\n", qal );
3386//              fprintf( stderr, "tal = %s\n", tal );
3387
3388                qal2 = cutal( qal, qal_display_start, qstart, qend );
3389                tal2 = cutal( tal, tal_display_start, tstart, tend );
3390
3391//              fprintf( stderr, "qal2 = %s\n", qal2 );
3392//              fprintf( stderr, "tal2 = %s\n", tal2 );
3393
3394//              fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3395                putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3396        }
3397//      fprintf( stderr, "count = %d\n", count );
3398    return count;
3399}
3400int ReadFasta34m10( FILE *fp, double *dis, int qmem, char **name, LocalHom *localhomlist )
3401{
3402    int count=0;
3403    char b[B];
3404        char *pt;
3405    static int junban[M];
3406        int overlapaa;
3407        int opt, qstart, qend, tstart, tend;
3408        double z, bits;
3409        int qal_display_start, tal_display_start;
3410        static char qal[N], tal[N];
3411        char *qal2, *tal2;
3412        int c;
3413
3414
3415    count = 0;
3416#if 0
3417    for( i=0; i<10000000 && count<nseq; i++ )
3418#else
3419    while( !feof( fp ) )
3420#endif
3421    {
3422        fgets_IGNORE_RESULT( b, B-1, fp );
3423        if( !strncmp( "+==========+", b, 12 ) )
3424        {
3425            junban[count] = atoi( b+12 );
3426
3427                        pt = strchr( b, ')' ) + 1;
3428                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3429            dis[junban[count]] = (double)opt;
3430            count++;
3431
3432        }
3433                else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3434                {
3435                        break;
3436                }
3437
3438    }
3439        if( !count ) return -1;
3440
3441        count = 0;
3442    while( 1 )
3443        {
3444                if( strncmp( ">>+==========+", b, 14 ) )
3445                {
3446                        fgets_IGNORE_RESULT( b, B-1, fp );
3447                        if( feof( fp ) ) break;
3448                        continue;
3449                }
3450                junban[count++] = atoi( b+14 );
3451//              fprintf( stderr, "t = %d\n", atoi( b+14 ) );
3452                while( fgets( b, B-1, fp ) )
3453                        if( !strncmp( "; fa_opt:", b, 9 ) || !strncmp( "; sw_s-w opt:", b, 13 ) ) break;
3454                pt = strstr( b, ":" ) +1;
3455                opt = atoi( pt );
3456
3457
3458                while( fgets( b, B-1, fp ) )
3459                        if( !strncmp( "_overlap:", b+4, 9 ) ) break;
3460                pt = strstr( b, ":" ) +1;
3461                overlapaa = atoi( pt );
3462
3463                while( fgets( b, B-1, fp ) )
3464                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3465                pt = strstr( b, ":" ) +1;
3466                qstart = atoi( pt ) - 1;
3467
3468                while( fgets( b, B-1, fp ) )
3469                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3470                pt = strstr( b, ":" ) +1;
3471                qend = atoi( pt ) - 1;
3472
3473                while( fgets( b, B-1, fp ) )
3474                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3475                pt = strstr( b, ":" ) +1;
3476                qal_display_start = atoi( pt ) - 1;
3477
3478                pt = qal;
3479                while( (c = fgetc( fp )) )
3480                {
3481                        if( c == '>' ) 
3482                        {
3483                                ungetc( c, fp );
3484                                break;
3485                        }
3486                        if( isalpha( c ) || c == '-' ) 
3487                        *pt++ = c;
3488                }
3489                *pt = 0;
3490
3491                while( fgets( b, B-1, fp ) )
3492                        if( !strncmp( "_start:", b+4, 7 ) ) break;
3493                pt = strstr( b, ":" ) + 1;
3494                tstart = atoi( pt ) - 1;
3495
3496                while( fgets( b, B-1, fp ) )
3497                        if( !strncmp( "_stop:", b+4, 6 ) ) break;
3498                pt = strstr( b, ":" ) + 1;
3499                tend = atoi( pt ) - 1;
3500
3501                while( fgets( b, B-1, fp ) )
3502                        if( !strncmp( "_display_start:", b+4, 15 ) ) break;
3503                pt = strstr( b, ":" ) + 1;
3504                tal_display_start = atoi( pt ) - 1;
3505
3506                pt = tal;
3507                while( ( c = fgetc( fp ) ) )
3508                {
3509                        if( c == '>' ) 
3510                        {
3511                                ungetc( c, fp );
3512                                break;
3513                        }
3514                        if( isalpha( c ) || c == '-' ) 
3515                        *pt++ = c;
3516                }
3517                *pt = 0;
3518
3519//              fprintf( stderr, "(%d-%d:%d-%d)\n", qstart, qend, tstart, tend );
3520//              fprintf( stderr, "qal_display_start = %d, tal_display_start = %d\n", qal_display_start, tal_display_start );
3521
3522//              fprintf( stderr, "qal = %s\n", qal );
3523//              fprintf( stderr, "tal = %s\n", tal );
3524
3525                qal2 = cutal( qal, qal_display_start, qstart, qend );
3526                tal2 = cutal( tal, tal_display_start, tstart, tend );
3527
3528//              fprintf( stderr, "qal2 = %s\n", qal2 );
3529//              fprintf( stderr, "tal2 = %s\n", tal2 );
3530
3531//              fprintf( stderr, "putting   %d - %d, opt = %d\n", qmem, junban[count-1], opt );
3532                putlocalhom( qal2, tal2, localhomlist+junban[count-1], qstart, tstart, opt, overlapaa );
3533        }
3534//      fprintf( stderr, "count = %d\n", count );
3535    return count;
3536}
3537int ReadFasta34m10_scoreonly_nucbk( FILE *fp, double *dis, int nin )
3538{
3539    int count=0;
3540    char b[B];
3541        char *pt;
3542    int pos;
3543        int opt;
3544        double z, bits;
3545
3546    count = 0;
3547    while( !feof( fp ) )
3548    {
3549        fgets_IGNORE_RESULT( b, B-1, fp );
3550        if( !strncmp( "+===========+", b, 13 ) )
3551        {
3552            pos = atoi( b+13 );
3553
3554                        if( strchr( b, 'r' ) ) continue;
3555
3556//                      pt = strchr( b, ')' ) + 1;
3557                        pt = strchr( b, ']' ) + 1;
3558                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3559            dis[pos] += (double)opt;
3560            count++;
3561#if 0
3562                        fprintf( stderr, "b=%s\n", b );
3563                        fprintf( stderr, "opt=%d\n", opt );
3564                        fprintf( stderr, "pos=%d\n", pos );
3565                        fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3566#endif
3567
3568        }
3569                else if( 0 == strncmp( ">>><<<", b, 6 ) )
3570                {
3571                        break;
3572                }
3573
3574    }
3575        if( !count ) return -1;
3576
3577    return count;
3578}
3579
3580int ReadFasta34m10_scoreonly_nuc( FILE *fp, double *dis, int nin )
3581{
3582    int count=0;
3583    char b[B];
3584        char *pt;
3585    int pos;
3586        int opt;
3587        double z, bits;
3588        int c;
3589        int *yonda;
3590
3591
3592        yonda = AllocateIntVec( nin );
3593        for( c=0; c<nin; c++ ) yonda[c] = 0;
3594        for( c=0; c<nin; c++ ) dis[c] = 0.0;
3595
3596    count = 0;
3597    while( !feof( fp ) )
3598    {
3599        fgets_IGNORE_RESULT( b, B-1, fp );
3600        if( !strncmp( "+===========+", b, 13 ) )
3601        {
3602            pos = atoi( b+13 );
3603
3604                        if( strchr( b, 'r' ) ) continue;
3605
3606//                      pt = strchr( b, ')' ) + 1;
3607                        pt = strchr( b, ']' ) + 1;
3608                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3609                        if( yonda[pos] == 0 )
3610                        {
3611                    dis[pos] += (double)opt;
3612                                yonda[pos] = 1;
3613                        }
3614            count++;
3615#if 0
3616                        fprintf( stderr, "b=%s\n", b );
3617                        fprintf( stderr, "opt=%d\n", opt );
3618                        fprintf( stderr, "pos=%d\n", pos );
3619                        fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3620#endif
3621
3622        }
3623        else if( !strncmp( ">>>", b, 3 ) )
3624                {
3625                        for( c=0; c<nin; c++ ) yonda[c] = 0;
3626                }
3627                else if( 0 == strncmp( ">>><<<", b, 6 ) )
3628                {
3629                        break;
3630                }
3631
3632    }
3633
3634        free( yonda );
3635
3636        if( !count ) return -1;
3637
3638    return count;
3639}
3640
3641int ReadFasta34m10_scoreonly( FILE *fp, double *dis, int nin )
3642{
3643    int count=0;
3644    char b[B];
3645        char *pt;
3646    int pos;
3647        int opt;
3648        double z, bits;
3649        int c;
3650        int *yonda;
3651
3652
3653        yonda = AllocateIntVec( nin );
3654        for( c=0; c<nin; c++ ) yonda[c] = 0;
3655        for( c=0; c<nin; c++ ) dis[c] = 0.0;
3656
3657    count = 0;
3658    while( !feof( fp ) )
3659    {
3660        fgets_IGNORE_RESULT( b, B-1, fp );
3661        if( !strncmp( "+===========+", b, 13 ) )
3662        {
3663            pos = atoi( b+13 );
3664
3665                        pt = strchr( b, ')' ) + 1;
3666                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3667                        if( yonda[pos] == 0 )
3668                        {
3669                    dis[pos] += (double)opt;
3670                                yonda[pos] = 1;
3671                        }
3672            count++;
3673#if 0
3674                        fprintf( stderr, "b=%s\n", b );
3675                        fprintf( stderr, "opt=%d\n", opt );
3676                        fprintf( stderr, "pos=%d\n", pos );
3677                        fprintf( stderr, "dis[pos]=%f\n", dis[pos] );
3678#endif
3679
3680        }
3681        else if( !strncmp( ">>>", b, 3 ) )
3682                {
3683                        for( c=0; c<nin; c++ ) yonda[c] = 0;
3684                }
3685                else if( 0 == strncmp( ">>><<<", b, 6 ) )
3686                {
3687                        break;
3688                }
3689
3690    }
3691
3692        free( yonda );
3693
3694        if( !count ) return -1;
3695
3696    return count;
3697}
3698int ReadFasta34( FILE *fp, double *dis, int nseq, char name[M][B], LocalHom *localhomlist )
3699{
3700    int count=0;
3701    char b[B];
3702        char *pt;
3703    static int junban[M];
3704        int overlapaa;
3705        int opt, qstart, qend, tstart, tend;
3706        double z, bits;
3707
3708
3709    count = 0;
3710#if 0
3711    for( i=0; i<10000000 && count<nseq; i++ )
3712#else
3713    while( !feof( fp ) )
3714#endif
3715    {
3716        fgets_IGNORE_RESULT( b, B-1, fp );
3717        if( !strncmp( "+==========+", b, 12 ) )
3718        {
3719            junban[count] = atoi( b+12 );
3720
3721                        pt = strchr( b, ')' ) + 1;
3722                        sscanf( pt, "%d %lf %lf",  &opt, &bits, &z ); 
3723            dis[junban[count]] = (double)opt;
3724            count++;
3725
3726        }
3727                else if( 0 == strncmp( ">>+==========+", b, 14 ) )
3728                {
3729                        break;
3730                }
3731
3732    }
3733        if( !count ) return -1;
3734
3735        count = 0;
3736    while( !feof( fp ) )
3737        {
3738                if( !strncmp(">>+==========+", b, 14 ) )
3739                {
3740            junban[count] = atoi( b+14 );
3741            count++;
3742                fgets_IGNORE_RESULT( b, B-1, fp ); // initn:
3743                        pt = strstr( b, "opt: " ) + 5;
3744                        localhomlist[junban[count-1]].opt = atof( pt );
3745                fgets_IGNORE_RESULT( b, B-1, fp ); // Smith-Waterman score
3746                        pt = strstr( b, "ungapped) in " ) + 13;
3747                        sscanf( pt, "%d", &overlapaa ); 
3748                        fprintf( stderr, "pt = %s, overlapaa = %d\n", pt, overlapaa );
3749                        pt = strstr( b, "overlap (" ) + 8;
3750                        sscanf( pt, "(%d-%d:%d-%d)", &qstart, &qend, &tstart, &tend ); 
3751                        localhomlist[junban[count-1]].overlapaa = overlapaa;
3752                        localhomlist[junban[count-1]].start1 = qstart-1;
3753                        localhomlist[junban[count-1]].end1   = qend-1;
3754                        localhomlist[junban[count-1]].start2 = tstart-1;
3755                        localhomlist[junban[count-1]].end2   = tend-1;
3756                }
3757        fgets_IGNORE_RESULT( b, B-1, fp );
3758        }
3759        fprintf( stderr, "count = %d\n", count );
3760    return count;
3761}
3762
3763int ReadFasta3( FILE *fp, double *dis, int nseq, char name[M][B] )
3764{
3765    int count=0;
3766    char b[B];
3767        char *pt;
3768    int junban[M];
3769        int initn, init1, opt;
3770        double z;
3771
3772    count = 0;
3773#if 0
3774    for( i=0; i<10000000 && count<nseq; i++ )
3775#else
3776    while( !feof( fp ) )
3777#endif
3778    {
3779        fgets_IGNORE_RESULT( b, B-1, fp );
3780        if( !strncmp( "+==========+", b, 12 ) )
3781        {
3782            junban[count] = atoi( b+12 );
3783
3784                        pt = strchr( b, ')' ) + 1;
3785                        sscanf( pt, "%d %d %d %lf", &initn, &init1, &opt, &z ); 
3786            dis[junban[count]] = (double)opt;
3787            count++;
3788        }
3789    }
3790    return 0;
3791}
3792
3793int ReadFasta( FILE *fp, double *dis, int nseq, char name[M][B] )
3794{
3795    int i, count=0;
3796    char b[B];
3797    int junban[M];
3798        int initn, init1, opt;
3799
3800    count = 0;
3801        for( i=0; i<nseq; i++ ) dis[i] = 0.0;
3802    for( i=0; !feof( fp ) && count<nseq; i++ )
3803    {
3804        fgets_IGNORE_RESULT( b, B-1, fp );
3805        if( !strncmp( "+==========+", b, 12 ) )
3806        {
3807            junban[count] = atoi( b+12 );
3808
3809                        sscanf( b+50, "%d %d %d", &initn, &init1, &opt ); 
3810            dis[junban[count]] = (double)opt;
3811            count++;
3812        }
3813    }
3814
3815/*
3816    count = 0;
3817    for( i=0; i<100000 && count<nseq; i++ )
3818    {
3819        fgets( b, B-1, fp );
3820        if( !strncmp( name[junban[count]], b, 20  ) )
3821        {
3822            dis[junban[count]] = atof( b+65 );
3823            count++;
3824        }
3825    }
3826*/
3827    return 0;
3828}
3829
3830
3831int ReadOpt( FILE *fp, int opt[M], int nseq, char name[M][B] )
3832{
3833    int i, count=0;
3834    char b[B];
3835    int junban[M];
3836        int optt, initn, init1;
3837
3838    count = 0;
3839    for( i=0; i<10000000 && count<nseq; i++ )
3840    {
3841        fgets_IGNORE_RESULT( b, B-1, fp );
3842        if( !strncmp( "+==========+", b, 12 ) )
3843        {
3844            junban[count] = atoi( b+12 );
3845                        sscanf( b+50, "%d %d %d", &initn, &init1, &optt ); 
3846            opt[junban[count]] = (double)optt;
3847            count++;
3848        }
3849    }
3850    return 0;
3851}
3852
3853int ReadOpt2( FILE *fp, int opt[M], int nseq, char name[M][B] )
3854{
3855    int i, count=0;
3856    char b[B];
3857    int junban[M];
3858
3859    count = 0;
3860    for( i=0; i<10000000 && count<nseq; i++ )
3861    {
3862        fgets_IGNORE_RESULT( b, B-1, fp );
3863        if( !strncmp( "+==========+", b, 12 ) )
3864        {
3865            junban[count] = atoi( b+12 );
3866            opt[junban[count]] = atoi( b+65 );
3867            count++;
3868        }
3869    }
3870    return 0;
3871}
3872
3873
3874
3875int writePre( int nseq, char **name, int nlen[M], char **aseq, int force )
3876{
3877#if USE_XCED
3878        int i, value;
3879        if( !signalSM )
3880        {
3881                if( force ) 
3882                {
3883                        rewind( prep_g );
3884                        if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3885#if 0
3886                        else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3887#else
3888                        else    writeData( prep_g, nseq, name, nlen, aseq );
3889#endif
3890                }
3891                return( 0 );
3892        }
3893        for( i=0; i<10; i++ )
3894        {
3895#if IODEBUG
3896                fprintf( stderr, "SEMAPHORE = %d\n", signalSM[SEMAPHORE] );
3897#endif
3898                if( signalSM[SEMAPHORE]-- > 0 )
3899                {
3900#if 0 /* /tmp/pre €ÎŽØ·ž€Ç€Ï€º€·€¿ */
3901                        if( ferror( prep_g ) ) prep_g = fopen( "pre", "w" );
3902                        if( !prep_g ) ErrorExit( "Cannot re-open pre." );
3903#endif
3904                        rewind( prep_g );
3905                        signalSM[STATUS] = IMA_KAITERU;
3906#if IODEBUG
3907                        if( force ) fprintf( stderr, "FINAL " );
3908#endif
3909                        if( devide ) dvWrite( prep_g, nseq, name, nlen, aseq );
3910                        else    WriteGapFill( prep_g, nseq, name, nlen, aseq );
3911                        /*
3912                        fprintf( prep_g, '\EOF' );
3913                        */
3914                        fflush( prep_g );
3915                        if( force ) signalSM[STATUS] = OSHIMAI;
3916                        else        signalSM[STATUS] = KAKIOWATTA;
3917                        value = 1;
3918                        signalSM[SEMAPHORE]++;
3919#if IODEBUG
3920                        fprintf( stderr, "signalSM[STATUS] = %c\n", signalSM[STATUS] );
3921#endif
3922                        break;
3923                }
3924                else
3925                {
3926#if IODEBUG
3927                        fprintf( stderr, "YONDERUKARA_AKIRAMERU\n" );
3928#endif
3929                        value = 0;
3930                        signalSM[SEMAPHORE]++;
3931                        if( !force ) break;
3932#if IODEBUG
3933                        fprintf( stderr, "MATSU\n" );
3934#endif
3935                        sleep( 1 );
3936                }
3937        }
3938        if( force && !value ) ErrorExit( "xced ga pre wo hanasanai \n" );
3939        return( value );
3940#else
3941        if( force ) 
3942        {
3943                rewind( prep_g );
3944                writeData_pointer( prep_g, nseq, name, nlen, aseq );
3945        }
3946#endif
3947        return( 0 );
3948}
3949
3950
3951void readOtherOptions( int *ppidptr, int *fftThresholdptr, int *fftWinSizeptr )
3952{
3953        if( calledByXced )
3954        {
3955                FILE *fp = fopen( "pre", "r" );
3956                char b[B];
3957                if( !fp ) ErrorExit( "Cannot open pre.\n" );
3958                fgets_IGNORE_RESULT( b, B-1, fp );
3959                sscanf( b, "%d %d %d", ppidptr, fftThresholdptr, fftWinSizeptr );
3960                fclose( fp );
3961#if IODEBUG
3962        fprintf( stderr, "b = %s\n", b );
3963        fprintf( stderr, "ppid = %d\n", ppid );
3964        fprintf( stderr, "fftThreshold = %d\n", fftThreshold );
3965        fprintf( stderr, "fftWinSize = %d\n", fftWinSize );
3966#endif
3967        }
3968        else
3969        {
3970                *ppidptr = 0;
3971                *fftThresholdptr = FFT_THRESHOLD;
3972                if( dorp == 'd' )
3973                        *fftWinSizeptr = FFT_WINSIZE_D;
3974                else
3975                        *fftWinSizeptr = FFT_WINSIZE_P;
3976        }
3977#if 0
3978        fprintf( stderr, "fftThresholdptr=%d\n", *fftThresholdptr );
3979        fprintf( stderr, "fftWinSizeptr=%d\n", *fftWinSizeptr );
3980#endif
3981}
3982
3983void initSignalSM( void )
3984{
3985//      int signalsmid;
3986
3987#if IODEBUG
3988        if( ppid ) fprintf( stderr, "PID of xced = %d\n", ppid );
3989#endif
3990        if( !ppid )
3991        {
3992                signalSM = NULL;
3993                return;
3994        }
3995
3996#if 0
3997        signalsmid = shmget( (key_t)ppid, 3, IPC_ALLOC | 0666 );
3998        if( signalsmid == -1 ) ErrorExit( "Cannot get Shared memory for signal.\n" );
3999        signalSM = shmat( signalsmid, 0, 0 );
4000        if( (int)signalSM == -1 ) ErrorExit( "Cannot attatch Shared Memory for signal!\n" );
4001        signalSM[STATUS] = IMA_KAITERU;
4002        signalSM[SEMAPHORE] = 1;
4003#endif
4004}
4005
4006void initFiles( void )
4007{
4008        char pname[100];
4009        if( ppid )
4010                sprintf( pname, "/tmp/pre.%d", ppid );
4011        else
4012                sprintf( pname, "pre" );
4013        prep_g = fopen( pname, "w" );
4014        if( !prep_g ) ErrorExit( "Cannot open pre" );
4015
4016        trap_g = fopen( "trace", "w" );
4017        if( !trap_g ) ErrorExit( "cannot open trace" );
4018        fprintf( trap_g, "PID = %d\n", getpid() );
4019        fflush( trap_g );
4020}
4021
4022
4023void WriteForFasta( FILE *fp, int locnjob, char **name, int nlen[M], char **aseq )
4024{
4025    static char b[N];
4026    int i, j;
4027    int nalen[M];
4028
4029    for( i=0; i<locnjob; i++ )
4030    {
4031        nalen[i] = strlen( aseq[i] );
4032        fprintf( fp, ">%s\n", name[i] );
4033        for( j=0; j<nalen[i]; j=j+C ) 
4034        {
4035            strncpy( b, aseq[i]+j, C ); b[C] = 0;
4036            fprintf( fp, "%s\n",b );
4037        }
4038    }
4039}
4040
4041void readlocalhomtable2( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
4042{
4043        double opt;
4044        static char buff[B];
4045        char infor[100];
4046        int i, j, overlapaa, start1, end1, start2, end2;
4047        LocalHom *tmpptr1, *tmpptr2;
4048
4049//      for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
4050
4051        while ( NULL != fgets( buff, B-1, fp ) )
4052        {
4053//              fprintf( stderr, "\n" );
4054                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4055                if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
4056
4057#if 0
4058                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4059#endif
4060
4061//              if( i < j )
4062                {
4063                        if( localhomtable[i][j].nokori++ > 0 )
4064                        {
4065                                tmpptr1 = localhomtable[i][j].last;
4066//                              fprintf( stderr, "reallocating, localhomtable[%d][%d].nokori = %d\n", i, j, localhomtable[i][j].nokori );
4067                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4068                                tmpptr1 = tmpptr1->next;
4069                                tmpptr1->extended = -1;
4070                                tmpptr1->next = NULL;
4071                                localhomtable[i][j].last = tmpptr1;
4072//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
4073                        }
4074                        else
4075                        {
4076                                tmpptr1 = localhomtable[i]+j;
4077//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", i, j, localhomtable[i][j].nokori );
4078                        }
4079       
4080                        tmpptr1->start1 = start1;
4081                        tmpptr1->start2 = start2;
4082                        tmpptr1->end1 = end1;
4083                        tmpptr1->end2 = end2;
4084//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4085//                      tmpptr1->opt = opt;
4086                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4087                        tmpptr1->overlapaa = overlapaa;
4088                        tmpptr1->korh = *infor;
4089                }
4090//              else
4091                {
4092                        if( localhomtable[j][i].nokori++ > 0 )
4093                        {
4094                                tmpptr2 = localhomtable[j][i].last;
4095                                tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4096                                tmpptr2 = tmpptr2->next;
4097                                tmpptr2->extended = -1;
4098                                tmpptr2->next = NULL;
4099                                localhomtable[j][i].last = tmpptr2;
4100//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
4101                        }
4102                        else
4103                        {
4104                                tmpptr2 = localhomtable[j]+i;
4105//                              fprintf( stderr, "### i,j = %d,%d, nokori=%d\n", j, i, localhomtable[j][i].nokori );
4106                        }
4107       
4108                        tmpptr2->start2 = start1;
4109                        tmpptr2->start1 = start2;
4110                        tmpptr2->end2 = end1;
4111                        tmpptr2->end1 = end2;
4112//                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4113//                      tmpptr2->opt = opt;
4114                        tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4115                        tmpptr2->overlapaa = overlapaa;
4116                        tmpptr2->korh = *infor;
4117       
4118//                      fprintf( stderr, "i=%d, j=%d, st1=%d, en1=%d, opt = %f\n", i, j, tmpptr1->start1, tmpptr1->end1, opt );
4119                }
4120
4121        }
4122}
4123void readlocalhomtable( FILE*fp, int njob, LocalHom **localhomtable, char *kozoarivec )
4124{
4125        double opt;
4126        static char buff[B];
4127        char infor[100];
4128        int i, j, overlapaa, start1, end1, start2, end2;
4129        int **nlocalhom = NULL;
4130        LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
4131
4132        nlocalhom = AllocateIntMtx( njob, njob );
4133        for( i=0; i<njob; i++ ) for( j=0; j<njob; j++ ) nlocalhom[i][j] = 0;
4134
4135        while ( NULL != fgets( buff, B-1, fp ) )
4136        {
4137//              fprintf( stderr, "\n" );
4138                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4139                if( *infor == 'k' ) kozoarivec[i] = kozoarivec[j] = 1;
4140
4141#if 0
4142                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4143#endif
4144
4145
4146//              if( i < j )
4147                {
4148                        if( nlocalhom[i][j]++ > 0 )
4149                        {
4150//                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4151                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4152                                tmpptr1 = tmpptr1->next;
4153                                tmpptr1->next = NULL;
4154                        }
4155                        else
4156                        {
4157                                tmpptr1 = localhomtable[i]+j;
4158//                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4159                        }
4160       
4161                        tmpptr1->start1 = start1; // CHUUI!!!!
4162                        tmpptr1->start2 = start2;
4163                        tmpptr1->end1 = end1; // CHUUI!!!!
4164                        tmpptr1->end2 = end2;
4165//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4166//                      tmpptr1->opt = opt;
4167                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4168                        tmpptr1->overlapaa = overlapaa;
4169                        tmpptr1->korh = *infor;
4170       
4171//                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4172                }
4173//              else
4174                {
4175                        if( nlocalhom[j][i]++ > 0 )
4176                        {
4177                                tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4178                                tmpptr2 = tmpptr2->next;
4179                                tmpptr2->next = NULL;
4180                        }
4181                        else
4182                                tmpptr2 = localhomtable[j]+i;
4183       
4184                        tmpptr2->start2 = start1; // CHUUI!!!!
4185                        tmpptr2->start1 = start2;
4186                        tmpptr2->end2 = end1; // CHUUI!!!!
4187                        tmpptr2->end1 = end2;
4188//                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4189//                      tmpptr2->opt = opt;
4190                        tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4191                        tmpptr2->overlapaa = overlapaa;
4192                        tmpptr2->korh = *infor;
4193
4194//                      fprintf( stderr, "j=%d, i=%d, opt = %f\n", j, i, opt );
4195                }
4196
4197        }
4198        FreeIntMtx( nlocalhom );
4199}
4200
4201
4202void readlocalhomtable_two( FILE*fp, int norg, int nadd, LocalHom **localhomtable, LocalHom **localhomtablex, char *kozoarivec ) // for test only
4203{
4204        double opt;
4205        static char buff[B];
4206        char infor[100];
4207        int i, j, overlapaa, start1, end1, start2, end2;
4208        int **nlocalhom = NULL;
4209        int **nlocalhomx = NULL;
4210        LocalHom *tmpptr1=NULL, *tmpptr2=NULL; // by D.Mathog, a guess
4211
4212        nlocalhom = AllocateIntMtx( norg, nadd );
4213        for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ ) nlocalhom[i][j] = 0;
4214        nlocalhomx = AllocateIntMtx( nadd, norg );
4215        for( i=0; i<nadd; i++ ) for( j=0; j<norg; j++ ) nlocalhomx[i][j] = 0;
4216
4217        while ( NULL != fgets( buff, B-1, fp ) )
4218        {
4219//              fprintf( stderr, "\n" );
4220                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4221                if( *infor == 'k' ) 
4222                {
4223                        fprintf( stderr, "Not supported!\n" );
4224                        exit( 1 );
4225                }
4226                j -= norg;
4227
4228#if 0
4229                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4230#endif
4231
4232
4233//              if( i < j )
4234                {
4235                        if( nlocalhom[i][j]++ > 0 )
4236                        {
4237//                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4238                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4239                                tmpptr1 = tmpptr1->next;
4240                                tmpptr1->next = NULL;
4241                        }
4242                        else
4243                        {
4244                                tmpptr1 = localhomtable[i]+j;
4245//                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4246                        }
4247       
4248                        tmpptr1->start1 = start1; // CHUUI!!!!
4249                        tmpptr1->start2 = start2;
4250                        tmpptr1->end1 = end1; // CHUUI!!!!
4251                        tmpptr1->end2 = end2;
4252//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4253//                      tmpptr1->opt = opt;
4254                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4255                        tmpptr1->overlapaa = overlapaa;
4256                        tmpptr1->korh = *infor;
4257       
4258//                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4259                }
4260
4261                {
4262                        if( nlocalhomx[j][i]++ > 0 )
4263                        {
4264                                tmpptr2->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4265                                tmpptr2 = tmpptr2->next;
4266                                tmpptr2->next = NULL;
4267                        }
4268                        else
4269                                tmpptr2 = localhomtablex[j]+i;
4270       
4271                        tmpptr2->start2 = start1+1; // CHUUI!!!!
4272                        tmpptr2->start1 = start2;
4273                        tmpptr2->end2 = end1+1; // CHUUI!!!!
4274                        tmpptr2->end1 = end2;
4275//                      tmpptr2->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4276//                      tmpptr2->opt = opt;
4277                        tmpptr2->opt = ( opt + 0.00 ) / 5.8  * 600;
4278                        tmpptr2->overlapaa = overlapaa;
4279                        tmpptr2->korh = *infor;
4280
4281//                      fprintf( stderr, "j=%d, i=%d, opt = %f\n", j, i, opt );
4282                }
4283
4284        }
4285        FreeIntMtx( nlocalhom );
4286        FreeIntMtx( nlocalhomx );
4287}
4288
4289void readlocalhomtable_one( FILE*fp, int norg, int nadd, LocalHom **localhomtable, char *kozoarivec ) // for test only
4290{
4291        double opt;
4292        static char buff[B];
4293        char infor[100];
4294        int i, j, overlapaa, start1, end1, start2, end2;
4295        int **nlocalhom = NULL;
4296        LocalHom *tmpptr1=NULL; // by D.Mathog, a guess
4297
4298        nlocalhom = AllocateIntMtx( norg, nadd );
4299        for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ ) nlocalhom[i][j] = 0;
4300
4301        while ( NULL != fgets( buff, B-1, fp ) )
4302        {
4303//              fprintf( stderr, "\n" );
4304                sscanf( buff, "%d %d %d %lf %d %d %d %d %s",  &i, &j, &overlapaa, &opt, &start1, &end1, &start2, &end2, infor );
4305                if( *infor == 'k' ) 
4306                {
4307                        fprintf( stderr, "Not supported!\n" );
4308                        exit( 1 );
4309                }
4310                j -= norg;
4311
4312#if 0
4313                if( start1 == end1 || start2 == end2 ) continue; //mondai ari
4314#endif
4315
4316
4317//              if( i < j )
4318                {
4319                        if( nlocalhom[i][j]++ > 0 )
4320                        {
4321//                              fprintf( stderr, "reallocating, nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4322                                tmpptr1->next = (LocalHom *)calloc( 1, sizeof( LocalHom ) );
4323                                tmpptr1 = tmpptr1->next;
4324                                tmpptr1->next = NULL;
4325                        }
4326                        else
4327                        {
4328                                tmpptr1 = localhomtable[i]+j;
4329//                              fprintf( stderr, "nlocalhom[%d][%d] = %d\n", i, j, nlocalhom[i][j] );
4330                        }
4331       
4332                        tmpptr1->start1 = start1; // CHUUI!!!!
4333                        tmpptr1->start2 = start2;
4334                        tmpptr1->end1 = end1; // CHUUI!!!!
4335                        tmpptr1->end2 = end2;
4336//                      tmpptr1->opt = ( opt / overlapaa + 0.00 ) / 5.8  * 600;
4337//                      tmpptr1->opt = opt;
4338                        tmpptr1->opt = ( opt + 0.00 ) / 5.8  * 600;
4339                        tmpptr1->overlapaa = overlapaa;
4340                        tmpptr1->korh = *infor;
4341       
4342//                      fprintf( stderr, "i=%d, j=%d, opt = %f\n", i, j, opt );
4343                }
4344
4345        }
4346        FreeIntMtx( nlocalhom );
4347}
4348
4349void outlocalhom_part( LocalHom **localhom, int norg, int nadd )
4350{
4351        int i, j;
4352        LocalHom *tmpptr;
4353        for( i=0; i<norg; i++ ) for( j=0; j<nadd; j++ )
4354        {
4355                tmpptr = localhom[i]+j;
4356                fprintf( stderr, "%d-%d\n", i, j+norg );
4357                do
4358                {
4359                        fprintf( stderr, "reg1=%d-%d, reg2=%d-%d, imp=%f, opt=%f\n", tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->importance, tmpptr->opt );
4360                }
4361                while( (tmpptr=tmpptr->next) );
4362        }
4363}
4364
4365void outlocalhom( LocalHom **localhom, int nseq )
4366{
4367        int i, j;
4368        LocalHom *tmpptr;
4369        for( i=0; i<nseq; i++ ) for( j=0; j<nseq; j++ )
4370        {
4371                tmpptr = localhom[i]+j;
4372                fprintf( stderr, "%d-%d\n", i, j );
4373                do
4374                {
4375                        fprintf( stderr, "reg1=%d-%d, reg2=%d-%d, imp=%f, opt=%f\n", tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->importance, tmpptr->opt );
4376                }
4377                while( (tmpptr=tmpptr->next) );
4378        }
4379}
4380
4381void outlocalhompt( LocalHom ***localhom, int n1, int n2 )
4382{
4383        int i, j;
4384        LocalHom *tmpptr;
4385        for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ )
4386        {
4387                tmpptr = localhom[i][j];
4388//              fprintf( stdout, "%d-%d\n", i, j );
4389                do
4390                {
4391                        fprintf( stdout, "%d-%d, reg1=%d-%d, reg2=%d-%d, imp=%f, opt=%f, wimp=%f\n", i, j, tmpptr->start1, tmpptr->end1, tmpptr->start2, tmpptr->end2, tmpptr->importance, tmpptr->opt, tmpptr->wimportance );
4392                }
4393                while( (tmpptr=tmpptr->next) );
4394        }
4395}
4396
4397void FreeLocalHomTable_part( LocalHom **localhomtable, int n, int m ) 
4398{
4399        int i, j;
4400        LocalHom *ppp, *tmpptr;
4401        for( i=0; i<n; i++ ) 
4402        {
4403                for( j=0; j<m; j++ )
4404                {
4405                        tmpptr=localhomtable[i]+j;
4406                        ppp = tmpptr->next;
4407                        for( ; tmpptr; tmpptr=ppp )
4408                        {
4409#if DEBUG
4410                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4411#endif
4412                                ppp = tmpptr->next;
4413                                if( tmpptr!=localhomtable[i]+j ) 
4414                                {
4415#if DEBUG
4416                                        fprintf( stderr, "freeing %p\n", tmpptr );
4417#endif
4418                                        free( tmpptr );
4419                                }
4420                        }
4421                }
4422#if DEBUG
4423                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4424#endif
4425                free( localhomtable[i] );
4426        }
4427#if DEBUG
4428        fprintf( stderr, "freeing localhomtable\n" );
4429#endif
4430        free( localhomtable );
4431#if DEBUG
4432        fprintf( stderr, "freed\n" );
4433#endif
4434}
4435
4436void FreeLocalHomTable_two( LocalHom **localhomtable, int n, int m ) 
4437{
4438        int i, j;
4439        LocalHom *ppp, *tmpptr;
4440        for( i=0; i<n; i++ ) 
4441        {
4442                for( j=0; j<m; j++ )
4443                {
4444                        tmpptr=localhomtable[i]+j;
4445                        ppp = tmpptr->next;
4446                        for( ; tmpptr; tmpptr=ppp )
4447                        {
4448#if DEBUG
4449                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4450#endif
4451                                ppp = tmpptr->next;
4452                                if( tmpptr!=localhomtable[i]+j ) 
4453                                {
4454#if DEBUG
4455                                        fprintf( stderr, "freeing %p\n", tmpptr );
4456#endif
4457                                        free( tmpptr );
4458                                }
4459                        }
4460                }
4461#if DEBUG
4462                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4463#endif
4464                free( localhomtable[i] );
4465        }
4466
4467        for( i=n; i<n+m; i++ ) 
4468        {
4469                for( j=0; j<n; j++ )
4470                {
4471                        tmpptr=localhomtable[i]+j;
4472                        ppp = tmpptr->next;
4473                        for( ; tmpptr; tmpptr=ppp )
4474                        {
4475#if DEBUG
4476                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4477#endif
4478                                ppp = tmpptr->next;
4479                                if( tmpptr!=localhomtable[i]+j ) 
4480                                {
4481#if DEBUG
4482                                        fprintf( stderr, "freeing %p\n", tmpptr );
4483#endif
4484                                        free( tmpptr );
4485                                }
4486                        }
4487                }
4488#if DEBUG
4489                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4490#endif
4491                free( localhomtable[i] );
4492        }
4493#if DEBUG
4494        fprintf( stderr, "freeing localhomtable\n" );
4495#endif
4496        free( localhomtable );
4497#if DEBUG
4498        fprintf( stderr, "freed\n" );
4499#endif
4500}
4501
4502void FreeLocalHomTable_one( LocalHom **localhomtable, int n, int m ) 
4503{
4504        int i, j;
4505        LocalHom *ppp, *tmpptr;
4506        for( i=0; i<n; i++ ) 
4507        {
4508                for( j=0; j<m; j++ )
4509                {
4510                        tmpptr=localhomtable[i]+j;
4511                        ppp = tmpptr->next;
4512                        for( ; tmpptr; tmpptr=ppp )
4513                        {
4514#if DEBUG
4515                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4516#endif
4517                                ppp = tmpptr->next;
4518                                if( tmpptr!=localhomtable[i]+j ) 
4519                                {
4520#if DEBUG
4521                                        fprintf( stderr, "freeing %p\n", tmpptr );
4522#endif
4523                                        free( tmpptr );
4524                                }
4525                        }
4526                }
4527#if DEBUG
4528                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4529#endif
4530                free( localhomtable[i] );
4531        }
4532
4533#if DEBUG
4534        fprintf( stderr, "freeing localhomtable\n" );
4535#endif
4536        free( localhomtable );
4537#if DEBUG
4538        fprintf( stderr, "freed\n" );
4539#endif
4540}
4541
4542void FreeLocalHomTable( LocalHom **localhomtable, int n ) 
4543{
4544        int i, j;
4545        LocalHom *ppp, *tmpptr;
4546        for( i=0; i<n; i++ ) 
4547        {
4548                for( j=0; j<n; j++ )
4549                {
4550                        tmpptr=localhomtable[i]+j;
4551                        ppp = tmpptr->next;
4552                        for( ; tmpptr; tmpptr=ppp )
4553                        {
4554#if DEBUG
4555                                fprintf( stderr, "i=%d, j=%d\n", i, j ); 
4556#endif
4557                                ppp = tmpptr->next;
4558                                if( tmpptr!=localhomtable[i]+j ) 
4559                                {
4560#if DEBUG
4561                                        fprintf( stderr, "freeing %p\n", tmpptr );
4562#endif
4563                                        free( tmpptr );
4564                                }
4565                        }
4566                }
4567#if DEBUG
4568                fprintf( stderr, "freeing localhomtable[%d]\n", i );
4569#endif
4570                free( localhomtable[i] );
4571        }
4572#if DEBUG
4573        fprintf( stderr, "freeing localhomtable\n" );
4574#endif
4575        free( localhomtable );
4576#if DEBUG
4577        fprintf( stderr, "freed\n" );
4578#endif
4579}
4580
4581char *progName( char *str )
4582{
4583    char *value; 
4584    if( ( value = strrchr( str, '/' ) ) != NULL )
4585        return( value+1 );
4586    else   
4587        return( str );
4588}
4589
4590static void tabtospace( char *str )
4591{
4592        char *p;
4593//      fprintf( stderr, "before = %s\n", str );
4594        while( NULL != ( p = strchr( str , '\t' ) ) )
4595        {
4596                *p = ' ';
4597        }
4598//      fprintf( stderr, "after = %s\n", str );
4599}
4600
4601static char *extractfirstword( char *str )
4602{
4603        char *val = str;
4604
4605        tabtospace( str );
4606        while( *str )
4607        {
4608                if( val == str && *str == ' ' )
4609                {
4610                        val++; str++;
4611                }
4612                else if( *str != ' ' )
4613                {
4614                        str++;
4615                }
4616                else if( *str == ' ' )
4617                {
4618                        *str = 0;
4619                }
4620        }
4621        return( val );
4622}
4623
4624void phylipout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, int *order, int namelen )
4625{
4626        int pos, pos2, j;
4627        if( namelen == -1 ) namelen = 10;
4628        pos = 0;
4629
4630        fprintf( fp, " %d %d\n", nseq, maxlen );
4631       
4632        while( pos < maxlen )
4633        {
4634                for( j=0; j<nseq; j++ )
4635                {
4636                        if( pos == 0 )
4637                                fprintf( fp, "%-*.*s", namelen, namelen, extractfirstword( name[order[j]]+1 ) );
4638                        else
4639                                fprintf( fp, "%-*.*s", namelen, namelen, "" );
4640
4641                        pos2 = pos;
4642                        while( pos2 < pos+41 && pos2 < maxlen )
4643                        {
4644                                fprintf( fp, " %.10s", seq[order[j]]+pos2 );
4645                                pos2 += 10;
4646                        }
4647                        fprintf( fp, "\n" );
4648                }
4649                fprintf( fp, "\n" );
4650                pos += 50;
4651        }
4652}
4653
4654void clustalout_pointer( FILE *fp, int nseq, int maxlen, char **seq, char **name, char *mark, char *comment, int *order, int namelen )
4655{
4656        int pos, j;
4657        if( namelen == -1 ) namelen = 15;
4658        pos = 0;
4659        if( comment == NULL )
4660                fprintf( fp, "CLUSTAL format alignment by MAFFT (v%s)\n\n", VERSION );
4661        else
4662                fprintf( fp, "CLUSTAL format alignment by MAFFT %s (v%s)\n\n", comment, VERSION );
4663       
4664        while( pos < maxlen )
4665        {
4666                fprintf( fp, "\n" );
4667                for( j=0; j<nseq; j++ )
4668                {
4669                        fprintf( fp, "%-*.*s ", namelen, namelen, extractfirstword( name[order[j]]+1 ) );
4670                        fprintf( fp, "%.60s\n", seq[order[j]]+pos ); // Ĺ€µ€¬°ã€Š€È€À€á¡£
4671                }
4672                if( mark )
4673                {
4674                        fprintf( fp, "%-*.*s ", namelen, namelen, "" );
4675                        fprintf( fp, "%.60s\n", mark + pos ); // Ĺ€µ€¬°ã€Š€È€À€á¡£
4676                }
4677                pos += 60;
4678        }
4679}
4680
4681
4682void writeData_reorder_pointer( FILE *fp, int locnjob, char **name, int *nlen, char **aseq, int *order )
4683{
4684        int i, j, k;
4685        int nalen;
4686
4687        for( i=0; i<locnjob; i++ )
4688        {
4689                k = order[i];
4690#if DEBUG
4691                fprintf( stderr, "i = %d in writeData\n", i );
4692#endif
4693                nalen = strlen( aseq[k] );
4694                fprintf( fp, ">%s\n", name[k]+1 );
4695                for( j=0; j<nalen; j=j+C )
4696                {
4697#if 0
4698                        strncpy( b, aseq[k]+j, C ); b[C] = 0;
4699                        fprintf( fp, "%s\n",b );
4700#else
4701                        fprintf( fp, "%.*s\n", C, aseq[k]+j );
4702#endif
4703                }
4704        }
4705}
4706void writeData_reorder( FILE *fp, int locnjob, char name[][B], int nlen[], char **aseq, int *order )
4707{
4708        int i, j, k;
4709        int nalen;
4710
4711        for( i=0; i<locnjob; i++ )
4712        {
4713                k = order[i];
4714#if DEBUG
4715                fprintf( stderr, "i = %d in writeData\n", i );
4716#endif
4717                nalen = strlen( aseq[k] );
4718                fprintf( fp, ">%s\n", name[k]+1 );
4719                for( j=0; j<nalen; j=j+C )
4720                {
4721#if 0
4722                        strncpy( b, aseq[k]+j, C ); b[C] = 0;
4723                        fprintf( fp, "%s\n",b );
4724#else
4725                        fprintf( fp, "%.*s\n", C, aseq[k]+j );
4726#endif
4727                }
4728        }
4729}
4730static void showaamtxexample()
4731{
4732        fprintf( stderr, "Format error in aa matrix\n" );
4733        fprintf( stderr, "# Example:\n" );
4734        fprintf( stderr, "# comment\n" );
4735        fprintf( stderr, "   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V\n" );
4736        fprintf( stderr, "A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0\n" );
4737        fprintf( stderr, "R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3\n" );
4738        fprintf( stderr, "...\n" );
4739        fprintf( stderr, "V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4\n" );
4740        fprintf( stderr, "frequency 0.07 0.05 0.04 0.05 0.02 .. \n" );
4741        fprintf( stderr, "# Example end\n" );
4742        fprintf( stderr, "Only the lower half is loaded\n" );
4743        fprintf( stderr, "The last line (frequency) is optional.\n" );
4744        exit( 1 );
4745}
4746
4747double *loadaamtx( void )
4748{
4749        int i, j, k, ii, jj;
4750        double *val;
4751        double **raw;
4752        int *map;
4753        char *aaorder = "ARNDCQEGHILKMFPSTWYV";
4754        char *inorder;
4755        char *line;
4756        char *ptr1;
4757        char *ptr2;
4758        char *mtxfname = "_aamtx";
4759        FILE *mf;
4760
4761        raw = AllocateDoubleMtx( 21, 20 );
4762        val = AllocateDoubleVec( 420 );
4763        map = AllocateIntVec( 20 );
4764
4765        if( dorp != 'p' )
4766        {
4767                fprintf( stderr, "User-defined matrix is not supported for DNA\n" );
4768                exit( 1 );
4769        }
4770
4771        mf = fopen( mtxfname, "r" );
4772        if( mf == NULL )
4773        {
4774                fprintf( stderr, "Cannot open the _aamtx file\n" );
4775                exit( 1 );
4776        }
4777
4778        inorder = calloc( 1000, sizeof( char ) );
4779        line = calloc( 1000, sizeof( char ) );
4780       
4781
4782        while( !feof( mf ) )
4783        {
4784                fgets_IGNORE_RESULT( inorder, 999, mf );
4785                if( inorder[0] != '#' ) break;
4786        }
4787        ptr1 = ptr2 = inorder;
4788        while( *ptr2 )
4789        {
4790                if( isalpha( *ptr2 ) )
4791                {
4792                        *ptr1 = toupper( *ptr2 );
4793                        ptr1++;
4794                }
4795                ptr2++;
4796        }
4797        inorder[20] = 0;
4798
4799        for( i=0; i<20; i++ )
4800        {
4801                ptr2 = strchr( inorder, aaorder[i] );
4802                if( ptr2 == NULL )
4803                {
4804                        fprintf( stderr, "%c: not found in the first 20 letters.\n", aaorder[i] );
4805                        showaamtxexample();
4806                }
4807                else
4808                {
4809                        map[i] = ptr2 - inorder;
4810                }
4811        }
4812
4813        i = 0;
4814        while( !feof( mf ) )
4815        {
4816                fgets_IGNORE_RESULT( line, 999, mf );
4817//              fprintf( stderr, "line = %s\n", line );
4818                if( line[0] == '#' ) continue;
4819                ptr1 = line;
4820//              fprintf( stderr, "line = %s\n", line );
4821                for( j=0; j<=i; j++ )
4822                {
4823                        while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4824                                ptr1++;
4825
4826                        raw[i][j] = atof( ptr1 );
4827//                      fprintf( stderr, "raw[][]=%f, %c-%c %d-%d\n", raw[i][j], inorder[i], inorder[j], i, j );
4828                        ptr1 = strchr( ptr1, ' ' );
4829                        if( ptr1 == NULL && j<i) showaamtxexample();
4830                }
4831                i++;
4832                if( i > 19 ) break;
4833        }
4834
4835        for( i=0; i<20; i++ ) raw[20][i] = -1.0;
4836        while( !feof( mf ) )
4837        {
4838                fgets_IGNORE_RESULT( line, 999, mf );
4839                if( line[0] == 'f' )
4840                {
4841//                      fprintf( stderr, "line = %s\n", line );
4842                        ptr1 = line;
4843                        for( j=0; j<20; j++ )
4844                        {
4845                                while( !isdigit( *ptr1 ) && *ptr1 != '-' && *ptr1 != '.' )
4846                                        ptr1++;
4847       
4848                                raw[20][j] = atof( ptr1 );
4849//                              fprintf( stderr, "raw[20][]=%f, %c %d\n", raw[20][j], inorder[i], j );
4850                                ptr1 = strchr( ptr1, ' ' );
4851                                if( ptr1 == NULL && j<19) showaamtxexample();
4852                        }
4853                        break;
4854                }
4855        }
4856
4857        k = 0;
4858        for( i=0; i<20; i++ )
4859        {
4860                for( j=0; j<=i; j++ )
4861                {
4862                        if( i != j )
4863                        {
4864                                ii = MAX( map[i], map[j] );
4865                                jj = MIN( map[i], map[j] );
4866                        }
4867                        else ii = jj = map[i];
4868                        val[k++] = raw[ii][jj];
4869//                      fprintf( stderr, "%c-%c, %f\n", aaorder[i], aaorder[j], val[k-1] );
4870                }
4871        }
4872        for( i=0; i<20; i++ ) val[400+i] = raw[20][map[i]];
4873
4874        fprintf( stderr, "inorder = %s\n", inorder );
4875        fclose( mf );
4876        free( inorder );
4877        free( line );
4878        FreeDoubleMtx( raw );
4879        free( map );
4880        return( val );
4881}
4882
4883static void tab2space( char *s ) // nen no tame
4884{
4885        while( *s )
4886        {
4887                if( *s == '\t' ) *s = ' ';
4888                s++;
4889        }
4890}
4891
4892static int readasubalignment( char *s, int *t, int *preservegaps )
4893{
4894        int v = 0;
4895        char status = 's';
4896        char *pt = s;
4897        *preservegaps = 0;
4898        tab2space( s );
4899        while( *pt )
4900        {
4901                if( *pt == ' ' )
4902                {
4903                        status = 's';
4904                }
4905                else
4906                {
4907                        if( status == 's' ) 
4908                        {
4909                                if( *pt == '\n' || *pt == '#' ) break;
4910                                status = 'n';
4911                                t[v] = atoi( pt );
4912                                if( t[v] == 0 )
4913                                {
4914                                        fprintf( stderr, "Format error? Sequences must be specified as 1, 2, 3...\n" );
4915                                        exit( 1 );
4916                                }
4917                                if( t[v] < 0 ) *preservegaps = 1;
4918                                t[v] = abs( t[v] );
4919                                t[v] -= 1;
4920                                v++;
4921                        }
4922                }
4923                pt++;
4924        }
4925        t[v] = -1;
4926        return( v );
4927}
4928
4929static int countspace( char *s )
4930{
4931        int v = 0;
4932        char status = 's';
4933        char *pt = s;
4934        tab2space( s );
4935        while( *pt )
4936        {
4937                if( *pt == ' ' )
4938                {
4939                        status = 's';
4940                }
4941                else
4942                {
4943                        if( status == 's' ) 
4944                        {
4945                                if( *pt == '\n' || *pt == '#' ) break;
4946                                v++;
4947                                status = 'n';
4948                                if( atoi( pt ) == 0 )
4949                                {
4950                                        fprintf( stderr, "Format error? Sequences should be specified as 1, 2, 3...\n" );
4951                                        exit( 1 );
4952                                }
4953                        }
4954                }
4955                pt++;
4956        }
4957        return( v );
4958}
4959
4960
4961void readsubalignmentstable( int nseq, int **table, int *preservegaps, int *nsubpt, int *maxmempt )
4962{
4963        FILE *fp;
4964        char *line;
4965        int linelen = 1000000;
4966        int nmem;
4967        int lpos;
4968        int i, p;
4969        int *tab01;
4970
4971        line = calloc( linelen, sizeof( char ) );
4972        fp = fopen( "_subalignmentstable", "r" );
4973        if( !fp )
4974        {
4975                fprintf( stderr, "Cannot open _subalignmentstable\n" );
4976                exit( 1 );
4977        }
4978        if( table == NULL )
4979        {
4980                *nsubpt = 0;
4981                *maxmempt = 0;
4982                while( 1 )
4983                {
4984                        fgets_IGNORE_RESULT( line, linelen-1, fp );
4985                        if( feof( fp ) ) break;
4986                        if( line[strlen(line)-1] != '\n' )
4987                        {
4988                                fprintf( stderr, "too long line? \n" );
4989                                exit( 1 );
4990                        }
4991                        if( line[0] == '#' ) continue;
4992                        if( atoi( line ) == 0 ) continue;
4993                        nmem = countspace( line );
4994                        if( nmem > *maxmempt ) *maxmempt = nmem;
4995                        (*nsubpt)++;
4996                }
4997        }
4998        else
4999        {
5000                tab01 = calloc( nseq, sizeof( int ) );
5001                for( i=0; i<nseq; i++ ) tab01[i] = 0;
5002                lpos = 0;
5003                while( 1 )
5004                {
5005                        fgets_IGNORE_RESULT( line, linelen-1, fp );
5006                        if( feof( fp ) ) break;
5007                        if( line[strlen(line)-1] != '\n' )
5008                        {
5009                                fprintf( stderr, "too long line? \n" );
5010                                exit( 1 );
5011                        }
5012                        if( line[0] == '#' ) continue;
5013                        if( atoi( line ) == 0 ) continue;
5014                        readasubalignment( line, table[lpos], preservegaps+lpos );
5015                        for( i=0; (p=table[lpos][i])!=-1; i++ )
5016                        {
5017                                if( tab01[p] )
5018                                {
5019                                        fprintf( stderr, "\nSequence %d appears in different groups.\n", p+1 );
5020                                        fprintf( stderr, "Hierarchical grouping is not supported.\n\n" );
5021                                        exit( 1 );
5022                                }
5023                                tab01[p] = 1;
5024                                if( p > nseq-1 )
5025                                {
5026                                        fprintf( stderr, "Sequence %d does not exist in the input sequence file.\n", p+1 );
5027                                        exit( 1 );
5028                                }
5029                        }
5030                        lpos++;
5031                }
5032                free( tab01 );
5033        }
5034        fclose( fp );
5035        free( line );
5036}
5037
5038
5039void readmccaskill( FILE *fp, RNApair **pairprob, int length )
5040{
5041        char gett[1000];
5042        int *pairnum;
5043        int i;
5044        int left, right;
5045        float prob;
5046        int c;
5047
5048        pairnum = (int *)calloc( length, sizeof( int ) );
5049        for( i=0; i<length; i++ ) pairnum[i] = 0;
5050
5051        c = getc( fp );
5052        {
5053                if( c != '>' )
5054                {
5055                        fprintf( stderr, "format error in hat4 - 1\n" );
5056                        exit( 1 );
5057                }
5058        }
5059        fgets_IGNORE_RESULT( gett, 999, fp );
5060        while( 1 )
5061        {
5062                if( feof( fp ) ) break;
5063                c = getc( fp );
5064                ungetc( c, fp );
5065                if( c == '>' || c == EOF )
5066                {
5067                        break;
5068                }
5069                fgets_IGNORE_RESULT( gett, 999, fp );
5070//              fprintf( stderr, "gett = %s\n", gett );
5071                sscanf( gett, "%d %d %f", &left, &right, &prob );
5072
5073                if( left >= length || right >= length )
5074                {
5075                        fprintf( stderr, "format error in hat4 - 2\n" );
5076                        exit( 1 );
5077                }
5078
5079                if( prob < 0.01 ) continue; // 080607, mafft ni dake eikyou
5080
5081                if( left != right && prob > 0.0 )
5082                {
5083                        pairprob[left] = (RNApair *)realloc( pairprob[left], (pairnum[left]+2) * sizeof( RNApair ) );
5084                        pairprob[left][pairnum[left]].bestscore = prob;
5085                        pairprob[left][pairnum[left]].bestpos = right;
5086                        pairnum[left]++;
5087                        pairprob[left][pairnum[left]].bestscore = -1.0;
5088                        pairprob[left][pairnum[left]].bestpos = -1;
5089//                      fprintf( stderr, "%d-%d, %f\n", left, right, prob );
5090
5091                        pairprob[right] = (RNApair *)realloc( pairprob[right], (pairnum[right]+2) * sizeof( RNApair ) );
5092                        pairprob[right][pairnum[right]].bestscore = prob;
5093                        pairprob[right][pairnum[right]].bestpos = left;
5094                        pairnum[right]++;
5095                        pairprob[right][pairnum[right]].bestscore = -1.0;
5096                        pairprob[right][pairnum[right]].bestpos = -1;
5097//                      fprintf( stderr, "%d-%d, %f\n", right, left, prob );
5098                }
5099        }
5100        free( pairnum );
5101}
5102
5103void readpairfoldalign( FILE *fp, char *s1, char *s2, char *aln1, char *aln2, int q1, int q2, int *of1, int *of2, int sumlen )
5104{
5105        char gett[1000];
5106        int *maptoseq1;
5107        int *maptoseq2;
5108        char dumc;
5109        int dumi;
5110        char sinseq[100], sinaln[100];
5111        int posinseq, posinaln;
5112        int alnlen;
5113        int i;
5114        int pos1, pos2;
5115        char *pa1, *pa2;
5116        char qstr[1000];
5117
5118        *of1 = -1;
5119        *of2 = -1;
5120
5121        maptoseq1 = AllocateIntVec( sumlen+1 );
5122        maptoseq2 = AllocateIntVec( sumlen+1 );
5123
5124        posinaln = 0; // foldalign ga alingment wo kaesanaitok no tame.
5125
5126        while( !feof( fp ) )
5127        {
5128                fgets_IGNORE_RESULT( gett, 999, fp );
5129                if( !strncmp( gett, "; ALIGNING", 10 ) ) break;
5130        }
5131        sprintf( qstr, "; ALIGNING            %d against %d\n", q1+1, q2+1 );
5132        if( strcmp( gett, qstr ) )
5133        {
5134                fprintf( stderr, "Error in FOLDALIGN\n" );
5135                fprintf( stderr, "qstr = %s, but gett = %s\n", qstr, gett );
5136                exit( 1 );
5137        }
5138
5139        while( !feof( fp ) )
5140        {
5141                fgets_IGNORE_RESULT( gett, 999, fp );
5142                if( !strncmp( gett, "; --------", 10 ) ) break;
5143        }
5144
5145
5146        while( !feof( fp ) )
5147        {
5148                fgets_IGNORE_RESULT( gett, 999, fp );
5149                if( !strncmp( gett, "; ********", 10 ) ) break;
5150//              fprintf( stderr, "gett = %s\n", gett );
5151                sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
5152                posinaln = atoi( sinaln );
5153                posinseq = atoi( sinseq );
5154//              fprintf( stderr, "posinseq = %d\n", posinseq );
5155//              fprintf( stderr, "posinaln = %d\n", posinaln );
5156                maptoseq1[posinaln-1] = posinseq-1;
5157        }
5158        alnlen = posinaln;
5159
5160        while( !feof( fp ) )
5161        {
5162                fgets_IGNORE_RESULT( gett, 999, fp );
5163                if( !strncmp( gett, "; --------", 10 ) ) break;
5164        }
5165
5166        while( !feof( fp ) )
5167        {
5168                fgets_IGNORE_RESULT( gett, 999, fp );
5169                if( !strncmp( gett, "; ********", 10 ) ) break;
5170//              fprintf( stderr, "gett = %s\n", gett );
5171                sscanf( gett, "%c %c %s %s %d %d", &dumc, &dumc, sinseq, sinaln, &dumi, &dumi );
5172                posinaln = atof( sinaln );
5173                posinseq = atof( sinseq );
5174//              fprintf( stderr, "posinseq = %d\n", posinseq );
5175//              fprintf( stderr, "posinaln = %d\n", posinaln );
5176                maptoseq2[posinaln-1] = posinseq-1;
5177        }
5178        if( alnlen != posinaln )
5179        {
5180                fprintf( stderr, "Error in foldalign?\n" );
5181                exit( 1 );
5182        }
5183
5184        pa1 = aln1;
5185        pa2 = aln2;
5186        for( i=0; i<alnlen; i++ )
5187        {
5188                pos1 = maptoseq1[i];
5189                pos2 = maptoseq2[i];
5190
5191                if( pos1 > -1 )
5192                        *pa1++ = s1[pos1];
5193                else
5194                        *pa1++ = '-';
5195
5196                if( pos2 > -1 )
5197                        *pa2++ = s2[pos2];
5198                else
5199                        *pa2++ = '-';
5200        }
5201        *pa1 = 0;
5202        *pa2 = 0;
5203
5204        *of1 = 0;
5205        for( i=0; i<alnlen; i++ )
5206        {
5207                *of1 = maptoseq1[i];
5208                if( *of1 > -1 ) break;
5209        }
5210        *of2 = 0;
5211        for( i=0; i<alnlen; i++ )
5212        {
5213                *of2 = maptoseq2[i];
5214                if( *of2 > -1 ) break;
5215        }
5216
5217//      fprintf( stderr, "*of1=%d, aln1 = :%s:\n", *of1, aln1 );
5218//      fprintf( stderr, "*of2=%d, aln2 = :%s:\n", *of2, aln2 );
5219
5220        free( maptoseq1 );
5221        free( maptoseq2 );
5222}
Note: See TracBrowser for help on using the repository browser.