/*****************************************************************
 * SQUID - a library of functions for biological sequence analysis
 * Copyright (C) 1992-2002 Washington University School of Medicine
 * 
 *     This source code is freely distributed under the terms of the
 *     GNU General Public License. See the files COPYRIGHT and LICENSE
 *     for details.
 *****************************************************************/

/* weight.c
 * SRE, Thu Mar  3 07:56:01 1994
 * 
 * Calculate weights for sequences in an alignment.
 * RCS $Id: weight.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: weight.c,v 1.9 2002/10/09 14:26:09 eddy Exp)
 */

#include <ctype.h>
#include <string.h>
#include "squid.h"
#include "sre_random.h"

static void upweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, int node);
static void downweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, 
		       float *fwt, int node);
static float simple_distance(char *s1, char *s2);
static int    simple_diffmx(char **aseqs,int num, float ***ret_dmx);

/* Function: GSCWeights()
 * 
 * Purpose:  Use Erik's tree-based algorithm to set weights for
 *           sequences in an alignment. upweight() and downweight()
 *           are derived from Graeme Mitchison's code.
 *           
 * Args:     aseq        - array of (0..nseq-1) aligned sequences
 *           nseq        - number of seqs in alignment  
 *           alen        - length of alignment
 *           wgt         - allocated [0..nseq-1] array of weights to be returned
 *           
 * Return:   (void)
 *           wgt is filled in.
 */
void
GSCWeights(char **aseq, int nseq, int alen, float *wgt)
{
  float **dmx;                 /* distance (difference) matrix */
  struct phylo_s *tree;
  float  *lwt, *rwt;           /* weight on left, right of this tree node */
  float  *fwt;                 /* final weight assigned to this node */
  int      i;
  
  /* Sanity check first
   */
  if (nseq == 1) { wgt[0] = 1.0; return; }

  /* I use a simple fractional difference matrix derived by
   * pairwise identity. Perhaps I should include a Poisson
   * distance correction.
   */
  MakeDiffMx(aseq, nseq, &dmx);
  if (! Cluster(dmx, nseq, CLUSTER_MIN, &tree))  Die("Cluster() failed");
  
  /* Allocations
   */
  lwt = MallocOrDie (sizeof(float) * (2 * nseq - 1));
  rwt = MallocOrDie (sizeof(float) * (2 * nseq - 1));
  fwt = MallocOrDie (sizeof(float) * (2 * nseq - 1));
  
  /* lwt and rwt are the total branch weight to the left and
   * right of a node or sequence. They are 0..2N-2.  0..N-1 are 
   * the sequences; these have weight 0. N..2N-2 are the actual
   * tree nodes.
   */
  for (i = 0; i < nseq; i++)
    lwt[i] = rwt[i] = 0.0;
				/* recursively calculate rwt, lwt, starting
				   at node nseq (the root) */
  upweight(tree, nseq, lwt, rwt, nseq);

				/* recursively distribute weight across the
				   tree */
  fwt[nseq] = nseq;
  downweight(tree, nseq, lwt, rwt, fwt, nseq);
				/* collect the weights */
  for (i = 0; i < nseq; i++)
    wgt[i]  = fwt[i];

  FMX2Free(dmx);
  FreePhylo(tree, nseq);
  free(lwt); free(rwt); free(fwt);
}

static void 
upweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, int node)
{
  int ld,rd;

  ld = tree[node-nseq].left;
  if (ld >= nseq) upweight(tree, nseq, lwt, rwt, ld);
  rd = tree[node-nseq].right;
  if (rd >= nseq) upweight(tree, nseq, lwt, rwt, rd);
  lwt[node] = lwt[ld] + rwt[ld] + tree[node-nseq].lblen;
  rwt[node] = lwt[rd] + rwt[rd] + tree[node-nseq].rblen;
}


static void 
downweight(struct phylo_s *tree, int nseq, float *lwt, float *rwt, float *fwt, int node)
{
  int ld,rd;
  float lnum, rnum;

  ld = tree[node-nseq].left;
  rd = tree[node-nseq].right;
  if (lwt[node] + rwt[node] > 0.0)
    {
      fwt[ld] = fwt[node] * (lwt[node] / (lwt[node] + rwt[node]));
      fwt[rd] = fwt[node] * (rwt[node] / (lwt[node] + rwt[node]));
    }
  else
    {
      lnum = (ld >= nseq) ? tree[ld-nseq].incnum : 1.0;
      rnum = (rd >= nseq) ? tree[rd-nseq].incnum : 1.0;
      fwt[ld] = fwt[node] * lnum / (lnum + rnum);
      fwt[rd] = fwt[node] * rnum / (lnum + rnum);
    }

  if (ld >= nseq) downweight(tree, nseq, lwt, rwt, fwt, ld);
  if (rd >= nseq) downweight(tree, nseq, lwt, rwt, fwt, rd);
}




