--- /dev/null
+/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
+
+/*********************************************************************
+ * Clustal Omega - Multiple sequence alignment
+ *
+ * Copyright (C) 2010 University College Dublin
+ *
+ * Clustal-Omega is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * This file is part of Clustal-Omega.
+ *
+ ********************************************************************/
+
+/*
+ * RCS $Id: mbed.c 254 2011-06-21 13:07:50Z andreas $
+ *
+ *
+ * Reimplementation from scratch of mBed:
+ * Blackshields et al. (2010); PMID 20470396
+ *
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <assert.h>
+#include <float.h>
+
+#include "mbed.h"
+#include "pair_dist.h"
+#include "symmatrix.h"
+#include "ktuple_pair.h"
+#include "tree.h"
+#include "util.h"
+#include "progress.h"
+#include "queue.h"
+
+#include "log.h"
+#include "kmpp/KMeans.h"
+#include "mbed.h"
+
+#define TIMING 0
+#if TIMING
+#include "squid/stopwatch.h"
+#endif
+
+
+/* If FULL_WITHIN_CLUSTER_DISTANCES is not 0, distances within each
+ * bisecting kmeans subcluster are not estimated using the vectors,
+ * but calculated normally (using ktuple or kimura). Surprisingly this
+ * results in 3% loss on a Homfam p24-h2010-08-09 subset (100-5000
+ * seqs in test, at least 5 ref seqs; MAX_SEQ 100 vs 10000; NUM_SEEDS
+ * log2 instead of log2^2). And of course it slows things down.
+ */
+#define FULL_WITHIN_CLUSTER_DISTANCES 1
+
+#define COMPUTE_WITHIN_SUBCLUSTER_AVERAGE 0
+
+/* Cluster size limits. Maximum is a soft limit, which might be
+ * exceeded if a K-Means split was unsuccesful, where unsuccesful
+ * might also mean that the minimum required number seqs. was not
+ * reached */
+#if FULL_WITHIN_CLUSTER_DISTANCES
+static const int MAX_ALLOWED_SEQ_PER_PRECLUSTER = 100;
+static const int MIN_REQUIRED_SEQ_PER_PRECLUSTER = 1;
+#else
+static const int MAX_ALLOWED_SEQ_PER_PRECLUSTER = 10000;
+static const int MIN_REQUIRED_SEQ_PER_PRECLUSTER = 100;
+#endif
+
+/* How many restarts per bisecting kmeans split. 10 give 0.5% increase
+ * in quality on original HOMFAM over just 1. It also increases kmeans
+ * time by a factor of ~3, but overall time is insignificant
+ * compared to pairdist/progressive alignment part.
+ */
+static const int RESTARTS_PER_SPLIT = 10;
+
+
+/* Use standard kmeans (lloyds) or kmeans++. Both are almost
+ * indistinguishable here, although kmeans++ is slightly ahead the
+ * fewer seeds you pick (and it should be better in theory)
+ */
+#define USE_KMEANS_LLOYDS 0
+
+
+#ifndef HAVE_LOG2
+#define log2(x) (log(x) / 0.69314718055994530942)
+#endif
+#define NUMBER_OF_SEEDS(n) pow(log2(((double)n)), 2)
+
+
+/* Seed selection method: SAMPLE_SEEDS_BY_LENGTH is the original mBed
+ * approach: Sample iSeeds sequence with constant stride from length-sorted X.
+ * It might be better to pick the seeds randomly, because length sorting might
+ * be more prone to including outliers (e.g. very long and very short seqs).
+ * However, we could never observer such a thing. So just stick to the
+ * original version as this also removes the random element.
+ */
+enum SEED_SELECTION_TYPE {
+ SELECT_SEEDS_RANDOMLY,
+ SELECT_SEEDS_BY_LENGTH
+};
+#define SEED_SELECTION SELECT_SEEDS_BY_LENGTH
+
+
+/* Tests on BAliBase (RV11,12,20,30,40,50; 10 runs each) show there is
+ * no difference between mbed-trees created from cosine or euclidean
+ * distances (simple version, just using disparities).
+ */
+#define USE_EUCLIDEAN_DISTANCE 1
+
+
+/* print some mbed pre-cluster usage to screen */
+#define PRINT_CLUSTER_DISTRIBUTION 0
+
+
+#define TRACE 0
+
+
+typedef struct {
+ /* Number of final clusters
+ */
+ int iNClusters;
+
+ /* Coordinates (columns) for each cluster (rows)
+ * valid indices: [0...iNClusters][0...dim-1]
+ */
+ double **ppdClusterCenters;
+
+ /* Dimensionality of input data and cluster center coordinates
+ */
+ int iDim;
+
+ /* Number of objects per cluster
+ */
+ int *piNObjsPerCluster;
+
+ /* Objects (indices) for each cluster. [i][j] will point to (index
+ * of) object j in cluster i
+ */
+ int **ppiObjIndicesPerCluster;
+} bisecting_kmeans_result_t;
+
+
+
+/**
+ * @brief Free KMeans result structure.
+ *
+ * @param[out] prResult_p
+ * K-Means result to free
+ *
+ * @see NewKMeansResult()
+ */
+void
+FreeKMeansResult(bisecting_kmeans_result_t **prResult_p)
+{
+ int iAux;
+ CKFREE((*prResult_p)->piNObjsPerCluster);
+ for (iAux=0; iAux<(*prResult_p)->iNClusters; iAux++) {
+ CKFREE((*prResult_p)->ppiObjIndicesPerCluster[iAux]);
+ CKFREE((*prResult_p)->ppdClusterCenters[iAux]);
+ }
+ CKFREE((*prResult_p)->ppiObjIndicesPerCluster);
+ CKFREE((*prResult_p)->ppdClusterCenters);
+ (*prResult_p)->iNClusters = 0;
+ (*prResult_p)->iDim = 0;
+ CKFREE(*prResult_p);
+}
+/*** end: FreeKMeansResult() ***/
+
+
+
+/**
+ * @brief Allocate new KMeans result
+ *
+ * @param[out] prKMeansResult_p
+ * K-Means result to free
+ *
+ * @see NewKMeansResult()
+ */
+void
+NewKMeansResult(bisecting_kmeans_result_t **prKMeansResult_p)
+{
+ (*prKMeansResult_p) = (bisecting_kmeans_result_t *)
+ CKCALLOC(1, sizeof(bisecting_kmeans_result_t));
+ (*prKMeansResult_p)->iNClusters = 0;
+ (*prKMeansResult_p)->iDim = 0;
+ (*prKMeansResult_p)->ppdClusterCenters = NULL;
+ (*prKMeansResult_p)->piNObjsPerCluster = NULL;
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster = NULL;
+
+}
+/*** end: NewKMeansResult() ***/
+
+
+
+/**
+ * @brief Calculate the euclidean distance between two vectors
+ *
+ * @param[in] v1
+ * First vector with dim dimensions
+ * @param[in] v2
+ * Second vector with dim dimensions
+ * @param[in] dim
+ * Dimensionality of v1 and v2
+ *
+ * @return euclidean distance as double
+ *
+ * @note Could probably be inlined. Also perfect for SSE
+ */
+double
+EuclDist(const double *v1, const double *v2, const int dim)
+{
+ int i; /* aux */
+ double dist; /* returned distance */
+
+ assert(v1!=NULL);
+ assert(v2!=NULL);
+
+ dist = 0.0;
+ for (i=0; i<dim; i++) {
+ dist += pow(v1[i]-v2[i], 2.0);
+ }
+
+ return sqrt(dist);
+}
+/*** end: EuclDist() ***/
+
+
+
+/**
+ * @brief Calculate the cosine distance between two vectors
+ *
+ * @param[in] v1
+ * First vector with dim dimensions
+ * @param[in] v2
+ * Second vector with dim dimensions
+ * @param[in] dim
+ * Dimensionality of v1 and v2
+ *
+ * @return cosine distance as double
+ *
+ * @note Could probably be inlined. Also perfect for SSE
+ */
+double
+CosDist(const double *v1, const double *v2, const int dim)
+{
+ int i; /* aux */
+ double dist; /* returned distance */
+ double s, sq1, sq2;
+
+ assert(v1!=NULL);
+ assert(v2!=NULL);
+
+ s = 0.0;
+ sq1 = 0.0;
+ sq2 = 0.0;
+ for (i=0; i<dim; i++) {
+ s += (v1[i]*v2[i]);
+ sq1 += pow(v1[i], 2.0);
+ sq2 += pow(v2[i], 2.0);
+ }
+ sq1 = sqrt(sq1);
+ sq2 = sqrt(sq2);
+
+ if ((sq1 * sq2) < DBL_EPSILON) {
+ dist = 1 - s / (sq1 * sq2);
+ } else {
+ /* if the denominator gets small, the fraction gets big, hence dist
+ gets small: */
+ dist = 0.0;
+ }
+ return dist;
+}
+/*** end: CosDist() ***/
+
+
+
+/**
+ * @brief convert sequences into mbed-like (distance) vector
+ * representation. Seeds (prMSeq sequence indices) have to be picked before
+ *
+ * @param[out] ppdSeqVec
+ * Pointer to preallocated matrix of size nseqs x iSeeds
+ * @param[in] prMSeq
+ * Sequences which are to be converted
+ * @param[in] piSeeds
+ * Array of sequences in prMSeq which are to be used as seeds.
+ * @param[in] iNumSeeds
+ * Number of seeds/elements in piSeeds
+ * @param[in] iPairDistType
+ * Type of pairwise distance comparison
+ *
+ */
+int
+SeqToVec(double **ppdSeqVec, mseq_t *prMSeq,
+ int *piSeeds, const int iNumSeeds,
+ const int iPairDistType)
+{
+ symmatrix_t *prDistmat = NULL;
+ int iSeqIndex;
+ int iSeedIndex;
+ /* indices for restoring order */
+ int *restore;
+ /* sorted copy of piSeeds */
+ int *piSortedSeeds;
+
+#if TIMING
+ Stopwatch_t *stopwatch = StopwatchCreate();
+ StopwatchZero(stopwatch);
+ StopwatchStart(stopwatch);
+#endif
+
+ assert(NULL != ppdSeqVec);
+ assert(iNumSeeds>0 && iNumSeeds<=prMSeq->nseqs);
+
+ piSortedSeeds = CKMALLOC(iNumSeeds * sizeof(int));
+ memcpy(piSortedSeeds, piSeeds, iNumSeeds*sizeof(int));
+
+ /* need them sorted, otherwise the swapping approach below might
+ * break
+ */
+ qsort(piSortedSeeds, iNumSeeds, sizeof(int), IntCmp);
+
+
+ /* rearrange mseq order so that we can use ktuple_pairdist code as
+ * is. That is, swap seeds and non-seeds such that the seeds end
+ * up in front of mseq. This way we can use the KTuplePairDist()
+ * code, without making a copy of mseq. Also, keep an array of
+ * indices which serves to restore the original sequence order
+ * after putting the seeds in front
+ *
+ */
+ restore = (int *) CKMALLOC(prMSeq->nseqs * sizeof(int));
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ restore[iSeqIndex] = iSeqIndex;
+ }
+ for (iSeedIndex=0; iSeedIndex<iNumSeeds; iSeedIndex++) {
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "Swapping seqs %d and %u", piSortedSeeds[iSeedIndex], iSeedIndex);
+#endif
+ if (piSortedSeeds[iSeedIndex]!=iSeedIndex) {
+ int swap;
+ SeqSwap(prMSeq, piSortedSeeds[iSeedIndex], iSeedIndex);
+
+ swap = restore[iSeedIndex];
+ restore[iSeedIndex] = restore[piSortedSeeds[iSeedIndex]];
+ restore[piSortedSeeds[iSeedIndex]] = swap;
+ }
+ }
+#if TRACE
+ printf("DEBUG(%s|%s():%d): restore indices =",
+ __FILE__, __FUNCTION__, __LINE__);
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ printf(" %u:%u", iSeqIndex, restore[iSeqIndex]);
+ }
+ printf("\n");
+
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "swapped seq no %d now: seq = %s",
+ iSeqIndex, prMSeq->sqinfo[iSeqIndex].name);
+ }
+#endif
+
+
+ /* convert sequences into vectors of distances by calling pairwise
+ * distance function for each seq/seed pair
+ *
+ */
+ if (0 != PairDistances(&prDistmat, prMSeq, iPairDistType,
+ 0, iNumSeeds, 0, prMSeq->nseqs,
+ NULL, NULL)) {
+ Log(&rLog, LOG_ERROR, "Could not compute pairwise distances for mbed.");
+ FreeSymMatrix(&prDistmat);
+ CKFREE(piSortedSeeds);
+ CKFREE(restore);
+ return -1;
+ }
+#if TRACE
+ {
+ char **labels;
+ labels = (char **) CKMALLOC(prMSeq->nseqs * sizeof(char *));
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ labels[iSeqIndex] = prMSeq->sqinfo[iSeqIndex].name;
+ }
+ SymMatrixPrint(prDistmat, labels, NULL);
+ CKFREE(labels);
+ }
+#endif
+
+
+ /* update ppdSeqVec according to prDistmat. keep in mind that we
+ * changed mseq order
+ *
+ */
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ for (iSeedIndex=0; iSeedIndex<iNumSeeds; iSeedIndex++) {
+ ppdSeqVec[restore[iSeqIndex]][iSeedIndex] =
+ SymMatrixGetValue(prDistmat, iSeqIndex, iSeedIndex);
+ }
+ }
+
+
+ /* restore mseq order by reverse swapping
+ *
+ * need reverse order now, so start from top index
+ */
+ iSeedIndex=iNumSeeds-1;
+ do {
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "Swapping seqs %d and %d", piSortedSeeds[iSeedIndex], iSeedIndex);
+#endif
+ if (piSortedSeeds[iSeedIndex]!=iSeedIndex) {
+ int swap;
+
+ SeqSwap(prMSeq, piSortedSeeds[iSeedIndex], iSeedIndex);
+
+ swap = restore[iSeedIndex];
+ restore[iSeedIndex] = restore[piSortedSeeds[iSeedIndex]];
+ restore[piSortedSeeds[iSeedIndex]] = swap;
+ }
+ } while (0 != iSeedIndex--);
+#if TRACE
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "restored seq no %d: seq = %s %s",
+ iSeqIndex, prMSeq->sqinfo[iSeqIndex].name, prMSeq->seq[iSeqIndex]);
+ }
+#endif
+
+#if TRACE
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ printf("DEBUG: seq %-4u as distance vector =", iSeqIndex);
+ for (iSeedIndex=0; iSeedIndex<iNumSeeds; iSeedIndex++) {
+ printf(" %f", ppdSeqVec[iSeqIndex][iSeedIndex]);
+ }
+ printf("\n");
+ }
+#endif
+
+
+ FreeSymMatrix(&prDistmat);
+ CKFREE(restore);
+#if TIMING
+ StopwatchStop(stopwatch);
+ StopwatchDisplay(stdout, "Total time for SeqToVec(): ", stopwatch);
+ StopwatchFree(stopwatch);
+#endif
+
+
+ return 0;
+}
+/*** end: SeqToVec() ***/
+
+
+
+/**
+ * @brief Select seeds to be used from an prMSeq
+ *
+ * @param[out] piSeeds
+ * Will store the indices of prMSeqs seqs used to be as seeds here. Must be preallocated.
+ * @param[in] iNumSeeds
+ * Number of seeds to be picked
+ * @param[in] iSelectionMethod
+ * Seed selection method to be used
+ * @param[in] prMSeq
+ * The prMSeq structure to pick sequences from
+ *
+ * @return: Non-zero on failure
+ *
+ */
+int
+SeedSelection(int *piSeeds, int iNumSeeds, int iSelectionMethod, mseq_t *prMSeq)
+{
+ /* seed iterator */
+ int iSeedIndex;
+ /* sequence iterator */
+ int iSeqIndex;
+
+ if (SELECT_SEEDS_RANDOMLY == iSelectionMethod) {
+ int *piPermArray;
+
+ Log(&rLog, LOG_INFO,
+ "Using %d seeds (randomly chosen) for mBed (from a total of %d sequences)",
+ iNumSeeds, prMSeq->nseqs);
+
+ PermutationArray(&piPermArray, iNumSeeds);
+ /* copy to piSeeds */
+ for (iSeedIndex=0; iSeedIndex<iNumSeeds; iSeedIndex++) {
+ piSeeds[iSeedIndex] = piPermArray[iSeedIndex];
+ }
+ CKFREE(piPermArray);
+
+ } else if (SELECT_SEEDS_BY_LENGTH == iSelectionMethod) {
+
+ /* step size for picking with constant stride */
+ int iStep;
+ int *piSeqLen = CKMALLOC(prMSeq->nseqs * sizeof(int));
+ int *piOrder = CKMALLOC(prMSeq->nseqs * sizeof(int));
+
+ Log(&rLog, LOG_INFO,
+ "Using %d seeds (chosen with constant stride from length sorted seqs) for mBed (from a total of %d sequences)",
+ iNumSeeds, prMSeq->nseqs);
+
+ iStep = prMSeq->nseqs / iNumSeeds; /* iStep will never get too big
+ * due to rounding */
+ /* first get an array of seq indices order according to
+ corresponding sequence length: piOrder */
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ piSeqLen[iSeqIndex] = prMSeq->sqinfo[iSeqIndex].len;
+ }
+ QSortAndTrackIndex(piOrder, piSeqLen, prMSeq->nseqs, 'd', FALSE);
+#if 0
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "Descending order (no %d): seq %d has len %d)",
+ iSeqIndex, piOrder[iSeqIndex], piSeqLen[piOrder[iSeqIndex]]);
+ }
+#endif
+ CKFREE(piSeqLen);
+ iSeedIndex = 0;
+ for (iSeqIndex=0; iSeedIndex<iNumSeeds; iSeqIndex+=iStep) {
+ piSeeds[iSeedIndex++] = piOrder[iSeqIndex];
+ }
+ CKFREE(piOrder);
+
+ } else {
+
+ Log(&rLog, LOG_ERROR, "Internal error: unknown seed selection type");
+ return -1;
+ }
+
+
+#ifdef PICK_SEEDS_FROM_KMEANS_CLUSTERS
+ /* do initial mbedding and kmeans. then pick seeds from each cluster and
+ do full mbed with these seeds. idea is that those seeds represent the
+ sequence space better than random seeds. however, this reduces the
+ quality slightly. random is almost always better.
+ */
+
+ if (1) {
+ /* seqs represented as distance vectors */
+ double **ppdSeqVec;
+ double dCost = -1.0;
+ /* assignments for each seq to the K newly created clusters */
+ int *piKMeansClusterAssignments = CKMALLOC(prMSeq->nseqs * sizeof(int));
+ double *pdKMeansClusterCenters = CKCALLOC(iNumSeeds * iNumSeeds, sizeof(double));
+ /* create a copy of ppdSeqVec suitable for KMeans */
+ double *pdKMeansVectors = CKMALLOC(prMSeq->nseqs * iNumSeeds * sizeof(double));
+ bool *pbClusterUsed = CKCALLOC(iNumSeeds, sizeof(bool));
+
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "FIXME Experimental feature: K-means on first round embedding to get better seeds");
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "FIXME Reuse seeds from clusters");
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "FIXME Could try log(n) instead of log(n)^2 in first round");
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "FIXME hardcoded iPairDistType PAIRDIST_KTUPLE");
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "FIXME Show that this is fast *and* good");
+
+ /*
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "Overriding iNumSeeds");
+ iNumSeeds = (int) pow(log2((double)prMSeq->nseqs), 2);
+ */
+
+ ppdSeqVec = (double **) CKMALLOC(prMSeq->nseqs * sizeof(double *));
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ ppdSeqVec[iSeqIndex] = (double *) CKMALLOC(iNumSeeds * sizeof(double));
+ }
+ if (0 != SeqToVec(ppdSeqVec, prMSeq, piSeeds, iNumSeeds, PAIRDIST_KTUPLE)) {
+ Log(&rLog, LOG_ERROR, "Could not convert sequences into vectors for mbed");
+ return -1;
+ }
+
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ int iDim;
+ for (iDim=0; iDim<iNumSeeds; ++iDim) {
+ pdKMeansVectors[iDim*iSeqIndex + iDim] = ppdSeqVec[iSeqIndex][iDim];
+ }
+ }
+
+
+ Log(&rLog, LOG_FORCED_DEBUG, "%s\n", "FIXME hardcoded RESTARTS_PER_SPLIT");
+ dCost = KMeans(prMSeq->nseqs, iNumSeeds, iNumSeeds,
+ pdKMeansVectors,
+ RESTARTS_PER_SPLIT, USE_KMEANS_LLOYDS,
+ pdKMeansClusterCenters, piKMeansClusterAssignments);
+ Log(&rLog, LOG_FORCED_DEBUG, "Best split cost = %f", dCost);
+
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "FIXME Check for Nan in cluster centers");
+
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "K-Means output:");
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ Log(&rLog, LOG_FORCED_DEBUG, " Raw assignment: Seq %u (%s) to cluster %d",
+ iSeqIndex, prMSeq->sqinfo[iSeqIndex].name,
+ piKMeansClusterAssignments[iSeqIndex]);
+ }
+#endif
+
+ Log(&rLog, LOG_FORCED_DEBUG, "FIXME %s", "proof of concept implementation: Pick first sequences from each clusters instead of reusing seeds");
+ iSeedIndex = 0;
+ for (iSeqIndex=0; iSeqIndex<prMSeq->nseqs; iSeqIndex++) {
+ int iAssignedCluster = piKMeansClusterAssignments[iSeqIndex];
+ if (pbClusterUsed[iAssignedCluster]) {
+ continue;
+ } else {
+ /*LOG_DEBUG("Picked seed %d from cluster %d", iSeqIndex,iAssignedCluster);*/
+ piSeeds[iSeedIndex++] = iSeqIndex;
+ pbClusterUsed[iAssignedCluster] = TRUE;
+ }
+ }
+ CKFREE(pbClusterUsed);
+
+ CKFREE(pdKMeansVectors);
+ CKFREE(pdKMeansClusterCenters);
+ CKFREE(piKMeansClusterAssignments);
+
+ }
+ /* if (1) */
+#endif
+
+
+ if (LOG_DEBUG <= rLog.iLogLevelEnabled) {
+ for (iSeedIndex=0; iSeedIndex<iNumSeeds; iSeedIndex++) {
+ Log(&rLog, LOG_DEBUG, "Picked sequence %d (%s) as seed no %d",
+ piSeeds[iSeedIndex], prMSeq->sqinfo[piSeeds[iSeedIndex]].name, iSeedIndex);
+ }
+ }
+
+ return 0;
+}
+/* end of SeedSelection() */
+
+
+
+
+/**
+ * @brief Bisecting K-Means clustering. Repeatedly calls K-Means with
+ * a K of 2 until no cluster has more than iMaxAllowedObjsPerCluster.
+ *
+ * @param[out] prKMeansResult_p
+ * Result of Bisecting KMeans. Will be allocated here.
+ * Caller has to free. See @see FreeKMeansResult()
+ * @param[in] iNObjs
+ * Number of objects/sequences to cluster
+ * @param[in] iDim
+ * Dimensionality of input data
+ * @param[in] ppdVectors
+ * each row holds iDim points for this object's coordinates
+ * @param[in] iMinRequiredObjsPerCluster
+ * Minimum number of objects per Cluster (inclusive)/
+ * @param[in] iMaxAllowedObjsPerCluster
+ * Maximum number of objects per Cluster (inclusive). Function returns once no
+ * cluster contains more then this number of objects. Soft limit!
+ *
+ * @note Convoluted code. Could use some restructuring. My apologies.
+ * AW
+ *
+ */
+void
+BisectingKmeans(bisecting_kmeans_result_t **prKMeansResult_p,
+ const int iNObjs, const int iDim, double **ppdVectors,
+ const int iMinRequiredObjsPerCluster,
+ const int iMaxAllowedObjsPerCluster)
+{
+ int iN, iD;
+ /* cluster centers for each cluster created at each split */
+ double *pdKClusterCenters;
+ /* keep track of updated object indices per newly created
+ * cluster */
+ int *piCurObjToUpdate;
+ /* number of times we called 2-means */
+ int iNRounds;
+ /* flag for unsuccessful cluster split */
+ bool bNaNDetected = FALSE;
+ /* flag for detected small cluster after split */
+ bool bSmallClusterDetected;
+ /* queue of clusters which are to be split */
+ int_queue_t rClusterSplitQueue;
+
+#if TIMING
+ Stopwatch_t *stopwatch = StopwatchCreate();
+#endif
+
+
+ piCurObjToUpdate = (int *) CKMALLOC(2 * sizeof(int));
+
+ NewKMeansResult(prKMeansResult_p);
+
+ /* new cluster centers created at each split/KMeans run
+ */
+ pdKClusterCenters = (double *) CKCALLOC(2 * iDim, sizeof(double));
+
+ /* init results by setting a first cluster that contains all objects
+ */
+ (*prKMeansResult_p)->iNClusters = 1;
+ (*prKMeansResult_p)->iDim = iDim;
+ /* fake center coordinates of first cluster */
+ (*prKMeansResult_p)->ppdClusterCenters =
+ (double **) CKMALLOC(1 * sizeof(double *));
+ (*prKMeansResult_p)->ppdClusterCenters[0] =
+ (double *) CKMALLOC(iDim * sizeof(double));
+ /* objects per cluster */
+ (*prKMeansResult_p)->piNObjsPerCluster =
+ (int *) CKMALLOC(1 * sizeof(int));
+ (*prKMeansResult_p)->piNObjsPerCluster[0] = iNObjs;
+ /* object indices per cluster */
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster =
+ (int **) CKMALLOC(1 * sizeof(int *));
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[0] = (int *)
+ CKMALLOC(iNObjs * sizeof(int));
+ for (iN=0; iN<iNObjs; iN++) {
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[0][iN] = iN;
+ }
+
+
+
+ /* Starting with the first cluster that now contains all the
+ * sequences keep splitting until no cluster contains more than
+ * iMaxAllowedObjsPerCluster
+ *
+ * Keep a queue of clusters (rClusterSplitQueue) to split
+ *
+ * At each split values/memory of the just split cluster will be
+ * reused and exactly one new only allocated.
+ *
+ * Example:
+ * Given the following cluster assignments
+ * 0 0 0 0 1 1 2 2 2 3 3 3 3 3 3 3 4 4
+ * and a K-Means split in cluster 3 at |:
+ * 0 0 0 0 1 1 2 2 2 3 3 3 | 3 3 3 3 4 4
+ * The result should be this:
+ * 0 0 0 0 1 1 2 2 2 3 3 3 | 5 5 5 5 4 4
+ *
+ *
+ */
+ INT_QUEUE_INIT(&rClusterSplitQueue);
+ if (iNObjs>iMaxAllowedObjsPerCluster) {
+ /* pish fake first cluster index */
+ INT_QUEUE_PUSH(&rClusterSplitQueue, 0);
+ }
+ iNRounds = 0;
+ while (! INT_QUEUE_EMPTY(&rClusterSplitQueue)) {
+ /* assignments for each seq to the K newly created clusters */
+ int *piKClusterAssignments;
+ /* number of objects in cluster that is to be split */
+ int iNObjsInClusterToSplit;
+ /* coordinates of points in cluster that is to be split
+ * array of size n*d where [d*i + j] gives coordinate j of point i
+ */
+ double *pdVectorsInClusterToSplit;
+ /* copy of object indices in split cluster */
+ int *piObjIndicesOfSplitCluster;
+ /* best cost of kmeans rounds */
+ double dCost = -1.0;
+
+ /* indices for the two created clusters
+ */
+ /* index of cluster to split */
+ int iClusterToSplot;
+ /* index of newly created cluster at each round. data for the
+ other created cluster goes to the one just split,
+ i.e. (*piClusterToSplot) */
+ int iNewClusterIdx;
+
+#if TIMING
+ StopwatchZero(stopwatch);
+ StopwatchStart(stopwatch);
+#endif
+
+ INT_QUEUE_POP(&rClusterSplitQueue, &iClusterToSplot);
+
+ iNObjsInClusterToSplit = (*prKMeansResult_p)->piNObjsPerCluster[iClusterToSplot];
+ piKClusterAssignments = (int *)
+ CKMALLOC(iNObjsInClusterToSplit * sizeof(int));
+
+ pdVectorsInClusterToSplit = (double *)
+ CKMALLOC(iNObjsInClusterToSplit * iDim * sizeof(double));
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ for (iD=0; iD<iDim; ++iD) {
+ int iThisObjIdx =
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iClusterToSplot][iN];
+ pdVectorsInClusterToSplit[iDim*iN + iD] = ppdVectors[iThisObjIdx][iD];
+ }
+ }
+
+#if TRACE
+ Log(&rLog, LOG_FOCRED_DEBUG, "Round %d: Will split cluster %d which has %d objects",
+ iNRounds, iClusterToSplot, iNObjsInClusterToSplit);
+ fprintf(stderr, "DEBUG(%s|%s():%d): Object indices in cluster to split are:",
+ __FILE__, __FUNCTION__, __LINE__);
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ fprintf(stderr, " %u",
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iClusterToSplot][iN]);
+ }
+ fprintf(stderr, "\n");
+ (void) fflush(stderr);
+#endif
+
+#if TRACE
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ fprintf(stderr,
+ "DEBUG(%s|%s():%d): Coordinate of object %u (real index %u) in cluster to split:",
+ __FILE__, __FUNCTION__, __LINE__, iN,
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iClusterToSplot][iN]);
+ for (iD=0; iD<iDim; iD++) {
+ fprintf(stderr, " %f", pdVectorsInClusterToSplit[iDim*iN + iD]);
+ }
+ fprintf(stderr, "\n");
+ (void) fflush(stderr);
+ }
+#endif
+
+ /* KMeans(1 "The number of points in the data set",
+ * 2 "The number of clusters to look for",
+ * 3 "The number of dimensions that the data set lives in",
+ * 4 "points: An array of size n*d where points[d*i + j] gives coordinate j of point i",
+ * 5 "attempts: The number of times to independently run k-means",
+ * 6 "use_lloyds_method: uses kmpp if false, otherwise lloyds method",
+ * 7 "centers: This can either be null or an array of size k*d.
+ * In the latter case, centers[d*i + j] will give coordinate j of center i.
+ * If the cluster is unused, it will contain NaN instead.",
+ * 8 "assignments: This can either be null or an array of size n.
+ * In the latter case, it will be filled with the cluster that each point is assigned to
+ * (an integer between 0 and k-1 inclusive).");
+ */
+ dCost = KMeans(iNObjsInClusterToSplit, 2, iDim,
+ pdVectorsInClusterToSplit,
+ RESTARTS_PER_SPLIT, USE_KMEANS_LLOYDS,
+ pdKClusterCenters, piKClusterAssignments);
+
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "Raw K-Means output:");
+ for (iN=0; iN<2; iN++) {
+ fprintf(stderr, "DEBUG(%s|%s():%d): Cluster Center %u =",
+ __FILE__, __FUNCTION__, __LINE__, iN);
+ for (iD=0; iD<iDim; iD++) {
+ fprintf(stderr, " %f", pdKClusterCenters[iN*iDim+iD]);
+ }
+ fprintf(stderr, "\n");
+ (void) fflush(stderr);
+ }
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ Log(&rLog, LOG_FORCED_DEBUG, " Raw assignment: Seq %u to cluster %d (of #%u)",
+ iN, piKClusterAssignments[iN], 2);
+ }
+#endif
+
+ /* real index of one of the newly created clusters. the other
+ * one is iClusterToSplot */
+ iNewClusterIdx = (*prKMeansResult_p)->iNClusters;
+
+
+ /* We don't want Clusters which are too small. Check here if a
+ * split created a small cluster and if yes, discard the
+ * solution. Because the cluster has already been removed from
+ * the queue, we can just continue.
+ */
+ bSmallClusterDetected = FALSE;
+ if (iMinRequiredObjsPerCluster>1) {
+ int iNObjsCluster[2];
+ iNObjsCluster[0] = 0; /* first cluster */
+ iNObjsCluster[1] = 0; /* second cluster */
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ iNObjsCluster[piKClusterAssignments[iN]]+=1;
+ }
+
+ if (iNObjsCluster[0]<iMinRequiredObjsPerCluster
+ ||
+ iNObjsCluster[1]<iMinRequiredObjsPerCluster) {
+ bSmallClusterDetected = TRUE;
+ Log(&rLog, LOG_FORCED_DEBUG, "Skipping this split because objs in 1st/2nd cluster = %d/%d < %d",
+ iNObjsCluster[0], iNObjsCluster[1], iMinRequiredObjsPerCluster);
+ }
+ }
+
+ /* Same logic as for small clusters applies if KMeans couldn't
+ * split the cluster. In this case one of its center
+ * coordinates will be NaN, which we check for in the
+ * following.
+ */
+ if (! bSmallClusterDetected) {
+ for (iN=0; iN<2; iN++) {
+ bNaNDetected = FALSE;
+ for (iD=0; iD<iDim; iD++) {
+ if (0 != isnan(pdKClusterCenters[iN*iDim+iN])) {
+ /* Got NaN as coordinate after splitting */
+ Log(&rLog, LOG_WARN, "%s(): Can't split cluster no. %d which has %d objects any further. %s",
+ __FUNCTION__,
+ iClusterToSplot, iNObjsInClusterToSplit,
+ "Hope it's not too big and doesn't slow things down."
+ );
+ bNaNDetected = TRUE;
+ break;
+ }
+ }
+ if (bNaNDetected) {
+ break;
+ }
+ }
+ }
+
+ /* Discarding split of this cluster. It has been removed from the
+ * queue already. so just continue
+ */
+ if (bNaNDetected || bSmallClusterDetected) {
+ CKFREE(piKClusterAssignments);
+ CKFREE(pdVectorsInClusterToSplit);
+ continue;
+ }
+
+
+ /* update cluster centers: pdClusterCenters
+ *
+ */
+ /* reuse memory of existing/old/split cluster
+ */
+ for (iN=0; iN<iDim; iN++) {
+ double dCoord = pdKClusterCenters[0*iDim+iN];
+ (*prKMeansResult_p)->ppdClusterCenters[iClusterToSplot][iN] = dCoord;
+ }
+ /* realloc and set new one
+ */
+ (*prKMeansResult_p)->ppdClusterCenters = (double **)
+ CKREALLOC((*prKMeansResult_p)->ppdClusterCenters,
+ ((*prKMeansResult_p)->iNClusters+1) * sizeof(double *));
+ (*prKMeansResult_p)->ppdClusterCenters[iNewClusterIdx] = (double *)
+ CKMALLOC(iDim * sizeof(double));
+ for (iD=0; iD<iDim; iD++) {
+ double dCoord = pdKClusterCenters[1*iDim+iD];
+ (*prKMeansResult_p)->ppdClusterCenters[iNewClusterIdx][iD] = dCoord;
+ }
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "* Update of cluster centers done. Cluster centers so far:");
+ for (iN=0; iN<(*prKMeansResult_p)->iNClusters+1; iN++) {
+ fprintf(stderr, "DEBUG(%s|%s():%d): Center %u =",
+ __FILE__, __FUNCTION__, __LINE__, iN);
+ for (iD=0; iD<iDim; iD++) {
+ fprintf(stderr, " %f", (*prKMeansResult_p)->ppdClusterCenters[iN][iD]);
+ }
+ fprintf(stderr, "\n");
+ (void) fflush(stderr);
+ }
+#endif
+
+
+ /* update #seqs per cluster: piNObjsPerCluster
+ *
+ */
+ (*prKMeansResult_p)->piNObjsPerCluster = (int *)
+ CKREALLOC((*prKMeansResult_p)->piNObjsPerCluster,
+ ((*prKMeansResult_p)->iNClusters+1) * sizeof(int));
+ /* init new and old one to zero */
+ (*prKMeansResult_p)->piNObjsPerCluster[iClusterToSplot] = 0;
+ (*prKMeansResult_p)->piNObjsPerCluster[iNewClusterIdx] = 0;
+ /* now update values */
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ if (0 == piKClusterAssignments[iN]) {
+ (*prKMeansResult_p)->piNObjsPerCluster[iClusterToSplot] += 1;
+ } else if (1 == piKClusterAssignments[iN]) {
+ (*prKMeansResult_p)->piNObjsPerCluster[iNewClusterIdx] += 1;
+ } else {
+ /* there used to be code for iK>=2 in r101 */
+ Log(&rLog, LOG_FATAL, "Internal error: split into more than two clusters (got assignment %d)",
+ piKClusterAssignments[iN]);
+ }
+ }
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "* Update of NObjs per cluster done:");
+ for (iN=0; iN<(*prKMeansResult_p)->iNClusters+1; iN++) {
+ Log(&rLog, LOG_FORCED_DEBUG, "Cluster %d contains %d sequences",
+ iN, (*prKMeansResult_p)->piNObjsPerCluster[iN]);
+ }
+#endif
+ /* queue clusters if still they are still too big
+ */
+ if ((*prKMeansResult_p)->piNObjsPerCluster[iClusterToSplot] > iMaxAllowedObjsPerCluster) {
+ INT_QUEUE_PUSH(&rClusterSplitQueue, iClusterToSplot);
+ }
+ if ((*prKMeansResult_p)->piNObjsPerCluster[iNewClusterIdx] > iMaxAllowedObjsPerCluster) {
+ INT_QUEUE_PUSH(&rClusterSplitQueue, iNewClusterIdx);
+ }
+
+
+
+ /* update which objects belong to which cluster:
+ *
+ * note: piNObjsPerCluster needs to be already updated
+ *
+ */
+ /* create a copy of the object indices in the split cluster
+ * (original will be overwritten)
+ */
+ piObjIndicesOfSplitCluster = (int *) CKMALLOC(iNObjsInClusterToSplit * sizeof(int));
+ memcpy(piObjIndicesOfSplitCluster,
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iClusterToSplot],
+ iNObjsInClusterToSplit * sizeof(int));
+
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster = (int **)
+ CKREALLOC((*prKMeansResult_p)->ppiObjIndicesPerCluster,
+ ((*prKMeansResult_p)->iNClusters+1) * sizeof(int *));
+
+
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iClusterToSplot] =
+ (int *) CKREALLOC((*prKMeansResult_p)->ppiObjIndicesPerCluster[iClusterToSplot],
+ (*prKMeansResult_p)->piNObjsPerCluster[iClusterToSplot] * sizeof(int));
+
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iNewClusterIdx] =
+ (int *) CKMALLOC((*prKMeansResult_p)->piNObjsPerCluster[iNewClusterIdx] * sizeof(int));
+
+
+ /* now reassign the object indices to the assigned cluster
+ */
+ piCurObjToUpdate[0] = 0;
+ piCurObjToUpdate[1] = 0;
+ for (iN=0; iN<iNObjsInClusterToSplit; iN++) {
+ int iThisObjIdx = piObjIndicesOfSplitCluster[iN];
+ int iThisClusterAssignment = piKClusterAssignments[iN];
+ int iThisOffset = piCurObjToUpdate[iThisClusterAssignment];
+ int iThisClusterIndex = 0;
+
+ if (0 == iThisClusterAssignment) {
+ iThisClusterIndex = iClusterToSplot;
+ } else if (1 == iThisClusterAssignment) {
+ iThisClusterIndex = iNewClusterIdx;
+ } else {
+ /* there used to be code for iK>=2 in r101 */
+ Log(&rLog, LOG_FATAL, "Internal error: split into more than two clusters (got assignment %d)",
+ piKClusterAssignments[iN]);
+ }
+#if 0
+ Log(&rLog, LOG_FORCED_DEBUG, "Setting (*prKMeansResult_p)->ppiObjIndicesPerCluster[%d][%d] = %d",
+ iThisClusterIndex, iThisOffset, iThisObjIdx);
+#endif
+ (*prKMeansResult_p)->ppiObjIndicesPerCluster[iThisClusterIndex][iThisOffset] = iThisObjIdx;
+ piCurObjToUpdate[iThisClusterAssignment]+=1;
+ }
+ CKFREE(piObjIndicesOfSplitCluster);
+#if TRACE
+ for (iN=0; iN<(*prKMeansResult_p)->iNClusters+1; iN++) {
+ int iObj;
+ fprintf(stderr, "DEBUG(%s|%s():%d): Objects in cluster %u: ",
+ __FILE__, __FUNCTION__, __LINE__, iN);
+ for (iObj=0; iObj<(*prKMeansResult_p)->piNObjsPerCluster[iN]; iObj++) {
+ fprintf(stderr, " %u", (*prKMeansResult_p)->ppiObjIndicesPerCluster[iN][iObj]);
+ }
+ fprintf(stderr, "\n");
+ (void) fflush(stderr);
+ }
+#endif
+
+
+#if TIMING
+ StopwatchStop(stopwatch);
+ StopwatchDisplay(stdout, "Total time after next round in Bisecting-KMeans: ", stopwatch);
+#endif
+
+
+ /* finally: increase number of clusters
+ */
+ (*prKMeansResult_p)->iNClusters += 1;
+ iNRounds += 1;
+ CKFREE(piKClusterAssignments);
+ CKFREE(pdVectorsInClusterToSplit);
+
+ } /* while */
+ INT_QUEUE_DESTROY(&rClusterSplitQueue);
+
+
+ Log(&rLog, LOG_DEBUG,
+ "Bisecting K-means finished after %d rounds (no more clusters to split)",
+ iNRounds);
+
+#if TIMING
+ StopwatchFree(stopwatch);
+#endif
+
+ /* @note could use progress/timer */
+
+ CKFREE(pdKClusterCenters);
+ CKFREE(piCurObjToUpdate);
+
+ return;
+}
+/*** end: BisectingKmeans() ***/
+
+
+
+/**
+ *
+ * @brief From scratch reimplementation of mBed: Blackshields et al.
+ * (2010); PMID 20470396.
+ *
+ * Idea is a follows:
+ * - convert sequences into vectors of distances
+ * - cluster the vectors using k-means
+ * - cluster each of the k clusters using upgma (used cached distances
+ * from above?)
+ * - join the sub-clusters to create on tree (use UPGMA on k-means
+ * medoids)
+ *
+ *
+ * @param[out] prMbedTree_p
+ * Created upgma tree. will be allocated here. use FreeMuscleTree()
+ * to free
+ * @param[in] prMSeq
+ * Multiple sequences
+ * @param[in] iPairDistType
+ * Distance measure for pairwise alignments
+ * @param[in] pcGuidetreeOut
+ * Passed down to GuideTreeUpgma()
+ *
+ * @return Zero on success, non-zero on error
+ *
+ */
+int
+Mbed(tree_t **prMbedTree_p, mseq_t *prMSeq, const int iPairDistType,
+ const char *pcGuidetreeOut)
+{
+ /* number of seeds */
+ int iNumSeeds;
+ /* seed indices matching prMSeq */
+ int *piSeeds;
+ /* seqs represented as distance vectors */
+ double **ppdSeqVec;
+ /* kmeans result */
+ bisecting_kmeans_result_t *prKMeansResult = NULL;
+ /* distance matrix of kmeans (pre-)cluster centers */
+ symmatrix_t *prPreClusterDistmat = NULL;
+ /* auxiliary for symmetric matrix output; tree routines etc */
+ char **ppcLabels;
+ int iNodeIndex;
+ /* mapping of cluster-center tree node indices to corresponding
+ cluster */
+ int *piClusterToTreeNode;
+ int iClusterIndex;
+ int iI, iJ;
+ FILE *pfOut;
+ progress_t *prSubClusterDistanceProgress;
+ bool bPrintCR = (rLog.iLogLevelEnabled<=LOG_VERBOSE) ? FALSE : TRUE;
+
+#if FULL_WITHIN_CLUSTER_DISTANCES
+ Log(&rLog, LOG_DEBUG, "Computing real distances within subclusters for mBed.");
+#else
+ Log(&rLog, LOG_DEBUG, "Computing vector distances within subclusters for mBed.");
+#endif
+
+#if MBED_TIMING
+ Stopwatch_t *stopwatch = StopwatchCreate();
+#endif
+
+ assert(NULL != prMbedTree_p);
+ assert(NULL != prMSeq);
+
+
+ iNumSeeds = (int) NUMBER_OF_SEEDS(prMSeq->nseqs);
+ if (iNumSeeds >= prMSeq->nseqs) {
+ /* -1 is condition for RandomUniqueIntArray */
+ iNumSeeds = prMSeq->nseqs-1;
+ Log(&rLog, LOG_DEBUG,
+ "Automatically determined number of seeds is bigger (or equal) the number of sequences. Will set it to %d",
+ iNumSeeds);
+ }
+
+
+ /* Turn sequences into vectors of distances to the seeds
+ *
+ */
+ piSeeds = (int *) CKMALLOC(iNumSeeds * sizeof(int));
+ if (0 != SeedSelection(piSeeds, iNumSeeds, SEED_SELECTION, prMSeq)) {
+ Log(&rLog, LOG_ERROR, "Something went wrong during seed selection for mbed");
+ return -1;
+ }
+ ppdSeqVec = (double **) CKMALLOC(prMSeq->nseqs * sizeof(double *));
+ for (iI=0; iI<prMSeq->nseqs; iI++) {
+ ppdSeqVec[iI] = (double *) CKMALLOC(iNumSeeds * sizeof(double));
+ }
+ if (0 != SeqToVec(ppdSeqVec, prMSeq, piSeeds, iNumSeeds, iPairDistType)) {
+ Log(&rLog, LOG_ERROR, "Could not convert sequences into vectors for mbed");
+ return -1;
+ }
+ CKFREE(piSeeds);
+
+
+ /* Calculate (pre-)clusters of sequence vectors by applying
+ * bisecting kmeans
+ *
+ */
+#if MBED_TIMING
+ /* start clock only here, to make sure we don't include pairwise
+ * distance computation */
+ StopwatchZero(stopwatch);
+ StopwatchStart(stopwatch);
+#endif
+
+ BisectingKmeans(&prKMeansResult, prMSeq->nseqs, iNumSeeds, ppdSeqVec,
+ MIN_REQUIRED_SEQ_PER_PRECLUSTER,
+ MAX_ALLOWED_SEQ_PER_PRECLUSTER);
+ Log(&rLog, LOG_INFO,
+ "mBed created %u cluster/s (with a minimum of %d and a soft maximum of %d sequences each)",
+ prKMeansResult->iNClusters,
+ MIN_REQUIRED_SEQ_PER_PRECLUSTER,
+ MAX_ALLOWED_SEQ_PER_PRECLUSTER);
+
+
+#if PRINT_CLUSTER_DISTRIBUTION
+ Log(&rLog, LOG_FORCED_DEBUG, "Bisecting Kmeans returned %d clusters", prKMeansResult->iNClusters);
+ for (iI=0; iI<prKMeansResult->iNClusters; iI++) {
+#if TRACE
+ int iD;
+ Log(&rLog, LOG_FORCED_DEBUG, "Diagnostic output for cluster %d follows:", iI);
+ fprintf(stderr, "DEBUG(%s|%s():%d): center coordinates =",
+ __FILE__, __FUNCTION__, __LINE__);
+ for (iD=0; iD<iNumSeeds; iD++) {
+ fprintf(stderr, " %f", prKMeansResult->ppdClusterCenters[iI][iD]);
+ }
+ fprintf(stderr, "\n");
+ fflush(stderr);
+#endif
+ Log(&rLog, LOG_FORCED_DEBUG, "Cluster %d has %d objects assigned",
+ iI, prKMeansResult->piNObjsPerCluster[iI]);
+#if TRACE
+ for (iJ=0; iJ<prKMeansResult->piNObjsPerCluster[iI]; iJ++) {
+ int iRealIndex = prKMeansResult->ppiObjIndicesPerCluster[iI][iJ];
+
+ Log(&rLog, LOG_FORCED_DEBUG, "Cluster %u: object %u has index %u (= seq %s)",
+ iI, iJ, iRealIndex, prMSeq->sqinfo[iRealIndex].name);
+ }
+#endif
+ }
+#endif
+
+
+ /* Cluster pre-clusters produced by k-means.
+ *
+ * Do this by calculating the vector distances of the cluster
+ * centers and applying UPGMA.
+ *
+ * @note could try to force-balance the tree here
+ *
+ */
+ if (0 != NewSymMatrix(&prPreClusterDistmat,
+ prKMeansResult->iNClusters, prKMeansResult->iNClusters)) {
+ Log(&rLog, LOG_FATAL, "%s", "Memory allocation for pre-cluster distance-matrix failed");
+ }
+ for (iI=0; iI<prKMeansResult->iNClusters; iI++) {
+ for (iJ=iI+1; iJ<prKMeansResult->iNClusters; iJ++) {
+ double dDist;
+ dDist = EuclDist(prKMeansResult->ppdClusterCenters[iI],
+ prKMeansResult->ppdClusterCenters[iJ],
+ iNumSeeds);
+ SymMatrixSetValue(prPreClusterDistmat, iI, iJ, dDist);
+ /* Log(&rLog, LOG_FORCED_DEBUG, "Euclidean distance between clusters %d and %d = %f",
+ iI, iJ, dDist); */
+ }
+ }
+
+
+ /* labels needed for the guide tree building routine only */
+ ppcLabels = (char **) CKMALLOC(prKMeansResult->iNClusters * sizeof(char*));
+ for (iI=0; iI<prKMeansResult->iNClusters; iI++) {
+ ppcLabels[iI] = (char *) CKMALLOC(32 * sizeof(char));
+ (void) snprintf(ppcLabels[iI], 32, "Subcluster-%u", iI);
+ }
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "Distance matrix for pre-cluster centers:");
+ SymMatrixPrint(prPreClusterDistmat, ppcLabels, NULL);
+#endif
+
+ GuideTreeUpgma(prMbedTree_p,
+ ppcLabels, prPreClusterDistmat, NULL);
+
+ for (iI=0; iI<prKMeansResult->iNClusters; iI++) {
+ CKFREE(ppcLabels[iI]);
+ }
+ CKFREE(ppcLabels);
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "Cluster-center guide-tree:");
+ LogTree(*prMbedTree_p);
+#endif
+
+
+
+ /* Now replace each leaf in the pre-cluster-center tree
+ * appropriately, i.e. with the corresponding sub-cluster.
+ *
+ * For each leaf/sub-cluster, create a distance matrix for the
+ * corresponding sequences. Use distances between vectors as
+ * approximated distances between sequences. Then create a
+ * guide-tree by applying UPGMA.
+ *
+ */
+
+ /* Get a mapping of (pre)cluster number and leaf-node indices in
+ * the cluster-center tree. We can add trees to prMbedTree_p
+ * because AppendTrees() guarantees that no other than the node to
+ * append to changes.
+ */
+ piClusterToTreeNode = (int*)
+ CKMALLOC(prKMeansResult->iNClusters * sizeof(int));
+ iNodeIndex = FirstDepthFirstNode(*prMbedTree_p);
+ do {
+ if (IsLeaf(iNodeIndex, *prMbedTree_p)) {
+ int iLeafId = GetLeafId(iNodeIndex, *prMbedTree_p);
+ piClusterToTreeNode[iLeafId] = iNodeIndex;
+ }
+ iNodeIndex = NextDepthFirstNode(iNodeIndex, *prMbedTree_p);
+ } while (NULL_NEIGHBOR != iNodeIndex);
+
+
+
+
+ /* Now step through all the leafs and replace them with the
+ * corresponding sub-trees
+ */
+ NewProgress(&prSubClusterDistanceProgress, LogGetFP(&rLog, LOG_INFO),
+ "Distance calculation within sub-clusters", bPrintCR);
+ /* for each cluster */
+ for (iClusterIndex=0;
+ iClusterIndex < prKMeansResult->iNClusters; iClusterIndex++) {
+ /* distance matrix for the sub-cluster */
+ symmatrix_t *prWithinClusterDistances = NULL;
+ int iNSeqInCluster;
+ tree_t *prSubClusterTree = NULL;
+
+ ProgressLog(prSubClusterDistanceProgress,
+ iClusterIndex, prKMeansResult->iNClusters, FALSE);
+
+#if FULL_WITHIN_CLUSTER_DISTANCES
+ mseq_t *prSubClusterMSeq;
+ int iPairDistType;
+ int iSeqIndex;
+ int iOldLogLevel;
+
+ Log(&rLog, LOG_DEBUG,
+ "%s\n", "Calling new Mbed use makes only sense if nseq>MAX_ALLOWED_SEQ_PER_PRECLUSTER");
+
+ if (TRUE == prMSeq->aligned) {
+ iPairDistType = PAIRDIST_SQUIDID_KIMURA;
+ } else {
+ iPairDistType = PAIRDIST_KTUPLE;
+ }
+#endif
+
+ iNSeqInCluster = prKMeansResult->piNObjsPerCluster[iClusterIndex];
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "#seqs in subcluster no %d = %d",
+ iClusterIndex, iNSeqInCluster);
+#endif
+
+#if FULL_WITHIN_CLUSTER_DISTANCES
+ /* create an mseq structure for sequences in this cluster
+ * don't need most of the members (e.g. orig_seq) but copy
+ * them anyway for the sake of completeness
+ */
+ NewMSeq(&prSubClusterMSeq);
+
+ prSubClusterMSeq->nseqs = iNSeqInCluster;
+ prSubClusterMSeq->seqtype = prMSeq->seqtype;
+ if (NULL!=prMSeq->filename) {
+ prSubClusterMSeq->filename = CkStrdup(prMSeq->filename);
+ }
+ prSubClusterMSeq->aligned = prMSeq->aligned; /* FS, r252 */
+ prSubClusterMSeq->seq = (char **)
+ CKMALLOC(prSubClusterMSeq->nseqs * sizeof(char *));
+ prSubClusterMSeq->orig_seq = (char **)
+ CKMALLOC(prSubClusterMSeq->nseqs * sizeof(char *));
+ prSubClusterMSeq->sqinfo = (SQINFO *)
+ CKMALLOC(prSubClusterMSeq->nseqs * sizeof(SQINFO));
+
+ for (iSeqIndex=0; iSeqIndex<iNSeqInCluster; iSeqIndex++) {
+ int iRealSeqIndex = prKMeansResult->ppiObjIndicesPerCluster[iClusterIndex][iSeqIndex];
+ prSubClusterMSeq->seq[iSeqIndex] = CkStrdup(prMSeq->seq[iRealSeqIndex]);
+ prSubClusterMSeq->orig_seq[iSeqIndex] = CkStrdup(prMSeq->orig_seq[iRealSeqIndex]);
+ SeqinfoCopy(&prSubClusterMSeq->sqinfo[iSeqIndex], &prMSeq->sqinfo[iRealSeqIndex]);
+#if TRACE
+ Log(&rLog, LOG_DEBUG, "seq no %d in cluster %d is %s (real index = %d)",
+ iSeqIndex, iClusterIndex, prSubClusterMSeq->sqinfo[iSeqIndex].name,
+ iRealSeqIndex);
+#endif
+ }
+#endif
+
+
+ /* Create a distance matrix for this sub-cluster
+ * (prWithinClusterDistances) by using the vector distances or
+ * ktuple distances.
+ *
+ * Then apply UPGMA to get a subcluster tree
+ * (prSubClusterTree) and append created tree to the
+ * pre-cluster-tree (prMbedTree_p)
+ *
+ */
+#if FULL_WITHIN_CLUSTER_DISTANCES
+
+ iOldLogLevel = rLog.iLogLevelEnabled;
+ rLog.iLogLevelEnabled = LOG_WARN;
+ /* compute distances, but be quiet */
+ if (PairDistances(&prWithinClusterDistances, prSubClusterMSeq, iPairDistType,
+ 0, prSubClusterMSeq->nseqs, 0, prSubClusterMSeq->nseqs,
+ NULL, NULL)) {
+ Log(&rLog, LOG_ERROR, "Couldn't compute pair distances");
+ return -1;
+ }
+ rLog.iLogLevelEnabled = iOldLogLevel;
+
+#if COMPUTE_WITHIN_SUBCLUSTER_AVERAGE
+ {
+ double dSum = 0.0;
+ for (iI=0; iI<prSubClusterMSeq->nseqs; iI++) {
+ for (iJ=iI+1; iJ<prSubClusterMSeq->nseqs; iJ++) {
+ dSum += SymMatrixGetValue(prWithinClusterDistances, iI, iJ);
+ }
+ }
+ Log(&rLog, LOG_FORCED_DEBUG,
+ "mean pair-wise distance within subcluster %d of %d = %f",
+ iClusterIndex, prKMeansResult->iNClusters, dSum/prSubClusterMSeq->nseqs);
+ }
+#endif
+
+#else
+ if (NewSymMatrix(&prWithinClusterDistances, iNSeqInCluster, iNSeqInCluster)!=0) {
+ Log(&rLog, LOG_FATAL, "%s", "Memory allocation for disparity matrix failed");
+ }
+ for (iI=0; iI<iNSeqInCluster; iI++) {
+ int iRealIndexI = prKMeansResult->ppiObjIndicesPerCluster[iClusterIndex][iI];
+
+ for (iJ=iI+1; iJ<iNSeqInCluster; iJ++) {
+ int iRealIndexJ = prKMeansResult->ppiObjIndicesPerCluster[iClusterIndex][iJ];
+ double dist;
+
+ /* Log(&rLog, LOG_FORCED_DEBUG, "Cluster %d: compute distance between %d:%s and %d:%s",
+ iClusterIndex, i, prMSeq->sqinfo[iRealIndexI].name,
+ iJ, prMSeq->sqinfo[iRealIndexJ].name); */
+
+ if (1 == USE_EUCLIDEAN_DISTANCE) {
+ dist = EuclDist(ppdSeqVec[iRealIndexI],
+ ppdSeqVec[iRealIndexJ], iNumSeeds);
+ } else {
+ dist = CosDist(ppdSeqVec[iRealIndexI],
+ ppdSeqVec[iRealIndexJ], iNumSeeds);
+ }
+ SymMatrixSetValue(prWithinClusterDistances, iI, iJ, dist);
+ }
+ }
+#endif
+
+ /* labels needed for the guide tree building routine only */
+ ppcLabels = (char**) CKMALLOC(iNSeqInCluster * sizeof(char*));
+ for (iI=0; iI<iNSeqInCluster; iI++) {
+#if FULL_WITHIN_CLUSTER_DISTANCES
+ ppcLabels[iI] = prSubClusterMSeq->sqinfo[iI].name;
+#else
+ int iRealIndex = prKMeansResult->ppiObjIndicesPerCluster[iClusterIndex][iI];
+ ppcLabels[iI] = prMSeq->sqinfo[iRealIndex].name;
+#endif
+ }
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "Distance matrix for seqs within sub cluster %d/%d",
+ iClusterIndex, prKMeansResult->iNClusters);
+ SymMatrixPrint(prWithinClusterDistances, ppcLabels, NULL);
+#endif
+
+
+ GuideTreeUpgma(&prSubClusterTree, ppcLabels,
+ prWithinClusterDistances, NULL);
+
+ CKFREE(ppcLabels); /* don't free members, they just point */
+#if 0
+ Log(&rLog, LOG_FORCED_DEBUG, "Cluster %d guide-tree:", iClusterIndex);
+ LogTree(prSubClusterTree);
+#endif
+
+
+ /* The guide tree id's (that point to the sequences) now start
+ * from 0, i.e. the association with the prMSeq numbering is
+ * broken and fixed in the following
+ */
+ for (iNodeIndex = 0; iNodeIndex < (int)GetNodeCount(prSubClusterTree); iNodeIndex++) {
+ if (IsLeaf(iNodeIndex, prSubClusterTree)) {
+ int iLeafId = GetLeafId(iNodeIndex, prSubClusterTree);
+ int iRealId = prKMeansResult->ppiObjIndicesPerCluster[iClusterIndex][iLeafId];
+#if 0
+ Log(&rLog, LOG_FORCED_DEBUG, "Correcting leaf node %d which has (wrong) id %d and name %s to id %d (prMSeq name %s)",
+ iNodeIndex, iLeafId,
+ GetLeafName(iNodeIndex, prSubClusterTree),
+ iRealId, prMSeq->sqinfo[iRealId].name);
+#endif
+ SetLeafId(prSubClusterTree, iNodeIndex, iRealId);
+ }
+ }
+
+
+ /* Append the newly created tree (prSubClusterTree) to the
+ * corresponding node index of prMbedTree_p.
+ */
+#if TRACE
+ Log(&rLog, LOG_FORCED_DEBUG, "Will join trees at leaf node %d = %s",
+ piClusterToTreeNode[iClusterIndex],
+ GetLeafName(piClusterToTreeNode[iClusterIndex], *prMbedTree_p));
+#endif
+
+ AppendTree(*prMbedTree_p,
+ piClusterToTreeNode[iClusterIndex],
+ prSubClusterTree);
+ /* Note: piClusterToTreeNode is still valid, because
+ * AppendTrees() guarantees that no other than the node to
+ * append to changes. */
+
+#if 0
+ Log(&rLog, LOG_FORCED_DEBUG, "%s", "prMbedTree_p after cluster %d has appended:", iClusterIndex);
+ LogTree(*prMbedTree_p);
+
+ if (0) {
+ char fname[] = "mbed-joined-tree.dnd";
+ FILE *pfOut;
+ if (NULL == (pfOut = fopen(fname, "w"))) {
+ Log(&rLog, LOG_FATAL, "Couldn't open %s for writing", fname);
+ }
+ MuscleTreeToFile(pfOut, *prMbedTree_p);
+ Log(&rLog, LOG_FORCED_DEBUG, "Joined tree written to %s", fname);
+ fclose(pfOut);
+ Log(&rLog, LOG_FATAL, "DEBUG EXIT");
+ }
+#endif
+
+ /* cleanup
+ */
+ FreeMuscleTree(prSubClusterTree);
+ FreeSymMatrix(&prWithinClusterDistances);
+#if FULL_WITHIN_CLUSTER_DISTANCES
+ FreeMSeq(&prSubClusterMSeq);
+#endif
+ } /* end for each cluster */
+ ProgressDone(prSubClusterDistanceProgress);
+ FreeProgress(&prSubClusterDistanceProgress);
+
+
+ if (NULL != pcGuidetreeOut) {
+ if (NULL == (pfOut = fopen(pcGuidetreeOut, "w"))) {
+ Log(&rLog, LOG_ERROR, "Couldn't open %s for writing", pcGuidetreeOut);
+ } else {
+ MuscleTreeToFile(pfOut, *prMbedTree_p);
+ Log(&rLog, LOG_INFO, "Guide tree written to %s", pcGuidetreeOut);
+ (void) fclose(pfOut);
+ }
+ }
+
+
+ /* cleanup
+ *
+ */
+#if MBED_TIMING
+ StopwatchStop(stopwatch);
+ StopwatchDisplay(stdout, "mBed time (without pairwise distance computation): ", stopwatch);
+ StopwatchFree(stopwatch);
+#endif
+
+ FreeKMeansResult(&prKMeansResult);
+ FreeSymMatrix(&prPreClusterDistmat);
+ for (iI=0; iI<prMSeq->nseqs; iI++) {
+ CKFREE(ppdSeqVec[iI]);
+ }
+ CKFREE(ppdSeqVec);
+ CKFREE(piClusterToTreeNode);
+
+#ifndef NDEBUG
+ TreeValidate(*prMbedTree_p);
+#endif
+
+ return 0;
+}
+/*** end: Mbed() ***/
git://source.jalview.org / jabaws.git / blobdiff