src/treespace.c

   1 /* treespace.c\r
   2    collection of tree perturbation routines\r
   3 */\r
   4 \r
   5 #include "paml.h"\r
   6 \r
   7 int MakeTreeIb (int ns, int Ib[], int rooted)\r
   8 {\r
   9 /* construct tree from Ib[] using the algorithm of adding species\r
  10    Ib[k] marks the branch to which the (k+3)_th species (or the root) \r
  11    is added.  Ib[k] should be in the range [0,k+3]\r
  12 */\r
  13    int i,j,k, is,it;\r
  14 \r
  15    tree.nbranch=3;\r
  16    for (i=0; i<3; i++)  { tree.branches[i][0]=3;  tree.branches[i][1]=i; }\r
  17    for (k=0; k<ns-3; k++) {\r
  18       is=k+3;       /* add species (is) to branch Ib[k] */\r
  19 \r
  20       for (i=0; i<tree.nbranch; i++)  for (j=0; j<2; j++)\r
  21          if (tree.branches[i][j]>=is) tree.branches[i][j]+=2;\r
  22       it=tree.branches[Ib[k]][1];\r
  23       tree.branches[Ib[k]][1]=is+1;\r
  24       tree.branches[tree.nbranch][0]=is+1;\r
  25       tree.branches[tree.nbranch++][1]=it;\r
  26       tree.branches[tree.nbranch][0]=is+1;\r
  27       tree.branches[tree.nbranch++][1]=is;\r
  28    }\r
  29    tree.root=tree.branches[0][0];\r
  30    BranchToNode ();\r
  31    \r
  32    if (rooted) {\r
  33       it=tree.branches[Ib[k]][0];\r
  34       tree.branches[Ib[k]][0]=tree.branches[tree.nbranch][0]=2*com.ns-2;\r
  35       tree.branches[tree.nbranch][1]=it;\r
  36       for (; it!=tree.root;  it=nodes[it].father) {\r
  37          tree.branches[nodes[it].ibranch][0]=it;\r
  38          tree.branches[nodes[it].ibranch][1]=nodes[it].father;\r
  39       }\r
  40       tree.root=2*com.ns-2;  tree.nbranch++;\r
  41       BranchToNode ();\r
  42    }\r
  43    return (0);\r
  44 }\r
  45 \r
  46 int GetTreeI (int itree, int ns, int rooted)\r
  47 {\r
  48 /* get the i_th tree.  Trees are ordered according to the algorithm of \r
  49    adding species.\r
  50    returns a random tree if itree==-1, in which case ns can be large\r
  51 */\r
  52    int i, M[NS-2], nM=ns-3+rooted, Ib[NS-2];\r
  53 \r
  54    for (i=0; i<nM-1; i++) M[i]=2*i+5;\r
  55    for (i=0,M[nM-1]=1; i<nM-2; i++) M[nM-1-i-2]*=M[nM-1-i-1];\r
  56 \r
  57    if (itree==-1)  for (i=0; i<nM; i++) Ib[i]=(int)((2*i+3)*rndu());\r
  58    else            for (i=0; i<nM; i++) {Ib[i]=itree/M[i]; itree%=M[i]; } \r
  59 /*\r
  60    if (noisy>3) {\r
  61       FOR (i, nM) printf ("%5d ", M[i]);   FPN (F0);\r
  62       FOR (i, nM) printf ("%5d ", Ib[i]);  FPN (F0);\r
  63    }\r
  64 */\r
  65    MakeTreeIb (ns, Ib, rooted);\r
  66    return (0);\r
  67 }\r
  68 \r
  69 \r
  70 /*\r
  71 int NumberofTrees (int ns, int rooted)\r
  72 {\r
  73    int i, ntree=1;\r
  74 \r
  75    if (ns>15) error2 ("ns too large in NumberofTrees().");\r
  76    for (i=4; i<=ns; i++)  ntree *= 2*i-5;\r
  77    if (rooted) ntree *= 2*i-3;\r
  78    return (ntree);\r
  79 }\r
  80 */\r
  81 \r
  82 \r
  83 int CountLHistories (void)\r
  84 {\r
  85 /* This counts the number of labeled histories for a given rooted tree.\r
  86 */\r
  87    int i,k, nLH, nLR[NS-1][2], change, *sons, j;\r
  88    double y=0;\r
  89 \r
  90    for(i=com.ns; i<tree.nnode; i++) \r
  91      if(nodes[i].nson!=2) error2("this works for rooted trees only");\r
  92    for(i=com.ns; i<tree.nnode; i++) \r
  93       nLR[i-com.ns][0] = nLR[i-com.ns][1] = -1;\r
  94    for(k=0; k<com.ns; k++) {\r
  95       for(i=tree.nnode-1,change=0; i>=com.ns; i--) {\r
  96          sons = nodes[i].sons;\r
  97          for(j=0; j<2; j++) {\r
  98             if(nLR[i-com.ns][j] != -1) continue;\r
  99             if(sons[j] < com.ns) {\r
 100                nLR[i-com.ns][j] = 0;\r
 101                change = 1;\r
 102             }\r
 103             else if(nLR[sons[j]-com.ns][0] != -1 && nLR[sons[j]-com.ns][1] != -1) {\r
 104                nLR[i-com.ns][j] = nLR[sons[j]-com.ns][0] + nLR[sons[j]-com.ns][1] + 1;\r
 105                change = 1;\r
 106             }\r
 107          }\r
 108       }\r
 109       if(!change) break;\r
 110    }\r
 111    for(i=0,nLH=1; i<tree.nnode-com.ns; i++) {\r
 112       /*\r
 113       printf("\nnode %2d (%2d %2d): %2d %2d ", i+com.ns, nodes[i+com.ns].sons[0], nodes[i+com.ns].sons[1], nLR[i][0], nLR[i][1]);\r
 114       */\r
 115       if(nLR[i][0]==-1 || nLR[i][1]==-1)\r
 116          error2("nLR = -1");\r
 117       if(nLR[i][0] && nLR[i][1]) {\r
 118          nLH *= (int)Binomial((double)(nLR[i][0]+nLR[i][1]), min2(nLR[i][0], nLR[i][1]), &y);\r
 119          if(y) error2("y!=0 not considered");\r
 120       }\r
 121    }\r
 122    return(nLH);   \r
 123 }\r
 124 \r
 125 \r
 126 \r
 127 int ListTrees (FILE* fout, int ns, int rooted)\r
 128 {\r
 129 /* list trees by adding species, works fine with large ns\r
 130 */\r
 131    int NTrees, NTreeRoot=3;\r
 132    int i, Ib[NS-2], ns1=ns+rooted, nM=ns1-3, finish;\r
 133 \r
 134    if(com.ns<=12) {\r
 135       printf ("%20s%20s%20s\n", "Taxa", "Unrooted trees", "Rooted trees");\r
 136       for (i=4,NTrees=1; i<=com.ns; i++)  \r
 137          printf ("%20d%20d%20d\n", i, (NTrees*=2*i-5), (NTreeRoot*=2*i-3));\r
 138       fprintf (fout, "%10d %10d\n", com.ns, (!rooted?NTrees:NTreeRoot));\r
 139    }\r
 140 \r
 141    if(com.ns<=26) {\r
 142       for (i=0; i<com.ns; i++)\r
 143          sprintf(com.spname[i], "%d", i+1);\r
 144    }\r
 145 \r
 146    for (i=0;i<nM;i++) Ib[i]=0;\r
 147    for (NTrees=0; ; ) {\r
 148       MakeTreeIb(ns, Ib, rooted);\r
 149       OutTreeN(fout, (com.ns<=26), 0);\r
 150 \r
 151       if(rooted) fprintf(fout, " [%7d %6d LHs]\n", NTrees++, CountLHistories());\r
 152       else fprintf(fout, " [%7d]\n", NTrees++);\r
 153 \r
 154       for (i=nM-1,Ib[nM-1]++,finish=0; i>=0; i--) {\r
 155          if (Ib[i]<2*i+3) break;\r
 156          if (i==0) { \r
 157             finish=1; \r
 158             break;\r
 159          }\r
 160          Ib[i]=0; Ib[i-1]++; \r
 161       }\r
 162       if (finish) break;\r
 163    }\r
 164    FPN(fout);\r
 165 \r
 166    return (0);\r
 167 }\r
 168 \r
 169 int GetIofTree (int rooted, int keeptree, double space[])\r
 170 {\r
 171 /* Get the index of tree.\r
 172    tree.branches[] are destroyed for reconstruction, \r
 173    and some branches are reversed which may affect nodes[] also.\r
 174    Both are restored before return if keeptree=1.\r
 175    Works with binary trees only.\r
 176    bA[nbranch*(ns-2)]\r
 177 */\r
 178    int M[NS-2], nM=com.ns-3+rooted;\r
 179    int i,j,k,is,it, b=0,a1,a2,Ib[NS-1], nid=tree.nnode-com.ns;\r
 180    char ns2=(char)(com.ns-2), *bA=(char*)space;  /* bA[b*ns2+j]: ancestors on branch b */\r
 181    struct TREEB tree0=tree;\r
 182 \r
 183    if (tree.nnode-com.ns!=com.ns-1-!rooted) error2 ("GetIofTree");\r
 184    if (com.ns>15) error2("ns too large in GetIofTree");\r
 185 \r
 186    /* find new root. \r
 187       Ib[]: No. of times inner nodes are visited on paths 1-2, 1-3, 2-3 */\r
 188    for (i=0; i<nid; i++) Ib[i]=0;\r
 189    for (i=0; i<2; i++)  for (j=i+1; j<3; j++)  {\r
 190       for (a1=i; a1!=-1; a1=nodes[a1].father) {\r
 191          for (a2=j; a2!=-1; a2=nodes[a2].father)  if (a1==a2) break;\r
 192          if (a1==a2) break;\r
 193       }\r
 194       for (k=0,Ib[a1-com.ns]++; k<2; k++) {  /* a1 is ancestor of i and j */\r
 195          a2=k?nodes[i].father:nodes[j].father;\r
 196          for (; a2!=a1; a2=nodes[a2].father)  Ib[a2-com.ns]++; \r
 197       }\r
 198    }\r
 199    /* reset root of tree at it */\r
 200    for (it=com.ns; it<com.ns+nid; it++) if (Ib[it-com.ns]==3) break;\r
 201    ReRootTree (it);\r
 202    for (i=0,tree.nbranch=3; i<3; i++)  {\r
 203       tree.branches[i][0]=tree.root;  tree.branches[i][1]=i;\r
 204       for (it=nodes[i].father,k=0; it!=tree.root; it=nodes[it].father)\r
 205          bA[i*ns2+k++]=(char)it;\r
 206       bA[i*ns2+k]=0;\r
 207    }\r
 208    \r
 209    for (is=3; is<com.ns; is++) { /* add species (is) on branch b at node it */\r
 210       for (it=nodes[is].father; it!=tree.root; it=nodes[it].father) {\r
 211     for (b=0; b<tree.nbranch; b++) if (strchr(bA+b*ns2,it)) break;\r
 212     if (b<tree.nbranch) break;\r
 213       }\r
 214       Ib[is-3]=b;\r
 215       tree.branches[tree.nbranch][0]=it;\r
 216       tree.branches[tree.nbranch++][1]=tree.branches[b][1];\r
 217       tree.branches[tree.nbranch][0]=it;\r
 218       tree.branches[tree.nbranch++][1]=is;\r
 219       tree.branches[b][1]=it;\r
 220 \r
 221       if (is==com.ns-1) break;\r
 222       for (i=0; i<3; i++) {  /* modify bA[] for the 3 affected branches */\r
 223          if (i) b=tree.nbranch-3+i;\r
 224          it=nodes[tree.branches[b][1]].father; \r
 225          for (k=0; it!=tree.branches[b][0]; it=nodes[it].father) \r
 226             bA[b*ns2+k++]=(char)it;\r
 227          bA[b*ns2+k]=0;\r
 228       }\r
 229    }  /* for (is) */\r
 230    if (rooted) {\r
 231       a1=nodes[k=tree0.root].sons[0];  a2=nodes[tree0.root].sons[1];\r
 232       if (nodes[a1].father==k)      k=a1;\r
 233       else if (nodes[a2].father==k) k=a2;\r
 234       else error2 ("rooooot");\r
 235       for (b=0; b<tree.nbranch; b++) if (tree.branches[b][1]==k) break;\r
 236       Ib[nM-1]=b;\r
 237    }\r
 238    if (keeptree)  { tree=tree0;  BranchToNode (); }\r
 239 \r
 240    for (i=0; i<nM-1; i++) M[i]=2*i+5;\r
 241    for (i=0,M[nM-1]=1; i<nM-2; i++) M[nM-1-i-2]*=M[nM-1-i-1];\r
 242    for (i=0,it=0; i<nM; i++) it+=Ib[i]*M[i];\r
 243    return (it);\r
 244 }\r
 245 \r
 246 \r
 247 void ReRootTree (int newroot)\r
 248 {\r
 249 /* reroot tree at newroot.  oldroot forgotten\r
 250    The order of branches is not changed.  \r
 251    Branch lengths, and other parameters for branches are updated.\r
 252    Note that node inode needs to be updated if com.oldconP[inode] == 0.\r
 253 */\r
 254    int oldroot=tree.root, a,b;  /* a->b becomes b->a */\r
 255 \r
 256    if (newroot==oldroot) return;\r
 257    for (b=newroot,a=nodes[b].father; b!=oldroot; b=a,a=nodes[b].father) {\r
 258       tree.branches[nodes[b].ibranch][0]=b;\r
 259       tree.branches[nodes[b].ibranch][1]=a;\r
 260 #if (BASEML || CODEML)\r
 261       if(a>=com.ns /* && com.method==1 */) com.oldconP[a]=0;  /* update the node */\r
 262 #endif\r
 263    }\r
 264 \r
 265    tree.root=newroot;\r
 266    BranchToNode ();\r
 267    for (b=oldroot,a=nodes[b].father; b!=newroot; b=a,a=nodes[b].father) { \r
 268       nodes[b].branch = nodes[a].branch; \r
 269       nodes[b].label  = nodes[a].label;\r
 270    }\r
 271    nodes[newroot].branch = -1;\r
 272    nodes[newroot].label = -1;\r
 273 \r
 274 #if (CODEML)\r
 275    /* omega's are moved in updateconP for NSbranchsites models */\r
 276    if(com.model && com.NSsites==0) { \r
 277       for (b=oldroot,a=nodes[b].father; b!=newroot; b=a,a=nodes[b].father)\r
 278          nodes[b].omega = nodes[a].omega;\r
 279       nodes[newroot].omega = -1;\r
 280    }\r
 281 #endif\r
 282 #if (BASEML)\r
 283    if(com.nhomo==2) { \r
 284       for (b=oldroot,a=nodes[b].father; b!=newroot; b=a,a=nodes[b].father)\r
 285          nodes[b].pkappa = nodes[a].pkappa;\r
 286       nodes[newroot].pkappa = NULL;\r
 287    }\r
 288 #endif\r
 289 \r
 290 }\r
 291 \r
 292 \r
 293 \r
 294 int NeighborNNI (int i_tree)\r
 295 {\r
 296 /* get the i_tree'th neighboring tree of tree by the nearest neighbor \r
 297    interchange (NNI), each tree has 2*(# internal branches) neighbors.\r
 298    works with rooted and unrooted binary trees.\r
 299 \r
 300    Gives the ip_th neighbor for interior branch ib. \r
 301    Involved branches are a..b, a..c, b..d,\r
 302    with a..b to be the internal branch.\r
 303    swap c with d, with d to be the ip_th son of b\r
 304 */ \r
 305    int i, a,b,c,d, ib=i_tree/2, ip=i_tree%2;\r
 306 \r
 307    if (tree.nbranch!=com.ns*2-2-(nodes[tree.root].nson==3)) \r
 308       error2 ("err NeighborNNI: multificating tree.");\r
 309 \r
 310    /* locate a,b,c,d */\r
 311    for (i=0,a=0; i<tree.nbranch; i++)\r
 312       if (tree.branches[i][1]>=com.ns && a++==ib) break;\r
 313    ib=i;\r
 314    a=tree.branches[ib][0];  b=tree.branches[ib][1];\r
 315    c=nodes[a].sons[0];      if(c==b) c=nodes[a].sons[1];\r
 316    d=nodes[b].sons[ip];\r
 317 \r
 318    /* swap nodes c and d */\r
 319    tree.branches[nodes[c].ibranch][1]=d;\r
 320    tree.branches[nodes[d].ibranch][1]=c;\r
 321    BranchToNode ();\r
 322    return (0);\r
 323 }\r
 324 \r
 325 int GetLHistoryI (int iLH)\r
 326 {\r
 327 /* Get the ILH_th labelled history.  This function is rather similar to \r
 328    GetTreeI which returns the I_th rooted or unrooted tree topology.\r
 329    The labeled history is recorded in the numbering of nodes: \r
 330    node # increases as the node gets older: \r
 331    node d corresponds to time 2*ns-2-d; tree.root=ns*2-2;\r
 332    t0=1 > t1 > t2 > ... > t[ns-2]\r
 333    k ranges from 0 to i(i-1)/2 and indexes the pair (s1 & s2, with s1<s2),\r
 334    out of i lineages, that coalesce.\r
 335 */\r
 336    int i,k, inode, nodea[NS], s1, s2, it;\r
 337 \r
 338    tree.nnode=com.ns*2-1;\r
 339    for (i=0; i<tree.nnode; i++)  {\r
 340       nodes[i].father=nodes[i].nson=-1;  FOR (k,com.ns) nodes[i].sons[k]=-1;\r
 341    }\r
 342    for (i=0,inode=com.ns; i<com.ns; i++) nodea[i]=i;\r
 343    for (i=com.ns,it=iLH; i>=2; i--)  {\r
 344       k=it%(i*(i-1)/2);  it/=(i*(i-1)/2); \r
 345       s2=(int)(sqrt(1.+8*k)-1)/2+1;  s1=k-s2*(s2-1)/2; /* s1<s2, k=s2*(s2-1)/2+s1 */\r
 346 \r
 347       if (s1>=s2 || s1<0)  printf("\nijk%6d%6d%6d", s1, s2, k);\r
 348 \r
 349       nodes[nodea[s1]].father=nodes[nodea[s2]].father=inode;\r
 350       nodes[inode].nson=2;\r
 351       nodes[inode].sons[0]=nodea[s1];\r
 352       nodes[inode].sons[1]=nodea[s2];\r
 353       nodea[s1]=inode;  nodea[s2]=nodea[i-1]; \r
 354       inode++;\r
 355    }\r
 356    tree.root=inode-1;\r
 357    NodeToBranch();\r
 358    return (0);\r
 359 }\r
 360 \r
 361 int GetIofLHistory (void)\r
 362 {\r
 363 /* Get the index of the labelled history (rooted tree with nodes ordered\r
 364    according to time).  \r
 365    Numbering of nodes: node # increases as the node gets older:\r
 366    node d corresponds to time 2*ns-2-d; tree.root=ns*2-2;\r
 367    t0=1 > t1 > t2 > ... > t[ns-2]\r
 368 */\r
 369    int index, i,j,k[NS+1], inode,nnode, nodea[NS], s[2];\r
 370 \r
 371    if (nodes[tree.root].nson!=2 || tree.nnode!=com.ns*2-1\r
 372      || tree.root!=com.ns*2-2)  error2("IofLH");\r
 373    for (i=0; i<com.ns; i++) nodea[i]=i;\r
 374    for (inode=nnode=com.ns,index=0; inode<com.ns*2-1; inode++,nnode--) {\r
 375       FOR (i,2) FOR (j,nnode)  \r
 376          if (nodes[inode].sons[i]==nodea[j]) { s[i]=j; break; }\r
 377       k[nnode]=max2(s[0],s[1]); s[0]=min2(s[0],s[1]); s[1]=k[nnode];\r
 378       k[nnode]=s[1]*(s[1]-1)/2+s[0];\r
 379       nodea[s[0]]=inode; nodea[s[1]]=nodea[nnode-1];\r
 380    }\r
 381    for (nnode=2,index=0; nnode<=com.ns; nnode++)\r
 382       index=nnode*(nnode-1)/2*index+k[nnode];\r
 383    return (index);\r
 384 }\r
 385 \r
 386 \r
 387 int ReorderNodes (char LHistory[])\r
 388 {\r
 389 /* changes interior node numbers so that the topology represent a labeled\r
 390    history.  LHistory recordes the order of interior nodes with [0] to be \r
 391    root, and [ns-2] the youngest node.   \r
 392    Changes node LHistory[k] to com.ns*2-3-k.  \r
 393    Uses only tree but not nodes[] but transforms both tree and nodes into \r
 394    the labeled history.\r
 395 */\r
 396    int i, j, k;\r
 397 \r
 398    if (tree.root!=com.ns*2-2 || LHistory[0]!=com.ns*2-2) {\r
 399       tree.root=com.ns*2-2;\r
 400 /*      printf("\nRoot changed to %d in ReorderNodes..\n", com.ns*2-2+1); */\r
 401    }\r
 402    FOR (i, tree.nbranch) FOR (j,2) \r
 403       if (tree.branches[i][j]>=com.ns) \r
 404          FOR (k, com.ns-1)\r
 405             if (tree.branches[i][j]==LHistory[k]) \r
 406                { tree.branches[i][j]=com.ns*2-2-k;  break; }\r
 407    BranchToNode();\r
 408 \r
 409    return (0);\r
 410 }\r