/* Function: VoronoiWeights()
 * 
 * Purpose:  Calculate weights using the scheme of Sibbald &
 *           Argos (JMB 216:813-818 1990). The scheme is
 *           slightly modified because the original algorithm
 *           actually doesn't work on gapped alignments.
 *           The sequences are assumed to be protein.
 *           
 * Args:     aseq  - array of (0..nseq-1) aligned sequences
 *           nseq  - number of sequences
 *           alen  - length of alignment
 *           wgt   - allocated [0..nseq-1] array of weights to be returned
 *
 * Return:   void
 *           wgt is filled in.
 */
void
VoronoiWeights(char **aseq, int nseq, int alen, float *wgt)
{
  float **dmx;                 /* distance (difference) matrix    */
  float  *halfmin;             /* 1/2 minimum distance to other seqs */
  char   **psym;                /* symbols seen in each column     */
  int     *nsym;                /* # syms seen in each column      */
  int      symseen[27];         /* flags for observed syms         */
  char    *randseq;             /* randomly generated sequence     */
  int      acol;		/* pos in aligned columns          */
  int      idx;                 /* index in sequences              */
  int      symidx;              /* 0..25 index for symbol          */
  int      i;			/* generic counter                 */
  float   min;			/* minimum distance                */
  float   dist;		/* distance between random and real */
  float   challenge, champion; /* for resolving ties              */
  int     itscale;		/* how many iterations per seq     */
  int     iteration;           
  int     best;			/* index of nearest real sequence  */

  /* Sanity check first
   */
  if (nseq == 1) { wgt[0] = 1.0; return; }

  itscale = 50;

  /* Precalculate 1/2 minimum distance to other
   * sequences for each sequence
   */
  if (! simple_diffmx(aseq, nseq, &dmx)) 
    Die("simple_diffmx() failed");
  halfmin = MallocOrDie (sizeof(float) * nseq);
  for (idx = 0; idx < nseq; idx++)
    {
      for (min = 1.0, i = 0; i < nseq; i++)
	{
	  if (i == idx) continue;
	  if (dmx[idx][i] < min) min = dmx[idx][i];
	}
      halfmin[idx] = min / 2.0;
    }
  Free2DArray((void **) dmx, nseq);

  /* Set up the random sequence generating model.
   */
  psym = MallocOrDie (alen * sizeof(char *));
  nsym = MallocOrDie (alen * sizeof(int));
  for (acol = 0; acol < alen; acol++)
    psym[acol] = MallocOrDie (27 * sizeof(char));

/* #ifdef ORIGINAL_SIBBALD_ALGORITHM_IS_BROKEN */
  for (acol = 0; acol < alen; acol++)
    {
      memset(symseen, 0, sizeof(int) * 27);
      for (idx = 0; idx < nseq; idx++)
	if (! isgap(aseq[idx][acol]))
	  {
	    if (isupper((int) aseq[idx][acol])) 
	      symidx = aseq[idx][acol] - 'A';
	    else
	      symidx = aseq[idx][acol] - 'a';
	    if (symidx >= 0 && symidx < 26)
	      symseen[symidx] = 1;
	  }
	else
	  symseen[26] = 1;	/* a gap */

      for (nsym[acol] = 0, i = 0; i < 26; i++)
	if (symseen[i]) 
	  {
	    psym[acol][nsym[acol]] = 'A'+i;
	    nsym[acol]++;
	  }
      if (symseen[26]) { psym[acol][nsym[acol]] = ' '; nsym[acol]++; }
    }
/* #endif ORIGINAL_SIBBALD_ALGORITHM_IS_BROKEN */

  /* Note: the original Sibbald&Argos algorithm calls for
   * bounding the sampled space using a template-like random
   * sequence generator. However, this leads to one minor
   * and one major problem. The minor problem is that
   * exceptional amino acids in a column can have a
   * significant effect by altering the amount of sampled
   * sequence space; the larger the data set, the worse
   * this problem becomes. The major problem is that 
   * there is no reasonable way to deal with gaps.
   * Gapped sequences simply inhabit a different dimensionality
   * and it's pretty painful to imagine calculating Voronoi
   * volumes when the N in your N-space is varying.
   * Note that all the examples shown by Sibbald and Argos
   * are *ungapped* examples.
   * 
   * The best way I've found to circumvent this problem is
   * just not to bound the sampled space; count gaps as
   * symbols and generate completely random sequences.
   */
#ifdef ALL_SEQUENCE_SPACE
  for (acol = 0; acol < alen; acol++)
    {
      strcpy(psym[acol], "ACDEFGHIKLMNPQRSTVWY ");
      nsym[acol] = 21;
    }
#endif
  
  /* Sibbald and Argos algorithm:
   *   1) assign all seqs weight 0.
   *   2) generate a "random" sequence
   *   3) calculate distance to every other sequence
   *      (if we get a distance < 1/2 minimum distance
   *       to other real seqs, we can stop)
   *   4) if unique closest sequence, increment its weight 1.
   *      if multiple closest seq, choose one randomly    
   *   5) repeat 2-4 for lots of iterations
   *   6) normalize all weights to sum to nseq.
   */
  randseq = MallocOrDie ((alen+1) * sizeof(char));

  best = 42.;			/* solely to silence GCC uninit warnings. */
  FSet(wgt, nseq, 0.0);
  for (iteration = 0; iteration < itscale * nseq; iteration++)
    {
      for (acol = 0; acol < alen; acol++)
	randseq[acol] = (nsym[acol] == 0) ? ' ' : psym[acol][CHOOSE(nsym[acol])];
      randseq[acol] = '\0';

      champion = sre_random();
      for (min = 1.0, idx = 0; idx < nseq; idx++)
	{
	  dist = simple_distance(aseq[idx], randseq);
	  if (dist < halfmin[idx]) 
	    { 
	      best = idx; 
	      break;      
	    } 
	  if (dist < min)          
	    { champion = sre_random(); best = idx; min = dist; } 
	  else if (dist == min)    
	    { 
	      challenge = sre_random(); 
	      if (challenge > champion)
		{ champion = challenge; best = idx; min = dist; }
	    }
	}
      wgt[best] += 1.0;
    }

  for (idx = 0; idx < nseq; idx++)
    wgt[idx] = wgt[idx] / (float) itscale;

  free(randseq);
  free(nsym);
  free(halfmin);
  Free2DArray((void **) psym, alen);
}


