+++ /dev/null
-
-/*
- * clearcut.c
- *
- * $Id$
- *
- *****************************************************************************
- *
- * Copyright (c) 2004, Luke Sheneman
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- *
- * + Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * + Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in
- * the documentation and/or other materials provided with the
- * distribution.
- * + The names of its contributors may not be used to endorse or promote
- * products derived from this software without specific prior
- * written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- *
- *****************************************************************************
- *
- * An implementation of the Relaxed Neighbor-Joining algorithm
- * of Evans, J., Sheneman, L., and Foster, J.
- *
- *
- * AUTHOR:
- *
- * Luke Sheneman
- * sheneman@cs.uidaho.edu
- *
- */
-
-
-
-#include <stdio.h>
-#include <string.h>
-#include <limits.h>
-#include <stdlib.h>
-#include <math.h>
-#include <time.h>
-#include <sys/time.h>
-#include <float.h>
-
-#include "distclearcut.h"
-#include "dmat.h"
-#include "fasta.h"
-#include "cmdargs.h"
-#include "common.h"
-#include "clearcut.h"
-#include "prng.h"
-
-
-/*
- * main() -
- *
- * The entry point to the program.
- *
- */
-int clearcut_main(int argc, char *argv[]) {
-
- DMAT *dmat; /* The working distance matrix */
- DMAT *dmat_backup = NULL;/* A backup distance matrix */
- NJ_TREE *tree; /* The phylogenetic tree */
- NJ_ARGS *nj_args; /* Structure for holding command-line arguments */
- long int i;
-
- /* some variables for tracking time */
- struct timeval tv;
- unsigned long long startUs, endUs;
-
-
- /* check and parse supplied command-line arguments */
- nj_args = NJ_handle_args(argc, argv);
-
- if(!nj_args) {
- fprintf(stderr, "Clearcut: Error processing command-line arguments.\n");
- exit(-1);
- }
-
- /* for verbose reporting, print the random number seed to stdout */
- if(nj_args->verbose_flag) {
- printf("PRNG SEED: %d\n", nj_args->seed);
- }
-
- /* Initialize Mersenne Twister PRNG */
- init_genrand(nj_args->seed);
-
-
- switch(nj_args->input_mode) {
-
- /* If the input type is a distance matrix */
- case NJ_INPUT_MODE_DISTANCE:
-
- /* parse the distance matrix */
- dmat = NJ_parse_distance_matrix(nj_args);
- if(!dmat) {
- exit(-1);
- }
-
- break;
-
- /* If the input type is a multiple sequence alignment */
- case NJ_INPUT_MODE_ALIGNED_SEQUENCES:
-
- /* build a distance matrix from a multiple sequence alignment */
- dmat = NJ_build_distance_matrix(nj_args);
- if(!dmat) {
- fprintf(stderr, "Clearcut: Failed to build distance matrix from alignment.\n");
- exit(-1);
- }
-
- break;
-
- default:
-
- fprintf(stderr, "Clearcut: Could not determine how to process input\n");
- exit(-1);
- }
-
- /*
- * Output the computed distance matrix,
- * if the user specified one.
- */
- if(nj_args->matrixout) {
- NJ_output_matrix(nj_args, dmat);
- }
-
- /*
- * If we are going to generate multiple trees from
- * the same distance matrix, we need to make a backup
- * of the original distance matrix.
- */
- if(nj_args->ntrees > 1) {
- dmat_backup = NJ_dup_dmat(dmat);
- }
-
- /* process n trees */
- for(i=0;i<nj_args->ntrees;i++) {
-
- /*
- * If the user has specified matrix shuffling, we need
- * to randomize the distance matrix
- */
- if(nj_args->shuffle) {
- NJ_shuffle_distance_matrix(dmat);
- }
-
- /* RECORD THE PRECISE TIME OF THE START OF THE NEIGHBOR-JOINING */
- gettimeofday(&tv, NULL);
- startUs = ((unsigned long long) tv.tv_sec * 1000000ULL)
- + ((unsigned long long) tv.tv_usec);
-
-
- /*
- * Invoke either the Relaxed Neighbor-Joining algorithm (default)
- * or the "traditional" Neighbor-Joining algorithm
- */
- if(nj_args->neighbor) {
- tree = NJ_neighbor_joining(nj_args, dmat);
- } else {
- tree = NJ_relaxed_nj(nj_args, dmat);
- }
-
- if(!tree) {
- fprintf(stderr, "Clearcut: Failed to construct tree.\n");
- exit(0);
- }
-
- /* RECORD THE PRECISE TIME OF THE END OF THE NEIGHBOR-JOINING */
- gettimeofday(&tv, NULL);
- endUs = ((unsigned long long) tv.tv_sec * 1000000ULL)
- + ((unsigned long long) tv.tv_usec);
-
- /* print the time taken to perform the neighbor join */
- if(nj_args->verbose_flag) {
- if(nj_args->neighbor) {
- fprintf(stderr, "NJ tree built in %llu.%06llu secs\n",
- (endUs - startUs) / 1000000ULL,
- (endUs - startUs) % 1000000ULL);
- } else {
- fprintf(stderr, "RNJ tree built in %llu.%06llu secs\n",
- (endUs - startUs) / 1000000ULL,
- (endUs - startUs) % 1000000ULL);
- }
- }
-
- /* Output the neighbor joining tree here */
- NJ_output_tree(nj_args, tree, dmat, i);
-
- NJ_free_tree(tree); /* Free the tree */
- NJ_free_dmat(dmat); /* Free the working distance matrix */
-
- /*
- * If we need to do another iteration, lets re-initialize
- * our working distance matrix.
- */
- if(nj_args->ntrees > 1 && i<(nj_args->ntrees-1) ) {
- dmat = NJ_dup_dmat(dmat_backup);
- }
- }
-
- /* Free the backup distance matrix */
- if(nj_args->ntrees > 1) {
- NJ_free_dmat(dmat_backup);
- }
-
- /* If verbosity, describe where the tree output is */
- if(nj_args->verbose_flag) {
- if(nj_args->neighbor) {
- printf("NJ tree(s) in %s\n", nj_args->outfilename);
- } else {
- printf("Relaxed NJ tree(s) in %s\n", nj_args->outfilename);
- }
- }
-
- return 0;
-}
-
-
-
-
-
-/*
- * NJ_find_hmin() - Find minimum transformed values along horizontal
- *
- *
- * INPUTS:
- * -------
- * dmat -- The distance matrix
- * a -- The index of the specific taxon in the distance matrix
- *
- * RETURNS:
- * --------
- * <float> -- The value of the selected minimum
- * min -- Used to transport the index of the minima out
- * of the function (by reference)
- * hmincount -- Return the number of minima along the horizontal
- * (by reference)
- *
- *
- * DESCRIPTION:
- * ------------
- *
- * A fast, inline function to find the smallest transformed value
- * along the "horizontal" portion of an entry in a distance matrix.
- *
- * Distance matrices are stored internally as continguously-allocated
- * upper-diagonal structures. With the exception of the taxa at
- * row 0 of this upper-diagonal matrix, all taxa have both a horizontal
- * and vertical component in the distance matrix. This function
- * scans the horizonal portion of the entry in the distance matrix
- * for the specified taxon and finds the minimum transformed value
- * along that horizontal component.
- *
- * Since multiple minima can exist along the horizontal portion
- * of the entry, I consider all minima and break ties
- * stochastically to help avoid systematic bias.
- *
- * Just searching along the horizontal portion of a row is very fast
- * since the data is stored linearly and contiguously in memory and
- * cache locality is exploited in the distance matrix representation.
- *
- * Look at nj.h for more information on how the distance matrix
- * is architected.
- *
- */
-static inline
-float
-NJ_find_hmin(DMAT *dmat,
- long int a,
- long int *min,
- long int *hmincount) {
-
- long int i; /* index variable for looping */
- int size; /* current size of distance matrix */
- int mindex = 0; /* holds the current index to the chosen minimum */
- float curval; /* used to hold current transformed values */
- float hmin; /* the value of the transformed minimum */
-
- float *ptr, *r2, *val; /* pointers used to reduce dereferencing in inner loop */
-
- /* values used for stochastic selection among multiple minima */
- float p, x;
- long int smallcnt;
-
- /* initialize the min to something large */
- hmin = (float)HUGE_VAL;
-
- /* setup some pointers to limit dereferencing later */
- r2 = dmat->r2;
- val = dmat->val;
- size = dmat->size;
-
- /* initialize values associated with minima tie breaking */
- p = 1.0;
- smallcnt = 0;
-
-
- ptr = &(val[NJ_MAP(a, a+1, size)]); /* at the start of the horiz. part */
- for(i=a+1;i<size;i++) {
-
- curval = *(ptr++) - (r2[a] + r2[i]); /* compute transformed distance */
-
- if(NJ_FLT_EQ(curval, hmin)) { /* approx. equal */
-
- smallcnt++;
-
- p = 1.0/(float)smallcnt;
- x = genrand_real2();
-
- /* select this minimum in a way which is proportional to
- the number of minima found along the row so far */
- if( x < p ) {
- mindex = i;
- }
-
- } else if (curval < hmin) {
-
- smallcnt = 1;
- hmin = curval;
- mindex = i;
- }
- }
-
- /* save off the the minimum index to be returned via reference */
- *min = mindex;
-
- /* save off the number of minima */
- *hmincount = smallcnt;
-
- /* return the value of the smallest tranformed distance */
- return(hmin);
-}
-
-
-
-
-
-
-
-
-/*
- * NJ_find_vmin() - Find minimum transformed distance along vertical
- *
- *
- * INPUTS:
- * -------
- * dmat -- The distance matrix
- * a -- The index of the specific taxon in the distance matrix
- *
- *
- * RETURNS:
- * --------
- * <float> -- The value of the selected minimum
- * min -- Used to transport the index of the minima out
- * of the function (by reference)
- * vmincount -- The number of minima along the vertical
- * return by reference.
- *
- * DESCRIPTION:
- * ------------
- *
- * A fast, inline function to find the smallest transformed value
- * along the "vertical" portion of an entry in a distance matrix.
- *
- * Distance matrices are stored internally as continguously-allocated
- * upper-diagonal matrices. With the exception of the taxa at
- * row 0 of this upper-diagonal matrix, all taxa have both a horizontal
- * and vertical component in the distance matrix. This function
- * scans the vertical portion of the entry in the distance matrix
- * for the specified taxon and finds the minimum transformed value
- * along that vertical component.
- *
- * Since multiple minima can exist along the vertical portion
- * of the entry, I consider all minima and break ties
- * stochastically to help avoid systematic bias.
- *
- * Due to cache locality reasons, searching along the vertical
- * component is going to be considerably slower than searching
- * along the horizontal.
- *
- * Look at nj.h for more information on how the distance matrix
- * is architected.
- *
- */
-static inline
-float
-NJ_find_vmin(DMAT *dmat,
- long int a,
- long int *min,
- long int *vmincount) {
-
- long int i; /* index variable used for looping */
- long int size; /* track the size of the matrix */
- long int mindex = 0;/* track the index to the minimum */
- float curval; /* track value of current transformed distance */
- float vmin; /* the index to the smallest "vertical" minimum */
-
- /* pointers which are used to reduce pointer dereferencing in inner loop */
- float *ptr, *r2, *val;
-
- /* values used in stochastically breaking ties */
- float p, x;
- long int smallcnt;
-
- /* initialize the vertical min to something really big */
- vmin = (float)HUGE_VAL;
-
- /* save off some values to limit dereferencing later */
- r2 = dmat->r2;
- val = dmat->val;
- size = dmat->size;
-
- p = 1.0;
- smallcnt = 0;
-
- /* start on the first row and work down */
- ptr = &(val[NJ_MAP(0, a, size)]);
- for(i=0;i<a;i++) {
-
- curval = *ptr - (r2[i] + r2[a]); /* compute transformed distance */
-
- if(NJ_FLT_EQ(curval, vmin)) { /* approx. equal */
-
- smallcnt++;
-
- p = 1.0/(float)smallcnt;
- x = genrand_real2();
-
- /* break ties stochastically to avoid systematic bias */
- if( x < p ) {
- mindex = i;
- }
-
- } else if (curval < vmin) {
-
- smallcnt = 1;
- vmin = curval;
- mindex = i;
- }
-
- /* increment our working pointer to the next row down */
- ptr += size-i-1;
- }
-
- /* pass back the index to the minimum found so far (by reference) */
- *min = mindex;
-
- /* pass back the number of minima along the vertical */
- *vmincount = smallcnt;
-
- /* return the value of the smallest transformed distance */
- return(vmin);
-}
-
-
-
-
-/*
- * NJ_permute() - Generate random permutation using the provably unbiased
- * Fisher-Yates Shuffle.
- *
- * INPUTS:
- * -------
- * perm -- A pointer to the array of long ints which will be filled.
- * size -- the length of the permutation vector
- *
- *
- * OUTPUTS:
- * -------
- * NONE
- *
- *
- * DESCRIPTION:
- * ------------
- *
- * Return a permuted list of numbers from 0 through size.
- * This is accomplished by initializing the permutation
- * as an ordered list of integers and then iterating
- * through and swapping two values selected according to the
- * Fisher-Yates method.
- *
- * This unbiased method for random permutation generation is
- * discussed in:
- *
- * Donald E. Knuth, The Art of Computer Programming,
- * Addison-Wesley, Volumes 1, 2, and 3, 3rd edition, 1998
- *
- */
-static inline
-void
-NJ_permute(long int *perm,
- long int size) {
-
- long int i; /* index used for looping */
- long int swap; /* we swap values to generate permutation */
- long int tmp; /* used for swapping values */
-
-
- /* check to see if vector of long ints is valid */
- if(!perm) {
- fprintf(stderr, "Clearcut: NULL permutation pointer in NJ_permute()\n");
- exit(-1);
- }
-
- /* init permutation as an ordered list of integers */
- for(i=0;i<size;i++) {
- perm[i] = i;
- }
-
- /*
- * Iterate across the array from i = 0 to size -1, swapping ith element
- * with a randomly chosen element from a changing range of possible values
- */
- for(i=0;i<size;i++) {
-
- /* choose which element we will swap with */
- swap = i + NJ_genrand_int31_top(size-i);
-
- /* swap elements here */
- if(i != swap) {
- tmp = perm[swap];
- perm[swap] = perm[i];
- perm[i] = tmp;
- }
- }
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_compute_r() - Compute post-join changes to r-vector. In this
- * case, we decrement all of the accumulated distances
- * in the r-vector for the two nodes that were
- * recently joined (i.e. x, y)
- *
- * INPUTS:
- * -------
- * dmat -- The distance matrix
- * a -- The index of one of the taxa that were joined
- * b -- The index of the other taxa that was joined
- *
- * RETURNS:
- * --------
- * NONE
- *
- * DESCRIPTION:
- * ------------
- *
- * This vector of floats is used as a summary of overall distances from
- * each entry in the distance matrix to every other entry. These values
- * are then used when computing the transformed distances from which
- * decisions concerning joining are made.
- *
- * For speed, we don't recompute r from scratch. Instead, we decrement
- * all entries in r by the appropriate amount. That is, r[i] -= dist(i, a)
- * and r[i] -= dist(i, b).
- *
- * As a speed optimization, I process the rows altogether for cache locality
- * purposes, and then process columns.
- *
- * The processing of the scaled r matrix (r2) is handled on-the-fly elsewhere.
- *
- */
-static inline
-void
-NJ_compute_r(DMAT *dmat,
- long int a,
- long int b) {
-
- long int i; /* a variable used in indexing */
- float *ptrx, *ptry; /* pointers into the distance matrix */
-
- /* some variables to limit pointer dereferencing in loop */
- long int size;
- float *r, *val;
-
- /* to limit pointer dereferencing */
- size = dmat->size;
- val = dmat->val;
- r = dmat->r+a+1;
-
- /*
- * Loop through the rows and decrement the stored r values
- * by the distances stored in the rows and columns of the distance
- * matrix which are being removed post-join.
- *
- * We do the rows altogether in order to benefit from cache locality.
- */
- ptrx = &(val[NJ_MAP(a, a+1, size)]);
- ptry = &(val[NJ_MAP(b, b+1, size)]);
-
- for(i=a+1;i<size;i++) {
- *r -= *(ptrx++);
-
- if(i>b) {
- *r -= *(ptry++);
- }
-
- r++;
- }
-
- /* Similar to the above loop, we now do the columns */
- ptrx = &(val[NJ_MAP(0, a, size)]);
- ptry = &(val[NJ_MAP(0, b, size)]);
- r = dmat->r;
- for(i=0;i<b;i++) {
- if(i<a) {
- *r -= *ptrx;
- ptrx += size-i-1;
- }
-
- *r -= *ptry;
- ptry += size-i-1;
- r++;
- }
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_check_additivity() - Check to make sure that addivity preserved by join
- *
- *
- * INPUTS:
- * -------
- * dmat -- distance matrix
- * a -- index into dmat for one of the rows to be joined
- * b -- index into dmat for another row to be joined
- *
- * OUTPUTS:
- * --------
- * int 1 if join adheres to additivity constraint
- * 0 if join does breaks additivity
- *
- * DESCRIPTION:
- * ------------
- *
- * Here we perform the check to make sure that by joining a and b we do not
- * also break consistency (i.e. preserves additivity) with the distances between
- * the taxa in the new clade and other nodes in the tree. This is done quite
- * efficiently by looking up the untransformed distance between node b and
- * some other "target" taxa in the distance matrix (which is not a nor b) and
- * comparing that distance to the distance computed by finding the distance
- * from node a to the proposed internal node "x" which joins (a,b).
- *
- * If dist(x,b) + dist (b, target) == dist(b, target) then additivity is
- * preserved, otherwise, additivity is not preserved. If we are in
- * additivity mode, this join should be rejected.
- *
- */
-static inline
-int
-NJ_check_additivity(DMAT *dmat,
- long int a,
- long int b) {
-
- float a2clade, b2clade;
- float clade_dist;
- long int target;
-
-
- /* determine target taxon here */
- if(b == dmat->size-1) {
- /* if we can't do a row here, lets do a column */
- if(a==0) {
- if(b==1) {
- target = 2;
- } else {
- target = 1;
- }
- } else {
- target = 0;
- }
- } else {
- target = b+1;
- }
-
-
- /* distance between a and the root of clade (a,b) */
- a2clade =
- ( (dmat->val[NJ_MAP(a, b, dmat->size)]) +
- (dmat->r2[a] - dmat->r2[b]) ) / 2.0;
-
- /* distance between b and the root of clade (a,b) */
- b2clade =
- ( (dmat->val[NJ_MAP(a, b, dmat->size)]) +
- (dmat->r2[b] - dmat->r2[a]) ) / 2.0;
-
- /* distance between the clade (a,b) and the target taxon */
- if(b<target) {
-
- /* compute the distance from the clade root to the target */
- clade_dist =
- ( (dmat->val[NJ_MAP(a, target, dmat->size)] - a2clade) +
- (dmat->val[NJ_MAP(b, target, dmat->size)] - b2clade) ) / 2.0;
-
- /*
- * Check to see that distance from clade root to target + distance from
- * b to clade root are equal to the distance from b to the target
- */
- if(NJ_FLT_EQ(dmat->val[NJ_MAP(b, target, dmat->size)],
- (clade_dist + b2clade))) {
- return(1); /* join is legitimate */
- } else {
- return(0); /* join is illigitimate */
- }
-
- } else {
-
- /* compute the distance from the clade root to the target */
- clade_dist =
- ( (dmat->val[NJ_MAP(target, a, dmat->size)] - a2clade) +
- (dmat->val[NJ_MAP(target, b, dmat->size)] - b2clade) ) / 2.0;
-
- /*
- * Check to see that distance from clade root to target + distance from
- * b to clade root are equal to the distance from b to the target
- */
- if(NJ_FLT_EQ(dmat->val[NJ_MAP(target, b, dmat->size)],
- (clade_dist + b2clade))) {
- return(1); /* join is legitimate */
- } else {
- return(0); /* join is illegitimate */
- }
- }
-}
-
-
-
-
-
-
-
-/*
- * NJ_check() - Check to see if two taxa can be joined
-
- *
- * INPUTS:
- * -------
- * nj_args -- Pointer to the data structure holding command-line args
- * dmat -- distance matrix
- * a -- index into dmat for one of the rows to be joined
- * b -- index into dmat for another row to be joined
- * min -- the minimum value found
- * additivity -- a flag (0 = not additive mode, 1 = additive mode)
- *
- * OUTPUTS:
- * --------
- * int 1 if join is okay
- * 0 if join is not okay
- *
- * DESCRIPTION:
- * ------------
- *
- * This function ultimately takes two rows and makes sure that the
- * intersection of those two rows, which has a transformed distance of
- * "min", is actually the smallest (or equal to the smallest)
- * transformed distance for both rows (a, b). If so, it returns
- * 1, else it returns 0.
- *
- * Basically, we want to join two rows only if the minimum
- * transformed distance on either row is at the intersection of
- * those two rows.
- *
- */
-static inline
-int
-NJ_check(NJ_ARGS *nj_args,
- DMAT *dmat,
- long int a,
- long int b,
- float min,
- int additivity) {
-
-
- long int i, size;
- float *ptr, *val, *r2;
-
-
- /* some aliases for speed and readability reasons */
- val = dmat->val;
- r2 = dmat->r2;
- size = dmat->size;
-
-
- /* now determine if joining a, b will result in broken distances */
- if(additivity) {
- if(!NJ_check_additivity(dmat, a, b)) {
- return(0);
- }
- }
-
- /* scan the horizontal of row b, punt if anything < min */
- ptr = &(val[NJ_MAP(b, b+1, size)]);
- for(i=b+1;i<size;i++) {
- if( NJ_FLT_LT( (*ptr - (r2[b] + r2[i])), min) ) {
- return(0);
- }
- ptr++;
- }
-
- /* scan the vertical component of row a, punt if anything < min */
- if(nj_args->norandom) { /* if we are doing random joins, we checked this */
- ptr = val + a;
- for(i=0;i<a;i++) {
- if( NJ_FLT_LT( (*ptr - (r2[i] + r2[a])), min) ) {
- return(0);
- }
- ptr += size-i-1;
- }
- }
-
- /* scan the vertical component of row b, punt if anything < min */
- ptr = val + b;
- for(i=0;i<b;i++) {
- if( NJ_FLT_LT( (*ptr - (r2[i] + r2[b])), min) && i!=a) {
- return(0);
- }
- ptr += size-i-1;
- }
-
- return(1);
-}
-
-
-
-
-
-
-
-/*
- * NJ_collapse() - Collapse the distance matrix by removing
- * rows a and b from the distance matrix and
- * replacing them with a single new row which
- * represents the internal node joining a and b
- *
- *
- * INPUTS:
- * -------
- * dmat -- A pointer to the distance matrix
- * vertex -- A pointer to the vertex vector (vector of tree nodes)
- * which is used in constructing the tree
- * a -- An index to a row in the distance matrix from which we
- * joined. This row will be collapsed.
- * b -- An index to a row in the distance matrix from which we
- * joined. This row will be collapsed.
- *
- * RETURNS:
- * --------
- * NONE
- *
- *
- * DESCRIPTION:
- * ------------
- *
- * This function collapses the distance matrix in a way which optimizes
- * cache locality and ultimately gives us a speed improvement due to
- * cache. At this point, we've decided to join rows a and b from
- * the distance matrix. We will remove rows a and b from the distance
- * matrix and replace them with a new row which represents the internal
- * node which joins rows a and b together.
- *
- * We always keep the matrix as compact as possible in order to
- * get good performance from our cache in subsequent operations. Cache
- * is the key to good performance here.
- *
- * Key Steps:
- * ----------
- *
- * 1) Fill the "a" row with the new distances of the internal node
- * joining a and b to all other rows.
- * 2) Copy row 0 into what was row b
- * 3) Increment the pointer to the start of the distance matrix
- * by one row.
- * 4) Decrement the size of the matrix by one row.
- * 5) Do roughly the same thing to the r vector in order to
- * keep it in sync with the distance matrix.
- * 6) Compute the scaled r vector (r2) based on the updated
- * r vector
- *
- * This keeps the distance matrix as compact as possible in memory, and
- * is a relatively fast operation.
- *
- * This function requires that a < b
- *
- */
-static inline
-void
-NJ_collapse(DMAT *dmat,
- NJ_VERTEX *vertex,
- long int a,
- long int b) {
-
-
- long int i; /* index used for looping */
- long int size; /* size of dmat --> reduce pointer dereferencing */
- float a2clade; /* distance from a to the new node that joins a and b */
- float b2clade; /* distance from b to the new node that joins a and b */
- float cval; /* stores distance information during loop */
- float *vptr; /* pointer to elements in first row of dist matrix */
- float *ptra; /* pointer to elements in row a of distance matrix */
- float *ptrb; /* pointer to elements in row b of distance matrix */
-
- float *val, *r, *r2; /* simply used to limit pointer dereferencing */
-
-
- /* We must assume that a < b */
- if(a >= b) {
- fprintf(stderr, "Clearcut: (a<b) constraint check failed in NJ_collapse()\n");
- exit(0);
- }
-
- /* some shortcuts to help limit dereferencing */
- val = dmat->val;
- r = dmat->r;
- r2 = dmat->r2;
- size = dmat->size;
-
- /* compute the distance from the clade components (a, b) to the new node */
- a2clade =
- ( (val[NJ_MAP(a, b, size)]) + (dmat->r2[a] - dmat->r2[b]) ) / 2.0;
- b2clade =
- ( (val[NJ_MAP(a, b, size)]) + (dmat->r2[b] - dmat->r2[a]) ) / 2.0;
-
-
- r[a] = 0.0; /* we are removing row a, so clear dist. in r */
-
- /*
- * Fill the horizontal part of the "a" row and finish computing r and r2
- * we handle the horizontal component first to maximize cache locality
- */
- ptra = &(val[NJ_MAP(a, a+1, size)]); /* start ptra at the horiz. of a */
- ptrb = &(val[NJ_MAP(a+1, b, size)]); /* start ptrb at comparable place */
- for(i=a+1;i<size;i++) {
-
- /*
- * Compute distance from new internal node to others in
- * the distance matrix.
- */
- cval =
- ( (*ptra - a2clade) +
- (*ptrb - b2clade) ) / 2.0;
-
- /* incr. row b pointer differently depending on where i is in loop */
- if(i<b) {
- ptrb += size-i-1; /* traverse vertically by incrementing by row */
- } else {
- ptrb++; /* traverse horiz. by incrementing by column */
- }
-
- /* assign the newly computed distance and increment a ptr by a column */
- *(ptra++) = cval;
-
- /* accumulate the distance onto the r vector */
- r[a] += cval;
- r[i] += cval;
-
- /* scale r2 on the fly here */
- r2[i] = r[i]/(float)(size-3);
- }
-
- /* fill the vertical part of the "a" column and finish computing r and r2 */
- ptra = val + a; /* start at the top of the columb for "a" */
- ptrb = val + b; /* start at the top of the columb for "b" */
- for(i=0;i<a;i++) {
-
- /*
- * Compute distance from new internal node to others in
- * the distance matrix.
- */
- cval =
- ( (*ptra - a2clade) +
- (*ptrb - b2clade) ) / 2.0;
-
- /* assign the newly computed distance and increment a ptr by a column */
- *ptra = cval;
-
- /* accumulate the distance onto the r vector */
- r[a] += cval;
- r[i] += cval;
-
- /* scale r2 on the fly here */
- r2[i] = r[i]/(float)(size-3);
-
- /* here, always increment by an entire row */
- ptra += size-i-1;
- ptrb += size-i-1;
- }
-
-
- /* scale r2 on the fly here */
- r2[a] = r[a]/(float)(size-3);
-
-
-
- /*
- * Copy row 0 into row b. Again, the code is structured into two
- * loops to maximize cache locality for writes along the horizontal
- * component of row b.
- */
- vptr = val;
- ptrb = val + b;
- for(i=0;i<b;i++) {
- *ptrb = *(vptr++);
- ptrb += size-i-1;
- }
- vptr++; /* skip over the diagonal */
- ptrb = &(val[NJ_MAP(b, b+1, size)]);
- for(i=b+1;i<size;i++) {
- *(ptrb++) = *(vptr++);
- }
-
- /*
- * Collapse r here by copying contents of r[0] into r[b] and
- * incrementing pointer to the beginning of r by one row
- */
- r[b] = r[0];
- dmat->r = r+1;
-
-
- /*
- * Collapse r2 here by copying contents of r2[0] into r2[b] and
- * incrementing pointer to the beginning of r2 by one row
- */
- r2[b] = r2[0];
- dmat->r2 = r2+1;
-
- /* increment dmat pointer to next row */
- dmat->val += size;
-
- /* decrement the total size of the distance matrix by one row */
- dmat->size--;
-
- return;
-}
-
-
-
-
-
-
-
-
-
-/*
- * NJ_neighbor_joining() - Perform a traditional Neighbor-Joining
- *
- *
- * INPUTS:
- * -------
- * nj_args -- A pointer to a structure containing the command-line arguments
- * dmat -- A pointer to the distance matrix
- *
- * RETURNS:
- * --------
- * NJ_TREE * -- A pointer to the Neighbor-Joining tree.
- *
- * DESCRIPTION:
- * ------------
- *
- * This function performs a traditional Neighbor-Joining operation in which
- * the distance matrix is exhaustively searched for the global minimum
- * transformed distance. The two nodes which intersect at the global
- * minimum transformed distance are then joined and the distance
- * matrix is collapsed. This process continues until there are only
- * two nodes left, at which point those nodes are joined.
- *
- */
-NJ_TREE *
-NJ_neighbor_joining(NJ_ARGS *nj_args,
- DMAT *dmat) {
-
-
- NJ_TREE *tree = NULL;
- NJ_VERTEX *vertex = NULL;
-
- long int a, b;
- float min;
-
-
- /* initialize the r and r2 vectors */
- NJ_init_r(dmat);
-
- /* allocate and initialize our vertex vector used for tree construction */
- vertex = NJ_init_vertex(dmat);
- if(!vertex) {
- fprintf(stderr, "Clearcut: Could not initialize vertex in NJ_neighbor_joining()\n");
- return(NULL);
- }
-
- /* we iterate until the working distance matrix has only 2 entries */
- while(vertex->nactive > 2) {
-
- /*
- * Find the global minimum transformed distance from the distance matrix
- */
- min = NJ_min_transform(dmat, &a, &b);
-
- /*
- * Build the tree by removing nodes a and b from the vertex array
- * and inserting a new internal node which joins a and b. Collapse
- * the vertex array similarly to how the distance matrix and r and r2
- * are compacted.
- */
- NJ_decompose(dmat, vertex, a, b, 0);
-
- /* decrement the r and r2 vectors by the distances corresponding to a, b */
- NJ_compute_r(dmat, a, b);
-
- /* compact the distance matrix and the r and r2 vectors */
- NJ_collapse(dmat, vertex, a, b);
- }
-
- /* Properly join the last two nodes on the vertex list */
- tree = NJ_decompose(dmat, vertex, 0, 1, NJ_LAST);
-
- /* return the computed tree to the calling function */
- return(tree);
-}
-
-
-
-
-
-
-
-
-/*
- * NJ_relaxed_nj() - Construct a tree using the Relaxed Neighbor-Joining
- *
- * INPUTS:
- * -------
- * nj_args -- A pointer to a data structure containing the command-line args
- * dmat -- A pointer to the distance matrix
- *
- * RETURNS:
- * --------
- *
- * NJ_TREE * -- A pointer to a Relaxed Neighbor-Joining tree
- *
- * DESCRIPTION:
- * ------------
- *
- * This function implements the Relaxed Neighbor-Joining algorithm of
- * Evans, J., Sheneman, L., and Foster, J.
- *
- * Relaxed Neighbor-Joining works by choosing a local minimum transformed
- * distance when determining when to join two nodes. (Traditional
- * Neighbor-Joining chooses a global minimum transformed distance).
- *
- * The algorithm shares the property with traditional NJ that if the
- * input distances are additive (self-consistent), then the algorithm
- * will manage to construct the true tree consistent with the additive
- * distances. Additivity state is tracked and every proposed join is checked
- * to make sure it maintains additivity constraints. If no
- * additivity-preserving join is possible in a single pass, then the distance
- * matrix is non-additive, and additivity checking is abandoned.
- *
- * The algorithm will either attempt joins randomly, or it will perform joins
- * in a particular order. The default behavior is to perform joins randomly,
- * but this can be switched off with a command-line switch.
- *
- * For randomized joins, all attempts are made to alleviate systematic bias
- * for the choice of rows to joins. All tie breaking is done in a way which
- * is virtually free of bias.
- *
- * To perform randomized joins, a random permutation is constructed which
- * specifies the order in which to attempt joins. I iterate through the
- * random permutation, and for each row in the random permutation, I find
- * the minimum transformed distance for that row. If there are multiple
- * minima, I break ties evenly. For the row which intersects our
- * randomly chosen row at the chosen minimum, if we are are still in
- * additivity mode, I check to see if joining the two rows will break
- * our additivity constraints. If not, I check to see if there exists
- * a transformed distance which is smaller than the minimum found on the
- * original row. If there is, then we proceed through the random permutation
- * trying additional rows in the random order specified in the permutation.
- * If there is no smaller minimum transformed distance on either of the
- * two rows, then we join them, collapse the distance matrix, and compute
- * a new random permutation.
- *
- * If the entire random permutation is traversed and no joins are possible
- * due to additivity constraints, then the distance matrix is not
- * additive, and additivity constraint-checking is disabled.
- *
- */
-NJ_TREE *
-NJ_relaxed_nj(NJ_ARGS *nj_args,
- DMAT *dmat) {
-
-
- NJ_TREE *tree;
- NJ_VERTEX *vertex;
- long int a, b, t, bh, bv, i;
- float hmin, vmin, hvmin;
- float p, q, x;
- int join_flag;
- int additivity_mode;
- long int hmincount, vmincount;
- long int *permutation = NULL;
-
-
-
- /* initialize the r and r2 vectors */
- NJ_init_r(dmat);
-
- additivity_mode = 1;
-
- /* allocate the permutation vector, if we are in randomize mode */
- if(!nj_args->norandom) {
- permutation = (long int *)calloc(dmat->size, sizeof(long int));
- if(!permutation) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_relaxed_nj()\n");
- return(NULL);
- }
- }
-
- /* allocate and initialize our vertex vector used for tree construction */
- vertex = NJ_init_vertex(dmat);
-
- /* loop until there are only 2 nodes left to join */
- while(vertex->nactive > 2) {
-
- switch(nj_args->norandom) {
-
- /* RANDOMIZED JOINS */
- case 0:
-
- join_flag = 0;
-
- NJ_permute(permutation, dmat->size-1);
- for(i=0;i<dmat->size-1 && (vertex->nactive>2) ;i++) {
-
- a = permutation[i];
-
- /* find min trans dist along horiz. of row a */
- hmin = NJ_find_hmin(dmat, a, &bh, &hmincount);
- if(a) {
- /* find min trans dist along vert. of row a */
- vmin = NJ_find_vmin(dmat, a, &bv, &vmincount);
- } else {
- vmin = hmin;
- bv = bh;
- vmincount = 0;
- }
-
- if(NJ_FLT_EQ(hmin, vmin)) {
-
- /*
- * The minima along the vertical and horizontal are
- * the same. Compute the proportion of minima along
- * the horizonal (p) and the proportion of minima
- * along the vertical (q).
- *
- * If the same minima exist along the horizonal and
- * vertical, we break the tie in a way which is
- * non-biased. That is, we break the tie based on the
- * proportion of horiz. minima versus vertical minima.
- *
- */
- p = (float)hmincount / ((float)hmincount + (float)vmincount);
- q = 1.0 - p;
- x = genrand_real2();
-
- if(x < p) {
- hvmin = hmin;
- b = bh;
- } else {
- hvmin = vmin;
- b = bv;
- }
- } else if(NJ_FLT_LT(hmin, vmin) ) {
- hvmin = hmin;
- b = bh;
- } else {
- hvmin = vmin;
- b = bv;
- }
-
- if(NJ_check(nj_args, dmat, a, b, hvmin, additivity_mode)) {
-
- /* swap a and b, if necessary, to make sure a < b */
- if(b < a) {
- t = a;
- a = b;
- b = t;
- }
-
- join_flag = 1;
-
- /* join taxa from rows a and b */
- NJ_decompose(dmat, vertex, a, b, 0);
-
- /* collapse matrix */
- NJ_compute_r(dmat, a, b);
- NJ_collapse(dmat, vertex, a, b);
-
- NJ_permute(permutation, dmat->size-1);
- }
- }
-
- /* turn off additivity if go through an entire cycle without joining */
- if(!join_flag) {
- additivity_mode = 0;
- }
-
- break;
-
-
-
- /* DETERMINISTIC JOINS */
- case 1:
-
- join_flag = 0;
-
- for(a=0;a<dmat->size-1 && (vertex->nactive > 2) ;) {
-
- /* find the min along the horizontal of row a */
- hmin = NJ_find_hmin(dmat, a, &b, &hmincount);
-
- if(NJ_check(nj_args, dmat, a, b, hmin, additivity_mode)) {
-
- join_flag = 1;
-
- /* join taxa from rows a and b */
- NJ_decompose(dmat, vertex, a, b, 0);
-
- /* collapse matrix */
- NJ_compute_r(dmat, a, b);
- NJ_collapse(dmat, vertex, a, b);
-
- if(a) {
- a--;
- }
-
- } else {
- a++;
- }
- }
-
- /* turn off additivity if go through an entire cycle without joining */
- if(!join_flag) {
- additivity_mode = 0;
- }
-
- break;
- }
-
- } /* WHILE */
-
- /* Join the last two nodes on the vertex list */
- tree = NJ_decompose(dmat, vertex, 0, 1, NJ_LAST);
-
- if(nj_args->verbose_flag) {
- if(additivity_mode) {
- printf("Tree is additive\n");
- } else {
- printf("Tree is not additive\n");
- }
- }
-
- if(vertex) {
- NJ_free_vertex(vertex);
- }
-
- if(!nj_args->norandom && permutation) {
- free(permutation);
- }
-
- return(tree);
-}
-
-
-
-
-
-
-/*
- * NJ_print_distance_matrix() -
- *
- * Print a distance matrix
- *
- */
-void
-NJ_print_distance_matrix(DMAT *dmat) {
-
- long int i, j;
-
- printf("ntaxa: %ld\n", dmat->ntaxa);
- printf(" size: %ld\n", dmat->size);
-
- for(i=0;i<dmat->size;i++) {
-
- for(j=0;j<dmat->size;j++) {
- if(j>i) {
- printf(" %0.4f", dmat->val[NJ_MAP(i, j, dmat->size)]);
- } else {
- printf(" -");
- }
- }
-
-
- if(dmat->r && dmat->r2) {
- printf("\t\t%0.4f", dmat->r[i]);
- printf("\t%0.4f", dmat->r2[i]);
-
-
-
- printf("\n");
-
- for(j=0;j<dmat->size;j++) {
- if(j>i) {
- printf(" %0.4f", dmat->val[NJ_MAP(i, j, dmat->size)] - (dmat->r2[i] + dmat->r2[j]));
- } else {
- printf(" ");
- }
- }
-
- printf("\n");
- }
- }
-
- printf("\n\n");
-
- return;
-}
-
-
-
-
-
-
-
-/*
- * NJ_output_tree() -
- *
- * A wrapper for the function that really prints the tree,
- * basically to get a newline in there conveniently. :-)
- *
- * Print n trees, as specified in command-args
- * using "count" variable from 0 to (n-1)
- *
- */
-void
-NJ_output_tree(NJ_ARGS *nj_args,
- NJ_TREE *tree,
- DMAT *dmat,
- long int count) {
-
- FILE *fp;
-
- if(nj_args->stdout_flag) {
- fp = stdout;
- } else {
-
- if(count == 0) {
- fp = fopen(nj_args->outfilename, "w"); /* open for writing */
- } else {
- fp = fopen(nj_args->outfilename, "a"); /* open for appending */
- }
-
- if(!fp) {
- fprintf(stderr, "Clearcut: Failed to open outfile %s\n", nj_args->outfilename);
- exit(-1);
- }
- }
-
- NJ_output_tree2(fp, nj_args, tree, tree, dmat);
- fprintf(fp, ";\n");
-
- if(!nj_args->stdout_flag) {
- fclose(fp);
- }
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_output_tree2() -
- *
- *
- */
-void
-NJ_output_tree2(FILE *fp,
- NJ_ARGS *nj_args,
- NJ_TREE *tree,
- NJ_TREE *root,
- DMAT *dmat) {
-
- if(!tree) {
- return;
- }
-
- if(tree->taxa_index != NJ_INTERNAL_NODE) {
-
- if(nj_args->expblen) {
- fprintf(fp, "%s:%e",
- dmat->taxaname[tree->taxa_index],
- tree->dist);
- } else {
- fprintf(fp, "%s:%f",
- dmat->taxaname[tree->taxa_index],
- tree->dist);
- }
-
- } else {
-
-
- if(tree->left && tree->right) {
- fprintf(fp, "(");
- }
- if(tree->left) {
- NJ_output_tree2(fp, nj_args, tree->left, root, dmat);
- }
-
- if(tree->left && tree->right) {
- fprintf(fp, ",");
- }
- if(tree->right) {
- NJ_output_tree2(fp, nj_args, tree->right, root, dmat);
- }
-
- if(tree != root->left) {
- if(tree->left && tree->right) {
- if(tree != root) {
- if(nj_args->expblen) {
- fprintf(fp, "):%e", tree->dist);
- } else {
- fprintf(fp, "):%f", tree->dist);
- }
- } else {
- fprintf(fp, ")");
- }
- }
- } else {
- fprintf(fp, ")");
- }
- }
-
- return;
-}
-
-
-
-
-
-
-
-/*
- * NJ_init_r()
- *
- * This function computes the r column in our matrix
- *
- */
-void
-NJ_init_r(DMAT *dmat) {
-
- long int i, j, size;
- long int index;
- float *r, *r2, *val;
- long int size1;
- float size2;
-
- r = dmat->r;
- r2 = dmat->r2;
- val = dmat->val;
- size = dmat->size;
- size1 = size-1;
- size2 = (float)(size-2);
-
- index = 0;
- for(i=0;i<size1;i++) {
- index++;
- for(j=i+1;j<size;j++) {
- r[i] += val[index];
- r[j] += val[index];
- index++;
- }
-
- r2[i] = r[i]/size2;
- }
-
- return;
-}
-
-
-
-
-
-
-
-
-
-
-
-/*
- * NJ_init_vertex() -
- *
- * Construct a vertex, which we will use to construct our tree
- * in a true bottom-up approach. The vertex construct is
- * basically the center node in the initial star topology.
- *
- */
-NJ_VERTEX *
-NJ_init_vertex(DMAT *dmat) {
-
- long int i;
- NJ_VERTEX *vertex;
-
- /* allocate the vertex here */
- vertex = (NJ_VERTEX *)calloc(1, sizeof(NJ_VERTEX));
-
- /* allocate the nodes in the vertex */
- vertex->nodes = (NJ_TREE **)calloc(dmat->ntaxa, sizeof(NJ_TREE *));
- vertex->nodes_handle = vertex->nodes;
-
- /* initialize our size and active variables */
- vertex->nactive = dmat->ntaxa;
- vertex->size = dmat->ntaxa;
-
- /* initialize the nodes themselves */
- for(i=0;i<dmat->ntaxa;i++) {
-
- vertex->nodes[i] = (NJ_TREE *)calloc(1, sizeof(NJ_TREE));
-
- vertex->nodes[i]->left = NULL;
- vertex->nodes[i]->right = NULL;
-
- vertex->nodes[i]->taxa_index = i;
- }
-
- return(vertex);
-}
-
-
-
-
-
-/*
- * NJ_decompose() -
- *
- * This function decomposes the star by creating new internal nodes
- * and joining two existing tree nodes to it
- *
- */
-NJ_TREE *
-NJ_decompose(DMAT *dmat,
- NJ_VERTEX *vertex,
- long int x,
- long int y,
- int last_flag) {
-
- NJ_TREE *new_node;
- float x2clade, y2clade;
-
- /* compute the distance from the clade components to the new node */
- if(last_flag) {
- x2clade =
- (dmat->val[NJ_MAP(x, y, dmat->size)]);
- } else {
- x2clade =
- (dmat->val[NJ_MAP(x, y, dmat->size)])/2 +
- ((dmat->r2[x] - dmat->r2[y])/2);
- }
-
- vertex->nodes[x]->dist = x2clade;
-
- if(last_flag) {
- y2clade =
- (dmat->val[NJ_MAP(x, y, dmat->size)]);
- } else {
- y2clade =
- (dmat->val[NJ_MAP(x, y, dmat->size)])/2 +
- ((dmat->r2[y] - dmat->r2[x])/2);
- }
-
- vertex->nodes[y]->dist = y2clade;
-
- /* allocate new node to connect two sub-clades */
- new_node = (NJ_TREE *)calloc(1, sizeof(NJ_TREE));
-
- new_node->left = vertex->nodes[x];
- new_node->right = vertex->nodes[y];
- new_node->taxa_index = NJ_INTERNAL_NODE; /* this is not a terminal node, no taxa index */
-
- if(last_flag) {
- return(new_node);
- }
-
- vertex->nodes[x] = new_node;
- vertex->nodes[y] = vertex->nodes[0];
-
- vertex->nodes = &(vertex->nodes[1]);
-
- vertex->nactive--;
-
- return(new_node);
-}
-
-
-
-/*
- * NJ_print_vertex() -
- *
- * For debugging, print the contents of the vertex
- *
- */
-void
-NJ_print_vertex(NJ_VERTEX *vertex) {
-
- long int i;
-
- printf("Number of active nodes: %ld\n", vertex->nactive);
-
- for(i=0;i<vertex->nactive;i++) {
- printf("%ld ", vertex->nodes[i]->taxa_index);
- }
- printf("\n");
-
- return;
-}
-
-
-
-
-
-
-
-
-
-/*
- * NJ_print_r() -
- *
- */
-void
-NJ_print_r(DMAT *dmat) {
-
- long int i;
-
- printf("\n");
- for(i=0;i<dmat->size;i++) {
- printf("r[%ld] = %0.2f\n", i, dmat->r[i]);
- }
- printf("\n");
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_print_taxanames() -
- *
- * Print taxa names here
- *
- */
-void
-NJ_print_taxanames(DMAT *dmat) {
-
- long int i;
-
- printf("Number of taxa: %ld\n", dmat->ntaxa);
-
- for(i=0;i<dmat->ntaxa;i++) {
- printf("%ld) %s\n", i, dmat->taxaname[i]);
- }
-
- printf("\n");
-
- return;
-}
-
-
-
-
-/*
- * NJ_shuffle_distance_matrix() -
- *
- * Randomize a distance matrix here
- *
- */
-void
-NJ_shuffle_distance_matrix(DMAT *dmat) {
-
-
- long int *perm = NULL;
- char **tmp_taxaname = NULL;
- float *tmp_val = NULL;
- long int i, j;
-
-
- /* alloc the random permutation and a new matrix to hold the shuffled vals */
- perm = (long int *)calloc(dmat->size, sizeof(long int));
- tmp_taxaname = (char **)calloc(dmat->size, sizeof(char *));
- tmp_val = (float *)calloc(NJ_NCELLS(dmat->ntaxa), sizeof(float));
- if(!tmp_taxaname || !perm || !tmp_val) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_shuffle_distance_matrix()\n");
- exit(-1);
- }
-
- /* compute a permutation which will describe how to shuffle the matrix */
- NJ_permute(perm, dmat->size);
-
- for(i=0;i<dmat->size;i++) {
- for(j=i+1;j<dmat->size;j++) {
-
- if(perm[j] < perm[i]) {
- tmp_val[NJ_MAP(i, j, dmat->size)] = dmat->val[NJ_MAP(perm[j], perm[i], dmat->size)];
- } else {
- tmp_val[NJ_MAP(i, j, dmat->size)] = dmat->val[NJ_MAP(perm[i], perm[j], dmat->size)];
- }
-
- }
-
- tmp_taxaname[i] = dmat->taxaname[perm[i]];
- }
-
- /* free our random permutation */
- if(perm) {
- free(perm);
- }
-
- /* free the old value matrix */
- if(dmat->val) {
- free(dmat->val);
- }
-
- /* re-assign the value matrix pointers */
- dmat->val = tmp_val;
- dmat->valhandle = dmat->val;
-
- /*
- * Free our old taxaname with its particular ordering
- * and re-assign to the new.
- */
- if(dmat->taxaname) {
- free(dmat->taxaname);
- }
- dmat->taxaname = tmp_taxaname;
-
- return;
-}
-
-
-
-/*
- * NJ_free_tree() -
- *
- * Free a given NJ tree
- */
-void
-NJ_free_tree(NJ_TREE *node) {
-
- if(!node) {
- return;
- }
-
- if(node->left) {
- NJ_free_tree(node->left);
- }
-
- if(node->right) {
- NJ_free_tree(node->right);
- }
-
- free(node);
-
- return;
-}
-
-
-
-
-
-
-
-
-
-/*
- * NJ_print_permutation()
- *
- * Print a permutation
- *
- */
-void
-NJ_print_permutation(long int *perm,
- long int size) {
-
- long int i;
-
- for(i=0;i<size-1;i++) {
- printf("%ld,", perm[i]);
- }
- printf("%ld\n", perm[size-1]);
-
- return;
-}
-
-
-
-
-/*
- * NJ_dup_dmat() -
- *
- * Duplicate a distance matrix
- *
- */
-DMAT *
-NJ_dup_dmat(DMAT *src) {
-
- long int i;
- DMAT *dest;
-
- /* allocate the resulting distance matrix */
- dest = (DMAT *)calloc(1, sizeof(DMAT));
- if(!dest) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_dup_dmat()\n");
- goto XIT_BAD;
- }
-
- dest->ntaxa = src->ntaxa;
- dest->size = src->size;
-
- /* allocate space for array of pointers to taxanames */
- dest->taxaname = (char **)calloc(dest->ntaxa, sizeof(char *));
- if(!dest->taxaname) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_dup_dmat()\n");
- goto XIT_BAD;
- }
-
- /* allocate space for the taxanames themselves */
- for(i=0;i<src->ntaxa;i++) {
- dest->taxaname[i] = (char *)calloc(strlen(src->taxaname[i])+1, sizeof(char));
- if(!dest->taxaname[i]) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_dup_dmat()\n");
- goto XIT_BAD;
- }
- }
-
- /* allocate space for the distance values */
- dest->val = (float *)calloc(NJ_NCELLS(src->ntaxa), sizeof(float));
- if(!dest->val) {
- fprintf(stderr, "Clearcut: Memory allocation error in NJ_dup_dmat()\n");
- goto XIT_BAD;
- }
-
- /* allocate space for the r and r2 vectors */
- dest->r = (float *)calloc(src->ntaxa, sizeof(float));
- dest->r2 = (float *)calloc(src->ntaxa, sizeof(float));
-
- /* copy titles */
- for(i=0;i<src->ntaxa;i++) {
- strcpy(dest->taxaname[i], src->taxaname[i]);
- }
-
- /* copy values */
- memcpy(dest->val, src->valhandle, NJ_NCELLS(src->ntaxa)*sizeof(float));
-
- /* copy r and r2 */
- memcpy(dest->r, src->rhandle, src->ntaxa*sizeof(float));
- memcpy(dest->r2, src->r2handle, src->ntaxa*sizeof(float));
-
- /* track some memory addresses */
- dest->valhandle = dest->val;
- dest->rhandle = dest->r;
- dest->r2handle = dest->r2;
-
- return(dest);
-
- XIT_BAD:
-
- /* free what we may have allocated */
- NJ_free_dmat(dest);
-
- return(NULL);
-}
-
-
-
-
-/*
- * NJ_free_dmat() -
- */
-void
-NJ_free_dmat(DMAT *dmat) {
-
- long int i;
-
- if(dmat) {
-
- if(dmat->taxaname) {
-
- for(i=0;i<dmat->ntaxa;i++) {
- if(dmat->taxaname[i]) {
- free(dmat->taxaname[i]);
- }
- }
-
- free(dmat->taxaname);
- }
-
- if(dmat->valhandle) {
- free(dmat->valhandle);
- }
-
- if(dmat->rhandle) {
- free(dmat->rhandle);
- }
-
- if(dmat->r2handle) {
- free(dmat->r2handle);
- }
-
- free(dmat);
- }
-
- return;
-}
-
-
-
-
-
-/*
- * NJ_free_vertex() -
- *
- * Free the vertex data structure
- *
- */
-void
-NJ_free_vertex(NJ_VERTEX *vertex) {
-
- if(vertex) {
- if(vertex->nodes_handle) {
- free(vertex->nodes_handle);
- }
- free(vertex);
- }
-
- return;
-}
-
-
-
-
-
-
-
-
-
-/*
- *
- * NJ_min_transform() - Find the smallest transformed value to identify
- * which nodes to join.
- *
- * INPUTS:
- * -------
- * dmat -- The distance matrix
- *
- * RETURNS:
- * --------
- * <float> -- The minimimum transformed distance
- * ret_i -- The row of the smallest transformed distance (by reference)
- * ret_j -- The col of the smallest transformed distance (by reference)
- *
- *
- * DESCRIPTION:
- * ------------
- *
- * Used only with traditional Neighbor-Joining, this function checks the entire
- * working distance matrix and identifies the smallest transformed distance.
- * This requires traversing the entire diagonal matrix, which is itself a
- * O(N^2) operation.
- *
- */
-float
-NJ_min_transform(DMAT *dmat,
- long int *ret_i,
- long int *ret_j) {
-
- long int i, j; /* indices used for looping */
- long int tmp_i = 0;/* to limit pointer dereferencing */
- long int tmp_j = 0;/* to limit pointer dereferencing */
- float smallest; /* track the smallest trans. dist */
- float curval; /* the current trans. dist in loop */
-
- float *ptr; /* pointer into distance matrix */
- float *r2; /* pointer to r2 matrix for computing transformed dists */
-
- smallest = (float)HUGE_VAL;
-
- /* track these here to limit pointer dereferencing in inner loop */
- ptr = dmat->val;
- r2 = dmat->r2;
-
- /* for every row */
- for(i=0;i<dmat->size;i++) {
- ptr++; /* skip diagonal */
- for(j=i+1;j<dmat->size;j++) { /* for every column */
-
- /* find transformed distance in matrix at i, j */
- curval = *(ptr++) - (r2[i] + r2[j]);
-
- /* if the transformed distanance is less than the known minimum */
- if(curval < smallest) {
-
- smallest = curval;
- tmp_i = i;
- tmp_j = j;
- }
- }
- }
-
- /* pass back (by reference) the coords of the min. transformed distance */
- *ret_i = tmp_i;
- *ret_j = tmp_j;
-
- return(smallest); /* return the min transformed distance */
-}
-
-
-
-
-
-
-
projects / mothur.git / blobdiff