src/bionjs.c

   1 /* bionjs.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 bionjs(double *D, int *N, int *edge1, int *edge2, double *edge_length,int* fsS)
  11 {       //assume missing values are denoted by -1
  12         double *S,*R , *v,*new_v, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
  13         int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
  14         /*for(i=0;i<n*(n-1)/2;i++)
  15           {if(isNA(D[i])){D[i]=-1;}
  16           }*/
  17         int *s;//s contains |Sxy|, which is all we need for agglomeration
  18         double *newR;
  19         int *newS;
  20         int fS=*fsS;
  21
  22         R = &Sdist;
  23         new_dist = &Ndist;
  24         otu_label = &o_l;
  25         n = *N;
  26         cur_nod = 2*n - 2;
  27
  28         R = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  29         v = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  30         new_v = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  31         S = (double*)R_alloc(n + 1, sizeof(double));
  32         newR = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  33         new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  34         otu_label = (int*)R_alloc(n + 1, sizeof(int));
  35         s = (int*)R_alloc(n*(n - 1)/2, sizeof(int));
  36         newS = (int*)R_alloc(n*(n - 1)/2, sizeof(int));
  37
  38         for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
  39
  40         k = 0;
  41         //populate the v matrix
  42         for(i=1;i<n;i++)
  43          for(j=i+1;j<=n;j++)
  44           {
  45            v[give_index(i,j,n)]=D[give_index(i,j,n)];
  46           }
  47         //compute Sxy and Rxy
  48         for(i=0;i<n*(n-1)/2;i++)
  49           {newR[i]=0;
  50            newS[i]=0;
  51            s[i]=0;
  52            R[i]=0;
  53           }
  54
  55         for(i=1;i<n;i++)
  56          for(j=i+1;j<=n;j++)
  57          {//algorithm assumes i,j /in Sij, so skip pair if it is not known
  58           if(D[give_index(i,j,n)]==-1)
  59             {
  60               continue;
  61             }
  62           for(k=1;k<=n;k++)
  63            {//ij is the pair for which we compute
  64             //skip k if we do not know the distances between it and i AND j
  65
  66              if(k==i || k==j)
  67                {
  68                   /* added 2012-04-02: */
  69                   if(i!=k)R[give_index(i,j,n)]+=D[give_index(i,k,n)];
  70                   if(j!=k)R[give_index(i,j,n)]+=D[give_index(j,k,n)];
  71                   /* end of addition */
  72                   s[give_index(i,j,n)]++;
  73                   continue;
  74                }
  75               if(D[give_index(i,k,n)]==-1 || D[give_index(j,k,n)]==-1)continue;
  76               //Rprintf("%i\n",k);
  77               s[give_index(i,j,n)]++;
  78               R[give_index(i,j,n)]+=D[give_index(i,k,n)];
  79               R[give_index(i,j,n)]+=D[give_index(j,k,n)];
  80            }
  81          }
  82
  83         /*for(i=1;i<n;i++)
  84                   {
  85                     for(j=i+1;j<=n;j++)
  86                       {
  87                         Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
  88                       }
  89                     Rprintf("\n");
  90                   }
  91
  92                 for(i=1;i<n;i++)
  93                   {
  94                     for(j=i+1;j<=n;j++)
  95                       {
  96                         Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
  97                       }
  98                     Rprintf("\n");
  99                   }*/
 100
 101         k=0;
 102         int sw=1;//if 1 then incomplete
 103         while (n > 3) {
 104
 105                 ij = 0;
 106                 for(i=1;i<n;i++)
 107                  for(j=i+1;j<=n;j++)
 108                   {newR[give_index(i,j,n)]=0;
 109                    newS[give_index(i,j,n)]=0;
 110                   }
 111                 smallest_S = -1e50;
 112                 if(sw==0)
 113                     for(i=1;i<=n;i++)
 114                        {S[i]=0;
 115                        }
 116
 117                 B=n-2;
 118                 if(sw==1)
 119                      {
 120                       choosePair(D,n,R,s,&sw,&OTU1,&OTU2,fS);
 121                      }
 122                  else{ //Rprintf("distance matrix is now complete\n");
 123                         for (i=1;i<=n;i++)
 124                          for(j=1;j<=n;j++)
 125                            {if(i==j)continue;
 126                              //Rprintf("give_index(%i,%i)=%i\n",i,j,give_index(i,j,n));
 127                              //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
 128                              S[i]+=D[give_index(i,j,n)];
 129                            }
 130                         B=n-2;
 131                         //Rprintf("n=%i,B=%f",n,B);
 132                         for (i = 1; i < n; i++) {
 133                          for (j = i + 1; j <= n; j++) {
 134                              //Rprintf("S[%i]=%f, S[%i]=%f, D[%i,%i]=%f, B=%f",i,S[i],j,S[j],i,j,D[give_index(i,j,n)],B);
 135                                 A=S[i]+S[j]-B*D[give_index(i,j,n)];
 136                                 //Rprintf("Q[%i,%i]=%f\n",i,j,A);
 137                                 if (A > smallest_S) {
 138                                         OTU1 = i;
 139                                         OTU2 = j;
 140                                         smallest_S = A;
 141                                         smallest = ij;
 142                                 }
 143                                 ij++;
 144                         }
 145                         }
 146                      }
 147
 148                 //Rprintf("agglomerating %i and %i, Q=%f \n",OTU1,OTU2,smallest_S);
 149
 150                 /*for(i=1;i<n;i++)
 151                   {
 152                     for(j=i+1;j<=n;j++)
 153                       {
 154                         Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
 155                       }
 156                     Rprintf("\n");
 157                   }
 158
 159                 for(i=1;i<n;i++)
 160                   {
 161                     for(j=i+1;j<=n;j++)
 162                       {
 163                         Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
 164                       }
 165                     Rprintf("\n");
 166                   }
 167
 168                 for(i=1;i<n;i++)
 169                   {
 170                     for(j=i+1;j<=n;j++)
 171                       {
 172                         Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
 173                       }
 174                     Rprintf("\n");
 175                   }*/
 176                 //update R and S, only if matrix still incomplete
 177                 if(sw==1)
 178                 for(i=1;i<n;i++)
 179                 {if(i==OTU1 || i==OTU2)continue;
 180                  for(j=i+1;j<=n;j++)
 181                   {if(j==OTU1 || j==OTU2)continue;
 182                     if(D[give_index(i,j,n)]==-1)continue;
 183                      if(D[give_index(i,OTU1,n)]!=-1 && D[give_index(j,OTU1,n)]!=-1)
 184                       {//OTU1 was considered for Rij, so now subtract
 185                        R[give_index(i,j,n)]-=(D[give_index(i,OTU1,n)]+D[give_index(j,OTU1,n)]);
 186                        s[give_index(i,j,n)]--;
 187                       }
 188                      if(D[give_index(i,OTU2,n)]!=-1 && D[give_index(j,OTU2,n)]!=-1)
 189                       {//OTU2 was considered for Rij, so now subtract
 190                        R[give_index(i,j,n)]-=(D[give_index(i,OTU2,n)]+D[give_index(j,OTU2,n)]);
 191                        s[give_index(i,j,n)]--;
 192                       }
 193                   }
 194                 }
 195
 196                 edge2[k] = otu_label[OTU1];
 197                 edge2[k + 1] = otu_label[OTU2];
 198                 edge1[k] = edge1[k + 1] = cur_nod;
 199
 200                 double sum=0;
 201                 double lamb=0;//the parameter used for matrix reduction
 202                 double lambSum=0;
 203                 for(i=1;i<=n;i++)
 204                 {if(i==OTU1 || i==OTU2)continue;
 205                  if(D[give_index(OTU1,i,n)]==-1 || D[give_index(OTU2,i,n)]==-1)continue;
 206                  sum+=(D[give_index(OTU1,i,n)]-D[give_index(OTU2,i,n)]);
 207                  lambSum+=(v[give_index(OTU2,i,n)]-v[give_index(OTU1,i,n)]);
 208                 }
 209                 //if we stil have incomplete distances
 210                 if(sw==1)
 211                 {
 212                  lamb=0.5+(1/(2*(s[give_index(OTU1,OTU2,n)]-2)*v[give_index(OTU2,OTU1,n)]))*lambSum;
 213                 }else{
 214                  lamb=0.5+(1/(2*(n-2)*v[give_index(OTU2,OTU1,n)]))*lambSum;
 215                      }
 216
 217                 //although s was updated above, s[otu1,otu2] has remained unchanged
 218                 //so it is safe to use it here
 219                 //if complete distanes, use N-2, else use S
 220                 int down=B;
 221                 if(sw==1){down=s[give_index(OTU1,OTU2,n)]-2;}
 222                 if(down<=0)
 223                   {error("distance information insufficient to construct a tree, leaves %i and %i isolated from tree",OTU1,OTU2);
 224                   }
 225                 //Rprintf("down=%f\n",B);
 226                 sum*=(1.0/(2*(down)));
 227                 //Rprintf("sum=%f\n",sum);
 228                 double dxy=D[give_index(OTU1,OTU2,n)]/2;
 229
 230                 //Rprintf("R[%i,%i]:%f \n",OTU1,OTU2,sum);
 231                 edge_length[k] = dxy+sum;//OTU1
 232                 //Rprintf("l1:%f \n",edge_length[k]);
 233                 edge_length[k + 1] = dxy-sum;//OTU2
 234                 //Rprintf("l2:%f \n",edge_length[k+1]);
 235                //no need to change distance matrix update for complete distance
 236                //case, as pairs will automatically fall in the right cathegory
 237
 238                 //OTU1=x, OTU2=y from formulas
 239                 A = D[give_index(OTU1,OTU2,n)];
 240                 ij = 0;
 241                 for (i = 1; i <= n; i++) {
 242                         if (i == OTU1 || i == OTU2) continue;
 243                         if(D[give_index(OTU1,i,n)]!=-1 && D[give_index(OTU2,i,n)]!=-1)
 244                          {
 245                             new_dist[ij]= lamb*(D[give_index(OTU1,i,n)]-edge_length[k])+(1-lamb)*(D[give_index(OTU2,i,n)]-edge_length[k+1]);
 246                             new_v[ij]=lamb*v[give_index(OTU1,i,n)]+(1-lamb)*v[give_index(OTU2,i,n)]-lamb*(1-lamb)*v[give_index(OTU1,OTU2,n)];
 247                          }else{
 248                          if(D[give_index(OTU1,i,n)]!=-1)
 249                                 {
 250                                  new_dist[ij]=D[give_index(OTU1,i,n)]-edge_length[k];
 251                                  new_v[ij]=v[give_index(OTU1,i,n)];
 252                                 }else{
 253                                       if(D[give_index(OTU2,i,n)]!=-1)
 254                                         {
 255                                             new_dist[ij]=D[give_index(OTU2,i,n)]-edge_length[k+1];
 256                                             new_v[ij]=v[give_index(OTU2,i,n)];
 257                                         }else{new_dist[ij]=-1;new_v[ij]=-1;}
 258                                      }
 259                               }
 260
 261                         ij++;
 262                 }
 263
 264                 for (i = 1; i < n; i++) {
 265                         if (i == OTU1 || i == OTU2) continue;
 266                         for (j = i + 1; j <= n; j++) {
 267                                 if (j == OTU1 || j == OTU2) continue;
 268                                 new_dist[ij] = D[DINDEX(i, j)];
 269                                 new_v[ij]=v[give_index(i,j,n)];
 270                                 ij++;
 271                         }
 272                 }
 273
 274                 /*for(i=1;i<n-1;i++)
 275                 {
 276                   for(j=i+1;j<=n-1;j++)
 277                    {Rprintf("%f ",new_dist[give_index(i,j,n-1)]);
 278                    }
 279                   Rprintf("\n");
 280                 }*/
 281                 //compute Rui, only if distance matrix is still incomplete
 282                 ij=0;
 283                 if(sw==1)
 284                 for(i=2;i<n;i++)
 285                   {
 286                    ij++;
 287                    if(new_dist[give_index(i,1,n-1)]==-1)continue;
 288
 289                    for(j=1;j<n;j++)
 290                      {
 291                        if(j==1 || j==i)
 292                        {
 293                          /* added 2012-04-02 */
 294                          if(i!=j)newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)];
 295                          if(1!=j)newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)];
 296                          /* end of addition */
 297                          newS[give_index(1,i,n-1)]++;
 298                          continue;
 299                        }
 300                        if(new_dist[give_index(i,j,n-1)]!=-1 && new_dist[give_index(1,j,n-1)]!=-1)
 301                         {
 302                           newS[give_index(1,i,n-1)]++;
 303                           newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)];
 304                           newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)];
 305                         }
 306                      }
 307                   }
 308                 //fill in the rest of R and S, again only if distance matrix still
 309                 //incomplete
 310                 if(sw==1)
 311                 for(i=1;i<n;i++)
 312                 {if(i==OTU1 || i==OTU2)continue;
 313                  for(j=i+1;j<=n;j++)
 314                   {if(j==OTU1 || j==OTU2)continue;
 315                    newR[ij]=R[give_index(i,j,n)];
 316                    newS[ij]=s[give_index(i,j,n)];
 317                    ij++;
 318                   }
 319                 }
 320                 //update newR and newS with the new taxa, again only if distance
 321                 //matrix is still incomplete
 322                 if(sw==1)
 323                 for(i=2;i<n-1;i++)
 324                 {if(new_dist[give_index(1,i,n-1)]==-1)continue;
 325                  for(j=i+1;j<=n-1;j++)
 326                   {if(new_dist[give_index(1,j,n-1)]==-1)continue;
 327                    if(new_dist[give_index(i,j,n-1)]==-1)continue; /* added 2012-04-02 */
 328                    newR[give_index(i,j,n-1)]+=(new_dist[give_index(1,i,n-1)]+new_dist[give_index(1,j,n-1)]);
 329                    newS[give_index(i,j,n-1)]++;
 330                   }
 331                 }
 332                 /* compute the branch lengths */
 333
 334
 335
 336                 /* update before the next loop
 337                    (we are sure that OTU1 < OTU2) */
 338                 if (OTU1 != 1)
 339                         for (i = OTU1; i > 1; i--)
 340                                 otu_label[i] = otu_label[i - 1];
 341                 if (OTU2 != n)
 342                         for (i = OTU2; i < n; i++)
 343                                 otu_label[i] = otu_label[i + 1];
 344                 otu_label[1] = cur_nod;
 345
 346
 347
 348                 n--;
 349                 for (i = 0; i < n*(n - 1)/2; i++)
 350                   {
 351                     D[i] = new_dist[i];
 352                     v[i] = new_v[i];
 353                     if(sw==1)
 354                        {
 355                         R[i] = newR[i];
 356                         s[i] = newS[i];
 357                        }
 358                   }
 359                 cur_nod--;
 360                 k = k + 2;
 361         }
 362         int dK=0;//number of known distances in final distance matrix
 363         int iUK=-1;//index of unkown distance, if we have one missing distance
 364         int iK=-1;//index of only known distance, only needed if dK==1
 365         for (i = 0; i < 3; i++) {
 366                 edge1[*N*2 - 4 - i] = cur_nod;
 367                 edge2[*N*2 - 4 - i] = otu_label[i + 1];
 368                 if(D[i]!=-1){dK++;iK=i;}else{iUK=i;}
 369         }
 370         if(dK==2)
 371          {//if two distances are known: assume our leaves are x,y,z, d(x,z) unknown
 372           //and edge weights of three edges are a,b,c, then any b,c>0 that
 373           //satisfy c-b=d(y,z)-d(x,y) a+c=d(y,z) are good edge weights, but for
 374           //simplicity we assume a=c if d(yz)<d(xy) a=b otherwise, and after some
 375           //algebra we get that we can set the missing distance equal to the
 376           //maximum of the already present distances
 377             double max=-1e50;
 378           for(i=0;i<3;i++)
 379             {if(i==iUK)continue;
 380              if(D[i]>max)max=D[i];
 381             }
 382           D[iUK]=max;
 383          }
 384         if(dK==1)
 385          {//through similar motivation as above, if we have just one known distance
 386           //we set the other two distances equal to it
 387           for(i=0;i<3;i++)
 388             {if(i==iK)continue;
 389              D[i]=D[iK];
 390             }
 391          }
 392         if(dK==0)
 393          {//no distances are known, we just set them to 1
 394           for(i=0;i<3;i++)
 395            {D[i]=1;
 396            }
 397          }
 398         edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
 399         edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
 400         edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
 401 }
 402
 403