/* Function: simple_distance()
 * 
 * Purpose:  For two identical-length null-terminated strings, return
 *           the fractional difference between them. (0..1)
 *           (Gaps don't count toward anything.)
 */
static float
simple_distance(char *s1, char *s2)
{
  int diff  = 0;
  int valid = 0;

  for (; *s1 != '\0'; s1++, s2++)
    {
      if (isgap(*s1) || isgap(*s2)) continue;
      if (*s1 != *s2) diff++;
      valid++;
    }
  return (valid > 0 ? ((float) diff / (float) valid) : 0.0);
}
    
/* Function: simple_diffmx()
 * 
 * Purpose:  Given a set of flushed, aligned sequences, construct
 *           an NxN fractional difference matrix using the
 *           simple_distance rule.
 *           
 * Args:     aseqs        - flushed, aligned sequences
 *           num          - number of aseqs
 *           ret_dmx      - RETURN: difference matrix (caller must free)
 *           
 * Return:   1 on success, 0 on failure.
 */
static int
simple_diffmx(char    **aseqs,
	      int       num,
	      float ***ret_dmx)
{
  float **dmx;                 /* RETURN: distance matrix           */
  int      i,j;			/* counters over sequences           */

  /* Allocate
   */
  if ((dmx = (float **) malloc (sizeof(float *) * num)) == NULL)
    Die("malloc failed");
  for (i = 0; i < num; i++)
    if ((dmx[i] = (float *) malloc (sizeof(float) * num)) == NULL)
      Die("malloc failed");

  /* Calculate distances, symmetric matrix
   */
  for (i = 0; i < num; i++)
    for (j = i; j < num; j++)
      dmx[i][j] = dmx[j][i] = simple_distance(aseqs[i], aseqs[j]);

  /* Return
   */
  *ret_dmx = dmx;
  return 1;
}



