src/triangMtds.c

   1 /* triangMtds.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 #include "ape.h"
   9
  10 void triangMtds(double* d, int* np, int* ed1,int* ed2, double* edLen)
  11  {
  12     int n=*np;
  13     int k=0;
  14     int i=0;
  15     int j=0;
  16     int ij=-1;
  17     int m[n+1];
  18     int c[n+1];
  19     int s[n+1];
  20     int o[n+1];
  21     for(i=1;i<=n;i++)
  22     {m[i]=0;
  23      c[i]=i;
  24      s[i]=0;
  25      for(j=1;j<=n;j++)
  26       {
  27         if(i==j){m[i]++;continue;}
  28         if(d[give_index(i,j,n)]==-1)continue;
  29         m[i]++;
  30       }
  31     }
  32
  33    for(i=1;i<n;i++)
  34     for(j=i+1;j<=n;j++)
  35       {
  36         if(m[i]<m[j])
  37          {
  38             int aux=m[i];
  39             m[i]=m[j];
  40             m[j]=aux;
  41             aux=c[i];
  42             c[i]=c[j];
  43             c[j]=aux;
  44          }
  45       }
  46    ij=1;
  47    while(m[ij]==n){ij++;}
  48    for(i=ij;i<n;i++)
  49      {if(s[c[i]]==1)continue;
  50       for(j=i+1;j<=n;j++)
  51       {
  52        if(s[c[j]]==1)continue;
  53        if(d[give_index(c[i],c[j],n)]==-1)
  54          {s[c[j]]=1;
  55          }
  56       }
  57      }
  58
  59    for(i=1;i<=n;i++)
  60      {
  61        if(s[i]==0)
  62          {
  63           k++;
  64           o[k]=i;
  65          }
  66      }
  67    double* sub_d=(double*)R_alloc(k*(k - 1)/2, sizeof(double));
  68    ij=0;
  69    for(i=1;i<n;i++)
  70    {if(s[i]==1)continue;
  71     for(j=i+1;j<=n;j++)
  72      {
  73       if(s[j]==1)continue;
  74       sub_d[ij++]=d[give_index(i,j,n)];
  75      }
  76    }
  77
  78    triangMtd(sub_d,&k,ed1,ed2,edLen);
  79    for(i=0;i<2*k-3;i++)
  80      {
  81        if(ed1[i]>k)
  82         {
  83           ed1[i]+=(n-k);
  84         }
  85        if(ed2[i]>k)
  86         {
  87           ed2[i]+=(n-k);
  88         }
  89      }
  90    for(i=0;i<2*k-3;i++)
  91      {
  92       if(ed2[i]<=k)
  93        {
  94           ed2[i]=o[ed2[i]];
  95        }
  96      }
  97    for(i=1;i<=n;i++)
  98     {
  99       if(s[i]==0)continue;//take only leaves not in Y
 100       m[i]=0;
 101       for(j=1;j<=n;j++)
 102        {
 103          if(s[j]==1)continue;//take only leaves already in Y
 104          if(d[give_index(i,j,n)]==-1)continue;//igonore if distance unknown
 105          m[i]++;
 106        }
 107     }
 108  int numEdges=2*k-4;//0-based, so subtract 1
 109  //Rprintf("numEdge=%i",numEdges);
 110  int nv=(k-2)+n;
 111  while(k<n)
 112  {
 113    //now find element in X\Y such that it has most known distances to already
 114    //built tree until all leaves are added or we can not find a place to attach
 115    //the new element
 116    //s[i]=1 => i not added to tree
 117     int max=-1;
 118     int maxPos=-1;
 119     for(i=1;i<=n;i++)
 120     {
 121      if(s[i]==0)continue;
 122      if(m[i]>max)
 123       {
 124         max=m[i];
 125         maxPos=i;
 126       }
 127     }
 128    s[maxPos]=0;//mark maxPos as added
 129    //calculate new m values for leaves not added, i.e we just increment any
 130    //already present value by 1 if we know the distance between i and maxPos
 131    for(i=1;i<=n;i++)
 132     {
 133       if(s[i]==0)continue;
 134       if(d[give_index(i,maxPos,n)]==-1)continue;
 135       m[i]++;
 136     }
 137
 138    //find path to attach maxPos to, grow tree
 139         double minDist=1e50;
 140         int z=maxPos;
 141         int x=-1,y=-1;
 142         for(i=1;i<n;i++)
 143         {if(s[i]==1 || d[give_index(i,z,n)]==-1 || i==z)continue;
 144          for(j=i+1;j<=n;j++)
 145           {
 146             if(s[j]==1 || d[give_index(j,z,n)]==-1 || j==z)continue;
 147             double tDist=(d[give_index(i,z,n)]+d[give_index(j,z,n)]-d[give_index(i,j,n)])/2;
 148             if(tDist<minDist)
 149              {
 150                 minDist=tDist;
 151                 x=i;
 152                 y=j;
 153              }
 154           }
 155         }
 156         if(x==-1 || y==-1)
 157          {error("could not build tree from given distance matrix");
 158          }
 159         int* ord=getPathBetween(x,y,n,ed1,ed2,numEdges);
 160         /*for(i=1;i<=2*n-2;i++)
 161          {Rprintf("ord[%i]=%i ",i,ord[i]);
 162          }
 163         Rprintf("\n");*/
 164         //look for the edge on the path x to y to subdivide
 165
 166         int p=x;
 167         double sum=0;
 168         double prevSum=0;
 169         int aux=0;
 170
 171         int subdiv=-1;//index of edge to subdivide
 172         //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)]);
 173         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
 174        // 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);
 175         //point from x
 176         //Rprintf("Adding leaf %i, between %i and %i\n",z,x,y);
 177        // Rprintf("%i situated at a distance of %d from tree",z,minDist);
 178         int sw=0;
 179         //Rprintf("path between %i and %i\n",x,y);
 180         //int cc=0;
 181         while(p!=y && sum<lx)
 182           { //cc++;
 183             aux=p;
 184             //error("%i -> %i ",p,ord[p]);
 185             p=ord[p];
 186             prevSum=sum;
 187             for(i=0;i<=numEdges;i++)
 188               {
 189                 if((ed1[i]==aux && ed2[i]==p)||(ed2[i]==aux && ed1[i]==p))
 190                   {
 191                     if(ed1[i]==aux && ed2[i]==p){sw=1;}
 192                     subdiv=i;
 193                     sum+=edLen[i];
 194                   }
 195               }
 196             //if(cc>1000)error("failed to follow path between x=%i y=%i\n",x,y);
 197           }
 198
 199
 200         nv++;
 201         //subdivide subdiv with a node labelled nv
 202         //length calculation
 203
 204         int edd=ed2[subdiv];
 205         ed2[subdiv]=nv;
 206         edLen[subdiv]= (sw==1)?(lx-prevSum):(sum-lx);//check which 'half' of the
 207                                                      //path the leaf belongs to
 208                                                      //and updates accordingly
 209         //error("sum=%f, prevsum=%f\n",sum,prevSum);
 210         //error("lx-prevSum=%f, sum-lx=%f, minDist=%f",lx-prevSum,sum-lx,minDist);
 211         //Rprintf("adding %i on path %i %i, at distance %f from %i, and %f from tree\n",z,x,y,lx,x,minDist);
 212        // Rprintf("subdividing edge %i\n",subdiv);
 213         numEdges++;
 214         ed1[numEdges]=nv;
 215         ed2[numEdges]=edd;
 216         edLen[numEdges]= (sw==1)?(sum-lx):(lx-prevSum);
 217         numEdges++;
 218         edLen[numEdges]=minDist;
 219         ed1[numEdges]=nv;
 220         ed2[numEdges]=z;
 221    k++;
 222  }
 223  }