bme.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: bme.cpp,v $
 * PRODUCTION Revision 1000.1  2004/06/01 18:09:50  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: bme.cpp,v 1000.1 2004/06/01 18:09:50 gouriano Exp $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's official duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the author in any work or product based on this material.
*
* ===========================================================================
*
* Author:  Richard Desper
*
* File Description:  bme.cpp
*
*    A part of the Miminum Evolution algorithm
*
*/

#include <ncbi_pch.hpp>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "graph.h"
#include "fastme.h"

BEGIN_NCBI_SCOPE
BEGIN_SCOPE(fastme)

meEdge *siblingEdge(meEdge *e);
boolean leaf(meNode *v);
meEdge *depthFirstTraverse(meTree *T, meEdge *e);
meEdge *topFirstTraverse(meTree *T, meEdge *e);

void BalWFext(meEdge *e, double **A) /*works except when e is the one edge
				  inserted to new vertex v by firstInsert*/
{
  meEdge *f, *g;
  if ((leaf(e->head)) && (leaf(e->tail)))
    e->distance = A[e->head->index][e->head->index];
  else if (leaf(e->head))
    {
      f = e->tail->parentEdge;
      g = siblingEdge(e);
      e->distance = 0.5*(A[e->head->index][g->head->index]
			 + A[e->head->index][f->head->index]
			 - A[g->head->index][f->head->index]);
    }
  else
    {
      f = e->head->leftEdge;
      g = e->head->rightEdge;
      e->distance = 0.5*(A[g->head->index][e->head->index]
			 + A[f->head->index][e->head->index]
			 - A[f->head->index][g->head->index]);
    }
}

void BalWFint(meEdge *e, double **A)
{
  int up, down, left, right;
  up = e->tail->index;
  down = (siblingEdge(e))->head->index;
  left = e->head->leftEdge->head->index;
  right = e->head->rightEdge->head->index;
  e->distance = 0.25*(A[up][left] + A[up][right] + A[left][down] + A[right][down]) - 0.5*(A[down][up] + A[left][right]);
} 

  

void assignBMEWeights(meTree *T, double **A)
{
  meEdge *e;
  e = depthFirstTraverse(T,NULL);
  while (NULL != e) {
    if ((leaf(e->head)) || (leaf(e->tail)))
      BalWFext(e,A);
    else
      BalWFint(e,A);
    e = depthFirstTraverse(T,e);
  }
}      

void BMEcalcDownAverage(meTree *T, meNode *v, meEdge *e, double **D, double **A)
{
  meEdge  *left, *right;
  if (leaf(e->head))
    A[e->head->index][v->index] = D[v->index2][e->head->index2]; 
  else
    {
      left = e->head->leftEdge;
      right = e->head->rightEdge;
      A[e->head->index][v->index] = 0.5 * A[left->head->index][v->index] 
	+ 0.5 * A[right->head->index][v->index];
    }
}

void BMEcalcUpAverage(meTree *T, meNode *v, meEdge *e, double **D, double **A)
{
  meEdge *up,*down;
  if (T->root == e->tail)
    A[v->index][e->head->index] = D[v->index2][e->tail->index2];
  /*for now, use convention
    v->index first => looking up
    v->index second => looking down */
  else
    {
      up = e->tail->parentEdge;
      down = siblingEdge(e);
      A[v->index][e->head->index] = 0.5 * A[v->index][up->head->index]
	+0.5  * A[down->head->index][v->index];
    }
}


void BMEcalcNewvAverages(meTree *T, meNode *v, double **D, double **A)
{
  /*loop over edges*/
  /*depth-first search*/
  meEdge *e;
  e = NULL;
  e = depthFirstTraverse(T,e);  /*the downward averages need to be
				  calculated from bottom to top */
  while(NULL != e)
    {
      BMEcalcDownAverage(T,v,e,D,A);
      e = depthFirstTraverse(T,e);
    }
  
  e = topFirstTraverse(T,e);   /*the upward averages need to be calculated 
				 from top to bottom */
  while(NULL != e)
    {
      BMEcalcUpAverage(T,v,e,D,A);
      e = topFirstTraverse(T,e);
    }
}


/*update Pair updates A[nearEdge][farEdge] and makes recursive call to subtree
  beyond farEdge*/
