开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
gme.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：23k
源码类别：

生物技术

开发平台：

C/C++

gme.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: gme.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:09:56  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: gme.cpp,v 1000.1 2004/06/01 18:09:56 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Author:  Richard Desper
  36. *
  37. * File Description:  gme.cpp
  38. *
  39. *    A part of the Miminum Evolution algorithm
  40. *
  41. */
  43. #include <ncbi_pch.hpp>
  44. #include <stdio.h>
  45. #include <stdlib.h>
  46. #include <math.h>
  47. #include "graph.h"
  48. #include "fastme.h"
  50. BEGIN_NCBI_SCOPE
  51. BEGIN_SCOPE(fastme)
  53. meEdge *siblingEdge(meEdge *e);
  54. boolean leaf(meNode *v);
  55. meEdge *depthFirstTraverse(meTree *T, meEdge *e);
  56. meEdge *topFirstTraverse(meTree *T, meEdge *e);
  58. /*from NNI.c*/
  59. void fillTableUp(meEdge *e, meEdge *f, double **A, double **D, meTree *T);
  61. /*from graph.c*/
  62. void updateSizes(meEdge *e, int direction);
  64. /*OLSint and OLSext use the average distance array to calculate weights
  65.   instead of using the meEdge average weight fields*/
  67. void OLSext(meEdge *e, double **A)
  68. {
  69.   meEdge *f, *g;
  70.   if(leaf(e->head))
  71.     {
  72.       f = siblingEdge(e);
  73.       e->distance = 0.5*(A[e->head->index][e->tail->index] 
  74.  + A[e->head->index][f->head->index]
  75.  - A[f->head->index][e->tail->index]);
  76.     }
  77.   else
  78.     {
  79.       f = e->head->leftEdge;
  80.       g = e->head->rightEdge;
  81.       e->distance = 0.5*(A[e->head->index][f->head->index]
  82.  + A[e->head->index][g->head->index]
  83.  - A[f->head->index][g->head->index]);
  84.     }
  85. }
  87. double wf(double lambda, double D_LR, double D_LU, double D_LD, 
  88.    double D_RU, double D_RD, double D_DU)
  89. {
  90.   double weight;
  91.   weight = 0.5*(lambda*(D_LU  + D_RD) + (1 -lambda)*(D_LD + D_RU)
  92. - (D_LR + D_DU));  
  93.   return(weight);
  94. }
  96. void OLSint(meEdge *e, double **A)
  97. {
  98.   double lambda;
  99.   meEdge *left, *right, *sib;
  100.   left = e->head->leftEdge;
  101.   right = e->head->rightEdge;
  102.   sib = siblingEdge(e);
  103.   lambda = ((double) sib->bottomsize*left->bottomsize + 
  104.     right->bottomsize*e->tail->parentEdge->topsize) /
  105.     (e->bottomsize*e->topsize);
  106.   e->distance = wf(lambda,A[left->head->index][right->head->index],
  107.    A[left->head->index][e->tail->index],
  108.    A[left->head->index][sib->head->index],
  109.    A[right->head->index][e->tail->index],
  110.    A[right->head->index][sib->head->index],
  111.    A[sib->head->index][e->tail->index]);
  112. }
  115. void assignOLSWeights(meTree *T, double **A)
  116. {
  117.   meEdge *e;
  118.   e = depthFirstTraverse(T,NULL);
  119.   while (NULL != e) {
  120.     if ((leaf(e->head)) || (leaf(e->tail)))
  121.       OLSext(e,A);
  122.     else
  123.       OLSint(e,A);
  124.     e = depthFirstTraverse(T,e);
  125.   }
  126. }      
  128. /*makes table of average distances from scratch*/
  129. void makeOLSAveragesTable(meTree *T, double **D, double **A)
  130. {
  131.   meEdge *e, *f, *g, *h;
  132.   meEdge *exclude;
  133.   e = f = NULL;
  134.   e = depthFirstTraverse(T,e);
  135.   while (NULL != e)
  136.     {
  137.       f = e;
  138.       exclude = e->tail->parentEdge;
  139.       /*we want to calculate A[e->head][f->head] for all edges
  140. except those edges which are ancestral to e.  For those
  141. edges, we will calculate A[e->head][f->head] to have a
  142. different meaning, later*/
  143.       if(leaf(e->head))
  144. while (NULL != f)
  145.   {
  146.     if (exclude != f)    
  147.       {
  148. if (leaf(f->head))
  149.   A[e->head->index][f->head->index] = A[f->head->index][e->head->index] = D[e->head->index2][f->head->index2];
  150. else
  151.   {
  152.     g = f->head->leftEdge;
  153.     h = f->head->rightEdge;
  154.     A[e->head->index][f->head->index] = A[f->head->index][e->head->index] = (g->bottomsize*A[e->head->index][g->head->index] + h->bottomsize*A[e->head->index][h->head->index])/f->bottomsize; 
  155.   }
  156.       }
  157.     else /*exclude == f*/
  158.       exclude = exclude->tail->parentEdge; 
  159.     f = depthFirstTraverse(T,f);
  160.   }
  161.       else 
  162. /*e->head is not a leaf, so we go recursively to values calculated for
  163.   the nodes below e*/
  164. while(NULL !=f )
  165.   {
  166.     if (exclude != f)       
  167.       {
  168. g = e->head->leftEdge;
  169. h = e->head->rightEdge;
  170. A[e->head->index][f->head->index] = A[f->head->index][e->head->index] = (g->bottomsize*A[f->head->index][g->head->index] + h->bottomsize*A[f->head->index][h->head->index])/e->bottomsize;
  171.       }
  172.     else
  173.       exclude = exclude->tail->parentEdge;
  174.     f = depthFirstTraverse(T,f);
  175.   }
  177.   /*now we move to fill up the rest of the table: we want
  178.     A[e->head->index][f->head->index] for those cases where e is an
  179.     ancestor of f, or vice versa.  We'll do this by choosing e via a
  180.     depth first-search, and the recursing for f up the path to the
  181.     root*/
  182.       f = e->tail->parentEdge;
  183.       if (NULL != f)
  184. fillTableUp(e,f,A,D,T);    
  185.       e = depthFirstTraverse(T,e);
  186.     } 
  188.   /*we are indexing this table by vertex indices, but really the
  189.     underlying object is the meEdge set.  Thus, the array is one element
  190.     too big in each direction, but we'll ignore the entries involving the root,
  191.     and instead refer to each meEdge by the head of that edge.  The head of
  192.     the root points to the meEdge ancestral to the rest of the tree, so
  193.     we'll keep track of up distances by pointing to that head*/
  195.   /*10/13/2001: collapsed three depth-first searces into 1*/
  196. }
  199. void GMEcalcDownAverage(meNode *v, meEdge *e, double **D, double **A)
  200. {
  201.   meEdge *left, *right;
  202.   if (leaf(e->head))
  203.     A[e->head->index][v->index] = D[v->index2][e->head->index2]; 
  204.   else
  205.     {
  206.       left = e->head->leftEdge;
  207.       right = e->head->rightEdge;
  208.       A[e->head->index][v->index] = 
  209. ( left->bottomsize * A[left->head->index][v->index] + 
  210.   right->bottomsize * A[right->head->index][v->index]) 
  211. / e->bottomsize;
  212.     }
  213. }
  215. void GMEcalcUpAverage(meNode *v, meEdge *e, double **D, double **A)
  216. {
  217.   meEdge *up, *down;
  218.   if (NULL == e->tail->parentEdge)
  219.     A[v->index][e->head->index] =  D[v->index2][e->tail->index2];
  220.   else
  221.     {
  222.       up = e->tail->parentEdge;
  223.       down = siblingEdge(e);
  224.       A[v->index][e->head->index] = 
  225. (up->topsize * A[v->index][up->head->index] + 
  226.  down->bottomsize * A[down->head->index][v->index])
  227. / e->topsize;
  228.       }
  229. }
  231. /*this function calculates average distance D_Xv for each X which is
  232.   a meSet of leaves of an induced submeTree of T*/
  233. void GMEcalcNewvAverages(meTree *T, meNode *v, double **D, double **A)
  234. {
  235.   /*loop over edges*/
  236.   /*depth-first search*/
  237.   meEdge *e;
  238.   e = NULL;
  239.   e = depthFirstTraverse(T,e);  /*the downward averages need to be
  240.   calculated from bottom to top */
  241.   while(NULL != e)
  242.     {
  243.       GMEcalcDownAverage(v,e,D,A);
  244.       e = depthFirstTraverse(T,e);
  245.     }
  246.   
  247.   e = topFirstTraverse(T,e);   /*the upward averages need to be calculated 
  248.  from top to bottom */
  249.   while(NULL != e)
  250.     {
  251.       GMEcalcUpAverage(v,e,D,A);
  252.       e = topFirstTraverse(T,e);
  253.     }
  254. }
  256. double wf4(double lambda, double lambda2, double D_AB, double D_AC, 
  257.    double D_BC, double D_Av, double D_Bv, double D_Cv)
  258. {
  259.   return(((1 - lambda) * (D_AC + D_Bv) + (lambda2 - 1)*(D_AB + D_Cv)
  260.  + (lambda - lambda2)*(D_BC + D_Av)));
  261. }
  264. /*testEdge cacluates what the OLS weight would be if v were inserted into
  265.   T along e.  Compare against known values for inserting along 
  266.   f = e->parentEdge */
  267. /*edges are tested by a top-first, left-first scheme. we presume
  268.   all distances are fixed to the correct weight for 
  269.   e->parentEdge, if e is a left-oriented edge*/
  270. void testEdge(meEdge *e, meNode *v, double **A)
  271. {
  272.   double lambda, lambda2;
  273.   meEdge *par, *sib;
  274.   sib = siblingEdge(e);
  275.   par = e->tail->parentEdge;
  276.   /*C is meSet above e->tail, B is meSet below e, and A is meSet below sib*/
  277.   /*following the nomenclature of Desper & Gascuel*/
  278.   lambda =  (((double) (sib->bottomsize + e->bottomsize*par->topsize))
  279.      / ((1 + par->topsize)*(par->bottomsize)));
  280.   lambda2 = (((double) (sib->bottomsize + e->bottomsize*par->topsize))
  281.      / ((1 + e->bottomsize)*(e->topsize)));
  282.   e->totalweight = par->totalweight 
  283.     + wf4(lambda,lambda2,A[e->head->index][sib->head->index],
  284.   A[sib->head->index][e->tail->index],
  285.   A[e->head->index][e->tail->index],
  286.   A[sib->head->index][v->index],A[e->head->index][v->index],
  287.   A[v->index][e->tail->index]);  
  288. }
  290. void printDoubleTable(double **A, int d)
  291. {
  292.   int i,j;
  293.   for(i=0;i<d;i++)
  294.     {
  295.       for(j=0;j<d;j++)
  296. printf("%lf ", A[i][j]);
  297.       printf("n");
  298.     }
  299. }
  301. void GMEsplitEdge(meTree *T, meNode *v, meEdge *e, double **A);
  303. meTree *GMEaddSpecies(meTree *T,meNode *v, double **D, double **A) 
  304.      /*the key function of the program addSpeices inserts
  305.        the meNode v to the meTree T.  It uses testEdge to see what the
  306.        weight would be if v split a particular edge.  Weights are assigned by
  307.        OLS formula*/
  308. {
  309.   meTree *T_e;
  310.   meEdge *e; /*loop variable*/
  311.   meEdge *e_min; /*points to best meEdge seen thus far*/
  312.   double w_min = 0.0;   /*used to keep track of meTree weights*/
  314.   if (verbose)
  315.     printf("Adding %s.n",v->label);
  316.  
  317.   /*initialize variables as necessary*/
  318.   /*CASE 1: T is empty, v is the first node*/
  319.   if (NULL == T)  /*create a meTree with v as only vertex, no edges*/
  320.     {
  321.       T_e = newTree();
  322.       T_e->root = v;  
  323.       /*note that we are rooting T arbitrarily at a leaf.  
  324. T->root is not the phylogenetic root*/
  325.       v->index = 0;
  326.       T_e->size = 1;
  327.       return(T_e);      
  328.     }
  329.   /*CASE 2: T is a single-vertex tree*/
  330.   if (1 == T->size)
  331. {
  332.   v->index = 1;
  333.   e = makeEdge("",T->root,v,0.0);
  334.   sprintf(e->label,"E1");
  335.   e->topsize = 1;
  336.   e->bottomsize = 1;
  337.   A[v->index][v->index] = D[v->index2][T->root->index2];
  338.   T->root->leftEdge = v->parentEdge = e;
  339.   T->size = 2;
  340.   return(T); 
  341. }
  342.   /*CASE 3: T has at least two nodes and an edge.  Insert new node
  343.     by breaking one of the edges*/
  344.   
  345.   v->index = T->size;
  346.   /*if (!(T->size % 100))
  347.     printf("T->size is %dn",T->size);*/
  348.   GMEcalcNewvAverages(T,v,D,A);
  349.   /*calcNewvAverges will assign values to all the meEdge averages of T which
  350.     include the meNode v.  Will do so using pre-existing averages in T and
  351.     information from A,D*/
  352.   e_min = T->root->leftEdge;  
  353.   e = e_min->head->leftEdge;  
  354.   while (NULL != e)
  355.     {
  356.       testEdge(e,v,A); 
  357.       /*testEdge tests weight of meTree if loop variable 
  358. e is the meEdge split, places this weight in e->totalweight field */
  359.       if (e->totalweight < w_min)
  360. {
  361.   e_min = e;
  362.   w_min = e->totalweight;
  363. }
  364.       e = topFirstTraverse(T,e);
  365.     }
  366.   /*e_min now points at the meEdge we want to split*/
  367.   GMEsplitEdge(T,v,e_min,A);
  368.   return(T);
  369. }
  371. void updateSubTreeAverages(double **A, meEdge *e, meNode *v, int direction);
  373. /*the ME version of updateAveragesMatrix does not update the entire matrix
  374.   A, but updates A[v->index][w->index] whenever this represents an average
  375.   of 1-distant or 2-distant subtrees*/
  377. void GMEupdateAveragesMatrix(double **A, meEdge *e, meNode *v, meNode *newNode)
  378. {
  379.   meEdge *sib, *par, *left, *right;
  380.   sib = siblingEdge(e);
  381.   left = e->head->leftEdge;
  382.   right = e->head->rightEdge;
  383.   par = e->tail->parentEdge;
  385.   /*we need to update the matrix A so all 1-distant, 2-distant, and
  386.     3-distant averages are correct*/
  387.   
  388.   /*first, initialize the newNode entries*/
  389.   /*1-distant*/
  390.   A[newNode->index][newNode->index] = 
  391.     (e->bottomsize*A[e->head->index][e->head->index]
  392.      + A[v->index][e->head->index])
  393.     / (e->bottomsize + 1);
  394.   /*1-distant for v*/
  395.   A[v->index][v->index] = 
  396.     (e->bottomsize*A[e->head->index][v->index] 
  397.      + e->topsize*A[v->index][e->head->index])
  398.     / (e->bottomsize + e->topsize);
  400.   /*2-distant for v,newNode*/
  401.   A[v->index][newNode->index] = A[newNode->index][v->index] = 
  402.     A[v->index][e->head->index];
  403.   
  404.   /*second 2-distant for newNode*/
  405.   A[newNode->index][e->tail->index] = A[e->tail->index][newNode->index]
  406.     = (e->bottomsize*A[e->head->index][e->tail->index]
  407.        + A[v->index][e->tail->index])/(e->bottomsize + 1);
  408.   /*third 2-distant for newNode*/
  409.   A[newNode->index][e->head->index] = A[e->head->index][newNode->index]
  410.     = A[e->head->index][e->head->index];
  411.    
  412.   if (NULL != sib)
  413.     {
  414.       /*fourth and last 2-distant for newNode*/
  415.       A[newNode->index][sib->head->index] =
  416. A[sib->head->index][newNode->index] = 
  417. (e->bottomsize*A[sib->head->index][e->head->index]
  418.  + A[sib->head->index][v->index]) / (e->bottomsize + 1);
  419.       updateSubTreeAverages(A,sib,v,SKEW); /*updates sib and below*/
  420.     }
  421.   if (NULL != par)
  422.     {
  423.       if (e->tail->leftEdge == e)
  424. updateSubTreeAverages(A,par,v,LEFT); /*updates par and above*/
  425.       else
  426. updateSubTreeAverages(A,par,v,RIGHT);
  427.     }
  428.   if (NULL != left)
  429.     updateSubTreeAverages(A,left,v,UP); /*updates left and below*/
  430.   if (NULL != right)
  431.     updateSubTreeAverages(A,right,v,UP); /*updates right and below*/  
  433.   /*1-dist for e->head*/
  434.   A[e->head->index][e->head->index] = 
  435.     (e->topsize*A[e->head->index][e->head->index]
  436.      + A[e->head->index][v->index]) / (e->topsize+1);
  437.   /*2-dist for e->head (v,newNode,left,right)
  438.     taken care of elsewhere*/
  439.   /*3-dist with e->head either taken care of elsewhere (below)
  440.     or unchanged (sib,e->tail)*/
  441.   
  442.   /*symmetrize the matrix (at least for distant-2 subtrees) */
  443.   A[v->index][e->head->index] = A[e->head->index][v->index];
  444.   /*and distant-3 subtrees*/
  445.   A[e->tail->index][v->index] = A[v->index][e->tail->index];
  446.   if (NULL != left)
  447.     A[v->index][left->head->index] = A[left->head->index][v->index];
  448.   if (NULL != right)
  449.     A[v->index][right->head->index] = A[right->head->index][v->index];
  450.   if (NULL != sib)
  451.     A[v->index][sib->head->index] = A[sib->head->index][v->index];
  453. }      
  454.   
  455. void GMEsplitEdge(meTree *T, meNode *v, meEdge *e, double **A)
  456. {
  457.   char nodelabel[NODE_LABEL_LENGTH];
  458.   char edgelabel[EDGE_LABEL_LENGTH];
  459.   meEdge *newPendantEdge;
  460.   meEdge *newInternalEdge;
  461.   meNode *newNode;
  462.     
  463.   sprintf(nodelabel,"I%d",T->size+1); 
  464.   newNode = makeNewNode(nodelabel,T->size + 1);  
  465.   
  466.   sprintf(edgelabel,"E%d",T->size); 
  467.   newPendantEdge = makeEdge(edgelabel,newNode,v,0.0);   
  468.   
  469.   sprintf(edgelabel,"E%d",T->size+1); 
  470.   newInternalEdge = makeEdge(edgelabel,newNode,e->head,0.0);   
  471.   
  472.   if (verbose)
  473.     printf("Inserting meNode %s on meEdge %s between nodes %s and %s.n",
  474.    v->label,e->label,e->tail->label,e->head->label); 
  475.   /*update the matrix of average distances*/
  476.   /*also updates the bottomsize, topsize fields*/
  477.   
  478.   GMEupdateAveragesMatrix(A,e,v,newNode);
  480.   newNode->parentEdge = e;
  481.   e->head->parentEdge = newInternalEdge;
  482.   v->parentEdge = newPendantEdge;
  483.   e->head = newNode;
  484.   
  485.   T->size = T->size + 2;
  487.   if (e->tail->leftEdge == e) 
  488.     {
  489.       newNode->leftEdge = newInternalEdge;
  490.       newNode->rightEdge = newPendantEdge;
  491.     }
  492.   else
  493.     {
  494.       newNode->leftEdge = newInternalEdge;
  495.       newNode->rightEdge = newPendantEdge;
  496.     }
  497.   
  498.   /*assign proper topsize, bottomsize values to the two new Edges*/
  499.   
  500.   newPendantEdge->bottomsize = 1; 
  501.   newPendantEdge->topsize = e->bottomsize + e->topsize;
  502.   
  503.   newInternalEdge->bottomsize = e->bottomsize;
  504.   newInternalEdge->topsize = e->topsize;  /*off by one, but we adjust
  505.     that below*/
  506.   
  507.   /*and increment these fields for all other edges*/
  508.   updateSizes(newInternalEdge,UP);
  509.   updateSizes(e,DOWN);
  510. }
  512. void updateSubTreeAverages(double **A, meEdge *e, meNode *v, int direction)
  513.      /*the monster function of this program*/
  514. {
  515.   meEdge *sib, *left, *right, *par;
  516.   left = e->head->leftEdge;
  517.   right = e->head->rightEdge;
  518.   sib = siblingEdge(e);
  519.   par = e->tail->parentEdge;
  520.   switch(direction)
  521.     {
  522.       /*want to preserve correctness of 
  523. all 1-distant, 2-distant, and 3-distant averages*/
  524.       /*1-distant updated at meEdge splitting the two trees*/
  525.       /*2-distant updated:
  526. (left->head,right->head) and (head,tail) updated at
  527. a given edge.  Note, NOT updating (head,sib->head)!
  528. (That would lead to multiple updating).*/
  529.       /*3-distant updated: at meEdge in center of quartet*/
  530.     case UP: /*point of insertion is above e*/
  531.       /*1-distant average of nodes below e to 
  532.        nodes above e*/
  533.       A[e->head->index][e->head->index] = 
  534. (e->topsize*A[e->head->index][e->head->index] + 
  535.  A[e->head->index][v->index])/(e->topsize + 1);      
  536.       /*2-distant average of nodes below e to 
  537. nodes above parent of e*/
  538.       A[e->head->index][par->head->index] = 
  539. A[par->head->index][e->head->index] = 
  540. (par->topsize*A[par->head->index][e->head->index]
  541.  + A[e->head->index][v->index]) / (par->topsize + 1);
  542.       /*must do both 3-distant averages involving par*/
  543.       if (NULL != left)
  544. {
  545.   updateSubTreeAverages(A, left, v, UP); /*and recursive call*/
  546.   /*3-distant average*/
  547.   A[par->head->index][left->head->index]
  548.     = A[left->head->index][par->head->index]
  549.     = (par->topsize*A[par->head->index][left->head->index]
  550.        + A[left->head->index][v->index]) / (par->topsize + 1);
  551. }
  552.       if (NULL != right)
  553. {
  554.   updateSubTreeAverages(A, right, v, UP);
  555.   A[par->head->index][right->head->index]
  556.     = A[right->head->index][par->head->index]
  557.     = (par->topsize*A[par->head->index][right->head->index]
  558.        + A[right->head->index][v->index]) / (par->topsize + 1);
  559. }
  560.       break;
  561.     case SKEW: /*point of insertion is skew to e*/
  562.       /*1-distant average of nodes below e to 
  563. nodes above e*/
  564.       A[e->head->index][e->head->index] = 
  565. (e->topsize*A[e->head->index][e->head->index] + 
  566.  A[e->head->index][v->index])/(e->topsize + 1);      
  567.       /*no 2-distant averages to update in this case*/
  568.       /*updating 3-distant averages involving sib*/
  569.       if (NULL != left)
  570. {
  571.   updateSubTreeAverages(A, left, v, UP);
  572.   A[sib->head->index][left->head->index]
  573.     = A[left->head->index][sib->head->index]
  574.     = (sib->bottomsize*A[sib->head->index][left->head->index]
  575.        + A[left->head->index][v->index]) / (sib->bottomsize + 1);
  576. }
  577.       if (NULL != right)
  578. {
  579.   updateSubTreeAverages(A, right, v, UP);
  580.   A[sib->head->index][right->head->index]
  581.     = A[right->head->index][sib->head->index]
  582.     = (sib->bottomsize*A[par->head->index][right->head->index]
  583.        + A[right->head->index][v->index]) / (sib->bottomsize + 1);
  584. }
  585.       break;
  588.     case LEFT: /*point of insertion is below the meEdge left*/
  589.       /*1-distant average*/
  590.       A[e->head->index][e->head->index] = 
  591. (e->bottomsize*A[e->head->index][e->head->index] + 
  592.  A[v->index][e->head->index])/(e->bottomsize + 1);        
  593.       /*2-distant averages*/
  594.       A[e->head->index][e->tail->index] = 
  595. A[e->tail->index][e->head->index] = 
  596. (e->bottomsize*A[e->head->index][e->tail->index] + 
  597.  A[v->index][e->tail->index])/(e->bottomsize + 1);  
  598.       A[left->head->index][right->head->index] = 
  599. A[right->head->index][left->head->index] = 
  600. (left->bottomsize*A[right->head->index][left->head->index]
  601.  + A[right->head->index][v->index]) / (left->bottomsize+1);
  602.       /*3-distant avereages involving left*/
  603.       if (NULL != sib)
  604. {
  605.   updateSubTreeAverages(A, sib, v, SKEW);
  606.   A[left->head->index][sib->head->index]
  607.     = A[sib->head->index][left->head->index]
  608.     = (left->bottomsize*A[left->head->index][sib->head->index]
  609.        + A[sib->head->index][v->index]) / (left->bottomsize + 1);
  610. }
  611.       if (NULL != par)
  612. {
  613.   if (e->tail->leftEdge == e)
  614.     updateSubTreeAverages(A, par, v, LEFT);
  615.   else
  616.     updateSubTreeAverages(A, par, v, RIGHT);
  617.   A[left->head->index][par->head->index]
  618.     = A[par->head->index][left->head->index]
  619.     = (left->bottomsize*A[left->head->index][par->head->index]
  620.        + A[v->index][par->head->index]) / (left->bottomsize + 1);
  621. }
  622.       break;    
  623.     case RIGHT: /*point of insertion is below the meEdge right*/
  624.       /*1-distant average*/
  625.       A[e->head->index][e->head->index] = 
  626. (e->bottomsize*A[e->head->index][e->head->index] + 
  627.  A[v->index][e->head->index])/(e->bottomsize + 1);        
  628.       /*2-distant averages*/
  629.       A[e->head->index][e->tail->index] = 
  630. A[e->tail->index][e->head->index] = 
  631. (e->bottomsize*A[e->head->index][e->tail->index] + 
  632.  A[v->index][e->tail->index])/(e->bottomsize + 1);  
  633.       A[left->head->index][right->head->index] = 
  634. A[right->head->index][left->head->index] = 
  635. (right->bottomsize*A[right->head->index][left->head->index]
  636.  + A[left->head->index][v->index]) / (right->bottomsize+1);
  637.       /*3-distant avereages involving right*/
  638.       if (NULL != sib)
  639. {
  640.   updateSubTreeAverages(A, sib, v, SKEW);
  641.   A[right->head->index][sib->head->index]
  642.     = A[sib->head->index][right->head->index]
  643.     = (right->bottomsize*A[right->head->index][sib->head->index]
  644.        + A[sib->head->index][v->index]) / (right->bottomsize + 1);
  645. }
  646.       if (NULL != par)
  647. {
  648.   if (e->tail->leftEdge == e)
  649.     updateSubTreeAverages(A, par, v, LEFT);
  650.   else
  651.     updateSubTreeAverages(A, par, v, RIGHT);
  652.   A[right->head->index][par->head->index]
  653.     = A[par->head->index][right->head->index]
  654.     = (right->bottomsize*A[right->head->index][par->head->index]
  655.        + A[v->index][par->head->index]) / (right->bottomsize + 1);
  656. }
  658.       break;
  659.     }
  660. }
  662. void assignBottomsize(meEdge *e)
  663. {
  664.   if (leaf(e->head))
  665.     e->bottomsize = 1;
  666.   else
  667.     {
  668.       assignBottomsize(e->head->leftEdge);
  669.       assignBottomsize(e->head->rightEdge);
  670.       e->bottomsize = e->head->leftEdge->bottomsize
  671. + e->head->rightEdge->bottomsize;
  672.     }
  673. }
  675. void assignTopsize(meEdge *e, int numLeaves)
  676. {
  677.   if (NULL != e)
  678.     {
  679.       e->topsize = numLeaves - e->bottomsize;
  680.       assignTopsize(e->head->leftEdge,numLeaves);
  681.       assignTopsize(e->head->rightEdge,numLeaves);
  682.     }
  683. }
  685. void assignAllSizeFields(meTree *T)
  686. {
  687.   assignBottomsize(T->root->leftEdge);
  688.   assignTopsize(T->root->leftEdge,T->size/2 + 1);
  689. }
  692. END_SCOPE(fastme)
  693. END_NCBI_SCOPE
  695. /*
  696.  * ===========================================================================
  697.  * $Log: gme.cpp,v $
  698.  * Revision 1000.1  2004/06/01 18:09:56  gouriano
  699.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
  700.  *
  701.  * Revision 1.2  2004/05/21 21:41:03  gorelenk
  702.  * Added PCH ncbi_pch.hpp
  703.  *
  704.  * Revision 1.1  2004/02/10 15:16:02  jcherry
  705.  * Initial version
  706.  *
  707.  * ===========================================================================
  708.  */