src/ewLasso.c

   1 /* ewLasso.c    2013-03-30 */
   2
   3 /* Copyright 2013 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 int isTripletCover(int nmb, int n, int** s, int stat, int sSoFar[n], int* a)//number of sides, number of leaves, sides, side under consideration, set so far, lasso
  11 {
  12         int ret=0;
  13         if(stat==nmb)return 1;
  14         int i=0;
  15
  16         for(i=1;i<=n;i++)
  17          {
  18         if(!s[stat][i])continue;//not in set
  19                 int sw=1, j;
  20
  21                 for(j=1;j<=n;j++)//check if all distances to previous candidates are present
  22                  {
  23                          if(!sSoFar[j])continue;//not in set so far
  24                          if(!a[i*(n+1)+j]){//if not, then i is not a good candidate for this side
  25                                    //Rprintf("failed to find distance between %i and %i, a[%i][%i]=%i \n",i,j,i,j,a[i*(n+1)+j]);
  26                                    sw=0;
  27                                }
  28                  }
  29
  30                 if(!sw)continue;//not all required distances are present
  31
  32                 sSoFar[i]=1;//try choosing i as representative for the side
  33                 ret+=isTripletCover(nmb,n,s,stat+1,sSoFar,a)>0?1:0;//see if, with i chosen, we can find leaves in other sides to satisfy the triplet cover condition
  34                 sSoFar[i]=0;
  35          }
  36   return ret;
  37 }
  38
  39 void ewLasso(double *D, int *N, int *e1, int *e2)
  40 {
  41         double *S, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
  42         int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
  43
  44         n=*N;
  45         int tCov=1;
  46         int* a = (int*)R_alloc((n+1)*(n+1), sizeof(int));//adjacency matrix of G_{\cL} graph
  47
  48         ij=0;
  49
  50         for(i=1;i<=n;i++)
  51          {
  52           for(j=1;j<=n;j++)
  53            {
  54         if(D[give_index(i,j,n)]==-1)//if missing value then no edge between pair of taxa (i,j) in G
  55                  {
  56                   a[i*(n+1)+j]=a[j*(n+1)+i]=0;
  57                  }
  58                 else
  59                  {
  60                   a[i*(n+1)+j]=a[j*(n+1)+i]=1;// otherwise edge between pair of taxa (i,j) in G
  61                  }
  62            }
  63          }
  64    //check for connectedness of G
  65
  66    int *q = (int*)R_alloc(2*n-1, sizeof(int));//BFS queue
  67    int *v = (int*)R_alloc(2*n-1, sizeof(int));//visited?
  68
  69    int p=0,u=1;//p-head of queue, u- position after last loaded element
  70
  71    for(i=1;i<=n;i++)v[i]=-1;
  72
  73    int stNBipartite=1, con=1, comp=1;
  74    int ini=1;
  75    ij=0;
  76    /*for(i=1;i<=n;i++)
  77     {
  78          for(j=1;j<=n;j++)
  79           {
  80                 Rprintf("a[%i][%i]=%i ",i,j,a[i*(n+1)+j]);
  81           }
  82          Rprintf("\n");
  83     }*/
  84
  85    while(comp)
  86    {
  87    q[p]=ini;
  88    v[ini]=1;
  89    comp=0;
  90    int stNBipartiteLoc=0;//check if current connected component is bipartite
  91    while(p<u)//BFS
  92     {
  93      int head=q[p];
  94          //Rprintf("head: %i\n",head);
  95          //Rprintf("unvisited neighbours: \n");
  96          for(i=1;i<=n;i++)
  97           {
  98                 if(i==head)continue;
  99         if(a[i*(n+1)+head]==0)continue;
 100                 if(v[i]==v[head])//same color as vertex from which we visit-> not bipartite
 101                   {
 102                           stNBipartiteLoc=1;
 103                   }
 104                 if(v[i]!=-1)continue;
 105                 //Rprintf("vertex %i \n",i);
 106                 q[u++]=i;
 107                 v[i]=1-v[head];
 108           }
 109          p++;
 110     }
 111     stNBipartite*=stNBipartiteLoc;//anding strngly-non-bipartite over all connected components
 112
 113
 114
 115         //check if all vertices have been visited
 116         for(int i=1;i<=n;i++)
 117          {
 118           if(v[i]==-1)
 119            {
 120                    comp=1;
 121                    p=0;
 122                    u=1;
 123                    ini=i;
 124                    con=0;
 125                    break;
 126            }
 127          }
 128    }
 129
 130    Rprintf("connected: %i\n",con);
 131    Rprintf("strongly non-bipartite: %i\n",stNBipartite);
 132
 133    //finally check if \cL is triplet cover of T
 134
 135    //adjencency matrix of tree, 1 to n are leaves
 136
 137    int *at= R_alloc((2*n-1)*(2*n-1), sizeof(int));
 138
 139    for(i=1;i<=2*n-2;i++)
 140     {
 141           for(j=1;j<=2*n-2;j++)at[i*(2*n-1)+j]=0;
 142     }
 143
 144    for(i=0;i<2*n-3;i++)
 145     {
 146                 //Rprintf("e1[%i]=%i e2[%i]=%i \n",i,e1[i],i,e2[i]);
 147                 at[e1[i]*(2*n-1)+e2[i]]=at[e2[i]*(2*n-1)+e1[i]]=1;
 148     }
 149
 150   /*for(i=1;i<2*n-1;i++)
 151     {
 152          for(j=1;j<2*n-1;j++)
 153           {
 154                 Rprintf("at[%i][%i]=%i ",i,j,at[i*(2*n-1)+j]);
 155           }
 156          Rprintf("\n");
 157     }*/
 158
 159    for(i=n+1;i<=2*n-2;i++)//for each interior vertex
 160     {
 161      for(j=1;j<2*n-1;j++)//reset queue and visited veectors
 162        {
 163                 v[j]=-1;
 164                 q[j]=0;
 165        }
 166
 167          v[i]=1;//'disconnect' graph at i
 168          int *l=(int*)R_alloc(2*n-2, sizeof(int));//vertices adjacent to i
 169
 170      int nmb=0;//number of found adjacent vertices of i
 171
 172          for(j=1;j<=2*n-2;j++)//find adjacent vertices
 173                  {
 174                if(at[i*(2*n-1)+j]==1)
 175                     {
 176                           l[nmb++]=j;
 177                     }
 178                  }
 179
 180          int** s=(int**)R_alloc(nmb,sizeof(int*));//set of leaves in each side, stored as presence/absence
 181
 182          for(j=0;j<nmb;j++)s[j]=(int*)R_alloc(n+1,sizeof(int));
 183
 184          for(j=0;j<nmb;j++)
 185           {
 186        for(int k=1;k<=n;k++)s[j][k]=0;
 187           }
 188
 189          /*Rprintf("for %i \n",i);
 190
 191          for(j=0;j<nmb;j++)Rprintf("l[%i]= %i ",j,l[j]);
 192
 193          Rprintf("\n");*/
 194
 195          for(j=0;j<nmb;j++)//for each adjancet vertex, find its leafset
 196            {
 197                  p=0;
 198                  u=1;
 199                  ini=l[j];
 200                  q[p]=ini;
 201                  v[ini]=1;
 202                  if(ini<=n)s[j][ini]=1;
 203          comp=0;
 204                  int nbr=0;
 205          while(p<u)//BFS
 206                         {
 207                                 int head=q[p];
 208                                 //Rprintf("head: %i\n",head);
 209                                 //Rprintf("unvisited neighbours: \n");
 210                         for(nbr=1;nbr<=2*n-1;nbr++)
 211                           {
 212                                 if(nbr==head)continue;
 213                                 if(at[nbr*(2*n-1)+head]==0)continue;
 214                                 if(v[nbr]!=-1)continue;
 215                                 //Rprintf("vertex %i \n",nbr);
 216                                 if(nbr<=n)s[j][nbr]=1;//if leaf, then set visited to true
 217                                 q[u++]=nbr;
 218                                 v[nbr]=1;
 219                           }
 220                         p++;
 221                   }
 222            }
 223
 224          /*Rprintf("sides for %i\n",i);
 225          for(j=0;j<nmb;j++)
 226           {
 227        int ii;
 228            for(ii=1;ii<=n;ii++)
 229              {
 230                         if(s[j][ii])Rprintf("%i ",ii);
 231              }
 232            Rprintf("\n");
 233           }*/
 234
 235          int* sSoFar= (int*)R_alloc(n+1,sizeof(int));
 236
 237          for(j=1;j<=n;j++)sSoFar[j]=0;
 238
 239          tCov*=isTripletCover(nmb,n,s,0,sSoFar,a)>0?1:0;
 240     }
 241    Rprintf("is triplet cover? %i \n",tCov);
 242 }