NNI.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:13k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: NNI.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:09:44  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: NNI.cpp,v 1000.1 2004/06/01 18:09:44 gouriano Exp $
  10. * ===========================================================================
  11. *
  12. *                            PUBLIC DOMAIN NOTICE
  13. *               National Center for Biotechnology Information
  14. *
  15. *  This software/database is a "United States Government Work" under the
  16. *  terms of the United States Copyright Act.  It was written as part of
  17. *  the author's official duties as a United States Government employee and
  18. *  thus cannot be copyrighted.  This software/database is freely available
  19. *  to the public for use. The National Library of Medicine and the U.S.
  20. *  Government have not placed any restriction on its use or reproduction.
  21. *
  22. *  Although all reasonable efforts have been taken to ensure the accuracy
  23. *  and reliability of the software and data, the NLM and the U.S.
  24. *  Government do not and cannot warrant the performance or results that
  25. *  may be obtained by using this software or data. The NLM and the U.S.
  26. *  Government disclaim all warranties, express or implied, including
  27. *  warranties of performance, merchantability or fitness for any particular
  28. *  purpose.
  29. *
  30. *  Please cite the author in any work or product based on this material.
  31. *
  32. * ===========================================================================
  33. *
  34. * Author:  Richard Desper
  35. *
  36. * File Description:  NNI.cpp
  37. *
  38. *    A part of the Miminum Evolution algorithm
  39. *
  40. */
  41. #include <ncbi_pch.hpp>
  42. #include <stdio.h>
  43. #include <stdlib.h>
  44. #include <math.h>
  45. #include "graph.h"
  46. #include "fastme.h"
  47. BEGIN_NCBI_SCOPE
  48. BEGIN_SCOPE(fastme)
  49. boolean leaf(meNode *v);
  50. meEdge *siblingEdge(meEdge *e);
  51. meEdge *depthFirstTraverse(meTree *T, meEdge *e);
  52. meEdge *findBottomLeft(meEdge *e);
  53. meEdge *topFirstTraverse(meTree *T, meEdge *e);
  54. meEdge *moveUpRight(meEdge *e);
  55. double wf(double lambda, double D_LR, double D_LU, double D_LD, 
  56.   double D_RU, double D_RD, double D_DU);
  57. /*NNI functions for unweighted OLS topological switches*/
  58. /*fillTableUp fills all the entries in D associated with
  59.   e->head,f->head and those edges g->head above e->head*/
  60. void fillTableUp(meEdge *e, meEdge *f, double **A, double **D, meTree *T)
  61. {
  62.   meEdge *g,*h;
  63.   if (T->root == f->tail)
  64.     {
  65.       if (leaf(e->head))
  66. A[e->head->index][f->head->index] = 
  67.   A[f->head->index][e->head->index] = 
  68.   D[e->head->index2][f->tail->index2];
  69.       else
  70. {
  71.   g = e->head->leftEdge;
  72.   h = e->head->rightEdge;
  73.   A[e->head->index][f->head->index] = 
  74.     A[f->head->index][e->head->index] =  
  75.     (g->bottomsize*A[f->head->index][g->head->index]
  76.      + h->bottomsize*A[f->head->index][h->head->index])
  77.     /e->bottomsize;  
  78. }
  79.     }
  80.   else 
  81.     {
  82.       g = f->tail->parentEdge;
  83.       fillTableUp(e,g,A,D,T); /*recursive call*/
  84.       h = siblingEdge(f);
  85.       A[e->head->index][f->head->index] = 
  86. A[f->head->index][e->head->index] =  
  87. (g->topsize*A[e->head->index][g->head->index]
  88.  + h->bottomsize*A[e->head->index][h->head->index])/f->topsize;    
  89.     }
  90. }
  91. void makeOLSAveragesTable(meTree *T, double **D, double **A);
  92. double **buildAveragesTable(meTree *T, double **D)
  93. {
  94.   int i,j;
  95.   double **A;
  96.   A = (double **) malloc(T->size*sizeof(double *));
  97.   for(i = 0; i < T->size;i++)
  98.     {
  99.       A[i] = (double *) malloc(T->size*sizeof(double));
  100.       for(j=0;j<T->size;j++)
  101. A[i][j] = 0.0;
  102.     }
  103.   makeOLSAveragesTable(T,D,A);
  104.   return(A);
  105. }
  106. double wf2(double lambda, double D_AD, double D_BC, double D_AC, double D_BD,
  107.    double D_AB, double D_CD)
  108. {
  109.   double weight;
  110.   weight = 0.5*(lambda*(D_AC + D_BD) + (1 - lambda)*(D_AD + D_BC)
  111. + (D_AB + D_CD));
  112.   return(weight);
  113. }
  114. int NNIEdgeTest(meEdge *e, meTree *T, double **A, double *weight)
  115. {
  116.   int a,b,c,d;
  117.   meEdge *f;
  118.   double *lambda;
  119.   double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
  120.   double w1,w2,w0;
  121.   
  122.   if ((leaf(e->tail)) || (leaf(e->head)))
  123.     return(NONE);
  124.   lambda = (double *)malloc(3*sizeof(double));
  125.   a = e->tail->parentEdge->topsize;
  126.   f = siblingEdge(e);
  127.   b = f->bottomsize;  
  128.   c = e->head->leftEdge->bottomsize;
  129.   d = e->head->rightEdge->bottomsize;
  130.   lambda[0] = ((double) b*c + a*d)/((a + b)*(c+d));
  131.   lambda[1] = ((double) b*c + a*d)/((a + c)*(b+d));    
  132.   lambda[2] = ((double) c*d + a*b)/((a + d)*(b+c));
  133.   
  134.   D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
  135.   D_LU = A[e->head->leftEdge->head->index][e->tail->index];
  136.   D_LD = A[e->head->leftEdge->head->index][f->head->index];
  137.   D_RU = A[e->head->rightEdge->head->index][e->tail->index];
  138.   D_RD = A[e->head->rightEdge->head->index][f->head->index];
  139.   D_DU = A[e->tail->index][f->head->index];
  140.   w0 = wf2(lambda[0],D_RU,D_LD,D_LU,D_RD,D_DU,D_LR);
  141.   w1 = wf2(lambda[1],D_RU,D_LD,D_DU,D_LR,D_LU,D_RD);
  142.   w2 = wf2(lambda[2],D_DU,D_LR,D_LU,D_RD,D_RU,D_LD);
  143.   free(lambda);
  144.   if (w0 <= w1)
  145.     {
  146.       if (w0 <= w2) /*w0 <= w1,w2*/
  147. {
  148.   *weight = 0.0;
  149.   return(NONE);
  150. }
  151.       else /*w2 < w0 <= w1 */
  152. {
  153.   *weight = w2 - w0;
  154.   if (verbose)
  155.     {
  156.       printf("Swap across %s. ",e->label);
  157.       printf("Weight dropping by %lf.n",w0 - w2);
  158.       printf("New weight should be %lf.n",T->weight + w2 - w0);
  159.     }
  160.   return(RIGHT);
  161. }
  162.     }
  163.   else if (w2 <= w1) /*w2 <= w1 < w0*/
  164.     {
  165.       *weight = w2 - w0;
  166.       if (verbose)
  167. {
  168.   printf("Swap across %s. ",e->label);
  169.   printf("Weight dropping by %lf.n",w0 - w2);
  170.   printf("New weight should be %lf.n",T->weight + w2 - w0);
  171. }
  172.       return(RIGHT);
  173.     }
  174.   else /*w1 < w2, w0*/
  175.     {
  176.       *weight = w1 - w0;
  177.       if (verbose)
  178. {
  179.   printf("Swap across %s. ",e->label);
  180.   printf("Weight dropping by %lf.n",w0 - w1);
  181.   printf("New weight should be %lf.n",T->weight + w1 - w0);
  182. }
  183.       return(LEFT);
  184.     }
  185. }
  186.  
  187. int *initPerm(int size);
  188. void NNIupdateAverages(double **A, meEdge *e, meEdge *par, meEdge *skew, 
  189.        meEdge *swap, meEdge *fixed, meTree *T)
  190. {
  191.   meNode *v;
  192.   meEdge *elooper;
  193.   v = e->head;
  194.   /*first, v*/
  195.   A[e->head->index][e->head->index] =  
  196.     (swap->bottomsize* 
  197.      ((skew->bottomsize*A[skew->head->index][swap->head->index]
  198.        + fixed->bottomsize*A[fixed->head->index][swap->head->index]) 
  199.       / e->bottomsize) +
  200.      par->topsize*
  201.      ((skew->bottomsize*A[skew->head->index][par->head->index]
  202.        + fixed->bottomsize*A[fixed->head->index][par->head->index]) 
  203.       / e->bottomsize)
  204.      ) / e->topsize; 
  205.   
  206.   elooper = findBottomLeft(e); /*next, we loop over all the edges 
  207.  which are below e*/
  208.   while (e != elooper)  
  209.     {
  210.       A[e->head->index][elooper->head->index] = 
  211. A[elooper->head->index][v->index] 
  212. = (swap->bottomsize*A[elooper->head->index][swap->head->index] +
  213.    par->topsize*A[elooper->head->index][par->head->index]) 
  214. / e->topsize;
  215.       elooper = depthFirstTraverse(T,elooper);
  216.     }
  217.   elooper = findBottomLeft(swap); /*next we loop over all the edges below and
  218.     including swap*/  
  219.   while (swap != elooper)
  220.   {
  221.     A[e->head->index][elooper->head->index] = 
  222.       A[elooper->head->index][e->head->index]
  223.       = (skew->bottomsize * A[elooper->head->index][skew->head->index] + 
  224.  fixed->bottomsize*A[elooper->head->index][fixed->head->index]) 
  225.       / e->bottomsize;
  226.     elooper = depthFirstTraverse(T,elooper);
  227.   }
  228.   /*now elooper = skew */
  229.   A[e->head->index][elooper->head->index] = 
  230.     A[elooper->head->index][e->head->index]
  231.     = (skew->bottomsize * A[elooper->head->index][skew->head->index] + 
  232.        fixed->bottomsize* A[elooper->head->index][fixed->head->index])
  233.     / e->bottomsize;
  234.   
  235.   /*finally, we loop over all the edges in the meTree 
  236.     on the far side of parEdge*/ 
  237.   elooper = T->root->leftEdge;
  238.   while ((elooper != swap) && (elooper != e)) /*start a top-first traversal*/
  239.     {
  240.       A[e->head->index][elooper->head->index] = 
  241. A[elooper->head->index][e->head->index]
  242. = (skew->bottomsize * A[elooper->head->index][skew->head->index] 
  243.    + fixed->bottomsize* A[elooper->head->index][fixed->head->index]) 
  244. / e->bottomsize;
  245.       elooper = topFirstTraverse(T,elooper);
  246.     }
  247.   /*At this point, elooper = par.
  248.     We finish the top-first traversal, excluding the submeTree below par*/
  249.   elooper = moveUpRight(par);
  250.   while (NULL != elooper)
  251.     {
  252.       A[e->head->index][elooper->head->index] 
  253. = A[elooper->head->index][e->head->index]
  254. = (skew->bottomsize * A[elooper->head->index][skew->head->index] + 
  255.    fixed->bottomsize* A[elooper->head->index][fixed->head->index]) 
  256. / e->bottomsize;
  257.       elooper = topFirstTraverse(T,elooper);
  258.     }
  259.   
  260. }
  261. void NNItopSwitch(meTree *T, meEdge *e, int direction, double **A)
  262. {
  263.   meEdge *par, *fixed;
  264.   meEdge *skew, *swap;
  265.   
  266.   if (verbose)
  267.     printf("Branch swap across meEdge %s.n",e->label);
  268.   if (LEFT == direction)
  269.     swap = e->head->leftEdge;
  270.   else
  271.     swap = e->head->rightEdge;
  272.   skew = siblingEdge(e);
  273.   fixed = siblingEdge(swap);
  274.   par = e->tail->parentEdge;
  275.   
  276.   if (verbose)
  277.     {
  278.       printf("Branch swap: switching edges %s and %s.n",skew->label,swap->label);
  279.     }
  280.   /*perform topological switch by changing f from (u,b) to (v,b)
  281.     and g from (v,c) to (u,c), necessitatates also changing parent fields*/
  282.   
  283.   swap->tail = e->tail;
  284.   skew->tail = e->head;
  285.   
  286.   if (LEFT == direction)
  287.     e->head->leftEdge = skew;
  288.   else
  289.     e->head->rightEdge = skew;
  290.   if (skew == e->tail->rightEdge)
  291.     e->tail->rightEdge = swap;
  292.   else
  293.     e->tail->leftEdge = swap;
  294.   /*both topsize and bottomsize change for e, but nowhere else*/
  295.   e->topsize = par->topsize + swap->bottomsize;
  296.   e->bottomsize = fixed->bottomsize + skew->bottomsize;
  297.   NNIupdateAverages(A, e, par, skew, swap, fixed,T);
  298. } /*end NNItopSwitch*/
  299. void reHeapElement(int *p, int *q, double *v, int length, int i);
  300. void pushHeap(int *p, int *q, double *v, int length, int i);
  301. void popHeap(int *p, int *q, double *v, int length, int i);
  302. void NNIRetestEdge(int *p, int *q, meEdge *e,meTree *T, double **avgDistArray, 
  303. double *weights, int *location, int *possibleSwaps)
  304. {
  305.   int tloc;
  306.   tloc = location[e->head->index+1];
  307.   location[e->head->index+1] = 
  308.     NNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
  309.   if (NONE == location[e->head->index+1])
  310.     {
  311.       if (NONE != tloc)
  312. popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);   
  313.     }
  314.   else
  315.     {
  316.       if (NONE == tloc)
  317. pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
  318.       else
  319. reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
  320.     }
  321. }
  322. void permInverse(int *p, int *q, int length);
  323. int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
  324. void NNI(meTree *T, double **avgDistArray, int *count)
  325. {
  326.   meEdge *e, *centerEdge;
  327.   meEdge **edgeArray;
  328.   int *location;
  329.   int *p,*q;
  330.   int i;
  331.   int possibleSwaps;
  332.   double *weights;
  333.   p = initPerm(T->size+1);
  334.   q = initPerm(T->size+1);
  335.   edgeArray = (meEdge **) malloc((T->size+1)*sizeof(double));
  336.   weights = (double *) malloc((T->size+1)*sizeof(double));
  337.   location = (int *) malloc((T->size+1)*sizeof(int));
  338.   for (i=0;i<T->size+1;i++)
  339.     {
  340.       weights[i] = 0.0;
  341.       location[i] = NONE;
  342.     }
  343.   e = findBottomLeft(T->root->leftEdge); 
  344.   /* *count = 0; */
  345.   while (NULL != e)
  346.     {
  347.       edgeArray[e->head->index+1] = e;
  348.       location[e->head->index+1] = 
  349. NNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
  350.       e = depthFirstTraverse(T,e);
  351.     } 
  352.   possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
  353.   permInverse(p,q,T->size+1);
  354.   /*we put the negative values of weights into a heap, indexed by p
  355.     with the minimum value pointed to by p[1]*/
  356.   /*p[i] is index (in edgeArray) of meEdge with i-th position 
  357.     in the heap, q[j] is the position of meEdge j in the heap */
  358.   while (weights[p[1]] < 0)
  359.     {
  360.       centerEdge = edgeArray[p[1]];
  361.       (*count)++;
  362.       T->weight = T->weight + weights[p[1]];
  363.       NNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
  364.       location[p[1]] = NONE;
  365.       weights[p[1]] = 0.0;  /*after the NNI, this meEdge is in optimal
  366.       configuration*/
  367.       popHeap(p,q,weights,possibleSwaps--,1);
  368.       /*but we must retest the other four edges*/
  369.       e = centerEdge->head->leftEdge;
  370.       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
  371.       e = centerEdge->head->rightEdge;
  372.       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);     
  373.       e = siblingEdge(centerEdge);
  374.       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
  375.       e = centerEdge->tail->parentEdge;
  376.       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
  377.     }
  378.   free(p);
  379.   free(q);
  380.   free(location);
  381.   free(edgeArray);
  382. }
  383. void NNIwithoutMatrix(meTree *T, double **D, int *count)
  384. {
  385.   double **avgDistArray;
  386.   avgDistArray = buildAveragesTable(T,D);
  387.   NNI(T,avgDistArray,count);
  388. }
  389. void NNIWithPartialMatrix(meTree *T,double **D,double **A,int *count)
  390. {
  391.   makeOLSAveragesTable(T,D,A);
  392.   NNI(T,A,count);
  393. }
  394. END_SCOPE(fastme)
  395. END_NCBI_SCOPE
  396. /*
  397.  * ===========================================================================
  398.  * $Log: NNI.cpp,v $
  399.  * Revision 1000.1  2004/06/01 18:09:44  gouriano
  400.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
  401.  *
  402.  * Revision 1.2  2004/05/21 21:41:03  gorelenk
  403.  * Added PCH ncbi_pch.hpp
  404.  *
  405.  * Revision 1.1  2004/02/10 15:16:00  jcherry
  406.  * Initial version
  407.  *
  408.  * ===========================================================================
  409.  */