diff --git a/Bio/Phylogeny/PhyBin.hs b/Bio/Phylogeny/PhyBin.hs
--- a/Bio/Phylogeny/PhyBin.hs
+++ b/Bio/Phylogeny/PhyBin.hs
@@ -9,12 +9,12 @@
 
 module Bio.Phylogeny.PhyBin
        ( driver, binthem, normalize, annotateWLabLists, unitTests, acquireTreeFiles,
-         deAnnotate )
+         deAnnotate, retrieveHighlights, matchAnyHighlight )
        where
 
 import qualified Data.Foldable as F
 import           Data.Function       (on)
-import           Data.List           (delete, minimumBy, sortBy, foldl1')
+import           Data.List           (delete, minimumBy, sortBy, foldl1', foldl', intersperse, isPrefixOf)
 import           Data.Maybe          (fromMaybe)
 import           Data.Either         (partitionEithers)
 import           Data.Time.Clock
@@ -24,14 +24,15 @@
 import qualified Data.Vector                 as V
 import qualified Data.Vector.Unboxed         as U
 import           Control.Monad       (forM, forM_, filterM, when, unless)
-import           Control.Concurrent.Async
+import qualified Control.Concurrent.Async as Async
 import           Control.Exception   (evaluate)
 import           Control.Applicative ((<$>),(<*>))
 import           Control.Concurrent  (Chan)
 import           System.FilePath     (combine)
-import           System.Directory    (doesFileExist, doesDirectoryExist,
+import           System.Directory    (doesFileExist, doesDirectoryExist, createDirectoryIfMissing,
                                       getDirectoryContents, getCurrentDirectory)
 import           System.IO           (openFile, hClose, IOMode(..), stdout)
+import           System.Info         (os)
 import           System.Process      (system)
 import           System.Exit         (ExitCode(..))
 import           Test.HUnit          ((~:),(~=?),Test,test)
@@ -48,18 +49,34 @@
 import           Bio.Phylogeny.PhyBin.Util
 
 import Debug.Trace
+----------------------------------------------------------------------------------------------------
 
 -- Turn on for extra invariant checking:
 debug :: Bool
 debug = True
 
+#ifdef SEQUENTIALIZE
+#warning "SEQUENTIALIING execution.  Disabling all parallelism."
+type Async t = t
+async a = a
+wait  x = return x 
+#else
+type Async t = Async.Async t
+async = Async.async
+wait  = Async.wait
+#endif
+
+-- | A dendrogram PLUS consensus trees at the intermediate points.
+data DendroPlus a = DPLeaf (FullTree a)
+                  | DPBranch !Double (FullTree ()) (DendroPlus a) (DendroPlus a)
+
 ----------------------------------------------------------------------------------------------------
 
 -- | Driver to put all the pieces together (parse, normalize, bin)
 driver :: PhyBinConfig -> IO ()
-driver PBC{ verbose, num_taxa, name_hack, output_dir, inputs=files,
-            do_graph, branch_collapse_thresh,
-            dist_thresh, clust_mode, print_rfmatrix } =
+driver cfg@PBC{ verbose, num_taxa, name_hack, output_dir, inputs=files,
+                do_graph, branch_collapse_thresh, highlights, 
+                dist_thresh, clust_mode, use_hashrf, print_rfmatrix } =
    -- Unused: do_draw
  do 
     --------------------------------------------------------------------------------
@@ -69,17 +86,18 @@
     --putStrLn$ "PHYBIN RUNNING IN DIRECTORY: "++ cd
 
     bl <- doesDirectoryExist output_dir
-    unless bl $ do
-      c <- system$ "mkdir -p "++output_dir
-      case c of
-        ExitSuccess     -> return ()
-        ExitFailure cde -> error$"Could not create output directory. 'mkdir' command failed with: "++show cde
+    unless bl $ createDirectoryIfMissing True output_dir
 
-    putStrLn$ "Cleaning away previous phybin outputs..."
-    system$ "rm -f "++output_dir++"/dendrogram.*"
-    system$ "rm -f "++output_dir++"/cluster*"
-    system$ "rm -f "++output_dir++"/distance_matrix.txt"
-    system$ "rm -f "++output_dir++"/WARNINGS.txt"
+    if isPrefixOf "mingw" os then
+      -- TODO: Need a portable version of this.  'filemanip' would do:
+      putStrLn$ "Not cleaning away previous phybin outputs (TODO: port this to Windows)."
+     else do 
+      putStrLn$ "Cleaning away previous phybin outputs..."
+      system$ "rm -f "++output_dir++"/dendrogram.*"
+      system$ "rm -f "++output_dir++"/cluster*"
+      system$ "rm -f "++output_dir++"/distance_matrix.txt"
+      system$ "rm -f "++output_dir++"/WARNINGS.txt"
+      return ()
 
     putStrLn$ "Parsing "++show (length files)++" Newick tree files."
     --putStrLn$ "\nFirst ten \n"++ concat (map (++"\n") $ map show $ take 10 files)
@@ -97,7 +115,12 @@
     let num_files = length goodFiles
     bstrs <- mapM B.readFile goodFiles
     let (labelTab, fulltrees) = parseNewicks name_hack (zip files bstrs)
-
+          
+    highlightTrs <- retrieveHighlights name_hack labelTab highlights
+        
+    putStrLn$ "\nTotal unique taxa ("++ show (M.size labelTab) ++"):\n  "++
+	      (unwords$ M.elems labelTab)
+--	      show (nest 2 $ sep $ map text $ M.elems labelTab)    
     --------------------------------------------------------------------------------
 
     case branch_collapse_thresh of 
@@ -138,10 +161,11 @@
       putStrLn$ "Number of bad/unreadable input tree files: " ++ show (length warnings2)
     putStrLn$ "Number of VALID trees (correct # of leaves/taxa): " ++ show (length validtrees)
     putStrLn$ "Total tree nodes contained in valid trees: "++ show (sum counts)
-    let all_branches = concatMap (F.foldr' (\x ls -> getBranchLen x:ls) []) validtrees
-    putStrLn$ "Average branch len over valid trees: "++ show (avg all_branches)
-    putStrLn$ "Max/Min branch lengths: "++ show (foldl1' max all_branches,
-                                                 foldl1' min all_branches)
+    let all_branches = concatMap (F.foldr' (\x ls -> getBranchLen x:ls) []) validtrees    
+    unless (null all_branches) $ do
+      putStrLn$ "Average branch len over valid trees: "++ show (avg all_branches)
+      putStrLn$ "Max/Min branch lengths: "++ show (foldl1' max all_branches,
+                                                   foldl1' min all_branches)
 
     --------------------------------------------------------------------------------
     -- Next, dispatch on the mode and do the actual clustering or binning.
@@ -150,54 +174,48 @@
     (classes,binlist,asyncs) <- case clust_mode of
       BinThem         -> do x <- doBins validtrees
                             -- A list of bins, sorted by size:
-                            let binlist = reverse $ sortBy (compare `on` fst3) $
-                                          map (\ (tr,OneCluster ls) -> (length ls, tr, OneCluster ls)) $
+                            let binlist = reverse $ sortBy (compare `on` fst) $
+                                          map (\ (tr,OneCluster ls) -> (length ls, OneCluster ls)) $
                                           M.toList x
                             return (x,binlist,[])
-      ClusterThem lnk -> do
-        (mat, dendro) <- doCluster lnk validtrees        
-        case print_rfmatrix of
-          False -> return ()
-          True -> do -- treeFiles <- acquireTreeFiles files
-                     -- let fn f = do raw <- B.readFile f
-                     --               let ls = map (`B.append` (B.pack ";")) $ 
-                     --                        B.splitWith (== ';') raw
-                     --               return (map (f,) ls)
-                     -- trees0 <- concat <$> mapM fn treeFiles
-                     -- FIXME: no name_hack here:
-                     -- let (lbls, trees) = parseNewicks id trees0 
-                     -- putStrLn$ "Read trees! "++show (length trees)
-                     -- putStrLn$ "Taxa: "++show (pPrint lbls)
-                     -- putStrLn$ "First tree: "++show (displayDefaultTree (head trees))
-                     printDistMat stdout mat
-        writeFile (combine output_dir ("dendrogram.txt"))
-                  (show$ fmap treename dendro)
-        putStrLn "Wrote full dendrogram to file dendrogram.txt"
-
-        gvizAsync <- async $ do
-          t0 <- getCurrentTime
-          let dot = dotDendrogram "dendrogram" 1.0 dendro
-          _ <- dotToPDF dot (combine output_dir "dendrogram.pdf")
-          t1 <- getCurrentTime          
-          putStrLn$ "Wrote dendrogram diagram to file dendrogram.pdf ("++show(diffUTCTime t1 t0)++")"
-
+      ClusterThem{linkage} -> do
+        (mat, dendro) <- doCluster use_hashrf num_taxa linkage validtrees        
+        when print_rfmatrix $ printDistMat stdout mat
         hnd <- openFile  (combine output_dir ("distance_matrix.txt")) WriteMode
         printDistMat hnd mat
         hClose hnd
-
+        --------------------
+        writeFile (combine output_dir ("dendrogram.txt"))
+                  (show$ fmap treename dendro)
+        putStrLn " [finished] Wrote full dendrogram to file dendrogram.txt"
+        let plotIt mnameMap = 
+              if True -- do_graph
+              then async (do 
+                t0 <- getCurrentTime
+                let dot = dotDendrogram cfg "dendrogram" 1.0 dendro mnameMap highlightTrs
+                _ <- dotToPDF dot (combine output_dir "dendrogram.pdf") 
+                t1 <- getCurrentTime          
+                putStrLn$ " [finished] Wrote dendrogram diagram to file dendrogram.pdf ("
+                          ++show(diffUTCTime t1 t0)++")")
+              else async (return ())
         case dist_thresh of
-          Nothing -> error "Fully hierarchical cluster output is not finished!  Use --editdist."
-          Just dstThresh -> do
+          Nothing -> do a <- plotIt Nothing
+                        return (M.empty,[],[a])
+                     -- error "Fully hierarchical cluster output is not finished!  Use --editdist."
+          Just dstThresh -> do            
             putStrLn$ "Combining all clusters at distance less than or equal to "++show dstThresh
             let clusts = sliceDendro (fromIntegral dstThresh) dendro
-                wlens  = map (\ (OneCluster l) -> (length l, error "need consensus tree", OneCluster l)) clusts
-                sorted0 = reverse$ sortBy (compare `on` fst3) wlens
---                sorted  = map thd3 sorted0
-            putStrLn$ "After flattening, cluster sizes are: "++show (map fst3 sorted0)
+                classes = clustsToMap clusts
+                -- Cache the lengths:
+                wlens  = [ (length ls, OneCluster ls) | OneCluster ls <- clusts ]
+                sorted0 = reverse$ sortBy (compare `on` fst) wlens -- TODO: Parallel sorting?
+                nameMap = clustsToNameMap (map snd sorted0)
+            gvizAsync <- plotIt (Just nameMap)
+            putStrLn$ "After flattening, cluster sizes are: "++show (map fst sorted0)
             -- Flatten out the dendogram:
-            return (clustsToMap clusts, sorted0, [gvizAsync])
+            return (classes, sorted0, [gvizAsync])
 
-    reportClusts clust_mode binlist
+    unless (null binlist)$ reportClusts clust_mode binlist
         
     ----------------------------------------
     -- TEST, TEMPTOGGLE: print out edge weights :
@@ -212,10 +230,7 @@
     -- Finally, produce all the required outputs.
     --------------------------------------------------------------------------------
 
-    putStrLn$ "\nTotal unique taxa ("++ show (M.size labelTab) ++"):\n"++ 
-	      show (nest 2 $ sep $ map text $ M.elems labelTab)
-
-    putStrLn$ "Final number of tree bins: "++ show (M.size classes)
+--    putStrLn$ "Final number of tree bins: "++ show (M.size classes)
 
     unless (null warnings1 && null warnings2) $
 	writeFile (combine output_dir "WARNINGS.txt")
@@ -230,13 +245,13 @@
                                  "Wrong number of taxa ("++ show n ++"): "++ file++"\n")
 		           warnings2))
 
-    case clust_mode of
-      BinThem         -> outputBins binlist output_dir do_graph
-      ClusterThem lnk -> outputClusters binlist output_dir do_graph
+    async2 <- case clust_mode of
+                BinThem       -> outputBins binlist output_dir do_graph
+                ClusterThem{} -> outputClusters num_taxa binlist output_dir do_graph
 
     -- Wait on parallel tasks:
     putStrLn$ "Waiting for asynchronous tasks to finish..."
-    mapM_ wait asyncs 
+    mapM_ wait (async2:asyncs)
     putStrLn$ "Finished."
     --------------------------------------------------------------------------------
     -- End driver
@@ -257,12 +272,15 @@
 	          binthem validtrees
     return (classes)
 
-doCluster :: C.Linkage -> [FullTree a] -> IO (DistanceMatrix, C.Dendrogram (FullTree a))
-doCluster linkage validtrees = do
-  putStrLn$ "Clustering using method "++show linkage
+doCluster :: Bool -> Int -> C.Linkage -> [FullTree a] -> IO (DistanceMatrix, C.Dendrogram (FullTree a))
+doCluster use_hashrf num_taxa linkage validtrees = do
+  t0 <- getCurrentTime
+  when use_hashrf$ putStrLn " Using HashRF-style algorithm..."
   let nwtrees  = map nwtree validtrees
-      numtrees = length validtrees 
-      mat      = distanceMatrix nwtrees
+      numtrees = length validtrees
+      mat = if use_hashrf 
+            then hashRF num_taxa nwtrees
+            else fst (distanceMatrix nwtrees)
       ixtrees  = zip [0..] validtrees
       dist (i,t1) (j,t2) | j == i     = 0
 --                         | i == numtrees-1 = 0 
@@ -270,29 +288,29 @@
                          | otherwise  = fromIntegral ((mat V.! j) U.! i)
       dist1 a b = trace ("Taking distance between "++show (fst a, fst b)) $ dist a b
       dendro = fmap snd $ C.dendrogram linkage ixtrees dist
-  -- putStrLn$ "Got numtrees ...  "++show numtrees
-  -- putStrLn$ "Got the distance matrix ...  "++show (V.length mat)
+  -- Force the distance matrix:
+  V.mapM_ evaluate mat
+  t1 <- getCurrentTime
+  putStrLn$ "Time to compute distance matrix: "++show(diffUTCTime t1 t0)
+  putStrLn$ "Clustering using method "++show linkage
   return (mat,dendro)
-  
--- reportClusts :: ClustMode -> BinResults StandardDecor -> IO [(Int, StrippedTree, OneCluster StandardDecor)]
+
+
+reportClusts :: ClustMode -> [(Int, OneCluster StandardDecor)] -> IO ()
 reportClusts mode binlist = do 
-    let 
-        taxa :: S.Set Int
-	taxa = S.unions$ map (S.fromList . all_labels . snd3) binlist
-        binsizes = map fst3 binlist
+    let binsizes = map fst binlist
 
     putStrLn$ " Outcome: "++show (length binlist)++" clusters found, "++show (length$ takeWhile (>1) binsizes)
              ++" non-singleton, top bin sizes: "++show(take 10 binsizes)
-    putStrLn$"  First 50 bin sizes, excluding singletons:"
-    forM_ (zip [1..50] binlist) $ \ (ind, (len, tr, OneCluster ftrees)) -> do
+    putStrLn$"  Up to first 30 bin sizes, excluding singletons:"
+    forM_ (zip [1..30] binlist) $ \ (ind, (len, OneCluster ftrees)) -> do
        when (len > 1) $ -- Omit that long tail of single element classes...
           putStrLn$show$
            hcat [text ("  * cluster#"++show ind++", members "++ show len ++", "), 
                  case mode of
                    BinThem -> vcat [text ("avg bootstraps "++show (get_bootstraps$ avg_trees$ map nwtree ftrees)++", "),
                                     text "all: " <> pPrint (filter (not . null) $ map (get_bootstraps . nwtree) ftrees)]
-                   ClusterThem _ -> hcat [] ]
-    return binlist
+                   ClusterThem{} -> hcat [] ]
 
 -- | Convert a flat list of clusters into a map from individual trees to clusters.
 clustsToMap :: [OneCluster StandardDecor] -> BinResults StandardDecor
@@ -303,6 +321,14 @@
     fn2 theclust acc (FullTree{nwtree}) =
       M.insert (anonymize_annotated nwtree) theclust acc
 
+-- | Map each tree NAME onto the one-based index (in sorted order) of the cluster it
+-- comes from.
+clustsToNameMap :: [OneCluster StandardDecor] -> M.Map TreeName Int
+clustsToNameMap clusts = foldl' fn M.empty (zip [1..] clusts)
+  where
+    fn acc (ix,(OneCluster ftrs)) =
+      foldl' (\acc' t -> M.insert (treename t) ix acc') acc ftrs
+
 flattenDendro :: (C.Dendrogram (FullTree DefDecor)) -> OneCluster StandardDecor
 flattenDendro dendro =
   case dendro of
@@ -327,65 +353,72 @@
 
 --------------------------------------------------------------------------------
 
-outputClusters binlist output_dir do_graph = do
+-- outputClusters :: (Num a1, Ord a1, Show a1) => Int -> [(a1, t, OneCluster a)] -> String -> Bool -> IO ()
+outputClusters :: Int -> [(Int, OneCluster a)] -> String -> Bool -> IO (Async ())
+outputClusters num_taxa binlist output_dir do_graph = do
     let numbins = length binlist
     let base i size = combine output_dir (filePrefix ++ show i ++"_"++ show size) 
-
-    forM_ (zip [1::Int ..] binlist) $ \ (i, (size, _tr, OneCluster ftrees)) -> do
+    let consTrs = [ FullTree "consensus" (labelTable $ head ftrees) nwtr 
+                  | (_, OneCluster ftrees) <- binlist
+                  , let nwtr :: NewickTree StandardDecor
+                        nwtr = annotateWLabLists $ fmap (const (Nothing,0)) $
+                               consensusTree num_taxa $ map nwtree ftrees ]
+    
+    forM_ (zip3 [1::Int ..] binlist consTrs) $ \ (i, (size, OneCluster ftrees), fullConsTr) -> do
        writeFile (base i size ++".txt") (concat$ map ((++"\n") . treename) ftrees)
-       -- TODO: CONSENSUS TREE:
-       -- writeFile   (base i size ++".tr")  (show (displayDefaultTree$ deAnnotate fullAvgTr) ++ ";\n") -- FIXME
+       writeFile (base i size ++"_consensus.tr") (show (displayStrippedTree fullConsTr) ++ "\n")
+       writeFile (base i size ++"_alltrees.tr")
+           (unlines [ show (displayStrippedTree ft) | ft <- ftrees ])
 
-    putStrLn$ "[finished] Wrote contents of each cluster to cluster<N>_<size>.txt"
-    putStrLn$ "           Wrote representative trees to cluster<N>_<size>.tr"  
-    when (do_graph) $ do
-      putStrLn$ "Next do the time consuming operation of writing out graphviz visualizations:"
-      forM_ (zip [1::Int ..] binlist) $ \ (i, (size, _tr, OneCluster membs)) -> do
+    putStrLn$ " [finished] Wrote contents of each cluster to cluster<N>_<size>.txt"
+    putStrLn$ " [finished] Wrote representative (consensus) trees to cluster<N>_<size>_consensus.tr"
+    if do_graph then do
+      putStrLn$ "Next start the time consuming operation of writing out graphviz visualizations:"
+      asyncs <- forM (zip3 [1::Int ..] binlist consTrs) $
+                \ (i, (size, OneCluster membs), fullConsTr) -> async$ do
     	 when (size > 1 || numbins < 100) $ do
-           -- TODO CONSENSUS TREE
---           let dot = dotNewickTree ("cluster #"++ show i) (1.0 / avg_branchlen (map nwtree membs)) fullAvgTr
---	   _ <- dotToPDF dot (base i size ++ ".pdf")
+           let dot = dotNewickTree ("cluster #"++ show i) 1.0 fullConsTr
+	   _ <- dotToPDF dot (base i size ++ ".pdf")
 	   return ()
-      putStrLn$ "[finished] Wrote visual representations of trees to "++filePrefix++"<N>_<size>.pdf"
-
-    return ()
-
+      async $ do
+        mapM_ wait asyncs
+        putStrLn$ " [finished] Wrote visual representations of consensus trees to "++filePrefix++"<N>_<size>.pdf"
+     else async (return ())
 
+outputBins :: [(Int, OneCluster StandardDecor)] -> String -> Bool -> IO (Async ())
 outputBins binlist output_dir  do_graph = do
     let numbins = length binlist
     let base i size = combine output_dir (filePrefix ++ show i ++"_"++ show size) 
-    let avgs = map (avg_trees . map nwtree . clustMembers . thd3) binlist
-    forM_ (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, _tr, OneCluster ftrees), avgTree) -> do
+    let avgs = map (avg_trees . map nwtree . clustMembers . snd) binlist
+    forM_ (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, OneCluster ftrees), avgTree) -> do
        let FullTree fstName labs _ = head ftrees
            fullAvgTr = FullTree fstName labs avgTree
