dab89c1682e091dfa3780ce95359ebea9f6b7d52 braney Thu May 18 09:49:17 2017 -0700 add licensing info to optimalLeaf library diff --git src/optimalLeaf/Leaf.cc src/optimalLeaf/Leaf.cc index e3fcf05..f1d09c7 100644 --- src/optimalLeaf/Leaf.cc +++ src/optimalLeaf/Leaf.cc @@ -1,218 +1,224 @@ +/* This code provided by Ziv Bar-Joseph with the explicit understanding that + * no licensing is required for it to be used in the UCSC Genome Browser + * as long as reference is made to this paper: + * https://www.ncbi.nlm.nih.gov/pubmed/12801867 + */ + // -------------------------------------------------------------------- // Leaf.cc // This is where the optimal ordering for two subtrees // is comnputed. Uses the known optimal ordering for // each of the subtrees to generate the new ordering of the combined tree. // --------------------------------------------------------------------- //#include <iostream.h> #include <stdlib.h> #include "Leaf.hh" #include "general.hh" int compDist(const void *d1,const void *d2); // Pair for best tree construction LeafPair::LeafPair(int l,int r,LeafPair *pl,LeafPair *pr,int t1,int t2) { leftLeaf = l; rightLeaf = r; preLeft = pl; preRight = pr; n1 = t1; n2 = t2; } // Each gene is assigned a leaf, initially it does not have // any pairing gene. Leaf::Leaf(int num,float **mat) { index = num; distMat = mat; LeafPair *oneLeaf = new LeafPair(index,-1,0,0,1,0); curDist = new LeafDist*[1]; curDist[0] = new LeafDist(index,0,oneLeaf); listSize = 1; newDist = 0; newSize = 0; bestNew = maxInf; } // replace previous pair list (from previous subtree) with // a new one from the current subtree void Leaf::replace() { int i; // delete previous distance list for(i=0;i<listSize;i++) delete curDist[i]; delete []curDist; listSize = newSize; // bestNew becomes the maximum distance that may help in future // trees. bestNew += maxAdd; curDist = new LeafDist*[listSize]; for(i=0;i<newSize;i++) { curDist[i] = newDist[i]; } delete []newDist; qsort(curDist,listSize,sizeof(LeafDist*),compDist); newSize = 0; bestNew = maxInf; } // Adds corner to list of corners, and finds the best distance // with corner on the other side. void Leaf::addToNew(Leaf **corner1,Leaf **corner2,int n1,int n2,int c1,int c2,float max1,int tot) { int fromNum,fromIndex; LeafDist **fDist; Leaf *curLeaf; int i,j; float *bestC = new float[tot+1]; float curVal,bestD; float bestPos; int *bForC = new int[tot+1]; int cForD,dPlace; float maxAC = maxInf; for(j=0;j<c1;j++) { fromIndex = corner1[j]->giveIndex(); bestC[fromIndex] = maxInf; for(i=0;i<listSize;i++) { bestPos = curDist[i]->dist + max1; if(bestPos > bestC[fromIndex]) // optimization i = listSize; else { curVal = curDist[i]->dist + distMat[curDist[i]->n][fromIndex]; if(curVal < bestC[fromIndex]) { bestC[fromIndex] = curVal; bForC[fromIndex] = i; } } } if(bestC[fromIndex] < maxAC) maxAC = bestC[fromIndex]; } for(j=0;j<c2;j++) { bestD = maxInf; curLeaf = corner2[j]; fromNum = curLeaf->giveSize(); fromIndex = curLeaf->giveIndex(); fDist = curLeaf->giveList(); for(i=0;i<fromNum;i++) { bestPos = curLeaf->curDist[i]->dist + maxAC; if(bestPos > bestD) i = fromNum; // optimization else { curVal = bestC[fDist[i]->n] + curLeaf->curDist[i]->dist; if(curVal < bestD) { bestD = curVal; cForD = curLeaf->curDist[i]->n; dPlace = i; } } } LeafPair *newPair = new LeafPair(index,fromIndex, curDist[bForC[cForD]]->best,fDist[dPlace]->best,n1,n2); newDist[newSize] = new LeafDist(fromIndex,bestD,newPair); newSize++; LeafPair *cornerPair = new LeafPair(fromIndex,index, fDist[dPlace]->best,curDist[bForC[cForD]]->best,n2,n1); curLeaf->addNewDist(index,bestD,cornerPair); } delete []bForC; delete []bestC; } // adds a new pair to the dist list void Leaf::addNewDist(int n,float dist,LeafPair *p) { newDist[newSize] = new LeafDist(n,dist,p); newSize++; } // find the minimum distance pair (since the list are ordered, // this is always the first pair on the list). LeafPair *Leaf::findMin(float *min) { (*min) = curDist[0]->dist; return curDist[0]->best; } // used for sorting the lists int compDist(const void *d1,const void *d2) { LeafDist **l1 = (LeafDist**)d1; LeafDist **l2 = (LeafDist**)d2; return((*l1)->dist > (*l2)->dist); } // the optimzation for the last two subtrees, no need // to compute distance to all leaves, the best will suffice (see paper). int Leaf::findLast(Leaf **corner1,Leaf **corner2,int n1,int n2,int tot) { LeafDist **myDist; Leaf *curLeaf; int i,j; float curVal,best = maxInf; float myVal; int myInd,bestIndl,bestIndr,mBestl,mBestr; LeafPair *lpre,*rpre; LeafDist ***fDist = new LeafDist**[n2]; for(i=0;i<n2;i++) { fDist[i] = corner2[i]->giveList(); } for(j=0;j<n1;j++) { curLeaf = corner1[j]; myDist = curLeaf->giveList(); myVal = myDist[0]->dist; myInd = curLeaf->giveIndex(); for(i=0;i<n2;i++) { curVal = myVal + fDist[i][0]->dist + distMat[myInd][corner2[i]->giveIndex()]; if(best > curVal) { best = curVal; bestIndl = myDist[0]->n; bestIndr = fDist[i][0]->n; mBestl = myInd; mBestr = corner2[i]->giveIndex(); } } } LeafPair *newPair; int place,size; for(i=0;i<n2;i++) { if(corner2[i]->giveIndex() == bestIndr) { size = corner2[i]->giveSize(); for(j=0;j<size;j++) { if(fDist[i][j]->n == mBestr) { rpre = fDist[i][j]->best; j = size; } } i = n2; } } for(i=0;i<n1;i++) { if(corner1[i]->giveIndex() == bestIndl) { size = corner1[i]->giveSize(); myDist = corner1[i]->giveList(); for(j=0;j<size;j++) { if(myDist[j]->n == mBestl) { lpre = myDist[j]->best; j = size; } } newPair = new LeafPair(bestIndl,bestIndr,lpre,rpre,n1,n2); place = i; corner1[i]->initNewSize(1); corner1[i]->addNewDist(bestIndr,best,newPair); corner1[i]->replace(); i = n1; } } delete []fDist; return place; }