| 1 | /* |
|---|
| 2 | * qltree.c Adachi, J. 1994.01.11 |
|---|
| 3 | * Copyright (C) 1992, 1993 J. Adachi & M. Hasegawa, All rights reserved. |
|---|
| 4 | */ |
|---|
| 5 | |
|---|
| 6 | #include "protml.h" |
|---|
| 7 | |
|---|
| 8 | |
|---|
| 9 | Infoqltree * |
|---|
| 10 | newinfoqltrees(n, maxbrnch) |
|---|
| 11 | int n; |
|---|
| 12 | int maxbrnch; |
|---|
| 13 | { |
|---|
| 14 | int i; |
|---|
| 15 | Infoqltree *info; |
|---|
| 16 | double *v; |
|---|
| 17 | |
|---|
| 18 | info = (Infoqltree *)malloc((unsigned)n * sizeof(Infoqltree)); |
|---|
| 19 | if (info == NULL) maerror("1 in newinfoqltrees()."); |
|---|
| 20 | v = (double *)malloc((unsigned)(n * maxbrnch) * sizeof(double)); |
|---|
| 21 | if (v == NULL) maerror("2 in newinfoqltrees()."); |
|---|
| 22 | for (i = 0; i < n; i++, v+=maxbrnch) |
|---|
| 23 | info[i].lengths = v; |
|---|
| 24 | |
|---|
| 25 | return info; |
|---|
| 26 | } /*_ newinfoqltrees */ |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | Infoaddtree * |
|---|
| 30 | newinfoaddtree(buftree) |
|---|
| 31 | int buftree; |
|---|
| 32 | { |
|---|
| 33 | Infoaddtree *info; |
|---|
| 34 | cvector m; |
|---|
| 35 | |
|---|
| 36 | info = (Infoaddtree *) malloc(sizeof(Infoaddtree)); |
|---|
| 37 | if (info == NULL) maerror("1 in newinfoaddtree()."); |
|---|
| 38 | m = (cvector) malloc((unsigned)buftree * sizeof(char)); |
|---|
| 39 | if (m == NULL) maerror("2 in newinfoaddtree()."); |
|---|
| 40 | info->ltplgy = m; |
|---|
| 41 | info->lklaprox = 0.0; |
|---|
| 42 | info->frequency = 0; |
|---|
| 43 | info->dp = NULL; |
|---|
| 44 | |
|---|
| 45 | return info; |
|---|
| 46 | } /*_ newinfoaddtree */ |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | void initturn(tr) |
|---|
| 50 | Tree *tr; |
|---|
| 51 | { |
|---|
| 52 | int i, j, mini, n; |
|---|
| 53 | double dis, mindis; |
|---|
| 54 | ivector turn; |
|---|
| 55 | Node *dp, *ap; |
|---|
| 56 | |
|---|
| 57 | mindis = fabs(Distanmat[0][1] - Distanmat[Numspc-1][1]); |
|---|
| 58 | mini = 1; |
|---|
| 59 | for ( i = 2; i < Numspc - 1; i++) { |
|---|
| 60 | dis = fabs(Distanmat[0][i] - Distanmat[Numspc-1][i]); |
|---|
| 61 | if (dis < mindis) { |
|---|
| 62 | mindis = dis; |
|---|
| 63 | mini = i; |
|---|
| 64 | } |
|---|
| 65 | } |
|---|
| 66 | turn = tr->bturn; |
|---|
| 67 | turn[0] = 0; |
|---|
| 68 | turn[1] = mini; |
|---|
| 69 | turn[2] = Numspc - 1; |
|---|
| 70 | for (i = 1, j = 3; i < Numspc - 1; i++) { |
|---|
| 71 | if (i != mini) |
|---|
| 72 | turn[j++] = i; |
|---|
| 73 | } |
|---|
| 74 | |
|---|
| 75 | for (i = 0; i < Maxspc; i++) { |
|---|
| 76 | n = turn[i]; |
|---|
| 77 | dp = tr->ebrnchp[i]; |
|---|
| 78 | ap = dp->kinp; |
|---|
| 79 | dp->num = n; |
|---|
| 80 | ap->num = n; |
|---|
| 81 | dp->eprob = Seqconint[n]; |
|---|
| 82 | } |
|---|
| 83 | } /* initturn */ |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | void randturn(tr) |
|---|
| 87 | Tree *tr; |
|---|
| 88 | { |
|---|
| 89 | int i, j, temp, n; |
|---|
| 90 | ivector turn; |
|---|
| 91 | Node *dp, *ap; |
|---|
| 92 | |
|---|
| 93 | turn = tr->bturn; |
|---|
| 94 | for (i = Maxspc - 1; i > 0; i--) { |
|---|
| 95 | j = (int)(((i + 1) / (RANDOM_MAX + 1.0)) * rand()); |
|---|
| 96 | temp = turn[i]; |
|---|
| 97 | turn[i] = turn[j]; |
|---|
| 98 | turn[j] = temp; |
|---|
| 99 | } |
|---|
| 100 | if (turn[0] > turn[1]) { temp=turn[0]; turn[0]=turn[1]; turn[1]=temp; } |
|---|
| 101 | if (turn[0] > turn[2]) { temp=turn[0]; turn[0]=turn[2]; turn[2]=temp; } |
|---|
| 102 | if (turn[1] > turn[2]) { temp=turn[1]; turn[1]=turn[2]; turn[2]=temp; } |
|---|
| 103 | |
|---|
| 104 | for (i = 0; i < Maxspc; i++) { |
|---|
| 105 | n = turn[i]; |
|---|
| 106 | dp = tr->ebrnchp[i]; |
|---|
| 107 | ap = dp->kinp; |
|---|
| 108 | dp->num = n; |
|---|
| 109 | ap->num = n; |
|---|
| 110 | dp->eprob = Seqconint[n]; |
|---|
| 111 | } |
|---|
| 112 | } /* randturn */ |
|---|
| 113 | |
|---|
| 114 | |
|---|
| 115 | void |
|---|
| 116 | convertdistan(tr, numspc, distanmat, distanvec) |
|---|
| 117 | dmatrix distanmat; |
|---|
| 118 | dvector distanvec; |
|---|
| 119 | Tree *tr; |
|---|
| 120 | int numspc; |
|---|
| 121 | { |
|---|
| 122 | int i, j, k; |
|---|
| 123 | ivector turn; |
|---|
| 124 | |
|---|
| 125 | turn = tr->bturn; |
|---|
| 126 | for (i = 1, k = 0; i < numspc; i++) { |
|---|
| 127 | for (j = 0; j < i; j++, k++) |
|---|
| 128 | distanvec[k] = |
|---|
| 129 | distanmat[turn[i]][turn[j]]; |
|---|
| 130 | } |
|---|
| 131 | } |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | void |
|---|
| 135 | praproxlkl2(fp, tr) |
|---|
| 136 | FILE *fp; |
|---|
| 137 | Tree *tr; |
|---|
| 138 | { |
|---|
| 139 | fprintf(fp, "%.1f\t%.3f\t",tr->lklhd,tr->rssleast); |
|---|
| 140 | } /*_ praproxlkl2 */ |
|---|
| 141 | |
|---|
| 142 | |
|---|
| 143 | int addotu(tr, cp, np, ip, cnspc) |
|---|
| 144 | Tree *tr; |
|---|
| 145 | Node *cp, *np, *ip; |
|---|
| 146 | int cnspc; |
|---|
| 147 | { |
|---|
| 148 | Node *op; |
|---|
| 149 | |
|---|
| 150 | if (cp == tr->rootp) tr->rootp = ip; |
|---|
| 151 | for (op = cp->isop->isop; op->isop != cp; op = op->isop) |
|---|
| 152 | ; |
|---|
| 153 | ip->isop = cp->isop; |
|---|
| 154 | op->isop = ip; |
|---|
| 155 | ip->kinp->isop = cp; |
|---|
| 156 | cp->isop = np; |
|---|
| 157 | np->isop = ip->kinp; |
|---|
| 158 | |
|---|
| 159 | if (Debug) prcurtree(tr); |
|---|
| 160 | pathing(tr); |
|---|
| 161 | lslength(tr, Distanvec, cnspc+1); |
|---|
| 162 | tr->lklhd = 0.0; |
|---|
| 163 | /* initpartlkl(tr); */ |
|---|
| 164 | /* aproxlkl(tr); */ |
|---|
| 165 | /* if (Logfl_optn) { |
|---|
| 166 | praproxlkl2(Logfp, tr); |
|---|
| 167 | fputctopology(Logfp, tr); |
|---|
| 168 | } */ |
|---|
| 169 | |
|---|
| 170 | cp->isop = ip->isop; |
|---|
| 171 | op->isop = cp; |
|---|
| 172 | np->isop = NULL; |
|---|
| 173 | ip->isop = NULL; |
|---|
| 174 | ip->kinp->isop = NULL; |
|---|
| 175 | if (ip == tr->rootp) tr->rootp = cp; |
|---|
| 176 | } /* addotu */ |
|---|
| 177 | |
|---|
| 178 | |
|---|
| 179 | addotual(tr, cp, np, ip, lengs) |
|---|
| 180 | Tree *tr; |
|---|
| 181 | Node *cp, *np, *ip; |
|---|
| 182 | dvector lengs; |
|---|
| 183 | { |
|---|
| 184 | Node *op, **bp; |
|---|
| 185 | int i, j; |
|---|
| 186 | |
|---|
| 187 | if (cp == tr->rootp) tr->rootp = ip; |
|---|
| 188 | for (op = cp->isop->isop; op->isop != cp; op = op->isop) |
|---|
| 189 | ; |
|---|
| 190 | ip->isop = cp->isop; |
|---|
| 191 | op->isop = ip; |
|---|
| 192 | ip->kinp->isop = cp; |
|---|
| 193 | cp->isop = np; |
|---|
| 194 | np->isop = ip->kinp; |
|---|
| 195 | |
|---|
| 196 | bp = tr->ebrnchp; |
|---|
| 197 | for (i = 0; i < Numspc; i++) { |
|---|
| 198 | bp[i]->length = lengs[i]; |
|---|
| 199 | bp[i]->kinp->length = lengs[i]; |
|---|
| 200 | } |
|---|
| 201 | bp = tr->ibrnchp; |
|---|
| 202 | for (j = 0, i = Numspc; j < Numibrnch; j++, i++) { |
|---|
| 203 | bp[j]->length = lengs[i]; |
|---|
| 204 | bp[j]->kinp->length = lengs[i]; |
|---|
| 205 | } |
|---|
| 206 | /* for (i = 0; i < Numbrnch; i++) printf("%6.2f",lengs[i]); putchar('\n'); */ |
|---|
| 207 | if (Debug) prcurtree(tr); |
|---|
| 208 | initpartlkl(tr); |
|---|
| 209 | aproxlkl(tr); |
|---|
| 210 | if (Logfl_optn) { |
|---|
| 211 | praproxlkl2(Logfp, tr); |
|---|
| 212 | fputctopology(Logfp, tr); |
|---|
| 213 | } |
|---|
| 214 | |
|---|
| 215 | cp->isop = ip->isop; |
|---|
| 216 | op->isop = cp; |
|---|
| 217 | np->isop = NULL; |
|---|
| 218 | ip->isop = NULL; |
|---|
| 219 | ip->kinp->isop = NULL; |
|---|
| 220 | if (ip == tr->rootp) tr->rootp = cp; |
|---|
| 221 | } /* addotual */ |
|---|
| 222 | |
|---|
| 223 | |
|---|
| 224 | void |
|---|
| 225 | roundtree(tr, cnspc, infoqltrees, qhead, qtail) |
|---|
| 226 | Tree *tr; |
|---|
| 227 | int cnspc; |
|---|
| 228 | Infoqltree *infoqltrees, *qhead, *qtail; |
|---|
| 229 | { |
|---|
| 230 | Node *cp, *rp, *np, *ip, *op, *ap, *maxlkp; |
|---|
| 231 | Infoqltree *xp, *yp, *zp; |
|---|
| 232 | boolean added; |
|---|
| 233 | int i, ntree; |
|---|
| 234 | double res, minres, rate; |
|---|
| 235 | dvector lengs; |
|---|
| 236 | |
|---|
| 237 | minres = 1e+36; |
|---|
| 238 | qtail = qhead; |
|---|
| 239 | np = tr->ebrnchp[cnspc]->kinp; |
|---|
| 240 | ip = tr->ibrnchp[cnspc - 3]->kinp; |
|---|
| 241 | ntree = 0; |
|---|
| 242 | cp = rp = tr->rootp; |
|---|
| 243 | do { |
|---|
| 244 | added = FALSE; |
|---|
| 245 | cp = cp->isop->kinp; |
|---|
| 246 | if (cp->isop == NULL) { /* external node */ |
|---|
| 247 | addotu(tr, cp->kinp, np, ip, cnspc); |
|---|
| 248 | added = TRUE; |
|---|
| 249 | ap = cp->kinp; |
|---|
| 250 | cp = cp->kinp; /* not descen */ |
|---|
| 251 | } else { /* internal node */ |
|---|
| 252 | if (cp->descen) { |
|---|
| 253 | addotu(tr, cp->kinp, np, ip, cnspc); |
|---|
| 254 | ap = cp->kinp; |
|---|
| 255 | added = TRUE; |
|---|
| 256 | } |
|---|
| 257 | } |
|---|
| 258 | if (added) { |
|---|
| 259 | res = tr->rssleast; |
|---|
| 260 | if (res < qtail->residual || ntree < MAXQLBUF) { |
|---|
| 261 | if (res < minres) |
|---|
| 262 | minres = res; |
|---|
| 263 | if (ntree < MAXQLBUF) { |
|---|
| 264 | yp = &infoqltrees[ntree]; |
|---|
| 265 | } else if (res < qtail->residual) { |
|---|
| 266 | yp = qtail; |
|---|
| 267 | qtail = yp->up; |
|---|
| 268 | } |
|---|
| 269 | yp->residual = res; |
|---|
| 270 | yp->ap = ap; |
|---|
| 271 | lengs = yp->lengths; |
|---|
| 272 | for (xp = qtail; res < xp->residual; zp = xp, xp = xp->up) |
|---|
| 273 | ; |
|---|
| 274 | yp->up = xp; |
|---|
| 275 | if (qtail == xp) |
|---|
| 276 | qtail = yp; |
|---|
| 277 | else |
|---|
| 278 | zp->up = yp; |
|---|
| 279 | for (i = 0; i < Numbrnch; i++) |
|---|
| 280 | lengs[i] = Brnlength[i]; |
|---|
| 281 | /* for (i = 0; i < Numbrnch; i++) printf("%6.2f",Brnlength[i]); |
|---|
| 282 | putchar('\n'); |
|---|
| 283 | for (i = 0; i < Numbrnch; i++) printf("%6.2f",lengs[i]); |
|---|
| 284 | putchar('\n'); */ |
|---|
| 285 | } |
|---|
| 286 | ntree++; |
|---|
| 287 | } |
|---|
| 288 | } while (cp != rp); |
|---|
| 289 | |
|---|
| 290 | if (Logfl_optn) putc('\n', Logfp); |
|---|
| 291 | if ((cnspc * 2 - 3) > MAXQLBUF) |
|---|
| 292 | ntree = MAXQLBUF; |
|---|
| 293 | else |
|---|
| 294 | ntree = cnspc * 2 - 3; |
|---|
| 295 | tr->lklmean = -1e+36; |
|---|
| 296 | xp = qtail; |
|---|
| 297 | while (xp != qhead) { |
|---|
| 298 | ntree--; |
|---|
| 299 | rate = xp->residual / minres; /* !? */ |
|---|
| 300 | if (rate < QRSRATE || ntree < 3) { |
|---|
| 301 | /* printf("%.3f %d\n", rate, ntree); */ |
|---|
| 302 | tr->rssleast = rate; |
|---|
| 303 | addotual(tr, xp->ap, np, ip, xp->lengths); |
|---|
| 304 | if (tr->lklhd > tr->lklmean) { |
|---|
| 305 | maxlkp = xp->ap; |
|---|
| 306 | tr->lklmean = tr->lklhd; |
|---|
| 307 | } |
|---|
| 308 | } |
|---|
| 309 | xp = xp->up; |
|---|
| 310 | } |
|---|
| 311 | |
|---|
| 312 | if (maxlkp == tr->rootp) tr->rootp = ip; |
|---|
| 313 | for (op = maxlkp->isop->isop; op->isop != maxlkp; op = op->isop) |
|---|
| 314 | ; |
|---|
| 315 | ip->isop = maxlkp->isop; |
|---|
| 316 | op->isop = ip; |
|---|
| 317 | ip->kinp->isop = maxlkp; |
|---|
| 318 | maxlkp->isop = np; |
|---|
| 319 | np->isop = ip->kinp; |
|---|
| 320 | #if 0 |
|---|
| 321 | if (Debug) prcurtree(tr); |
|---|
| 322 | pathing(tr); |
|---|
| 323 | lslength(tr, Distanvec, cnspc+1); |
|---|
| 324 | if (Logfl_optn) fprintf(Logfp, "\t%.3f\n", tr->rssleast / minres); /* !? */ |
|---|
| 325 | initpartlkl(tr); |
|---|
| 326 | aproxlkl(tr); |
|---|
| 327 | if (Logfl_optn) praproxlkl2(Logfp, tr); |
|---|
| 328 | if (Logfl_optn) fputctopology(Logfp, tr); |
|---|
| 329 | #endif |
|---|
| 330 | } /*_ roundtree */ |
|---|
| 331 | |
|---|
| 332 | |
|---|
| 333 | void |
|---|
| 334 | qtreeinit(tr) |
|---|
| 335 | Tree *tr; |
|---|
| 336 | { |
|---|
| 337 | Node *sp0, *sp1, *sp2; |
|---|
| 338 | |
|---|
| 339 | sp0 = tr->ebrnchp[0]->kinp; |
|---|
| 340 | sp1 = tr->ebrnchp[1]->kinp; |
|---|
| 341 | sp2 = tr->ebrnchp[2]->kinp; |
|---|
| 342 | sp0->isop = sp1; |
|---|
| 343 | sp1->isop = sp2; |
|---|
| 344 | sp2->isop = sp0; |
|---|
| 345 | tr->rootp = sp2; |
|---|
| 346 | } /*_ qtreeinit */ |
|---|
| 347 | |
|---|
| 348 | |
|---|
| 349 | void |
|---|
| 350 | tableaddtree(head, numaddtree) |
|---|
| 351 | Infoaddtree *head; |
|---|
| 352 | int numaddtree; |
|---|
| 353 | { |
|---|
| 354 | Infoaddtree *cp, *tail; |
|---|
| 355 | |
|---|
| 356 | printf("%d / %d ", numaddtree, Cnotree); |
|---|
| 357 | fputs(Modelname, stdout); |
|---|
| 358 | fputs(" model", stdout); |
|---|
| 359 | for (cp = head; cp != NULL; cp = cp->dp) { |
|---|
| 360 | tail = cp; |
|---|
| 361 | } |
|---|
| 362 | printf(" approx ln L %.1f", head->lklaprox); |
|---|
| 363 | printf(" ... %.1f", tail->lklaprox); |
|---|
| 364 | printf(" diff %.1f\n", head->lklaprox - tail->lklaprox); |
|---|
| 365 | for (cp = head; cp != NULL; cp = cp->dp) { |
|---|
| 366 | if (Info_optn ) printf("%d\t%.1f ", cp->frequency, cp->lklaprox); |
|---|
| 367 | fputs(cp->ltplgy, stdout); |
|---|
| 368 | fputs(";\n", stdout); |
|---|
| 369 | } |
|---|
| 370 | } /*_ tableaddtree */ |
|---|