/*root is head or tail of meEdge being split, depending on direction toward
  v*/
void updatePair(double **A, meEdge *nearEdge, meEdge *farEdge, meNode *v,
		meNode *root, double dcoeff, int direction)
{
  meEdge *sib;
  switch(direction) /*the various cases refer to where the new vertex has
		      been inserted, in relation to the meEdge nearEdge*/
    {
    case UP: /*this case is called when v has been inserted above 
	       or skew to farEdge*/
      /*do recursive calls first!*/
      if (NULL != farEdge->head->leftEdge)
	updatePair(A,nearEdge,farEdge->head->leftEdge,v,root,dcoeff,UP);
      if (NULL != farEdge->head->rightEdge)
	updatePair(A,nearEdge,farEdge->head->rightEdge,v,root,dcoeff,UP);
      A[farEdge->head->index][nearEdge->head->index] =
	A[nearEdge->head->index][farEdge->head->index]
	= A[farEdge->head->index][nearEdge->head->index]
	+ dcoeff*A[farEdge->head->index][v->index]
	- dcoeff*A[farEdge->head->index][root->index];
      break; 
    case DOWN: /*called when v has been inserted below farEdge*/
      if (NULL != farEdge->tail->parentEdge)
	updatePair(A,nearEdge,farEdge->tail->parentEdge,v,root,dcoeff,DOWN);
      sib = siblingEdge(farEdge);
      if (NULL != sib)
	updatePair(A,nearEdge,sib,v,root,dcoeff,UP);
      A[farEdge->head->index][nearEdge->head->index] =
	A[nearEdge->head->index][farEdge->head->index]
	= A[farEdge->head->index][nearEdge->head->index]
	+ dcoeff*A[v->index][farEdge->head->index]
	- dcoeff*A[farEdge->head->index][root->index];	  
    }
}

void updateSubTree(double **A, meEdge *nearEdge, meNode *v, meNode *root,
		   meNode *newNode, double dcoeff, int direction)
{
  meEdge *sib;
  switch(direction)
    {
    case UP: /*newNode is above the meEdge nearEdge*/
      A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
      A[newNode->index][nearEdge->head->index] = 
	A[nearEdge->head->index][newNode->index] =
	A[nearEdge->head->index][root->index];	
      if (NULL != nearEdge->head->leftEdge)
	updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff, UP);
      if (NULL != nearEdge->head->rightEdge)
	updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff, UP);
      updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
      break;
    case DOWN: /*newNode is below the meEdge nearEdge*/
      A[nearEdge->head->index][v->index] = A[v->index][nearEdge->head->index];
      A[newNode->index][nearEdge->head->index] = 
	A[nearEdge->head->index][newNode->index] =
	0.5*(A[nearEdge->head->index][root->index] 
	     + A[v->index][nearEdge->head->index]);
      sib = siblingEdge(nearEdge);
      if (NULL != sib)
	updateSubTree(A, sib, v, root, newNode, 0.5*dcoeff, SKEW);
      if (NULL != nearEdge->tail->parentEdge)
	updateSubTree(A, nearEdge->tail->parentEdge, v, root, newNode, 0.5*dcoeff, DOWN);
      updatePair(A, nearEdge, nearEdge, v, root, dcoeff, DOWN);
      break;
    case SKEW: /*newNode is neither above nor below nearEdge*/
      A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
      A[newNode->index][nearEdge->head->index] = 
	A[nearEdge->head->index][newNode->index] =
	0.5*(A[nearEdge->head->index][root->index] + 
	     A[nearEdge->head->index][v->index]);	
      if (NULL != nearEdge->head->leftEdge)
	updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff,SKEW);
      if (NULL != nearEdge->head->rightEdge)
	updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff,SKEW);
      updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
    }
}


/*we update all the averages for nodes (u1,u2), where the insertion point of 
  v is in "direction" from both u1 and u2 */
/*The general idea is to proceed in a direction from those edges already corrected
 */

/*r is the root of the meTree relative to the inserted node*/

void BMEupdateAveragesMatrix(double **A, meEdge *e, meNode *v,meNode *newNode)
{
  meEdge *sib, *par, *left, *right;
  /*first, update the v,newNode entries*/
  A[newNode->index][newNode->index] = 0.5*(A[e->head->index][e->head->index]