/* Function: BlosumWeights()
 * Date:     SRE, Fri Jul 16 17:33:59 1999 (St. Louis)
 * 
 * Purpose:  Assign weights to a set of aligned sequences
 *           using the BLOSUM rule:
 *             - do single linkage clustering at some pairwise identity
 *             - in each cluster, give each sequence 1/clustsize
 *               total weight.
 *
 *           The clusters have no pairwise link >= maxid. 
 *           
 *           O(N) in memory. Probably ~O(NlogN) in time; O(N^2)
 *           in worst case, which is no links between sequences
 *           (e.g., values of maxid near 1.0).
 * 
 * Args:     aseqs - alignment
 *           nseq  - number of seqs in alignment
 *           alen  - # of columns in alignment
 *           maxid - fractional identity (e.g. 0.62 for BLOSUM62)
 *           wgt   - [0..nseq-1] array of weights to be returned
 */               
void
BlosumWeights(char **aseqs, int nseq, int alen, float maxid, float *wgt)
{
  int  *c, nc;
  int  *nmem;        /* number of seqs in each cluster */
  int   i;           /* loop counter */

  SingleLinkCluster(aseqs, nseq, alen, maxid, &c, &nc);

  FSet(wgt, nseq, 1.0);
  nmem = MallocOrDie(sizeof(int) * nc);

  for (i = 0; i < nc;   i++) nmem[i] = 0;
  for (i = 0; i < nseq; i++) nmem[c[i]]++;
  for (i = 0; i < nseq; i++) wgt[i] = 1. / (float) nmem[c[i]];

  free(nmem);
  free(c);
  return;
}


/* Function: PositionBasedWeights()
 * Date:     SRE, Fri Jul 16 17:47:22 1999 [St. Louis]
 *
 * Purpose:  Implementation of Henikoff and Henikoff position-based
 *           weights (JMB 243:574-578, 1994) [Henikoff94b].
 *           
 *           A significant advantage of this approach that Steve and Jorja
 *           don't point out is that it is O(N) in memory, unlike
 *           many other approaches like GSC weights or Voronoi.
 *           
 *           A potential disadvantage that they don't point out
 *           is that in the theoretical limit of infinite sequences
 *           in the alignment, weights go flat: eventually every
 *           column has at least one representative of each of 20 aa (or 4 nt)
 *           in it.
 *           
 *           They also don't give a rule for how to handle gaps.
 *           The rule used here seems the obvious and sensible one
 *           (ignore them). This means that longer sequences
 *           initially get more weight; hence a "double
 *           normalization" in which the weights are first divided
 *           by sequence length (to compensate for that effect),
 *           then normalized to sum to nseq.
 *           
 * Limitations:
 *           Implemented in a way that's alphabet-independent:
 *           it uses the 26 upper case letters as "residues".
 *           Any alphabetic character in aseq is interpreted as
 *           a unique "residue" (case insensitively; lower case
 *           mapped to upper case). All other characters are
 *           interpreted as gaps.
 *           
 *           This way, we don't have to pass around any alphabet
 *           type info (DNA vs. RNA vs. protein) and don't have
 *           to deal with remapping IUPAC degenerate codes
 *           probabilistically. However, on the down side,
 *           a sequence with a lot of degenerate IUPAC characters
 *           will get an artifactually high PB weight.
 *           
 * Args:     aseq   - sequence alignment to weight
 *           nseq   - number of sequences in alignment
 *           alen   - length of alignment
 *           wgt    - RETURN: weights filled in (pre-allocated 0..nseq-1)
 *
 * Returns:  (void)
 *           wgt is allocated (0..nseq-1) by caller, and filled in here.
 */
void
PositionBasedWeights(char **aseq, int nseq, int alen, float *wgt)
{
  int  rescount[26];		/* count of A-Z residues in a column   */
  int  nres;			/* number of different residues in col */
  int  idx, pos;                /* indices into aseq                   */
  int  x;
  float norm;

  FSet(wgt, nseq, 0.0);
  for (pos = 0; pos < alen; pos++)
    {
      for (x = 0; x < 26; x++) rescount[x] = 0;
      for (idx = 0; idx < nseq; idx++)
	if (isalpha(aseq[idx][pos]))
	  rescount[toupper(aseq[idx][pos]) - 'A'] ++;

      nres = 0;
      for (x = 0; x < 26; x++) 
	if (rescount[x] > 0) nres++;

      for (idx = 0; idx < nseq; idx++)
	if (isalpha(aseq[idx][pos]))
	  wgt[idx] += 1. / (float) (nres * rescount[toupper(aseq[idx][pos]) - 'A']);
    }

  for (idx = 0; idx < nseq; idx++)
    wgt[idx] /= (float) DealignedLength(aseq[idx]);
  norm = (float) nseq / FSum(wgt, nseq);
  FScale(wgt, nseq, norm);
  return;
}