-         
-       --putStrLn$ ("  WRITING " ++ combine output_dir (filePrefix ++ show i ++"_"++ show size ++".txt"))
        writeFile (base i size ++".txt") (concat$ map ((++"\n") . treename) ftrees)
-       -- writeFile (base i size ++".tr")  (show (pPrint tr) ++ ";\n")
-       -- Printing the average tree instead of the stripped cannonical one:
        when debug$ do
          writeFile (base i size ++".dbg") (show (pPrint avgTree) ++ "\n")
-       writeFile   (base i size ++".tr")  (show (displayDefaultTree$ deAnnotate fullAvgTr) ++ ";\n") -- FIXME
-
-    putStrLn$ "[finished] Wrote contents of each bin to "++filePrefix++"<N>_<binsize>.txt"
-    putStrLn$ "           Wrote representative trees to "++filePrefix++"<N>_<binsize>.tr" 
-    when (do_graph) $ do
-      putStrLn$ "Next do the time consuming operation of writing out graphviz visualizations:"
-      forM_ (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, _tr, OneCluster membs), avgTree) -> do
-         let FullTree fstName labs _ = head membs
-             fullAvgTr = FullTree fstName labs avgTree 
-    	 when (size > 1 || numbins < 100) $ do 
-           let dot = dotNewickTree ("cluster #"++ show i) (1.0 / avg_branchlen (map nwtree membs))
-                                   --(annotateWLabLists$ fmap (const 0) tr)
-                                   -- TEMP FIXME -- using just ONE representative tree:
-                                   ( --trace ("WEIGHTED: "++ show (head$ trees bentry)) $ 
-                                     --(head$ trees bentry) 
- 				    fullAvgTr)
-	   _ <- dotToPDF dot (base i size ++ ".pdf")
-	   return ()
-      putStrLn$ "[finished] Wrote visual representations of trees to bin<N>_<binsize>.pdf"
-
+       -- Printing the average tree instead of the stripped cannonical one:         
+       writeFile   (base i size ++"_avg.tr")  (show (displayDefaultTree$ deAnnotate fullAvgTr) ++ "\n") -- FIXME
 
+    putStrLn$ " [finished] Wrote contents of each bin to "++filePrefix++"<N>_<binsize>.txt"
+    putStrLn$ "            Wrote representative trees to "++filePrefix++"<N>_<binsize>_avg.tr" 
+    if do_graph then do
+       putStrLn$ "Next start the time consuming operation of writing out graphviz visualizations:"
+       asyncs <- forM (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, OneCluster membs), avgTree) -> async $ do
+          let FullTree fstName labs _ = head membs
+              fullAvgTr = FullTree fstName labs avgTree 
+          when (size > 1 || numbins < 100) $ do 
+            let dot = dotNewickTree ("cluster #"++ show i) (1.0 / avg_branchlen (map nwtree membs))
+                                    --(annotateWLabLists$ fmap (const 0) tr)
+                                    -- TEMP FIXME -- using just ONE representative tree:
+                                    ( --trace ("WEIGHTED: "++ show (head$ trees bentry)) $ 
+                                      --(head$ trees bentry) 
+                                     fullAvgTr)
+            _ <- dotToPDF dot (base i size ++ ".pdf")
+            return ()
+       async $ do
+         mapM_ wait asyncs
+         putStrLn$ " [finished] Wrote visual representations of trees to "++filePrefix++"<N>_<size>.pdf"
+     else async (return ())
 --------------------------------------------------------------------------------
 
 -- Monadic mapAccum
@@ -436,3 +469,30 @@
 type TempDecor = (Double, (Int, Int), Int, [Label])
 
 avg ls = sum ls / fromIntegral (length ls)
