+++ /dev/null
-/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
-
-/* This the fast UPGMA algorithm (O(N^2)) as implemented in Bob Edgar's
- * Muscle (UPGMA2.cpp; version 3.7) ported to pure C.
- *
- * Muscle's code is public domain and so is this code here.
- *
- * From http://www.drive5.com/muscle/license.htm:
- * """
- * MUSCLE is public domain software
- *
- * The MUSCLE software, including object and source code and
- * documentation, is hereby donated to the public domain.
- *
- * Disclaimer of warranty
- * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
- * EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION IMPLIED
- * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * """
- *
- */
-
-/*
- * RCS $Id: muscle_upgma.c 230 2011-04-09 15:37:50Z andreas $
- *
- *
- * Notes:
- * ------
- * LINKAGE become linkage_t here
- *
- * Replaced the the following member functions for DistCalc DC:
- * DC.GetId = sequence id as int
- * DC.GetName = sequence name
- * DC.GetCount = matrix dim
- * DC.DistRange = vector / matrix row for object i with index j<i
- *
- * Log() has been replaced with Clustal's Info(), Quiet() with Log(&rLog, LOG_FATAL)
- *
- * Made TriangleSubscript() and g_ulTriangleSize ulong to prevent overflow for many sequences
- */
-
-#ifndef ulint
-/* limit use of unsigned vars (see coding_style_guideline.txt) */
-typedef unsigned long int ulong;
-#endif
-
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <assert.h>
-
-#include "util.h"
-#include "log.h"
-#include "symmatrix.h"
-
-#include "muscle_tree.h"
-#include "muscle_upgma.h"
-
-/* from distcalc.h */
-typedef float dist_t;
-static const dist_t BIG_DIST = (dist_t) 1e29;
-/* from muscle.h */
-static const unsigned uInsane = 8888888;
-
-
-
-
-/*static inline*/
-ulong TriangleSubscript(uint uIndex1, uint uIndex2);
-
-
-
-
-#define TRACE 0
-
-#ifndef MIN
-#define MIN(x, y) ((x) < (y) ? (x) : (y))
-#endif
-#ifndef MIN
-#define MAX(x, y) ((x) > (y) ? (x) : (y))
-#endif
-#define AVG(x, y) (((x) + (y))/2)
-
-static uint g_uLeafCount;
-static ulong g_ulTriangleSize;
-static uint g_uInternalNodeCount;
-static uint g_uInternalNodeIndex;
-
-/* Triangular distance matrix is g_Dist, which is allocated
- * as a one-dimensional vector of length g_ulTriangleSize.
- * TriangleSubscript(i,j) maps row,column=i,j to the subscript
- * into this vector.
- * Row / column coordinates are a bit messy.
- * Initially they are leaf indexes 0..N-1.
- * But each time we create a new node (=new cluster, new subtree),
- * we re-use one of the two rows that become available (the children
- * of the new node). This saves memory.
- * We keep track of this through the g_uNodeIndex vector.
- */
-static dist_t *g_Dist;
-
-/* Distance to nearest neighbor in row i of distance matrix.
- * Subscript is distance matrix row.
- */
-static dist_t *g_MinDist;
-
-/* Nearest neighbor to row i of distance matrix.
- * Subscript is distance matrix row.
- */
-static uint *g_uNearestNeighbor;
-
-/* Node index of row i in distance matrix.
- * Node indexes are 0..N-1 for leaves, N..2N-2 for internal nodes.
- * Subscript is distance matrix row.
- */
-static uint *g_uNodeIndex;
-
-/* The following vectors are defined on internal nodes,
- * subscripts are internal node index 0..N-2.
- * For g_uLeft/Right, value is the node index 0 .. 2N-2
- * because a child can be internal or leaf.
- */
-static uint *g_uLeft;
-static uint *g_uRight;
-static dist_t *g_Height;
-static dist_t *g_LeftLength;
-static dist_t *g_RightLength;
-
-
-/*** CalcDistRange
- *
- * Imitation of DistCalc.DistRange
- *
- * Sets values of row (vector / matrix row) to distances for object i with index j<i
- *
- * row must be preallocated
- */
-void CalcDistRange(symmatrix_t *distmat, uint i, dist_t *row)
-{
- uint j;
- for (j = 0; j < i; ++j) {
- row[j] = SymMatrixGetValue(distmat, i, j);
- }
-}
-/* end of CalcDistRange */
-
-
-
-/*static inline*/
-ulong
-TriangleSubscript(uint uIndex1, uint uIndex2)
-{
- ulong v;
-#ifndef NDEBUG
- if (uIndex1 >= g_uLeafCount || uIndex2 >= g_uLeafCount)
- Log(&rLog, LOG_FATAL, "TriangleSubscript(%u,%u) %u", uIndex1, uIndex2, g_uLeafCount);
-#endif
- if (uIndex1 >= uIndex2)
- v = uIndex2 + (uIndex1*(uIndex1 - 1))/2;
- else
- v = uIndex1 + (uIndex2*(uIndex2 - 1))/2;
- assert(v < (g_uLeafCount*(g_uLeafCount - 1))/2);
- return v;
-}
-
-#ifdef UNUSED
-static void ListState()
-{
- uint i, j;
- Info("Dist matrix\n");
- Info(" ");
- for (i = 0; i < g_uLeafCount; ++i)
- {
- if (uInsane == g_uNodeIndex[i])
- continue;
- Info(" %5u", g_uNodeIndex[i]);
- }
- Info("\n");
-
- for (i = 0; i < g_uLeafCount; ++i)
- {
- if (uInsane == g_uNodeIndex[i])
- continue;
- Info("%5u ", g_uNodeIndex[i]);
- for (j = 0; j < g_uLeafCount; ++j)
- {
- if (uInsane == g_uNodeIndex[j])
- continue;
- if (i == j)
- Info(" ");
- else
- {
- ulong v = TriangleSubscript(i, j);
- Info("%5.2g ", g_Dist[v]);
- }
- }
- Info("\n");
- }
-
- Info("\n");
- Info(" i Node NrNb Dist\n");
- Info("----- ----- ----- --------\n");
- for (i = 0; i < g_uLeafCount; ++i)
- {
- if (uInsane == g_uNodeIndex[i])
- continue;
- Info("%5u %5u %5u %8.3f\n",
- i,
- g_uNodeIndex[i],
- g_uNearestNeighbor[i],
- g_MinDist[i]);
- }
-
- Info("\n");
- Info(" Node L R Height LLength RLength\n");
- Info("----- ----- ----- ------ ------- -------\n");
- for (i = 0; i <= g_uInternalNodeIndex; ++i)
- Info("%5u %5u %5u %6.2g %6.2g %6.2g\n",
- i,
- g_uLeft[i],
- g_uRight[i],
- g_Height[i],
- g_LeftLength[i],
- g_RightLength[i]);
-}
-#endif
-/* ifdef UNUSED */
-
-/**
- * @brief Creates a UPGMA in O(N^2) tree from given distmat
- *
- * @param[out] tree
- * newly created rooted UPGMA tree
- * @param[in] distmat
- * distance matrix to be clustered
- * @param[in] linkage
- * linkage type
- * @param[in] names
- * leaf names, will be copied
- *
- * @note called UPGMA2() in Muscle3.7.
- * caller has to free with FreeMuscleTree()
- *
- * @see FreeMuscleTree()
- */
-void
-MuscleUpgma2(tree_t *tree, symmatrix_t *distmat, linkage_t linkage, char **names)
-{
- int i, j;
- uint *Ids;
-
- /* only works on full symmetric matrices */
- assert (distmat->nrows==distmat->ncols);
-
- g_uLeafCount = distmat->ncols;
- g_ulTriangleSize = (g_uLeafCount*(g_uLeafCount - 1))/2;
- g_uInternalNodeCount = g_uLeafCount - 1;
-
- g_Dist = (dist_t *) CKMALLOC(g_ulTriangleSize * sizeof(dist_t));
-
- g_uNodeIndex = (uint*) CKMALLOC(sizeof(uint) * g_uLeafCount);
- g_uNearestNeighbor = (uint*) CKMALLOC(sizeof(uint) * g_uLeafCount);
- g_MinDist = (dist_t *) CKMALLOC(sizeof(dist_t) * g_uLeafCount);
- Ids = (uint*) CKMALLOC(sizeof(uint) * g_uLeafCount);
- /* NOTE: we replaced Names with argument names */
-
- /**
- * left and right node indices, as well as left and right
- * branch-length and height for for internal nodes
- */
- g_uLeft = (uint*) CKMALLOC(sizeof(uint) * g_uInternalNodeCount);
- g_uRight = (uint*) CKMALLOC(sizeof(uint) * g_uInternalNodeCount);
- g_Height = (dist_t*) CKMALLOC(sizeof(dist_t) * g_uInternalNodeCount);
- g_LeftLength = (dist_t*) CKMALLOC(sizeof(dist_t) * g_uInternalNodeCount);
- g_RightLength = (dist_t*) CKMALLOC(sizeof(dist_t) * g_uInternalNodeCount);
-
- for (i = 0; i < g_uLeafCount; ++i) {
- g_MinDist[i] = BIG_DIST;
- g_uNodeIndex[i] = i;
- g_uNearestNeighbor[i] = uInsane;
- Ids[i] = i;
- }
-
- for (i = 0; i < g_uInternalNodeCount; ++i) {
- g_uLeft[i] = uInsane;
- g_uRight[i] = uInsane;
- g_LeftLength[i] = BIG_DIST;
- g_RightLength[i] = BIG_DIST;
- g_Height[i] = BIG_DIST;
- }
-
-/* Compute initial NxN triangular distance matrix.
- * Store minimum distance for each full (not triangular) row.
- * Loop from 1, not 0, because "row" is 0, 1 ... i-1,
- * so nothing to do when i=0.
- */
- for (i = 1; i < g_uLeafCount; ++i) {
- dist_t *Row = g_Dist + TriangleSubscript(i, 0);
- CalcDistRange(distmat, i, Row);
- for (j = 0; j < i; ++j) {
- const dist_t d = Row[j];
- if (d < g_MinDist[i]) {
- g_MinDist[i] = d;
- g_uNearestNeighbor[i] = j;
- }
- if (d < g_MinDist[j]) {
- g_MinDist[j] = d;
- g_uNearestNeighbor[j] = i;
- }
- }
- }
-
-#if TRACE
- Info("Initial state:\n");
- ListState();
-#endif
-
- for (g_uInternalNodeIndex = 0;
- g_uInternalNodeIndex < g_uLeafCount - 1;
- ++g_uInternalNodeIndex) {
-
- dist_t dtNewMinDist = BIG_DIST;
- uint uNewNearestNeighbor = uInsane;
-
-#if TRACE
- Info("\n");
- Info("Internal node index %5u\n", g_uInternalNodeIndex);
- Info("-------------------------\n");
-#endif
-
- /* Find nearest neighbors */
- uint Lmin = uInsane;
- uint Rmin = uInsane;
- dist_t dtMinDist = BIG_DIST;
- for (j = 0; j < g_uLeafCount; ++j) {
- dist_t d;
- if (uInsane == g_uNodeIndex[j])
- continue;
-
- d = g_MinDist[j];
- if (d < dtMinDist) {
- dtMinDist = d;
- Lmin = j;
- Rmin = g_uNearestNeighbor[j];
- assert(uInsane != Rmin);
- assert(uInsane != g_uNodeIndex[Rmin]);
- }
- }
-
- assert(Lmin != uInsane);
- assert(Rmin != uInsane);
- assert(dtMinDist != BIG_DIST);
-
-#if TRACE
- Info("Nearest neighbors Lmin %u[=%u] Rmin %u[=%u] dist %.3g\n",
- Lmin,
- g_uNodeIndex[Lmin],
- Rmin,
- g_uNodeIndex[Rmin],
- dtMinDist);
-#endif
-
- /* Compute distances to new node
- * New node overwrites row currently assigned to Lmin
- */
- for ( j = 0; j < g_uLeafCount; ++j) {
- ulong vL, vR;
- dist_t dL, dR;
- dist_t dtNewDist;
-
- if (j == Lmin || j == Rmin)
- continue;
- if (uInsane == g_uNodeIndex[j])
- continue;
-
- vL = TriangleSubscript(Lmin, j);
- vR = TriangleSubscript(Rmin, j);
- dL = g_Dist[vL];
- dR = g_Dist[vR];
- dtNewDist = 0.0;
-
- switch (linkage) {
- case LINKAGE_AVG:
- dtNewDist = AVG(dL, dR);
- break;
-
- case LINKAGE_MIN:
- dtNewDist = MIN(dL, dR);
- break;
-
- case LINKAGE_MAX:
- dtNewDist = MAX(dL, dR);
- break;
-/* couldn't be arsed to figure out proper usage of g_dSUEFF */
-#if 0
- case LINKAGE_BIASED:
- dtNewDist = g_dSUEFF*AVG(dL, dR) + (1 - g_dSUEFF)*MIN(dL, dR);
- break;
-#endif
- default:
- Log(&rLog, LOG_FATAL, "UPGMA2: Invalid LINKAGE_%u", linkage);
- }
-
- /* Nasty special case.
- * If nearest neighbor of j is Lmin or Rmin, then make the new
- * node (which overwrites the row currently occupied by Lmin)
- * the nearest neighbor. This situation can occur when there are
- * equal distances in the matrix. If we don't make this fix,
- * the nearest neighbor pointer for j would become invalid.
- * (We don't need to test for == Lmin, because in that case
- * the net change needed is zero due to the change in row
- * numbering).
- */
- if (g_uNearestNeighbor[j] == Rmin)
- g_uNearestNeighbor[j] = Lmin;
-
-#if TRACE
- Info("New dist to %u = (%u/%.3g + %u/%.3g)/2 = %.3g\n",
- j, Lmin, dL, Rmin, dR, dtNewDist);
-#endif
- g_Dist[vL] = dtNewDist;
- if (dtNewDist < dtNewMinDist) {
- dtNewMinDist = dtNewDist;
- uNewNearestNeighbor = j;
- }
- }
-
- assert(g_uInternalNodeIndex < g_uLeafCount - 1 || BIG_DIST != dtNewMinDist);
- assert(g_uInternalNodeIndex < g_uLeafCount - 1 || uInsane != uNewNearestNeighbor);
-
- const ulong v = TriangleSubscript(Lmin, Rmin);
- const dist_t dLR = g_Dist[v];
- const dist_t dHeightNew = dLR/2;
- const uint uLeft = g_uNodeIndex[Lmin];
- const uint uRight = g_uNodeIndex[Rmin];
- const dist_t HeightLeft =
- uLeft < g_uLeafCount ? 0 : g_Height[uLeft - g_uLeafCount];
- const dist_t HeightRight =
- uRight < g_uLeafCount ? 0 : g_Height[uRight - g_uLeafCount];
-
- g_uLeft[g_uInternalNodeIndex] = uLeft;
- g_uRight[g_uInternalNodeIndex] = uRight;
- g_LeftLength[g_uInternalNodeIndex] = dHeightNew - HeightLeft;
- g_RightLength[g_uInternalNodeIndex] = dHeightNew - HeightRight;
- g_Height[g_uInternalNodeIndex] = dHeightNew;
-
- /* Row for left child overwritten by row for new node */
- g_uNodeIndex[Lmin] = g_uLeafCount + g_uInternalNodeIndex;
- g_uNearestNeighbor[Lmin] = uNewNearestNeighbor;
- g_MinDist[Lmin] = dtNewMinDist;
-
- /* Delete row for right child */
- g_uNodeIndex[Rmin] = uInsane;
-
-#if TRACE
- Info("\nInternalNodeIndex=%u Lmin=%u Rmin=%u\n",
- g_uInternalNodeIndex, Lmin, Rmin);
- ListState();
-#endif
- }
-
- uint uRoot = g_uLeafCount - 2;
-
-#if TRACE
- Log(&rLog, LOG_FORCED_DEBUG, "uRoot=%d g_uLeafCount=%d g_uInternalNodeCount=%d", uRoot, g_uLeafCount, g_uInternalNodeCount);
- for (i=0; i<g_uInternalNodeCount; i++) {
- Log(&rLog, LOG_FORCED_DEBUG, "internal node=%d: g_uLeft=%d g_uRight=%d g_LeftLength=%f g_RightLength=%f g_Height=%f",
- i, g_uLeft[i], g_uRight[i],
- g_LeftLength[i], g_RightLength[i],
- g_Height[i]);
- }
- for (i=0; i<g_uLeafCount; i++) {
- Log(&rLog, LOG_FORCED_DEBUG, "leaf node=%d: Ids=%d names=%s",
- i, Ids[i], names[i]);
- }
-#endif
-
- MuscleTreeCreate(tree, g_uLeafCount, uRoot,
- g_uLeft, g_uRight,
- g_LeftLength, g_RightLength,
- Ids, names);
-#if TRACE
- tree.LogMe();
-#endif
-
- free(g_Dist);
-
- free(g_uNodeIndex);
- free(g_uNearestNeighbor);
- free(g_MinDist);
- free(g_Height);
-
- free(g_uLeft);
- free(g_uRight);
- free(g_LeftLength);
- free(g_RightLength);
-
- /* NOTE: Muscle's "Names" variable is here the argument "names" */
- free(Ids);
-}
-/*** end of UPGMA2 ***/
git://source.jalview.org / jabaws.git / blobdiff