NNI.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:13k
- /*
- * ===========================================================================
- * PRODUCTION $Log: NNI.cpp,v $
- * PRODUCTION Revision 1000.1 2004/06/01 18:09:44 gouriano
- * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
- * PRODUCTION
- * ===========================================================================
- */
- /* $Id: NNI.cpp,v 1000.1 2004/06/01 18:09:44 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: NNI.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)
- boolean leaf(meNode *v);
- meEdge *siblingEdge(meEdge *e);
- meEdge *depthFirstTraverse(meTree *T, meEdge *e);
- meEdge *findBottomLeft(meEdge *e);
- meEdge *topFirstTraverse(meTree *T, meEdge *e);
- meEdge *moveUpRight(meEdge *e);
- double wf(double lambda, double D_LR, double D_LU, double D_LD,
- double D_RU, double D_RD, double D_DU);
- /*NNI functions for unweighted OLS topological switches*/
- /*fillTableUp fills all the entries in D associated with
- e->head,f->head and those edges g->head above e->head*/
- void fillTableUp(meEdge *e, meEdge *f, double **A, double **D, meTree *T)
- {
- meEdge *g,*h;
- if (T->root == f->tail)
- {
- if (leaf(e->head))
- A[e->head->index][f->head->index] =
- A[f->head->index][e->head->index] =
- D[e->head->index2][f->tail->index2];
- else
- {
- g = e->head->leftEdge;
- h = e->head->rightEdge;
- 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;
- }
- }
- else
- {
- g = f->tail->parentEdge;
- fillTableUp(e,g,A,D,T); /*recursive call*/
- h = siblingEdge(f);
- A[e->head->index][f->head->index] =
- A[f->head->index][e->head->index] =
- (g->topsize*A[e->head->index][g->head->index]
- + h->bottomsize*A[e->head->index][h->head->index])/f->topsize;
- }
- }
- void makeOLSAveragesTable(meTree *T, double **D, double **A);
- double **buildAveragesTable(meTree *T, double **D)
- {
- int i,j;
- double **A;
- A = (double **) malloc(T->size*sizeof(double *));
- for(i = 0; i < T->size;i++)
- {
- A[i] = (double *) malloc(T->size*sizeof(double));
- for(j=0;j<T->size;j++)
- A[i][j] = 0.0;
- }
- makeOLSAveragesTable(T,D,A);
- return(A);
- }
- double wf2(double lambda, double D_AD, double D_BC, double D_AC, double D_BD,
- double D_AB, double D_CD)
- {
- double weight;
- weight = 0.5*(lambda*(D_AC + D_BD) + (1 - lambda)*(D_AD + D_BC)
- + (D_AB + D_CD));
- return(weight);
- }
- int NNIEdgeTest(meEdge *e, meTree *T, double **A, double *weight)
- {
- int a,b,c,d;
- meEdge *f;
- double *lambda;
- double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
- double w1,w2,w0;
-
- if ((leaf(e->tail)) || (leaf(e->head)))
- return(NONE);
- lambda = (double *)malloc(3*sizeof(double));
- a = e->tail->parentEdge->topsize;
- f = siblingEdge(e);
- b = f->bottomsize;
- c = e->head->leftEdge->bottomsize;
- d = e->head->rightEdge->bottomsize;
- lambda[0] = ((double) b*c + a*d)/((a + b)*(c+d));
- lambda[1] = ((double) b*c + a*d)/((a + c)*(b+d));
- lambda[2] = ((double) c*d + a*b)/((a + d)*(b+c));
-
- D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
- D_LU = A[e->head->leftEdge->head->index][e->tail->index];
- D_LD = A[e->head->leftEdge->head->index][f->head->index];
- D_RU = A[e->head->rightEdge->head->index][e->tail->index];
- D_RD = A[e->head->rightEdge->head->index][f->head->index];
- D_DU = A[e->tail->index][f->head->index];
- w0 = wf2(lambda[0],D_RU,D_LD,D_LU,D_RD,D_DU,D_LR);
- w1 = wf2(lambda[1],D_RU,D_LD,D_DU,D_LR,D_LU,D_RD);
- w2 = wf2(lambda[2],D_DU,D_LR,D_LU,D_RD,D_RU,D_LD);
- free(lambda);
- if (w0 <= w1)
- {
- if (w0 <= w2) /*w0 <= w1,w2*/
- {
- *weight = 0.0;
- return(NONE);
- }
- else /*w2 < w0 <= w1 */
- {
- *weight = w2 - w0;
- if (verbose)
- {
- printf("Swap across %s. ",e->label);
- printf("Weight dropping by %lf.n",w0 - w2);
- printf("New weight should be %lf.n",T->weight + w2 - w0);
- }
- return(RIGHT);
- }
- }
- else if (w2 <= w1) /*w2 <= w1 < w0*/
- {
- *weight = w2 - w0;
- if (verbose)
- {
- printf("Swap across %s. ",e->label);
- printf("Weight dropping by %lf.n",w0 - w2);
- printf("New weight should be %lf.n",T->weight + w2 - w0);
- }
- return(RIGHT);
- }
- else /*w1 < w2, w0*/
- {
- *weight = w1 - w0;
- if (verbose)
- {
- printf("Swap across %s. ",e->label);
- printf("Weight dropping by %lf.n",w0 - w1);
- printf("New weight should be %lf.n",T->weight + w1 - w0);
- }
- return(LEFT);
- }
- }
-
- int *initPerm(int size);
- void NNIupdateAverages(double **A, meEdge *e, meEdge *par, meEdge *skew,
- meEdge *swap, meEdge *fixed, meTree *T)
- {
- meNode *v;
- meEdge *elooper;
- v = e->head;
- /*first, v*/
- A[e->head->index][e->head->index] =
- (swap->bottomsize*
- ((skew->bottomsize*A[skew->head->index][swap->head->index]
- + fixed->bottomsize*A[fixed->head->index][swap->head->index])
- / e->bottomsize) +
- par->topsize*
- ((skew->bottomsize*A[skew->head->index][par->head->index]
- + fixed->bottomsize*A[fixed->head->index][par->head->index])
- / e->bottomsize)
- ) / e->topsize;
-
- elooper = findBottomLeft(e); /*next, we loop over all the edges
- which are below e*/
- while (e != elooper)
- {
- A[e->head->index][elooper->head->index] =
- A[elooper->head->index][v->index]
- = (swap->bottomsize*A[elooper->head->index][swap->head->index] +
- par->topsize*A[elooper->head->index][par->head->index])
- / e->topsize;
- elooper = depthFirstTraverse(T,elooper);
- }
- elooper = findBottomLeft(swap); /*next we loop over all the edges below and
- including swap*/
- while (swap != elooper)
- {
- A[e->head->index][elooper->head->index] =
- A[elooper->head->index][e->head->index]
- = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
- fixed->bottomsize*A[elooper->head->index][fixed->head->index])
- / e->bottomsize;
- elooper = depthFirstTraverse(T,elooper);
- }
- /*now elooper = skew */
- A[e->head->index][elooper->head->index] =
- A[elooper->head->index][e->head->index]
- = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
- fixed->bottomsize* A[elooper->head->index][fixed->head->index])
- / e->bottomsize;
-
- /*finally, we loop over all the edges in the meTree
- on the far side of parEdge*/
- elooper = T->root->leftEdge;
- while ((elooper != swap) && (elooper != e)) /*start a top-first traversal*/
- {
- A[e->head->index][elooper->head->index] =
- A[elooper->head->index][e->head->index]
- = (skew->bottomsize * A[elooper->head->index][skew->head->index]
- + fixed->bottomsize* A[elooper->head->index][fixed->head->index])
- / e->bottomsize;
- elooper = topFirstTraverse(T,elooper);
- }
- /*At this point, elooper = par.
- We finish the top-first traversal, excluding the submeTree below par*/
- elooper = moveUpRight(par);
- while (NULL != elooper)
- {
- A[e->head->index][elooper->head->index]
- = A[elooper->head->index][e->head->index]
- = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
- fixed->bottomsize* A[elooper->head->index][fixed->head->index])
- / e->bottomsize;
- elooper = topFirstTraverse(T,elooper);
- }
-
- }
- void NNItopSwitch(meTree *T, meEdge *e, int direction, double **A)
- {
- meEdge *par, *fixed;
- meEdge *skew, *swap;
-
- if (verbose)
- printf("Branch swap across meEdge %s.n",e->label);
- if (LEFT == direction)
- swap = e->head->leftEdge;
- else
- swap = e->head->rightEdge;
- skew = siblingEdge(e);
- fixed = siblingEdge(swap);
- par = e->tail->parentEdge;
-
- if (verbose)
- {
- printf("Branch swap: switching edges %s and %s.n",skew->label,swap->label);
- }
- /*perform topological switch by changing f from (u,b) to (v,b)
- and g from (v,c) to (u,c), necessitatates also changing parent fields*/
-
- swap->tail = e->tail;
- skew->tail = e->head;
-
- if (LEFT == direction)
- e->head->leftEdge = skew;
- else
- e->head->rightEdge = skew;
- if (skew == e->tail->rightEdge)
- e->tail->rightEdge = swap;
- else
- e->tail->leftEdge = swap;
- /*both topsize and bottomsize change for e, but nowhere else*/
- e->topsize = par->topsize + swap->bottomsize;
- e->bottomsize = fixed->bottomsize + skew->bottomsize;
- NNIupdateAverages(A, e, par, skew, swap, fixed,T);
- } /*end NNItopSwitch*/
- void reHeapElement(int *p, int *q, double *v, int length, int i);
- void pushHeap(int *p, int *q, double *v, int length, int i);
- void popHeap(int *p, int *q, double *v, int length, int i);
- void NNIRetestEdge(int *p, int *q, meEdge *e,meTree *T, double **avgDistArray,
- double *weights, int *location, int *possibleSwaps)
- {
- int tloc;
- tloc = location[e->head->index+1];
- location[e->head->index+1] =
- NNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
- if (NONE == location[e->head->index+1])
- {
- if (NONE != tloc)
- popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
- }
- else
- {
- if (NONE == tloc)
- pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
- else
- reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
- }
- }
- void permInverse(int *p, int *q, int length);
- int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
- void NNI(meTree *T, double **avgDistArray, int *count)
- {
- meEdge *e, *centerEdge;
- meEdge **edgeArray;
- int *location;
- int *p,*q;
- int i;
- int possibleSwaps;
- double *weights;
- p = initPerm(T->size+1);
- q = initPerm(T->size+1);
- edgeArray = (meEdge **) malloc((T->size+1)*sizeof(double));
- weights = (double *) malloc((T->size+1)*sizeof(double));
- location = (int *) malloc((T->size+1)*sizeof(int));
- for (i=0;i<T->size+1;i++)
- {
- weights[i] = 0.0;
- location[i] = NONE;
- }
- e = findBottomLeft(T->root->leftEdge);
- /* *count = 0; */
- while (NULL != e)
- {
- edgeArray[e->head->index+1] = e;
- location[e->head->index+1] =
- NNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
- e = depthFirstTraverse(T,e);
- }
- possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
- permInverse(p,q,T->size+1);
- /*we put the negative values of weights into a heap, indexed by p
- with the minimum value pointed to by p[1]*/
- /*p[i] is index (in edgeArray) of meEdge with i-th position
- in the heap, q[j] is the position of meEdge j in the heap */
- while (weights[p[1]] < 0)
- {
- centerEdge = edgeArray[p[1]];
- (*count)++;
- T->weight = T->weight + weights[p[1]];
- NNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
- location[p[1]] = NONE;
- weights[p[1]] = 0.0; /*after the NNI, this meEdge is in optimal
- configuration*/
- popHeap(p,q,weights,possibleSwaps--,1);
- /*but we must retest the other four edges*/
- e = centerEdge->head->leftEdge;
- NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
- e = centerEdge->head->rightEdge;
- NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
- e = siblingEdge(centerEdge);
- NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
- e = centerEdge->tail->parentEdge;
- NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
- }
- free(p);
- free(q);
- free(location);
- free(edgeArray);
- }
- void NNIwithoutMatrix(meTree *T, double **D, int *count)
- {
- double **avgDistArray;
- avgDistArray = buildAveragesTable(T,D);
- NNI(T,avgDistArray,count);
- }
- void NNIWithPartialMatrix(meTree *T,double **D,double **A,int *count)
- {
- makeOLSAveragesTable(T,D,A);
- NNI(T,A,count);
- }
- END_SCOPE(fastme)
- END_NCBI_SCOPE
- /*
- * ===========================================================================
- * $Log: NNI.cpp,v $
- * Revision 1000.1 2004/06/01 18:09:44 gouriano
- * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
- *
- * Revision 1.2 2004/05/21 21:41:03 gorelenk
- * Added PCH ncbi_pch.hpp
- *
- * Revision 1.1 2004/02/10 15:16:00 jcherry
- * Initial version
- *
- * ===========================================================================
- */