src/triangMtd.c

   1 /* triangMtd.c    2012-04-02 */
   2
   3 /* Copyright 2011-2012 Andrei-Alin Popescu */
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 /*
   9  * leafs labelled 1 to n. root labelled n+1. other nodes labelled n+1 to m
  10  */
  11
  12 #include <R.h>
  13
  14 int give_indexx(int i, int j, int n)
  15 {       if (i == j) return 0;
  16         if (i > j) return(n*(j - 1) - j*(j - 1)/2 + i - j - 1);
  17         else return(n*(i - 1) - i*(i - 1)/2 + j - i - 1);
  18 }
  19
  20 int pred(int k, int* ed1, int* ed2, int numEdges)
  21 /* find the predecesor of vertex k */
  22 {
  23         int i = 0;
  24         for (i = 0; i <= numEdges; i++)
  25                 if (ed2[i] == k) return ed1[i];
  26         return -1;
  27 }
  28
  29 int* getPathBetween(int x, int y, int n, int* ed1, int* ed2, int numEdges)
  30 //get the path between vertices x and y in an array ord.
  31 //ord[i]=j means that we go between i and j on the path between x and y
  32  {  int i=0;
  33     int k=x;
  34     int ch[2*n-1];//ch[i]==1 implies {k,pred(k)} on path between x and y
  35     for(i=1;i<=2*n-2;i++)
  36       {ch[i]=0;
  37       }
  38
  39     while(k!=n+1)
  40       {
  41         ch[k]=1;
  42         k=pred(k,ed1,ed2,numEdges);
  43       }
  44     k=y;
  45
  46     while(k!=n+1)
  47       {
  48         ch[k]++;
  49         k=pred(k,ed1,ed2,numEdges);
  50       }
  51         int *ord=(int*)malloc(sizeof(int)*(2*n-1));
  52         //starting from x, fill ord
  53
  54
  55         int p=x;
  56
  57         while(ch[p]==1)
  58           {
  59             int aux=p;
  60             p=pred(p,ed1,ed2,numEdges);
  61             ord[aux]=p;
  62           }
  63         p=y;
  64
  65         while(ch[p]==1)
  66           {
  67             int aux=p;
  68             p=pred(p,ed1,ed2,numEdges);
  69             ord[p]=aux;//other way
  70           }
  71
  72       return ord;
  73  }
  74
  75 int getLength(int x, int y, int* ed1, int* ed2, int numEdges, int* edLen)
  76 /* get length of edge {x,y}, -1 if edge does not exist */
  77 {
  78         int i = 0;
  79         for (i = 0; i <= numEdges; i++)
  80                 if ((ed1[i] == x && ed2[i] == y) || (ed1[i] == y && ed2[i] == x))
  81                         return edLen[i];
  82         return -1;
  83 }
  84
  85 void triangMtd(double* d, int* np, int* ed1,int* ed2, double* edLen)
  86 {
  87     int n=*np;
  88     int i=0;
  89     int j=0;
  90     int ij=-1;
  91
  92     for(i=0;i<n-1;i++)
  93     {
  94      for(j=i+1;j<n;j++)
  95        {ij++;
  96          //error("%f ,giveindex= %f",d[ij],d[give_indexx(i,j,n)]);
  97        }
  98     }
  99
 100     /*double d[n*n];
 101
 102
 103     for(i=0;i<n;i++)
 104     {d[i*n+i]=0;
 105      for(j=i;j<n;j++)
 106        {d[i*n+j]=d[j*n+i];
 107        }
 108     }
 109     for(i=0;i<n;i++)
 110     {
 111      for(j=0;j<n;j++)
 112        {error("%f ",d[i*n+j]);
 113        }
 114      error("\n");
 115     }*/
 116
 117     double minW=0;
 118     int x=-1,y=-1,numEdges=0;
 119     int st[n+1];//array holding at position i the starting point of attachment path of leaf i
 120     int ed[n+1];//array holding at position i the ending point of attachment path of leaf i
 121     double l[n+1];//array holding at position i the distance of leaf i from the partially constructed tree
 122     int w[n+1];//presence array:w[i]=1 if leaf i is added to the tree, 0 otherwise
 123     int wSize=0;//number of leaves added so far
 124
 125     for(i=1;i<=n;i++){w[i]=0;}
 126     //find the minimum distance in our matrix
 127
 128     int sw=0;
 129     //choose two leaves x,y such that d[x][y] is minimal
 130     for(i=1;i<n;i++)
 131     {
 132      for(j=i+1;j<=n;j++)
 133        {if(d[give_indexx(i,j,n)]<=0)continue;
 134         minW=d[give_indexx(i,j,n)];
 135         x=i;
 136         y=j;
 137         sw=1;
 138         break;
 139        }
 140      if(sw)break;
 141     }
 142     for(i=1;i<n;i++)
 143      for(j=i+1;j<=n;j++)
 144        { if(i==j)continue;
 145          if(d[give_indexx(i,j,n)]<minW)
 146            {
 147              minW=d[give_indexx(i,j,n)];
 148              x=i;
 149              y=j;
 150            }
 151        }
 152
 153     w[x]=1; w[y]=1;//mark x and y as added
 154
 155     //calculate the distance between leaves not in w and leaves in w
 156     for(i=1;i<=n;i++)
 157       {
 158         if(w[i]){continue;}
 159         st[i]=x;ed[i]=y;
 160         l[i]=0.5*(d[give_indexx(i,x,n)]+d[give_indexx(i,y,n)]-d[give_indexx(x,y,n)]);//distance of leaf i from the
 161 //initial tree
 162       }
 163
 164
 165     wSize+=3;//since we construct a star tree on three leaves
 166     int nv=n+1;//since first n numbers are reserved for leaves
 167     double minDist=1000;
 168     int x3=1;
 169     //search for x3 that is closest to the tree
 170     for(i=1;i<=n;i++)
 171           {if(w[i])continue;//skip if leaf already added
 172             if(l[i]<minDist)
 173               {
 174                 minDist=l[i];
 175                 x3=i;
 176               }
 177           }
 178     //construct initial star tree on three leaves:x,y,x3. nv is the interior
 179     //vertex and also the root of the tree
 180     ed1[numEdges]=nv;
 181     ed2[numEdges]=x;
 182     edLen[numEdges]=0.5*(d[give_indexx(x3,x,n)]+d[give_indexx(y,x,n)]-d[give_indexx(y,x3,n)]);
 183     numEdges++;
 184     ed1[numEdges]=nv;
 185     ed2[numEdges]=y;
 186     edLen[numEdges]=0.5*(d[give_indexx(y,x3,n)]+d[give_indexx(y,x,n)]-d[give_indexx(x,x3,n)]);
 187     numEdges++;
 188     ed1[numEdges]=nv;
 189     ed2[numEdges]=x3;
 190     edLen[numEdges]=0.5*(d[give_indexx(x3,x,n)]+d[give_indexx(y,x3,n)]-d[give_indexx(y,x,n)]);
 191     w[x3]=1;
 192
 193         /*Rprintf("new tree\n");
 194         for(i=0;i<2*n-3;i++)
 195         {
 196         Rprintf("%i->%i of length:%f \n",ed1[i],ed2[i],edLen[i]);
 197         }
 198         Rprintf("end new tree\n");*/
 199
 200     //calculate distance of leaves not yet added to the star tree
 201     int s;
 202     for(s=1;s<=n;s++)
 203          {if(w[s])continue;
 204            for(i=1;i<=n;i++)
 205             {  if(i==x3)continue;
 206                if(!w[i])continue;
 207                double newL=0.5*(d[give_indexx(i,s,n)]+d[give_indexx(s,x3,n)]-d[give_indexx(i,x3,n)]);
 208                if(newL<l[s])
 209                   {
 210                    l[s]=newL;
 211                    st[s]=i;
 212                    ed[s]=x3;
 213                   }
 214             }
 215          }
 216
 217     //partial tree construction begins here
 218
 219     while(wSize<n)
 220       { double minDist=1e50;
 221         int z=1;
 222
 223         //search for leaf z which is closest to partial tree
 224         for(i=1;i<=n;i++)
 225           {if(w[i])continue;//skip if leaf already added
 226             if(l[i]<minDist)
 227               {
 228                 minDist=l[i];
 229                 z=i;
 230               }
 231           }
 232         //x and y are the leaves on the path between which z is attached
 233         x=st[z];y=ed[z];
 234
 235
 236         int* ord=getPathBetween(x,y,n,ed1,ed2,numEdges);
 237        /* Rprintf("ord\n");
 238         for(i=1;i<=2*n-2;i++)
 239           {Rprintf("ord[%i]=%i ",i,ord[i]);
 240           }
 241         Rprintf("\n");*/
 242         //order the path from x to y, in an array ord, ord[i]=j means vertex i comes before j
 243         //first count number of edges on path (i.e count all i s.t ch[i]==1)
 244
 245
 246         //look for the edge on the path x to y to subdivide
 247         int p=x;
 248         double sum=0;
 249         double prevSum=0;
 250         int aux=0;
 251         int subdiv=-1;//index of edge to subdivide
 252         //error("d[y,x]=%f,d[z,x]=%f,d[z,y]=%f\n",d[give_indexx(y,x,n)],d[give_indexx(z,x,n)],d[give_indexx(z,y,n)]);
 253         double lx=0.5*(d[give_indexx(y,x,n)]+d[give_indexx(z,x,n)]-d[give_indexx(z,y,n)]);//distance of attachment
 254        // Rprintf("adding %i on the path between %i and %i at a distance from x of %f and a distance of %f from tree",z,x,y,lx,minDist);
 255         //point from x
 256         //Rprintf("Adding leaf %i, between %i and %i\n",z,x,y);
 257        // Rprintf("%i situated at a distance of %d from tree",z,minDist);
 258         int sw=0;
 259         //Rprintf("path between %i and %i\n",x,y);
 260
 261         while(p!=y && sum<lx)
 262           {
 263             aux=p;
 264             //error("%i -> %i ",p,ord[p]);
 265             p=ord[p];
 266             prevSum=sum;
 267             for(i=0;i<=numEdges;i++)
 268               {
 269                 if((ed1[i]==aux && ed2[i]==p)||(ed2[i]==aux && ed1[i]==p))
 270                   {
 271                     if(ed1[i]==aux && ed2[i]==p){sw=1;}
 272                     subdiv=i;
 273                     sum+=edLen[i];
 274                   }
 275               }
 276           }
 277
 278         nv++;
 279         //subdivide subdiv with a node labelled nv
 280         //length calculation
 281         //multifurcating vertices
 282
 283         int edd=ed2[subdiv];
 284         ed2[subdiv]=nv;
 285         edLen[subdiv]= (sw==1)?(lx-prevSum):(sum-lx);//check which 'half' of the
 286                                                      //path the leaf belongs to
 287                                                      //and updates accordingly
 288         //error("sum=%f, prevsum=%f\n",sum,prevSum);
 289         //error("lx-prevSum=%f, sum-lx=%f, minDist=%f",lx-prevSum,sum-lx,minDist);
 290         numEdges++;
 291         ed1[numEdges]=nv;
 292         ed2[numEdges]=edd;
 293         edLen[numEdges]= (sw==1)?(sum-lx):(lx-prevSum);
 294         numEdges++;
 295         edLen[numEdges]=minDist;
 296         ed1[numEdges]=nv;
 297         ed2[numEdges]=z;
 298
 299         wSize++;
 300         w[z]=1;
 301
 302         /*update distance matrix, only needed for incomplete distances
 303         int ii;
 304       for(ii=0;ii<n;ii++)
 305       {if(!w[ii])continue;
 306        error("path between %i and %i\n",ii,z);
 307         int* ord=getPathBetween(ii,z,n,ed1,ed2,numEdges);
 308
 309         p=ii;
 310         int newDist=0;
 311         error("path");
 312         for(i=0;i<2*n-2;i++)
 313           {error("ord[%i]=%i ",i,ord[i]);
 314           }
 315        error("\n");
 316        error("innn\n");
 317         while(p!=z)
 318           { //error("p=%i\n",p);
 319             int aux=p;
 320             p=ord[p];
 321             newDist+=getLength(ii,z,ed1,ed2,numEdges,edLen);
 322           }
 323
 324         error("outt\n");
 325
 326         d[ii][z]=d[z][ii]=newDist;
 327        }*/
 328         //update l[s] for all s not yet added
 329         int s;
 330         for(s=1;s<=n;s++)
 331          {if(w[s])continue;
 332            for(i=1;i<=n;i++)
 333             {if(i==z)continue;
 334                if(i!=x && i!=y)continue;//we only consider x and y as being other leaf
 335                                         //and take the minimum of them as being new distance
 336                double newL=0.5*(d[give_indexx(i,s,n)]+d[give_indexx(z,s,n)]-d[give_indexx(i,z,n)]);//one of leaves is
 337                                                       //z, since
 338                                                       //all pairs not cotaining z
 339                                                       //will remain unchanged
 340                if(newL<l[s])
 341                   {
 342                    l[s]=newL;
 343                    st[s]=i;
 344                    ed[s]=z;
 345                   }
 346             }
 347          }
 348         free(ord);
 349         /*Rprintf("new tree\n");
 350         for(i=0;i<2*n-3;i++)
 351         {
 352         Rprintf("%i->%i of length:%f \n",ed1[i],ed2[i],edLen[i]);
 353         }
 354         Rprintf("end new tree\n");*/
 355       }
 356    //for(i=0;i<=numEdges;i++){ed1[i]++;ed2[i]++;}
 357  }