/* Function: FilterAlignment()
 * Date:     SRE, Wed Jun 30 09:19:30 1999 [St. Louis]
 * 
 * Purpose:  Constructs a new alignment by removing near-identical 
 *           sequences from a given alignment (where identity is 
 *           calculated *based on the alignment*).
 *           Does not affect the given alignment.
 *           Keeps earlier sequence, discards later one. 
 *           
 *           Usually called as an ad hoc sequence "weighting" mechanism.
 *           
 * Limitations:
 *           Unparsed Stockholm markup is not propagated into the
 *           new alignment.
 *           
 * Args:     msa      -- original alignment
 *           cutoff   -- fraction identity cutoff. 0.8 removes sequences > 80% id.
 *           ret_new  -- RETURN: new MSA, usually w/ fewer sequences
 *                         
 * Return:   (void)
 *           ret_new must be free'd by caller: MSAFree().
 */
void
FilterAlignment(MSA *msa, float cutoff, MSA **ret_new)
{
  int     nnew;			/* number of seqs in new alignment */
  int    *list;
  int    *useme;
  float   ident;
  int     i,j;
  int     remove;

				/* find which seqs to keep (list) */
				/* diff matrix; allow ragged ends */
  list  = MallocOrDie (sizeof(int) * msa->nseq);
  useme = MallocOrDie (sizeof(int) * msa->nseq);
  for (i = 0; i < msa->nseq; i++) useme[i] = FALSE;

  nnew = 0;
  for (i = 0; i < msa->nseq; i++)
    {
      remove = FALSE;
      for (j = 0; j < nnew; j++)
	{
	  ident = PairwiseIdentity(msa->aseq[i], msa->aseq[list[j]]);
	  if (ident > cutoff)
	    { 
	      remove = TRUE; 
	      printf("removing %12s -- fractional identity %.2f to %s\n", 
		     msa->sqname[i], ident,
		     msa->sqname[list[j]]);
	      break; 
	    }
	}
      if (remove == FALSE) {
	list[nnew++] = i;
	useme[i]     = TRUE;
      }
    }

  MSASmallerAlignment(msa, useme, ret_new);
  free(list);
  free(useme);
  return;
}


/* Function: SampleAlignment()
 * Date:     SRE, Wed Jun 30 10:13:56 1999 [St. Louis]
 * 
 * Purpose:  Constructs a new, smaller alignment by sampling a given
 *           number of sequences at random. Does not change the
 *           alignment nor the order of the sequences.
 *           
 *           If you ask for a sample that is larger than nseqs,
 *           it silently returns the original alignment.
 *           
 *           Not really a weighting method, but this is as good
 *           a place as any to keep it, since it's similar in
 *           construction to FilterAlignment().
 *           
 * Args:     msa      -- original alignment
 *           sample   -- number of sequences in new alignment (0 < sample <= nseq)
 *           ret_new  -- RETURN: new MSA 
 *                         
 * Return:   (void)
 *           ret_new must be free'd by caller: MSAFree().
 */
void
SampleAlignment(MSA *msa, int sample, MSA **ret_new)
{
  int    *list;                 /* array for random selection w/o replace */
  int    *useme;                /* array of flags 0..nseq-1: TRUE to use */
  int     i, idx;
  int     len;

  /* Allocations
   */
  list  = (int *) MallocOrDie (sizeof(int) * msa->nseq);
  useme = (int *) MallocOrDie (sizeof(int) * msa->nseq);
  for (i = 0; i < msa->nseq; i++)
    {
      list[i]  = i;
      useme[i] = FALSE;
    }

  /* Sanity check.
   */
  if (sample >= msa->nseq) sample = msa->nseq;

				/* random selection w/o replacement */
  for (len = msa->nseq, i = 0; i < sample; i++)
    {
      idx = CHOOSE(len);
      printf("chose %d: %s\n", list[idx], msa->sqname[list[idx]]);
      useme[list[idx]] = TRUE;
      list[idx] = list[--len];
    }

  MSASmallerAlignment(msa, useme, ret_new);
  free(list);
  free(useme);
  return;
}


