src/mvr.c

   1 /* mvr.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 void mvr(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_length)
  11 {
  12         double *S, Sdist, *new_v, 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
  15         S = &Sdist;
  16         new_dist = &Ndist;
  17         otu_label = &o_l;
  18         n = *N;
  19         cur_nod = 2*n - 2;
  20
  21         S = (double*)R_alloc(n + 1, sizeof(double));
  22         new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  23         new_v = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  24         otu_label = (int*)R_alloc(n + 1, sizeof(int));
  25
  26         for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
  27
  28         k = 0;
  29
  30         while (n > 3) {
  31
  32                 /*for(i=1;i<n;i++)
  33                   {
  34                     for(j=i+1;j<=n;j++)
  35                       {
  36                         Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
  37                       }
  38                     Rprintf("\n");
  39                   } */
  40
  41                 for (i = 1; i <= n; i++)
  42                 {double sum=0;
  43                   for( j=1;j<=n;j++)
  44                    {if(i==j)continue;
  45                      //Rprintf("index(%i,%i)=%i\n",i,j,give_index(i,j,n));
  46                      //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
  47                      sum+=D[give_index(i,j,n)];
  48                    }
  49                 S[i]=sum;
  50                 //Rprintf("\n");
  51                 //Rprintf("S[%i]=%f\n",i,S[i]);
  52                 //Rprintf("\n");
  53                 }
  54                 ij = 0;
  55                 smallest_S = 1e50;
  56                 B = n - 2;
  57                 for (i = 1; i < n; i++) {
  58                         for (j = i + 1; j <= n; j++) {
  59
  60                                 A = B*D[ij] - S[i] - S[j];
  61                                 /*Rprintf("D[ij]=%f\n",D[ij]);
  62                                 Rprintf("S[%i]=%f\n",i,S[i]);
  63                                 Rprintf("S[%i]=%f\n",j,S[j]);
  64                                 Rprintf("B=%f\n",B);
  65                                 Rprintf("A=%f\n",(B*D[ij] - S[i] - S[j]));
  66                                 Rprintf("Q[%i,%i]=%f\n",i,j,A);*/
  67                                 if (A < smallest_S) {
  68                                         OTU1 = i;
  69                                         OTU2 = j;
  70                                         smallest_S = A;
  71                                         smallest = ij;
  72                                 }
  73                                 ij++;
  74                         }
  75                 }
  76
  77                 /*Rprintf("agglomerating %i and %i, Q=%f \n",OTU1,OTU2,smallest_S);
  78
  79                 for(i=1;i<n;i++)
  80                   {
  81                     for(j=i+1;j<=n;j++)
  82                       {
  83                         Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
  84                       }
  85                     Rprintf("\n");
  86                   }
  87
  88                 for(i=1;i<n;i++)
  89                   {
  90                     for(j=i+1;j<=n;j++)
  91                       {
  92                         Rprintf("v[%i,%i]=%f ",i,j,v[give_index(i,j,n)]);
  93                       }
  94                     Rprintf("\n");
  95                   }*/
  96
  97                 edge2[k] = otu_label[OTU1];
  98                 edge2[k + 1] = otu_label[OTU2];
  99                 edge1[k] = edge1[k + 1] = cur_nod;
 100
 101                 /* get the distances between all OTUs but the 2 selected ones
 102                    and the latter:
 103                    a) get the sum for both
 104                    b) compute the distances for the new OTU */
 105                 double miu=0;
 106                 double miuSum=0;
 107                 for(i=1;i<=n;i++)
 108                  {
 109                    if(i == OTU1 || i==OTU2)continue;
 110                    //Rprintf("index(%i,%i)=%i index(%i,%i)=%i",i,OTU1,give_index(i,OTU1,n),i,OTU2,give_index(i,OTU2,n));
 111                    miuSum+=(1/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]));
 112                  }
 113                 miuSum=1/miuSum;
 114                 miu=miuSum/2;
 115
 116                 double eLenSum=0;
 117                 for(i=1;i<=n;i++)
 118                  {
 119                    if(i == OTU1 || i==OTU2)continue;
 120
 121                    double wi=miu/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]);
 122                    eLenSum+=wi*(D[give_index(i,OTU1,n)]-D[give_index(i,OTU2,n)]);
 123                  }
 124
 125                 edge_length[k]=D[give_index(OTU1,OTU2,n)]/2 + eLenSum;
 126
 127                 eLenSum=0;
 128                 for(i=1;i<=n;i++)
 129                  {
 130                    if(i == OTU1 || i==OTU2)continue;
 131
 132                    double wi=miu/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]);
 133                    eLenSum+=wi*(D[give_index(i,OTU2,n)]-D[give_index(i,OTU1,n)]);
 134                  }
 135
 136                 edge_length[k+1]=D[give_index(OTU1,OTU2,n)]/2 + eLenSum;
 137
 138                 A = D[smallest];
 139                 ij = 0;
 140                 for (i = 1; i <= n; i++) {
 141                         if (i == OTU1 || i == OTU2) continue;
 142                         double xi = D[give_index(i, OTU1, n)]; /* dist between OTU1 and i */
 143                         double yi = D[give_index(i, OTU2, n)]; /* dist between OTU2 and i */
 144                         double lamb=v[give_index(i,OTU2,n)]/(v[give_index(i,OTU2,n)]+v[give_index(i,OTU1,n)]);
 145                         new_dist[ij] = lamb*(xi-edge_length[k])+(1-lamb)*(yi-edge_length[k]);
 146                         new_v[ij]=(v[give_index(i,OTU2,n)]*v[give_index(i,OTU1,n)])/(v[give_index(i,OTU2,n)]+v[give_index(i,OTU1,n)]);
 147                         ij++;
 148                 }
 149                 /* compute the branch lengths */
 150                 //Rprintf("l2:%f \n",edge_length[k+1]);
 151                 /* update before the next loop
 152                    (we are sure that OTU1 < OTU2) */
 153                 if (OTU1 != 1)
 154                         for (i = OTU1; i > 1; i--)
 155                                 otu_label[i] = otu_label[i - 1];
 156                 if (OTU2 != n)
 157                         for (i = OTU2; i < n; i++)
 158                                 otu_label[i] = otu_label[i + 1];
 159                 otu_label[1] = cur_nod;
 160
 161                 for (i = 1; i < n; i++) {
 162                         if (i == OTU1 || i == OTU2) continue;
 163                         for (j = i + 1; j <= n; j++) {
 164                                 if (j == OTU1 || j == OTU2) continue;
 165                                 new_dist[ij] = D[DINDEX(i, j)];
 166                                 new_v[ij]=v[give_index(i,j,n)];
 167                                 ij++;
 168                         }
 169                 }
 170
 171                 n--;
 172                 for (i = 0; i < n*(n - 1)/2; i++)
 173                  {D[i] = new_dist[i];
 174                   v[i] = new_v[i];
 175                  }
 176                 cur_nod--;
 177                 k = k + 2;
 178         }
 179
 180         for (i = 0; i < 3; i++) {
 181                 edge1[*N*2 - 4 - i] = cur_nod;
 182                 edge2[*N*2 - 4 - i] = otu_label[i + 1];
 183         }
 184
 185         edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
 186         edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
 187         edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
 188 }
 189