--- /dev/null
+/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
+
+/* Module for deriving sequence weights from a tree. Largely based on
+ * Bob Edgar's Muscle (mainly clwwt.cpp; version 3.7). Ported to pure
+ * C. Most functions where apparently based on Clustal 1.8 anyway.
+ *
+ * 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: weights.c 231 2011-04-09 17:13:06Z andreas $
+ */
+
+/*
+ * Documentation from Muscle
+ *
+ * """
+ * Compute weights by the CLUSTALW method.
+ * Thompson, Higgins and Gibson (1994), CABIOS (10) 19-29;
+ * see also CLUSTALW paper.
+ *
+ * Weights are computed from the edge lengths of a rooted tree.
+ *
+ * Define the strength of an edge to be its length divided by the number
+ * of leaves under that edge. The weight of a sequence is then the sum
+ * of edge strengths on the path from the root to the leaf.
+ *
+ * Example.
+ *
+ * 0.2
+ * -----A 0.1
+ * -x ------- B 0.7
+ * --------y ----------- C
+ * 0.3 ----------z
+ * 0.4 -------------- D
+ * 0.8
+ *
+ * Edge Length Leaves Strength
+ * ---- ----- ------ --------
+ * xy 0.3 3 0.1
+ * xA 0.2 1 0.2
+ * yz 0.4 2 0.2
+ * yB 0.1 1 0.1
+ * zC 0.7 1 0.7
+ * zD 0.8 1 0.8
+ *
+ * Leaf Path Strengths Weight
+ * ---- ---- --------- ------
+ * A xA 0.2 0.2
+ * B xy-yB 0.1 + 0.1 0.2
+ * C xy-yz-zC 0.1 + 0.2 + 0.7 1.0
+ * D xy-yz-zD 0.1 + 0.2 + 0.8 1.1
+ * """
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <assert.h>
+#include "log.h"
+#include "muscle_tree.h"
+#include "weights.h"
+
+
+/*
+ #undef DEBUG
+*/
+
+
+/**
+ * @brief FIXME
+ *
+ * @param[out] puLeavesUnderNode
+ * FIXME
+ * @param[in] prTree
+ * FIXME
+ * @param[in] uNodeIndex
+ * FIXME
+ *
+ * @return The return value
+ *
+ * @note see Muscle3.7:clwwt.cpp
+ *
+ */
+uint
+CountLeaves(uint *puLeavesUnderNode, tree_t *prTree, uint uNodeIndex)
+{
+ uint uLeft;
+ uint uRight;
+ uint uRightCount;
+ uint uLeftCount;
+ uint uCount;
+
+
+ if (IsLeaf(uNodeIndex, prTree)) {
+ puLeavesUnderNode[uNodeIndex] = 1;
+ return 1;
+ }
+
+ uLeft = GetLeft(uNodeIndex, prTree);
+ uRight = GetRight(uNodeIndex, prTree);
+ uRightCount = CountLeaves(puLeavesUnderNode, prTree, uRight);
+ uLeftCount = CountLeaves(puLeavesUnderNode, prTree, uLeft);
+ uCount = uRightCount + uLeftCount;
+
+ puLeavesUnderNode[uNodeIndex] = uCount;
+
+ return uCount;
+}
+/*** end: CountLeaves() ***/
+
+
+
+
+/**
+ * @brief Normalise values in a double array to values between 0 and 1.
+ *
+ * @param[out] p
+ * double array with n elements
+ * @param[in] n
+ * number of elements in p
+ *
+ * @note From Muscle3.7: intmath.cpp:Normalize()
+ *
+ */
+void
+Normalise(double *p, uint n) {
+ unsigned i;
+ double dSum = 0.0;
+ for (i = 0; i < n; ++i) {
+ dSum += p[i];
+ }
+ if (0.0 == dSum) {
+ Log(&rLog, LOG_FATAL, "Normalise, sum=0");
+ }
+ for (i = 0; i < n; ++i) {
+ p[i] /= dSum;
+ }
+}
+/*** end: Normalise() ***/
+
+
+
+/**
+ * @brief Calculate "Clustal" weights from a tree.
+ *
+ * FIXME see doc in muscle:clwwt.cpp
+ *
+ * @param[out] pdWeights_p
+ * Will contain a weight for each leaf/sequence. Allocated here. User
+ * has to free
+ * @param[in] prTree
+ * Tree to derive weights from
+ *
+ * @return 0 on success, non-zero otherwise
+ *
+ * @note Largely copied from Muscle3.7: clwwt.cpp:CalcClustalWWeights()
+ *
+ * @warning FIXME Not sure if Muscle routines are most efficient here.
+ * Couldn't we do all this while traversing the tree and thereby safe
+ * time?
+ *
+ */
+int
+CalcClustalWeights(double **pdWeights_p, tree_t *prTree)
+{
+ int i; /* aux */
+ uint uLeafCount;
+ uint uNodeCount;
+ uint *puLeavesUnderNode;
+ uint uLeavesUnderRoot;
+ uint uRootNodeIndex;
+ double *pdStrengths;
+ uint uNodeIndex;
+ bool bLogWeights = FALSE; /* verbose output of weights */
+
+
+ assert(NULL != pdWeights_p);
+ assert(NULL != prTree);
+
+ if (rLog.iLogLevelEnabled <= LOG_DEBUG) {
+ bLogWeights = TRUE;
+ }
+
+ uLeafCount = GetLeafCount(prTree);
+ uNodeCount = GetNodeCount(prTree);
+
+
+ (*pdWeights_p) = (double *) CKMALLOC(uNodeCount * sizeof(double));
+
+ if (0 == uLeafCount) {
+ return 0;
+ } else if (1 == uLeafCount) {
+ (*pdWeights_p)[0] = 1.0;
+ return 0;
+ } else if (2 == uLeafCount) {
+ (*pdWeights_p)[0] = 0.5;
+ (*pdWeights_p)[1] = 0.5;
+ return 0;
+ }
+
+ if (!IsRooted(prTree)) {
+ Log(&rLog, LOG_ERROR, "Tree must be rooted to get weights");
+ CKFREE(pdWeights_p);
+ return -1;
+ }
+
+
+#ifdef TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "Weights follow");
+ fprintf(stderr, "Node Leaves Length Strength\n");
+ fprintf(stderr, "---- ------ -------- --------\n");
+ /* 1234 123456 12345678 12345678 */
+#endif
+
+ uRootNodeIndex = GetRootNodeIndex(prTree);
+ puLeavesUnderNode = (uint *) CKCALLOC(uNodeCount, sizeof(uint));
+
+ uLeavesUnderRoot = CountLeaves(puLeavesUnderNode, prTree, uRootNodeIndex);
+ if (uLeavesUnderRoot != uLeafCount) {
+ Log(&rLog, LOG_FATAL, "Internal error, root count %u %u",
+ uLeavesUnderRoot, uLeafCount);
+ }
+#if 0
+ for (uNodeIndex=0; uNodeIndex<uNodeCount; uNodeIndex++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "LeavesUnderNode[%d]=%d", uNodeIndex, puLeavesUnderNode[uNodeIndex]);
+ }
+#endif
+
+ pdStrengths = (double *) CKMALLOC(uNodeCount * sizeof(double));
+
+ for (uNodeIndex=0; uNodeIndex < uNodeCount; uNodeIndex++) {
+ uint uParent;
+ double dLength;
+ uint uLeaves;
+ double dStrength;
+
+ if (IsRoot(uNodeIndex, prTree)) {
+ pdStrengths[uNodeIndex] = 0.0;
+ continue;
+ }
+
+ uParent = GetParent(uNodeIndex, prTree);
+ dLength = GetEdgeLength(uNodeIndex, uParent, prTree);
+ uLeaves = puLeavesUnderNode[uNodeIndex];
+ dStrength = dLength / (double) uLeaves;
+ pdStrengths[uNodeIndex] = dStrength;
+
+#ifdef TRACE
+ fprintf(stderr, "%4u %6u %8g %8g\n", uNodeIndex, uLeaves, dLength, dStrength);
+#endif
+ }
+
+
+
+
+
+ if (bLogWeights){
+ fprintf(stderr, "\n");
+ fprintf(stderr, " Seq Path..Weight\n");
+ fprintf(stderr, "-------------------- ------------\n");
+ }
+ for (i=0; i<uLeafCount; i++) {
+ double dWeight = 0.0;
+ unsigned uLeafNodeIndex;
+ unsigned uNode;
+
+ uLeafNodeIndex = LeafIndexToNodeIndex(i, prTree);
+ uNode = uLeafNodeIndex;
+
+ if (bLogWeights){
+ fprintf(stderr, "%20.20s %4u ", GetLeafName(uLeafNodeIndex, prTree), uLeafNodeIndex);
+ }
+ if (! IsLeaf(uLeafNodeIndex, prTree)) {
+ Log(&rLog, LOG_FATAL,
+ "Internal error: non-leaf-node %d", uLeafNodeIndex);
+ }
+
+ /*LOG_DEBUG("dWeight = %f", dWeight);*/
+ while (! IsRoot(uNode, prTree)) {
+ dWeight += pdStrengths[uNode];
+ /*LOG_DEBUG("dWeight +== %f", pdStrengths[uNode]);*/
+ uNode = GetParent(uNode, prTree);
+ if (bLogWeights){
+ fprintf(stderr, "->%u(%g)", uNode, pdStrengths[uNode]);
+ }
+ }
+ /* AW: no idea what this is, but it's done like this in Muscle */
+ if (dWeight < 0.0001) {
+#ifdef TRACE
+ fprintf(stderr, "zero->one");
+#endif
+ dWeight = 1.0;
+ }
+
+ /* @note: the only difference to the muscle code is here: we
+ * use the input index for storing weights, instead of the
+ * tree leaf index
+ */
+ (*pdWeights_p)[GetLeafId(uLeafNodeIndex, prTree)] = dWeight;
+ if (bLogWeights){
+ fprintf(stderr, " = %g\n", dWeight);
+ }
+ }
+
+#if 0
+ for (i=0; i<uLeafCount; i++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "Weights before normalisation: pdWeights_p[%d]=%f", i, (*pdWeights_p)[i]);
+ /*LOG_DEBUG("Should be %d", GetLeafId(LeafIndexToNodeIndex(i, prTree), prTree));*/
+ }
+#endif
+
+ Normalise((*pdWeights_p), uLeafCount);
+
+
+ CKFREE(puLeavesUnderNode);
+ CKFREE(pdStrengths);
+
+ return 0;
+}
+/*** end: CalcWeights() ***/
git://source.jalview.org / jabaws.git / blobdiff