1 /* treePop.c 2011-10-11 */
3 /* Copyright 2011 Andrei-Alin Popescu */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
15 while (a[i] == 0) i++; /* count number of elements = 0 */
18 while ((a[i] & (1 << b)) == 0) b--;
23 short count_bits(uint8_t n)
25 short c; /* c accumulates the total bits set in v */
27 n &= n - 1; /* clear the least significant bit set */
31 /* Print n as a binary number */
32 /*void printbits(uint8_t n)
35 i = 1 << (sizeof(n) * 8 - 1);
38 if (n & i) Rprintf("1");
44 uint8_t* setdiff(uint8_t* x, uint8_t *y, int nrow) //x-y
47 uint8_t* ret = (uint8_t*)R_alloc(nrow, sizeof(uint8_t));
48 for (i = 0; i < nrow; i++) {
49 uint8_t tmp = (~y[i]);
51 /* Rprintf("x[%i]=",i);
61 ret[i] = (x[i] & tmp);
66 void treePop(int* splits, double* w,int* ncolp,int* np, int* ed1, int* ed2, double* edLen)
70 int nrow=ceil(n/(double)8);
71 uint8_t split[nrow][ncol];
74 /*Rprintf("n=%i nrow=%i ncol=%i\n",n,nrow,ncol);
80 Rprintf("%i ",splits[i*nrow+j]);
89 split[j][i]=(uint8_t)splits[i*nrow+j];
92 /*Rprintf("short-ed\n");
97 printbits(split[i][j]);
103 uint8_t vlabs[2*n-1][nrow];
104 for(i=0;i<nrow-1;i++)
108 vlabs[n+1][nrow-1]=~((uint8_t)(pow(2,8-(n%8))-1));
110 int bitt_count[ncol];
111 uint8_t msk=~((uint8_t)(pow(2,8-(n%8))-1));//mask out trailing bits
116 for(j=0;j<nrow-1;j++)
117 { //Rprintf("countbits(%i)=%i ",split[j][i],count_bits(split[j][i]));
118 sum+=count_bits(split[j][i]);
120 uint8_t bt=split[nrow-1][i];
122 //Rprintf("countbits(%i)=%i ",bt,count_bits(bt));
124 // Rprintf("bit_count[%i]=%i ",i,sum);
130 split[j][i]=~split[j][i];
132 split[nrow-1][i]&=msk;
143 for(i=0;i<ncol-1;i++)
145 for(j=i+1;j<ncol;j++)
147 if(bitt_count[i]<bitt_count[j])
150 bitt_count[i]=bitt_count[j];
162 //Rprintf("split %i\n",i);
164 for(ii=0;ii<nrow;ii++)
165 {//copy split into sp
166 sp[ii]=split[ii][ind[i]];
168 //search for node whose labellings are a superset of the current split
169 for(j=n+1;j<=nNodes;j++)
172 for(ii=0;ii<nrow;ii++)
173 {//copy vertex labeling into vl
176 int sw=0;//print current split
177 for(ii=0;ii<nrow;ii++)
179 //Rprintf("sp[%i]= ",ii);
182 //Rprintf("\n");//print current label
183 for(ii=0;ii<nrow;ii++)
185 // Rprintf("vl[%i]= ",ii);
186 // printbits(vl[ii]);
189 uint8_t* sd=setdiff(sp,vl,nrow);
190 //print the setdifference
191 for(ii=0;ii<nrow/*-1*/;ii++)
192 { //Rprintf("sd[%i]=%i ",ii,sd[ii]);
198 //Rprintf("sd[%i]=%i ",nrow-1,sd[nrow-1]);
199 if(sd[nrow-1]!=0)sw=1;
200 if(sw==0)//setdiff==0, so we split vertex j
201 { // Rprintf("vertex %i",j);
205 // Rprintf("bitt_count[%i]=%i\n",i,bitt_count[i]);
206 if(bitt_count[i]>=2)//if not trivial split
214 // Rprintf("gn=%i\n",gn);
216 edLen[numEdges]=w[ind[i]];
218 uint8_t* sdd=setdiff(vl,sp,nrow);
219 for(ii=0;ii<nrow;ii++)//label current vertex byset difference
220 //and new vertex by actual split
222 vlabs[j][ii]=sdd[ii];
223 vlabs[gn][ii]=sp[ii];
226 // Rprintf("new v\n");
228 for(ii=1;ii<=nNodes;ii++)
229 {//Rprintf("node %i : ",ii);
230 for(jj=0;jj<nrow;jj++)
232 //printbits(vlabs[ii][jj]);