85bd10da37f403d155c4434d90edcd146e682948
braney
  Mon May 15 13:03:32 2017 -0700
add sorting to composite wiggles.

diff --git src/optimalLeaf/Leaf.cc src/optimalLeaf/Leaf.cc
new file mode 100644
index 0000000..e3fcf05
--- /dev/null
+++ src/optimalLeaf/Leaf.cc
@@ -0,0 +1,218 @@
+// --------------------------------------------------------------------
+//                         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;
+}
+
+
+
+