1 /* ewLasso.c 2013-03-30 */
3 /* Copyright 2013 Andrei-Alin Popescu */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
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
13 if(stat==nmb)return 1;
18 if(!s[stat][i])continue;//not in set
21 for(j=1;j<=n;j++)//check if all distances to previous candidates are present
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]);
30 if(!sw)continue;//not all required distances are present
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
39 void ewLasso(double *D, int *N, int *e1, int *e2)
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;
46 int* a = (int*)R_alloc((n+1)*(n+1), sizeof(int));//adjacency matrix of G_{\cL} graph
54 if(D[give_index(i,j,n)]==-1)//if missing value then no edge between pair of taxa (i,j) in G
56 a[i*(n+1)+j]=a[j*(n+1)+i]=0;
60 a[i*(n+1)+j]=a[j*(n+1)+i]=1;// otherwise edge between pair of taxa (i,j) in G
64 //check for connectedness of G
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?
69 int p=0,u=1;//p-head of queue, u- position after last loaded element
71 for(i=1;i<=n;i++)v[i]=-1;
73 int stNBipartite=1, con=1, comp=1;
80 Rprintf("a[%i][%i]=%i ",i,j,a[i*(n+1)+j]);
90 int stNBipartiteLoc=0;//check if current connected component is bipartite
94 //Rprintf("head: %i\n",head);
95 //Rprintf("unvisited neighbours: \n");
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
104 if(v[i]!=-1)continue;
105 //Rprintf("vertex %i \n",i);
111 stNBipartite*=stNBipartiteLoc;//anding strngly-non-bipartite over all connected components
115 //check if all vertices have been visited
116 for(int i=1;i<=n;i++)
130 Rprintf("connected: %i\n",con);
131 Rprintf("strongly non-bipartite: %i\n",stNBipartite);
133 //finally check if \cL is triplet cover of T
135 //adjencency matrix of tree, 1 to n are leaves
137 int *at= R_alloc((2*n-1)*(2*n-1), sizeof(int));
139 for(i=1;i<=2*n-2;i++)
141 for(j=1;j<=2*n-2;j++)at[i*(2*n-1)+j]=0;
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;
150 /*for(i=1;i<2*n-1;i++)
154 Rprintf("at[%i][%i]=%i ",i,j,at[i*(2*n-1)+j]);
159 for(i=n+1;i<=2*n-2;i++)//for each interior vertex
161 for(j=1;j<2*n-1;j++)//reset queue and visited veectors
167 v[i]=1;//'disconnect' graph at i
168 int *l=(int*)R_alloc(2*n-2, sizeof(int));//vertices adjacent to i
170 int nmb=0;//number of found adjacent vertices of i
172 for(j=1;j<=2*n-2;j++)//find adjacent vertices
174 if(at[i*(2*n-1)+j]==1)
180 int** s=(int**)R_alloc(nmb,sizeof(int*));//set of leaves in each side, stored as presence/absence
182 for(j=0;j<nmb;j++)s[j]=(int*)R_alloc(n+1,sizeof(int));
186 for(int k=1;k<=n;k++)s[j][k]=0;
189 /*Rprintf("for %i \n",i);
191 for(j=0;j<nmb;j++)Rprintf("l[%i]= %i ",j,l[j]);
195 for(j=0;j<nmb;j++)//for each adjancet vertex, find its leafset
202 if(ini<=n)s[j][ini]=1;
208 //Rprintf("head: %i\n",head);
209 //Rprintf("unvisited neighbours: \n");
210 for(nbr=1;nbr<=2*n-1;nbr++)
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
224 /*Rprintf("sides for %i\n",i);
230 if(s[j][ii])Rprintf("%i ",ii);
235 int* sSoFar= (int*)R_alloc(n+1,sizeof(int));
237 for(j=1;j<=n;j++)sSoFar[j]=0;
239 tCov*=isTripletCover(nmb,n,s,0,sSoFar,a)>0?1:0;
241 Rprintf("is triplet cover? %i \n",tCov);