+
+-- | Parse trees in addition to the main inputs (for --highlight).
+retrieveHighlights :: (String->String) -> LabelTable -> [FilePath] -> IO [[NewickTree ()]]
+retrieveHighlights name_hack labelTab ls =
+  mapM parseHighlight ls
+  where
+    parseHighlight file = do 
+      bs <- B.readFile file
+      let (lt2,htr) = parseNewick labelTab name_hack file bs
+      unless (lt2 == labelTab) $
+        error$"Tree given as --highlight includes taxa not present in main tree set: "++
+              show(M.keys$ M.difference lt2 labelTab)            
+      return (map (fmap (const())) htr)
+
+
+-- | Create a predicate that tests trees for consistency with the set of --highlight
+-- (consensus) trees.
+matchAnyHighlight :: [[NewickTree ()]] -> NewickTree () -> Bool      
+-- matchAnyHighlight :: [[NewickTree ()]] -> NewickTree () -> Maybe Int
+-- If there is a match, return the index of the highlight that matched.
+matchAnyHighlight highlightTrs =
+   let matchers = map mkMatcher highlightTrs
+       mkMatcher ls = let fns = map compatibleWith ls -- Multiple predicate functions
+                      in \ tr -> -- Does it match any predicate?
+                          or$ map (\f -> f tr) fns   
+   in \ newtr -> 
+       any (\f -> f newtr) matchers
diff --git a/Bio/Phylogeny/PhyBin/Binning.hs b/Bio/Phylogeny/PhyBin/Binning.hs
--- a/Bio/Phylogeny/PhyBin/Binning.hs
+++ b/Bio/Phylogeny/PhyBin/Binning.hs
@@ -30,10 +30,8 @@
 import           System.Directory    (doesFileExist, doesDirectoryExist,
                                       getDirectoryContents, getCurrentDirectory)
 import           System.IO           (openFile, hClose, IOMode(ReadMode))
-import           System.Process      (system)
 import           System.Exit         (ExitCode(..))
 import           Test.HUnit          ((~:),(~=?),Test,test)
-import qualified HSH 
 
 -- For vizualization:
 import           Text.PrettyPrint.HughesPJClass hiding (char, Style)
@@ -179,13 +177,13 @@
 
 -- TODO: Salvage any of these tests that are worthwhile and get them into the unit tests:	        	
 tt :: AnnotatedTree
-tt = normalize $ annotateWLabLists $ snd$ parseNewick M.empty id "" "(A,(C,D,E),B);"
+tt = normalize $ annotateWLabLists $ head$ snd$ parseNewick M.empty id "" "(A,(C,D,E),B);"
 
 norm4 :: FullTree StandardDecor
 norm4 = norm "((C,D,E),B,A);"
 
 norm5 :: AnnotatedTree
-norm5 = normalize$ annotateWLabLists$ snd$ parseNewick M.empty id "" "(D,E,C,(B,A));"
+norm5 = normalize$ annotateWLabLists$ head$ snd$ parseNewick M.empty id "" "(D,E,C,(B,A));"
 
 
 ----------------------------------------------------------------------------------------------------
@@ -332,6 +330,8 @@
 annotateWLabLists :: NewickTree DefDecor -> AnnotatedTree
 annotateWLabLists tr = case tr of 
   NTLeaf (bs,bl) n      -> NTLeaf (StandardDecor bl bs 1 [n]) n
+  NTInterior (bs,bl) [] -> NTInterior (StandardDecor bl bs 0 []) []
+  -- NTInterior (bs,bl) [] -> error "annotateWLabLists: internal invariant broken.  Shouldn't have NTInterior with null children."
   NTInterior (bs,bl) ls -> 
       let children = map annotateWLabLists ls in 
       NTInterior (StandardDecor bl bs
@@ -372,7 +372,9 @@
 ----------------------------------------------------------------------------------------------------
 
 tre1 :: (LabelTable, NewickTree DefDecor)
-tre1 = parseNewick M.empty id "" "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
+tre1 = (x,head y)
+ where
+   (x,y) = parseNewick M.empty id "" "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
 
 tre1draw :: IO (Chan (), FullTree StandardDecor)
 tre1draw = viewNewickTree "tre1"$ (FullTree "" (fst tre1) (annotateWLabLists (snd tre1)))
@@ -384,7 +386,7 @@
 norm = norm2 . B.pack
 
 norm2 :: B.ByteString -> FullTree StandardDecor
-norm2 bstr = FullTree "" tbl (normalize $ annotateWLabLists tr)
+norm2 bstr = FullTree "" tbl (normalize $ annotateWLabLists$ head tr)
   where
     (tbl,tr) = parseNewick M.empty id "test" bstr
 
diff --git a/Bio/Phylogeny/PhyBin/CoreTypes.hs b/Bio/Phylogeny/PhyBin/CoreTypes.hs
--- a/Bio/Phylogeny/PhyBin/CoreTypes.hs
+++ b/Bio/Phylogeny/PhyBin/CoreTypes.hs
@@ -12,7 +12,7 @@
          ClustMode(..), TreeName,
          
          -- * Tree operations
-         displayDefaultTree,
+         displayDefaultTree, displayStrippedTree, 
          treeSize, numLeaves, liftFT,
          get_dec, set_dec, get_children, 
          map_labels, all_labels, foldIsomorphicTrees,
@@ -21,7 +21,7 @@
          avg_branchlen, get_bootstraps,
 
          -- * Command line config options
-         PhyBinConfig(..), default_phybin_config,
+         PhyBinConfig(..), default_phybin_config, 
 
          -- * General helpers
          Label, LabelTable,
@@ -85,6 +85,9 @@
 instance Foldable FullTree where
   foldMap f (FullTree _ _ tr) = foldMap f tr
 
+instance Functor FullTree where
+  fmap f (FullTree n l tr) = FullTree n l $ fmap f tr
+
 instance Pretty dec => Pretty (NewickTree dec) where 
  -- I'm using displayDefaultTree for the "prettiest" printing and
  -- replacing pPrint with a whitespace-improved version of show:
@@ -114,6 +117,13 @@
          Just val -> base <> text ":[" <> text (show val) <> text "]"
       where base = parens$ sep$ map_but_last (<>text",") $ map loop ls
 
+displayStrippedTree :: FullTree a -> Doc
+displayStrippedTree orig = loop tr <> ";"
+  where
+    (FullTree _ mp tr) = orig -- normalize orig
+    loop (NTLeaf _ name) = text (mp M.! name)
+    loop (NTInterior _ ls) = parens$ sep$ map_but_last (<>text",") $ map loop ls
+
 ----------------------------------------------------------------------------------------------------
 -- Labels
 ----------------------------------------------------------------------------------------------------
@@ -177,7 +187,7 @@
            , labelTable :: LabelTable
            , nwtree     :: NewickTree a 
            }
- deriving (Show)
+ deriving (Show, Ord, Eq)
 
 liftFT :: (NewickTree t -> NewickTree a) -> FullTree t -> FullTree a
 liftFT fn (FullTree nm labs x) = FullTree nm labs (fn x)
@@ -205,6 +215,10 @@
       , do_graph :: Bool
       , do_draw :: Bool
       , clust_mode  :: ClustMode
+      , highlights  :: [FilePath] -- [FullTree ()]
+      , show_trees_in_dendro :: Bool
+      , show_interior_consensus :: Bool
+      , use_hashrf  :: Bool  
       , print_rfmatrix :: Bool
       , dist_thresh :: Maybe Int
       , branch_collapse_thresh :: Maybe Double -- ^ Branches less than this length are collapsed.
@@ -214,19 +228,28 @@
 default_phybin_config :: PhyBinConfig
 default_phybin_config = 
  PBC { verbose = False
-      , num_taxa = error "must be able to determine the number of taxa expected in the dataset.  (Supply it manually.)"
+      , num_taxa = error "must be able to determine the number of taxa expected in the dataset.  (Supply it manually with -n.)"
       , name_hack = id -- Default, no transformation of leaf-labels
       , output_dir = "./phybin_out/"
       , inputs = []
       , do_graph = False
       , do_draw = False
-      , clust_mode = BinThem
+--      , clust_mode = BinThem
+      , clust_mode = ClusterThem C.UPGMA
+#ifdef USE_HASHRF
+      , use_hashrf = True
+#else
+      , use_hashrf = False
+#endif
+      , highlights     = []
+      , show_trees_in_dendro = False
+      , show_interior_consensus = False
       , print_rfmatrix = False
       , dist_thresh = Nothing
       , branch_collapse_thresh = Nothing
      }
 
-data ClustMode = BinThem | ClusterThem C.Linkage 
+data ClustMode = BinThem | ClusterThem { linkage :: C.Linkage }
 
 ----------------------------------------------------------------------------------------------------
 -- * Simple utility functions for the core types:
diff --git a/Bio/Phylogeny/PhyBin/Parser.hs b/Bio/Phylogeny/PhyBin/Parser.hs
--- a/Bio/Phylogeny/PhyBin/Parser.hs
+++ b/Bio/Phylogeny/PhyBin/Parser.hs
@@ -24,12 +24,24 @@
 -- | Parse a bytestring into a NewickTree with branch lengths.  The
 --   first argument is file from which the data came and is just for
 --   better error messages.
-parseNewick :: LabelTable -> NameHack -> String -> B.ByteString -> (LabelTable, NewickTree DefDecor)
+parseNewick :: LabelTable -> NameHack -> String -> B.ByteString -> (LabelTable, [NewickTree DefDecor])
 parseNewick tbl0 name_hack file input =
   extractLabelTable tbl0 $ 
-  runB file (newick_parser name_hack) $
+  runB file (many1$ newick_parser name_hack) $
   B.filter (not . isSpace) input
 
+-- treeFiles <- acquireTreeFiles files
+                     -- let fn f = do raw <- B.readFile f
+                     --               let ls = map (`B.append` (B.pack ";")) $ 
+                     --                        B.splitWith (== ';') raw
+                     --               return (map (f,) ls)
+                     -- trees0 <- concat <$> mapM fn treeFiles
+                     -- FIXME: no name_hack here:
+                     -- let (lbls, trees) = parseNewicks id trees0 
+                     -- putStrLn$ "Read trees! "++show (length trees)
+                     -- putStrLn$ "Taxa: "++show (pPrint lbls)
+                     -- putStrLn$ "First tree: "++show (displayDefaultTree (head trees))
+
 -- | Parse a list of trees, starting with an empty map of labels and accumulating a final map.
 parseNewickFiles :: NameHack -> [String] -> IO (LabelTable, [FullTree DefDecor])
 parseNewickFiles fn nms = do
@@ -41,34 +53,43 @@
 parseNewicks name_hack pairs = (labtbl, fullTrs)
  where
    fullTrs = [ FullTree (takeBaseName file) labtbl tr
-             | (file,_) <- pairs
-             | tr       <- trs ]
+             | (file,tr) <- trs ]
    (labtbl, trs) = P.foldr fn (M.empty,[]) pairs
    fn (file,bstr) (!acc,!ls) =
