src/njs.c

   1 /* njs.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 "ape.h"
   9
  10 int H(double t)
  11 {
  12         if (t >= 0) return 1;
  13         return 0;
  14 }
  15
  16 void choosePair(double* D,int n,double* R,int* s, int* sw, int* x, int* y,int fS)
  17 {
  18     int i=0,j=0,k=0;
  19     int sww=0;
  20     double cFS[fS];
  21     int iFS[fS];
  22     int jFS[fS];
  23     for(k=0;k<fS;k++)
  24               {cFS[k]=-1e50;
  25                iFS[k]=0;
  26                jFS[k]=0;
  27               }
  28     double max=-1e50;
  29     for (i = 1; i < n; i++)
  30       {
  31        for (j = i + 1; j <= n; j++)
  32          {if(D[give_index(i,j,n)]==-1){sww=1;continue;}
  33            //Rprintf("R[%i,%i]=%f\n",i,j,R[give_index(i,j,n)]);
  34            //Rprintf("s[%i,%i]=%i\n",i,j,s[give_index(i,j,n)]);
  35            //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
  36            int tr=0;
  37             double numb=((R[give_index(i,j,n)])/(s[give_index(i,j,n)]-2))-D[give_index(i,j,n)];
  38            //Rprintf("numb(%i,%i)=%f\n",i,j,numb);
  39            for(k=0;k<fS, cFS[k]>numb;k++);
  40            //Rprintf("k=%i ",k);
  41            for(tr=fS-1;tr>k;tr--)
  42              {cFS[tr]=cFS[tr-1];
  43               iFS[tr]=iFS[tr-1];
  44               jFS[tr]=jFS[tr-1];
  45              }
  46
  47             if(k<fS){cFS[k]=numb;
  48                      iFS[k]=i;
  49                      jFS[k]=j;
  50                     }
  51           }
  52       }
  53
  54     //no missing distances, just return the one with maximal Q-criterion
  55    if(sww==0){*x=iFS[0];*y=jFS[0];*sw=0;return;
  56              }
  57    /*for(k=0;k<fS;k++)
  58               {Rprintf("value[%i]=%f ",k,cFS[k]);
  59                Rprintf("i=%i ",iFS[k]);
  60                Rprintf("j=%i ",jFS[k]);
  61                Rprintf("\n");
  62               }*/
  63   //calculate N*(x,y)
  64   for(i=0;i<fS;i++)
  65    {
  66     double nb=nxy(iFS[i],jFS[i],n,D);
  67     if(nb>max){max=nb;}
  68     cFS[i]=nb;
  69    }
  70   /*Rprintf("all N*(x,y)\n");
  71   for(k=0;k<fS;k++)
  72               {Rprintf("value[%i]=%f ",k,cFS[k]);
  73                Rprintf("i=%i ",iFS[k]);
  74                Rprintf("j=%i ",jFS[k]);
  75                Rprintf("\n");
  76               }*/
  77   int dk=0;
  78   //shift the max N*xy to the front of the array
  79   for(i=0;i<fS;i++)
  80    {
  81       if(cFS[i]==max)
  82         {
  83           cFS[dk]=cFS[i];
  84           iFS[dk]=iFS[i];
  85           jFS[dk]=jFS[i];
  86           dk++;
  87         }
  88    }
  89   //if just one pair realises max N*xy, return it
  90   if(dk==1){*x=iFS[0];*y=jFS[0];return;}
  91   fS=dk;
  92  /*Rprintf("max n*xy values\n");
  93  for(k=0;k<fS;k++)
  94               {Rprintf("value[%i]=%f ",k,cFS[k]);
  95                Rprintf("i=%i ",iFS[k]);
  96                Rprintf("j=%i ",jFS[k]);
  97                Rprintf("\n");
  98               }*/
  99  max=-1e50;
 100  //on the fron of the array containing max N*xy values compute cxy
 101  for(i=0;i<fS;i++)
 102   {
 103     double nb=cxy(iFS[i],jFS[i],n,D);
 104     if(nb>max){max=nb;}
 105     cFS[i]=nb;
 106   }
 107  /*Rprintf("all c*xy values\n");
 108  for(k=0;k<fS;k++)
 109               {Rprintf("value[%i]=%f ",k,cFS[k]);
 110                Rprintf("i=%i ",iFS[k]);
 111                Rprintf("j=%i ",jFS[k]);
 112                Rprintf("\n");
 113               }*/
 114   //and again shift maximal C*xy values at the fron of the array
 115   dk=0;
 116   for(i=0;i<fS;i++)
 117    {
 118       if(cFS[i]==max)
 119         {
 120           cFS[dk]=cFS[i];
 121           iFS[dk]=iFS[i];
 122           jFS[dk]=jFS[i];
 123           dk++;
 124         }
 125    }
 126  //if just one C*xy with maximal value, return pair realising it
 127  if(dk==1){*x=iFS[0];*y=jFS[0];return;}
 128  fS=dk;
 129  /*Rprintf("max c*xy values\n");
 130  for(k=0;k<fS;k++)
 131               {Rprintf("value[%i]=%f ",k,cFS[k]);
 132                Rprintf("i=%i ",iFS[k]);
 133                Rprintf("j=%i ",jFS[k]);
 134                Rprintf("\n");
 135               }*/
 136  max=-1e50;
 137  //on the front of the array containing max C*xy compute m*xy
 138  for(i=0;i<fS;i++)
 139   {
 140     double nb=mxy(iFS[i],jFS[i],n,D);
 141     if(nb>max){max=nb;}
 142     cFS[i]=nb;
 143   }
 144  /*Rprintf("all m*xy values\n");
 145  for(k=0;k<fS;k++)
 146               {Rprintf("value[%i]=%f ",k,cFS[k]);
 147                Rprintf("i=%i ",iFS[k]);
 148                Rprintf("j=%i ",jFS[k]);
 149                Rprintf("\n");
 150               }*/
 151   //again shift maximal m*xy values to the fron of the array
 152   dk=0;
 153   for(i=0;i<fS;i++)
 154    {
 155       if(cFS[i]==max)
 156         {
 157           cFS[dk]=cFS[i];
 158           iFS[dk]=iFS[i];
 159           jFS[dk]=jFS[i];
 160           dk++;
 161         }
 162    }
 163  //if just one maximal value for m*xy return the pair realising it, found at 0
 164  if(dk==1){*x=iFS[0];*y=jFS[0];return;}
 165  fS=dk;
 166  /*Rprintf("max m*xy values\n");
 167  for(k=0;k<fS;k++)
 168               {Rprintf("value[%i]=%f ",k,cFS[k]);
 169                Rprintf("i=%i ",iFS[k]);
 170                Rprintf("j=%i ",jFS[k]);
 171                Rprintf("\n");
 172               }*/
 173  //and calculate cnxy on these values, but this time we do not shift, but simply
 174  //return the pair realising the maximum, stored at iPos
 175  max=-1e50;
 176  int iPos=-1;
 177  for(i=0;i<fS;i++)
 178   {
 179     double nb=cnxy(iFS[i],jFS[i],n,D);
 180     if(nb>max){max=nb;iPos=i;}
 181     cFS[i]=nb;
 182   }
 183  /*Rprintf("cN*xy\n");
 184  Rprintf("value[%i]=%f ",0,cFS[0]);
 185  Rprintf("i=%i ",iFS[0]);
 186  Rprintf("j=%i ",jFS[0]);
 187  Rprintf("\n");*/
 188  *x=iFS[iPos];*y=jFS[iPos];
 189 }
 190
 191 double cnxy(int x, int y, int n,double* D)
 192 {
 193     int i=0;
 194     int j=0;
 195     double nMeanXY=0;
 196     //Rprintf("cN[%i,%i]\n",x,y);
 197     for(i=1;i<=n;i++)
 198      {if(i==x || i==y)continue;
 199       for(j=1;j<=n;j++)
 200       {if(i==j)continue;
 201        if(j==y || j==x)continue;
 202        double n1=0;
 203        double n2=0;
 204        n1=D[give_index(i,x,n)];
 205        n2=D[give_index(j,y,n)];
 206        if(n1==-1 || n2==-1 || D[give_index(i,j,n)]==-1)continue;
 207        nMeanXY+=(n1+n2-D[give_index(x,y,n)]-D[give_index(i,j,n)]);
 208        //Rprintf("cnMeanXY after (%i,%i)=%f\n",i,j,nMeanXY);
 209       }
 210      }
 211     return nMeanXY;
 212 }
 213
 214 int mxy(int x,int y,int n,double* D)
 215 {
 216     int i=0;
 217     int mx[n+1];
 218     int my[n+1];
 219     for(i=1;i<=n;i++)
 220       {
 221         mx[i]=0;my[i]=0;
 222       }
 223     for(i=1;i<=n;i++)
 224       {
 225         if(i!=x && D[give_index(x,i,n)]==-1)
 226           {
 227             mx[i]=1;
 228           }
 229         if(i!=y && D[give_index(y,i,n)]==-1)
 230           {
 231             my[i]=1;
 232           }
 233       }
 234     /*for(i=1;i<=n;i++)
 235       {
 236         Rprintf("mx[%i]=%i ",i,mx[i]);
 237       }
 238     Rprintf("\n");
 239
 240     for(i=1;i<=n;i++)
 241       {
 242         Rprintf("my[%i]=%i ",i,my[i]);
 243       }
 244     Rprintf("\n");*/
 245
 246     int xmy=0;
 247     int ymx=0;
 248     for(i=1;i<=n;i++)
 249       {
 250         if(mx[i]==1 && my[i]==0)
 251           {
 252             xmy++;
 253           }
 254         if(my[i]==1 && mx[i]==0)
 255           {
 256             ymx++;
 257           }
 258       }
 259     //Rprintf("xmy=%i, ymx=%i, xmy+ymx=%i\n",xmy,ymx,xmy+ymx);
 260     return xmy+ymx;
 261 }
 262 double nxy(int x, int y, int n,double* D)
 263 {
 264     int sCXY=0;
 265     int i=0;
 266     int j=0;
 267     double nMeanXY=0;
 268     //Rprintf("N[%i,%i]\n",x,y);
 269     for(i=1;i<=n;i++)
 270      {if(i==x || i==y)continue;
 271       for(j=1;j<=n;j++)
 272       {if(i==j)continue;
 273        if(j==x || j==y)continue;
 274        double n1=0;
 275        double n2=0;
 276        n1=D[give_index(i,x,n)];
 277        n2=D[give_index(j,y,n)];
 278        if(n1==-1 || n2==-1 || D[give_index(i,j,n)]==-1)continue;
 279        sCXY++;
 280        //Rprintf("considered pair (%i,%i)\n",i,j);
 281        nMeanXY+=H(n1+n2-D[give_index(x,y,n)]-D[give_index(i,j,n)]);
 282        //Rprintf("nMeanXY after (%i,%i)=%f\n",i,j,nMeanXY);
 283       }
 284      }
 285     //Rprintf("sCXY=%i",sCXY);
 286     return nMeanXY/sCXY;
 287 }
 288
 289 int cxy(int x, int y, int n,double* D)
 290 {
 291     int sCXY=0;
 292     int i=0;
 293     int j=0;
 294
 295     for(i=1;i<=n;i++)
 296      {if(i==x || i==y)continue;
 297       for(j=1;j<=n;j++)
 298       {if(i==j)continue;
 299        if(j==x || j==y)continue;
 300        double n1=0;
 301        double n2=0;
 302        if(i!=x)n1=D[give_index(i,x,n)];
 303        if(j!=y)n2=D[give_index(j,y,n)];
 304        if(n1==-1 || n2==-1 || D[give_index(i,j,n)]==-1)continue;
 305        sCXY++;
 306       }
 307      }
 308     return sCXY;
 309 }
 310
 311 void njs(double *D, int *N, int *edge1, int *edge2, double *edge_length,int *fsS)
 312 {       //assume missing values are denoted by -1
 313         double *S,*R, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
 314         int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
 315         /*for(i=0;i<n*(n-1)/2;i++)
 316           {if(isNA(D[i])){D[i]=-1;}
 317           }*/
 318         int *s;//s contains |Sxy|, which is all we need for agglomeration
 319         double *newR;
 320         int *newS;
 321         int fS=*fsS;
 322
 323         R = &Sdist;
 324         new_dist = &Ndist;
 325         otu_label = &o_l;
 326         n = *N;
 327         cur_nod = 2*n - 2;
 328
 329         R = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
 330         S = (double*)R_alloc(n + 1, sizeof(double));
 331         newR = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
 332         new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
 333         otu_label = (int*)R_alloc(n + 1, sizeof(int));
 334         s = (int*)R_alloc(n*(n - 1)/2, sizeof(int));
 335         newS = (int*)R_alloc(n*(n - 1)/2, sizeof(int));
 336
 337         for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
 338
 339         k = 0;
 340
 341         //compute Sxy and Rxy
 342         for(i=0;i<n*(n-1)/2;i++)
 343           {newR[i]=0;
 344            newS[i]=0;
 345            s[i]=0;
 346            R[i]=0;
 347           }
 348
 349         for(i=1;i<n;i++)
 350          for(j=i+1;j<=n;j++)
 351          {//algorithm assumes i,j /in Sij, so skip pair if it is not known
 352           if(D[give_index(i,j,n)]==-1)
 353             {
 354               continue;
 355             }
 356           //Rprintf("for %i and %i :\n",i,j);
 357           for(k=1;k<=n;k++)
 358            {//ij is the pair for which we compute
 359             //skip k if we do not know the distances between it and i AND j
 360               if(k==i || k==j)
 361                {
 362                   s[give_index(i,j,n)]++;
 363                   continue;
 364                }
 365               if(D[give_index(i,k,n)]==-1 || D[give_index(j,k,n)]==-1)continue;
 366               //Rprintf("%i\n",k);
 367               s[give_index(i,j,n)]++;
 368               R[give_index(i,j,n)]+=D[give_index(i,k,n)];
 369               R[give_index(i,j,n)]+=D[give_index(j,k,n)];
 370            }
 371          }
 372
 373         /*for(i=1;i<n;i++)
 374                   {
 375                     for(j=i+1;j<=n;j++)
 376                       {
 377                         Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
 378                       }
 379                     Rprintf("\n");
 380                   }
 381
 382                 for(i=1;i<n;i++)
 383                   {
 384                     for(j=i+1;j<=n;j++)
 385                       {
 386                         Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
 387                       }
 388                     Rprintf("\n");
 389                   }*/
 390
 391         k=0;
 392         int sw=1;//if 1 then incomplete
 393         while (n > 3) {
 394                 ij = 0;
 395                 for(i=1;i<n;i++)
 396                  for(j=i+1;j<=n;j++)
 397                   {newR[give_index(i,j,n)]=0;
 398                    newS[give_index(i,j,n)]=0;
 399                   }
 400                 smallest_S = -1e50;
 401                 if(sw==0)
 402                     for(i=1;i<=n;i++)
 403                        {S[i]=0;
 404                        }
 405                 B=n-2;
 406                 if(sw==1)
 407                      {
 408                       choosePair(D,n,R,s,&sw,&OTU1,&OTU2,fS);
 409                      }
 410                  else{ //Rprintf("distance matrix is now complete\n");
 411                         for (i=1;i<=n;i++)
 412                          for(j=1;j<=n;j++)
 413                            {if(i==j)continue;
 414                              //Rprintf("give_index(%i,%i)=%i\n",i,j,give_index(i,j,n));
 415                              //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
 416                              S[i]+=D[give_index(i,j,n)];
 417                            }
 418                         B=n-2;
 419                         //Rprintf("n=%i,B=%f",n,B);
 420                         for (i = 1; i < n; i++) {
 421                          for (j = i + 1; j <= n; j++) {
 422                              //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);
 423                                 A=S[i]+S[j]-B*D[give_index(i,j,n)];
 424                                 //Rprintf("Q[%i,%i]=%f\n",i,j,A);
 425                                 if (A > smallest_S) {
 426                                         OTU1 = i;
 427                                         OTU2 = j;
 428                                         smallest_S = A;
 429                                         smallest = ij;
 430                                 }
 431                                 ij++;
 432                         }
 433                         }
 434                      }
 435
 436                 /*Rprintf("agglomerating %i and %i, Q=%f \n",OTU1,OTU2,smallest_S);
 437
 438                 for(i=1;i<n;i++)
 439                   {
 440                     for(j=i+1;j<=n;j++)
 441                       {
 442                         Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
 443                       }
 444                     Rprintf("\n");
 445                   }
 446
 447                 for(i=1;i<n;i++)
 448                   {
 449                     for(j=i+1;j<=n;j++)
 450                       {
 451                         Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
 452                       }
 453                     Rprintf("\n");
 454                   }
 455
 456                 for(i=1;i<n;i++)
 457                   {
 458                     for(j=i+1;j<=n;j++)
 459                       {
 460                         Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
 461                       }
 462                     Rprintf("\n");
 463                   }*/
 464                 //update Rxy and Sxy, only if matrix still incomplete
 465                 if(sw==1)
 466                 for(i=1;i<n;i++)
 467                 {if(i==OTU1 || i==OTU2)continue;
 468                  for(j=i+1;j<=n;j++)
 469                   {if(j==OTU1 || j==OTU2)continue;
 470                     if(D[give_index(i,j,n)]==-1)continue;
 471                      if(D[give_index(i,OTU1,n)]!=-1 && D[give_index(j,OTU1,n)]!=-1)
 472                       {//OTU1 was considered for Rij, so now subtract
 473                        R[give_index(i,j,n)]-=(D[give_index(i,OTU1,n)]+D[give_index(j,OTU1,n)]);
 474                        s[give_index(i,j,n)]--;
 475                       }
 476                      if(D[give_index(i,OTU2,n)]!=-1 && D[give_index(j,OTU2,n)]!=-1)
 477                       {//OTU2 was considered for Rij, so now subtract
 478                        R[give_index(i,j,n)]-=(D[give_index(i,OTU2,n)]+D[give_index(j,OTU2,n)]);
 479                        s[give_index(i,j,n)]--;
 480                       }
 481                   }
 482                 }
 483
 484                 edge2[k] = otu_label[OTU1];
 485                 edge2[k + 1] = otu_label[OTU2];
 486                 edge1[k] = edge1[k + 1] = cur_nod;
 487
 488                 /* get the distances between all OTUs but the 2 selected ones
 489                    and the latter:
 490                    a) get the sum for both
 491                    b) compute the distances for the new OTU */
 492
 493                 double sum=0;
 494
 495                 for(i=1;i<=n;i++)
 496                  {if(i==OTU1 || i==OTU2)continue;
 497                  if(D[give_index(OTU1,i,n)]==-1 || D[give_index(OTU2,i,n)]==-1)continue;
 498                  sum+=(D[give_index(OTU1,i,n)]-D[give_index(OTU2,i,n)]);
 499                  }
 500                 //although s was updated above, s[otu1,otu2] has remained unchanged
 501                 //so it is safe to use it here
 502                 //if complete distanes, use N-2, else use S
 503                 int down=B;
 504                 if(sw==1){down=s[give_index(OTU1,OTU2,n)]-2;}
 505                 if(down==0)
 506                   {error("distance information insufficient to construct a tree, leaves %i and %i isolated from tree",OTU1,OTU2);
 507                   }
 508                 //Rprintf("down=%i\n",down);
 509                 sum*=(1.0/(2*(down)));
 510                 //Rprintf("sum=%f\n",sum);
 511                 double dxy=D[give_index(OTU1,OTU2,n)]/2;
 512
 513                 //Rprintf("R[%i,%i]:%f \n",OTU1,OTU2,sum);
 514                 edge_length[k] = dxy+sum;//OTU1
 515                 //Rprintf("l1:%f \n",edge_length[k]);
 516                 edge_length[k + 1] = dxy-sum;//OTU2
 517                 //Rprintf("l2:%f \n",edge_length[k+1]);
 518                //no need to change distance matrix update for complete distance
 519                //case, as pairs will automatically fall in the right cathegory
 520                 A = D[give_index(OTU1,OTU2,n)];
 521                 ij = 0;
 522                 for (i = 1; i <= n; i++) {
 523                         if (i == OTU1 || i == OTU2) continue;
 524                         if(D[give_index(OTU1,i,n)]!=-1 && D[give_index(OTU2,i,n)]!=-1)
 525                          {
 526                             new_dist[ij]=0.5*(D[give_index(OTU1,i,n)]-edge_length[k]+D[give_index(OTU2,i,n)]-edge_length[k+1]);
 527                          }else{
 528                          if(D[give_index(OTU1,i,n)]!=-1)
 529                                 {
 530                                  new_dist[ij]=D[give_index(OTU1,i,n)]-edge_length[k];
 531                                 }else{
 532                                       if(D[give_index(OTU2,i,n)]!=-1)
 533                                         {
 534                                             new_dist[ij]=D[give_index(OTU2,i,n)]-edge_length[k+1];
 535                                         }else{new_dist[ij]=-1;}
 536                                      }
 537                               }
 538
 539                         ij++;
 540                 }
 541
 542                 for (i = 1; i < n; i++) {
 543                         if (i == OTU1 || i == OTU2) continue;
 544                         for (j = i + 1; j <= n; j++) {
 545                                 if (j == OTU1 || j == OTU2) continue;
 546                                 new_dist[ij] = D[DINDEX(i, j)];
 547                                 ij++;
 548                         }
 549                 }
 550
 551                 /*for(i=1;i<n-1;i++)
 552                 {
 553                   for(j=i+1;j<=n-1;j++)
 554                    {Rprintf("%f ",new_dist[give_index(i,j,n-1)]);
 555                    }
 556                   Rprintf("\n");
 557                 }*/
 558                 //compute Rui, only if distance matrix is still incomplete
 559                 ij=0;
 560                 if(sw==1)
 561                 for(i=2;i<n;i++)
 562                   {
 563                    ij++;
 564                    if(new_dist[give_index(i,1,n-1)]==-1)continue;
 565
 566                    for(j=1;j<n;j++)
 567                      { if(j==1 || j==i)
 568                        {
 569                          newS[give_index(1,i,n-1)]++;
 570                          continue;
 571                        }
 572                        if(new_dist[give_index(i,j,n-1)]!=-1 && new_dist[give_index(1,j,n-1)]!=-1)
 573                         {
 574                           newS[give_index(1,i,n-1)]++;
 575                           newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)];
 576                           newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)];
 577                         }
 578                      }
 579                   }
 580                 //fill in the rest of R and S, again only if distance matrix still
 581                 //incomplete
 582                 for(i=1;i<n;i++)
 583                 {if(i==OTU1 || i==OTU2)continue;
 584                  for(j=i+1;j<=n;j++)
 585                   {if(j==OTU1 || j==OTU2)continue;
 586                    newR[ij]=R[give_index(i,j,n)];
 587                    newS[ij]=s[give_index(i,j,n)];
 588                    ij++;
 589                   }
 590                 }
 591                 //update newR and newS with the new taxa, again only if distance
 592                 //matrix is still incomplete
 593                 if(sw==1)
 594                 for(i=2;i<n-1;i++)
 595                 {if(new_dist[give_index(1,i,n-1)]==-1)continue;
 596                  for(j=i+1;j<=n-1;j++)
 597                   {if(new_dist[give_index(1,j,n-1)]==-1)continue;
 598                     if(new_dist[give_index(i,j,n-1)]==-1)continue;
 599                    newR[give_index(i,j,n-1)]+=(new_dist[give_index(1,i,n-1)]+new_dist[give_index(1,j,n-1)]);
 600                    newS[give_index(i,j,n-1)]++;
 601                   }
 602                 }
 603
 604                 /* update before the next loop
 605                    (we are sure that OTU1 < OTU2) */
 606                 if (OTU1 != 1)
 607                         for (i = OTU1; i > 1; i--)
 608                                 otu_label[i] = otu_label[i - 1];
 609                 if (OTU2 != n)
 610                         for (i = OTU2; i < n; i++)
 611                                 otu_label[i] = otu_label[i + 1];
 612                 otu_label[1] = cur_nod;
 613
 614                 n--;
 615                 for (i = 0; i < n*(n - 1)/2; i++)
 616                   {
 617                     D[i] = new_dist[i];
 618                     if(sw==1)
 619                        {
 620                         R[i] = newR[i];
 621                         s[i] = newS[i];
 622                        }
 623                   }
 624                 cur_nod--;
 625                 k = k + 2;
 626         }
 627         int dK=0;//number of known distances in final distance matrix
 628         int iUK=-1;//index of unkown distance, if we have one missing distance
 629         int iK=-1;//index of only known distance, only needed if dK==1
 630         for (i = 0; i < 3; i++) {
 631                 edge1[*N*2 - 4 - i] = cur_nod;
 632                 edge2[*N*2 - 4 - i] = otu_label[i + 1];
 633                 if(D[i]!=-1){dK++;iK=i;}else{iUK=i;}
 634         }
 635         if(dK==2)
 636          {//if two distances are known: assume our leaves are x,y,z, d(x,z) unknown
 637           //and edge weights of three edges are a,b,c, then any b,c>0 that
 638           //satisfy c-b=d(y,z)-d(x,y) a+c=d(y,z) are good edge weights, but for
 639           //simplicity we assume a=c if d(yz)<d(xy) a=b otherwise, and after some
 640           //algebra we get that we can set the missing distance equal to the
 641           //maximum of the already present distances
 642             double max=-1e50;
 643           for(i=0;i<3;i++)
 644             {if(i==iUK)continue;
 645              if(D[i]>max)max=D[i];
 646             }
 647           D[iUK]=max;
 648          }
 649         if(dK==1)
 650          {//through similar motivation as above, if we have just one known distance
 651           //we set the other two distances equal to it
 652           for(i=0;i<3;i++)
 653             {if(i==iK)continue;
 654              D[i]=D[iK];
 655             }
 656          }
 657         if(dK==0)
 658          {//no distances are known, we just set them to 1
 659           for(i=0;i<3;i++)
 660            {D[i]=1;
 661            }
 662          }
 663         edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
 664         edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
 665         edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
 666 }
 667