source: tags/arb_5.1/GDE/MOLPHY/mltree.c

Last change on this file was 6143, checked in by westram, 15 years ago
  • backport [6141] (parts changing code, but only strings and comments)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.6 KB
Line 
1/*
2 * mltree.c   Adachi, J.   1996.03.01
3 * Copyright (C) 1992-1996 J. Adachi & M. Hasegawa, All rights reserved.
4 */
5
6#include "protml.h"
7
8#define SLSDEBUG 0
9#define MLSITE 0
10#define USERTREEDEBUG 0
11
12
13Node **
14new_npvector(n)
15int n;
16/* memory allocate a node pointer vector */
17{
18        Node **v;
19
20        v = (Node **) malloc((unsigned)n * sizeof(Node *));
21        if (v == NULL) maerror("in new_npvector().");
22        return v;
23}
24
25void
26free_npvector(v)
27Node **v;
28{
29        free((char *) v);
30}
31
32
33Tree *
34new_tree(maxspc, maxibrnch, numptrn, seqconint)
35int maxspc, maxibrnch, numptrn;
36imatrix seqconint;
37{
38        int n, i;
39        Tree *tr;
40        Node *dp, *up;
41
42        tr = (Tree *) malloc(sizeof(Tree));
43        if (tr == NULL) maerror("tr in new_tree().");
44        tr->ebrnchp = (Node **) malloc((unsigned)maxspc * sizeof(Node *));
45        if (tr->ebrnchp == NULL) maerror("ebrnchp in new_tree().");
46        tr->ibrnchp = (Node **) malloc((unsigned)maxibrnch * sizeof(Node *));
47        if (tr->ibrnchp == NULL) maerror("ibrnchp in new_tree().");
48        tr->bturn = new_ivector(maxspc);
49        for (n = 0; n < maxspc; n++) {
50                tr->bturn[n] = n;
51                dp = (Node *) malloc(sizeof(Node));
52                if (dp == NULL) maerror("dp in new_tree().");
53                up = (Node *) malloc(sizeof(Node));
54                if (up == NULL) maerror("up in new_tree().");
55                dp->isop = NULL;
56                up->isop = NULL;
57                dp->kinp = up;
58                up->kinp = dp;
59                dp->descen = TRUE;
60                up->descen = FALSE;
61                dp->num = n;
62                up->num = n;
63                dp->length = 0.0;
64                up->length = 0.0;
65                dp->lklhdl = 0.0;
66                up->lklhdl = 0.0;
67                dp->paths = new_ivector(maxspc);
68                up->paths = dp->paths;
69                for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
70                dp->paths[n] = 1;
71                dp->eprob = seqconint[n];
72                up->eprob = NULL;
73                dp->iprob = NULL;
74                up->iprob = new_dmatrix(numptrn, Tpmradix);
75                tr->ebrnchp[n] = dp;
76        }
77        for (n = 0; n < maxibrnch; n++) {
78                dp = (Node *) malloc(sizeof(Node));
79                if (dp == NULL) maerror("dp in new_tree().");
80                up = (Node *) malloc(sizeof(Node));
81                if (up == NULL) maerror("up in new_tree().");
82                dp->isop = NULL;
83                up->isop = NULL;
84                dp->kinp = up;
85                up->kinp = dp;
86                dp->descen = TRUE;
87                up->descen = FALSE;
88                dp->num = n + maxspc;
89                up->num = n + maxspc;
90                dp->length = 0.0;
91                up->length = 0.0;
92                dp->lklhdl = 0.0;
93                up->lklhdl = 0.0;
94                dp->paths = new_ivector(maxspc);
95                up->paths = dp->paths;
96                for (i = 0; i < maxspc; i++) dp->paths[i] = 0;
97                dp->eprob = NULL;
98                up->eprob = NULL;
99                dp->iprob = new_dmatrix(numptrn, Tpmradix);
100                up->iprob = new_dmatrix(numptrn, Tpmradix);
101                tr->ibrnchp[n] = dp;
102        }
103        tr->rootp = NULL;
104
105        return tr;
106} /*_ new_tree */
107
108
109int
110getbuftree(numspc, identif)
111int numspc;
112cmatrix identif;
113{
114        int i, buf;
115
116        buf = numspc * 3; /* - 3 */
117        for (i = 0; i < numspc; i++)
118                buf += strlen(identif[i]);
119        return buf;
120} /*_ getbuftree */
121
122
123void
124changedistan(distanmat, distanvec, numspc)
125dmatrix distanmat;
126dvector distanvec;
127int numspc;
128{
129        int i, j, k;
130
131        for (k = 0, i = 1; i < numspc; i++) {
132                for (j = 0; j < i; j++, k++)
133                        distanvec[k] = distanmat[i][j];
134        }
135}
136
137
138void
139getproportion(proportion, distanvec, maxpair)
140double *proportion;
141dvector distanvec;
142int maxpair;
143{
144        int i;
145        double maxdis;
146
147        maxdis = 0.0;
148        for (i = 0; i < maxpair; i++) {
149                if (distanvec[i] > maxdis) {
150                        maxdis = distanvec[i];
151                }
152        }
153        *proportion = (double)MAXCOLUMN / (maxdis * 3.0);
154        if (*proportion > 1.0) *proportion = 1.0;
155        if (Debug) printf("Proportion: %.5f   maxdis: %.5f\n", *proportion,maxdis);
156} /* getproportion */
157
158
159Infotree *
160newinfotrees(numtree, buftree)
161int numtree;
162int buftree;
163{
164        int i;
165        Infotree *info;
166        cvector m;
167
168        info = (Infotree *) malloc((unsigned)numtree * sizeof(Infotree));
169        if (info == NULL) maerror("in newinfotrees().");
170        m = (cvector) malloc((unsigned)(numtree * buftree) * sizeof(char));
171        if (m == NULL) maerror("2 in newinfotrees().");
172        for (i = 0; i < numtree; i++, m+=buftree)
173                info[i].ltplgy = m;
174
175        return info;
176} /*_ newinfotrees */
177
178
179Infoaltree *
180newinfoaltrees(numaltree, buftree)
181int numaltree;
182int buftree;
183{
184        int i;
185        Infoaltree *info;
186        cvector m;
187
188        info = (Infoaltree *) malloc((unsigned)numaltree * sizeof(Infoaltree));
189        if (info == NULL) maerror("in newinfoaltrees().");
190        m = (cvector) malloc((unsigned)(numaltree * buftree) * sizeof(char));
191        if (m == NULL) maerror("2 in newinfoaltrees().");
192        for (i = 0; i < numaltree; i++, m+=buftree)
193                info[i].ltplgy = m;
194
195        return info;
196} /*_ newinfoaltrees */
197
198
199void
200getnumtree(ifp, numtree)
201FILE *ifp;
202int *numtree;
203{
204        char *cp;
205        char line[BUFLINE];
206
207        while (fgets(line, BUFLINE, ifp) != NULL) {
208                /* fputs(line, ofp); */
209                if (sscanf(line, "%d", numtree) == 1) {
210                        if (Verbs_optn && User_optn && !Aprox_optn)
211                                fprintf(stderr,"%d trees\n", *numtree); /* stdout */
212                        return;
213                } else {
214                        if (*line == '#')
215                                continue;
216                        for (cp = line; isspace(*cp) && *cp != '\0'; cp++)
217                        if (*cp != '\0') {
218                                fputs(line, stderr);
219                                fprintf(stderr,
220                                        "\nCan't read number of user trees, bad format.\n");
221                                exit(1);
222                        } else {
223                                continue;
224                        }
225                }
226        }
227        fprintf(stderr, "\nCan't read number of user trees.\n");
228        exit(1);
229} /*_ getnumtree */
230
231
232void
233getusertree(ifp, strtree, buftree)
234FILE *ifp;
235cvector strtree;
236int buftree;
237{
238        char line[BUFLINE];
239        char *cp, *np, *xp, *ap, *bp;
240        int par, bra, nsp;
241        boolean word;
242
243        word = FALSE;
244        par = bra = nsp = 0;
245        np = strtree;
246        strtree[buftree - 1] = '\0';
247        while (fgets(line, BUFLINE, ifp) != NULL) {
248                if (!isalnum(*line)) word = FALSE;
249                /* printf("%s\n", line); */
250                if (*line == '#') continue; /* comment line */
251                for (cp = line; *cp != '\0'; ) {
252                        switch(*cp) {
253
254                        case ';':
255                                        if (np == strtree) goto next;
256                                        if (par == 0 && bra == 0) { /* tree */
257                                                *np = '\0';
258                                                if (nsp != Numspc) {
259                                                        fprintf(stderr,
260                                                                "ERROR: different number of OTU, %d/%d. \n",
261                                                                nsp, Numspc);
262                                                        fprintf(stderr, "%s\n", strtree);
263                                                        exit(1);
264                                                }
265                                                /* printf("%s\n", strtree); */
266                                                return;
267                                        } else { /* no tree */
268                                                if (par != 0) {
269                                                        fprintf(stderr,
270                                                                "ERROR, bad match parenthesis \"()\" !\n");
271                                                        if (par > 0)
272                                                                fprintf(stderr,"')' is %d less than '('.\n",
273                                                                        par);
274                                                        if (par < 0)
275                                                                fprintf(stderr,"')' is %d more than '('.\n",
276                                                                        abs(par));
277                                                }
278                                                if (bra != 0) {
279                                                        fprintf(stderr,
280                                                                "ERROR, bad match brace \"{}\" !\n");
281                                                        if (bra > 0)
282                                                                fprintf(stderr,"'}' is %d less than '{'.\n",
283                                                                        par);
284                                                        if (par < 0)
285                                                                fprintf(stderr,"'}' is %d more than '{'.\n",
286                                                                        abs(par));
287                                                }
288                                                fprintf(stderr, "%s\n", strtree);
289                                                exit(1);
290                                        }
291                                        break;
292                        case '(': *np++ = *cp++; par++; break;
293                        case ')': *np++ = *cp++; par--; break;
294                        case '{': *np++ = *cp++; bra++; break;
295                        case '}': *np++ = *cp++; bra--; break;
296                        default:
297                                if (isspace(*cp) || *cp == ',') {
298                                        cp++;
299                                } else { /* identifier */
300                                        if (isalnum(*cp) || *cp == '_' || *cp == '-') {
301                                                if (!word) {
302                                                        np--;
303                                                        if (isalnum(*np) || *np == '_' || *np == '-') {
304                                                                np++;
305                                                                *np = ' ';
306                                                        }
307                                                        np++;
308                                                }
309                                                ap = cp; bp = np;
310                                                while (isalnum(*cp) || *cp == '_' || *cp == '-') {
311                                                        *np++ = *cp++;
312                                                }
313                                                if (!word) nsp++;
314                                                *cp == '\0' ? (word = TRUE) : (word = FALSE);
315                                                if (Debug) printf("%3d %1d %-10.10s %-20.20s\n",
316                                                        nsp,word,bp,ap);
317                                        } else {
318                                                fprintf(stderr, "bad name: %s\n", cp);
319                                                exit(1);
320                                        }
321                                }
322
323                        } /* switch */
324                } /* while */
325                if (par == 0 && bra == 0 && np != strtree) {
326                        *np = '\0';
327                        /* printf("%s\n", strtree); */
328                        return;
329                }
330                if (np > (strtree + buftree - 2)) {
331                        fprintf(stderr, "Too long, users tree.\n");
332                        fprintf(stderr, "You may forget tree delimiter \";\".\n");
333                        fprintf(stderr, "\"%s\"\n", strtree);
334                        exit(1);
335                }
336#if USERTREEDEBUG
337                printf("%s\n", strtree);
338#endif
339                next: ;
340        }
341        fprintf(stderr, "Can't read users trees of input file.\n");
342        exit(1);
343} /*_ getusertree */
344
345
346Node *
347internalnode(tr, cpp, ninode, st)
348Tree *tr;
349char **cpp;
350int *ninode;
351char *st;
352{
353        Node *xp, *np, *rp;
354        int i, j, dvg;
355        char ident[MAXWORD];
356        char *idp;
357
358        /* fprintf(stderr, "in1:%s\n", *cpp); */
359
360        if (**cpp == ' ') (*cpp)++;
361        if (**cpp == '(') {
362                if (Debug) printf("internal1: %c\n", **cpp);
363                (*cpp)++;
364                xp = internalnode(tr, cpp, ninode, st);
365                xp->isop = xp;
366                dvg = 1;
367                while (**cpp != ')') {
368                        dvg++;
369                        np = internalnode(tr, cpp, ninode, st);
370                        np->isop = xp->isop;
371                        xp->isop = np;
372                        xp = np;
373                }
374                if (Debug) printf("internal2: %c\n", **cpp);
375                (*cpp)++;
376                if (dvg < 2) {
377                        fputs("ERROR, redundancy or unnecessary \"()\" !\n", stderr);
378                        fprintf(stderr, "%s\n", st);
379                        exit(1);
380                }
381                /* printf("ninode %d\n", *ninode + 1); */
382                rp = tr->ibrnchp[*ninode];
383                rp->isop = xp->isop;
384                xp->isop = rp;
385
386                for (j = 0; j < Numspc; j++) rp->paths[j] = 0;
387                for (xp = rp->isop; xp != rp; xp = xp->isop) {
388                        for (j = 0; j < Numspc; j++) {
389                                if (xp->paths[j] == 1) rp->paths[j] = 1;
390                        }
391                }
392                /*
393                if (Debug) {
394                        for (j = 0; j < Numspc; j++) printf("%2d",rp->paths[j]);
395                        putchar('\n');
396                }
397                */
398                (*ninode)++;
399                if (*ninode > Maxibrnch) {
400                        fputs("ERROR, too many number of internal branch!\n", stderr);
401                        fputs("This tree may be rooted tree.\n", stderr);
402                        fprintf(stderr, "%s\n", st);
403                        exit(1);
404                }
405
406                /* fprintf(stderr, "in2:%s\n", *cpp); */
407
408                return rp->kinp;
409        } else if (isalnum(**cpp)) {
410                if (Debug) printf("external: %c\n", **cpp);
411                for (idp = ident; **cpp!=' ' && **cpp!='(' && **cpp!=')'; (*cpp)++) {
412                        *idp++ = **cpp;
413                        if (Debug) putchar(**cpp);
414                }
415                *idp = '\0';
416                if (Debug) putchar('\n');
417
418                for (i = 0; i < Numspc; i++) {
419                        /* puts(Identif[i]); */
420                        if (!strcmp(ident, Identif[i])) {
421                                return tr->ebrnchp[i]->kinp;
422                        }
423                }
424                fprintf(stderr, "ERROR, abnormal identifier(OTU name): %s !\n", ident);
425                exit(1);
426        } else {
427                fprintf(stderr, "ERROR, abnormal tree topology character: %s\n", *cpp);
428                exit(1);
429        }
430        return NULL;
431} /*_ internalnode */
432
433
434void
435constructtree(tr, strtree)
436Tree *tr;
437cvector strtree;
438{
439        char *cp;
440        int ninode, dvg;
441        Node *xp, *np;
442
443#if USERTREEDEBUG
444        printf("%s\n", strtree);
445#endif
446        ninode = 0;
447        cp = strtree;
448        if (*cp == '(') {
449                if (Debug) printf("roottre0: %c\n", *cp);
450                cp++;
451                xp = internalnode(tr, &cp, &ninode, strtree);
452                xp->isop = xp;
453                dvg = 1;
454                while (*cp != ')') {
455                        dvg++;
456                        np = internalnode(tr, &cp, &ninode, strtree);
457                        np->isop = xp->isop;
458                        xp->isop = np;
459                        xp = np;
460                }
461                if (Debug) printf("roottre0: %c\n", *cp);
462                if (dvg < 3) {
463                        fputs("ERROR: this tree is not unroot tree!\n", stderr);
464                        fprintf(stderr, "%s\n", strtree);
465                        exit(1);
466                }
467                tr->rootp = xp;
468                Numibrnch = ninode;
469                Numbrnch = Numspc + ninode;
470        } else {
471                fprintf(stderr, "ERROR users tree:\n");
472                fprintf(stderr, "%s\n", strtree);
473                exit(1);
474        }
475} /*_ constructtree */
476
477
478void
479prcurtree(tr)
480Tree *tr;
481{
482        Node *dp, *up, *cp;
483        int i, j;
484        double sum;
485
486/*      printf("\nStructure of Tree\n"); */
487        printf("\n%4s%5s%5s%5s%7s%7s%37s\n",
488        "num","kinp","isop","isop","descen","length", "namesp");
489        for (i = 0; i < Numspc; i++) {
490                dp = tr->ebrnchp[i];
491                up = dp->kinp;
492                printf("%4d", dp->num+1);
493                printf("%5d", up->num+1);
494                if (dp->isop == NULL) printf("%5s", "null");
495                else                  printf("%5d", dp->isop->num+1);
496                if (up->isop == NULL) printf("%5s", "null");
497                else                  printf("%5d", up->isop->num+1);
498                printf("%3d", dp->descen);
499                printf("%3d", up->descen);
500                printf("%8.3f", dp->length);
501                printf("%8.3f", up->length);
502                printf("%5d%5d", dp->eprob[0],dp->eprob[Numptrn-1]);
503                for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[0][j];
504                printf("%5.0f", sum * 1000.0);
505                for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[Numptrn-1][j];
506                printf("%5.0f", sum * 1000.0);
507                printf("   %s", Identif[dp->num]);
508                putchar('\n');
509        }
510        for (i = 0; i < Numibrnch; i++) {
511                dp = tr->ibrnchp[i];
512                up = dp->kinp;
513                printf("%4d", dp->num+1);
514                printf("%5d", up->num+1);
515                if (dp->isop == NULL) printf("%5s", "null");
516                else                  printf("%5d", dp->isop->num+1);
517                if (up->isop == NULL) printf("%5s", "null");
518                else                  printf("%5d", up->isop->num+1);
519                printf("%3d", dp->descen);
520                printf("%3d", up->descen);
521                printf("%8.3f", dp->length);
522                printf("%8.3f", up->length);
523                for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += dp->iprob[0][j];
524                printf("%5.0f", sum * 1000.0);
525                for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += dp->iprob[Numptrn-1][j];
526                printf("%5.0f", sum * 1000.0);
527                for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[0][j];
528                printf("%5.0f", sum * 1000.0);
529                for (sum = 0.0, j = 0; j < Tpmradix; j++) sum += up->iprob[Numptrn-1][j];
530                printf("%5.0f", sum * 1000.0);
531                for (cp = dp->isop; cp != dp; cp = cp->isop) {
532                        printf("%5d", cp->num+1);
533                }
534                putchar('\n');
535        }
536                dp = tr->rootp->isop;
537                printf("%4d", dp->num+1);
538                for (cp = dp->isop; cp != dp; cp = cp->isop) {
539                        printf("%5d", cp->num+1);
540                }
541                putchar('\n');
542} /*_ prcurtree */
543
544
545void
546pathing(tr)
547Tree *tr;
548{
549        int i, j;
550        ivector pthsd, pthsa;
551        Node *rp, *cp, *xp, *dp;
552
553        cp = rp = tr->rootp;
554        do {
555                cp = cp->isop->kinp;
556                if (cp->isop == NULL) { /* external node */
557                        /*
558                        pthsd = cp->paths;
559                        */
560                        cp = cp->kinp;
561                        /*
562                        pthsa = cp->paths;
563                        for (i = 0; i < Numspc; i++) {
564                                pthsd[i] == 1 ? (pthsa[i] = 0) : (pthsa[i] = 1);
565                        }
566                        */
567                } else { /* internal node */
568                        if (!cp->descen) { /* ascent */
569                                dp = cp->kinp;
570                                pthsd = dp->paths;
571                                for (i = 0; i < Numspc; i++) pthsd[i] = 0;
572                                for (xp = dp->isop; xp != dp; xp = xp->isop) {
573                                        pthsa = xp->kinp->paths;
574                                        for (i = 0; i < Numspc; i++) {
575                                                if (pthsa[i]) pthsd[i] = 1;
576                                        }
577                                }
578                                /*
579                                pthsa = cp->paths;
580                                for (i = 0; i < Numspc; i++) {
581                                        pthsd[i] == 1 ? (pthsa[i] = 0) : (pthsa[i] = 1);
582                                }
583                                */
584                        }
585                }
586        } while (cp != rp);
587/* */
588        if (Debug_optn) {
589        puts("\npath of internal nodes");
590        for (i = 0; i < Maxibrnch; i++) {
591                pthsd = tr->ibrnchp[i]->paths;
592                for (j = 0; j < Numspc; j++) printf("%2d", pthsd[j]); putchar('\n');
593        }
594        }
595/* */
596} /*_ pathing */
597
598
599#if 0
600static void
601leastsquares(am, ym)
602dmatrix am;
603dvector ym;
604{
605        int i, j, k;
606        double pivot, element;
607        dmatrix im;
608
609        im = new_dmatrix(Numbrnch, Numbrnch);
610        for (i = 0; i < Numbrnch; i++) {
611                for (j = 0; j < Numbrnch; j++)
612                        im[i][j] = 0.0;
613                im[i][i] = 1.0;
614        }
615        for (k = 0; k < Numbrnch; k++) {
616                pivot = am[k][k];
617                ym[k] /= pivot;
618                for (j = 0; j < Numbrnch; j++) {
619                        am[k][j] /= pivot;
620                        im[k][j] /= pivot;
621                }
622                for (i = 0; i < Numbrnch; i++) {
623                        if (k != i) {
624                                element = am[i][k];
625                                ym[i] -= element * ym[k];
626                                for (j = 0; j < Numbrnch; j++) {
627                                        am[i][j] -= element * am[k][j];
628                                        im[i][j] -= element * im[k][j];
629                                }
630                        }
631                }
632        }
633/*
634        putchar('\n');
635        for (i = 0; i < Numbrnch; i++) {
636                for (j = 0; j < Numbrnch; j++)
637                        printf("%10.6f", im[i][j]);
638                putchar('\n');
639        }
640*/
641        free_dmatrix(im);
642} /*_ leastsquares */
643#endif
644
645
646static void
647luequation(amat, yvec, size)
648dmatrix amat;
649dvector yvec;
650int size;
651{
652        /* SOLVE THE LINEAR SET OF EQUATIONS ON LU DECOMPOSITION */
653    double eps = 1.0e-20; /* ! */
654        int i, j, k, l, maxi, idx;
655        double sum, tmp, maxb, aw;
656        dvector wk;
657        ivector index;
658
659        wk = new_dvector(size);
660        index = new_ivector(size);
661        aw = 1.0;
662        for (i = 0; i < size; i++) {
663                maxb = 0.0;
664                for (j = 0; j < size; j++) {
665                        if (fabs(amat[i][j]) > maxb)
666                                maxb = fabs(amat[i][j]);
667                }
668                if (maxb == 0.0) {
669                        fprintf(stderr, "luequation: singular matrix\n");
670                        exit(1);
671                }
672                wk[i] = 1.0 / maxb;
673        }
674        for (j = 0; j < size; j++) {
675                for (i = 0; i < j; i++) {
676                        sum = amat[i][j];
677                        for (k = 0; k < i; k++)
678                                sum -= amat[i][k] * amat[k][j];
679                        amat[i][j] = sum;
680                }
681                maxb = 0.0;
682                for (i = j; i < size; i++) {
683                        sum = amat[i][j];
684                        for (k = 0; k < j; k++)
685                                sum -= amat[i][k] * amat[k][j];
686                        amat[i][j] = sum;
687                        tmp = wk[i] * fabs(sum);
688                        if (tmp >= maxb) {
689                                maxb = tmp;
690                                maxi = i;
691                        }
692                }
693                if (j != maxi) {
694                        for (k = 0; k < size; k++) {
695                                tmp = amat[maxi][k];
696                                amat[maxi][k] = amat[j][k];
697                                amat[j][k] = tmp;
698                        }
699                        aw = -aw;
700                        wk[maxi] = wk[j];
701                }
702                index[j] = maxi;
703                if (amat[j][j] == 0.0)
704                        amat[j][j] = eps;
705                if (j != size - 1) {
706                        tmp = 1.0 / amat[j][j];
707                        for (i = j + 1; i < size; i++)
708                                amat[i][j] *= tmp;
709                }
710        }
711
712        l = -1;
713        for (i = 0; i < size; i++) {
714                idx = index[i];
715                sum = yvec[idx];
716                yvec[idx] = yvec[i];
717                if (l != -1) {
718                        for (j = l; j < i; j++)
719                                sum -= amat[i][j] * yvec[j];
720                } else if (sum != 0.0)
721                        l = i;
722                yvec[i] = sum;
723        }
724        for (i = size - 1; i >= 0; i--) {
725                sum = yvec[i];
726                for (j = i + 1; j < size; j++)
727                        sum -= amat[i][j] * yvec[j];
728                yvec[i] = sum / amat[i][i];
729        }
730        free_ivector(index);
731        free_dvector(wk);
732} /*_ luequation */
733
734
735void
736lslength(tr, distanvec, numspc)
737Tree *tr;
738dvector distanvec;
739int numspc;
740{
741        int i, i1, j, j1, j2, k, numibrnch, numbrnch, numpair;
742        double sum, leng, alllen, rss;
743        ivector pths;
744        dmatrix atmt, atamt;
745        Node **ebp, **ibp;
746
747        if (Debug) printf("numspc = %d\n", numspc);
748        numibrnch = Numibrnch;
749        numbrnch = numspc + numibrnch;
750        numpair = (numspc * (numspc - 1)) / 2;
751        atmt = new_dmatrix(numbrnch, numpair);
752        atamt = new_dmatrix(numbrnch, numbrnch);
753        ebp = tr->ebrnchp;
754        ibp = tr->ibrnchp;
755        if (Debug) {
756                putchar('\n');
757                for (i = 0; i < numspc; i++) {
758                        for (j = 0; j < numspc; j++) printf("%2d",ebp[i]->paths[j]);
759                        putchar('\n');
760                }
761                for (i = 0; i < numibrnch; i++) {
762                        for (j = 0; j < numspc; j++) printf("%2d",ibp[i]->paths[j]);
763                        putchar('\n');
764                }
765        }
766
767        for (i = 0; i < numspc; i++) {
768                for (j1 = 1, j = 0; j1 < numspc; j1++) {
769                        if (j1 == i) {
770                                for (j2 = 0; j2 < j1; j2++, j++) {
771                                        atmt[i][j] = 1.0;
772                                }
773                        } else {
774                                for (j2 = 0; j2 < j1; j2++, j++) {
775                                        if (j2 == i)
776                                                atmt[i][j] = 1.0;
777                                        else
778                                                atmt[i][j] = 0.0;
779                                }
780                        }
781                }
782        }
783        for (i1 = 0, i = numspc; i1 < numibrnch; i1++, i++) {
784                pths = ibp[i1]->paths;
785                for (j1 = 1, j = 0; j1 < numspc; j1++) {
786                        for (j2 = 0; j2 < j1; j2++, j++) {
787                                if (pths[j1] != pths[j2])
788                                        atmt[i][j] = 1.0;
789                                else
790                                        atmt[i][j] = 0.0;
791                        }
792                }
793        }
794
795        if (Debug) {
796                putchar('\n');
797                for (i = 0; i < numpair; i++) {
798                        for (j = 0; j < numbrnch; j++)
799                                printf("%2.0f",atmt[j][i]); /* !? */
800                        printf("%6.1f\n",distanvec[i]);
801                }
802        }
803#ifdef DIST
804                putchar('\n');
805                for (i = 0; i < numpair; i++)
806                        printf("%5.1f",distanvec[i]);
807                putchar('\n');
808#endif
809        for (i = 0; i < numbrnch; i++) {
810                for (j = 0; j <= i; j++) {
811                        for (k = 0, sum = 0.0; k < numpair; k++)
812                                sum += atmt[i][k] * atmt[j][k];
813                        atamt[i][j] = sum;
814                        atamt[j][i] = sum;
815                }
816        }
817        for (i = 0; i < numbrnch; i++) {
818                for (k = 0, sum = 0.0; k < numpair; k++)
819                        sum += atmt[i][k] * distanvec[k];
820                Brnlength[i] = sum;
821        }
822#if 0
823        putchar('\n');
824        for (i = 0; i < numbrnch; i++) {
825                for (j = 0; j < numbrnch; j++)
826                        printf("%5.1f", atamt[i][j]);
827                printf("%7.3f",Brnlength[i]);
828                putchar('\n');
829        }
830#endif
831        luequation(atamt, Brnlength, numbrnch);
832#if 0
833        putchar('\n');
834        for (i = 0; i < numbrnch; i++) {
835                for (j = 0; j < numbrnch; j++)
836                        printf("%5.1f", atamt[i][j]);
837                printf("%7.3f",Brnlength[i]);
838                putchar('\n');
839        }
840#endif
841        if (Debug) {
842                for (i = 0; i < numbrnch; i++) {
843                        printf("%7.3f",Brnlength[i]);
844                } putchar('\n');
845        }
846        for (i = 0, rss = 0.0; i < numpair; i++) {
847                sum = distanvec[i];
848                for (j = 0; j < numbrnch; j++) {
849                        if (atmt[j][i] == 1.0 && Brnlength[j] > 0.0) /* caution */
850                                sum -= Brnlength[j];
851                }
852                rss += sum * sum;
853        }
854        tr->rssleast = sqrt(rss); /* / numpair !? */
855        alllen = 0.0;
856        for (i = 0; i < numspc; i++) {
857                leng = Brnlength[i];
858        /*      if (leng > 0.0) leng = 0.0   caution */
859                alllen += leng;
860                if (leng < Llimit) leng = Llimit;
861                if (leng > Ulimit) leng = Ulimit;
862                ebp[i]->length = leng;
863                ebp[i]->kinp->length = leng;
864                Brnlength[i] = leng;
865        }
866        for (i = 0, j = numspc; i < numibrnch; i++, j++) {
867                leng = Brnlength[j];
868        /*      if (leng > 0.0) leng = 0.0   caution */
869                alllen += leng;
870                if (leng < Llimit) leng = Llimit;
871                if (leng > Ulimit) leng = Ulimit;
872                ibp[i]->length = leng;
873                ibp[i]->kinp->length = leng;
874                Brnlength[j] = leng;
875        }
876        tr->tbldis = alllen;
877        free_dmatrix(atmt);
878        free_dmatrix(atamt);
879} /*_ lslength */
880
881
882void
883slslength(tr, dmat, ns) /* simplify least square method (not equal LS) */
884Tree *tr;
885dmatrix dmat;
886int ns;
887{
888        int i, j, k, m, ne, ni, l, r, ll, rr;
889        int numibrnch, numbrnch;
890        double sumc, suml, sumr, leng, alllen, rss, coef;
891        ivector pths, lpaths, rpaths, npaths;
892        Node **ebp, **ibp, *dp, *ap, *xp;
893
894        lpaths = new_ivector(ns);
895        rpaths = new_ivector(ns);
896        npaths = new_ivector(ns);
897        numibrnch = Numibrnch;
898        numbrnch = ns + numibrnch;
899        ebp = tr->ebrnchp;
900        ibp = tr->ibrnchp;
901        alllen = 0.0;
902
903#if SLSDEBUG
904        for (i = 0; i < ns; i++) {
905                for (j = 0; j < ns; j++) printf("%6.0f", dmat[i][j]*1000);
906                putchar('\n');
907        } putchar('\n');
908        for (i = 0; i < ns; i++) {
909                for (j = 0, sumc = 0.0; j < ns; j++) sumc += dmat[i][j];
910                printf("%6.0f", sumc*1000);
911        } putchar('\n'); putchar('\n');
912#endif /* SLSDEBUG */
913
914        for ( ne = 0; ne < numbrnch; ne++) {
915
916                for (k = 0; k < ns; k++) {
917                        npaths[k] = rpaths[k] = 0;
918                }
919                if (ne < ns) { /* internal */
920                        l = 1;
921                        ll = 1;
922                        dp = ebp[ne];
923                        ap = ebp[ne]->kinp;
924                        pths = dp->paths;
925                        for (k = 0; k < ns; k++) lpaths[k] = pths[k];
926                        npaths[ne] = 1;
927                } else { /* external */
928                        ni = ne - ns;
929                        dp = ibp[ni];
930                        ap = ibp[ni]->kinp;
931                        pths = dp->paths;
932                        for (k = 0, l = 0; k < ns; k++) {
933                                lpaths[k] = 0;
934                                if (pths[k]) l++;
935                        }
936                        for (xp = dp->isop, j = 1; xp != dp; xp = xp->isop, j++) {
937                                pths = xp->paths;
938                                for (k = 0, m = 0; k < ns; k++) if (pths[k]) m++;
939                                if (xp->descen) {
940                                        m = ns - m;
941                                        for (k = 0; k < ns; k++) {
942                                                if (!pths[k]) {
943                                                        lpaths[k] = j;
944                                                        npaths[k] = m;
945                                                }
946                                        }
947                                } else {
948                                        for (k = 0; k < ns; k++) {
949                                                if (pths[k]) {
950                                                        lpaths[k] = j;
951                                                        npaths[k] = m;
952                                                }
953                                        }
954                                }
955                        }
956                        ll = j - 1;
957                }
958                r = ns - l;
959                for (xp = ap->isop, j = 1; xp != ap; xp = xp->isop, j++) {
960                        pths = xp->paths;
961                        for (k = 0, m = 0; k < ns; k++) if (pths[k]) m++;
962                        if (xp->descen) {
963                                m = ns - m;
964                                for (k = 0; k < ns; k++) {
965                                        if (!pths[k]) {
966                                                rpaths[k] = j;
967                                                npaths[k] = m;
968                                        }
969                                }
970                        } else {
971                                for (k = 0; k < ns; k++) {
972                                        if (pths[k]) {
973                                                rpaths[k] = j;
974                                                npaths[k] = m;
975                                        }
976                                }
977                        }
978                }
979                rr = j - 1;
980
981#if SLSDEBUG
982                for (k = 0; k < ns; k++) printf("%2d", lpaths[k]); putchar('\n');
983                for (k = 0; k < ns; k++) printf("%2d", rpaths[k]); putchar('\n');
984                for (k = 0; k < ns; k++) printf("%2d", npaths[k]); putchar('\n');
985#endif /* SLSDEBUG */
986
987                sumc = suml = sumr = 0.0;
988                for (i = 0; i < ns - 1; i++) {
989                        for (j = i + 1; j < ns; j++) {
990                                coef = (double)npaths[i] * (double)npaths[j];
991                                if (lpaths[i]) {
992                                        if (rpaths[j]) {
993                                                sumc += dmat[i][j] / coef;
994#if                                             SLSDEBUG
995                                                printf("c1:%3d%3d%8.0f%7.0f\n",
996                                                        i+1,j+1,sumc*1000,dmat[i][j]/coef*1000);
997#endif                                  /* SLSDEBUG */
998                                        } else if (lpaths[i] != lpaths[j]) {
999                                                suml += dmat[i][j] / coef;
1000#if                                             SLSDEBUG
1001                                                printf("l :%3d%3d%8.0f%7.0f\n",
1002                                                        i+1,j+1,suml*1000,dmat[i][j]/coef*1000);
1003#endif                                  /* SLSDEBUG */
1004                                        }
1005                                } else {
1006                                        if (!rpaths[j]) {
1007                                                sumc += dmat[i][j] / coef;
1008#if                                             SLSDEBUG
1009                                                printf("c2:%3d%3d%8.0f%7.0f\n",
1010                                                        i+1,j+1,sumc*1000,dmat[i][j]/coef*1000);
1011#endif                                  /* SLSDEBUG */
1012                                        } else if (rpaths[i] != rpaths[j]) {
1013                                                sumr += dmat[i][j] / coef;
1014#if                                             SLSDEBUG
1015                                                printf("r :%3d%3d%8.0f%7.0f\n",
1016                                                        i+1,j+1,sumr*1000,dmat[i][j]/coef*1000);
1017#endif                                  /* SLSDEBUG */
1018                                        }
1019                                }
1020                        }
1021                }
1022
1023#if             SLSDEBUG
1024                printf("%3s:%3d",(ne < ns ? "ext" : "int"), ne+1);
1025                printf(" %3d%3d%3d%3d", l, r, ll, rr);
1026                printf(" %5.2f%5.2f", (ne < ns ? 0 : (double)rr / (double)(ll-1)),
1027                        (double)ll / (double)(rr-1));
1028                printf("%7.0f%7.0f%7.0f", sumc*1000, suml*1000, sumr*1000);
1029#endif  /* SLSDEBUG */
1030                if (ne >= ns) suml *= (double)rr / (double)(ll-1);
1031                sumr *= (double)ll / (double)(rr-1);
1032                leng = ( sumc - sumr - suml ) / (double)(ll * rr);
1033#if             SLSDEBUG
1034                printf("   leng:%9.5f\n", leng);
1035#endif  /* SLSDEBUG */
1036
1037                alllen += leng;
1038
1039                if (leng < 0.0) leng = 0.0;  /* caution!? */
1040                if (ne < ns) /* external */
1041                        pths = ebp[ne]->paths;
1042                else /* internal */
1043                        pths = ibp[ni]->paths;
1044                for (i = 1; i < ns ; i++) {
1045                        for (j = 0; j < i; j++) {
1046                                if (pths[i] != pths[j]) dmat[i][j] -= leng;
1047                        }
1048                }
1049
1050                if (leng > Ulimit) leng = Ulimit;
1051                if (ne < ns)  { /* external */
1052                        if (leng < LOWERLIMIT) leng = LOWERLIMIT;
1053                        ebp[ne]->length = ebp[ne]->kinp->length = leng;
1054                } else { /* internal */
1055                        if (leng < Llimit) leng = Llimit;
1056                        ibp[ni]->length = ibp[ni]->kinp->length = leng;
1057                }
1058                Brnlength[ne] = leng; /* !? */
1059
1060        } /* for ne */
1061
1062        for (i = 1, rss = 0.0; i < ns ; i++) {
1063                for (j = 0; j < i; j++) {
1064                        rss += dmat[i][j] * dmat[i][j];
1065#if                     SLSDEBUG
1066                        printf("%3d%3d%9.3f\n", i+1, j+1, dmat[i][j]);
1067#endif          /* SLSDEBUG */
1068                }
1069        }
1070        tr->rssleast = sqrt(rss); /* / numpair !? */
1071#if SLSDEBUG
1072        printf("rss: %9.3f%9.3f\n", rss, tr->rssleast);
1073#endif /* SLSDEBUG */
1074        for (i = 1; i < ns ; i++) {
1075                for (j = 0; j < i; j++) dmat[i][j] = dmat[j][i];
1076        }
1077
1078        tr->tbldis = alllen;
1079        free_ivector(lpaths);
1080        free_ivector(rpaths);
1081        free_ivector(npaths);
1082} /*_ slslength */
1083
1084
1085#define DEBUGFM 0
1086#define LSLEN   0
1087void
1088fmlength(tr, dmat, ns)
1089Tree *tr;
1090dmatrix dmat;
1091int ns;
1092{
1093        int i, j, k, no, nx, ny, m;
1094        double sy, x, xav, alllen, rss;
1095        ivector otux, path;
1096        dvector xv, yv;
1097        dmatrix dism;
1098        Node **otup, *cp, *rp, *dp, *xp;
1099#if     LSLEN
1100        ivector otun;
1101#endif  /* LSLEN */
1102
1103        dism = new_dmatrix(ns, ns);
1104        for (i = 0; i < ns; i++) {
1105                for (j = 0; j < ns; j++) dism[i][j] = dmat[i][j];
1106        }
1107        otux = new_ivector(ns);
1108        xv = new_dvector(ns);
1109        yv = new_dvector(ns);
1110        otup = (Node **)new_npvector(ns);
1111#if     LSLEN
1112        otun = new_ivector(ns);
1113#endif  /* LSLEN */
1114        for (i = 0; i < ns; i++) {
1115                otux[i] = 0;
1116                otup[i] = tr->ebrnchp[i]->kinp;
1117#if             LSLEN
1118                otun[i] = 1;
1119#endif  /* LSLEN */
1120        }
1121        alllen = 0.0;
1122
1123        no = ns;
1124        cp = rp = tr->rootp;
1125        do {
1126                cp = cp->isop->kinp;
1127                if (cp->isop == NULL) { /* external node */
1128                        cp = cp->kinp; /* not descen */
1129                } else { /* internal node */
1130                        if (!cp->descen) {
1131                                dp = cp->kinp;
1132                                m = ns;
1133                                for (xp = dp->isop, nx = 0; xp != dp; xp = xp->isop, nx++) {
1134                                        i = xp->num;
1135                                        otux[i] = i;
1136                                        if (i < m) m = i;
1137                                }
1138                                ny = no - nx;
1139                                for (xp = dp->isop, sy = 0; xp != dp; xp = xp->isop) {
1140                                        i = xp->num;
1141                                        for (j = 0, xv[i] = 0, yv[i] = 0; j < ns; j++) {
1142                                                if (otup[j] != NULL) {
1143                                                        if (otux[j])
1144                                                                xv[i] += dism[i][j];
1145                                                        else
1146                                                                yv[i] += dism[i][j];
1147                                                }
1148                                        }
1149                                        yv[i] /= ny;
1150                                        sy += yv[i];
1151                                }
1152#if                             DEBUGFM
1153                                printf("\nint. baranch %3d", cp->num+1);
1154                                printf("   m:%3d  nx:%3d  ny:%3d  no:%3d\n", m+1, nx, ny, no);
1155                                printf("%3s", "");
1156                                for (j = 0; j < ns; j++) printf("%4d", j+1); putchar('\n');
1157                                for (i = 0; i < ns; i++) {
1158                                        printf("%3d", i+1);
1159                                        for (j = 0; j < ns; j++) printf("%4.0f", dism[i][j] * 10);
1160                                        putchar('\n');
1161                                }
1162                                printf("\n%3s", "otp");
1163                                for (j = 0; j < ns; j++) printf("%4d", otup[j] != NULL ? 1 : 0);
1164                                putchar('\n');
1165                                printf("%3s", "otx");
1166                                for (j = 0; j < ns; j++) printf("%4d", otux[j]); putchar('\n');
1167                                printf("%3s", "xv");
1168                                for (j = 0; j < ns; j++) printf("%4.0f", xv[j]*10);
1169                                putchar('\n');
1170                                printf("%3s", "yv");
1171                                for (j = 0; j < ns; j++) printf("%4.0f", yv[j]*10);
1172                                putchar('\n');
1173#endif                  /* DEBUGFM */
1174                                for (xp = dp->isop, xav = 0; xp != dp; xp = xp->isop) {
1175                                        i = xp->num;
1176                                        x = (xv[i]+yv[i])/nx - (sy-yv[i])/(nx*(nx-1));
1177                                        alllen += x; /* caution!? */
1178                                        if (x < Llimit) x = Llimit;
1179                                        xp->length = xp->kinp->length = x;
1180                                        xav += x;
1181                                        Brnlength[xp->kinp->num] = x;
1182#if                                     DEBUGFM
1183                                        printf("b-len %3d%9.4f%9.4f%9.4f\n", xp->kinp->num+1
1184                                                x, (xv[i]+yv[i])/nx, (sy-yv[i])/(nx*(nx-1)));
1185#endif                          /* DEBUGFM */
1186                                        if (i != m) {
1187                                                for (j = 0; j < ns; j++) {
1188                                                        if (otup[j] != NULL && j != m) {
1189                                                                dism[m][j] += dism[i][j];
1190                                                                dism[j][m] += dism[j][i];
1191                                                        }
1192#if                                                     DEBUGFM
1193                                                        dism[i][j] = dism[j][i] = 0.0;
1194#endif                                          /* DEBUGFM */
1195                                                }
1196                                                otup[i] = NULL;
1197                                        }
1198                                        otux[i] = 0;
1199                                        if (xp->kinp->isop != NULL) xp->num = xp->kinp->num;
1200                                }
1201
1202                                for (j = 0; j < ns; j++) {
1203                                        if (otup[j] != NULL && j != m) {
1204                                                dism[m][j] = (dism[m][j] - xav)/nx;
1205                                                dism[j][m] = (dism[j][m] - xav)/nx;
1206                                        }
1207                                        xv[j] = yv[j] = 0.0;
1208                                }
1209                                cp->num = m;
1210                                no -= (nx - 1);
1211
1212                        }
1213                }
1214        } while (cp != rp);
1215
1216        xp = rp; nx = 0;
1217        do {
1218                xp = xp->isop;
1219                i = xp->num;
1220                otux[i] = 1;
1221                nx++;
1222        } while (xp != rp);
1223        ny = no - nx;
1224        xp = rp; sy = 0.0;
1225        do {
1226                xp = xp->isop;
1227                i = xp->num;
1228                for (j = 0, xv[i] = 0; j < ns; j++) {
1229                        if (otup[j] != NULL && otux[j]) xv[i] += dism[i][j];
1230                }
1231                sy += xv[i];
1232        } while (xp != rp);
1233        sy /= 2;
1234#if     DEBUGFM
1235        printf("root   nx:%3d  ny:%3d  no:%3d\n", nx, ny, no);
1236        printf("%3s", "");
1237        for (j = 0; j < ns; j++) printf("%4d", j+1); putchar('\n');
1238        for (i = 0; i < ns; i++) {
1239                printf("%3d", i+1);
1240                for (j = 0; j < ns; j++) printf("%4.0f", dism[i][j] * 10);
1241                putchar('\n');
1242        }
1243        printf("\n%3s", "otp");
1244        for (j = 0; j < ns; j++) printf("%4d", otup[j] != NULL ? 1 : 0);
1245        putchar('\n');
1246        printf("%3s", "otx");
1247        for (j = 0; j < ns; j++) printf("%4d", otux[j]); putchar('\n');
1248        printf("%3s", "xv");
1249        for (j = 0; j < ns; j++) printf("%4.0f", xv[j] * 10); putchar('\n');
1250        printf("sy:%9.5f\n", sy);
1251#endif /* DEBUGFM */
1252        xp = rp;
1253        do {
1254                xp = xp->isop;
1255                i = xp->num;
1256                x = ( xv[i] - (sy-xv[i])/(nx-2) ) / (nx-1);
1257                if (xp->kinp->isop == NULL) {
1258                        if (x < LOWERLIMIT) x = LOWERLIMIT;
1259                } else {
1260                        if (x < Llimit) x = Llimit;
1261                }
1262                xp->length = xp->kinp->length = x;
1263                alllen += x;
1264                Brnlength[xp->kinp->num] = x;
1265#if             DEBUGFM
1266                printf("b-len %3d%9.4f%9.4f%9.4f%9.4f\n",
1267                        xp->kinp->unm+1, x, xv[i], sy-xv[i], (sy-xv[i])/(nx-2) );
1268#endif  /* DEBUGFM */
1269                if (xp->kinp->isop != NULL) xp->num = xp->kinp->num;
1270        } while (xp != rp);
1271
1272#if     DEBUGFM
1273        for (j = 0; j < Numbrnch; j++) printf("%3d %9.4f\n", j+1, Brnlength[j]);
1274#endif  /* DEBUGFM */
1275
1276        for (i = 0; i < ns-1; i++) {
1277                for (j = i+1; j < ns; j++) dism[i][j] = dmat[i][j];
1278        }
1279        for (i = 0; i < ns; i++) {
1280                x = tr->ebrnchp[i]->length;
1281                path = tr->ebrnchp[i]->paths;
1282                for (j = 0; j < ns; j++) {
1283                        if (path[j]) {
1284                                for (k = 0; k < ns; k++) {
1285                                        if (path[k] == 0) {
1286                                                if (j < k) {
1287                                                        dism[j][k] -= x;
1288#if 0
1289                                                        printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1290                                                                i+1,j+1,k+1,x,dism[j][k],dmat[j][k]);
1291#endif
1292                                                } else {
1293                                                        dism[k][j] -= x;
1294#if 0
1295                                                        printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1296                                                                i+1,j+1,k+1,x,dism[k][j],dmat[k][j]);
1297#endif
1298                                                }
1299                                        }
1300                                }
1301                        }
1302                }
1303        }
1304        for (i = 0; i < Numibrnch; i++) {
1305                x = tr->ibrnchp[i]->length;
1306                if (x == Llimit) x = 0.0; /* caution */
1307                path = tr->ibrnchp[i]->paths;
1308                for (j = 0; j < ns; j++) {
1309                        if (path[j]) {
1310                                for (k = 0; k < ns; k++) {
1311                                        if (path[k] == 0) {
1312                                                if (j < k) {
1313                                                        dism[j][k] -= x;
1314#if 0
1315                                                        printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1316                                                                i+ns+1,j+1,k+1,x,dism[j][k],dmat[j][k]);
1317#endif
1318                                                } else {
1319                                                        dism[k][j] -= x;
1320#if 0
1321                                                        printf("%3d%3d%3d%9.4f%9.4f%9.4f\n",
1322                                                                i+ns+1,j+1,k+1,x,dism[k][j],dmat[k][j]);
1323#endif
1324                                                }
1325                                        }
1326                                }
1327                        }
1328                }
1329        }
1330        for (i = 0, rss = 0.0; i < ns-1; i++) {
1331                for (j = i+1; j < ns; j++) {
1332                        rss += dism[i][j] * dism[i][j];
1333#if 0
1334                        printf("%3d%3d%9.4f%9.4f%9.4f\n",
1335                                i+1, j+1, dism[i][j], dmat[i][j], rss);
1336#endif
1337                }
1338        }
1339
1340        tr->rssleast = rss;
1341        tr->tbldis = alllen;
1342        free_dmatrix(dism);
1343        free_ivector(otux);
1344        free_dvector(xv);
1345        free_dvector(yv);
1346        free_npvector(otup);
1347#if     LSLEN
1348        free_ivector(otun);
1349#endif  /* LSLEN */
1350} /*_ fmlength */
1351
1352
1353extern double probnormal();
1354
1355void
1356resulttree(tr)
1357Tree *tr;
1358{
1359        boolean llimit, rtrif;
1360        int n, ne;
1361        double reli;
1362        Node *ep, *ip;
1363
1364        if (Relitrif != NULL) rtrif = TRUE; else rtrif = FALSE;
1365        printf("\nNo.%-3d", Cnotree + 1); /* offset */
1366        printf("%9s%7s%5s", "ext.","branch","S.E.");
1367        printf("%6s%7s%5s", "int.","branch","S.E.");
1368        if (Relia_optn) printf("  %5s%7s%6s", "LBP ", "2nd ", "pair");
1369        /* if (Info_optn) printf("%14s", "Likelihood"); */
1370        putchar('\n');
1371#if 0
1372        for (n = 0; n < 28; n++) putchar('-');
1373        for (n = 0; n <  5; n++) putchar(' ');
1374        for (n = 0; n < 16; n++) putchar('-'); putchar('\n');
1375#endif
1376        llimit = FALSE;
1377        for (n = 0; n < Numspc; n++) {
1378                ep = tr->ebrnchp[n];
1379                ne = ep->num;
1380                fputid(stdout, Identif[ne], 10);
1381                fputs(" ", stdout);
1382                printf("%3d", ne + 1); /* offset */
1383                if (ep->length == LOWERLIMIT)
1384                        /* printf("%7s %-5s", "lower", "limit"); */
1385                        printf("%7.2f%6.4s", ep->length, "----");
1386                else if (ep->length == Ulimit)
1387                        printf("%7s %-5s", "UPPER", "LIMIT");
1388                else
1389                        printf("%7.2f%6.2f", ep->length, sqrt(ep->kinp->lklhdl));
1390                if (n < Numibrnch) {
1391                        ip = tr->ibrnchp[n];
1392                        printf("%5d", n + 1 + Numspc); /* offset */
1393                        if (ip->length == Llimit) {
1394                                printf("%7s %-5s", "lower", "limit");
1395                                llimit = TRUE;
1396                        } else if (ip->length == Ulimit)
1397                                printf("%7s %-5s", "UPPER", "LIMIT");
1398                        else
1399                                printf("%7.2f%6.2f", ip->length, sqrt(ip->kinp->lklhdl));
1400                        if (Relia_optn && Relistat[n] >= 0) {
1401                                if (Reliprob[n][0] != 1.0)
1402                                        printf("  %.3f", Reliprob[n][0]);
1403                                else
1404                                        printf("  %.1f  ", Reliprob[n][0]);
1405                                (Relistat[n] > 0) ? putchar('*') : putchar(' ');
1406                                if (Reliprob[n][0] != 1.0)
1407                                        printf(" %.3f", Reliprob[n][1]);
1408                                else
1409                                        printf(" %.1f  ", Reliprob[n][1]);
1410                                printf(" %3d&%-3d", Relinum[n][0]+1, Relinum[n][1]+1);
1411                        }
1412#if 0
1413                        reli = probnormal( ip->length / sqrt(ip->kinp->lklhdl) );
1414                        if (ip->length == Llimit)
1415                                printf("%6s", "----");
1416                        else
1417                                printf("%6.2f", reli);
1418#endif
1419                /*      if (Info_optn) printf("%14.3f", ip->lklhdl); */
1420                        putchar('\n');
1421                } else {
1422                        if (n == Numspc - 3) {
1423                                printf("%7s%11.2f", "TBL :",tr->tblength * 1.0);
1424                                printf("%7s %d", "iter:", Numit);
1425                                if (!Converg)
1426                                        printf("%s", "        non convergence!");
1427                                else if (Converg == 2)
1428                                        printf("%s", " just before convergence");
1429                                putchar('\n');
1430                        } else if (n == Numspc - 2) {
1431                                printf("%7s%11.2f +- %.2f\n", "ln L:",
1432                                        tr->lklhd, sqrt(tr->varilkl));
1433                        } else if (n == Numspc - 1) {
1434                                printf("%7s%11.2f", "AIC :",tr->aic);
1435                                if (llimit) printf("%14s %.3f", "lower limit:",Llimit);
1436                                putchar('\n');
1437                        } else {
1438                                putchar('\n');
1439                        }
1440                }
1441        }
1442} /*_ resulttree */
1443
1444
1445void
1446bootstrap(infotrs, lklptrn)
1447Infotree *infotrs;
1448LPMATRIX lklptrn;
1449{
1450        int i, j, k, n, imax, maxtree, nsite, nptrn, same, allsame;
1451        double lklmax, coefrand;
1452        ivector addweight;
1453        dvector boots;
1454
1455        addweight = new_ivector(Numsite);
1456        for ( j = 0, k = 0; k < Numptrn; k++ ) {
1457                for ( i = 0, imax = Weight[k]; i < imax; i++ )
1458                        addweight[j++] = k;
1459        }
1460        if (j != Numsite) {
1461                fputs("ERROR in bootstrap(), sum Weight != Numsite\n", stderr);
1462                exit(1);
1463        }
1464
1465        coefrand = (double)Numsite / ((double)RANDOM_MAX + 1.0);
1466        boots = new_dvector(Numtree);
1467        allsame = 0;
1468        for (n = 0; n < Numtree; n++) infotrs[n].bsprob = 0;
1469
1470#if MLSITE
1471        for (n = 0; n < Numtree; n++) infotrs[n].mlsite = 0;
1472        for (k = 0; k < Numptrn; k++) {
1473                for (n = 1, maxtree = 0, lklmax = lklptrn[0][k]; n < Numtree; n++) {
1474                        if (lklptrn[n][k] > lklmax) {
1475                                lklmax = lklptrn[n][k];
1476                                maxtree = n;
1477                        }
1478                }
1479                infotrs[maxtree].mlsite += Weight[k];
1480        }
1481#endif
1482
1483        for (i = 0; i < NUMBOOTS; i++) {
1484                if (Verbs_optn && i % 100 == 99) fprintf(stderr, " %d", i+1);
1485                for (n = 0; n < Numtree; n++)
1486                        boots[n] = 0.0;
1487                for (k = 0; k < Numsite; k++) {
1488                        nsite = (int)( coefrand * (double)rand() ); /* RANDOM */
1489                        nptrn = addweight[nsite];
1490                        for (n = 0; n < Numtree; n++)
1491                                boots[n] += lklptrn[n][nptrn];
1492                }
1493                maxtree = 0;
1494                lklmax = boots[0];
1495                same = 0;
1496                for (n = 1; n < Numtree; n++) {
1497                        if (boots[n] >= lklmax) {
1498                                if (boots[n] > lklmax) {
1499                                        maxtree = n;
1500                                        lklmax = boots[n];
1501                                        same = 0;
1502                                } else {
1503                                        same++;
1504                                }
1505                        }
1506                }
1507                allsame += same;
1508                infotrs[maxtree].bsprob++;
1509                if (Debug_optn && i < 10)
1510                printf("%3d%4d%10.2f%6d%6d\n", i, maxtree, lklmax, nptrn, nsite);
1511        }
1512        if (allsame > 0)
1513                printf("\nsame bootstrap likelihood occurred %d times\n", allsame);
1514        free_dvector(boots);
1515        free_ivector(addweight);
1516} /*_ bootstrap */
1517
1518
1519extern double probnormal();
1520extern double uprobnormal();
1521
1522void
1523tabletree(infotrs, lklptrn)
1524Infotree *infotrs;
1525LPMATRIX lklptrn;
1526{
1527        int i, k, maxi, n;
1528        double ldiff, suml1, suml2, sdlkl, nn1, z;
1529        Infotree *info;
1530        LPVECTOR mlklptrn, alklptrn;
1531
1532        maxi = 52;
1533        printf("\n%s %s %s", Prog_name, VERSION, Modelname);
1534        printf(" %d trees %d OTUs %d sites. %s\n\n",
1535                Numtree, Numspc, Numsite, Comment);
1536        printf("%4s%9s%11s%6s%6s%6s%10s",
1537        "Tree","ln L","Diff ln L","S.E.","#Para","AIC","Diff AIC");
1538        /* if (Mevol_optn) */
1539        printf("%7s", "   TBL ");
1540        if (Boots_optn) printf("%8s", "RELL-BP");
1541        putchar('\n');
1542        for (i = 0; i < maxi; i++) putchar('-');
1543        /* if (Mevol_optn) */
1544        { for (i = 0; i < 7; i++) putchar('-'); }
1545        if (Boots_optn) { for (i = 0; i < 8; i++) putchar('-'); }
1546        putchar('\n');
1547        mlklptrn = lklptrn[Maxlkltree];
1548        nn1 = (double)Numsite / (double)(Numsite-1);
1549        for (n = 0; n < Numtree; n++) {
1550                info = &infotrs[n];
1551                alklptrn = lklptrn[n];
1552                printf("%-4d%11.1f%8.1f", n+1, info->lklhd, info->lklhd - Maxlkl); /* offset */
1553                if(n == Maxlkltree) {
1554                        fputs(" <-best",stdout);
1555                } else {
1556                        for (suml1 = suml2 = 0.0, k = 0; k < Numptrn; k++) {
1557                                ldiff = alklptrn[k] - mlklptrn[k];
1558                                suml1 += ldiff * Weight[k];
1559                                suml2 += ldiff * ldiff * Weight[k];
1560                        }
1561                        suml1 /= Numsite;
1562                        sdlkl = sqrt( nn1 * (suml2 - suml1*suml1*Numsite) );
1563                        printf("%7.1f", sdlkl);
1564#if 0
1565                        for (suml2 = 0.0, k = 0; k < Numptrn; k++) {
1566                                ldiff = alklptrn[k] - mlklptrn[k] - suml1;
1567                                suml2 += ldiff * ldiff * Weight[k];
1568                        }
1569                        sdlkl = sqrt( nn1 * suml2 );
1570                        printf("%7.1f", sdlkl);
1571#endif
1572                }
1573                printf("%5d%10.1f%7.1f", info->npara, info->aic, info->aic - Minaic);
1574                /* if (Mevol_optn) { */
1575                        if (n == Mintbltree) {
1576                                printf("%6s ", "ME");
1577                        } else {
1578                                printf("%7.1f", info->tblength - Mintbl);
1579                        }
1580                /* } */
1581                if (Boots_optn) printf("%8.4f", (double)info->bsprob / NUMBOOTS);
1582#if MLSITE
1583                if (Boots_optn) printf("%4d", info->mlsite);
1584#endif
1585#if 0
1586                if(n == Maxlkltree) {
1587                        fputs("  Base",stdout);
1588                } else {
1589                        z = (info->lklhd - Maxlkl) / sdlkl;
1590                        printf(" %.3f", probnormal(z)); /* NORMAL */
1591                }
1592#endif
1593                putchar('\n');
1594        }
1595} /*_ tabletree */
1596
1597
1598void
1599tableinfo(infotrs)
1600Infotree *infotrs;
1601{
1602        int n;
1603        Infotree *info;
1604
1605        putchar('\n');
1606        for (n = 0; n < Numtree; n++) {
1607                info = &infotrs[n];
1608                if (Boots_optn) printf("%.4f\t", (double) info->bsprob / NUMBOOTS);
1609                printf("%.1f\t%d\t", info->lklhd - Maxlkl, n+1); /* offset */
1610                puts(info->ltplgy);
1611        }
1612} /*_ tableinfo */
1613
1614#if 0
1615void
1616rerootq(tr, numspc)
1617Tree *tr;
1618int numspc;
1619{
1620        Node *cp, *rp, *op, *xp, *yp;
1621        int i;
1622
1623        for (i = 0; i < numspc; i++) {
1624                if (tr->bturn[i] == numspc - 1)
1625                         break;
1626        }
1627        rp = tr->ebrnchp[i]->kinp;
1628        tr->rootp = rp;
1629        cp = rp;
1630        do {
1631                cp = cp->isop->kinp;
1632                if (cp->isop == NULL) { /* external node */
1633                        cp->descen = TRUE;
1634                        cp = cp->kinp;
1635                        cp->descen = FALSE;
1636                } else { /* internal node */
1637                        if (cp->descen != -1) {
1638                                cp->descen = TRUE;
1639                                cp->kinp->descen = -1;
1640                                tr->ibrnchp[cp->num - Maxspc] = cp;
1641                        } else {
1642                                cp->descen = FALSE;
1643                        }
1644                        if (!cp->descen) {
1645                                op = cp->kinp;
1646                                xp = op->isop;
1647                                yp = xp->isop;
1648                                if (xp->num > yp->num) {
1649                                        op->isop = yp;
1650                                        yp->isop = xp;
1651                                        xp->isop = op;
1652                                }
1653                                cp->num = op->isop->num;
1654                                xp->num = xp->kinp->num;
1655                                yp->num = yp->kinp->num;
1656                        }
1657                }
1658        } while (cp != rp);
1659        op = cp;
1660        xp = op->isop;
1661        yp = xp->isop;
1662        if (xp->num > yp->num) {
1663                op->isop = yp;
1664                yp->isop = xp;
1665                xp->isop = op;
1666        }
1667        xp->num = xp->kinp->num;
1668        yp->num = yp->kinp->num;
1669        op->num = op->kinp->num;
1670} /* rerootq */
1671#endif
1672
1673
1674void
1675outlklhd(lklptrn)
1676LPMATRIX lklptrn;
1677{
1678        int i, j, k, l, maxk;
1679
1680#if 1
1681        fprintf(Lklfp, "%d %d ", Numtree, Numsite);
1682        header(Lklfp, &Maxspc, &Maxsite, &Comment);
1683        for (i = 0; i < Numtree; i++) {
1684                fprintf(Lklfp, "# %d %.1f %s;\n",
1685                        i+1, Infotrees[i].lklhd, Infotrees[i].ltplgy);
1686                for (j = 0, l = 0; j < Numptrn; j++) {
1687                        for (k = 0, maxk = Weight[j]; k < maxk; k++, l++) {
1688                                fprintf(Lklfp, "%16.8e", lklptrn[i][j]);
1689                                if ((l+1) % 5 == 0) fputc('\n', Lklfp);
1690                        }
1691                }
1692                if (l % 5 != 0) fputc('\n', Lklfp);
1693        }
1694#else
1695        for (j = 0, l = 0; j < Numptrn; j++) {
1696                for (k = 0, maxk = Weight[j]; k < maxk; k++, l++) {
1697                        for (i = 1; i < Numtree; i++) {
1698                                fprintf(Lklfp, "%15.8f", lklptrn[0][j] - lklptrn[i][j]);
1699                        }
1700                        fputc('\n', Lklfp);
1701                }
1702        }
1703#endif
1704} /* outlklhd */
1705
1706
1707void
1708putsortseq(tr)
1709Tree *tr;
1710{
1711        Node *rp, *cp;
1712        int i, j, k;
1713
1714        printf("%d %d %s\n", Maxspc, Maxsite, Comment);
1715        cp = rp = tr->rootp;
1716        do {
1717                cp = cp->isop->kinp;
1718                if (cp->isop == NULL) { /* external node */
1719                        i = cp->num;
1720                        printf("%s", Identif[i]);
1721                        if (Sciname[i]) printf(" %s", Sciname[i]);
1722                        putchar('\n');
1723                        for (k = 0; k < Maxsite; k += 60) {
1724                                for (j = k; j < k + 60 && j < Maxsite; j++)
1725                                        putchar(int2acid(Seqchar[i][j]));
1726                                putchar('\n');
1727                        }
1728                        cp = cp->kinp;
1729                }
1730        } while (cp != rp);
1731} /* putsortseq */
Note: See TracBrowser for help on using the repository browser.