src/treePop.c

   1 /* treePop.c    2011-10-11 */
   2
   3 /* Copyright 2011 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 <math.h>
   9 #include <R.h>
  10 #include <stdint.h>
  11
  12 int lsb(char* a)
  13 {
  14         int i = 0;
  15         while (a[i] == 0) i++; /* count number of elements = 0 */
  16
  17         int b = 7;
  18         while ((a[i] & (1 << b)) == 0) b--;
  19
  20         return i*8 + (8 - b);
  21 }
  22
  23 short count_bits(uint8_t n)
  24 {
  25         short c; /* c accumulates the total bits set in v */
  26         for (c = 0; n; c++)
  27                 n &= n - 1; /* clear the least significant bit set */
  28         return c;
  29 }
  30
  31 /* Print n as a binary number */
  32 /*void printbits(uint8_t n)
  33 {
  34         uint8_t i;
  35         i = 1 << (sizeof(n) * 8 - 1);
  36
  37         while (i > 0) {
  38                 if (n & i) Rprintf("1");
  39                 else Rprintf("0");
  40                 i >>= 1;
  41         }
  42 }*/
  43
  44 uint8_t* setdiff(uint8_t* x, uint8_t *y, int nrow) //x-y
  45 {
  46         int i = 0;
  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]);
  50
  51 /*              Rprintf("x[%i]=",i);
  52                 printbits(x[i]);
  53                 Rprintf("\n");
  54                 Rprintf("y[%i]=",i);
  55                 printbits(y[i]);
  56                 Rprintf("\n");
  57                 Rprintf("tmp=\n");
  58                 printbits(tmp);
  59                 Rprintf("\n"); */
  60
  61                 ret[i] = (x[i] & tmp);
  62         }
  63         return ret;
  64 }
  65
  66 void treePop(int* splits, double* w,int* ncolp,int* np, int* ed1, int* ed2, double* edLen)
  67   {
  68     int n=*np;
  69     int ncol=*ncolp;
  70     int nrow=ceil(n/(double)8);
  71     uint8_t split[nrow][ncol];
  72     int i=0,j=0;
  73
  74     /*Rprintf("n=%i nrow=%i ncol=%i\n",n,nrow,ncol);
  75     Rprintf("got\n");
  76     for(i=0;i<ncol;i++)
  77     {
  78      for(j=0;j<nrow;j++)
  79        {
  80           Rprintf("%i ",splits[i*nrow+j]);
  81        }
  82      Rprintf("\n");
  83     }*/
  84
  85     for(i=0;i<ncol;i++)
  86     {
  87      for(j=0;j<nrow;j++)
  88        {
  89          split[j][i]=(uint8_t)splits[i*nrow+j];
  90        }
  91     }
  92     /*Rprintf("short-ed\n");
  93     for(i=0;i<nrow;i++)
  94     {
  95       for(j=0;j<ncol;j++)
  96        {
  97          printbits(split[i][j]);
  98          Rprintf("\n");
  99        }
 100       Rprintf("\n");
 101     }*/
 102
 103    uint8_t vlabs[2*n-1][nrow];
 104    for(i=0;i<nrow-1;i++)
 105     {
 106        vlabs[n+1][i]=255;
 107     }
 108    vlabs[n+1][nrow-1]=~((uint8_t)(pow(2,8-(n%8))-1));
 109
 110    int bitt_count[ncol];
 111    uint8_t msk=~((uint8_t)(pow(2,8-(n%8))-1));//mask out trailing bits
 112    //printbits(msk);
 113    for(i=0;i<ncol;i++)
 114     {
 115       int sum=0;
 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]);
 119        }
 120       uint8_t bt=split[nrow-1][i];
 121       bt&=msk;
 122       //Rprintf("countbits(%i)=%i ",bt,count_bits(bt));
 123       sum+=count_bits(bt);
 124      // Rprintf("bit_count[%i]=%i ",i,sum);
 125       //Rprintf("\n");
 126       if(sum>n/2)
 127        {
 128          for(j=0;j<nrow;j++)
 129           {
 130              split[j][i]=~split[j][i];
 131           }
 132          split[nrow-1][i]&=msk;
 133          sum=n-sum;
 134        }
 135       bitt_count[i]=sum;
 136     }
 137    int ind[ncol];
 138    for(i=0;i<ncol;i++)
 139      {
 140        ind[i]=i;
 141      }
 142
 143    for(i=0;i<ncol-1;i++)
 144      {
 145        for(j=i+1;j<ncol;j++)
 146         {
 147           if(bitt_count[i]<bitt_count[j])
 148             { int aux;
 149               aux=bitt_count[i];
 150               bitt_count[i]=bitt_count[j];
 151               bitt_count[j]=aux;
 152               aux=ind[i];
 153               ind[i]=ind[j];
 154               ind[j]=aux;
 155             }
 156         }
 157      }
 158    int nNodes=n+1;
 159    int numEdges=0;
 160    for(i=0;i<ncol;i++)
 161     {  int ii=0;
 162        //Rprintf("split %i\n",i);
 163        uint8_t sp[nrow];
 164        for(ii=0;ii<nrow;ii++)
 165          {//copy split into sp
 166           sp[ii]=split[ii][ind[i]];
 167          }
 168       //search for node whose labellings are a superset of the current split
 169       for(j=n+1;j<=nNodes;j++)
 170        {
 171           uint8_t vl[nrow];
 172           for(ii=0;ii<nrow;ii++)
 173             {//copy vertex labeling into vl
 174               vl[ii]=vlabs[j][ii];
 175             }
 176          int sw=0;//print current split
 177          for(ii=0;ii<nrow;ii++)
 178            {
 179              //Rprintf("sp[%i]= ",ii);
 180             //printbits(sp[ii]);
 181            }
 182          //Rprintf("\n");//print current label
 183          for(ii=0;ii<nrow;ii++)
 184            {
 185             // Rprintf("vl[%i]= ",ii);
 186            //  printbits(vl[ii]);
 187            }
 188          //Rprintf("\n");
 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]);
 193              if(sd[ii]!=0)sw=1;
 194            }
 195         // Rprintf("\n");
 196
 197          sd[nrow-1]&=msk;
 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);
 202              ed1[numEdges]=j;
 203
 204              int gn=-1;
 205             // Rprintf("bitt_count[%i]=%i\n",i,bitt_count[i]);
 206              if(bitt_count[i]>=2)//if not trivial split
 207              {nNodes++;
 208               gn=nNodes;
 209              }
 210              else
 211              {
 212                gn=lsb(sp);
 213              }
 214             // Rprintf("gn=%i\n",gn);
 215              ed2[numEdges]=gn;
 216              edLen[numEdges]=w[ind[i]];
 217              numEdges++;
 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
 221                {
 222                  vlabs[j][ii]=sdd[ii];
 223                  vlabs[gn][ii]=sp[ii];
 224                }
 225              //print new labels
 226             // Rprintf("new v\n");
 227              int jj=0;
 228              for(ii=1;ii<=nNodes;ii++)
 229               {//Rprintf("node %i : ",ii);
 230                for(jj=0;jj<nrow;jj++)
 231                 {
 232                   //printbits(vlabs[ii][jj]);
 233                  // Rprintf(" ");
 234                 }
 235               // Rprintf("\n");
 236               }
 237              break;
 238           }
 239
 240        }
 241     }
 242   }