- /* find the 'C' index of the smallest value of S */
- smallest = 0;
- for (i = 1; i < n*(n - 1)/2; i++)
- if (S[smallest] > S[i]) smallest = i;
-
- GET_I_AND_J
-
- SET_CLADE
-
- /* get the distances between all OTUs but the 2 selected ones
- and the latter:
- a) get the sum for both
- b) compute the distances for the new OTU */
- A = B = ij = 0;
- for (i = 1; i <= n; i++) {
- if (i == OTU1 || i == OTU2) continue;
- x = D[give_index(i, OTU1, n)]; /* distance between OTU1 and i */
- y = D[give_index(i, OTU2, n)]; /* distance between OTU2 and i */
- new_dist[ij] = (x + y)/2;
- A += x;
- B += y;
- ij++;
- }
- /* compute the branch lengths */
- A /= n - 2;
- B /= n - 2;
- edge_length[k] = (D[smallest] + A - B)/2;
- edge_length[k + 1] = (D[smallest] + B - A)/2;
- DI[cur_nod - *N - 1] = D[smallest];
-
- /* update before the next loop */
- if (OTU1 > OTU2) { /* make sure that OTU1 < OTU2 */
- i = OTU1;
- OTU1 = OTU2;
- OTU2 = i;
- }
- if (OTU1 != 1)
- for (i = OTU1 - 1; i > 0; i--) otu_label[i] = otu_label[i - 1];
- if (OTU2 != n)
- for (i = OTU2; i <= n; i++) otu_label[i - 1] = otu_label[i];
- otu_label[0] = cur_nod;
-
- for (i = 1; i < n; i++) {
- if (i == OTU1 || i == OTU2) continue;
- for (j = i + 1; j <= n; j++) {
- if (j == OTU1 || j == OTU2) continue;
- new_dist[ij] = D[DINDEX(i, j)];
- ij++;
- }
- }
-
- n--;
- for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
-
- cur_nod--;
- k = k + 2;
- }
-
- for (i = 0; i < 3; i++) {
- edge1[*N*2 - 4 - i] = cur_nod;
- edge2[*N*2 - 4 - i] = otu_label[i];
- }
-
- edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
- edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
- edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
-
- for (i = 0; i < *N*2 - 3; i++) {
- if (edge2[i] <= *N) continue;
- /* In case there are zero branch lengths (see above): */
- if (DI[edge2[i] - *N - 1] == 0) continue;
- edge_length[i] -= DI[edge2[i] - *N - 1]/2;
- }