- POINTERS *trees; /* list of subtrees */
- tree *T; /* the returned tree */
- char *chain1; /* stringized branch-length */
-// char *chain2; /* idem */
- char *str; /* the string containing final tree */
- int *a, *b; /* pair to be agglomerated */
- float **delta; /* delta matrix */
- float la; /* first taxon branch-length */
- float lb; /* second taxon branch-length */
- float vab; /* variance of Dab */
- float lamda = 0.5;
- int i;
-// int ok;
- int r; /* number of subtrees */
- int n; /* number of taxa */
- int x, y;
-// double t;
- a=(int*)calloc(1,sizeof(int));
- b=(int*)calloc(1,sizeof(int));
- chain1=(char *)calloc(MAX_LABEL_LENGTH,sizeof(char));
- str = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
- /* added by EP */
- if (strlen(str))
- strncpy(str, "", strlen(str));
- /* end */
-
-// fscanf(input,"%d",&n);
- n = *N;
- /* Create the delta matrix */
- delta=(float **)calloc(n+1,sizeof(float*));
- for(i=1; i<= n; i++)
- {
- delta[i]=(float *)calloc(n+1, sizeof(float));
- if(delta[i] == NULL)
- {
- Rprintf("Out of memories!!");
- exit(0);
- }
- }
- trees=(POINTERS *)calloc(n+1,sizeof(POINTERS));
- if(trees == NULL)
- {
- Rprintf("Out of memories!!");
- exit(0);
- }
- /* initialise and symmetrize the running delta matrix */
- r=n;
- *a=0;
- *b=0;
- Initialize(delta, X, labels, n, trees);
-// ok=Symmetrize(delta, n);
-
-// if(!ok)
-// Rprintf("\n The matrix is not symmetric.\n ");
- while (r > 3) /* until r=3 */
- {
- Compute_sums_Sx(delta, n); /* compute the sum Sx */
- Best_pair(delta, r, a, b, n); /* find the best pair by */
- vab=Variance(*a, *b, delta); /* minimizing (1) */
- la=Branch_length(*a, *b, delta, r); /* compute branch-lengths */
- lb=Branch_length(*b, *a, delta, r); /* using formula (2) */
-// if (!isNJ)
- lamda=Lamda(*a, *b, vab, delta, n, r); /* compute lambda* using (9)*/
- for(i=1; i <= n; i++)
- {
- if(!Emptied(i,delta) && (i != *a) && (i != *b))
- {
- if(*a > i)
- {
- x=*a;
- y=i;
- }
- else
- {
- x=i;
- y=*a; /* apply reduction formulae */
- } /* 4 and 10 to delta */
- delta[x][y]=Reduction4(*a, la, *b, lb, i, lamda, delta);
- delta[y][x]=Reduction10(*a, *b, i, lamda, vab, delta);
- }
- }
- strncpy(chain1,"",1); /* agglomerate the subtrees */
- strncat(chain1,"(",1); /* a and b together with the*/
- Concatenate(chain1, *a, trees, 0); /* branch-lengths according */
- strncpy(chain1,"",1); /* to the NEWWICK format */
- strncat(chain1,":",1);
- snprintf(chain1+strlen(chain1),MAX_LABEL_LENGTH,"%f",la);
-
- strncat(chain1,",",1);
- Concatenate(chain1,*a, trees, 1);
- trees[*a].tail->suiv=trees[*b].head;
- trees[*a].tail=trees[*b].tail;
- strncpy(chain1,"",1);
- strncat(chain1,":",1);
- snprintf(chain1+strlen(chain1),MAX_LABEL_LENGTH,"%f",lb);
-
- strncat(chain1,")",1);
- Concatenate(chain1, *a, trees, 1);
- delta[*b][0]=1.0; /* make the b line empty */
- trees[*b].head=NULL;
- trees[*b].tail=NULL;
- r=r-1;
- }
-
- FinishStr (delta, n, trees, str); /* compute the branch-lengths*/
- /* of the last three subtrees*/
- /* and print the tree in the */
- /* output-file */
- T = readNewickString (str, n);
- T = detrifurcate(T);
-// compareSets(T,species);
- partitionSizes(T);
- *treeStr = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
- /* added by EP */
- if (strlen(*treeStr))
- strncpy(*treeStr, "", strlen(*treeStr));
- /* end */
- NewickPrintTreeStr (T, *treeStr);
-
- for(i=1; i<=n; i++)
- {
- delta[i][0]=0.0;
- trees[i].head=NULL;
- trees[i].tail=NULL;
- }
- free(delta);
- free(trees);
- /* free (str); */
- free (chain1);
- free (a);
- free (b);
- freeTree(T);
- T = NULL;
-// return (ret);
- return;