source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/addsingle.c

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)

File size: 76.3 KB
Line 
1#include "mltaln.h"
2
3#define SMALLMEMORY 1
4
5#define DEBUG 0
6#define IODEBUG 0
7#define SCOREOUT 0
8
9static int nadd;
10static int treein;
11static int topin;
12static int treeout;
13static int distout;
14static int noalign;
15static int multidist;
16static int maxdist = 1;
17
18static float lenfaca, lenfacb, lenfacc, lenfacd;
19static int tuplesize;
20
21#define PLENFACA 0.01
22#define PLENFACB 10000
23#define PLENFACC 10000
24#define PLENFACD 0.1
25#define D6LENFACA 0.01
26#define D6LENFACB 2500
27#define D6LENFACC 2500
28#define D6LENFACD 0.1
29#define D10LENFACA 0.01
30#define D10LENFACB 1000000
31#define D10LENFACC 1000000
32#define D10LENFACD 0.0
33
34typedef struct _thread_arg
35{
36        int njob; 
37        int nadd; 
38        int *nlen; 
39        int *follows; 
40        char **name; 
41        char **seq; 
42        LocalHom **localhomtable; 
43        float **iscore; 
44        float **nscore; 
45        int *istherenewgap; 
46        int **newgaplist;
47        RNApair ***singlerna; 
48        double *eff_kozo_mapped; 
49        int alloclen;
50        Treedep *dep;
51        int  ***topol;
52        float  **len;
53        Addtree *addtree;
54#ifdef enablemultithread
55        int *iaddshare;
56        int thread_no;
57        pthread_mutex_t *mutex_counter;
58#endif
59} thread_arg_t;
60
61
62#ifdef enablemultithread
63typedef struct _gaplist2alnxthread_arg
64{
65//      int thread_no;
66        int ncycle;
67        int *jobpospt;
68        int tmpseqlen;
69        int lenfull;
70        char **seq;
71        int *newgaplist;
72        int *posmap;
73        pthread_mutex_t *mutex;
74} gaplist2alnxthread_arg_t;
75
76typedef struct _distancematrixthread_arg
77{
78        int thread_no;
79        int njob;
80        int norg;
81        int *jobpospt;
82        int **pointt;
83        int *nogaplen;
84        float **imtx;
85        float **nmtx;
86        float *selfscore;
87        pthread_mutex_t *mutex;
88} distancematrixthread_arg_t;
89
90typedef struct _jobtable2d
91{
92    int i; 
93    int j; 
94} Jobtable2d;
95
96typedef struct _dndprethread_arg
97{
98        int njob;
99        int thread_no;
100        float *selfscore;
101        float **mtx;
102        char **seq;
103        Jobtable2d *jobpospt;
104        pthread_mutex_t *mutex;
105} dndprethread_arg_t;
106
107#endif
108
109typedef struct _blocktorealign
110{
111        int start;
112        int end;
113        int nnewres;
114} Blocktorealign;
115
116static void cnctintvec( int *res, int *o1, int *o2 )
117{
118        while( *o1 != -1 ) *res++ = *o1++;
119        while( *o2 != -1 ) *res++ = *o2++;
120        *res = -1;
121}
122
123static void countnewres( int len, Blocktorealign *realign, int *posmap, int *gaplist )
124{
125        int i, regstart, regend, len1;
126        regstart = 0;
127        len1 = len+1;
128        for( i=0; i<len1; i++ )
129        {
130                if( realign[i].nnewres || gaplist[i] )
131                {
132                        regend = posmap[i]-1;
133                        realign[i].start = regstart;
134                        realign[i].end = regend;
135                }
136                if( gaplist[i] )
137                {
138                        realign[i].nnewres++;
139//                      fprintf( stderr, "hit? reg = %d-%d\n", regstart, regend );
140                }
141                regstart = posmap[i]+1;
142        }
143}
144static void fillgap( char *s, int len )
145{
146        int orilen = strlen( s );
147        s += orilen;
148        len -= orilen;
149        while( len-- )
150                *s++ = '-';
151        *s = 0;
152}
153
154static int lencomp( const void *a, const void *b ) // osoikamo
155{
156        char **ast = (char **)a;
157        char **bst = (char **)b;
158        int lena = strlen( *ast );
159        int lenb = strlen( *bst );
160//      fprintf( stderr, "a=%s, b=%s\n", *ast, *bst );
161//      fprintf( stderr, "lena=%d, lenb=%d\n", lena, lenb );
162        if( lena > lenb ) return -1;
163        else if( lena < lenb ) return 1;
164        else return( 0 );
165}
166
167static int dorealignment_tree( Blocktorealign *block, char **fullseq, int *fullseqlenpt, int norg, int ***topol, int *follows )
168{
169        int i, j, k, posinold, newlen, *nmem;
170        int n0, n1, localloclen, nhit, hit1, hit2;
171        int *pickhistory;
172        int nprof1, nprof2, pos, zure;
173        char **prof1, **prof2;
174        int *iinf0, *iinf1;
175        int *group, *nearest, *g2n, ngroup;
176        char ***mem;
177        static char **tmpaln0 = NULL;
178        static char **tmpaln1 = NULL;
179        static char **tmpseq;
180        int ***topolpick;
181        int *tmpint;
182        int *intptr, *intptrx;
183        char *tmpseq0, *cptr, **cptrptr;
184
185
186        localloclen = 4 * ( block->end - block->start + 1 );     // ookisugi?
187        tmpaln0 = AllocateCharMtx( njob, localloclen );
188        tmpaln1 = AllocateCharMtx( njob, localloclen );
189        tmpseq = AllocateCharMtx( 1, *fullseqlenpt * 4 );
190        iinf0 = AllocateIntVec( njob );
191        iinf1 = AllocateIntVec( njob );
192        nearest = AllocateIntVec( njob ); // oosugi
193
194        posinold = block->start;
195
196        n0 = 0;
197        n1 = 0;
198        for( i=0; i<njob; i++ )
199        {
200                strncpy( tmpseq[0], fullseq[i] + block->start, block->end - block->start + 1 );
201                tmpseq[0][block->end - block->start + 1] = 0;
202                commongappick( 1, tmpseq );
203                if( tmpseq[0][0] != 0 )
204                {
205                        if( i < norg )
206                        {
207                                fprintf( stderr, "BUG!!!!\n" );
208                                exit( 1 );
209                        }
210                        strcpy( tmpaln0[n0], tmpseq[0] );
211                        iinf0[n0] = i;
212                        nearest[n0] = follows[i-norg];
213                        n0++;
214                }
215                else
216                {
217                        strcpy( tmpaln1[n0], "" );
218                        iinf1[n1] = i;
219                        n1++;
220                }
221        }
222        mem = AllocateCharCub( n0, n0+1, 0 ); // oosugi
223        nmem = AllocateIntVec( n0 ); // oosugi
224        g2n = AllocateIntVec( n0 ); // oosugi
225        group = AllocateIntVec( n0 ); // oosugi
226        for( i=0; i<n0; i++ ) mem[i][0] = NULL;
227        for( i=0; i<n0; i++ ) nmem[i] = 0;
228        ngroup = 0;
229        for( i=0; i<n0; i++ )
230        {
231                for( j=0; j<i; j++ ) if( nearest[j] == nearest[i] ) break;
232                if( j == i ) group[i] = ngroup++;
233                else group[i] = group[j];
234
235                for( j=0; mem[group[i]][j]; j++ )
236                        ;
237                mem[group[i]][j] = tmpaln0[i];
238                mem[group[i]][j+1] = NULL;
239                nmem[group[i]]++;
240                g2n[group[i]] = nearest[i];
241//              fprintf( stderr, "%d -> %d -> group%d\n", i, nearest[i], group[i] );
242//              fprintf( stderr, "mem[%d][%d] = %s\n", group[i], j, mem[group[i]][j] );
243        }
244
245        for( i=0; i<ngroup; i++ )
246        {
247//              fprintf( stderr, "before sort:\n" );
248//              for( j=0; j<nmem[i]; j++ ) fprintf( stderr, "%s\n", mem[i][j] );
249//              fprintf( stderr, "\n" );
250                qsort( mem[i], nmem[i], sizeof( char * ), lencomp );
251//              fprintf( stderr, "after sort:\n" );
252//              for( j=0; j<nmem[i]; j++ ) fprintf( stderr, "%s\n", mem[i][j] );
253//              fprintf( stderr, "\n" );
254        }
255
256#if 0
257        for( i=1; i<n0; i++ )
258        {
259                profilealignment( 1, n1, i, tmpaln0+i, tmpaln1, tmpaln0, localloclen, alg );
260        }
261        newlen = strlen( tmpaln0[0] );
262        for( i=0; i<n1; i++ ) eq2dash( tmpaln1[i] );
263#else
264//      newlen = 0;
265        for( i=0; i<ngroup; i++ )
266        {
267//              for( k=0; mem[i][k]; k++ ) fprintf( stderr, "mem[%d][%d] = %s\n", i, j, mem[i][k] );
268
269                for( j=1; j<nmem[i]; j++ )
270                {
271                        profilealignment2( 1, j, mem[i]+j, mem[i], localloclen, alg );
272                }
273//              for( j=0; j<nmem[i]; j++ ) fprintf( stderr, "j=%d, %s\n", j, mem[i][j] );
274
275#if 0 // iru
276                if( ( j = strlen( mem[i][0] ) ) > newlen ) newlen = j;
277                for( j=0; j<=i; j++ )
278                {
279                        for( k=0; mem[j][k]; k++ )
280                                fillgap( mem[j][k], newlen );
281                }
282#endif
283
284        }
285#if 0
286        fprintf( stderr, "After ingroupalignment (original order):\n" );
287        for( i=0; i<n0; i++ ) fprintf( stderr, "%s\n", tmpaln0[i] );
288#endif
289#endif
290
291        topolpick = AllocateIntCub( ngroup, 2, ngroup );
292        pickhistory = AllocateIntVec( ngroup );
293        tmpint = AllocateIntVec( 2 );
294        prof1 = AllocateCharMtx( n0, 0 );
295        prof2 = AllocateCharMtx( n0, 0 );
296        for( i=0; i<ngroup; i++ )
297        {
298                topolpick[i][0][0] = -1;
299                topolpick[i][1][0] = -1;
300                pickhistory[i] = -1;
301        }
302
303        nhit = 0;
304        for( i=0; i<norg-1; i++ )
305        {
306                for( intptr=topol[i][0]; *intptr>-1; intptr++ )
307                {
308                        for( intptrx=g2n,k=0; k<ngroup; intptrx++,k++ ) 
309                        {
310                                if( *intptr == *intptrx )
311                                {
312                                        hit1 = k;
313                                        goto exitloop1;
314                                }
315                        }
316                }
317                continue; 
318                exitloop1:
319//              fprintf( stderr, "hit1! group%d -> %d\n", k, topol[i][0][j] );
320
321                for( intptr=topol[i][1]; *intptr>-1; intptr++ )
322                {
323                        for( intptrx=g2n,k=0; k<ngroup; intptrx++,k++ ) 
324                        {
325                                if( *intptr == *intptrx )
326                                {
327                                        hit2 = k;
328                                        goto exitloop2;
329                                }
330                        }
331                }
332                continue; 
333                exitloop2:
334
335                if( pickhistory[hit1] == -1 )
336                {
337                        topolpick[nhit][0][0] = hit1;
338                        topolpick[nhit][0][1] = -1;
339                }
340                else
341                {
342                        intcpy( topolpick[nhit][0], topolpick[pickhistory[hit1]][0] );
343                        intcat( topolpick[nhit][0], topolpick[pickhistory[hit1]][1] );
344                }
345                if( pickhistory[hit2] == -1 )
346                {
347                        topolpick[nhit][1][0] = hit2;
348                        topolpick[nhit][1][1] = -1;
349                }
350                else
351                {
352                        intcpy( topolpick[nhit][1], topolpick[pickhistory[hit2]][0] );
353                        intcat( topolpick[nhit][1], topolpick[pickhistory[hit2]][1] );
354                }
355
356                pickhistory[hit1] = nhit;
357                pickhistory[hit2] = nhit;
358                nhit++;
359//              g2n[hit1] = -1;
360//              g2n[hit2] = -1;
361
362//              fprintf( stderr, "hit2! group%d -> %d\n", k, topol[i][1][j] );
363
364#if 0
365                fprintf( stderr, "\nHIT!!! \n" );
366                fprintf( stderr, "\nSTEP %d\n", i );
367                for( j=0; topol[i][0][j]>-1; j++ ) fprintf( stderr, "%3d ", topol[i][0][j] );
368                fprintf( stderr, "\n" );
369                for( j=0; topol[i][1][j]>-1; j++ ) fprintf( stderr, "%3d ", topol[i][1][j] );
370                fprintf( stderr, "\n" );
371#endif
372        }
373
374        for( i=0; i<ngroup-1; i++ )
375        {
376#if 0
377                fprintf( stderr, "\npickSTEP %d\n", i );
378                for( j=0; topolpick[i][0][j]>-1; j++ ) fprintf( stderr, "%3d ", topolpick[i][0][j] );
379                fprintf( stderr, "\n" );
380                for( j=0; topolpick[i][1][j]>-1; j++ ) fprintf( stderr, "%3d ", topolpick[i][1][j] );
381                fprintf( stderr, "\n" );
382#endif
383
384                pos = 0;
385//              for( j=0; topolpick[i][0][j]>-1; j++ ) for( k=0; (cptr=mem[topolpick[i][0][j]][k]); k++ ) prof1[pos++] = cptr;
386                for( intptr=topolpick[i][0]; *intptr>-1; intptr++ ) 
387                        for( cptrptr=mem[*intptr]; (cptr=*cptrptr); cptrptr++ ) 
388                                prof1[pos++] = cptr;
389                nprof1 = pos;
390                pos = 0;
391//              for( j=0; topolpick[i][1][j]>-1; j++ ) for( k=0; (cptr=mem[topolpick[i][1][j]][k]); k++ ) prof2[pos++] = cptr;
392                for( intptr=topolpick[i][1]; *intptr>-1; intptr++ ) 
393                        for( cptrptr=mem[*intptr]; (cptr=*cptrptr); cptrptr++ ) 
394                                prof2[pos++] = cptr;
395                nprof2 = pos;
396
397
398                profilealignment2( nprof1, nprof2, prof1, prof2, localloclen, alg );
399#if 0
400                for( j=0; j<nprof1; j++ ) fprintf( stderr, "prof1[%d] = %s\n", j, prof1[j] );
401                for( j=0; j<nprof2; j++ ) fprintf( stderr, "prof2[%d] = %s\n", j, prof2[j] );
402                fprintf( stderr, "done.\n" );
403#endif
404        }
405        newlen = strlen( tmpaln0[0] );
406        for( j=0; j<n1; j++ ) fillgap( tmpaln1[j], newlen );
407
408#if 0
409        fprintf( stderr, "After rerealignment (original order):\n" );
410        for( i=0; i<n0; i++ ) fprintf( stderr, "%s\n", tmpaln0[i] );
411#endif
412
413//      newlen = strlen( tmpaln0[0] );
414        zure = ( block->end - block->start + 1 - newlen );
415//      fprintf( stderr, "zure = %d, localloclen=%d, newlen=%d\n", zure, localloclen, newlen );
416
417
418        if( *fullseqlenpt < strlen( fullseq[0] ) - (block->end-block->start+1)  + newlen + 1 )
419        {
420                *fullseqlenpt = strlen( fullseq[0] ) * 2;
421                fprintf( stderr, "reallocating..." );
422                for( i=0; i<njob; i++ )
423                {
424                        fullseq[i] = realloc( fullseq[i], *fullseqlenpt * sizeof( char ) );
425                        if( !fullseq[i] )
426                        {
427                                fprintf( stderr, "Cannot reallocate seq[][]\n" );
428                                exit( 1 );
429                        }
430                }
431                fprintf( stderr, "done.\n" );
432        }
433
434
435        tmpseq0 = tmpseq[0];
436        posinold = block->end+1;
437        for( i=0; i<n0; i++ ) 
438        {
439                strncpy( tmpseq0, tmpaln0[i], newlen );
440                strcpy( tmpseq0+newlen, fullseq[iinf0[i]] + posinold );
441                strcpy( fullseq[iinf0[i]]+block->start, tmpseq0 );
442        }
443        for( i=0; i<n1; i++ ) 
444        {
445//              eq2dash( tmpaln1[i] );
446                strncpy( tmpseq0, tmpaln1[i], newlen );
447                strcpy( tmpseq0+newlen, fullseq[iinf1[i]] + posinold );
448                strcpy( fullseq[iinf1[i]]+block->start, tmpseq0 );
449        }
450        FreeCharMtx( tmpaln0 );
451        FreeCharMtx( tmpaln1 );
452        FreeCharMtx( tmpseq );
453        for( i=0; i<n0; i++ ) 
454        {
455//              for( j=0; j<njob; j++ ) free( mem[i][j] );
456                free( mem[i] );
457        }
458        free( mem );
459        free( nmem );
460        free( iinf0 );
461        free( iinf1 );
462        free( group );
463        free( g2n );
464        free( prof1 );
465        free( prof2 );
466        free( nearest );
467        FreeIntCub( topolpick );
468        free( pickhistory );
469        free( tmpint );
470
471        return( zure );
472}
473
474
475#if 0
476static int dorealignment( Blocktorealign *block, char **fullseq, int alloclen, int fullseqlen, int norg )
477{
478        int i, posinnew, posinold, newlen;
479        int n0, n1;
480        int zure;
481        static int *iinf0, *iinf1;
482        static char **tmpaln0 = NULL;
483        static char **tmpaln1 = NULL;
484        static char **tmpseq;
485        char *opt, *npt;
486
487        if( tmpaln0 == NULL )
488        {
489                tmpaln0 = AllocateCharMtx( njob, alloclen );
490                tmpaln1 = AllocateCharMtx( njob, alloclen );
491                tmpseq = AllocateCharMtx( 1, fullseqlen );
492                iinf0 = AllocateIntVec( njob );
493                iinf1 = AllocateIntVec( njob );
494        }
495        posinold = block->start;
496
497
498        n0 = 0;
499        n1 = 0;
500        for( i=0; i<njob; i++ )
501        {
502                strncpy( tmpseq[0], fullseq[i] + block->start, block->end - block->start + 1 );
503                tmpseq[0][block->end - block->start + 1] = 0;
504                commongappick( 1, tmpseq );
505//              if( strlen( tmpseq[0] ) > 0 )
506                if( tmpseq[0][0] != 0 )
507                {
508                        if( i < norg )
509                        {
510                                fprintf( stderr, "BUG!!!!\n" );
511                                exit( 1 );
512                        }
513                        strcpy( tmpaln0[n0], tmpseq[0] );
514                        iinf0[n0] = i;
515                        n0++;
516                }
517                else
518                {
519                        strcpy( tmpaln1[n0], "" );
520                        iinf1[n1] = i;
521                        n1++;
522                }
523        }
524
525
526        for( i=1; i<n0; i++ )
527        {
528                profilealignment( 1, n1, i, tmpaln0+i, tmpaln1, tmpaln0, alloclen, alg ); // n1 ha allgap
529        }
530
531#if 1
532        fprintf( stderr, "After realignment:\n" );
533        for( i=0; i<n0; i++ ) fprintf( stderr, "%s\n", tmpaln0[i] );
534//      for( i=0; i<n1; i++ ) fprintf( stderr, "%s\n", tmpaln1[i] );
535#endif
536
537        newlen = strlen( tmpaln0[0] );
538        for( i=0; i<n0; i++ ) strncpy( fullseq[iinf0[i]]+block->start, tmpaln0[i], newlen );
539        for( i=0; i<n1; i++ )
540        {
541                eq2dash( tmpaln1[i] );
542                strncpy( fullseq[iinf1[i]] + block->start, tmpaln1[i], newlen );
543        }
544
545        posinold = block->end+1;
546        posinnew = block->start + newlen;
547
548
549        zure = ( block->end - block->start + 1 - strlen( tmpaln0[0] ) );
550
551        for( i=0; i<njob; i++ )
552        {
553#if 0
554                strcpy( fullseq[i]+posinnew, fullseq[i]+posinold ); // ??
555#else
556                opt = fullseq[i] + posinold;
557                npt = fullseq[i] + posinnew;
558                while( ( *npt++ = *opt++ ) );
559                *npt = 0;
560#endif
561        }
562
563        return( zure );
564}
565#endif
566
567static void adjustposmap( int len, int *posmap, int *gaplist )
568{
569        int *newposmap;
570        int *mpt1, *mpt2;
571        int lenbk, zure;
572        newposmap = calloc( len+2, sizeof( int ) );
573        lenbk = len;
574        zure = 0;
575        mpt1 = newposmap;
576        mpt2 = posmap;
577
578#if 0
579        int i;
580        fprintf( stderr, "posmapa = " );
581        for( i=0; i<len+2; i++ )
582        {
583                fprintf( stderr, "%3d ", posmap[i] );
584        }
585        fprintf( stderr, "\n" );
586#endif
587
588        while( len-- ) 
589        {
590                zure += *gaplist++;
591                *mpt1++ = *mpt2++ + zure;
592        }
593        zure += *gaplist++;
594        *mpt1 = *mpt2 + zure;
595
596        mpt1 = newposmap;
597        mpt2 = posmap;
598        len = lenbk;
599        while( len-- ) *mpt2++ = *mpt1++;
600        *mpt2 = *mpt1;
601        free( newposmap );
602#if 0
603        fprintf( stderr, "posmapm = " );
604        for( i=0; i<lenbk+2; i++ )
605        {
606                fprintf( stderr, "%3d ", posmap[i] );
607        }
608        fprintf( stderr, "\n" );
609#endif
610}
611
612static int insertgapsbyotherfragments_compact( int len, char *a, char *s, int *l, int *p )
613{
614        int gaplen;
615        int i, pos, pi;
616        int prevp = -1;
617        int realignment = 0;
618//      fprintf( stderr, "### insertgapsbyotherfragments\n" );
619        for( i=0; i<len; i++ )
620        {
621                gaplen = l[i];
622                pi = p[i];
623                pos = prevp + 1;
624//              fprintf( stderr, "gaplen = %d\n", gaplen );
625                while( gaplen-- )
626                {
627                        pos++;
628                        *a++ = *s++;
629                }
630//              fprintf( stderr, "pos = %d, pi = %d\n", pos, pi );
631                while( pos++ < pi )
632                {
633                        *a++ = '=';
634                        realignment = 1;
635                }
636                *a++ = *s++;
637                prevp = pi;
638        }
639        gaplen = l[i];
640        pi = p[i];
641        pos = prevp + 1;
642        while( gaplen-- )
643        {
644                pos++;
645                *a++ = *s++;
646        }
647        while( pos++ < pi )
648        {
649                *a++ = '=';
650                realignment = 1;
651        }
652        *a = 0;
653        return( realignment );
654}
655
656void makegaplistcompact( int len, int *p, int *c, int *l )
657{
658        int i;
659        int pg;
660        int prep = -1;
661        for( i=0; i<len+2; i++ )
662        {
663                if( ( pg = p[i]-prep-1) > 0 && l[i] > 0 )
664                {
665                        if( pg < l[i] ) 
666                        {
667                                c[i] = l[i] - pg;
668                        }
669                        else
670                        {
671                                c[i] = 0;
672                        }
673                }
674                else
675                {
676                        c[i] = l[i];
677                }
678                prep = p[i];
679        }
680}
681
682
683void gaplist2alnx( int len, char *a, char *s, int *l, int *p, int lenlimit )
684{
685        int gaplen;
686        int pos, pi, posl;
687        int prevp = -1;
688        int reslen = 0;
689        char *sp;
690//      char *abk = a;
691
692#if 0
693        int i;
694        char *abk = a;
695        fprintf( stderr, "s = %s\n", s );
696        fprintf( stderr, "posmap  = " );
697        for( i=0; i<len+2; i++ )
698        {
699                fprintf( stderr, "%3d ", p[i] );
700        }
701        fprintf( stderr, "\n" );
702        fprintf( stderr, "gaplist = " );
703        for( i=0; i<len+2; i++ )
704        {
705                fprintf( stderr, "%3d ", l[i] );
706        }
707        fprintf( stderr, "\n" );
708#endif
709        while( len-- )
710        {
711                gaplen = *l++;
712                pi = *p++;
713
714                if( (reslen+=gaplen) > lenlimit )
715                {
716                        fprintf( stderr, "Length over. Please recompile!\n" );
717                        exit( 1 );
718                }
719                while( gaplen-- ) *a++ = '-';
720
721                pos = prevp + 1;
722                sp = s + pos;
723                if( ( posl = pi - pos ) )
724                {
725                        if( ( reslen += posl ) > lenlimit )
726                        {
727                                fprintf( stderr, "Length over. Please recompile\n" );
728                                exit( 1 );
729                        }
730                        while( posl-- ) *a++ = *sp++;
731                }
732
733                if( reslen++ > lenlimit )
734                {
735                        fprintf( stderr, "Length over. Please recompile\n" );
736                        exit( 1 );
737                }
738                *a++ = *sp;
739                prevp = pi;
740        }
741
742        gaplen = *l;
743        pi = *p;
744        if( (reslen+=gaplen) > lenlimit )
745        {
746                fprintf( stderr, "Length over. Please recompile\n" );
747                exit( 1 );
748        }
749        while( gaplen-- ) *a++ = '-';
750
751        pos = prevp + 1;
752        sp = s + pos;
753        if( ( posl = pi - pos ) )
754        {
755                if( ( reslen += posl ) > lenlimit )
756                {
757                        fprintf( stderr, "Length over. Please recompile\n" );
758                        exit( 1 );
759                }
760                while( posl-- ) *a++ = *sp++;
761        }
762        *a = 0;
763//      fprintf( stderr, "reslen = %d, strlen(a) = %d\n", reslen, strlen( abk ) );
764//      fprintf( stderr, "a = %s\n", abk );
765}
766
767static void makenewgaplist( int *l, char *a )
768{
769        while( 1 )
770        {
771                while( *a == '=' )
772                {
773                        a++;
774                        (*l)++;
775//                      fprintf( stderr, "a[] (i) = %s, *l=%d\n", a, *(l) );
776                }
777                *++l = 0;
778                if( *a == 0 ) break;
779                a++;
780        }
781        *l = -1;
782}
783
784
785void arguments( int argc, char *argv[] )
786{
787    int c;
788
789        nthread = 1;
790        outnumber = 0;
791        scoreout = 0;
792        treein = 0;
793        topin = 0;
794        rnaprediction = 'm';
795        rnakozo = 0;
796        nevermemsave = 0;
797        inputfile = NULL;
798        addfile = NULL;
799        addprofile = 1;
800        fftkeika = 0;
801        constraint = 0;
802        nblosum = 62;
803        fmodel = 0;
804        calledByXced = 0;
805        devide = 0;
806        use_fft = 0; // chuui
807        force_fft = 0;
808        fftscore = 1;
809        fftRepeatStop = 0;
810        fftNoAnchStop = 0;
811    weight = 3;
812    utree = 1;
813        tbutree = 1;
814    refine = 0;
815    check = 1;
816    cut = 0.0;
817    disp = 0;
818    outgap = 1;
819    alg = 'A';
820    mix = 0;
821        tbitr = 0;
822        scmtd = 5;
823        tbweight = 0;
824        tbrweight = 3;
825        checkC = 0;
826        treemethod = 'X';
827        contin = 0;
828        scoremtx = 1;
829        kobetsubunkatsu = 0;
830        dorp = NOTSPECIFIED;
831        ppenalty = NOTSPECIFIED;
832        ppenalty_ex = NOTSPECIFIED;
833        poffset = NOTSPECIFIED;
834        kimuraR = NOTSPECIFIED;
835        pamN = NOTSPECIFIED;
836        geta2 = GETA2;
837        fftWinSize = NOTSPECIFIED;
838        fftThreshold = NOTSPECIFIED;
839        RNAppenalty = NOTSPECIFIED;
840        RNAppenalty_ex = NOTSPECIFIED;
841        RNApthr = NOTSPECIFIED;
842        TMorJTT = JTT;
843        consweight_multi = 1.0;
844        consweight_rna = 0.0;
845        nadd = 0;
846        multidist = 0;
847        tuplesize = -1;
848
849    while( --argc > 0 && (*++argv)[0] == '-' )
850        {
851        while ( ( c = *++argv[0] ) )
852                {
853            switch( c )
854            {
855                                case 'i':
856                                        inputfile = *++argv;
857                                        fprintf( stderr, "inputfile = %s\n", inputfile );
858                                        --argc;
859                    goto nextoption;
860                                case 'I':
861                                        nadd = atoi( *++argv );
862                                        fprintf( stderr, "nadd = %d\n", nadd );
863                                        --argc;
864                                        goto nextoption;
865                                case 'e':
866                                        RNApthr = (int)( atof( *++argv ) * 1000 - 0.5 );
867                                        --argc;
868                                        goto nextoption;
869                                case 'o':
870                                        RNAppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
871                                        --argc;
872                                        goto nextoption;
873                                case 'f':
874                                        ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
875//                                      fprintf( stderr, "ppenalty = %d\n", ppenalty );
876                                        --argc;
877                                        goto nextoption;
878                                case 'g':
879                                        ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
880                                        fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
881                                        --argc;
882                                        goto nextoption;
883                                case 'h':
884                                        poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
885//                                      fprintf( stderr, "poffset = %d\n", poffset );
886                                        --argc;
887                                        goto nextoption;
888                                case 'k':
889                                        kimuraR = atoi( *++argv );
890                                        fprintf( stderr, "kappa = %d\n", kimuraR );
891                                        --argc;
892                                        goto nextoption;
893                                case 'b':
894                                        nblosum = atoi( *++argv );
895                                        scoremtx = 1;
896                                        fprintf( stderr, "blosum %d / kimura 200\n", nblosum );
897                                        --argc;
898                                        goto nextoption;
899                                case 'j':
900                                        pamN = atoi( *++argv );
901                                        scoremtx = 0;
902                                        TMorJTT = JTT;
903                                        fprintf( stderr, "jtt/kimura %d\n", pamN );
904                                        --argc;
905                                        goto nextoption;
906                                case 'm':
907                                        pamN = atoi( *++argv );
908                                        scoremtx = 0;
909                                        TMorJTT = TM;
910                                        fprintf( stderr, "tm %d\n", pamN );
911                                        --argc;
912                                        goto nextoption;
913                                case 'l':
914                                        fastathreshold = atof( *++argv );
915                                        constraint = 2;
916                                        --argc;
917                                        goto nextoption;
918                                case 'r':
919                                        consweight_rna = atof( *++argv );
920                                        rnakozo = 1;
921                                        --argc;
922                                        goto nextoption;
923                                case 'c':
924                                        consweight_multi = atof( *++argv );
925                                        --argc;
926                                        goto nextoption;
927                                case 'C':
928                                        nthread = atoi( *++argv );
929                                        fprintf( stderr, "nthread = %d\n", nthread );
930                                        --argc; 
931                                        goto nextoption;
932                                case 'R':
933                                        rnaprediction = 'r';
934                                        break;
935                                case 's':
936                                        RNAscoremtx = 'r';
937                                        break;
938#if 1
939                                case 'a':
940                                        fmodel = 1;
941                                        break;
942#endif
943                                case 'K':
944                                        addprofile = 0;
945                                        break;
946                                case 'y':
947                                        distout = 1;
948                                        break;
949                                case 't':
950                                        treeout = 1;
951                                        break;
952                                case 'T':
953                                        noalign = 1;
954                                        break;
955                                case 'D':
956                                        dorp = 'd';
957                                        break;
958                                case 'P':
959                                        dorp = 'p';
960                                        break;
961#if 1
962                                case 'O':
963                                        outgap = 0;
964                                        break;
965#else
966                                case 'O':
967                                        fftNoAnchStop = 1;
968                                        break;
969#endif
970                                case 'S':
971                                        scoreout = 1;
972                                        break;
973#if 0
974                                case 'e':
975                                        fftscore = 0;
976                                        break;
977                                case 'r':
978                                        fmodel = -1;
979                                        break;
980                                case 'R':
981                                        fftRepeatStop = 1;
982                                        break;
983                                case 's':
984                                        treemethod = 's';
985                                        break;
986#endif
987                                case 'X':
988                                        treemethod = 'X';
989                                        break;
990                                case 'E':
991                                        treemethod = 'E';
992                                        break;
993                                case 'q':
994                                        treemethod = 'q';
995                                        break;
996                                case 'n' :
997                                        outnumber = 1;
998                                        break;
999#if 0
1000                                case 'a':
1001                                        alg = 'a';
1002                                        break;
1003#endif
1004                                case 'Q':
1005                                        alg = 'Q';
1006                                        break;
1007                                case 'H':
1008                                        alg = 'H';
1009                                        break;
1010                                case 'A':
1011                                        alg = 'A';
1012                                        break;
1013                                case 'M':
1014                                        alg = 'M';
1015                                        break;
1016                                case 'N':
1017                                        nevermemsave = 1;
1018                                        break;
1019                                case 'B':
1020                                        break;
1021                                case 'F':
1022                                        use_fft = 1;
1023                                        break;
1024                                case 'G':
1025                                        force_fft = 1;
1026                                        use_fft = 1;
1027                                        break;
1028                                case 'U':
1029                                        treein = 1;
1030                                        break;
1031                                case 'V':
1032                                        topin = 1;
1033                                        break;
1034                                case 'u':
1035                                        tbrweight = 0;
1036                                        weight = 0;
1037                                        break;
1038                                case 'v':
1039                                        tbrweight = 3;
1040                                        break;
1041                                case 'd':
1042                                        multidist = 1;
1043                                        break;
1044                                case 'W':
1045                                        tuplesize = atoi( *++argv );
1046                                        --argc;
1047                                        goto nextoption;
1048#if 0
1049                                case 'd':
1050                                        disp = 1;
1051                                        break;
1052#endif
1053/* Modified 01/08/27, default: user tree */
1054                                case 'J':
1055                                        tbutree = 0;
1056                                        break;
1057/* modification end. */
1058                                case 'z':
1059                                        fftThreshold = atoi( *++argv );
1060                                        --argc; 
1061                                        goto nextoption;
1062                                case 'w':
1063                                        fftWinSize = atoi( *++argv );
1064                                        --argc;
1065                                        goto nextoption;
1066                                case 'Z':
1067                                        checkC = 1;
1068                                        break;
1069                default:
1070                    fprintf( stderr, "illegal option %c\n", c );
1071                    argc = 0;
1072                    break;
1073            }
1074                }
1075                nextoption:
1076                        ;
1077        }
1078    if( argc == 1 )
1079    {
1080        cut = atof( (*argv) );
1081        argc--;
1082    }
1083    if( argc != 0 ) 
1084    {
1085        fprintf( stderr, "options: Check source file !\n" );
1086        exit( 1 );
1087    }
1088        if( tbitr == 1 && outgap == 0 )
1089        {
1090                fprintf( stderr, "conflicting options : o, m or u\n" );
1091                exit( 1 );
1092        }
1093        if( alg == 'C' && outgap == 0 )
1094        {
1095                fprintf( stderr, "conflicting options : C, o\n" );
1096                exit( 1 );
1097        }
1098}
1099
1100
1101static float treebase( int nseq, int *nlen, char **aseq, int nadd, char *mergeoralign, char **mseq1, char **mseq2, int ***topol, double *effarr, int *alloclen, LocalHom **localhomtable, RNApair ***singlerna, double *effarr_kozo )
1102{
1103        int i, l, m;
1104        int len1nocommongap, len2nocommongap;
1105        int len1, len2;
1106        int clus1, clus2;
1107        float pscore, tscore;
1108        char *indication1, *indication2;
1109        double *effarr1 = NULL;
1110        double *effarr2 = NULL;
1111        double *effarr1_kozo = NULL;
1112        double *effarr2_kozo = NULL;
1113        LocalHom ***localhomshrink = NULL;
1114        int *fftlog;
1115        int m1, m2;
1116        int *gaplen;
1117        int *gapmap;
1118        int *alreadyaligned;
1119        float dumfl = 0.0;
1120        int ffttry;
1121        RNApair ***grouprna1, ***grouprna2;
1122
1123        if( rnakozo && rnaprediction == 'm' )
1124        {
1125                grouprna1 = (RNApair ***)calloc( nseq, sizeof( RNApair ** ) );
1126                grouprna2 = (RNApair ***)calloc( nseq, sizeof( RNApair ** ) );
1127        }
1128        else
1129        {
1130                grouprna1 = grouprna2 = NULL;
1131        }
1132
1133        fftlog = AllocateIntVec( nseq );
1134        effarr1 = AllocateDoubleVec( nseq );
1135        effarr2 = AllocateDoubleVec( nseq );
1136        indication1 = AllocateCharVec( 150 );
1137        indication2 = AllocateCharVec( 150 );
1138        alreadyaligned = AllocateIntVec( nseq );
1139        if( constraint )
1140        {
1141                localhomshrink = (LocalHom ***)calloc( nseq, sizeof( LocalHom ** ) );
1142#if SMALLMEMORY
1143                if( multidist )
1144                {
1145                        for( i=0; i<nseq; i++) localhomshrink[i] = (LocalHom **)calloc( 1, sizeof( LocalHom *) );
1146                }
1147                else
1148#endif
1149                {
1150                        for( i=0; i<nseq; i++) localhomshrink[i] = (LocalHom **)calloc( nseq, sizeof( LocalHom *) );
1151                }
1152        }
1153        effarr1_kozo = AllocateDoubleVec( nseq ); //tsuneni allocate sareru.
1154        effarr2_kozo = AllocateDoubleVec( nseq ); //tsuneni allocate sareru.
1155        for( i=0; i<nseq; i++ ) effarr1_kozo[i] = 0.0;
1156        for( i=0; i<nseq; i++ ) effarr2_kozo[i] = 0.0;
1157
1158        gaplen = AllocateIntVec( *alloclen+10 ); // maikai shokika
1159        gapmap = AllocateIntVec( *alloclen+10 ); // maikai shokika
1160        for( i=0; i<nseq-1; i++ ) alreadyaligned[i] = 1;
1161        alreadyaligned[nseq-1] = 0;
1162
1163        for( l=0; l<nseq; l++ ) fftlog[l] = 1;
1164
1165
1166        if( constraint ) 
1167        {
1168#if SMALLMEMORY
1169                if( multidist )
1170                        dontcalcimportance_firstone( nseq, effarr, aseq, localhomtable );
1171                else
1172                        calcimportance( nseq, effarr, aseq, localhomtable );
1173#else
1174                calcimportance( nseq, effarr, aseq, localhomtable );
1175#endif
1176        }
1177
1178        tscore = 0.0;
1179        for( l=0; l<nseq-1; l++ )
1180        {
1181                if( mergeoralign[l] == 'n' )
1182                {
1183//                      fprintf( stderr, "SKIP!\n" );
1184#if 0
1185                        free( topol[l][0] );
1186                        free( topol[l][1] );
1187                        free( topol[l] );
1188#endif
1189                        continue;
1190                }
1191
1192                m1 = topol[l][0][0];
1193                m2 = topol[l][1][0];
1194        len1 = strlen( aseq[m1] );
1195        len2 = strlen( aseq[m2] );
1196        if( *alloclen < len1 + len2 )
1197        {
1198#if 0
1199                        fprintf( stderr, "\nReallocating.." );
1200                        *alloclen = ( len1 + len2 ) + 1000;
1201                        ReallocateCharMtx( aseq, nseq, *alloclen + 10  );
1202                        gaplen = realloc( gaplen, ( *alloclen + 10 ) * sizeof( int ) );
1203                        if( gaplen == NULL )
1204                        {
1205                                fprintf( stderr, "Cannot realloc gaplen\n" );
1206                                exit( 1 );
1207                        }
1208                        gapmap = realloc( gapmap, ( *alloclen + 10 ) * sizeof( int ) );
1209                        if( gapmap == NULL )
1210                        {
1211                                fprintf( stderr, "Cannot realloc gapmap\n" );
1212                                exit( 1 );
1213                        }
1214                        fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
1215#else
1216                        fprintf( stderr, "Length over!\n" );
1217                        exit( 1 );
1218#endif
1219                }
1220
1221                if( effarr_kozo )
1222                {
1223                        clus1 = fastconjuction_noname_kozo( topol[l][0], aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1224                        clus2 = fastconjuction_noname_kozo( topol[l][1], aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1225                }
1226                else
1227                {
1228                        clus1 = fastconjuction_noname( topol[l][0], aseq, mseq1, effarr1, effarr, indication1 );
1229                        clus2 = fastconjuction_noname( topol[l][1], aseq, mseq2, effarr2, effarr, indication2 );
1230                }
1231
1232                if( mergeoralign[l] == '1' || mergeoralign[l] == '2' )
1233                {
1234                        newgapstr = "=";
1235                }
1236                else
1237                        newgapstr = "-";
1238
1239
1240                len1nocommongap = len1;
1241                len2nocommongap = len2;
1242                if( mergeoralign[l] == '1' ) // nai
1243                {
1244                        findcommongaps( clus2, mseq2, gapmap );
1245                        commongappick( clus2, mseq2 );
1246                        len2nocommongap = strlen( mseq2[0] );
1247                }
1248                else if( mergeoralign[l] == '2' )
1249                {
1250                        findcommongaps( clus1, mseq1, gapmap );
1251                        commongappick( clus1, mseq1 );
1252                        len1nocommongap = strlen( mseq1[0] );
1253                }
1254               
1255
1256//              fprintf( trap_g, "\nSTEP-%d\n", l );
1257//              fprintf( trap_g, "group1 = %s\n", indication1 );
1258//              fprintf( trap_g, "group2 = %s\n", indication2 );
1259//
1260#if 1
1261//              fprintf( stderr, "\rSTEP % 5d /%d ", l+1, nseq-1 );
1262//              fflush( stderr );
1263#else
1264                fprintf( stdout, "STEP %d /%d\n", l+1, nseq-1 );
1265                fprintf( stderr, "STEP %d /%d\n", l+1, nseq-1 );
1266                fprintf( stderr, "group1 = %.66s", indication1 );
1267                if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
1268                fprintf( stderr, "\n" );
1269                fprintf( stderr, "group2 = %.66s", indication2 );
1270                if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
1271                fprintf( stderr, "\n" );
1272#endif
1273
1274
1275
1276//              for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
1277
1278                if( constraint )
1279                {
1280#if SMALLMEMORY
1281                        if( multidist )
1282                        {
1283                                fastshrinklocalhom_one( topol[l][0], topol[l][1], nseq-1, localhomtable, localhomshrink );
1284                        }
1285                        else
1286#endif
1287                        {
1288                                fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
1289                        }
1290
1291//                      msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
1292//                      fprintf( stdout, "localhomshrink =\n" );
1293//                      outlocalhompt( localhomshrink, clus1, clus2 );
1294//                      weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
1295//                      fprintf( stderr, "after weight =\n" );
1296//                      outlocalhompt( localhomshrink, clus1, clus2 );
1297                }
1298                if( rnakozo && rnaprediction == 'm' )
1299                {
1300                        makegrouprna( grouprna1, singlerna, topol[l][0] );
1301                        makegrouprna( grouprna2, singlerna, topol[l][1] );
1302                }
1303
1304
1305/*
1306                fprintf( stderr, "before align all\n" );
1307                display( aseq, nseq );
1308                fprintf( stderr, "\n" );
1309                fprintf( stderr, "before align 1 %s \n", indication1 );
1310                display( mseq1, clus1 );
1311                fprintf( stderr, "\n" );
1312                fprintf( stderr, "before align 2 %s \n", indication2 );
1313                display( mseq2, clus2 );
1314                fprintf( stderr, "\n" );
1315*/
1316
1317
1318                if( !nevermemsave && ( constraint != 2  && alg != 'M'  && ( len1 > 30000 || len2 > 30000 ) ) )
1319                {
1320                        fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
1321                        alg = 'M';
1322                        if( commonIP ) FreeIntMtx( commonIP );
1323                        commonIP = NULL; // 2013/Jul17
1324                        commonAlloc1 = 0;
1325                        commonAlloc2 = 0;
1326                }
1327
1328
1329//              if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
1330                if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
1331                else                                               ffttry = 0;
1332//              ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
1333//              fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
1334//              fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
1335                if( constraint == 2 )
1336                {
1337                        if( alg == 'M' )
1338                        {
1339                                fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
1340                                exit( 1 );
1341                        }
1342                        fprintf( stderr, "c" );
1343                        if( alg == 'A' )
1344                        {
1345                                imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1346                                if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1347                                pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1348                        }
1349                        else if( alg == 'H' )
1350                        {
1351                                imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1352                                pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1353                        }
1354                        else if( alg == 'Q' )
1355                        {
1356                                imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1357                                if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1358                                pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1359                        }
1360                        else if( alg == 'R' )
1361                        {
1362                                imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1363                                pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1364                        }
1365                }
1366                else if( force_fft || ( use_fft && ffttry ) )
1367                {
1368                        fprintf( stderr, "f" );
1369                        if( alg == 'M' )
1370                        {
1371                                fprintf( stderr, "m" );
1372                                pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
1373                        }
1374                        else
1375                                pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
1376                }
1377                else
1378                {
1379                        fprintf( stderr, "d" );
1380                        fftlog[m1] = 0;
1381                        switch( alg )
1382                        {
1383                                case( 'a' ):
1384                                        pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
1385                                        break;
1386                                case( 'M' ):
1387                                        fprintf( stderr, "m" );
1388                                        pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1389                                        break;
1390                                case( 'A' ):
1391                                        pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1392                                        break;
1393                                case( 'Q' ):
1394                                        pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1395                                        break;
1396                                case( 'R' ):
1397                                        pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1398                                        break;
1399                                case( 'H' ):
1400                                        pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1401                                        break;
1402                                default:
1403                                        ErrorExit( "ERROR IN SOURCE FILE" );
1404                        }
1405                }
1406
1407
1408                nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
1409
1410//              fprintf( stderr, "aseq[last] = %s\n", aseq[nseq-1] );
1411
1412#if SCOREOUT
1413                fprintf( stderr, "score = %10.2f\n", pscore );
1414#endif
1415                tscore += pscore;
1416/*
1417                fprintf( stderr, "after align 1 %s \n", indication1 );
1418                display( mseq1, clus1 );
1419                fprintf( stderr, "\n" );
1420                fprintf( stderr, "after align 2 %s \n", indication2 );
1421                display( mseq2, clus2 );
1422                fprintf( stderr, "\n" );
1423*/
1424
1425//              writePre( nseq, name, nlen, aseq, 0 );
1426
1427                if( disp ) display( aseq, nseq );
1428
1429                if( mergeoralign[l] == '1' ) // jissainiha nai. atarashii hairetsu ha saigo dakara.
1430                {
1431                        adjustgapmap( strlen( mseq2[0] )-len2nocommongap+len2, gapmap, mseq2[0] );
1432                        restorecommongaps( nseq, aseq, topol[l][0], topol[l][1], gapmap, *alloclen, '-' );
1433                        findnewgaps( clus2, 0, mseq2, gaplen );
1434                        insertnewgaps( nseq, alreadyaligned, aseq, topol[l][1], topol[l][0], gaplen, gapmap, *alloclen, alg, '-' );
1435//                      for( i=0; i<nseq; i++ ) eq2dash( aseq[i] );
1436                        for( i=0; (m=topol[l][0][i])>-1; i++ ) alreadyaligned[m] = 1;
1437                }
1438                if( mergeoralign[l] == '2' )
1439                {
1440                        adjustgapmap( strlen( mseq1[0] )-len1nocommongap+len1, gapmap, mseq1[0] );
1441                        restorecommongaps( nseq, aseq, topol[l][0], topol[l][1], gapmap, *alloclen, '-' );
1442                        findnewgaps( clus1, 0, mseq1, gaplen );
1443                        insertnewgaps( nseq, alreadyaligned, aseq, topol[l][0], topol[l][1], gaplen, gapmap, *alloclen, alg, '-' );
1444//                      for( i=0; i<nseq; i++ ) eq2dash( aseq[i] );
1445                        for( i=0; (m=topol[l][1][i])>-1; i++ ) alreadyaligned[m] = 1;
1446                }
1447
1448#if 0
1449                free( topol[l][0] );
1450                free( topol[l][1] );
1451                free( topol[l] );
1452#endif
1453        }
1454
1455#if SCOREOUT
1456        fprintf( stderr, "totalscore = %10.2f\n\n", tscore );
1457#endif
1458        free( gaplen );
1459        free( gapmap );
1460        if( rnakozo && rnaprediction == 'm' )
1461        {
1462                free( grouprna1 );
1463                free( grouprna2 );
1464        }
1465        free( fftlog ); // iranai
1466        free( effarr1 );
1467        free( effarr2 );
1468        free( indication1 );
1469        free( indication2 );
1470        free( alreadyaligned );
1471        if( constraint )
1472        {
1473                for( i=0; i<nseq; i++ ) free( localhomshrink[i] ); // ??
1474                free( localhomshrink );
1475        }
1476        free( effarr1_kozo );
1477        free( effarr2_kozo );
1478        return( pscore );
1479}
1480
1481
1482
1483
1484static void mtxcpy( int norg, int njobc, float ***iscorec, float **iscore )
1485{
1486        int i, nlim, n;
1487        float *fpt, *fptc;
1488       
1489        *iscorec = AllocateFloatHalfMtx( njobc );
1490        nlim = norg-1;
1491        for( i=0; i<nlim; i++ )
1492        {
1493                fptc = (*iscorec)[i]+1;
1494                fpt  = iscore[i]+1;
1495                n = norg-i-1;
1496                while( n-- )
1497                        *fptc++ = *fpt++;
1498//              for( j=i+1; j<norg; j++ )       
1499//                      (*iscorec)[i][j-i] = iscore[i][j-i];
1500        }
1501}
1502
1503
1504static void     *addsinglethread( void *arg )
1505        {
1506                thread_arg_t *targ = (thread_arg_t *)arg;
1507                int *nlenc;
1508                char **namec;
1509                Treedep *depc;
1510                char **mseq1, **mseq2;
1511                float **iscorec;
1512//              float **iscorecbk; // to speedup
1513                double *effc;
1514                int ***topolc;
1515                float **lenc;
1516                LocalHom **localhomtablec = NULL;
1517                int *memlist0;
1518                int *memlist1;
1519                int *addmem;
1520                int njobc, norg;
1521                char **bseq;
1522                int i, j, iadd, rep, neighbor;
1523                char *mergeoralign;
1524
1525#ifdef enablemultithread
1526                int thread_no = targ->thread_no;
1527                int *iaddshare = targ->iaddshare; 
1528#endif
1529                int njob = targ->njob; 
1530                int *follows = targ->follows;
1531                int nadd = targ->nadd;
1532                int *nlen = targ->nlen;
1533                char **name = targ->name;
1534                char **seq = targ->seq;
1535                LocalHom **localhomtable = targ->localhomtable;
1536                float **iscore = targ->iscore;
1537                float **nscore = targ->nscore;
1538                int *istherenewgap = targ->istherenewgap;
1539                int **newgaplist = targ->newgaplist;
1540                RNApair ***singlerna = targ->singlerna;
1541                double *eff_kozo_mapped = targ->eff_kozo_mapped;
1542                int alloclen = targ->alloclen;
1543                Treedep *dep = targ->dep;
1544                int ***topol = targ->topol;
1545                float **len = targ->len;
1546                Addtree *addtree = targ->addtree;
1547                float pscore;
1548
1549
1550//              fprintf( stderr, "\nPreparing thread %d\n", thread_no );
1551                norg = njob - nadd;
1552                njobc = norg+1;
1553
1554                addmem = AllocateIntVec( nadd+1 );
1555                depc = (Treedep *)calloc( njobc, sizeof( Treedep ) );
1556                mseq1 = AllocateCharMtx( njob, 0 );
1557                mseq2 = AllocateCharMtx( njob, 0 );
1558                bseq = AllocateCharMtx( njobc, alloclen );
1559                namec = AllocateCharMtx( njob, 0 );
1560                nlenc = AllocateIntVec( njob );
1561                mergeoralign = AllocateCharVec( njob );
1562                if( constraint )
1563                {
1564                        localhomtablec = (LocalHom **)calloc( njobc, sizeof( LocalHom *) ); // motto chiisaku dekiru.
1565#if SMALLMEMORY
1566                        if( multidist )
1567                        {
1568                                for( i=0; i<njobc; i++) localhomtablec[i] = (LocalHom *)calloc( 1, sizeof( LocalHom ) ); // motto chiisaku dekiru.
1569                        }
1570                        else
1571#endif
1572                        {
1573                                for( i=0; i<njobc; i++) localhomtablec[i] = (LocalHom *)calloc( njobc, sizeof( LocalHom ) ); // motto chiisaku dekiru.
1574                                for( i=0; i<norg; i++ ) for( j=0; j<norg; j++ ) localhomtablec[i][j] = localhomtable[i][j]; // iru!
1575                        }
1576                }
1577
1578
1579                topolc = AllocateIntCub( njobc, 2, 0 );
1580                lenc = AllocateFloatMtx( njobc, 2 );
1581                effc = AllocateDoubleVec( njobc );
1582//              for( i=0; i<norg; i++ ) nlenc[i] = strlen( seq[i] );
1583                for( i=0; i<norg; i++ ) nlenc[i] = nlen[i];
1584                for( i=0; i<norg; i++ ) namec[i] = name[i];
1585                memlist0 = AllocateIntVec( norg+1 );
1586                memlist1 = AllocateIntVec( 2 );
1587                for( i=0; i<norg; i++ ) memlist0[i] = i;
1588                memlist0[norg] = -1;
1589
1590//              fprintf( stderr, "\ndone. %d\n", thread_no );
1591
1592//              mtxcpy( norg, norg, &iscorecbk, iscore ); // to speedup?
1593
1594
1595                iadd = -1;
1596                while( 1 )
1597                {
1598#ifdef enablemultithread
1599                        if( nthread )
1600                        {
1601                                pthread_mutex_lock( targ->mutex_counter );
1602                                iadd = *iaddshare;
1603                                if( iadd == nadd )
1604                                {
1605                                        pthread_mutex_unlock( targ->mutex_counter );
1606                                        break;
1607                                }
1608                                fprintf( stderr, "\r%d / %d (thread %d)                    \r", iadd, nadd, thread_no );
1609                                ++(*iaddshare);
1610                                pthread_mutex_unlock( targ->mutex_counter );
1611                        }
1612                        else
1613#endif
1614                        {
1615                                iadd++;
1616                                if( iadd == nadd ) break;
1617                                fprintf( stderr, "\r%d / %d                    \r", iadd, nadd );
1618                        }
1619
1620                        for( i=0; i<norg; i++ ) strcpy( bseq[i], seq[i] );
1621                        gappick0( bseq[norg], seq[norg+iadd] );
1622
1623                        mtxcpy( norg, njobc, &iscorec, iscore );
1624       
1625                        if( multidist || tuplesize > 0 )
1626                        {
1627                                for( i=0; i<norg; i++ ) iscorec[i][norg-i] = nscore[i][iadd];
1628                        }
1629                        else
1630                        {
1631                                for( i=0; i<norg; i++ ) iscorec[i][norg-i] = iscore[i][norg+iadd-i];
1632                        }
1633
1634
1635#if 0
1636                        for( i=0; i<njobc-1; i++ )
1637                        {
1638                                fprintf( stderr, "i=%d\n", i );
1639                                for( j=i+1; j<njobc; j++ )
1640                                {
1641                                        fprintf( stderr, "%d-%d, %f\n", i, j, iscorec[i][j-i] );
1642                                }
1643                        }
1644#endif
1645                        nlenc[norg] = nlen[norg+iadd];
1646                        namec[norg] = name[norg+iadd];
1647                        if( constraint) 
1648                        {
1649                                for( i=0; i<norg; i++ )
1650                                {
1651#if SMALLMEMORY
1652                                        if( multidist )
1653                                        {
1654                                                localhomtablec[i][0] = localhomtable[i][iadd];
1655//                                              localhomtablec[norg][i] = localhomtable[norg+iadd][i];
1656                                        }
1657                                        else
1658#endif
1659                                        {
1660                                                localhomtablec[i][norg] = localhomtable[i][norg+iadd];
1661                                                localhomtablec[norg][i] = localhomtable[norg+iadd][i];
1662                                        }
1663                                }
1664//                              localhomtablec[norg][norg] = localhomtable[norg+iadd][norg+iadd]; // iranai!!
1665                        }
1666       
1667//                      fprintf( stderr, "Constructing a UPGMA tree %d ... ", iadd );
1668//                      fflush( stderr );
1669       
1670
1671//                      if( iadd == 0 )
1672//                      {
1673//                      }
1674//                      fixed_musclesupg_float_realloc_nobk_halfmtx( njobc, iscorec, topolc, lenc, depc, 0 );
1675                        neighbor = addonetip( njobc, topolc, lenc, iscorec, topol, len, dep, treeout, addtree, iadd, name );
1676
1677                        if( noalign ) 
1678                        {
1679                                FreeFloatHalfMtx( iscorec, njobc );
1680                                continue;
1681                        }
1682       
1683                        if( tbrweight )
1684                        {
1685                                weight = 3; 
1686                                counteff_simple_float_nostatic( njobc, topolc, lenc, effc );
1687                        }
1688                        else
1689                        {
1690                                for( i=0; i<njobc; i++ ) effc[i] = 1.0;
1691                        }
1692               
1693                        FreeFloatHalfMtx( iscorec, njobc );
1694//                      FreeFloatMtx( lenc );
1695               
1696#if 0
1697                        for( i=0; i<njobc-1; i++ )
1698                        {
1699                                fprintf( stderr, "topol[%d] = \n", i );
1700                                for( j=0; topolc[i][0][j]!=-1; j++ ) fprintf( stderr, "%d ", topolc[i][0][j] );
1701                                fprintf( stderr, "\n" );
1702                                for( j=0; topolc[i][1][j]!=-1; j++ ) fprintf( stderr, "%d ", topolc[i][1][j] );
1703                                fprintf( stderr, "\n" );
1704                        }
1705
1706                        fprintf( stderr, "\nneighbor = %d, iadd = %d\n", neighbor, iadd );
1707#endif
1708                        follows[iadd] = neighbor;
1709
1710                        for( i=0; i<njobc-1; i++ ) mergeoralign[i] = 'n';
1711                        for( j=njobc-1; j<njobc; j++ )
1712                        {
1713                                addmem[0] = j;
1714                                addmem[1] = -1;
1715                                for( i=0; i<njobc-1; i++ )
1716                                {
1717                                        if( samemember( topolc[i][0], addmem ) ) // arieru
1718                                        {
1719//                                              fprintf( stderr, "HIT!\n" );
1720                                                if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1721                                                else mergeoralign[i] = '1';
1722                                        }
1723                                        else if( samemember( topolc[i][1], addmem ) )
1724                                        {
1725//                                              fprintf( stderr, "HIT!\n" );
1726                                                if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1727                                                else mergeoralign[i] = '2';
1728                                        }
1729                                }
1730                        }
1731       
1732//                      for( i=0; i<1; i++ ) addmem[i] = njobc-1+i;
1733                        addmem[0] = njobc-1;
1734                        addmem[1] = -1;
1735                        for( i=0; i<njobc-1; i++ )
1736                        {
1737                                if( includemember( topolc[i][0], addmem ) && includemember( topolc[i][1], addmem ) )
1738                                {
1739                                        mergeoralign[i] = 'w';
1740                                }
1741                                else if( includemember( topolc[i][0], addmem ) )
1742                                {
1743                                        mergeoralign[i] = '1';
1744//                                      fprintf( stderr, "HIT 1! iadd=%d", iadd );
1745                                }
1746                                else if( includemember( topolc[i][1], addmem ) )
1747                                {
1748                                        mergeoralign[i] = '2';
1749//                                      fprintf( stderr, "HIT 2! iadd=%d", iadd );
1750                                }
1751                        }
1752#if 0
1753                        for( i=0; i<njob-1; i++ )
1754                        {
1755                                fprintf( stderr, "mem0 = " );
1756                                for( j=0; topol[i][0][j]>-1; j++ )      fprintf( stderr, "%d ", topol[i][0][j] );
1757                                fprintf( stderr, "\n" );
1758                                fprintf( stderr, "mem1 = " );
1759                                for( j=0; topol[i][1][j]>-1; j++ )      fprintf( stderr, "%d ", topol[i][1][j] );
1760                                fprintf( stderr, "\n" );
1761                                fprintf( stderr, "i=%d, mergeoralign[] = %c\n", i, mergeoralign[i] );
1762                        }
1763#endif
1764
1765
1766#if 0
1767                        for( i=0; i<norg; i++ ) fprintf( stderr, "seq[%d, iadd=%d] = \n%s\n", i, iadd, seq[i] );
1768                        fprintf( stderr, "gapmapS (iadd=%d) = \n", iadd );
1769                        for( i=0; i<lennocommongap; i++ ) fprintf( stderr, "%d\n", gapmapS[i] );
1770#endif
1771
1772
1773//                      fprintf( stderr, "Progressive alignment ... \r" );
1774               
1775#if 0
1776                        pthread_mutex_lock( targ->mutex_counter );
1777                        fprintf( stdout, "\nmergeoralign (iadd=%d) = ", iadd );
1778                        for( i=0; i<njobc-1; i++ ) fprintf( stdout, "%c", mergeoralign[i] );
1779                        fprintf( stdout, "\n" );
1780                        pthread_mutex_unlock( targ->mutex_counter );
1781#endif
1782                        singlerna = NULL;
1783                        pscore = treebase( njobc, nlenc, bseq, 1, mergeoralign, mseq1, mseq2, topolc, effc, &alloclen, localhomtablec, singlerna, eff_kozo_mapped );
1784#if 0
1785                        pthread_mutex_lock( targ->mutex_counter );
1786//                      fprintf( stdout, "res (iadd=%d) = %s, pscore=%f\n", iadd, bseq[norg], pscore );
1787//                      fprintf( stdout, "effc (iadd=%d) = ", iadd );
1788//                      for( i=0; i<njobc; i++ ) fprintf( stdout, "%f ", effc[i] );
1789//                      fprintf( stdout, "\n" );
1790                        pthread_mutex_unlock( targ->mutex_counter );
1791#endif
1792       
1793       
1794#if 0
1795                        fprintf( trap_g, "done.\n" );
1796                        fclose( trap_g );
1797#endif
1798//                      fprintf( stdout, "\n>seq[%d, iadd=%d] = \n%s\n", norg+iadd, iadd, seq[norg+iadd] );
1799//                      fprintf( stdout, "\n>bseq[%d, iadd=%d] = \n%s\n", norg, iadd, bseq[norg] );
1800       
1801                        strcpy( seq[norg+iadd], bseq[norg] );
1802       
1803                        rep = -1;
1804                        for( i=0; i<norg; i++ )
1805                        {
1806//                              fprintf( stderr, "Checking %d/%d\n", i, norg );
1807                                if( strchr( bseq[i], '=' ) ) break;
1808                        }
1809                        if( i == norg ) 
1810                                istherenewgap[iadd] = 0;
1811                        else
1812                        {
1813                                rep = i;
1814                                istherenewgap[iadd] = 1;
1815                                makenewgaplist( newgaplist[iadd], bseq[rep] );
1816//                              for( i=0; newgaplist[iadd][i]!=-1; i++ ) fprintf( stderr, "%d: %d\n", i, newgaplist[iadd][i] );
1817                        }
1818                        eq2dash( seq[norg+iadd] );
1819
1820                }
1821
1822
1823                if( constraint )
1824                {
1825                        for( i=0; i<njobc; i++ ) free( localhomtablec[i] );
1826                        free( localhomtablec );
1827                }
1828                free( mergeoralign );
1829                FreeCharMtx( bseq );
1830                free( namec );
1831                free( nlenc  );
1832                free( depc );
1833                FreeIntCub( topolc );
1834                FreeFloatMtx( lenc );
1835                FreeDoubleVec( effc );
1836                free( memlist0 );
1837                free( memlist1 );
1838                free( addmem );
1839                free( mseq1 );
1840                free( mseq2 );
1841                Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
1842                A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
1843                if( commonIP ) FreeIntMtx( commonIP );
1844                commonIP = NULL;
1845                commonAlloc1 = commonAlloc2 = 0;
1846//              FreeFloatHalfMtx( iscorecbk, norg );
1847
1848                return( NULL );
1849        }
1850
1851static int maxl;
1852static int tsize;
1853static int nunknown = 0;
1854
1855void seq_grp_nuc( int *grp, char *seq )
1856{
1857        int tmp;
1858        int *grpbk = grp;
1859        while( *seq )
1860        {
1861                tmp = amino_grp[(int)*seq++];
1862                if( tmp < 4 )
1863                        *grp++ = tmp;
1864                else
1865                        nunknown++;
1866        }
1867        *grp = END_OF_VEC;
1868        if( grp - grpbk < tuplesize )
1869        {
1870//              fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
1871//              exit( 1 );
1872                *grpbk = -1;
1873        }
1874}
1875
1876void seq_grp( int *grp, char *seq )
1877{
1878        int tmp;
1879        int *grpbk = grp;
1880        while( *seq )
1881        {
1882                tmp = amino_grp[(int)*seq++];
1883                if( tmp < 6 )
1884                        *grp++ = tmp;
1885                else
1886                        nunknown++;
1887        }
1888        *grp = END_OF_VEC;
1889        if( grp - grpbk < 6 )
1890        {
1891//              fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
1892//              exit( 1 );
1893                *grpbk = -1;
1894        }
1895}
1896
1897void makecompositiontable_p( int *table, int *pointt )
1898{
1899        int point;
1900
1901        while( ( point = *pointt++ ) != END_OF_VEC )
1902                table[point]++;
1903}
1904
1905int commonsextet_p( int *table, int *pointt )
1906{
1907        int value = 0;
1908        int tmp;
1909        int point;
1910        static TLS int *memo = NULL;
1911        static TLS int *ct = NULL;
1912        int *cp;
1913
1914        if( table == NULL )
1915        {
1916                if( memo ) free( memo );
1917                if( ct ) free( ct );
1918                return( 0 );
1919        }
1920
1921        if( *pointt == -1 )
1922                return( 0 );
1923
1924        if( !memo )
1925        {
1926                memo = (int *)calloc( tsize, sizeof( int ) );
1927                if( !memo ) ErrorExit( "Cannot allocate memo\n" );
1928                ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) ); // chuui!!
1929                if( !ct ) ErrorExit( "Cannot allocate ct\n" );
1930        }
1931
1932        cp = ct;
1933        while( ( point = *pointt++ ) != END_OF_VEC )
1934        {
1935                tmp = memo[point]++;
1936                if( tmp < table[point] )
1937                        value++;
1938                if( tmp == 0 ) *cp++ = point;
1939        }
1940        *cp = END_OF_VEC;
1941       
1942        cp =  ct;
1943        while( *cp != END_OF_VEC )
1944                memo[*cp++] = 0;
1945
1946        return( value );
1947}
1948
1949void makepointtable_nuc_dectet( int *pointt, int *n )
1950{
1951        int point;
1952        register int *p;
1953
1954        if( *n == -1 )
1955        {
1956                *pointt = -1;
1957                return;
1958        }
1959
1960        p = n;
1961        point  = *n++ *262144;
1962        point += *n++ * 65536;
1963        point += *n++ * 16384;
1964        point += *n++ *  4096;
1965        point += *n++ *  1024;
1966        point += *n++ *   256;
1967        point += *n++ *    64;
1968        point += *n++ *    16;
1969        point += *n++ *     4;
1970        point += *n++;
1971        *pointt++ = point;
1972
1973        while( *n != END_OF_VEC )
1974        {
1975                point -= *p++ *262144;
1976                point *= 4;
1977                point += *n++;
1978                *pointt++ = point;
1979
1980        }
1981        *pointt = END_OF_VEC;
1982}
1983
1984void makepointtable_nuc_octet( int *pointt, int *n )
1985{
1986        int point;
1987        register int *p;
1988
1989        if( *n == -1 )
1990        {
1991                *pointt = -1;
1992                return;
1993        }
1994
1995        p = n;
1996        point  = *n++ * 16384;
1997        point += *n++ *  4096;
1998        point += *n++ *  1024;
1999        point += *n++ *   256;
2000        point += *n++ *    64;
2001        point += *n++ *    16;
2002        point += *n++ *     4;
2003        point += *n++;
2004        *pointt++ = point;
2005
2006        while( *n != END_OF_VEC )
2007        {
2008                point -= *p++ * 16384;
2009                point *= 4;
2010                point += *n++;
2011                *pointt++ = point;
2012        }
2013        *pointt = END_OF_VEC;
2014}
2015
2016void makepointtable_nuc( int *pointt, int *n )
2017{
2018        int point;
2019        register int *p;
2020
2021        if( *n == -1 )
2022        {
2023                *pointt = -1;
2024                return;
2025        }
2026
2027        p = n;
2028        point  = *n++ *  1024;
2029        point += *n++ *   256;
2030        point += *n++ *    64;
2031        point += *n++ *    16;
2032        point += *n++ *     4;
2033        point += *n++;
2034        *pointt++ = point;
2035
2036        while( *n != END_OF_VEC )
2037        {
2038                point -= *p++ * 1024;
2039                point *= 4;
2040                point += *n++;
2041                *pointt++ = point;
2042        }
2043        *pointt = END_OF_VEC;
2044}
2045
2046void makepointtable( int *pointt, int *n )
2047{
2048        int point;
2049        register int *p;
2050
2051        if( *n == -1 )
2052        {
2053                *pointt = -1;
2054                return;
2055        }
2056
2057        p = n;
2058        point  = *n++ *  7776;
2059        point += *n++ *  1296;
2060        point += *n++ *   216;
2061        point += *n++ *    36;
2062        point += *n++ *     6;
2063        point += *n++;
2064        *pointt++ = point;
2065
2066        while( *n != END_OF_VEC )
2067        {
2068                point -= *p++ * 7776;
2069                point *= 6;
2070                point += *n++;
2071                *pointt++ = point;
2072        }
2073        *pointt = END_OF_VEC;
2074}
2075
2076#ifdef enablemultithread
2077
2078void *dndprethread( void *arg )
2079{
2080        dndprethread_arg_t *targ = (dndprethread_arg_t *)arg;
2081        int njob = targ->njob;
2082        int thread_no = targ->thread_no;
2083        float *selfscore = targ->selfscore;
2084        float **mtx = targ->mtx;
2085        char **seq = targ->seq;
2086        Jobtable2d *jobpospt = targ->jobpospt;
2087
2088        int i, j;
2089        float ssi, ssj, bunbo;
2090        float mtxv;
2091
2092        if( njob == 1 ) return( NULL );
2093       
2094        while( 1 )
2095        {
2096                pthread_mutex_lock( targ->mutex );
2097                j = jobpospt->j;
2098                i = jobpospt->i;
2099                j++;
2100//              fprintf( stderr, "\n i=%d, j=%d before check\n", i, j );
2101                if( j == njob )
2102                {
2103//                      fprintf( stderr, "\n j = %d, i = %d, njob = %d\n", j, i, njob );
2104                        fprintf( stderr, "%4d/%4d (thread %4d), dndpre\r", i+1, njob, thread_no );
2105                        i++;
2106                        j = i + 1;
2107                        if( i == njob-1 )
2108                        {
2109//                              fprintf( stderr, "\n i=%d, njob-1=%d\n", i, njob-1 );
2110                                pthread_mutex_unlock( targ->mutex );
2111                                return( NULL );
2112                        }
2113                }
2114//              fprintf( stderr, "\n i=%d, j=%d after check\n", i, j );
2115                jobpospt->j = j;
2116                jobpospt->i = i;
2117                pthread_mutex_unlock( targ->mutex );
2118
2119                ssi = selfscore[i];
2120                ssj = selfscore[j];
2121
2122                bunbo = MIN( ssi, ssj );
2123                if( bunbo == 0.0 )
2124                        mtxv = maxdist;
2125                else
2126                        mtxv = maxdist * ( 1.0 - (float)naivepairscore11( seq[i], seq[j], penalty * 10  ) / bunbo );
2127#if 1
2128                if( mtxv > 9.0 || mtxv < 0.0 )
2129                {
2130                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
2131                        exit( 1 );
2132                }
2133#else // CHUUI!!!  2012/05/16
2134                if( mtxv > 2.0 )
2135                {
2136                        mtxv = 2.0;
2137                }
2138                if( mtxv < 0.0 )
2139                {
2140                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
2141                        exit( 1 );
2142                }
2143#endif
2144                mtx[i][j-i] = mtxv;
2145        }
2146}
2147
2148static void *gaplist2alnxthread( void *arg )
2149{
2150        gaplist2alnxthread_arg_t *targ = (gaplist2alnxthread_arg_t *)arg;
2151//      int thread_no = targ->thread_no;
2152        int ncycle = targ->ncycle;
2153        char **seq = targ->seq;
2154        int *newgaplist = targ->newgaplist;
2155        int *posmap = targ->posmap;
2156        int *jobpospt = targ->jobpospt;
2157        int tmpseqlen = targ->tmpseqlen;
2158        int lenfull = targ->lenfull;
2159        char *tmpseq1;
2160        int i;
2161
2162        tmpseq1 = AllocateCharVec( tmpseqlen );
2163
2164        while( 1 )
2165        {
2166                pthread_mutex_lock( targ->mutex );
2167                i = *jobpospt;
2168                if( i == ncycle )
2169                {
2170                        pthread_mutex_unlock( targ->mutex );
2171                        free( tmpseq1 );
2172                        return( NULL );
2173                }
2174                *jobpospt = i+1;
2175                pthread_mutex_unlock( targ->mutex );
2176
2177                gaplist2alnx( lenfull, tmpseq1, seq[i], newgaplist, posmap, tmpseqlen  );
2178//              fprintf( stderr, ">%s (iadd=%d)\n%s\n", name[i], iadd, tmpseq1 );
2179                strcpy( seq[i], tmpseq1 );
2180        }
2181}
2182
2183static void *distancematrixthread( void *arg )
2184{
2185        distancematrixthread_arg_t *targ = (distancematrixthread_arg_t *)arg;
2186        int thread_no = targ->thread_no;
2187        int njob = targ->njob;
2188        int norg = targ->norg;
2189        int *jobpospt = targ->jobpospt;
2190        int **pointt = targ->pointt;
2191        float **imtx = targ->imtx;
2192        float **nmtx = targ->nmtx;
2193        float *selfscore = targ->selfscore;
2194        int *nogaplen = targ->nogaplen;
2195
2196        float lenfac, bunbo, longer, shorter, mtxv;
2197        int *table1;
2198        int i, j;
2199
2200        while( 1 )
2201        {
2202                pthread_mutex_lock( targ->mutex );
2203                i = *jobpospt;
2204                if( i == norg )
2205                {
2206                        pthread_mutex_unlock( targ->mutex );
2207                        commonsextet_p( NULL, NULL );
2208                        return( NULL );
2209                }
2210                *jobpospt = i+1;
2211                pthread_mutex_unlock( targ->mutex );
2212
2213                table1 = (int *)calloc( tsize, sizeof( int ) );
2214                if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
2215                if( i % 100 == 0 )
2216                {
2217                        fprintf( stderr, "\r% 5d / %d (thread %4d)", i+1, norg, thread_no );
2218                }
2219                makecompositiontable_p( table1, pointt[i] );
2220       
2221                for( j=i+1; j<njob; j++ ) 
2222                {
2223                        mtxv = (float)commonsextet_p( table1, pointt[j] );
2224                        if( nogaplen[i] > nogaplen[j] )
2225                        {
2226                                longer=(float)nogaplen[i];
2227                                shorter=(float)nogaplen[j];
2228                        }
2229                        else
2230                        {
2231                                longer=(float)nogaplen[j];
2232                                shorter=(float)nogaplen[i];
2233                        }
2234                        lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
2235                        bunbo = MIN( selfscore[i], selfscore[j] );
2236
2237                        if( j < norg )
2238                        {
2239                                if( bunbo == 0.0 )
2240                                        imtx[i][j-i] = 1.0;
2241                                else
2242                                        imtx[i][j-i] = ( 1.0 - mtxv / bunbo ) * lenfac;
2243
2244                        }
2245                        else
2246                        {
2247                                if( bunbo == 0.0 )
2248                                        nmtx[i][j-norg] = 1.0;
2249                                else
2250                                        nmtx[i][j-norg] = ( 1.0 - mtxv / bunbo ) * lenfac;
2251                        }
2252                } 
2253                free( table1 );
2254
2255//              for( j=i+1; j<norg; j++ )
2256//                      imtx[i][j-i] = (float)commonsextet_p( table1, pointt[j] );
2257//              for( j=norg; j<njob; j++ )
2258//                      nmtx[i][j-norg] = (float)commonsextet_p( table1, pointt[j] );
2259//              free( table1 );
2260        }
2261}
2262#endif
2263
2264
2265void ktupledistancematrix( int nseq, int norg, int nlenmax, char **seq, char **name, float **imtx, float **nmtx )
2266{
2267        char *tmpseq;
2268        int *grpseq;
2269        int **pointt;
2270        int i, j;
2271        int *nogaplen;
2272        int *table1;
2273        float lenfac, bunbo, longer, shorter, mtxv;
2274        float *selfscore;
2275        selfscore = AllocateFloatVec( nseq );
2276
2277        fprintf( stderr, "\n\nMaking a distance matrix ..\n" );
2278        fflush( stderr );
2279
2280        tmpseq = AllocateCharVec( nlenmax+1 );
2281        grpseq = AllocateIntVec( nlenmax+1 );
2282        pointt = AllocateIntMtx( nseq, nlenmax+1 );
2283        nogaplen = AllocateIntVec( nseq ); 
2284
2285        if( dorp == 'd' ) tsize = (int)pow( 4, tuplesize );
2286        else              tsize = (int)pow( 6, 6 );
2287
2288        if( dorp == 'd' && tuplesize == 6 )
2289        {
2290                lenfaca = D6LENFACA;
2291                lenfacb = D6LENFACB;
2292                lenfacc = D6LENFACC;
2293                lenfacd = D6LENFACD;
2294        }
2295        else if( dorp == 'd' && tuplesize == 10 )
2296        {
2297                lenfaca = D10LENFACA;
2298                lenfacb = D10LENFACB;
2299                lenfacc = D10LENFACC;
2300                lenfacd = D10LENFACD;
2301        }
2302        else   
2303        {
2304                lenfaca = PLENFACA;
2305                lenfacb = PLENFACB;
2306                lenfacc = PLENFACC;
2307                lenfacd = PLENFACD;
2308        }
2309
2310        maxl = 0;
2311        for( i=0; i<nseq; i++ ) 
2312        {
2313                gappick0( tmpseq, seq[i] );
2314                nogaplen[i] = strlen( tmpseq );
2315                if( nogaplen[i] < 6 )
2316                {
2317//                      fprintf( stderr, "Seq %d, too short, %d characters\n", i+1, nogaplen[i] );
2318//                      fprintf( stderr, "Please use mafft-ginsi, mafft-linsi or mafft-ginsi\n\n\n" );
2319//                      exit( 1 );
2320                }
2321                if( nogaplen[i] > maxl ) maxl = nogaplen[i];
2322                if( dorp == 'd' ) /* nuc */
2323                {
2324                        seq_grp_nuc( grpseq, tmpseq );
2325//                      makepointtable_nuc( pointt[i], grpseq );
2326//                      makepointtable_nuc_octet( pointt[i], grpseq );
2327                        if( tuplesize == 10 )
2328                                makepointtable_nuc_dectet( pointt[i], grpseq );
2329                        else if( tuplesize == 6 )
2330                                makepointtable_nuc( pointt[i], grpseq );
2331                        else
2332                        {
2333                                fprintf( stderr, "tuplesize=%d: not supported\n", tuplesize );
2334                                exit( 1 );
2335                        }
2336                }
2337                else                 /* amino */
2338                {
2339                        seq_grp( grpseq, tmpseq );
2340                        makepointtable( pointt[i], grpseq );
2341                }
2342
2343        }
2344        if( nunknown ) fprintf( stderr, "\nWARNING : %d unknown characters\n", nunknown );
2345
2346        for( i=0; i<nseq; i++ ) // serial de jubun
2347        {
2348                table1 = (int *)calloc( tsize, sizeof( int ) );
2349                if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
2350                makecompositiontable_p( table1, pointt[i] );
2351
2352                selfscore[i] = (float)commonsextet_p( table1, pointt[i] );
2353                free( table1 );
2354        }
2355
2356#ifdef enablemultithread
2357        if( nthread > 0 )
2358        {
2359                distancematrixthread_arg_t *targ; 
2360                int jobpos;
2361                pthread_t *handle;
2362                pthread_mutex_t mutex;
2363
2364                jobpos = 0; 
2365                targ = calloc( nthread, sizeof( distancematrixthread_arg_t ) ); 
2366                handle = calloc( nthread, sizeof( pthread_t ) ); 
2367                pthread_mutex_init( &mutex, NULL );
2368
2369                for( i=0; i<nthread; i++ )
2370                {
2371                        targ[i].thread_no = i;
2372                        targ[i].njob = nseq;
2373                        targ[i].norg = norg;
2374                        targ[i].jobpospt = &jobpos;
2375                        targ[i].pointt = pointt;
2376                        targ[i].imtx = imtx;
2377                        targ[i].nmtx = nmtx;
2378                        targ[i].selfscore = selfscore;
2379                        targ[i].nogaplen = nogaplen;
2380                        targ[i].mutex = &mutex;
2381
2382                        pthread_create( handle+i, NULL, distancematrixthread, (void *)(targ+i) );
2383                }
2384       
2385                for( i=0; i<nthread; i++ )
2386                {
2387                        pthread_join( handle[i], NULL );
2388                }
2389                pthread_mutex_destroy( &mutex );
2390                free( handle );
2391                free( targ );
2392
2393        }
2394        else
2395#endif
2396        {
2397                for( i=0; i<norg; i++ )
2398                {
2399                        table1 = (int *)calloc( tsize, sizeof( int ) );
2400                        if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
2401                        if( i % 100 == 0 )
2402                        {
2403                                fprintf( stderr, "\r% 5d / %d", i+1, norg );
2404                                fflush( stderr );
2405                        }
2406                        makecompositiontable_p( table1, pointt[i] );
2407       
2408                        for( j=i+1; j<nseq; j++ ) 
2409                        {
2410                                mtxv = (float)commonsextet_p( table1, pointt[j] );
2411                                if( nogaplen[i] > nogaplen[j] )
2412                                {
2413                                        longer=(float)nogaplen[i];
2414                                        shorter=(float)nogaplen[j];
2415                                }
2416                                else
2417                                {
2418                                        longer=(float)nogaplen[j];
2419                                        shorter=(float)nogaplen[i];
2420                                }
2421                                lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
2422                                bunbo = MIN( selfscore[i], selfscore[j] );
2423
2424                                if( j < norg )
2425                                {
2426                                        if( bunbo == 0.0 )
2427                                                imtx[i][j-i] = 1.0;
2428                                        else
2429                                                imtx[i][j-i] = ( 1.0 - mtxv / bunbo ) * lenfac;
2430                                }
2431                                else
2432                                {
2433                                        if( bunbo == 0.0 )
2434                                                nmtx[i][j-norg] = 1.0;
2435                                        else
2436                                                nmtx[i][j-norg] = ( 1.0 - mtxv / bunbo ) * lenfac;
2437                                }
2438                        } 
2439                        free( table1 );
2440                }
2441        }
2442
2443        fprintf( stderr, "\ndone.\n\n" );
2444        fflush( stderr );
2445
2446        for( i=0; i<norg; i++ )
2447        {
2448                for( j=i+1; j<norg; j++ ) 
2449                {
2450
2451                }
2452                for( j=norg; j<nseq; j++ ) 
2453                {
2454                }
2455        }
2456        free( grpseq );
2457        free( tmpseq );
2458        FreeIntMtx( pointt );
2459    free( nogaplen );
2460    free( selfscore );
2461
2462#if 0 // writehat2 wo kakinaosu
2463        if( distout )
2464        {
2465                hat2p = fopen( "hat2", "w" );
2466                WriteFloatHat2_pointer_halfmtx( hat2p, nseq, name, mtx );
2467                fclose( hat2p );
2468        }
2469#endif
2470}
2471
2472void dndpre( int nseq, char **seq, float **mtx ) // not used yet
2473{
2474        int i, j, ilim;
2475        float *selfscore;
2476        float mtxv;
2477        float ssi, ssj, bunbo;
2478
2479        selfscore = AllocateFloatVec( nseq );
2480
2481        for( i=0; i<nseq; i++ )
2482        {
2483                selfscore[i] = (float)naivepairscore11( seq[i], seq[i], 0 );
2484        }
2485#ifdef enablemultithread
2486        if( nthread > 0 )
2487        {
2488                dndprethread_arg_t *targ;
2489                Jobtable2d jobpos;
2490                pthread_t *handle;
2491                pthread_mutex_t mutex;
2492
2493                jobpos.i = 0;
2494                jobpos.j = 0;
2495
2496                targ = calloc( nthread, sizeof( dndprethread_arg_t ) );
2497                handle = calloc( nthread, sizeof( pthread_t ) );
2498                pthread_mutex_init( &mutex, NULL );
2499
2500                for( i=0; i<nthread; i++ )
2501                {
2502                        targ[i].thread_no = i;
2503                        targ[i].njob = nseq;
2504                        targ[i].selfscore = selfscore;
2505                        targ[i].mtx = mtx;
2506                        targ[i].seq = seq;
2507                        targ[i].jobpospt = &jobpos;
2508                        targ[i].mutex = &mutex;
2509
2510                        pthread_create( handle+i, NULL, dndprethread, (void *)(targ+i) );
2511                }
2512
2513                for( i=0; i<nthread; i++ )
2514                {
2515                        pthread_join( handle[i], NULL );
2516                }
2517                pthread_mutex_destroy( &mutex );
2518
2519        }
2520        else
2521#endif
2522        {
2523                ilim = nseq-1;
2524                for( i=0; i<ilim; i++ )
2525                {
2526                        ssi = selfscore[i];
2527                        fprintf( stderr, "%4d/%4d\r", i+1, nseq );
2528
2529                        for( j=i+1; j<nseq; j++ )
2530                        {
2531                                ssj = selfscore[j];
2532                                bunbo = MIN( ssi, ssj );
2533                                if( bunbo == 0.0 )
2534                                        mtxv = maxdist;
2535                                else
2536                                        mtxv = maxdist * ( 1.0 - (float)naivepairscore11( seq[i], seq[j], penalty * 10 ) / bunbo );
2537
2538#if 1
2539                                if( mtxv > 9.0 || mtxv < 0.0 )
2540                                {
2541                                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
2542                                        exit( 1 );
2543                                }
2544#else // CHUUI!!!  2012/05/16
2545                                if( mtxv > 2.0 )
2546                                {
2547                                        mtxv = 2.0;
2548                                }
2549                                if( mtxv < 0.0 )
2550                                {
2551                                        fprintf( stderr, "Distance %d-%d is strange, %f.\n", i, j, mtxv );
2552                                        exit( 1 );
2553                                }
2554#endif
2555                                mtx[i][j-i] = mtxv;
2556                        }
2557                }
2558        }
2559       
2560#if TEST
2561        for( i=0; i<nseq-1; i++ ) for( j=i+1; j<nseq; j++ ) 
2562                fprintf( stdout, "i=%d, j=%d, mtx[][] = %f\n", i, j, mtx[i][j] );
2563#endif
2564        free( selfscore );
2565
2566}
2567
2568int main( int argc, char *argv[] )
2569{
2570        static int  *nlen;     
2571        static char **name, **seq;
2572        static char  **tmpseq;
2573        static char  *tmpseq1;
2574//      static char  *check1, *check2;
2575        static float **iscore, **iscore_kozo;
2576        static double *eff_kozo, *eff_kozo_mapped = NULL;
2577        int i, j, f, ien;
2578        int iadd;
2579        static int ***topol_kozo;
2580        Treedep *dep;
2581        static float **len_kozo;
2582        FILE *prep;
2583        FILE *infp;
2584        FILE *hat2p;
2585        int alignmentlength;
2586        char c;
2587        int alloclen, fullseqlen, tmplen;
2588        LocalHom **localhomtable = NULL;
2589        static char *kozoarivec;
2590        int nkozo;
2591        int njobc, norg, lenfull;
2592        int **newgaplist_o;
2593        int *newgaplist_compact;
2594        int **follower;
2595        int *follows;
2596        int *istherenewgap;
2597        int zure;
2598        int *posmap;
2599        int *ordertable;
2600        FILE *orderfp;
2601        int tmpseqlen;
2602        Blocktorealign *realign;
2603        RNApair ***singlerna;
2604        int ***topol;
2605        float **len;
2606        float **iscoreo, **nscore;
2607        FILE *fp;
2608        Addtree *addtree;
2609
2610
2611        arguments( argc, argv );
2612#ifndef enablemultithread
2613        nthread = 0;
2614#endif
2615
2616        if( fastathreshold < 0.0001 ) constraint = 0;
2617
2618        if( inputfile )
2619        {
2620                infp = fopen( inputfile, "r" );
2621                if( !infp ) 
2622                {
2623                        fprintf( stderr, "Cannot open %s\n", inputfile );
2624                        exit( 1 );
2625                }
2626        }
2627        else   
2628                infp = stdin;
2629
2630        getnumlen( infp );
2631        rewind( infp );
2632
2633
2634        nkozo = 0;
2635
2636        if( njob < 2 )
2637        {
2638                fprintf( stderr, "At least 2 sequences should be input!\n"
2639                                                 "Only %d sequence found.\n", njob ); 
2640                exit( 1 );
2641        }
2642
2643        norg = njob-nadd;
2644        njobc = norg+1;
2645        fprintf( stderr, "norg = %d\n", norg );
2646        fprintf( stderr, "njobc = %d\n", njobc );
2647        if( norg > 1000 || nadd > 1000 ) use_fft = 0;
2648
2649        fullseqlen = alloclen = nlenmax*4+1; //chuui!
2650        seq = AllocateCharMtx( njob, alloclen );
2651
2652        name = AllocateCharMtx( njob, B+1 );
2653        nlen = AllocateIntVec( njob );
2654
2655
2656        if( multidist || tuplesize > 0 )
2657        {
2658                iscore = AllocateFloatHalfMtx( norg );
2659                nscore = AllocateFloatMtx( norg, nadd );
2660        }
2661        else
2662        {
2663                iscore = AllocateFloatHalfMtx( njob );
2664                nscore = NULL;
2665        }
2666
2667        kozoarivec = AllocateCharVec( njob );
2668
2669
2670        ordertable = AllocateIntVec( norg+1 );
2671
2672
2673        if( constraint )
2674        {
2675#if SMALLMEMORY
2676                if( multidist )
2677                {
2678                        localhomtable = (LocalHom **)calloc( norg, sizeof( LocalHom *) );
2679                        for( i=0; i<norg; i++)
2680                        {
2681                                localhomtable[i] = (LocalHom *)calloc( nadd, sizeof( LocalHom ) );
2682                                for( j=0; j<nadd; j++)
2683                                {
2684                                        localhomtable[i][j].start1 = -1;
2685                                        localhomtable[i][j].end1 = -1;
2686                                        localhomtable[i][j].start2 = -1;
2687                                        localhomtable[i][j].end2 = -1;
2688                                        localhomtable[i][j].overlapaa = -1.0;
2689                                        localhomtable[i][j].opt = -1.0;
2690                                        localhomtable[i][j].importance = -1.0;
2691                                        localhomtable[i][j].next = NULL;
2692                                        localhomtable[i][j].korh = 'h';
2693                                }
2694                        }
2695//                      localhomtable = (LocalHom **)calloc( norg+nadd, sizeof( LocalHom *) );
2696//                      for( i=norg; i<norg+nadd; i++) // hontou ha iranai
2697//                      {
2698//                              localhomtable[i] = (LocalHom *)calloc( norg, sizeof( LocalHom ) );
2699//                              for( j=0; j<norg; j++)
2700//                              {
2701//                                      localhomtable[i][j].start1 = -1;
2702//                                      localhomtable[i][j].end1 = -1;
2703//                                      localhomtable[i][j].start2 = -1;
2704//                                      localhomtable[i][j].end2 = -1;
2705//                                      localhomtable[i][j].overlapaa = -1.0;
2706//                                      localhomtable[i][j].opt = -1.0;
2707//                                      localhomtable[i][j].importance = -1.0;
2708//                                      localhomtable[i][j].next = NULL;
2709//                                      localhomtable[i][j].korh = 'h';
2710//                              }
2711//                      }
2712                }
2713                else
2714#endif
2715                {
2716                        localhomtable = (LocalHom **)calloc( njob, sizeof( LocalHom *) );
2717                        for( i=0; i<njob; i++)
2718                        {
2719                                localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
2720                                for( j=0; j<njob; j++)
2721                                {
2722                                        localhomtable[i][j].start1 = -1;
2723                                        localhomtable[i][j].end1 = -1;
2724                                        localhomtable[i][j].start2 = -1;
2725                                        localhomtable[i][j].end2 = -1;
2726                                        localhomtable[i][j].overlapaa = -1.0;
2727                                        localhomtable[i][j].opt = -1.0;
2728                                        localhomtable[i][j].importance = -1.0;
2729                                        localhomtable[i][j].next = NULL;
2730                                        localhomtable[i][j].korh = 'h';
2731                                }
2732                        }
2733                }
2734
2735                fprintf( stderr, "Loading 'hat3' ... " );
2736                prep = fopen( "hat3", "r" );
2737                if( prep == NULL ) ErrorExit( "Make hat3." );
2738#if SMALLMEMORY
2739                if( multidist )
2740                {
2741//                      readlocalhomtable_two( prep, norg, nadd, localhomtable, localhomtable+norg, kozoarivec );
2742                        readlocalhomtable_one( prep, norg, nadd, localhomtable, kozoarivec );
2743                }
2744                else
2745#endif
2746                {
2747                        readlocalhomtable( prep, njob, localhomtable, kozoarivec );
2748                }
2749
2750                fclose( prep );
2751                fprintf( stderr, "\ndone.\n" );
2752
2753
2754                nkozo = 0;
2755                for( i=0; i<njob; i++ ) 
2756                {
2757//                      fprintf( stderr, "kozoarivec[%d] = %d\n", i, kozoarivec[i] );
2758                        if( kozoarivec[i] ) nkozo++;
2759                }
2760                if( nkozo )
2761                {
2762                        topol_kozo = AllocateIntCub( nkozo, 2, 0 );
2763                        len_kozo = AllocateFloatMtx( nkozo, 2 );
2764                        iscore_kozo = AllocateFloatHalfMtx( nkozo );
2765                        eff_kozo = AllocateDoubleVec( nkozo );
2766                        eff_kozo_mapped = AllocateDoubleVec( njob );
2767                }
2768
2769
2770#if SMALLMEMORY
2771//              outlocalhom_part( localhomtable, norg, nadd );
2772#else
2773//              outlocalhom( localhomtable, njob );
2774#endif
2775
2776#if 0
2777                fprintf( stderr, "Extending localhom ... " );
2778                extendlocalhom2( njob, localhomtable );
2779                fprintf( stderr, "done.\n" );
2780#endif
2781        }
2782
2783#if 0
2784        readData( infp, name, nlen, seq );
2785#else
2786        readData_pointer( infp, name, nlen, seq );
2787        fclose( infp );
2788#endif
2789
2790        constants( njob, seq );
2791
2792#if 0
2793        fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
2794#endif
2795
2796        initSignalSM();
2797
2798        initFiles();
2799
2800//      WriteOptions( trap_g );
2801
2802        c = seqcheck( seq );
2803        if( c )
2804        {
2805                fprintf( stderr, "Illegal character %c\n", c );
2806                exit( 1 );
2807        }
2808
2809        alignmentlength = strlen( seq[0] );
2810        for( i=0; i<norg; i++ )
2811        {
2812                if( alignmentlength != strlen( seq[i] ) )
2813                {
2814                        fprintf( stderr, "#################################################################################\n" );
2815                        fprintf( stderr, "# ERROR!                                                                        #\n" );
2816                        fprintf( stderr, "# The original%4d sequences must be aligned                                    #\n", njob-nadd );
2817                        fprintf( stderr, "#################################################################################\n" );
2818                        exit( 1 );
2819                }
2820        }
2821        if( addprofile )
2822        {
2823                fprintf( stderr, "Not supported!\n" );
2824                exit( 1 );
2825        }
2826
2827        if( tuplesize > 0 ) // if mtx is internally computed
2828        {
2829                if( multidist == 1 )
2830                {
2831                        ktupledistancematrix( njob, norg, nlenmax, seq, name, iscore, nscore ); // iscore ha muda.
2832
2833//                      hat2p = fopen( "hat2-1", "w" );
2834//                      WriteFloatHat2_pointer_halfmtx( hat2p, njob, name, iscore );
2835//                      fclose( hat2p );
2836
2837                        dndpre( norg, seq, iscore );
2838//                      fprintf( stderr, "Loading 'hat2i' (aligned sequences) ... " );
2839//                      prep = fopen( "hat2i", "r" );
2840//                      if( prep == NULL ) ErrorExit( "Make hat2i." );
2841//                      readhat2_floathalf_pointer( prep, njob-nadd, name, iscore );
2842//                      fclose( prep );
2843//                      fprintf( stderr, "done.\n" );
2844
2845//                      hat2p = fopen( "hat2-2", "w" );
2846//                      WriteFloatHat2_pointer_halfmtx( hat2p, norg, name, iscore );
2847//                      fclose( hat2p );
2848                }
2849                else
2850                {
2851                        ktupledistancematrix( njob, norg, nlenmax, seq, name, iscore, nscore );
2852                }
2853        }
2854        else
2855        {
2856                if( multidist == 1 )
2857                {
2858                        fprintf( stderr, "Loading 'hat2n' (aligned sequences - new sequences) ... " );
2859                        prep = fopen( "hat2n", "r" );
2860                        if( prep == NULL ) ErrorExit( "Make hat2n." );
2861                        readhat2_floathalf_part_pointer( prep, njob, nadd, name, nscore );
2862                        fclose( prep );
2863                        fprintf( stderr, "done.\n" );
2864               
2865                        fprintf( stderr, "Loading 'hat2i' (aligned sequences) ... " );
2866                        prep = fopen( "hat2i", "r" );
2867                        if( prep == NULL ) ErrorExit( "Make hat2i." );
2868                        readhat2_floathalf_pointer( prep, njob-nadd, name, iscore );
2869                        fclose( prep );
2870                        fprintf( stderr, "done.\n" );
2871                }
2872                else
2873                {
2874                        fprintf( stderr, "Loading 'hat2' ... " );
2875                        prep = fopen( "hat2", "r" );
2876                        if( prep == NULL ) ErrorExit( "Make hat2." );
2877                        readhat2_floathalf_pointer( prep, njob, name, iscore );
2878                        fclose( prep );
2879                        fprintf( stderr, "done.\n" );
2880                }
2881        }
2882
2883#if 1
2884        if( distout )
2885        {
2886                fprintf( stderr, "Error in v6.936!! Please contact kazutaka.katoh@aist.go.jp\n" );
2887                exit( 1 );
2888                hat2p = fopen( "hat2", "w" );
2889                WriteFloatHat2_pointer_halfmtx( hat2p, norg, name, iscore );
2890                fclose( hat2p );
2891                exit( 1 );
2892        }
2893#endif
2894
2895
2896        singlerna = NULL;
2897
2898        commongappick( norg, seq );
2899
2900        lenfull = strlen( seq[0] );
2901
2902//      newgaplist_o = AllocateIntMtx( nadd, alloclen ); //ookisugi
2903        newgaplist_o = AllocateIntMtx( nadd, lenfull*2 );
2904        newgaplist_compact = AllocateIntVec( lenfull*2 );
2905        istherenewgap = AllocateIntVec( nadd );
2906        follower = AllocateIntMtx( norg, 1 );
2907        for( i=0; i<norg; i++ ) follower[i][0] = -1;
2908        follows = AllocateIntVec( nadd );
2909
2910        dep = (Treedep *)calloc( norg, sizeof( Treedep ) );
2911        topol = AllocateIntCub( norg, 2, 0 );
2912        len = AllocateFloatMtx( norg, 2 );
2913//      iscoreo = AllocateFloatHalfMtx( norg );
2914        mtxcpy( norg, norg, &iscoreo, iscore );
2915
2916        if( treeout )
2917        {
2918                addtree = (Addtree *)calloc( nadd, sizeof( Addtree ) );
2919                if( !addtree )
2920                {
2921                        fprintf( stderr, "Cannot allocate addtree\n" );
2922                        exit( 1 );
2923                }
2924        }
2925
2926
2927//      nlim = norg-1;
2928//      for( i=0; i<nlim; i++ )
2929//      {
2930//              fptc = iscoreo[i]+1;
2931//              fpt  = iscore[i]+1;
2932//              j = norg-i-1;
2933//              while( j-- )
2934//                      *fptc++ = *fpt++;
2935////    for( j=i+1; j<norg; j++ )       
2936////            iscoreo[i][j-i] = iscore[i][j-i];
2937//      }
2938
2939//      fprintf( stderr, "building a tree.." );
2940        if( treeout )
2941                fixed_musclesupg_float_realloc_nobk_halfmtx_treeout( norg, iscoreo, topol, len, name, nlen, dep );
2942        else
2943                fixed_musclesupg_float_realloc_nobk_halfmtx( norg, iscoreo, topol, len, dep, 0 );
2944//      fprintf( stderr, "done.\n" );
2945
2946        if( norg > 1 ) 
2947                cnctintvec( ordertable, topol[norg-2][0], topol[norg-2][1] );
2948        else
2949        {
2950                ordertable[0] = 0; ordertable[1] = -1;
2951        }
2952        FreeFloatHalfMtx( iscoreo, norg );
2953
2954#ifdef enablemultithread
2955        if( nthread )
2956        {
2957                pthread_t *handle;
2958                pthread_mutex_t mutex_counter;
2959                thread_arg_t *targ;
2960                int *iaddsharept;
2961
2962                targ = calloc( nthread, sizeof( thread_arg_t ) );
2963                handle = calloc( nthread, sizeof( pthread_t ) );
2964                pthread_mutex_init( &mutex_counter, NULL );
2965                iaddsharept = calloc( 1, sizeof(int) );
2966                *iaddsharept = 0;
2967
2968                for( i=0; i<nthread; i++ )
2969                {
2970                        targ[i].thread_no = i;
2971                        targ[i].follows = follows;
2972                        targ[i].njob = njob; 
2973                        targ[i].nadd = nadd; 
2974                        targ[i].nlen = nlen; 
2975                        targ[i].name = name;
2976                        targ[i].seq = seq;
2977                        targ[i].localhomtable = localhomtable;
2978                        targ[i].iscore = iscore;
2979                        targ[i].nscore = nscore;
2980                        targ[i].istherenewgap = istherenewgap;
2981                        targ[i].newgaplist = newgaplist_o;
2982                        targ[i].singlerna = singlerna;
2983                        targ[i].eff_kozo_mapped = eff_kozo_mapped;
2984                        targ[i].alloclen = alloclen;
2985                        targ[i].iaddshare = iaddsharept;
2986                        targ[i].dep = dep;
2987                        targ[i].topol = topol;
2988                        targ[i].len = len;
2989                        targ[i].addtree = addtree;
2990                        targ[i].mutex_counter = &mutex_counter;
2991                        pthread_create( handle+i, NULL, addsinglethread, (void *)(targ+i) );
2992                }
2993                for( i=0; i<nthread; i++ )
2994                {
2995                        pthread_join( handle[i], NULL );
2996                }
2997                pthread_mutex_destroy( &mutex_counter );
2998                free( handle );
2999                free( targ );
3000                free( iaddsharept );
3001        }
3002        else
3003#endif
3004        {
3005                thread_arg_t *targ;
3006                targ = calloc( 1, sizeof( thread_arg_t ) );
3007                targ[0].follows = follows;
3008                targ[0].njob = njob; 
3009                targ[0].nadd = nadd; 
3010                targ[0].nlen = nlen; 
3011                targ[0].name = name;
3012                targ[0].seq = seq;
3013                targ[0].localhomtable = localhomtable;
3014                targ[0].iscore = iscore;
3015                targ[0].nscore = nscore;
3016                targ[0].istherenewgap = istherenewgap;
3017                targ[0].newgaplist = newgaplist_o;
3018                targ[0].singlerna = singlerna;
3019                targ[0].eff_kozo_mapped = eff_kozo_mapped;
3020                targ[0].alloclen = alloclen;
3021                targ[0].dep = dep;
3022                targ[0].topol = topol;
3023                targ[0].len = len;
3024                targ[0].addtree = addtree;
3025                addsinglethread( targ );
3026                free( targ );
3027        }
3028        free( dep );
3029        FreeFloatMtx( len );
3030        if( multidist || tuplesize > 0 ) FreeFloatMtx( nscore );
3031
3032//      for( i=0; i<nadd; i++ ) fprintf( stdout, ">%s (%d) \n%s\n", name[norg+i], norg+i, seq[norg+i] );
3033
3034        if( treeout )
3035        {
3036                fp = fopen( "infile.tree", "a" );
3037                if( fp == 0 )
3038                {
3039                        fprintf( stderr, "File error!\n" );
3040                        exit( 1 );
3041                }
3042                for( i=0; i<nadd; i++ )
3043                {
3044                        fprintf( fp, "\n" );
3045                        fprintf( fp, "%8d: %s\n", norg+i+1, name[norg+i]+1 );
3046                        fprintf( fp, "          nearest sequence: %d\n", addtree[i].nearest + 1 );
3047                        fprintf( fp, "          approximate distance: %f\n", addtree[i].dist1 );
3048                        fprintf( fp, "          sister group: %s\n", addtree[i].neighbors );
3049                        fprintf( fp, "          approximate distance: %f\n", addtree[i].dist2 );
3050                        free( addtree[i].neighbors );
3051                }
3052                fclose( fp );
3053                free( addtree );
3054        }
3055
3056        for( iadd=0; iadd<nadd; iadd++ )
3057        {
3058                f = follows[iadd];
3059                for( i=0; follower[f][i]!=-1; i++ )
3060                        ;
3061                if( !(follower[f] = realloc( follower[f], (i+2)*sizeof(int) ) ) )
3062                {
3063                        fprintf( stderr, "Cannot reallocate follower[]" );
3064                        exit( 1 );
3065                }
3066                follower[f][i] = iadd;
3067                follower[f][i+1] = -1;
3068#if 0
3069                fprintf( stderr, "\nfollowers of %d = ", f );
3070                for( i=0; follower[f][i]!=-1; i++ )
3071                        fprintf( stderr, "%d ", follower[f][i]  );
3072                fprintf( stderr, "\n" );
3073#endif
3074        }
3075
3076        orderfp = fopen( "order", "w" );
3077        if( !orderfp )
3078        {
3079                fprintf( stderr, "Cannot open 'order'\n" );
3080                exit( 1 );
3081        }
3082        for( i=0; ordertable[i]!=-1; i++ )
3083        {
3084                fprintf( orderfp, "%d\n", ordertable[i] );
3085//              for( j=0; follower[i][j]!=-1; j++ )
3086//                      fprintf( orderfp, "%d\n", follower[i][j]+norg );
3087                for( j=0; follower[ordertable[i]][j]!=-1; j++ )
3088                        fprintf( orderfp, "%d\n", follower[ordertable[i]][j]+norg );
3089//                      fprintf( orderfp, "%d -> %d\n", follower[i][j]+norg, i );
3090        }
3091        fclose( orderfp );
3092
3093        posmap = AllocateIntVec( lenfull+2 );
3094        realign = calloc( lenfull+2, sizeof( Blocktorealign ) );
3095        for( i=0; i<lenfull+1; i++ ) posmap[i] = i;
3096        for( i=0; i<lenfull+1; i++ )
3097        {
3098                realign[i].nnewres = 0;
3099                realign[i].start = 0;
3100                realign[i].end = 0;
3101        }
3102
3103        fprintf( stderr, "\n\nCombining ..\n" );
3104        fflush( stderr );
3105        tmpseqlen = alloclen * 100;
3106        tmpseq = AllocateCharMtx( 1, tmpseqlen );
3107//      check1 = AllocateCharVec( tmpseqlen );
3108//      check2 = AllocateCharVec( tmpseqlen );
3109//      gappick0( check2, seq[0] );
3110        for( iadd=0; iadd<nadd; iadd++ )
3111        {
3112//              fprintf( stderr, "%d / %d\r", iadd, nadd );
3113                fflush( stderr );
3114
3115//              fprintf( stderr, "\niadd == %d\n", iadd );
3116                makegaplistcompact( lenfull, posmap, newgaplist_compact, newgaplist_o[iadd] );
3117                if( iadd == 0 || istherenewgap[iadd] )
3118                {
3119                        tmpseq1 = tmpseq[0];
3120//                      gaplist2alnx( lenfull, tmpseq1, seq[0], newgaplist_o[iadd], posmap, tmpseqlen );
3121                        gaplist2alnx( lenfull, tmpseq1, seq[0], newgaplist_compact, posmap, tmpseqlen );
3122//                      fprintf( stderr, "len = %d ? %d\n", strlen( tmpseq1 ), alloclen );
3123                        if( ( tmplen = strlen( tmpseq1 ) ) >= fullseqlen )
3124                        {
3125                                fullseqlen = tmplen * 2+1;
3126//                              fprintf( stderr, "Length over!\n" );
3127//                              fprintf( stderr, "strlen(tmpseq1)=%d\n", (int)strlen( tmpseq1 ) );
3128                                fprintf( stderr, "reallocating..." );
3129//                              fprintf( stderr, "alloclen=%d\n", alloclen );
3130//                              fprintf( stderr, "Please recompile!\n" );
3131//                              exit( 1 );
3132                                for( i=0; i<njob; i++ )
3133                                {
3134                                        seq[i] = realloc( seq[i], fullseqlen * sizeof( char ) );
3135                                        if( !seq[i] )
3136                                        {
3137                                                fprintf( stderr, "Cannot reallocate seq[][]\n" );
3138                                                exit( 1 );
3139                                        }
3140                                }
3141                                fprintf( stderr, "done.\n" );
3142                        }
3143                        strcpy( seq[0], tmpseq1 );
3144
3145                        ien = norg+iadd;
3146#ifdef enablemultithread
3147                        if( nthread > 0 && ien > 500 )
3148                        {
3149                                gaplist2alnxthread_arg_t *targ;
3150                                int jobpos;
3151                                pthread_t *handle;
3152                                pthread_mutex_t mutex;
3153                                fprintf( stderr, "%d / %d (threads %d-%d)\r", iadd, nadd, 0, nthread );
3154
3155                                targ = calloc( nthread, sizeof( gaplist2alnxthread_arg_t ) );
3156                                handle = calloc( nthread, sizeof( pthread_t ) );
3157                                pthread_mutex_init( &mutex, NULL );
3158                                jobpos = 1;
3159                                for( i=0; i<nthread; i++ )
3160                                {
3161//                                      targ[i].thread_no = i;
3162                                        targ[i].ncycle = ien;
3163                                        targ[i].jobpospt = &jobpos;
3164                                        targ[i].tmpseqlen = tmpseqlen;
3165                                        targ[i].lenfull = lenfull;
3166                                        targ[i].seq = seq;
3167//                                      targ[i].newgaplist = newgaplist_o[iadd];
3168                                        targ[i].newgaplist = newgaplist_compact;
3169                                        targ[i].posmap = posmap;
3170                                        targ[i].mutex = &mutex;
3171
3172                                        pthread_create( handle+i, NULL, gaplist2alnxthread, (void *)(targ+i) );
3173                                }
3174                                for( i=0; i<nthread; i++ )
3175                                {
3176                                        pthread_join( handle[i], NULL );
3177                                }
3178                                pthread_mutex_destroy( &mutex );
3179                                free( handle );
3180                                free( targ );
3181                        }
3182                        else
3183#endif
3184                        {
3185                                fprintf( stderr, "%d / %d\r", iadd, nadd );
3186                                for( i=1; i<ien; i++ )
3187                                {
3188                                        tmpseq1 = tmpseq[0];
3189                                        if( i == 1 ) fprintf( stderr, " %d / %d\r", iadd, nadd );
3190//                                      gaplist2alnx( lenfull, tmpseq1, seq[i], newgaplist_o[iadd], posmap, tmpseqlen  );
3191                                        gaplist2alnx( lenfull, tmpseq1, seq[i], newgaplist_compact, posmap, tmpseqlen  );
3192//                                      fprintf( stderr, ">%s (iadd=%d)\n%s\n", name[i], iadd, tmpseq1 );
3193                                        strcpy( seq[i], tmpseq1 );
3194                                }
3195                        }
3196                }
3197                tmpseq1 = tmpseq[0];
3198//              insertgapsbyotherfragments_simple( lenfull, tmpseq1, seq[norg+iadd], newgaplist_o[iadd], posmap );
3199                insertgapsbyotherfragments_compact( lenfull, tmpseq1, seq[norg+iadd], newgaplist_o[iadd], posmap );
3200//              fprintf( stderr, "%d = %s\n", iadd, tmpseq1 );
3201                eq2dash( tmpseq1 );
3202                strcpy( seq[norg+iadd], tmpseq1 );
3203
3204//              adjustposmap( lenfull, posmap, newgaplist_o[iadd] );
3205                adjustposmap( lenfull, posmap, newgaplist_compact );
3206                countnewres( lenfull, realign, posmap, newgaplist_o[iadd] ); // muda?
3207//              countnewres( lenfull, realign, posmap, newgaplist_compact ); // muda?
3208
3209        }
3210        fprintf( stderr, "\r   done.                      \n\n" );
3211
3212#if 0
3213        for( i=0; i<njob; i++ )
3214        {
3215                fprintf( stdout, ">%s\n", name[i] );
3216                fprintf( stdout, "%s\n", seq[i] );
3217        }
3218#endif
3219
3220#if 0
3221        fprintf( stderr, "realign[].nnewres = " );
3222        for( i=0; i<lenfull+1; i++ )
3223        {
3224                fprintf( stderr, "%d ", realign[i].nnewres );
3225        }
3226        fprintf( stderr, "\n" );
3227#endif
3228
3229        for( i=0; i<lenfull+1; i++ )
3230        {
3231        }
3232
3233        for( i=0; i<lenfull+1; i++ )
3234        {
3235                if( realign[i].nnewres > 1 ) 
3236                {
3237//                      fprintf( stderr, "i=%d: %d-%d\n", i, realign[i].start, realign[i].end );
3238                        fprintf( stderr, "\rRealigning %d/%d           \r", i, lenfull );
3239//                      zure = dorealignment_compact( realign+i, seq, &fullseqlen, norg );
3240//                      zure = dorealignment_order( realign+i, seq, &fullseqlen, norg, ordertable, follows );
3241                        zure = dorealignment_tree( realign+i, seq, &fullseqlen, norg, topol, follows );
3242#if 0
3243                        gappick0( check1, seq[0] );
3244                        fprintf( stderr, "check1 = %s\n", check1 );
3245                        if( strcmp( check1, check2 ) )
3246                        {
3247                                fprintf( stderr, "CHANGED!!!!!\n" );
3248                                exit( 1 );
3249                        }
3250#endif
3251                        for( j=i+1; j<lenfull+1; j++ )
3252                        {
3253                                if( realign[j].nnewres )
3254                                {
3255                                        realign[j].start -= zure;
3256                                        realign[j].end -= zure;
3257                                }
3258                        }
3259                }
3260        }
3261        FreeIntCub( topol );
3262        fprintf( stderr, "\r   done.                      \n\n" );
3263
3264        fflush( stderr );
3265
3266
3267        FreeIntMtx( newgaplist_o );
3268        FreeIntVec( newgaplist_compact );
3269        FreeIntVec( posmap );
3270        free( realign );
3271        free( istherenewgap );
3272        FreeIntMtx( follower );
3273        free( follows );
3274        free( ordertable );
3275        FreeCharMtx( tmpseq );
3276
3277
3278        writeData_pointer( prep_g, njob, name, nlen, seq );
3279#if 0
3280        writeData( stdout, njob, name, nlen, bseq );
3281        writePre( njob, name, nlen, bseq, !contin );
3282        writeData_pointer( prep_g, njob, name, nlen, aseq );
3283#endif
3284#if IODEBUG
3285        fprintf( stderr, "OSHIMAI\n" );
3286#endif
3287
3288#if SMALLMEMORY
3289        if( multidist )
3290        {
3291//              if( constraint ) FreeLocalHomTable_two( localhomtable, norg, nadd );
3292                if( constraint ) FreeLocalHomTable_one( localhomtable, norg, nadd );
3293        }
3294        else
3295#endif
3296        {
3297                if( constraint ) FreeLocalHomTable( localhomtable, njob );
3298        }
3299
3300        SHOWVERSION;
3301        return( 0 );
3302}
Note: See TracBrowser for help on using the repository browser.