X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=clearcut.cpp;fp=clearcut.cpp;h=524683819e3adca6ffd8ffc96afdab27be8d3f39;hb=0caf3fbabaa3ece404f8ce77f4c883dc5b1bf1dc;hp=0000000000000000000000000000000000000000;hpb=1b73ff67c83892a025e597dabd9df6fe7b58206a;p=mothur.git

diff --git a/clearcut.cpp b/clearcut.cpp
new file mode 100644
index 0000000..5246838
--- /dev/null
+++ b/clearcut.cpp
@@ -0,0 +1,2159 @@
+
+/*
+ * 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 */
+}
+
+
+
+
+
+
+