X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Fmvrs.c;h=fcfdc2520ff140c02e5ad5d59c2109678f587f4d;hb=dff741171e7afe3f9aaa2d9cb19c2f91995e8623;hp=d2aa342135a68536e407bfdeb62b3446a6f849e4;hpb=ef34642fe85694df0a80ef534137b8f6b427b67e;p=ape.git diff --git a/src/mvrs.c b/src/mvrs.c index d2aa342..fcfdc25 100644 --- a/src/mvrs.c +++ b/src/mvrs.c @@ -1,6 +1,6 @@ -/* mvrs.c 2011-10-11 */ +/* mvrs.c 2012-04-02 */ -/* Copyright 2011 Andrei-Alin Popescu */ +/* Copyright 2011-2012 Andrei-Alin Popescu */ /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -60,6 +60,8 @@ void mvrs(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_leng if(k==i || k==j) { + if(i!=k)R[give_index(i,j,n)]+=D[give_index(i,k,n)]; + if(j!=k)R[give_index(i,j,n)]+=D[give_index(j,k,n)]; s[give_index(i,j,n)]++; //Rprintf("%i",s[give_index(i,j,n)]); @@ -129,7 +131,7 @@ void mvrs(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_leng for (j = i + 1; j <= n; j++) { //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); A=S[i]+S[j]-B*D[give_index(i,j,n)]; - Rprintf("Q[%i,%i]=%f\n",i,j,A); + //Rprintf("Q[%i,%i]=%f\n",i,j,A); if (A > smallest_S) { OTU1 = i; OTU2 = j; @@ -288,6 +290,8 @@ void mvrs(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_leng { if(j==1 || j==i) { + if(i!=j)newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)]; + if(1!=j)newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)]; newS[give_index(1,i,n-1)]++; continue; } @@ -318,6 +322,7 @@ void mvrs(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_leng {if(new_dist[give_index(1,i,n-1)]==-1)continue; for(j=i+1;j<=n-1;j++) {if(new_dist[give_index(1,j,n-1)]==-1)continue; + if(new_dist[give_index(i,j,n-1)]==-1)continue; newR[give_index(i,j,n-1)]+=(new_dist[give_index(1,i,n-1)]+new_dist[give_index(1,j,n-1)]); newS[give_index(i,j,n-1)]++; } @@ -392,6 +397,3 @@ void mvrs(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_leng edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2; edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2; } - - -