/* Function: SingleLinkCluster()
 * Date:     SRE, Fri Jul 16 15:02:57 1999 [St. Louis]
 *
 * Purpose:  Perform simple single link clustering of seqs in a
 *           sequence alignment. A pairwise identity threshold
 *           defines whether two sequences are linked or not.
 *           
 *           Important: runs in O(N) memory, unlike standard
 *           graph decomposition algorithms that use O(N^2)
 *           adjacency matrices or adjacency lists. Requires
 *           O(N^2) time in worst case (which is when you have
 *           no links at all), O(NlogN) in "average"
 *           case, and O(N) in best case (when there is just
 *           one cluster in a completely connected graph.
 *           
 *           (Developed because hmmbuild could no longer deal
 *           with GP120, a 16,013 sequence alignment.)
 *           
 * Limitations: 
 *           CASE-SENSITIVE. Assumes aseq have been put into
 *           either all lower or all upper case; or at least,
 *           within a column, there's no mixed case.
 *           
 * Algorithm: 
 *           I don't know if this algorithm is published. I 
 *           haven't seen it in graph theory books, but that might
 *           be because it's so obvious that nobody's bothered.
 *           
 *           In brief, we're going to do a breadth-first search
 *           of the graph, and we're going to calculate links
 *           on the fly rather than precalculating them into
 *           some sort of standard adjacency structure.
 *           
 *           While working, we keep two stacks of maximum length N:
 *                a : list of vertices that are still unconnected.
 *                b : list of vertices that we've connected to 
 *                    in our current breadth level, but we haven't
 *                    yet tested for other connections to a.
 *           The current length (number of elements in) a and b are
 *           kept in na, nb.
 *                    
 *           We store our results in an array of length N:
 *                c : assigns each vertex to a component. for example
 *                    c[4] = 1 means that vertex 4 is in component 1.
 *                    nc is the number of components. Components
 *                    are numbered from 0 to nc-1. We return c and nc
 *                    to our caller.
 *                    
 *           The algorithm is:
 *           
 *           Initialisation: 
 *                a  <-- all the vertices
 *                na <-- N
 *                b  <-- empty set
 *                nb <-- 0
 *                nc <-- 0
 *                
 *           Then:
 *                while (a is not empty)
 *                  pop a vertex off a, push onto b
 *                  while (b is not empty)
 *                    pop vertex v off b
 *                    assign c[v] = nc
 *                    for each vertex w in a:
 *                       compare v,w. If w is linked to v, remove w
 *                       from a, push onto b.
 *                  nc++     
 *           q.e.d. :)       
 *
 * Args:     aseq   - aligned sequences
 *           nseq   - number of sequences in aseq
 *           alen   - alignment length
 *           maxid  - fractional identity threshold 0..1. if id >= maxid, seqs linked
 *           ret_c  - RETURN: 0..nseq-1 assignments of seqs to components (clusters)
 *           ret_nc - RETURN: number of components
 *
 * Returns:  void.
 *           ret_c is allocated here. Caller free's with free(*ret_c)
 */
void
SingleLinkCluster(char **aseq, int nseq, int alen, float maxid, 
		  int **ret_c, int *ret_nc)
{
  int *a, na;                   /* stack of available vertices */
  int *b, nb;                   /* stack of working vertices   */
  int *c;                       /* array of results            */
  int  nc;			/* total number of components  */
  int  v,w;			/* index of a working vertices */
  int  i;			/* loop counter */

  /* allocations and initializations
   */
  a = MallocOrDie (sizeof(int) * nseq);
  b = MallocOrDie (sizeof(int) * nseq);
  c = MallocOrDie (sizeof(int) * nseq);
  for (i = 0; i < nseq; i++) a[i] = i;
  na = nseq;
  nb = 0;
  nc = 0;

  /* Main algorithm
   */
  while (na > 0)
    {
      v = a[na-1]; na--;	/* pop a vertex off a, */
      b[nb] = v;   nb++;	/* and push onto b     */
      while (nb > 0)
	{
	  v    = b[nb-1]; nb--;	/* pop vertex off b          */
	  c[v] = nc;		/* assign it to component nc */
	  for (i = na-1; i >= 0; i--)/* backwards, becase of deletion/swapping we do*/
	    if (simple_distance(aseq[v], aseq[a[i]]) <= 1. - maxid) /* linked? */
	      {			
		w = a[i]; a[i] = a[na-1]; na--;	/* delete w from a (note swap) */
		b[nb] = w; nb++;                /* push w onto b */
	      }
	}
      nc++;
    }

  /* Cleanup and return
   */
  free(a);
  free(b);
  *ret_c  = c;
  *ret_nc = nc;
  return;
}