-     let (acc',tr) = parseNewick acc name_hack file bstr
-     in (acc', tr:ls)
+     let (acc',trs) = parseNewick acc name_hack file bstr
+         names      = if length trs > 1
+                      then zipWith (\x y -> x++"_"++show y) (repeat file) [0..]
+                      else [file]
+     in (acc', (zip names trs) ++ ls)
 
 runB :: Show a => String -> Parser a -> B.ByteString -> a
 runB file p input = case (parse p "" input) of
 	         Left err -> error ("parse error in file "++ show file ++" at "++ show err)
 		 Right x  -> x
 
-extractLabelTable :: LabelTable -> TempTree -> (LabelTable, NewickTree DefDecor)
-extractLabelTable tbl0 tr = (finMap,finTree)
+extractLabelTable :: LabelTable -> [TempTree] -> (LabelTable, [NewickTree DefDecor])
+extractLabelTable tbl0 trs = (finMap,finTree)
  where
    flipped = M.fromList $ L.map (\(x,y)->(y,x)) $ M.toList tbl0
-   -- (_,finMap,finTree) = loop (S.fromList (M.elems tbl0)) tbl0 tr
-   (_,finMap,finTree) = loop flipped tbl0 tr
-   
-   loop seen acc (NTLeaf (d,Just nm) _)
+   -- (_,finMap,finTree) = loop2 (S.fromList (M.elems tbl0)) tbl0 tr
+   (_,finMap,finTree) = loop1 flipped tbl0 trs
+
+   -- Mere plumbing:
+   loop1 seen acc [] = (seen, acc, [])
+   loop1 seen acc (hd:tl) =
+     let (seen2,acc2,hd') = loop2 seen acc hd 
+         (seen3,acc3,tl') = loop1 seen2 acc2 tl in
+     (seen3,acc3, hd':tl')
+
+   loop2 seen acc (NTLeaf (d,Just nm) _)
      | M.member nm seen = (seen, acc, NTLeaf d (seen M.! nm))
      | otherwise = let nxt = M.size acc in
                    (M.insert nm nxt seen,
                     M.insert nxt nm acc,  NTLeaf d nxt)
-   loop seen1 acc1 (NTInterior (d,Nothing) chlds) =
+   loop2 seen1 acc1 (NTInterior (d,Nothing) chlds) =
      let (seen',acc',ls') = 
           P.foldr (\ x (seen2,acc2,ls) ->
-                   let (seen3,acc3,x') = loop seen2 acc2 x in
+                   let (seen3,acc3,x') = loop2 seen2 acc2 x in
                    (seen3, acc3, x':ls))
                   (seen1,acc1,[])
                   chlds
@@ -173,7 +194,7 @@
 --------------------------------------------------------------------------------
 
 tre1 :: NewickTree DefDecor
-tre1 = snd $ parseNewick M.empty id "" $ B.pack "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
+tre1 = head $ snd $ parseNewick M.empty id "" $ B.pack "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
 
 unitTests :: Test
 unitTests = test
diff --git a/Bio/Phylogeny/PhyBin/RFDistance.hs b/Bio/Phylogeny/PhyBin/RFDistance.hs
--- a/Bio/Phylogeny/PhyBin/RFDistance.hs
+++ b/Bio/Phylogeny/PhyBin/RFDistance.hs
@@ -1,20 +1,36 @@
-{-# LANGUAGE ScopedTypeVariables, CPP #-}
+{-# LANGUAGE ScopedTypeVariables, CPP, BangPatterns #-}
 
 module Bio.Phylogeny.PhyBin.RFDistance
-       (DenseLabelSet, DistanceMatrix, 
-        allBips, foldBips,
-        distanceMatrix, printDistMat)
+       (
+         -- * Types
+         DenseLabelSet, DistanceMatrix,
+
+         -- * Bipartition (Bip) utilities
+         allBips, foldBips, dispBip,
+         consensusTree, bipsToTree, filterCompatible, compatibleWith,
+
+         -- * ADT for dense sets
+         mkSingleDense, mkEmptyDense, bipSize,
+         denseUnions, denseDiff, invertDense, markLabel,
+         
+        -- * Methods for computing distance matrices
+        distanceMatrix, hashRF, 
+
+        -- * Output
+        printDistMat)
        where
 
 import           Control.Monad
+import           Control.Monad.ST
+import           Data.Function       (on)
 import           Data.Word
 import qualified Data.Vector                 as V
-import qualified Data.Vector.Unboxed.Mutable as MV
+import qualified Data.Vector.Mutable         as MV
+import qualified Data.Vector.Unboxed.Mutable as MU
 import qualified Data.Vector.Unboxed         as U
-import qualified Data.Vector.Unboxed.Bit     as UB
-import qualified Data.Bit                    as B
 import           Text.PrettyPrint.HughesPJClass hiding (char, Style)
 import           System.IO      (hPutStrLn, hPutStr, Handle)
+import           System.IO.Unsafe
 
 -- import           Control.LVish
 -- import qualified Data.LVar.Set   as IS
@@ -26,17 +42,33 @@
 import           Bio.Phylogeny.PhyBin.CoreTypes
 -- import           Data.BitList
 import qualified Data.Set as S
+import qualified Data.List as L
 import qualified Data.IntSet as SI
 import qualified Data.Map.Strict as M
 import qualified Data.Foldable as F
+import qualified Data.Traversable as T
 import           Data.Monoid
 import           Prelude as P
 import           Debug.Trace
 
+#ifdef BITVEC_BIPS
+import qualified Data.Vector.Unboxed.Bit     as UB
+import qualified Data.Bit                    as B
+#endif
+
+-- I don't understand WHY, but I seem to get the same answers WITHOUT this.
+-- Normalization and symmetric difference do make things somewhat slower (e.g. 1.8
+-- seconds vs. 2.2 seconds for 150 taxa / 100 trees)
+#define NORMALIZATION
+-- define BITVEC_BIPS
+
 --------------------------------------------------------------------------------
 -- A data structure choice
 --------------------------------------------------------------------------------
 
+-- type DenseLabelSet s = BitList
+
+
 -- | Dense sets of taxa, aka Bipartitions or BiPs
 --   We assume that taxa labels have been mapped onto a dense, contiguous range of integers [0,N).
 -- 
@@ -49,27 +81,72 @@
 -- 
 --   A set that is more than a majority of the taxa can be normalized by "flipping",
 --   i.e. taking the taxa that are NOT in that set.
-type DenseLabelSet = SI.IntSet 
--- type DenseLabelSet s = BitList
--- type DenseLabelSet = UB.Vector B.Bit
+#ifdef BITVEC_BIPS
 
--- M.write vec lab (B.fromBool True)
--- mkEmptyDense size = U.replicate size (B.fromBool False)    
+#  if 1
+type DenseLabelSet = UB.Vector B.Bit
+markLabel lab = UB.modify (\vec -> MU.write vec lab (B.fromBool True)) 
+mkEmptyDense  size = UB.replicate size (B.fromBool False)
+mkSingleDense size ind = markLabel ind (mkEmptyDense size)
+denseUnions        = UB.unions
+bipSize            = UB.countBits
+denseDiff          = UB.difference
+invertDense size bip = UB.invert bip
+dispBip labs bip = show$ map (\(ix,_) -> (labs M.! ix)) $
+                        filter (\(_,bit) -> B.toBool bit) $
+                        zip [0..] (UB.toList bip)
+denseIsSubset a b = UB.or (UB.difference b a)
+traverseDense_ fn bip =
+  U.ifoldr' (\ ix bit acc ->
+              (if B.toBool bit
+               then fn ix
+               else return ()) >> acc)
+        (return ()) bip
 
--- markLabel lab set = IS.putInSet lab set 
--- mkEmptyDense _size = IS.newEmptySet
+#  else
+-- TODO: try tracking the size:
+data DenseLabelSet = DLS {-# UNPACK #-} !Int (UB.Vector B.Bit)
+markLabel lab (DLS _ vec)= DLS (UB.modify (\vec -> return (MU.write vec lab (B.fromBool True))) ) vec
+-- ....
+#  endif
 
+#else
+type DenseLabelSet = SI.IntSet
+markLabel lab set   = SI.insert lab set 
+mkEmptyDense _size  = SI.empty
+mkSingleDense _size = SI.singleton
+denseUnions _size   = SI.unions 
+bipSize             = SI.size
+denseDiff           = SI.difference
+denseIsSubset       = SI.isSubsetOf
+
+dispBip labs bip = "[" ++ unwords strs ++ "]"
+  where strs = map (labs M.!) $ SI.toList bip
+invertDense size bip = loop SI.empty (size-1)
+  where -- There's nothing for it but to iterate and test for membership:
+    loop !acc ix | ix < 0           = acc
+                 | SI.member ix bip = loop acc (ix-1)
+                 | otherwise        = loop (SI.insert ix acc) (ix-1)
+traverseDense_ fn bip =
+  -- FIXME: need guaranteed non-allocating way to do this.
+  SI.foldr' (\ix acc ->  fn ix >> acc) (return ()) bip
+#endif
+
 markLabel    :: Label -> DenseLabelSet -> DenseLabelSet
 mkEmptyDense :: Int -> DenseLabelSet
+mkSingleDense :: Int -> Label -> DenseLabelSet
 denseUnions  :: Int -> [DenseLabelSet] -> DenseLabelSet
 bipSize      :: DenseLabelSet -> Int
 
-markLabel lab set  = SI.insert lab set 
-mkEmptyDense _size = SI.empty
-denseUnions _size  = SI.unions 
-bipSize            = SI.size
+-- | Print a BiPartition in a pretty form
+dispBip      :: LabelTable -> DenseLabelSet -> String
 
+-- | Assume that total taxa are 0..N-1, invert membership:
+invertDense  :: Int -> DenseLabelSet -> DenseLabelSet
 
+traverseDense_ :: Monad m => (Int -> m ()) -> DenseLabelSet -> m ()
+
+
 --------------------------------------------------------------------------------
 -- Dirt-simple reference implementation
 --------------------------------------------------------------------------------
@@ -77,20 +154,26 @@
 type DistanceMatrix = V.Vector (U.Vector Int)
 
 -- | Returns a triangular distance matrix encoded as a vector.
-distanceMatrix :: [NewickTree a] -> DistanceMatrix
+--   Also return the set-of-BIPs representation for each tree.
+distanceMatrix :: [NewickTree a] -> (DistanceMatrix, V.Vector (S.Set DenseLabelSet))
 distanceMatrix lst = 
    let sz = P.length lst
        eachbips = V.fromList $ map allBips lst
---   in V.generate (sz-1) $ \ i ->
-   in V.generate sz $ \ i ->        
-      U.generate i  $ \ j ->
-      S.size (S.difference (eachbips V.! i) (eachbips V.! j))
-  
+       mat = V.generate sz $ \ i ->        
+             U.generate i  $ \ j ->
+             let diff1 = S.size (S.difference (eachbips V.! i) (eachbips V.! j))
+                 diff2 = S.size (S.difference (eachbips V.! j) (eachbips V.! i))
+             in diff1 + diff2
+   in (mat, eachbips)
+
 -- | The number of bipartitions implied by a tree is one per EDGE in the tree.  Thus
 -- each interior node carries a list of BiPs the same length as its list of children.
 labelBips :: NewickTree a -> NewickTree (a, [DenseLabelSet])
 labelBips tr =
-    trace ("labelbips "++show allLeaves++" "++show size) $
+--    trace ("labelbips "++show allLeaves++" "++show size) $
+#ifdef NORMALIZATION  
+    fmap (\(a,ls) -> (a,map (normBip size) ls)) $
+#endif
     loop tr
   where    
     size = numLeaves tr
@@ -98,22 +181,28 @@
     loop (NTLeaf dec lab) = NTLeaf (dec, [markLabel lab zero]) lab      
     loop (NTInterior dec chlds) =
       let chlds' = map loop chlds
-          sets   = map (normBip . denseUnions size . snd . get_dec) chlds' in
+          sets   = map (denseUnions size . snd . get_dec) chlds' in
       NTInterior (dec, sets) chlds'
 
-    halfSize = size `quot` 2
-    normBip bip =
-      let flipped = SI.difference allLeaves bip in
-      case compare (SI.size bip) halfSize of
-        LT -> bip 
-        GT -> flipped -- Flip it
-        EQ -> -- This is a painful case, we need a tie-breaker
-              min bip flipped
-           
     allLeaves = leafSet tr
-    leafSet (NTLeaf _ lab)    = SI.singleton lab
+    leafSet (NTLeaf _ lab)    = mkSingleDense size lab
     leafSet (NTInterior _ ls) = denseUnions size $ map leafSet ls
 
+-- normBip :: DenseLabelSet -> DenseLabelSet -> DenseLabelSet
+--    normBip allLeaves bip =
+normBip :: Int -> DenseLabelSet -> DenseLabelSet    
+normBip totsize bip =
+  let -- size     = bipSize allLeaves
+      halfSize = totsize `quot` 2
+--      flipped  = denseDiff allLeaves bip
+      flipped  = invertDense totsize bip 
+  in 
+  case compare (bipSize bip) halfSize of
+    LT -> bip 
+    GT -> flipped -- Flip it
+    EQ -> min bip flipped -- This is a painful case, we need a tie-breaker
+    
+
 foldBips :: Monoid m => (DenseLabelSet -> m) -> NewickTree a -> m
 foldBips f tr = F.foldMap f' (labelBips tr)
  where f' (_,bips) = F.foldMap f bips
@@ -127,6 +216,7 @@
 --------------------------------------------------------------------------------
 -- First, necessary types:
 
+-- UNFINISHED:
 #if 0
 -- | A collection of all observed bipartitons (bips) with a mapping of which trees
 -- contain which Bips.
@@ -143,12 +233,78 @@
 -- | Tree's are identified simply by their order within the list of input trees.
 -- type TreeID = Int
 #endif
+
 --------------------------------------------------------------------------------
+-- Alternate way of slicing the problem: HashRF
+--------------------------------------------------------------------------------
 
 -- The distance matrix is an atomically-bumped matrix of numbers.
 -- type DistanceMat s = NA.NatArray s Word32
 -- Except... bump isn't supported by our idempotent impl.
 
+-- | This version slices the problem a different way.  A single pass over the trees
+-- populates the table of bipartitions.  Then the table can be processed (locally) to
+-- produce (non-localized) increments to a distance matrix.
+hashRF :: Int -> [NewickTree a] -> DistanceMatrix
+hashRF num_taxa trees = build M.empty (zip [0..] trees)
+  where
+    num_trees = length trees
+    -- First build the table:
+    build acc [] = ingest acc
+    build acc ((ix,hd):tl) =
+      let bips = allBips hd
+          acc' = S.foldl' fn acc bips
+          fn acc bip = M.alter fn2 bip acc
+          fn2 (Just membs) = Just (markLabel ix membs)
+          fn2 Nothing      = Just (mkSingleDense num_taxa ix)
+      in      
+      build acc' tl
+
+    -- Second, ingest the table to construct the distance matrix:
+    ingest :: M.Map DenseLabelSet DenseLabelSet -> DistanceMatrix
+    ingest bipTable = runST theST
+      where
+       theST :: forall s0 . ST s0 DistanceMatrix
+       theST = do 
+        -- Triangular matrix, starting narrow and widening:
+        matr <- MV.new num_trees
+        -- Too bad MV.replicateM is insufficient.  It should pass index.  
+        -- Instead we write this C-style:
+        for_ (0,num_trees) $ \ ix -> do 
+          row <- MU.replicate ix (0::Int)
+          MV.write matr ix row
+          return ()
+
+        unsafeIOToST$ putStrLn$" Built matrix for dim "++show num_trees
+
+        let bumpMatr i j | j < i     = incr i j
+                         | otherwise = incr j i
+            incr :: Int -> Int -> ST s0 ()
+            incr i j = do -- Not concurrency safe yet:
+--                          unsafeIOToST$ putStrLn$" Reading at position "++show(i,j)
+                          row <- MV.read matr i
+                          elm <- MU.read row j
+                          MU.write row j (elm+1)
+                          return ()
+            fn bipMembs =
+              -- Here we quadratically consider all pairs of trees and ask whether
+              -- their edit distance is increased based on this particular BiP.
+              -- Actually, as an optimization, it is sufficient to consider only the
+              -- cartesian product of those that have and those that don't.
+              let haveIt   = bipMembs
+                  -- Depending on how invertDense is written, it could be useful to
+                  -- fuse this in and deforest "dontHave".
+                  dontHave = invertDense num_trees bipMembs
+                  fn1 trId = traverseDense_ (fn2 trId) dontHave
+                  fn2 trId1 trId2 = bumpMatr trId1 trId2
+              in
+--                 trace ("Computed donthave "++ show dontHave) $ 
+                 traverseDense_ fn1 haveIt
+        F.traverse_ fn bipTable
+        v1 <- V.unsafeFreeze matr
+        T.traverse (U.unsafeFreeze) v1
+
+
 #if 0
 -- | Returns a (square) distance matrix encoded as a vector.
 distanceMatrix :: [AnnotatedTree] -> IO (U.Vector Word)
@@ -157,7 +313,7 @@
 --   IM.IMapSnap (table :: M.Map DenseLabelSet (Snapshot IS.ISet TreeID)) <- runParThenFreeze par
    IM.IMapSnap table <- runParThenFreeze par
    let sz = P.length lst
-   v <- MV.replicate (sz*sz) (0::Word)
+   v <- MU.replicate (sz*sz) (0::Word)
    let fn set () =
          
    F.foldrM 
@@ -181,6 +337,8 @@
 #endif
 
 --------------------------------------------------------------------------------
+-- Miscellaneous Helpers
+--------------------------------------------------------------------------------
 
 instance Pretty a => Pretty (S.Set a) where
  pPrint s = pPrint (S.toList s)
@@ -196,3 +354,82 @@
       hPutStr h " "
     hPutStr h "0\n"          
   hPutStrLn h "-----------------------------------------"
+
+-- My own forM for numeric ranges (not requiring deforestation optimizations).
+-- Inclusive start, exclusive end.
+{-# INLINE for_ #-}
+for_ :: Monad m => (Int, Int) -> (Int -> m ()) -> m ()
+for_ (start, end) _fn | start > end = error "for_: start is greater than end"
+for_ (start, end) fn = loop start
+  where
+   loop !i | i == end  = return ()
+           | otherwise = do fn i; loop (i+1)
+
+-- | Which of a set of trees are compatible with a consensus?
+filterCompatible :: NewickTree a -> [NewickTree b] -> [NewickTree b]
+filterCompatible consensus trees =
+    let cbips = allBips consensus in
+    [ tr | tr <- trees
+         , cbips `S.isSubsetOf` allBips tr ]
+
+-- | Is a tree compatible with a consensus?
+--   This is more efficient if partially applied then used repeatedly.
+compatibleWith :: NewickTree a -> NewickTree b -> Bool
+compatibleWith consensus newTr =
+  S.isSubsetOf (allBips consensus) (allBips newTr)
+
+-- Consensus between two trees, which may even have different label maps.
+consensusTreeFull (FullTree n1 l1 t1) (FullTree n2 l2 t2) =
+  error "FINISHME - consensusTreeFull"
+
+-- | Take only the bipartitions that are agreed on by all trees.
+consensusTree :: Int -> [NewickTree a] -> NewickTree ()
+consensusTree _ [] = error "Cannot take the consensusTree of the empty list"
+consensusTree num_taxa (hd:tl) = bipsToTree num_taxa intersection
+  where
+    intersection = L.foldl' S.intersection (allBips hd) (map allBips tl)
+--     intersection = loop (allBips hd) tl
+--     loop :: S.Set DenseLabelSet -> [NewickTree a] -> S.Set DenseLabelSet
+--     loop !remain []      = remain
+--     -- Was attempting to use foldBips here as an optimization:
+-- --     loop !remain (hd:tl) = loop (foldBips S.delete hd remain) tl
+--     loop !remain (hd:tl) = loop (S.difference remain (allBips hd)) tl    
+      
+-- | Convert from bipartitions BACK to a single tree.
+bipsToTree :: Int -> S.Set DenseLabelSet -> NewickTree ()
+bipsToTree num_taxa origbip =
+--  trace ("Doing bips in order: "++show sorted++"\n") $ 
+  loop lvl0 sorted
+  where
+    -- We consider each subset in increasing size order.
+    -- FIXME: If we tweak the order on BIPs, then we can just use S.toAscList here:
+    sorted = L.sortBy (compare `on` bipSize) (S.toList origbip)
+
+    lvl0 = [ (mkSingleDense num_taxa ix, NTLeaf () ix)
+           | ix <- [0..num_taxa-1] ]
+
+    -- VERY expensive!  However, due to normalization issues this is necessary for now:
+    -- TODO: in the future make it possible to definitively denormalize.
+    -- isMatch bip x = denseIsSubset x bip || denseIsSubset x (invertDense num_taxa bip)
+    isMatch bip x = denseIsSubset x bip 
+
+    -- We recursively glom together subtrees until we have a complete tree.
+    -- We only process larger subtrees after we have processed all the smaller ones.
+    loop !subtrees [] =
+      case subtrees of
+        []    -> error "bipsToTree: internal error"
+        [(_,one)] -> one
+        lst   -> NTInterior () (map snd lst)
+    loop !subtrees (bip:tl) =
+--      trace (" -> looping, subtrees "++show subtrees) $ 
+      let (in_,out) = L.partition (isMatch bip. fst) subtrees in
+      case in_ of
+        [] -> error $"bipsToTree: Internal error!  No match for bip: "++show bip
+              ++" out is\n "++show out++"\n and remaining bips "++show (length tl)
+              ++"\n when processing orig bip set:\n  "++show origbip
+          -- loop out tl
+        _ -> 
+         -- Here all subtrees that match the current bip get merged:
+         loop ((denseUnions num_taxa (map fst in_),
+                NTInterior ()        (map snd in_)) : out) tl
+
diff --git a/Bio/Phylogeny/PhyBin/Util.hs b/Bio/Phylogeny/PhyBin/Util.hs
--- a/Bio/Phylogeny/PhyBin/Util.hs
+++ b/Bio/Phylogeny/PhyBin/Util.hs
@@ -27,11 +27,10 @@
 import           System.FilePath     (combine)
 import           System.Directory    (doesFileExist, doesDirectoryExist,
                                       getDirectoryContents, getCurrentDirectory)
-import           System.IO           (openFile, hClose, IOMode(ReadMode))
-import           System.Process      (system)
+import           System.IO           (openFile, hClose, IOMode(ReadMode), stderr,
+                                      hPutStr, hPutStrLn)
 import           System.Exit         (ExitCode(..))
 import           Test.HUnit          ((~:),(~=?),Test,test)
-import qualified HSH 
 
 -- For vizualization:
 import           Text.PrettyPrint.HughesPJClass hiding (char, Style)
@@ -89,19 +88,20 @@
 
       --stat   <- if exists then getFileStatus path else return (error "internal invariant")
       -- [2010.09.23] This is no longer really necessary:
-      if not exists then do 
-	 putStr$ "Input not a file/directory, assuming wildcard, using 'find' for expansion"
-	 entries <- HSH.run$ "find " ++ path	 
-	 putStrLn$ "("++show (length entries)++" files found):  "++ show path
-	 return entries
+      if not exists then do
+        error$ "No file or directory found at this path!: "++path
+	 -- hPutStr stderr$ "Input not a file/directory, assuming wildcard, using 'find' for expansion"
+	 -- entries <- HSH.run$ "find " ++ path	 
+	 -- hPutStrLn stderr$ "("++show (length entries)++" files found):  "++ show path
+	 -- return entries
        else do
 	 isdir <- is_directory path
 	 reg  <- is_regular_file path
 	 if isdir then do 
-	    putStr$ "Input is a directory, reading all regular files contained "
+	    hPutStr stderr$ "Input is a directory, reading all regular files contained "
 	    children <- getDirectoryContents path
 	    filtered <- filterM is_regular_file $ map (combine path) children
-	    putStrLn$ "("++show (length filtered)++" regular files found):  "++ show path
+	    hPutStrLn stderr$ "("++show (length filtered)++" regular files found):  "++ show path
 	    return$ filtered
           else if reg then do 
 	    return [path]
diff --git a/Bio/Phylogeny/PhyBin/Visualize.hs b/Bio/Phylogeny/PhyBin/Visualize.hs
--- a/Bio/Phylogeny/PhyBin/Visualize.hs
+++ b/Bio/Phylogeny/PhyBin/Visualize.hs
@@ -10,8 +10,11 @@
        where
 import           Text.Printf        (printf)
 import           Data.List          (elemIndex, isPrefixOf)
-import           Data.Maybe         (fromJust)
-import           Data.Map           ((!))
+import           Data.List.Split    (chunksOf)
+import           Data.Maybe         (fromJust, isJust)
+import qualified Data.Map           as M
+import qualified Data.Set           as S
+import qualified Data.Vector        as V
 import           Data.Text.Lazy     (pack)
 import           Control.Monad      (void)
 import           Control.Concurrent  (Chan, newChan, writeChan, forkIO)
@@ -19,11 +22,13 @@
 import qualified Data.GraphViz        as Gv hiding (parse, toLabel)
 import qualified Data.GraphViz.Attributes.Complete as GA
 import qualified Data.GraphViz.Attributes.Colors   as GC
+import           Data.GraphViz.Attributes.Colors   (Color(RGB))
 -- import           Test.HUnit          ((~:),(~=?),Test,test)
 
 import qualified Data.Clustering.Hierarchical as C
 
 import           Bio.Phylogeny.PhyBin.CoreTypes
+import           Bio.Phylogeny.PhyBin.RFDistance (filterCompatible, compatibleWith, consensusTree)
 
 ----------------------------------------------------------------------------------------------------
 -- Visualization with GraphViz and FGL:
@@ -33,7 +38,7 @@
 toGraph :: FullTree StandardDecor -> G.Gr String Double
 toGraph (FullTree _ tbl tree) = G.run_ G.empty $ loop tree
   where
- fromLabel ix = tbl ! ix
+ fromLabel ix = tbl M.! ix
  loop (NTLeaf _ name) = 
     do let str = fromLabel name
        _ <- G.insMapNodeM str
@@ -59,24 +64,171 @@
        mapM_ (\x -> G.insMapEdgeM (node, x, branchLen$ get_dec x)) ls
        return ()
 
--- dendrogramToGraph :: C.Dendrogram (FullTree a) -> G.Gr (Either String String) Double
-dendrogramToGraph :: C.Dendrogram (FullTree a) -> G.Gr String Double
-dendrogramToGraph x = G.run_ G.empty $ void$ loop x
+----------------------------------------------------------------------------------------------------
+-- Dendrograms
+ ----------------------------------------------------------------------------------------------------
+
+-- | Some duplicated code with dotNewickTree.
+dotDendrogram :: PhyBinConfig -> String -> Double -> C.Dendrogram (FullTree a) ->
+                 Maybe (M.Map TreeName Int) -> [[NewickTree ()]] -> Gv.DotGraph G.Node
+dotDendrogram PBC{show_trees_in_dendro, show_interior_consensus}
+              title edge_scale origDendro mNameMap highlightTrs =
+  Gv.graphToDot myparams (G.nmap uid graph)
+ where
+  (charsDropped, dendro) = truncateNames origDendro
+  -- This is ugly, but we modify the name map to match:
+  nameMap' = fmap (M.mapKeys (drop charsDropped)) mNameMap  
+  graph :: DendroGraph
+  graph = dendrogramToGraph num_taxa dendro
+
+  num_taxa = numLeaves $ nwtree fstLeaf
+  FullTree{labelTable} = fstLeaf
+  fstLeaf = getFstLeaf dendro
+  getFstLeaf (C.Leaf ft)      = ft
+  getFstLeaf (C.Branch _ l _) = getFstLeaf l 
+
+  uidsToNames = M.fromList $
+                map (\nd@NdLabel{uid} -> (uid,nd)) $
+--                filter (isJust . tre) $ 
+                map (fromJust . G.lab graph) $  
+                G.nodes graph
+
+  matchers = map mkMatcher highlightTrs
+  mkMatcher ls = let fns = map compatibleWith ls -- Multiple predicate functions
+                 in \ tr -> -- Does it match any predicate?
+                    or$ map (\f -> f tr) fns 
+  
+  wcolors = zip matchers defaultPalette
+  findColor tr = loop wcolors
+    where loop [] = Nothing
+          loop ((f,col):rst) | f tr = Just col
+                             | otherwise = loop rst
+  
+  -- Map uids to highlight color
+  highlightMap = M.map fromJust $
+                 M.filter isJust $
+                 M.map (\ (NdLabel _ _ _ ct) -> findColor ct)
+                 uidsToNames
+
+  myparams :: Gv.GraphvizParams G.Node String Double () String -- (Either String String)
+  myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]],
+                                Gv.fmtNode= nodeAttrs,
+                                Gv.fmtEdge= edgeAttrs
+                              }
+  nodeAttrs :: (Int, UniqueNodeName) -> [GA.Attribute]
+  nodeAttrs (_num, uid) =
+    let uid' = if show_trees_in_dendro
+               then printed_tree++"\n"++uid
+               else uid
+        printed_tree =
+          case M.lookup uid uidsToNames of
+            Nothing -> ""
+            Just NdLabel{clumpSize,consensus} ->
+              (if clumpSize>1 then "size "++show clumpSize++"\n" else "")
+              ++ show (displayStrippedTree (FullTree "" labelTable consensus))
+        (tag,shp,styl) = -- case eith of
+          if isPrefixOf "DUMMY_" uid
+          then (if show_trees_in_dendro && show_interior_consensus
+                then printed_tree else "",
+                if show_interior_consensus
+                then GA.BoxShape -- GA.PlainText
+                else GA.PointShape,
+                [ GA.Color [weighted$ GA.X11Color Gv.Transparent] ]
+               )
+          else (uid', GA.Ellipse, [ GA.Style [GA.SItem GA.Filled []]])
+        highlightColor = 
+          case M.lookup uid highlightMap of
+            Nothing  -> []
+            Just col -> [ GA.Color [weighted col ] ]
+        clustColor | not (null highlightTrs) = []
+                   | otherwise = 
+          case (nameMap', M.lookup uid uidsToNames) of
+            (Just nm, Just NdLabel{tre=Just FullTree{treename}}) ->
+              case M.lookup treename nm of
+                Nothing -> []
+                -- Here we color the top TEN clusters:
+                Just ind | ind <= 10 -> [ GA.Color [weighted$ defaultPaletteV V.! (ind-1) ] ]
+                         | otherwise -> []
+            -- TODO: shall we color intermediate nodes?
+            _ -> []
+    in 
+    [ GA.Label$ GA.StrLabel$ pack tag
+    , GA.Shape shp
+    ] ++ styl ++ highlightColor ++ clustColor
+
+  edgeAttrs = getEdgeAttrs edge_scale
+
+
+type UniqueNodeName = String
+-- type DendroGraph = G.Gr (Maybe TreeName,UniqueNodeName) Double
+type DendroGraph = G.Gr NdLabel Double
+
+-- | When we first convert to a graph representation, there is a bunch of information
+-- hanging off of each node.
+data NdLabel =
+  NdLabel
+  { uid       :: UniqueNodeName
+  , tre       :: Maybe (FullTree ())
+  , clumpSize :: !Int
+  , consensus :: NewickTree ()
+  }
+  deriving (Show, Ord, Eq)
+
+-- | Create a graph using TreeNames for node labels and edit-distance for edge labels.
+dendrogramToGraph :: Int -> C.Dendrogram (FullTree a) -> DendroGraph
+dendrogramToGraph num_taxa orig =
+  G.run_ G.empty $ void$
+    loop (fmap (fmap (const())) orig)
   where
- -- deEither (Left s)  = s
- -- deEither (Right s) = s
- loop node@(C.Leaf FullTree{treename}) = G.insMapNodeM (treename)
+-- loop node@(C.Leaf ft) = G.insMapNodeM (NdLabel (treename ft) (Just$ fmap (const()) ft) 1 ft)
+ loop node@(C.Leaf ft) =
+   let stripped = fmap (const()) ft in
+   G.insMapNodeM (NdLabel (treename ft) (Just stripped) 1 (nwtree stripped))
+ loop node@(C.Branch 0 left right) = do
+   -- As a preprocessing step we collapse clusters that are separated by zero edit distance.
+   ----------------------------------------
+   let lvs = collapseZeroes left ++ collapseZeroes right
+       nms = map (treename) lvs
+       lens = map length nms
+       total = sum lens
+       avg   = total `quot` length nms
+       -- The goal here is to make an approximately square arrangement:
+       -- goal: avg * perline == total / (avg * perline)
+       perline = ceiling$ sqrt (fromIntegral total / ((fromIntegral avg)^2))
+       chunked = chunksOf perline nms
+       fatname = unlines (map unwords chunked)
+   G.insMapNodeM (NdLabel fatname (Just (head lvs))
+                  (length lvs) (consensusTree num_taxa (map nwtree lvs)))
+   ----------------------------------------   
  loop node@(C.Branch dist left right) =
-    do (_,l) <- loop left
-       (_,r) <- loop right
-       let ndname = "DUMMY_"++(l++"_"++r)
-                    -- Right (deEither l ++"_"++ deEither r)
-       (midN,mid) <- G.insMapNodeM ndname
-       G.insMapEdgeM (l, mid, dist)
-       G.insMapEdgeM (r, mid, dist)
+    do (_,ll@(NdLabel lid _ s1 c1)) <- loop left
+       (_,rr@(NdLabel rid _ s2 c2)) <- loop right
+       -- Interior nodes do NOT have their names drawn:
+       let ndname = "DUMMY_"++(lid++"_"++rid)  -- HACK!
+       (midN,mid) <- G.insMapNodeM (NdLabel ndname Nothing (s1+s2) (consensusTree num_taxa [c1,c2]))
+       G.insMapEdgeM (ll, mid, dist)
+       G.insMapEdgeM (rr, mid, dist)
        return (midN,mid)
 
+ collapseZeroes (C.Leaf tr)      = [tr]
+ collapseZeroes (C.Branch 0 l r) = collapseZeroes l ++ collapseZeroes r
+ collapseZeroes oth = error "dendrogramToGraph: internal error.  Not expecting non-zero branch length here."
 
+
+-- | The plot looks nicer when the names aren't bloated with repeated stuff. This
+-- replaces all tree names with potentially shorter names by removing the common prefix.
+-- Returns how many characters were dropped.
+truncateNames :: C.Dendrogram (FullTree a) -> (Int, C.Dendrogram (FullTree a))
+truncateNames dendro = (prefChars, fmap chopName dendro)
+ where 
+  chopName ft = ft{ treename= drop prefChars (treename ft) }
+  prefChars = length$ commonPrefix$ S.toList$ allNames dendro
+  allNames (C.Leaf tr)      = S.singleton (treename tr)
+  allNames (C.Branch _ l r) = S.union (allNames l) (allNames r)
+
+----------------------------------------------------------------------------------------------------
+
+
 -- | Open a GUI window to displaya tree.
 --
 --   Fork a thread that then runs graphviz.
@@ -121,7 +273,7 @@
     Gv.graphToDot myparams graph
  where 
   graph = toGraph2 atree
-  fromLabel ix = tbl ! ix  
+  fromLabel ix = tbl M.! ix  
   myparams :: Gv.GraphvizParams G.Node (NewickTree StandardDecor) Double () (NewickTree StandardDecor)
   myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]],
                                 Gv.fmtNode= nodeAttrs, Gv.fmtEdge= edgeAttrs }
@@ -137,6 +289,36 @@
   edgeAttrs = getEdgeAttrs edge_scale
 
 
+-- | Some arbitrarily chosen colors from the X11 set:
+defaultPaletteV :: V.Vector GA.Color
+defaultPaletteV = V.fromList defaultPalette
+
+defaultPalette :: [GA.Color]
+defaultPalette = concat$ replicate 4 $ map GA.X11Color 
+  [ Gv.Aquamarine
+  , Gv.PaleVioletRed
+  , Gv.MediumPurple
+  , Gv.PaleGreen
+  , Gv.PapayaWhip
+  , Gv.SkyBlue
+  , Gv.Yellow
+  , Gv.Crimson
+  , Gv.Gray
+  , Gv.PaleGoldenrod
+  ]
+
+-- Grabbed the first couple palettes off a website:
+altPalette :: V.Vector GA.Color
+altPalette = V.fromList $ concat $ replicate 3 $ 
+  -- http://www.colourlovers.com/palette/2962435/Autumn_Roses
+  [ RGB 159 74 81, RGB 217 183 173, RGB 149 91 116, RGB 185 138 148
+  --  , RGB 101 69 82 -- too dark
+  ] ++
+  -- http://www.colourlovers.com/palette/2962434/Earthy_warm
+  [ RGB 108 74 39, RGB 207 179 83, RGB 180 149 60, RGB 244 242 185
+  -- , RGB 61 63 39
+  ]
+  
 getEdgeAttrs :: Double -> (t, t1, Double) -> [GA.Attribute]
 getEdgeAttrs edge_scale = edgeAttrs
  where 
@@ -152,8 +334,6 @@
 			    then [GA.Color [weighted$ GA.X11Color Gv.Red],
                                   GA.LWidth 3.0, GA.Len minlen]
 			    else [GA.Len draw_weight]
-
-  weighted c = GC.WC {GC.wColor=c, GC.weighting=Nothing}
   minlen = 0.7
   maxlen = 3.0
   compute_draw_weight w scale = 
@@ -161,42 +341,14 @@
     -- Don't draw them too big or it gets annoying:
     (min scaled maxlen)
 
--- | Some duplicated code with dotNewickTree.
-dotDendrogram :: String -> Double -> C.Dendrogram (FullTree a) -> Gv.DotGraph G.Node
-dotDendrogram title edge_scale dendro =
-  Gv.graphToDot myparams graph
- where
---  graph :: G.Gr (Either String String) Double
-  graph :: G.Gr String Double   
-  graph = dendrogramToGraph dendro
-  myparams :: Gv.GraphvizParams G.Node String Double () String -- (Either String String)
-  myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]],
-                                Gv.fmtNode= nodeAttrs,
-                                Gv.fmtEdge= edgeAttrs
-                              }
---  nodeAttrs :: (Int, C.Dendrogram(FullTree StandardDecor)) -> [GA.Attribute]
-  nodeAttrs :: (Int, String) -> [GA.Attribute]
-  nodeAttrs (_num, eith) =
-    let (tag,shp) = -- case eith of
-          -- Left treename -> (take 60 treename, GA.Ellipse)
-          -- Right _       -> ("", GA.PointShape)
-          if isPrefixOf "DUMMY_" eith
-          then ("", GA.PointShape)          
-          else (take 60 eith, GA.Ellipse)
-    in 
-    [ GA.Label$ GA.StrLabel$ pack tag
-    , GA.Shape shp
-    , GA.Style [GA.SItem GA.Filled []]
-    ]
-  edgeAttrs = getEdgeAttrs edge_scale
-
+weighted c = GC.WC {GC.wColor=c, GC.weighting=Nothing}
 
 -- | This version shows the ordered/rooted structure of the normalized tree.
 dotNewickTree_debug :: String -> FullTree StandardDecor -> Gv.DotGraph G.Node
 dotNewickTree_debug title atree@(FullTree _ tbl tree) = Gv.graphToDot myparams graph
  where 
   graph = toGraph2 atree
-  fromLabel ix = tbl ! ix    
+  fromLabel ix = tbl M.! ix    
   myparams :: Gv.GraphvizParams G.Node (NewickTree StandardDecor) Double () (NewickTree StandardDecor)
   myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]],
 			        Gv.fmtNode= nodeAttrs, Gv.fmtEdge= edgeAttrs }
@@ -262,3 +414,13 @@
 -- TEMP / HACK:
 prettyPrint' :: Show a => a -> String
 prettyPrint' = show
+
+-- | Common prefix of a list of lists.
+commonPrefix :: Eq a => [[a]] -> [a]
+commonPrefix [] = []
+commonPrefix ls@(hd:tl)
+  | any null ls = []
+  | otherwise =
+    if   all ((== (head hd)) . head) tl
+    then head hd : commonPrefix (map tail ls)
+    else           commonPrefix (map tail ls)
diff --git a/Main.hs b/Main.hs
--- a/Main.hs
+++ b/Main.hs
@@ -1,8 +1,8 @@
-{-# LANGUAGE RecordWildCards, TupleSections #-}
+{-# LANGUAGE RecordWildCards, TupleSections, NamedFieldPuns #-}
 {-# OPTIONS_GHC -fwarn-unused-imports -fwarn-incomplete-patterns #-}
 
 module Main where
-import           Data.List (sort)
+import           Data.List (sort, intersperse, foldl')
 import qualified Data.ByteString.Lazy.Char8 as B
 import qualified Data.Map as M
 import qualified Data.Set as S
@@ -16,15 +16,19 @@
 import           Test.HUnit            (runTestTT, Test, test, (~:))
 
 import Control.Applicative ((<$>))
+import qualified Data.Vector                 as V
+import           Test.HUnit                  as HU
+
 import Data.GraphViz (runGraphvizCanvas,GraphvizCommand(Dot),GraphvizCanvas(Xlib))
 import Text.PrettyPrint.HughesPJClass hiding (char, Style)
 
 import Bio.Phylogeny.PhyBin.CoreTypes          
 import Bio.Phylogeny.PhyBin           (driver, binthem, normalize, annotateWLabLists,
-                                       unitTests, acquireTreeFiles, deAnnotate)
+                                       unitTests, acquireTreeFiles, deAnnotate,
+                                       retrieveHighlights, matchAnyHighlight)
 import Bio.Phylogeny.PhyBin.Parser    (parseNewick, parseNewicks, parseNewickFiles, unitTests)
 import Bio.Phylogeny.PhyBin.Visualize (viewNewickTree, dotNewickTree_debug)
-import Bio.Phylogeny.PhyBin.RFDistance (distanceMatrix, printDistMat, allBips)
+import Bio.Phylogeny.PhyBin.RFDistance 
 import Bio.Phylogeny.PhyBin.PreProcessor
 
 import qualified Data.Clustering.Hierarchical as C
@@ -49,8 +53,13 @@
     | View
     | TabDelimited Int Int
 
+    | Highlight FilePath
+    | ShowTreesInDendro | ShowInterior
+
+    | HashRF Bool
     | SelfTest
-    | RFMatrix | LineSetDiffMode | PrintNorms | PrintReg
+    | RFMatrix | LineSetDiffMode
+    | PrintNorms | PrintReg | PrintConsensus | PrintMatching
     | Cluster C.Linkage
     | BinningMode
     | EditDistThresh Int
@@ -77,29 +86,43 @@
        
 {- -- TODO: FIXME: IMPLEMENT THIS:
      , Option []        []          (NoArg NullOpt)  ""
-     , Option ['t']     ["tabbed"]  (ReqArg parseTabDelim "NUM1:NUM2")$  "assume the input is a tab-delimited file with gene names \n"++
+     , Option ['t']     ["tabbed"]  (ReqArg parseTabDelim "NUM1:NUM2")$  "assume the input is a ab-delimited file with gene names \n"++
 		                                                        "in column NUM1 and Newick trees in NUM2"
 -}       
      , Option []        []          (NoArg NullOpt)  ""
      , Option []        []  (NoArg$ error "internal problem")  "----------------------------- Clustering Options ------------------------------"
 
-     , Option []    ["bin"]      (NoArg BinningMode)$  "Use simple binning, the cheapest form of 'clustering'"
+     , Option []    ["bin"]      (NoArg BinningMode)$                 "Use simple binning, the cheapest form of 'clustering'"
      , Option []    ["single"]   (NoArg$ Cluster C.SingleLinkage)  $  "Use single-linkage clustering (nearest neighbor)"
      , Option []    ["complete"] (NoArg$ Cluster C.CompleteLinkage)$  "Use complete-linkage clustering (furthest neighbor)"
-     , Option []    ["UPGMA"]    (NoArg$ Cluster C.UPGMA)          $  "Use Unweighted Pair Group Method (average linkage)"
+     , Option []    ["UPGMA"]    (NoArg$ Cluster C.UPGMA)          $  "Use Unweighted Pair Group Method (average linkage) - DEFAULT mode"
 
      , Option []    ["editdist"]  (ReqArg (EditDistThresh . read) "DIST")$
-                                  "Combine all clusters separated by DIST or less.  Report a flat list of clusters."
-     , Option []    ["dendogram"] (NoArg DendogramOnly)$ "Report a hierarchical clustering (default)"
+                                  "Combine all clusters separated by DIST or less.  Report a flat list of clusters.\n"++
+                                  "Irrespective of whether this is activated, a hierarchical clustering (dendogram.pdf) is produced."  
+--     , Option []    ["dendogram"] (NoArg DendogramOnly)$ "Report a hierarchical clustering (default)"
+
+     , Option []        []     (NoArg$ error "internal problem")  "  Select Robinson-Foulds (symmetric difference) distance algorithm:"
+     , Option []    ["simple"] (NoArg$ HashRF False)
+       ((if hashRF then "" else "(default) ")++ "use direct all-to-all comparison")
+     , Option []    ["hashrf"] (NoArg$ HashRF True)
+       ((if hashRF then "(default) " else "")++"use a variant of the HashRF algorithm for the distance matrix")
        
      , Option []        []          (NoArg NullOpt)  ""
      , Option []        []  (NoArg$ error "internal problem")  "----------------------------- Visualization --------------------------------"
-     , Option ['g']     ["graphbins"] (NoArg Graph)  "use graphviz to produce .dot and .pdf output files named cluster1.*, cluster2.*, etc"
--- TODO: Produce the consensus tree as well as the individual trees.
-     , Option ['d']     ["drawbins"]  (NoArg Draw)   "like -g, but open GUI windows to show a tree for each bin"
+     , Option ['g']     ["graphbins"] (NoArg Graph)  "use graphviz to produce .dot and .pdf output files"
+      -- TODO: Produce the consensus tree as well as the individual trees.
+     , Option ['d']     ["drawbins"]  (NoArg Draw)   "like -g, but open GUI windows to show each bin's tree"
 
      , Option ['w']     ["view"]    (NoArg View)$  "for convenience, \"view mode\" simply displays input Newick files without binning" 
 
+     , Option [] ["showtrees"] (NoArg ShowTreesInDendro)  "Print (textual) tree topology inside the nodes of the dendrogram"
+     , Option [] ["highlight"] (ReqArg Highlight "FILE") $ 
+           "Highlight nodes in the tree-of-trees (dendrogram) consistent with the.\n"++
+           "given tree file.  Multiple highlights are permitted and use different colors."
+     , Option [] ["interior"] (NoArg ShowInterior)
+          "Show the consensus trees for interior nodes in the dendogram, rather than just points."
+           
      , Option []        []          (NoArg NullOpt)  ""
      , Option []        []  (NoArg$ error "internal problem")  "---------------------------- Tree pre-processing -----------------------------"
 
@@ -141,7 +164,11 @@
      , Option [] ["setdiff"] (NoArg LineSetDiffMode) "for convenience, print the set difference between cluster*.txt files"
      , Option [] ["print"]      (NoArg PrintReg)     "simply print out a concise form of each input tree"       
      , Option [] ["printnorms"] (NoArg PrintNorms)   "simply print out a concise and NORMALIZED form of each input tree"
+     , Option [] ["consensus"]  (NoArg PrintConsensus) "print a strict consensus tree for the inputs, then exit"
+     , Option [] ["matching"] (NoArg PrintMatching) "print a list of tree names that match any --highlight argument"
      ]
+ where
+   hashRF = use_hashrf default_phybin_config
 
 usage :: String
 usage = "\nUsage: phybin [OPTION...] files or directories...\n\n"++
@@ -181,47 +208,7 @@
 defaultErr errs = error $ "ERROR!\n" ++ (concat errs ++ usageInfo usage options)
 
 --------------------------------------------------------------------------------
--- Aggregated Unit Tests
---------------------------------------------------------------------------------
 
-allUnitTests :: Test
--- allUnitTests = unitTests ++
-allUnitTests = test 
-  [ Bio.Phylogeny.PhyBin.unitTests
-  , Bio.Phylogeny.PhyBin.Parser.unitTests
-  , "norm/Bip1" ~: (testNorm prob1)
-  ]
---  Bio.Phylogeny.PhyBin.Parser.unitTests
-
---    "(5_, (19, ((3_, 14), ((2_, 1_), (7_, (6_, 18))))), 13)"
--- [2013.07.23]      
--- This was INCORRECTLY normalizing to:
---     ((1_, 2_), (7_, (18, 6_)), ((14, 3_), (19, (13, 5_))))
-prob1 = "(5_, (19, ((3_, 14), ((2_, 1_), (7_, (6_, 18))))), 13);"
-
--- | Make sure that the normalized version of a tree yields the same bipartitions as
--- the unnormalized one.
-testNorm :: String -> IO ()
-testNorm str = do
-  let (labs,parsed) = parseNewick M.empty id "test" (B.pack str)
-      normed = normalize $ annotateWLabLists parsed
-      bips1  = allBips parsed
-      bips2  = allBips normed
-      added   = S.difference bips2 bips1
-      removed = S.difference bips1 bips2
-      dispBips bip = show$
-        map (map (labs M.!)) $ 
-        map IS.toList$ S.toList bip
-  unless (bips1 == bips2) $ do
-    putStrLn$ "Normalized this: "++show (displayDefaultTree $ FullTree "" labs parsed)
-    putStrLn$ "To this        : "++show (displayDefaultTree $ deAnnotate $ FullTree "" labs normed)
-    error$ "Normalization added and removed these bipartitions, respectively:\n  "
-           ++dispBips added ++"\n  "
-           ++dispBips removed
-
-
---------------------------------------------------------------------------------
-
 main :: IO ()
 main = 
   do argv <- getArgs 
@@ -241,6 +228,10 @@
 
      	   RFMatrix -> return cfg { print_rfmatrix= True }
 
+           ShowTreesInDendro -> return cfg { show_trees_in_dendro = True }
+           ShowInterior      -> return cfg { show_interior_consensus = True }
+           Highlight path    -> return cfg { highlights = path : highlights cfg }
+     
            LineSetDiffMode -> do
              bss <- mapM B.readFile files
              case map (S.fromList . B.lines) bss of
@@ -251,10 +242,13 @@
                oth -> error $"Line set difference mode expects two files as input, got "++show(length oth)
              exitSuccess
 
-           PrintNorms -> return cfg
-           PrintReg   -> return cfg           
-     
+           PrintNorms     -> return cfg
+           PrintReg       -> return cfg
+           PrintConsensus -> return cfg
+           PrintMatching  -> return cfg
+           
            Cluster lnk -> return cfg { clust_mode = ClusterThem lnk }
+           HashRF  bl  -> return cfg { use_hashrf = bl }
            BinningMode -> return cfg { clust_mode = BinThem }
            EditDistThresh n -> return cfg { dist_thresh = Just n }
            DendogramOnly    -> return cfg { dist_thresh = Nothing }
@@ -302,6 +296,25 @@
            liftFT (normalize . annotateWLabLists) ft
        exitSuccess
      ------------------------------------------------------------
+     when (elem PrintConsensus opts) $ do 
+       (_,fts) <- parseNewickFiles (name_hack config) all_inputs
+       putStrLn $ "Strict Consensus Tree of "++show (length fts)++" trees:"
+       when (null fts) $ error "No trees provided!"
+       let ctree = consensusTree (num_taxa config) (map nwtree fts)
+           FullTree{labelTable} = head fts
+       print$ displayStrippedTree$ FullTree "" labelTable ctree
+       exitSuccess       
+     ------------------------------------------------------------
+     when (elem PrintMatching opts) $ do 
+       (labs,fts) <- parseNewickFiles (name_hack config) all_inputs
+       when (null$ highlights config) $ error "No --highlight given, so --matching makes no sense.  Matching what?"
+       htrs <- retrieveHighlights (name_hack config) labs (highlights config)
+       let isMatch = matchAnyHighlight htrs
+       forM_ fts $ \ FullTree{treename,nwtree} ->
+         when (isMatch$ fmap (const()) nwtree) $
+          putStrLn treename
+       exitSuccess       
+     ------------------------------------------------------------
      when (View `elem` opts) $ do 
        view_graphs config
        exitSuccess
@@ -316,7 +329,7 @@
 		str <- B.readFile file
 		putStrLn$ "Parsed: " ++ (B.unpack str)
                 let (tbl,tr) = parseNewick M.empty name_hack file str
- 	        (chan, _tr) <- viewNewickTree file (FullTree "" tbl (annotateWLabLists tr))
+ 	        (chan, _tr) <- viewNewickTree file (FullTree "" tbl (annotateWLabLists (head tr)))
 	        return chan
 	      forM_ chans readChan 
 	      return ()
@@ -354,7 +367,50 @@
 temp = driver default_phybin_config{ num_taxa=7, inputs=["../datasets/test.tr"] }
 
 
+--------------------------------------------------------------------------------
+-- Aggregated Unit Tests
+--------------------------------------------------------------------------------
+
+allUnitTests :: Test
+-- allUnitTests = unitTests ++
+allUnitTests = test 
+  [ Bio.Phylogeny.PhyBin.unitTests
+  , Bio.Phylogeny.PhyBin.Parser.unitTests
+  , "norm/Bip1" ~: (testNorm prob1)
+  , "bipTreeConversion" ~: testBipConversion
+  , "t3" ~: t3_consensusTest, "t4" ~: t4_consensusTest, "t5" ~: t5_consensusTest
+  ]
+
+-- [2013.07.23]      
+-- This was INCORRECTLY normalizing to:
+--     ((1_, 2_), (7_, (18, 6_)), ((14, 3_), (19, (13, 5_))))
+prob1 :: String
+prob1 = "(5_, (19, ((3_, 14), ((2_, 1_), (7_, (6_, 18))))), 13);"
+
+-- | Make sure that the normalized version of a tree yields the same bipartitions as
+-- the unnormalized one.
+testNorm :: String -> IO ()
+testNorm str = do
+  let (labs,parseds) = parseNewick M.empty id "test" (B.pack str)
+      parsed = head parseds
+      normed = normalize $ annotateWLabLists parsed
+      bips1  = allBips parsed
+      bips2  = allBips normed
+      added   = S.difference bips2 bips1
+      removed = S.difference bips1 bips2
+      -- dispBips bip = show$
+      --   map (map (labs M.!)) $ 
+      --   map IS.toList$ S.toList bip
+  unless (bips1 == bips2) $ do
+    putStrLn$ "Normalized this: "++show (displayDefaultTree $ FullTree "" labs parsed)
+    putStrLn$ "To this        : "++show (displayDefaultTree $ deAnnotate $ FullTree "" labs normed)
+    error$ "Normalization added and removed these bipartitions, respectively:\n  "
+           ++concat (intersperse " " (map (dispBip labs) (S.toList added))) ++"\n  "
+           ++concat (intersperse " " (map (dispBip labs) (S.toList removed)))
+
+
 -- 112 and 13
+rftest :: IO ()
 rftest = do 
   (mp,[t1,t2]) <- parseNewickFiles (take 2) ["tests/13.tr", "tests/112.tr"]
   putStrLn$ "Tree 13           : " ++ show (displayDefaultTree t1)
@@ -369,7 +425,7 @@
   putStrLn$ "13  collapsed 0.03: " ++show (disp$ liftFT (collapseBranchLenThresh 0.03) t1)
   putStrLn$ "112 collapsed 0.03: " ++show (disp$ liftFT (collapseBranchLenThresh 0.03) t2)  
 
-  let mat = distanceMatrix [nwtree t1, nwtree t2]
+  let (mat,_) = distanceMatrix [nwtree t1, nwtree t2]
   printDistMat stdout mat
   return ()
  where
@@ -377,6 +433,83 @@
     let collapsed :: AnnotatedTree 
         collapsed = normalize$ annotateWLabLists tr
     in displayDefaultTree$ deAnnotate $  FullTree nm labs collapsed
+
+-- | This test was done with --editdist 4 --complete
+t3_consensusTest :: IO ()
+t3_consensusTest = consensusTest "./tests/t3_consensus/cluster1_284_alltrees.tr"
+                                 "./tests/t3_consensus/cluster1_284_consensus.tr"
+
+-- | This test was done with --editdist 0 --complete
+t4_consensusTest :: IO ()
+t4_consensusTest = consensusTest "./tests/t4_consensus/cluster1_16_alltrees.tr"
+                                 "./tests/t4_consensus/cluster1_16_consensus.tr"
+
+-- | This test was done with --editdist 1 --complete
+t5_consensusTest :: IO ()
+t5_consensusTest = consensusTest "./tests/t5_consensus/cluster1_35_alltrees.tr"
+                                 "./tests/t5_consensus/cluster1_35_consensus.tr"
+
+consensusTest :: String -> String -> IO ()
+consensusTest alltrees consensus = do  
+  (_,ctree:ftrees)  <- parseNewickFiles id [consensus,alltrees]
+  let num_taxa      = numLeaves (nwtree ctree)
+      plainTrs      = map nwtree  ftrees 
+      eachbips      = map allBips plainTrs
+      totalBips     = foldl' S.union        S.empty eachbips
+      intersectBips = foldl' S.intersection S.empty eachbips
+      FullTree _ labs _ = ctree
+      linesPrnt x = unlines (map (("  "++) . dispBip labs) $ S.toList x)
+  putStrLn$ "Bips in each: "++     show (map S.size eachbips)
+  putStrLn$ "Total bips in all: "++show (S.size totalBips)  
+  putStrLn$ "Bips in common: "++   show (S.size intersectBips)
+  putStrLn$ "Bips of first member:\n" ++ linesPrnt (head eachbips)
+  putStrLn$ "Some bips in the union that are NOT in the intersection:\n" ++
+     linesPrnt (S.fromList$ take 20$ S.toList$ S.difference totalBips intersectBips)
+  let cbips = allBips $ nwtree ctree
+  putStrLn$ "ConsensusBips ("++show (S.size cbips)++"):\n"++linesPrnt cbips
+  putStrLn$"Things in the consensus that should NOT be:\n"++linesPrnt (S.difference cbips intersectBips)
+  putStrLn$"Things not in the consensus that SHOULD be:\n"++linesPrnt (S.difference intersectBips cbips)
+
+  putStrLn$ "Now recomputing consensus tree for "++show num_taxa++" taxa"
+  let ctree2 = consensusTree num_taxa plainTrs
+      cbips2 = allBips ctree2
+  putStrLn$ "Freshly recomputed consensusBips ("++show (S.size cbips2)++"):\n"++linesPrnt cbips2
+  HU.assertEqual "Consensus tree on disk should match computed one:"
+         cbips cbips2 -- (allBips$ fmap (const ()) $ nwtree ctree)         
+  
+  putStrLn " Partial distance matrix WITHIN this cluster:"
+  let (mat,_) = distanceMatrix (map nwtree ftrees)
+  printDistMat stdout (V.take 30 mat)
+  HU.assertBool "Consensus should only include bicbips2ps in the members" (S.isSubsetOf cbips totalBips)
+  HU.assertEqual "Consensus tree matches intersected bips" cbips intersectBips 
+  return ()
+
+testBipConversion :: IO ()
+testBipConversion = 
+  do (_,trs) <- allTestTrees
+     mapM_ testOne trs
+     putStrLn "All trees passed bip conversion."
+  where
+    testOne (FullTree{nwtree}) = do
+      let sz    = numLeaves nwtree   
+          bips  = allBips nwtree
+          tr2   = bipsToTree sz bips
+          bips2 = allBips tr2
+      assertEqual "round trip bips->tree->bips" bips bips2
+
+-- | Read in all test trees which we happen to have put in the repository for testing
+-- purposes.
+allTestTrees :: IO (LabelTable, [FullTree DefDecor])
+allTestTrees =
+  parseNewickFiles id $
+  [ "./tests/t3_consensus/cluster1_284_alltrees.tr"
+  , "./tests/t3_consensus/cluster1_284_consensus.tr"
+  , "./tests/t4_consensus/cluster1_16_alltrees.tr"
+  , "./tests/t4_consensus/cluster1_16_consensus.tr"
+  , "./tests/t5_consensus/cluster1_35_alltrees.tr"
+  , "./tests/t5_consensus/cluster1_35_consensus.tr"
+  ]
+  
 
 ----------------------------------------------------------------------------------------------------
 -- TODO: expose a command line argument for testing.
diff --git a/phybin.cabal b/phybin.cabal
--- a/phybin.cabal
+++ b/phybin.cabal
@@ -1,5 +1,5 @@
 Name:           phybin
-Version:        0.2.2
+Version:        0.2.11
 License: BSD3
 License-file:   LICENSE
 Stability: Beta
@@ -14,9 +14,17 @@
 -- 0.1.2.5 -- bump for graphviz API changes
 -- 0.2 -- Add Robinson-Foulds distance, use Int labels.
 -- 0.2.2 -- misc changes and expose library
+-- 0.2.4 -- add consensus trees
+-- 0.2.6 -- colorization, hashrf, misc improvements
+-- 0.2.7 -- Add command line opt --showtrees
+-- 0.2.8 -- Add command line opt --highlight
+-- 0.2.9 -- Add command line opt --interior
+-- 0.2.10 -- Add command line opt --matching
+-- 0.2.11 -- Cleanup and windows compatibility.
 
 -- homepage: http://code.haskell.org/phybin
-homepage: http://people.csail.mit.edu/newton/phybin/
+-- homepage: http://people.csail.mit.edu/newton/phybin/
+homepage: http://www.cs.indiana.edu/~rrnewton/projects/phybin/
 
 Copyright: Copyright (c) 2010 Ryan Newton
 Synopsis: Utility for clustering phylogenetic trees in Newick format based on Robinson-Foulds distance.
@@ -26,7 +34,7 @@
    It can do simple binning of identical trees or more complex clustering based on 
    an all-to-all Robinson-Foulds distance matrix.
    .
-   Output files that characterize the size and contents of each bin or cluster (including
+   phybin produces output files that characterize the size and contents of each bin or cluster (including
    generating GraphViz-based visual representations of the tree topologies).
 
 Category: Bioinformatics
@@ -40,6 +48,18 @@
   type:     git
   location: git://github.com/rrnewton/PhyBin.git
 
+Flag hashrf
+  description: Use the HashRF algorithm by default instead of the naive one.
+  default: True
+
+Flag bitvec
+  description: Use bitvectors rather than IntSets for bipartitions.
+  default: False
+
+Flag sequential
+  description: Don't use any parallelism at all.
+  default: False
+
 Library
   Exposed-modules:   Bio.Phylogeny.PhyBin
                      Bio.Phylogeny.PhyBin.RFDistance
@@ -49,34 +69,50 @@
                      Bio.Phylogeny.PhyBin.Util
                      Bio.Phylogeny.PhyBin.Visualize
 --  Other-modules:     Data.BitList
-  Build-Depends:     base >= 3 && < 5, directory, process, containers, unix, 
+  Build-Depends:     base >= 3 && < 5, directory, process, containers, 
                      async, time,
                      filepath, 
                      graphviz >= 2999.16, 
                      text >= 0.11 && < 0.12,
                      prettyclass, fgl,
-                     HSH, HUnit, bytestring, 
+                     HUnit, bytestring, 
                      -- For bytestring.lazy support:
                      parsec >= 3.1.0, 
-                     bitvec >= 0.1, vector >= 0.10,
-                     hierarchical-clustering >= 0.4
+                     vector >= 0.10,
+                     hierarchical-clustering >= 0.4, split >= 0.2
 --                     lattice-par, 
   GHC-Options: -O2 -funbox-strict-fields -rtsopts
 
+  if flag(hashrf)
+     CPP-Options: -DUSE_HASHRF
+  if flag(bitvec)
+     CPP-Options: -DBITVEC_BIPS
+     Build-Depends:  bitvec >= 0.1
+  if flag(sequential)
+     CPP-Options: -DSEQUENTIALIZE
+
 Executable phybin
   Main-is:           Main.hs
   Build-Depends:     phybin
   -- DUPLICATE:
-  Build-Depends:     base >= 3 && < 5, directory, process, containers, unix, 
+  Build-Depends:     base >= 3 && < 5, directory, process, containers, 
                      async, time,
                      filepath, 
                      graphviz >= 2999.16, 
                      text >= 0.11 && < 0.12,
                      prettyclass, fgl,
-                     HSH, HUnit, bytestring, 
+                     HUnit, bytestring, 
                      -- For bytestring.lazy support:
                      parsec >= 3.1.0, 
-                     bitvec >= 0.1, vector >= 0.10,
-                     hierarchical-clustering >= 0.4
+                     vector >= 0.10,
+                     hierarchical-clustering >= 0.4, split >= 0.2
 --                     lattice-par, 
   GHC-Options: -O2 -funbox-strict-fields -rtsopts
+
+  if flag(hashrf)
+     CPP-Options: -DUSE_HASHRF
+  if flag(bitvec)
+     CPP-Options: -DBITVEC_BIPS
+     Build-Depends:  bitvec >= 0.1
+  if flag(sequential)
+     CPP-Options: -DSEQUENTIALIZE