diff --git a/Dvda.hs b/Dvda.hs
--- a/Dvda.hs
+++ b/Dvda.hs
@@ -16,22 +16,20 @@
               -- * symbolic expression type
             , Expr
               -- * construct FunGraphs
-            , toFunGraph
-            , cse
+--            , cse
               -- * show/summarize FunGraphs
-            , previewGraph
-            , previewGraph'
+--            , previewGraph
+--            , previewGraph'
               -- * compile and link function
 --            , buildHSFunction
 --            , buildHSFunctionPure
 --            , buildHSFunctionFromGraph
-              -- * Heterogenous inputs/outputs
-            , (:*)(..)
+            , module Dvda.Algorithm
             ) where
 
 import Dvda.AD ( rad )
-import Dvda.CSE ( cse )
+--import Dvda.CSE ( cse )
 import Dvda.Expr ( Expr, sym, symDependent, symDependentN )
-import Dvda.FunGraph ( toFunGraph, (:*)(..) )
-import Dvda.Vis ( previewGraph, previewGraph' )
+import Dvda.Algorithm
+--import Dvda.Vis ( previewGraph, previewGraph' )
 --import Dvda.HSBuilder
diff --git a/Dvda/AD.hs b/Dvda/AD.hs
--- a/Dvda/AD.hs
+++ b/Dvda/AD.hs
@@ -34,7 +34,7 @@
     dfdg = dualPerturbation $ unop (Dual g 1)
 
 backpropNode :: (Ord a, Num a) => Expr a -> Expr a -> [(Expr a, Expr a)]
-backpropNode sens e@(ESym (SymDependent name k dep_)) = (e,sens):(backpropNode (sens*primal') dep)
+backpropNode sens e@(ESym (SymDependent name k dep_)) = (e,sens):backpropNode (sens*primal') dep
   where
     primal' = ESym (SymDependent name (k+1) dep_)
     dep = ESym dep_
@@ -55,6 +55,7 @@
 backpropNode sens (EFloating (Log x))   = bpUnary sens x log
 backpropNode sens (EFloating (Sin x))   = bpUnary sens x sin
 backpropNode sens (EFloating (Cos x))   = bpUnary sens x cos
+backpropNode sens (EFloating (Tan x))   = bpUnary sens x tan
 backpropNode sens (EFloating (ASin x))  = bpUnary sens x asin
 backpropNode sens (EFloating (ATan x))  = bpUnary sens x atan
 backpropNode sens (EFloating (ACos x))  = bpUnary sens x acos
@@ -68,7 +69,7 @@
 backprop :: (Num a, Ord a, Hashable a) => Expr a -> HashMap (Expr a) (Expr a)
 backprop x = HM.fromListWith (+) (backpropNode 1 x)
 
-rad :: (Num a, Ord a, Hashable a) => Expr a -> [Expr a] -> [Expr a]
-rad x args = map (\arg -> HM.lookupDefault 0 arg sensitivities) args
+rad :: (Num a, Ord a, Hashable a, Functor f) => Expr a -> f (Expr a) -> f (Expr a)
+rad x = fmap (\arg -> HM.lookupDefault 0 arg sensitivities)
   where
     sensitivities = backprop x
diff --git a/Dvda/Algorithm.hs b/Dvda/Algorithm.hs
new file mode 100644
--- /dev/null
+++ b/Dvda/Algorithm.hs
@@ -0,0 +1,37 @@
+{-# OPTIONS_GHC -Wall #-}
+
+module Dvda.Algorithm
+       ( Algorithm
+       , constructAlgorithm
+       , runAlgorithm
+       , runAlgorithm'
+       , toSymbolicAlg
+       , squashIsSame
+       ) where
+
+import qualified Data.Vector.Generic as G
+
+import Dvda.Algorithm.Construct ( Algorithm(..), AlgOp(..), constructAlgorithm, squashWorkVector )
+import Dvda.Algorithm.Eval ( runAlgorithm, runAlgorithm' )
+import Dvda.Expr ( Expr(..), GExpr(..) )
+
+-- | test to see of SSA algorithm works the same as Live variables version
+squashIsSame :: (Eq (v a), G.Vector v a) => v a -> Algorithm a -> Bool
+squashIsSame x alg = runAlgorithm alg x == runAlgorithm (squashWorkVector alg) x
+
+-- | Convert an algorithm into a symbolic algorithm.
+--   Is there any reason to keep this when we have toFloatingAlg?
+toSymbolicAlg :: Eq a => Algorithm a -> Algorithm (Expr a)
+toSymbolicAlg (Algorithm ind outd ops ws) = Algorithm ind outd (map opToExpr ops) ws
+  where
+    opToExpr :: Eq a => AlgOp a -> AlgOp (Expr a)
+    opToExpr (InputOp k idx) = InputOp k idx
+    opToExpr (OutputOp k idx) = OutputOp k idx
+    opToExpr (NormalOp k gexpr) = NormalOp k (g2e gexpr)
+      where
+        g2e :: Eq a => GExpr a b -> GExpr (Expr a) b
+        g2e (GSym x) = GSym x
+        g2e (GConst x) = GConst (EConst x)
+        g2e (GNum x) = GNum x
+        g2e (GFractional x) = GFractional x
+        g2e (GFloating x) = GFloating x
diff --git a/Dvda/Algorithm/Construct.hs b/Dvda/Algorithm/Construct.hs
new file mode 100644
--- /dev/null
+++ b/Dvda/Algorithm/Construct.hs
@@ -0,0 +1,153 @@
+{-# OPTIONS_GHC -Wall #-}
+
+module Dvda.Algorithm.Construct
+       ( Algorithm(..)
+       , AlgOp(..)
+       , Node(..)
+       , InputIdx(..)
+       , OutputIdx(..)
+       , constructAlgorithm
+       , squashWorkVector
+       ) where
+
+import qualified Data.Foldable as F
+import Data.Maybe ( fromMaybe )
+import qualified Data.Traversable as T
+import qualified Data.IntMap as IM
+import qualified Data.Vector as V
+import qualified Data.HashMap.Lazy as HM
+
+import Dvda.Expr
+import Dvda.Algorithm.FunGraph ( FunGraph(..), Node(..), toFunGraph )
+
+newtype InputIdx = InputIdx Int deriving Show
+newtype OutputIdx = OutputIdx Int deriving Show
+
+data AlgOp a = InputOp  Node InputIdx
+             | OutputOp Node OutputIdx
+             | NormalOp Node (GExpr a Node)
+             deriving Show
+
+data Algorithm a = Algorithm { algInDims :: Int
+                             , algOutDims :: Int
+                             , algOps :: [AlgOp a]
+                             , algWorkSize :: Int
+                             }
+
+newtype LiveNode = LiveNode Int deriving Show
+newtype NodeMap a = NodeMap (IM.IntMap a) deriving Show
+nmEmpty :: NodeMap a
+nmEmpty = NodeMap IM.empty
+
+nmInsertWith :: (a -> a -> a) -> Node -> a -> NodeMap a -> NodeMap a
+nmInsertWith f (Node k) v (NodeMap im) = NodeMap (IM.insertWith f k v im)
+
+nmLookup :: Node -> NodeMap a -> Maybe a
+nmLookup (Node k) (NodeMap im) = IM.lookup k im
+
+nmInsert :: Node -> a -> NodeMap a -> NodeMap a
+nmInsert (Node k) val (NodeMap im) = NodeMap $ IM.insert k val im
+
+squashWorkVec' :: NodeMap Int -> NodeMap LiveNode -> [LiveNode] -> [AlgOp a] -> [AlgOp a]
+squashWorkVec' accessMap liveMap0 (LiveNode pool0:pools) (InputOp k inIdx:xs) =
+  InputOp (Node pool0) inIdx : squashWorkVec' accessMap liveMap pools xs
+  where
+    -- input to the the first element of the live pool
+    -- update the liveMap to reflect this
+    -- update the pool to reflect this
+    liveMap = nmInsertWith err k (LiveNode pool0) liveMap0
+    err = error "SSA node written to more than once"
+squashWorkVec' accessMap0 liveMap0 pool0 (OutputOp k outIdx:xs) =
+  OutputOp (Node lk) outIdx : squashWorkVec' accessMap liveMap0 pool xs
+  where
+    -- output from node looked up from live variables
+    (LiveNode lk) = fromMaybe noLiveErr (nmLookup k liveMap0)
+      where noLiveErr = error "OutputOp couldn't find node in live map"
+    -- decrement access map, if references are now zero, add live node back to pool
+    (accessMap, pool) = case nmLookup k accessMap0 of
+      Just 0 -> error "squashWorkVec': accessed something with 0 references"
+      Just 1 -> (nmInsert k 0 accessMap0, LiveNode lk:pool0)
+      Just n -> (nmInsert k (n-1) accessMap0, pool0)
+      Nothing -> error "squashWorkVec': node not in access map"
+squashWorkVec' accessMap0 liveMap0 pool0 (NormalOp k gexpr0:xs) =
+  NormalOp (Node retLiveK) gexpr : squashWorkVec' accessMap liveMap pool xs
+  where
+    decrement (am0, p0) depk = case nmLookup depk am0 of
+      Just 0 -> error "squashWorkVec': accessed something with 0 references"
+      Just 1 -> ((nmInsert depk     0 am0, LiveNode lk:p0), Node lk)
+      Just n -> ((nmInsert depk (n-1) am0,             p0), Node lk)
+      Nothing -> error "squashWorkVec': node not in access map"
+      where
+        LiveNode lk = fromMaybe (error "depsLiveKs missing") (nmLookup depk liveMap0)
+
+    ((accessMap, LiveNode retLiveK:pool), gexpr) =
+      T.mapAccumL decrement (accessMap0, pool0) gexpr0
+    liveMap = nmInsert k (LiveNode retLiveK) liveMap0
+squashWorkVec' _ _ _ [] = []
+squashWorkVec' _ _ [] _ = error "squashWorkVec': empty pool"
+
+-- | Converts SSA to live variables.
+--   This reduces the size of the work vector by re-using dead registers.
+--   Does this break if it's called more than once?
+--   Maybe these should have different types
+squashWorkVector :: Algorithm a -> Algorithm a
+squashWorkVector alg =
+  Algorithm { algOps = newAlgOps
+            , algInDims = algInDims alg
+            , algOutDims = algOutDims alg
+            , algWorkSize = workVectorSize newAlgOps
+            }
+  where
+    addOne k = nmInsertWith (+) k (1::Int)
+    countAccesses accMap  (InputOp _ _:xs) = countAccesses accMap xs
+    countAccesses accMap  (OutputOp k _:xs) = countAccesses (addOne k accMap) xs
+    countAccesses accMap0 (NormalOp _ gexpr:xs) = countAccesses accMap xs
+      where
+        accMap = F.foldr addOne accMap0 gexpr
+    countAccesses accMap [] = accMap
+    accesses = countAccesses nmEmpty (algOps alg)
+
+    newAlgOps = squashWorkVec' accesses nmEmpty (map LiveNode [0..]) (algOps alg)
+
+
+graphToAlg :: [(Node,GExpr a Node)] -> V.Vector (Sym,InputIdx) -> V.Vector (Node,OutputIdx)
+              -> [AlgOp a]
+graphToAlg rgr0 inSyms outIdxs = f rgr0
+  where
+    inSymMap = HM.fromList (F.toList inSyms :: [(Sym,InputIdx)])
+    outIdxMap = IM.fromList (map (\(Node k, x) -> (k, x)) (F.toList outIdxs :: [(Node,OutputIdx)]))
+
+    f ((k@(Node k'),GSym s):xs) = case HM.lookup s inSymMap of
+      Nothing -> error "toAlg: symbolic is not in inputs"
+      Just inIdx -> case IM.lookup k' outIdxMap of
+        -- sym is an input
+        Nothing -> InputOp k inIdx : f xs
+        -- sym is an input and an output
+        Just outIdx -> InputOp k inIdx : OutputOp k outIdx : f xs
+    f ((k@(Node k'),x):xs) = case IM.lookup k' outIdxMap of
+      -- no input or output
+      Nothing -> NormalOp k x : f xs
+      -- output only
+      Just outIdx -> NormalOp k x : OutputOp k outIdx : f xs
+    f [] = []
+
+workVectorSize :: [AlgOp a] -> Int
+workVectorSize = workVectorSize' (-1)
+  where
+    workVectorSize' n (NormalOp (Node m) _:xs) = workVectorSize' (max n m) xs
+    workVectorSize' n (InputOp  (Node m) _:xs) = workVectorSize' (max n m) xs
+    workVectorSize' n (OutputOp (Node m) _:xs) = workVectorSize' (max n m) xs
+    workVectorSize' n [] = n+1
+
+-- | create a SSA algorithm from a vector of symbolic inputs and outputs
+constructAlgorithm :: V.Vector (Expr a) -> V.Vector (Expr a) -> IO (Algorithm a)
+constructAlgorithm inputVecs outputVecs = do
+  fg <- toFunGraph inputVecs outputVecs
+  let inputIdxs  = V.map (\(k,x) -> (x,  InputIdx k)) (V.indexed ( fgInputs fg))
+      outputIdxs = V.map (\(k,x) -> (x, OutputIdx k)) (V.indexed (fgOutputs fg))
+      ops = graphToAlg (fgReified fg) inputIdxs outputIdxs
+  return Algorithm { algInDims  = V.length inputIdxs
+                   , algOutDims = V.length outputIdxs
+                   , algOps = ops
+                   , algWorkSize = workVectorSize ops
+                   }
diff --git a/Dvda/Algorithm/Eval.hs b/Dvda/Algorithm/Eval.hs
new file mode 100644
--- /dev/null
+++ b/Dvda/Algorithm/Eval.hs
@@ -0,0 +1,99 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# Language Rank2Types #-}
+{-# Language FlexibleContexts #-}
+
+module Dvda.Algorithm.Eval
+       ( runAlgorithm
+       , runAlgorithm'
+       ) where
+
+import Control.Monad.ST ( ST, runST )
+import Data.Vector.Generic ( (!) )
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as GM
+
+import Dvda.Expr
+import Dvda.Algorithm.Construct ( Algorithm(..), AlgOp(..), InputIdx(..), OutputIdx(..) )
+import Dvda.Algorithm.FunGraph ( Node(..) )
+
+newtype RtOp v a = RtOp (forall s. (G.Mutable v) s a -> v a -> (G.Mutable v) s a -> ST s ())
+
+-- | purely run an algoritm
+runAlgorithm :: G.Vector v a => Algorithm a -> v a -> Either String (v a)
+runAlgorithm alg =
+  runAlg'' (algInDims alg) (algOutDims alg) (algWorkSize alg) (map toRtOp (algOps alg))
+  where
+    runAlg'' :: G.Vector v a => Int -> Int -> Int -> [RtOp v a] -> v a -> Either String (v a)
+    runAlg'' inSize outSize workSize ops inputVec
+      | G.length inputVec /= inSize =
+        Left $ "runAlg: input dimension mismatch, given: " ++ show (G.length inputVec) ++
+               ", expected: " ++ show inSize
+      | otherwise = Right $ runST $ do
+        workVec <- GM.new workSize
+        outputVec <- GM.new outSize
+        mapM_ (\(RtOp op) -> op workVec inputVec outputVec) ops
+        G.freeze outputVec
+
+-- | run an algoritm in the ST monad, mutating a user-provided output vector
+runAlgorithm' :: G.Vector v a => Algorithm a -> v a -> G.Mutable v s a -> ST s (Maybe String)
+runAlgorithm' alg =
+  runAlg'' (algInDims alg) (algOutDims alg) (algWorkSize alg) (map toRtOp (algOps alg))
+  where
+    runAlg'' :: G.Vector v a => Int -> Int -> Int -> [RtOp v a] -> v a -> G.Mutable v s a -> ST s (Maybe String)
+    runAlg'' inSize outSize workSize ops inputVec outputVec
+      | G.length inputVec /= inSize =
+        return $ Just $ "runAlg': input dimension mismatch, given: " ++ show (G.length inputVec) ++
+                        ", expected: " ++ show inSize
+      | GM.length outputVec /= outSize =
+        return $ Just $ "runAlg': output dimension mismatch, given: " ++ show (GM.length outputVec) ++
+                        ", expected: " ++ show outSize
+      | otherwise = do
+        workVec <- GM.new workSize
+        mapM_ (\(RtOp op) -> op workVec inputVec outputVec) ops
+        return Nothing
+
+bin :: GM.MVector (G.Mutable v) a => Node -> Node -> Node -> (a -> a -> a) -> RtOp v a
+bin (Node k) (Node kx) (Node ky) f = RtOp $ \work _ _ -> do
+  x <- GM.read work kx
+  y <- GM.read work ky
+  GM.write work k (f x y)
+
+un :: GM.MVector (G.Mutable v) a => Node -> Node -> (a -> a) -> RtOp v a
+un (Node k) (Node kx) f = RtOp $ \work _ _ -> GM.read work kx >>= GM.write work k . f
+
+toRtOp :: G.Vector v a => AlgOp a -> RtOp v a
+toRtOp (InputOp (Node k) (InputIdx i)) = RtOp $ \work input _ -> GM.write work k (input ! i)
+toRtOp (OutputOp (Node k) (OutputIdx i)) = RtOp $ \work _ output -> do
+  GM.read work k >>= GM.write output i
+toRtOp (NormalOp (Node k) (GConst c)) =
+  RtOp $ \work _ _ -> GM.write work k c
+toRtOp (NormalOp (Node k) (GNum (FromInteger x))) =
+  RtOp $ \work _ _ -> GM.write work k (fromIntegral x)
+toRtOp (NormalOp (Node k) (GFractional (FromRational x))) =
+  RtOp $ \work _ _ -> GM.write work k (fromRational x)
+
+toRtOp (NormalOp k (GNum (Mul x y)))  = bin k x y (*)
+toRtOp (NormalOp k (GNum (Add x y)))  = bin k x y (+)
+toRtOp (NormalOp k (GNum (Sub x y)))  = bin k x y (-)
+toRtOp (NormalOp k (GNum (Negate x))) = un k x negate
+toRtOp (NormalOp k (GFractional (Div x y)))   = bin k x y (/)
+
+toRtOp (NormalOp k (GNum (Abs x)))            = un k x abs
+toRtOp (NormalOp k (GNum (Signum x)))         = un k x signum
+toRtOp (NormalOp k (GFloating (Pow x y)))     = bin k x y (**)
+toRtOp (NormalOp k (GFloating (LogBase x y))) = bin k x y logBase
+toRtOp (NormalOp k (GFloating (Exp x)))       = un k x exp
+toRtOp (NormalOp k (GFloating (Log x)))       = un k x log
+toRtOp (NormalOp k (GFloating (Sin x)))       = un k x sin
+toRtOp (NormalOp k (GFloating (Cos x)))       = un k x cos
+toRtOp (NormalOp k (GFloating (Tan x)))       = un k x tan
+toRtOp (NormalOp k (GFloating (ASin x)))      = un k x asin
+toRtOp (NormalOp k (GFloating (ATan x)))      = un k x atan
+toRtOp (NormalOp k (GFloating (ACos x)))      = un k x acos
+toRtOp (NormalOp k (GFloating (Sinh x)))      = un k x sinh
+toRtOp (NormalOp k (GFloating (Cosh x)))      = un k x cosh
+toRtOp (NormalOp k (GFloating (Tanh x)))      = un k x tanh
+toRtOp (NormalOp k (GFloating (ASinh x)))     = un k x asinh
+toRtOp (NormalOp k (GFloating (ATanh x)))     = un k x atanh
+toRtOp (NormalOp k (GFloating (ACosh x)))     = un k x acosh
+toRtOp (NormalOp _ (GSym _)) = error "runAlg: there's symbol in my algorithm"
diff --git a/Dvda/Algorithm/FunGraph.hs b/Dvda/Algorithm/FunGraph.hs
new file mode 100644
--- /dev/null
+++ b/Dvda/Algorithm/FunGraph.hs
@@ -0,0 +1,84 @@
+{-# OPTIONS_GHC -Wall #-}
+
+module Dvda.Algorithm.FunGraph
+       ( FunGraph(..)
+       , Node(..)
+       , toFunGraph
+       ) where
+
+import Control.Applicative ( (<$>) )
+import Data.Foldable ( Foldable )
+import qualified Data.Foldable as F
+import qualified Data.Graph as Graph
+import qualified Data.HashSet as HS
+import Data.Traversable ( Traversable )
+
+import Dvda.Expr
+import Dvda.Algorithm.Reify ( ReifyGraph(..), Node(..), reifyGraph )
+
+data FunGraph f g a = FunGraph { fgInputs :: f Sym
+                               , fgOutputs :: g Node
+                               , fgReified :: [(Node, GExpr a Node)]
+                               , fgTopSort :: [(Node, GExpr a Node)]
+                               }
+
+-- | find any symbols which are parents of outputs, but are not supplied by the user
+detectMissingInputs :: Foldable f => f (Expr a) -> [(Node, GExpr a Node)] -> [Sym]
+detectMissingInputs exprs gr = HS.toList $ HS.difference allGraphInputs allUserInputs
+  where
+    allUserInputs =
+      let f (ESym name) acc = name : acc
+          f _ _ = error $ "detectMissingInputs given non-ESym input" -- \"" ++ show e ++ "\""
+      in HS.fromList $ F.foldr f [] exprs
+
+    allGraphInputs =
+      let f (_, GSym name) acc = name : acc
+          f _ acc = acc
+      in HS.fromList $ foldr f [] gr
+
+-- | if the same input symbol (like ESym "x") is given at two different places throw an exception
+findConflictingInputs :: Foldable f => f Sym -> [Sym]
+findConflictingInputs syms = HS.toList redundant
+  where
+    redundant = snd $ F.foldl f (HS.empty, HS.empty) syms
+      where
+        f (knownExprs, redundantExprs) s
+          | HS.member s knownExprs = (knownExprs, HS.insert s redundantExprs)
+          | otherwise = (HS.insert s knownExprs, redundantExprs)
+
+-- | Take inputs and outputs and traverse the outputs reifying all expressions
+--   and creating a hashmap of StableNames. Once the hashmap is created,
+--   lookup the provided inputs and return a FunGraph which contains an
+--   expression graph, input/output indices, and other useful functions.
+--   StableNames may be non-deterministic so this function may return graphs
+--   with greater or fewer CSE's eliminated.
+--   If CSE is then performed on the graph, the result is deterministic.
+toFunGraph :: (Functor f, Foldable f, Traversable g) =>
+              f (Expr a) -> g (Expr a) -> IO (FunGraph f g a)
+toFunGraph inputExprs outputExprs = do
+  -- reify the outputs
+  (ReifyGraph rgr, outputIndices) <- reifyGraph outputExprs
+  let userInputSyms = fmap f inputExprs
+        where
+          f (ESym s) = s
+          f _ = error $ "ERROR: toFunGraph given non-ESym input" -- \"" ++ show x ++ "\""
+      fg = FunGraph { fgInputs = userInputSyms
+                    , fgOutputs = outputIndices
+                    , fgReified = reverse rgr
+                    , fgTopSort = topSort
+                    }
+
+      -- make sure all the inputs are symbolic, and find their indices in the Expr graph
+      (gr, lookupVertex, lookupKey) =
+        Graph.graphFromEdges $ map (\(k,gexpr) -> (gexpr, k, F.toList gexpr)) rgr
+      lookupG k = (\(g,_,_) -> g) <$> lookupVertex <$> lookupKey k
+
+      topSort = map lookup' $ reverse $ map ((\(_,k,_) -> k) . lookupVertex) $ Graph.topSort gr
+
+      lookup' k = case lookupG k of
+        Nothing -> error "DVDA internal error"
+        Just g -> (k,g)
+  return $ case (detectMissingInputs inputExprs rgr, findConflictingInputs userInputSyms) of
+    ([],[]) -> fg
+    (xs,[]) -> error $ "toFunGraph found inputs that were not provided by the user: " ++ show xs
+    ( _,xs) -> error $ "toFunGraph found conflicting inputs: " ++ show xs
diff --git a/Dvda/Algorithm/Reify.hs b/Dvda/Algorithm/Reify.hs
new file mode 100644
--- /dev/null
+++ b/Dvda/Algorithm/Reify.hs
@@ -0,0 +1,170 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# Language BangPatterns #-}
+{-# Language ScopedTypeVariables #-}
+
+-- | This file is a modified version from Andy Gill's data-reify package
+--   It is modified to use Data.HashTable.IO, which gives a speed improvement
+--   at the expense of portability. This also gives me a more convenient
+--   sandbox to investigate other performance tweaks, though it is unclear
+--   if I have made anything any faster.
+
+module Dvda.Algorithm.Reify
+       ( ReifyGraph(..)
+       , Node(..)
+       , reifyGraph
+       ) where
+
+import Control.Monad.State.Strict ( StateT(..), runStateT )
+import Data.Hashable ( Hashable(..) )
+import Control.Applicative ( pure )
+import Data.Traversable ( Traversable )
+import qualified Data.Traversable as T
+import System.Mem.StableName ( StableName, makeStableName, hashStableName )
+import Unsafe.Coerce ( unsafeCoerce )
+
+import Dvda.Expr
+
+import qualified Data.HashTable.IO as H
+type HashTable k v = H.CuckooHashTable k v
+
+newtype Node = Node Int deriving (Ord, Eq)
+
+instance Show Node where
+  show (Node k) = '@' : show k
+
+data ReifyGraph e = ReifyGraph [(Node,e Node)]
+
+mapAccumM' :: (Monad m, Functor m, Traversable t) =>
+             (a -> b -> m (c, a)) -> a -> t b -> m (t c, a)
+mapAccumM' f = flip (runStateT . T.traverse (StateT . flip f))
+--{-# INLINE mapAccumM' #-}
+
+mapAccumM :: (Monad m, Functor m, Traversable t) =>
+             (a -> b -> m (a, c)) -> a -> t b -> m (t c, a)
+mapAccumM f' = mapAccumM' f
+  where
+    f acc z = do
+      (x,y) <- f' acc z
+      return (y,x)
+--{-# INLINE mapAccumM #-}
+
+mapDeRef :: (acc -> Expr a -> IO (acc, Node)) -> acc -> Expr a -> IO (acc, GExpr a Node)
+mapDeRef _ acc0 (ESym name) = pure (acc0, GSym name)
+mapDeRef _ acc0 (EConst c)  = pure (acc0, GConst c)
+mapDeRef f acc0 (ENum (Mul x y)) = do
+  (acc1, fx) <- f acc0 x
+  (acc2, fy) <- f acc1 y
+  return (acc2, GNum (Mul fx fy))
+mapDeRef f acc0 (ENum (Add x y)) = do
+  (acc1, fx) <- f acc0 x
+  (acc2, fy) <- f acc1 y
+  return (acc2, GNum (Add fx fy))
+mapDeRef f acc0 (ENum (Sub x y)) = do
+  (acc1, fx) <- f acc0 x
+  (acc2, fy) <- f acc1 y
+  return (acc2, GNum (Sub fx fy))
+mapDeRef f acc0 (ENum (Negate x)) = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GNum (Negate fx))
+mapDeRef f acc0 (ENum (Abs x)) = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GNum (Abs fx))
+mapDeRef f acc0 (ENum (Signum x)) = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GNum (Signum fx))
+mapDeRef _ acc0 (ENum (FromInteger k)) = pure (acc0, GNum (FromInteger k))
+mapDeRef f acc0 (EFractional (Div x y)) = do
+  (acc1, fx) <- f acc0 x
+  (acc2, fy) <- f acc1 y
+  return (acc2, GFractional (Div fx fy))
+mapDeRef _ acc0 (EFractional (FromRational x)) = pure (acc0, GFractional (FromRational x))
+mapDeRef f acc0 (EFloating (Pow x y))     = do
+  (acc1, fx) <- f acc0 x
+  (acc2, fy) <- f acc1 y
+  return (acc2, GFloating (Pow fx fy))
+mapDeRef f acc0 (EFloating (LogBase x y)) = do
+  (acc1, fx) <- f acc0 x
+  (acc2, fy) <- f acc1 y
+  return (acc2, GFloating (LogBase fx fy))
+mapDeRef f acc0 (EFloating (Exp   x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Exp fx))
+mapDeRef f acc0 (EFloating (Log   x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Log   fx))
+mapDeRef f acc0 (EFloating (Sin   x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Sin   fx))
+mapDeRef f acc0 (EFloating (Cos   x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Cos   fx))
+mapDeRef f acc0 (EFloating (Tan   x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Tan   fx))
+mapDeRef f acc0 (EFloating (ASin  x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (ASin  fx))
+mapDeRef f acc0 (EFloating (ATan  x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (ATan  fx))
+mapDeRef f acc0 (EFloating (ACos  x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (ACos  fx))
+mapDeRef f acc0 (EFloating (Sinh  x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Sinh  fx))
+mapDeRef f acc0 (EFloating (Cosh  x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Cosh  fx))
+mapDeRef f acc0 (EFloating (Tanh  x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (Tanh  fx))
+mapDeRef f acc0 (EFloating (ASinh x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (ASinh fx))
+mapDeRef f acc0 (EFloating (ATanh x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (ATanh fx))
+mapDeRef f acc0 (EFloating (ACosh x))     = do
+  (acc1, fx) <- f acc0 x
+  return (acc1, GFloating (ACosh fx))
+-- {-# INLINE mapDeRef #-}
+
+
+reifyGraph :: forall t a . Traversable t => t (Expr a) -> IO (ReifyGraph (GExpr a), t Node)
+reifyGraph m = do
+  ht <- H.new :: IO (HashTable DynStableName Node)
+  let findNodes :: ([(Node, GExpr a Node)],Node) -> Expr a ->
+                    IO (([(Node, GExpr a Node)],Node), Node)
+      findNodes !(!tab0, nextUnique@(Node nextUnique')) expr = do
+        stableName <- makeDynStableName expr
+        lu <- H.lookup ht stableName
+        case lu of
+          Just var -> return ((tab0,nextUnique), var)
+          Nothing -> do
+            let var = nextUnique
+            H.insert ht stableName var
+            ((tab1,nextNextUnique), res) <- mapDeRef findNodes (tab0, Node (nextUnique' + 1)) expr
+            let tab2 :: [(Node,GExpr a Node)]
+                tab2 = (var,res) : tab1
+            return ((tab2,nextNextUnique), var)
+      -- {-# INLINE findNodes #-}
+
+  (root, (pairs,_)) <- mapAccumM findNodes ([], Node 0) m
+  return (ReifyGraph pairs, root)
+
+
+-- Stable names that not use phantom types.
+-- As suggested by Ganesh Sittampalam.
+newtype DynStableName = DynStableName (StableName ()) deriving Eq
+
+instance Hashable DynStableName where
+  hashWithSalt salt = (salt `hashWithSalt`) . hashDynStableName
+hashDynStableName :: DynStableName -> Int
+hashDynStableName (DynStableName sn) = hashStableName sn
+
+makeDynStableName :: a -> IO DynStableName
+makeDynStableName !a = do
+  st <- makeStableName a
+  return $ DynStableName (unsafeCoerce st)
+--{-# INLINE makeDynStableName #-}
diff --git a/Dvda/CGen.hs b/Dvda/CGen.hs
deleted file mode 100644
--- a/Dvda/CGen.hs
+++ /dev/null
@@ -1,295 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language TemplateHaskell #-}
-{-# Language TypeFamilies #-}
-{-# Language FlexibleContexts #-}
-
-module Dvda.CGen ( showC
-                 , showMex
-                 , MatrixStorageOrder(..)
-                 ) where
-
-
-import Data.Hashable ( Hashable )
-import Data.List ( intercalate )
-import FileLocation ( err )
-import Text.Printf ( printf )
-
-import Dvda.Expr ( GExpr(..), Floatings(..), Nums(..), Fractionals(..) )
-import Dvda.FunGraph ( FunGraph, MVS(..), topSort, fgInputs, fgOutputs, fgLookupGExpr )
-import Dvda.HashMap ( HashMap )
-import qualified Dvda.HashMap as HM
-
-data MatrixStorageOrder = RowMajor | ColMajor
-
--- | take a list of pair of inputs to indices which reference them
---  create a hashmap from GSyms to strings which hold the declaration
-makeInputMap :: (Eq a, Hashable a, Show a)
-                => MatrixStorageOrder -> [MVS (GExpr a Int)] -> HashMap (GExpr a Int) String
-makeInputMap matStorageOrder ins = HM.fromList $ concat $ zipWith writeInput [(0::Int)..] ins
-  where
-    writeInput inputK (Sca g) = [(g, printf "*input%d; /* %s */" inputK (show g))]
-    writeInput inputK (Vec gs) = zipWith f [(0::Int)..] gs
-      where
-        f inIdx g = (g, printf "input%d[%d]; /* %s */" inputK inIdx (show g))
-    writeInput inputK (Mat gs)
-      | any ((ncols /=) . length) gs =
-          error $ "writeInputs [[GraphRef]] matrix got inconsistent column dimensions: "++ show (map length gs)
-      | otherwise = zipWith f [(r,c) | r <- [0..(nrows-1)], c <- [0..(ncols-1)]] (concat gs)
-      where
-        nrows = length gs
-        ncols = if nrows == 0 then 0 else length (head gs)
-        f (rowIdx,colIdx) g = (g,printf "input%d[%d][%d]; /* %s */" inputK fstIdx sndIdx (show g))
-          where
-            (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (rowIdx,colIdx)
-                                                      ColMajor -> (colIdx,rowIdx)
-
-writeInputPrototypes :: MatrixStorageOrder -> [MVS a] -> [String]
-writeInputPrototypes matStorageOrder ins = concat $ zipWith inputPrototype [(0::Int)..] ins
-  where
-    inputPrototype inputK (Sca _) = ["const double * input" ++ show inputK]
-    inputPrototype inputK (Vec gs) = ["const double input" ++ show inputK ++ "[" ++ show (length gs) ++ "]"]
-    inputPrototype inputK (Mat gs)
-      | any ((ncols /=) . length) gs =
-          error $ "writeInputs [[GraphRef]] matrix got inconsistent column dimensions: "++ show (map length gs)
-      | otherwise = ["const double input" ++ show inputK ++ "[" ++ show fstIdx ++ "][" ++ show sndIdx ++ "]"]
-      where
-        nrows = length gs
-        ncols = if nrows == 0 then 0 else length (head gs)
-        (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (nrows,ncols)
-                                                  ColMajor -> (ncols,nrows)
-
-writeOutputs :: MatrixStorageOrder -> [MVS Int] -> ([String], [String])
-writeOutputs matStorageOrder ins = (concatMap fst dcs, concatMap snd dcs)
-  where
-    dcs :: [([String],[String])]
-    dcs = zipWith writeOutput ins [0..]
-
-    writeOutput :: MVS Int -> Int -> ([String], [String])
-    writeOutput (Sca gref) outputK = (decls, prototype)
-      where
-        decls = [printf "/* output %d */" outputK, printf "(*output%d) = %s;" outputK (nameNode gref)]
-        prototype = ["double * const output" ++ show outputK]
-    writeOutput (Vec grefs) outputK = (decls, prototype)
-      where
-        prototype = ["double output" ++ show outputK ++ "[" ++ show (length grefs) ++ "]"]
-        decls = (printf "/* output %d */" outputK):
-                zipWith f [(0::Int)..] grefs
-          where
-            f outIdx gref = printf "output%d[%d] = %s;" outputK outIdx (nameNode gref)
-    writeOutput (Mat grefs) outputK
-      | any ((ncols /=) . length) grefs =
-          error $ "writeOutputs [[GraphRef]] matrix got inconsistent column dimensions: "++ show (map length grefs)
-      | otherwise = (decls, prototype)
-      where
-        nrows = length grefs
-        ncols = if nrows == 0 then 0 else length (head grefs)
-        prototype = ["double output" ++ show outputK ++ "[" ++ show fstIdx ++ "][" ++ show sndIdx ++ "]"]
-          where
-            (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (nrows,ncols)
-                                                      ColMajor -> (ncols,nrows)
-        decls = (printf "/* output %d */" outputK):
-                zipWith f [(r,c) | r <- [0..(nrows-1)], c <- [0..(ncols-1)]] (concat grefs)
-          where
-            f (rowIdx,colIdx) gref = printf "output%d[%d][%d] = %s;" outputK fstIdx sndIdx (nameNode gref)
-              where
-                (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (rowIdx,colIdx)
-                                                          ColMajor -> (colIdx,rowIdx)
-
-
-createMxOutputs :: [MVS Int] -> [String]
-createMxOutputs xs = concat $ zipWith createMxOutput xs [0..]
-  where
-    createMxOutput :: MVS Int -> Int -> [String]
-    createMxOutput (Sca _) outputK =
-      [ "    if ( " ++ show outputK ++ " < nlhs ) {"
-      , "        plhs[" ++ show outputK ++ "] = mxCreateDoubleScalar( 0 );"
-      , "        outputs[" ++ show outputK ++ "] = mxGetPr( plhs[" ++ show outputK ++ "] );"
-      , "    } else"
-      , "        outputs[" ++ show outputK ++ "] = (double*)malloc( sizeof(double) );"
-      ]
-    createMxOutput (Vec grefs) outputK =
-      [ "    if ( " ++ show outputK ++ " < nlhs ) {"
-      , "        plhs[" ++ show outputK ++ "] = mxCreateDoubleMatrix( " ++ show (length grefs) ++ ", 1, mxREAL );"
-      , "        outputs[" ++ show outputK ++ "] = mxGetPr( plhs[" ++ show outputK ++ "] );"
-      , "    } else"
-      , "        outputs[" ++ show outputK ++ "] = (double*)malloc( " ++ show (length grefs) ++ "*sizeof(double) );"
-      ]
-    createMxOutput (Mat grefs) outputK =
-      [ "    if ( " ++ show outputK ++ " < nlhs ) {"
-      , "        plhs[" ++ show outputK ++ "] = mxCreateDoubleMatrix( " ++ show nrows++ ", " ++ show ncols ++ ", mxREAL );"
-      , "        outputs[" ++ show outputK ++ "] = mxGetPr( plhs[" ++ show outputK ++ "] );"
-      , "    } else"
-      , "        outputs[" ++ show outputK ++ "] = (double*)malloc( " ++ show (nrows*ncols) ++ "*sizeof(double) );"
-      ]
-      where
-        nrows = length grefs
-        ncols = if nrows == 0 then 0 else length (head grefs)
-
-
-checkMxInputDims :: MVS a -> String -> Int -> [String]
-checkMxInputDims (Sca _) functionName inputK =
-  [ "    if ( 1 != mxGetM( prhs[" ++ show inputK ++ "] ) || 1 != mxGetN( prhs[" ++ show inputK ++ "] ) ) {"
-  , "        char errMsg[200];"
-  , "        sprintf(errMsg,"
-  , "                \"mex function '" ++ functionName ++ "' got incorrect dimensions for input " ++ show (1+inputK) ++ "\\n\""
-  , "                \"expected dimensions: (1, 1) but got (%zu, %zu)\","
-  , "                mxGetM( prhs[" ++ show inputK ++ "] ),"
-  , "                mxGetN( prhs[" ++ show inputK ++ "] ) );"
-  , "        mexErrMsgTxt(errMsg);"
-  , "    }"
-  ]
-checkMxInputDims (Vec grefs) functionName inputK =
-  [ "    if ( !( " ++ show nrows ++ " == mxGetM( prhs[" ++ show inputK ++ "] ) && 1 == mxGetN( prhs[" ++ show inputK ++ "] ) ) && !( " ++ show nrows ++ " == mxGetN( prhs[" ++ show inputK ++ "] ) && 1 == mxGetM( prhs[" ++ show inputK ++ "] ) ) ) {"
-  , "        char errMsg[200];"
-  , "        sprintf(errMsg,"
-  , "                \"mex function '" ++ functionName ++ "' got incorrect dimensions for input " ++ show (1+inputK) ++ "\\n\""
-  , "                \"expected dimensions: (" ++ show nrows ++ ", 1) or (1, " ++ show nrows ++ ") but got (%zu, %zu)\","
-  , "                mxGetM( prhs[" ++ show inputK ++ "] ),"
-  , "                mxGetN( prhs[" ++ show inputK ++ "] ) );"
-  , "        mexErrMsgTxt(errMsg);"
-  , "    }"
-  ]
-  where
-    nrows = length grefs
-checkMxInputDims (Mat grefs) functionName inputK =
-  [ "    if ( " ++ show nrows ++ " != mxGetM( prhs[" ++ show inputK ++ "] ) || " ++ show ncols ++ " != mxGetN( prhs[" ++ show inputK ++ "] ) ) {"
-  , "        char errMsg[200];"
-  , "        sprintf(errMsg,"
-  , "                \"mex function '" ++ functionName ++ "' got incorrect dimensions for input " ++ show (1+inputK) ++ "\\n\""
-  , "                \"expected dimensions: (" ++ show nrows ++ ", " ++ show ncols ++ ") but got (%zu, %zu)\","
-  , "                mxGetM( prhs[" ++ show inputK ++ "] ),"
-  , "                mxGetN( prhs[" ++ show inputK ++ "] ) );"
-  , "        mexErrMsgTxt(errMsg);"
-  , "    }"
-  ]
-  where
-    nrows = length grefs
-    ncols = if nrows == 0 then 0 else length (head grefs)
-
-
--- | Turns a FunGraph into a string containing C code
-showC :: (Eq a, Show a, Hashable a) => MatrixStorageOrder -> String -> FunGraph a -> String
-showC matStorageOrder functionName fg = txt
-  where
-    inPrototypes = writeInputPrototypes matStorageOrder (fgInputs fg)
-    (outDecls, outPrototypes) = writeOutputs matStorageOrder (fgOutputs fg)
-    inputMap = makeInputMap matStorageOrder (fgInputs fg)
-    mainDecls = let f k = case fgLookupGExpr fg k of
-                      Just v -> cAssignment inputMap k v
-                      Nothing -> error $ "couldn't find node " ++ show k ++ " in fungraph :("
-                in map f $ reverse $ topSort fg
-  
-    body = unlines $ map ("    "++) $
-           mainDecls ++ [""] ++
-           outDecls
-  
-    txt = "#include <math.h>\n\n" ++
-          "void " ++ functionName ++ " ( " ++ (intercalate ", " (inPrototypes++outPrototypes)) ++ " )\n{\n" ++
-          body ++ "}\n"
-
-nameNode :: Int -> String
-nameNode k = "v_" ++ show k
-
-cAssignment :: (Eq a, Hashable a, Show a) => HashMap (GExpr a Int) String -> Int -> GExpr a Int -> String
-cAssignment inputMap k g@(GSym _) = case HM.lookup g inputMap of
-  Nothing -> error $ "cAssignment: couldn't find " ++ show g ++ " in the input map"
-  Just str -> "const double " ++ nameNode k ++ " = " ++ str
-cAssignment inputMap k gexpr = "const double " ++ nameNode k ++ " = " ++ toCOp gexpr ++ ";"
-  where
-    bin :: Int -> Int -> String -> String
-    bin x y op = nameNode x ++ " " ++ op ++ " " ++ nameNode y
-    
-    un :: Int -> String -> String
-    un x op = op ++ "( " ++ nameNode x ++ " )"
-
-    asTypeOfG :: a -> GExpr a b -> a
-    asTypeOfG x _ = x
-    
-    toCOp (GSym _)                       = $(err "This should be impossible")
-    toCOp (GConst c)                     = show c
-    toCOp (GNum (Mul x y))               = bin x y "*"
-    toCOp (GNum (Add x y))               = bin x y "+"
-    toCOp (GNum (Sub x y))               = bin x y "-"
-    toCOp (GNum (Negate x))              = un x "-"
-    toCOp (GNum (Abs x))                 = un x "abs"
-    toCOp (GNum (Signum x))              = un x "sign"
-    toCOp (GNum (FromInteger x))         = show x
-    toCOp (GFractional (Div x y))        = bin x y "/"
-    toCOp (GFractional (FromRational x)) = show (fromRational x `asTypeOfG` gexpr)
-    toCOp (GFloating (Pow x y))          = "pow( " ++ nameNode x ++ ", " ++ nameNode y ++ " )"
-    toCOp (GFloating (LogBase x y))      = "log( " ++ nameNode y ++ ") / log( " ++ nameNode x ++ " )"
-    toCOp (GFloating (Exp x))            = un x "exp"
-    toCOp (GFloating (Log x))            = un x "log"
-    toCOp (GFloating (Sin x))            = un x "sin"
-    toCOp (GFloating (Cos x))            = un x "cos"
-    toCOp (GFloating (ASin x))           = un x "asin"
-    toCOp (GFloating (ATan x))           = un x "atan"
-    toCOp (GFloating (ACos x))           = un x "acos"
-    toCOp (GFloating (Sinh x))           = un x "sinh"
-    toCOp (GFloating (Cosh x))           = un x "cosh"
-    toCOp (GFloating (Tanh x))           = un x "tanh"
-    toCOp (GFloating (ASinh _))          = error "C generation doesn't support ASinh"
-    toCOp (GFloating (ATanh _))          = error "C generation doesn't support ATanh"
-    toCOp (GFloating (ACosh _))          = error "C generation doesn't support ACosh"
-
-
-showMex :: (Eq a, Show a, Hashable a) => String -> FunGraph a -> String
-showMex functionName fg = cText ++ "\n\n\n" ++ mexFun functionName (fgInputs fg) (fgOutputs fg)
-  where
-    cText = showC ColMajor functionName fg -- matlab is column major >_<
-
-mexFun :: String -> [MVS a] -> [MVS Int] -> String
-mexFun functionName ins outs =
-  unlines $
-  [ "#include \"mex.h\""
-  , []
-  , "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])"
-  , "{"
-  , "    /* check number of inputs  */"
-  , "    if ( " ++ show nrhs ++ " != nrhs ) {"
-  , "        char errMsg[200];"
-  , "        sprintf(errMsg,"
-  , "                \"mex function '" ++ functionName ++ "' given incorrect number of inputs\\n\""
-  , "                \"expected: " ++ show nrhs ++ " but got %d\","
-  , "                nrhs);"
-  , "        mexErrMsgTxt(errMsg);"
-  , "    }"
-  , []
-  , "    /* check the dimensions of the input arrays */"
-  ] ++ concat (zipWith (\x -> checkMxInputDims x functionName) ins [0..]) ++
-  [ []
-  , "    /* check number of outputs  */"
-  , "    if ( " ++ show nlhs ++ " < nlhs ) {"
-  , "        char errMsg[200];"
-  , "        sprintf(errMsg,"
-  , "                \"mex function '" ++ functionName ++ "' saw too many outputs\\n\""
-  , "                \"expected <= " ++ show nlhs ++ " but got %d\","
-  , "                nlhs);"
-  , "        mexErrMsgTxt(errMsg);"
-  , "    }"
-  , []
-  , "    /* create the output arrays, if no output is provided by user create a dummy output */"
-  , "    double * outputs[" ++ show nlhs ++ "];"
-  ] ++ createMxOutputs outs ++ -- e.g.: plhs[0] = mxCreateDoubleMatrix(1,ncols,mxREAL);
-  [ []
-  , "    /* call the c function */"
-  , "    " ++ functionName ++ "( " ++ intercalate ", " (inputPtrs ++ outputPtrs) ++ " );"
-  , []
-  , "    /* free the unused dummy outputs */"
-  , "    int k;"
-  , "    for ( k = " ++ show (nlhs - 1) ++ "; nlhs <= k; k-- )"
-  , "        free( outputs[k] );"
-  , "}"
-  ]
-  where
-    nlhs = length outs
-    nrhs = length ins
-    inputPtrs  = zipWith (\i k -> cast i "const " ++ "mxGetPr(prhs[" ++ show k ++ "])") ins  [(0::Int)..]
-    outputPtrs = zipWith (\o k -> cast o    ""    ++ "(outputs[" ++ show k ++ "])")     outs [(0::Int)..]
-
-    cast :: MVS a -> String -> String
-    cast (Sca _) _ = ""
-    cast (Vec _) _ = ""
-    cast (Mat xs) cnst = "(" ++ cnst ++ "double (*)[" ++ show nrows ++ "])" -- column major order
-      where
-        nrows = length xs
diff --git a/Dvda/CSE.hs b/Dvda/CSE.hs
deleted file mode 100644
--- a/Dvda/CSE.hs
+++ /dev/null
@@ -1,153 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.CSE ( cse
-                ) where
-
-import Control.Monad.ST ( ST, runST )
-import Data.Foldable ( toList )
-import Data.Hashable ( Hashable )
-import Data.IntMap ( IntMap )
-import qualified Data.IntMap as IM
-import Data.Tuple ( swap )
-
-import Dvda.Expr ( GExpr(..), Floatings(..), Fractionals(..), Nums(..) )
-import Dvda.FunGraph
-
-import qualified Data.HashTable.Class as HT
-import qualified Data.HashTable.ST.Cuckoo as C
-type HashTable s v k = C.HashTable s v k
-
-cse :: (Eq a, Hashable a) => FunGraph a -> FunGraph a
-cse fg = nodelistToFunGraph (map swap htList) (fgInputs fg) outputIndices
-  where
-    (htList, im) = cse' (fgLookupGExpr fg) (fgOutputs fg)
-    -- since the fgInputs are all symbolic (GSym _) there is no need for mapping old inputs to new inputs
-    outputIndices = let
-      oldIndexToNewIndex k = case IM.lookup k im of
-        Just k' -> k'
-        Nothing -> error $
-                   "CSE error, in mapping old output indices to new, found an old one which was missing from" ++
-                   "the old --> new Int mapping"
-      in map (fmap oldIndexToNewIndex) (fgOutputs fg)
-
-cse' ::
-  (Eq a, Hashable a)
-  => (Int -> Maybe (GExpr a Int))
-  -> [MVS Int]
-  -> ([(GExpr a Int, Int)], IntMap Int)
-cse' lookupFun outputIndices = runST $ do
-  ht <- HT.new
-  let -- folding function
-      f (im,n) [] = return (im,n)
-      f (im0,n0) (k:ks) = do
-        (_,im,n) <- insertOldNode k lookupFun ht im0 n0
-        f (im,n) ks
-  -- outputs
-  (oldToNewIdx,_) <- f (IM.empty,0) (concatMap toList outputIndices)
-  htList <- HT.toList ht
-  return (htList, oldToNewIdx)
-
-  
----- | take in an Int that represents a node in the original graph
----- see if that int has been inserted in the new graph
-insertOldNode ::
-  (Eq a, Hashable a)
-  => Int -- ^ Int to be inserted
-  -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference
-  -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references
-  -> IntMap Int -- ^ intmap of old int reference to new int references
-  -> Int -- ^ next free index
-  -> ST s (Int, IntMap Int, Int)
-insertOldNode kOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0 =
-  case IM.lookup kOld oldNodeToNewNode0 of
-    -- if the int has already been inserted in the new graph, return it
-    Just k -> return (k, oldNodeToNewNode0, nextFreeInt0)
-    -- if the int has not yet been inserted, then insert it
-    -- get the old GExpr to which this node corresponds
-    Nothing ->  case lookupOldGExpr kOld of
-      Nothing -> error $ "in CSE, insertOldNode got an old key \"" ++ show kOld ++
-                 "\" with was not found in the old graph"
-      -- insert this old GExpr
-      Just oldGExpr -> do
-        (k, oldNodeToNewNode1, nextFreeInt1) <- insertOldGExpr oldGExpr lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0
-        return (k, IM.insert kOld k oldNodeToNewNode1, nextFreeInt1)
-
-insertOldGExpr ::
-  (Eq a, Hashable a)
-  => GExpr a Int -- ^ GExpr to be inserted
-  -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference
-  -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references
-  -> IntMap Int -- ^ intmap of old int reference to new int references
-  -> Int -- ^ next free index
-  -> ST s (Int, IntMap Int, Int)
-
-insertOldGExpr g@(GSym _)                       = \_ ->  cseInsert g
-insertOldGExpr g@(GConst _)                     = \_ ->  cseInsert g
-insertOldGExpr g@(GNum (FromInteger _))         = \_ ->  cseInsert g
-insertOldGExpr g@(GFractional (FromRational _)) = \_ ->  cseInsert g
-
-insertOldGExpr (GNum (Mul x y))          = insertOldGExprBinary GNum Mul x y
-insertOldGExpr (GNum (Add x y))          = insertOldGExprBinary GNum Add x y
-insertOldGExpr (GNum (Sub x y))          = insertOldGExprBinary GNum Sub x y
-insertOldGExpr (GFractional (Div x y))   = insertOldGExprBinary GFractional Div x y
-insertOldGExpr (GFloating (Pow x y))     = insertOldGExprBinary GFloating Pow x y
-insertOldGExpr (GFloating (LogBase x y)) = insertOldGExprBinary GFloating LogBase x y
-                                         
-insertOldGExpr (GNum (Negate x))         = insertOldGExprUnary  GNum Negate x
-insertOldGExpr (GNum (Abs x))            = insertOldGExprUnary  GNum Abs x
-insertOldGExpr (GNum (Signum x))         = insertOldGExprUnary  GNum Signum x
-insertOldGExpr (GFloating (Exp x))       = insertOldGExprUnary  GFloating Exp x
-insertOldGExpr (GFloating (Log x))       = insertOldGExprUnary  GFloating Log x
-insertOldGExpr (GFloating (Sin x))       = insertOldGExprUnary  GFloating Sin x
-insertOldGExpr (GFloating (Cos x))       = insertOldGExprUnary  GFloating Cos x
-insertOldGExpr (GFloating (ASin x))      = insertOldGExprUnary  GFloating ASin x
-insertOldGExpr (GFloating (ATan x))      = insertOldGExprUnary  GFloating ATan x
-insertOldGExpr (GFloating (ACos x))      = insertOldGExprUnary  GFloating ACos x
-insertOldGExpr (GFloating (Sinh x))      = insertOldGExprUnary  GFloating Sinh x
-insertOldGExpr (GFloating (Cosh x))      = insertOldGExprUnary  GFloating Cosh x
-insertOldGExpr (GFloating (Tanh x))      = insertOldGExprUnary  GFloating Tanh x
-insertOldGExpr (GFloating (ASinh x))     = insertOldGExprUnary  GFloating ASinh x
-insertOldGExpr (GFloating (ATanh x))     = insertOldGExprUnary  GFloating ATanh x
-insertOldGExpr (GFloating (ACosh x))     = insertOldGExprUnary  GFloating ACosh x
-
-insertOldGExprBinary ::
-  (Eq a, Hashable a)
-  => (f -> GExpr a Int)
-  -> (Int -> Int -> f)
-  -> Int -> Int
-  -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference
-  -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references
-  -> IntMap Int -- ^ intmap of old int reference to new int references
-  -> Int -- ^ next free index
-  -> ST s (Int, IntMap Int, Int)
-insertOldGExprBinary gnum mul kxOld kyOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0 = do
-  (kx, oldNodeToNewNode1,nextFreeInt1) <- insertOldNode kxOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0
-  (ky, oldNodeToNewNode2,nextFreeInt2) <- insertOldNode kyOld lookupOldGExpr ht oldNodeToNewNode1 nextFreeInt1
-  let newGExpr = gnum (mul kx ky)
-  cseInsert newGExpr ht oldNodeToNewNode2 nextFreeInt2
-
-insertOldGExprUnary ::
-  (Eq a, Hashable a)
-  => (f -> GExpr a Int)
-  -> (Int -> f)
-  -> Int
-  -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference
-  -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references
-  -> IntMap Int -- ^ intmap of old int reference to new int references
-  -> Int -- ^ next free index
-  -> ST s (Int, IntMap Int, Int)
-insertOldGExprUnary gnum neg kxOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0 = do
-  (kx, oldNodeToNewNode1,nextFreeInt1) <- insertOldNode kxOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0
-  let newGExpr = gnum (neg kx)
-  cseInsert newGExpr ht oldNodeToNewNode1 nextFreeInt1
-
-cseInsert :: (Eq a, Hashable a) => GExpr a Int -> HashTable s (GExpr a Int) Int -> IntMap Int -> Int
-             -> ST s (Int, IntMap Int, Int)
-cseInsert gexpr ht oldNodeToNewNode0 nextFreeInt0 = do
-  lu <- HT.lookup ht gexpr
-  case lu of
-    Just k -> return (k, oldNodeToNewNode0, nextFreeInt0)
-    Nothing -> do
-      HT.insert ht gexpr nextFreeInt0
-      return (nextFreeInt0, oldNodeToNewNode0, nextFreeInt0+1)
-        
diff --git a/Dvda/Codegen/Gcc.hs b/Dvda/Codegen/Gcc.hs
deleted file mode 100644
--- a/Dvda/Codegen/Gcc.hs
+++ /dev/null
@@ -1,32 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.Codegen.Gcc ( compileWithGcc
-                        ) where
-
-import System.Process(runCommand, waitForProcess)
-import System.Exit(ExitCode(ExitSuccess))
-import Control.Monad(when)
-
--- | whether to print the gcc call when generating code
-spewGccCall :: Bool
-spewGccCall = True
-
--- | take in source file and object, return string suitible for calling to compile
-gccString :: FilePath -> FilePath -> String
-gccString src obj = "gcc -O2 -std=gnu99 -fPIC -shared -Wall -Wextra -Werror " ++ src ++ " -o " ++ obj
-
--- | take in name of source and future object, compile object
-compileWithGcc :: FilePath -> FilePath -> IO ()
-compileWithGcc srcname objname = do
-  -- compile new object
-  let compileString = gccString srcname objname
-
-  -- print compilation string
-  when spewGccCall $ putStrLn compileString
-  
-  -- run compilation string
-  p <- runCommand compileString
-  
-  -- check for errors
-  exitCode <- waitForProcess p
-  when (exitCode /= ExitSuccess) $ error $ "failed compiling " ++ srcname
diff --git a/Dvda/Codegen/WriteFile.hs b/Dvda/Codegen/WriteFile.hs
deleted file mode 100644
--- a/Dvda/Codegen/WriteFile.hs
+++ /dev/null
@@ -1,32 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.Codegen.WriteFile ( writeSourceFile
-                              ) where
-
-import System.Directory
-import Control.Monad ( when )
-
----- | return directory to use for temp files
----- | create this directory and print message if it doesn't exist
---dvdaDir :: IO FilePath
---dvdaDir = do
---  dir <- getAppUserDataDirectory "dvda"
-
-writeSourceFile :: String -> FilePath -> FilePath -> IO FilePath
-writeSourceFile source functionDir sourceName = do
-  -- make function directory if it doesn't exist
-  createDirectoryIfMissing False functionDir
-  
-  -- filenames
-  let sourcePath  = functionDir ++ "/" ++ sourceName
-      
-  -- if the source already exists, make sure it matches the old source
-  srcExists <- doesFileExist sourcePath
-  when srcExists $ do
-    putStrLn $ "file \"" ++ sourcePath ++ "\" already exists, overwriting"
-  
-  -- write  source
-  putStrLn $ "writing " ++ sourcePath
-  writeFile sourcePath source
-
-  return sourcePath
diff --git a/Dvda/Examples.hs b/Dvda/Examples.hs
deleted file mode 100644
--- a/Dvda/Examples.hs
+++ /dev/null
@@ -1,54 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.Examples ( doCse
-                     , showFg
-                     , cgen
-                     , mexgen
-                     ) where
-
-import Dvda.Expr
-import Dvda.FunGraph
-import Dvda.CGen
-import Dvda.Vis ( previewGraph )
-import Dvda.CSE ( cse )
-import Dvda.AD ( rad )
-
--- a random function to use in different examples
-someFunGraph :: IO (FunGraph Double)
-someFunGraph = toFunGraph inputs outputs
-  where
-    x = sym "x" :: Expr Double
-    y = sym "y"
-    z = sym "z"
-    w = sym "w"
-    w1 = sym "w1"
-    w2 = sym "w2"
-    w3 = sym "w3"
-    f0 = x*y + z + w1 + w2
-    f2 = f0 * w2/w3
-    
-    f1 = [f0/2, f0*y, w, 0.0, 0]
-    boo = x
-
-    inputs = boo :* [y]:*[[z]] :* [w3,w1,w2,w]
-    outputs = f0:*f1:*f2:*[[f0*f0]]:*(rad f2 [x,y,z,w,w1,w2,w3])
-
--- | do cse on a fungraph and count nodes
-doCse :: IO ()
-doCse = do
-  fg' <- someFunGraph
-  putStrLn $ "fungraph has " ++ show (countNodes fg') ++ " nodes"
-  let fg = cse fg'
-  putStrLn $ "fungraph has " ++ show (countNodes fg) ++ " nodes after cse"
-
--- | show a fungraph
-showFg :: IO ()
-showFg = someFunGraph >>= previewGraph
-
--- | c code generation
-cgen :: IO ()
-cgen = fmap (showC RowMajor "foo") someFunGraph >>= putStrLn
-
--- | mex function generation
-mexgen :: IO ()
-mexgen = fmap (showMex "foo") someFunGraph >>= putStrLn
diff --git a/Dvda/Expr.hs b/Dvda/Expr.hs
--- a/Dvda/Expr.hs
+++ b/Dvda/Expr.hs
@@ -1,11 +1,12 @@
 {-# OPTIONS_GHC -Wall #-}
+{-# Language StandaloneDeriving #-}
 {-# Language GADTs #-}
-{-# Language TemplateHaskell #-}
 {-# Language TypeFamilies #-}
-{-# Language StandaloneDeriving #-}
+{-# Language DeriveGeneric #-}
 {-# Language DeriveDataTypeable #-}
-{-# Language FlexibleInstances #-}
-{-# Language FlexibleContexts #-}
+{-# Language DeriveFunctor #-}
+{-# Language DeriveFoldable #-}
+{-# Language DeriveTraversable #-}
 
 module Dvda.Expr ( Expr(..)
                  , GExpr(..)
@@ -13,12 +14,10 @@
                  , Fractionals(..)
                  , Floatings(..)
                  , Sym(..)
-                 , isVal
                  , sym
                  , symDependent
                  , symDependentN
                  , const'
-                 , getParents
                  , extractLinearPart
                  , getConst
                  , substitute
@@ -27,15 +26,15 @@
                  , fromNeg
                  ) where
 
-import Control.Applicative ( (<$>), (<*>), pure )
-import Data.Data ( Data, Typeable, Typeable1, Typeable2 )
-import Data.Hashable ( Hashable, hash, combine )
+import Control.Applicative ( (<$>), pure )
+import Data.Hashable ( Hashable(..), hash )
 import Data.Ratio ( (%) )
-
---import Test.QuickCheck -- ( Arbitrary(..) )
+import GHC.Generics ( Generic )
+import Data.Monoid ( mempty )
+import qualified Data.Foldable as F
+import qualified Data.Traversable as T
 
 import qualified Dvda.HashMap as HM
-import Dvda.Reify ( MuRef(..) )
 
 commutativeMul :: Bool
 commutativeMul = True
@@ -43,13 +42,16 @@
 commutativeAdd :: Bool
 commutativeAdd = True
 
-data Sym = Sym String                  -- doesn't depend on independent variable, or is an independent variable
-         | SymDependent String Int Sym -- depends on independent variable, Int specifies the nth derivative
-           deriving (Eq, Ord)
+data Sym =
+    Sym String -- doesn't depend on independent variable, or is an independent variable
+  | SymDependent String Int Sym -- depends on independent variable, Int specifies the nth derivative
+  deriving (Eq, Ord, Generic)
 
 instance Show Sym where
-  show (Sym name) = name
-  show (SymDependent name k s) = name ++ replicate k '\'' ++ "(" ++ show s ++ ")"
+  showsPrec d (Sym name) = showParen (d >= 9) $ showString name
+  showsPrec d (SymDependent name k s) =
+    showParen (d >= 9) $
+    showString $ name ++ replicate k '\'' ++ "(" ++ show s ++ ")"
 
 data Expr a where
   ESym :: Sym -> Expr a
@@ -64,10 +66,12 @@
             | Negate a
             | Abs a
             | Signum a
-            | FromInteger Integer deriving Ord
+            | FromInteger Integer
+            deriving (Ord, Generic, Functor, F.Foldable, T.Traversable)
 
 data Fractionals a = Div a a
-                   | FromRational Rational deriving (Eq, Ord)
+                   | FromRational Rational
+                   deriving (Eq, Ord, Generic, Functor, F.Foldable, T.Traversable)
 
 data Floatings a = Pow a a
                  | LogBase a a
@@ -75,6 +79,7 @@
                  | Log a
                  | Sin a
                  | Cos a
+                 | Tan a
                  | ASin a
                  | ATan a
                  | ACos a
@@ -83,31 +88,18 @@
                  | Tanh a
                  | ASinh a
                  | ATanh a
-                 | ACosh a deriving (Eq, Ord)
-
-deriving instance Data Sym
-deriving instance Data a => Data (Nums a)
-deriving instance Data a => Data (Fractionals a)
-deriving instance Data a => Data (Floatings a)
-deriving instance (Data a, Floating a) => Data (Expr a)
-deriving instance (Data a, Data b, Floating a) => Data (GExpr a b)
-
-deriving instance Typeable Sym
-deriving instance Typeable1 Nums
-deriving instance Typeable1 Fractionals
-deriving instance Typeable1 Floatings
-deriving instance Typeable1 Expr
-deriving instance Typeable2 GExpr
+                 | ACosh a
+                 deriving (Eq, Ord, Generic, Functor, F.Foldable, T.Traversable)
 
 ----------------------- Show instances -------------------------
 showsInfixBinary :: (Show a, Show b) => Int -> Int -> String -> a -> b -> ShowS
-showsInfixBinary d prec op u v = showParen (d > prec) $
+showsInfixBinary d prec op u v = showParen (d >= prec) $
                                  showsPrec prec u .
                                  showString op .
                                  showsPrec prec v
 
 showsUnary :: Show a => Int -> Int -> String -> a -> ShowS
-showsUnary d prec op u = showParen (d > prec) $
+showsUnary d prec op u = showParen (d >= prec) $
                          showString op .
                          showsPrec prec u
 
@@ -118,11 +110,11 @@
   showsPrec d (Negate x) = showsUnary d 7 "-" x
   showsPrec d (Abs x) = showsUnary d 10 "abs" x
   showsPrec d (Signum x) = showsUnary d 10 "signum" x
-  showsPrec _ (FromInteger k) = showString (show k)
+  showsPrec d (FromInteger k) = showParen (d >= 9) $ showString (show k)
 
 instance Show a => Show (Fractionals a) where
   showsPrec d (Div x y) = showsInfixBinary d 7 " / " x y
-  showsPrec _ (FromRational r) = showString $ show (fromRational r :: Double)
+  showsPrec d (FromRational r) = showParen (d >= 9) $ showString $ show (fromRational r :: Double)
 
 instance Show a => Show (Floatings a) where
   showsPrec d (Pow x y) = showsInfixBinary d 8 " ** " x y
@@ -131,6 +123,7 @@
   showsPrec d (Log x)   = showsUnary d 10 "log" x
   showsPrec d (Sin x)   = showsUnary d 10 "sin" x
   showsPrec d (Cos x)   = showsUnary d 10 "cos" x
+  showsPrec d (Tan x)   = showsUnary d 10 "tan" x
   showsPrec d (ASin x)  = showsUnary d 10 "asin" x
   showsPrec d (ATan x)  = showsUnary d 10 "atan" x
   showsPrec d (ACos x)  = showsUnary d 10 "acos" x
@@ -142,8 +135,9 @@
   showsPrec d (ACosh x) = showsUnary d 10 "acosh" x
 
 instance Show a => Show (Expr a) where
-  showsPrec _ (ESym s) = showString (show s)
-  showsPrec _ (EConst x) = showString (show x)
+  showsPrec d (ESym s) = showParen (d > 9) $
+                         showString (show s)
+  showsPrec d (EConst x) = showParen (d >= 9) $ showString (show x)
   showsPrec d (ENum x) = showsPrec d x
   showsPrec d (EFractional x) = showsPrec d x
   showsPrec d (EFloating x) = showsPrec d x
@@ -171,83 +165,57 @@
   (Signum x) == (Signum y) = x == y
   (FromInteger x) == (FromInteger y) = x == y
   _ == _ = False
-  
 
+
 ----------------------------- hashable instances --------------------------
-instance Hashable Sym where
-  hash (Sym name) = hash "Sym" `combine` hash name
-  hash (SymDependent name k s) = hash ("SymDependent", name, k, s)
+instance Hashable Sym
 
 instance Hashable a => Hashable (Nums a) where
-  hash (Mul x y)  = hash "Mul" `combine` hx `combine` hy
+  hashWithSalt s (Mul x y)  = s `hashWithSalt` "Mul" `hashWithSalt` hx `hashWithSalt` hy
     where
       hx' = hash x
       hy' = hash y
       (hx, hy)
         | commutativeMul = (min hx' hy', max hx' hy')
         | otherwise = (hx', hy')
-  hash (Add x y)  = hash "Add" `combine` hx `combine` hy
+  hashWithSalt s (Add x y)  = s `hashWithSalt` "Add" `hashWithSalt` hx `hashWithSalt` hy
     where
       hx' = hash x
       hy' = hash y
       (hx, hy)
         | commutativeAdd = (min hx' hy', max hx' hy')
         | otherwise = (hx', hy')
-  hash (Sub x y)  = hash "Sub" `combine` hash x `combine` hash y
-  hash (Negate x)      = hash "Negate"      `combine` hash x
-  hash (Abs x)         = hash "Abs"         `combine` hash x
-  hash (Signum x)      = hash "Signum"      `combine` hash x
-  hash (FromInteger x) = hash "FromInteger" `combine` hash x
-
-instance Hashable a => Hashable (Fractionals a) where
-  hash (Div x y)  = hash "Div" `combine` hash x `combine` hash y
-  hash (FromRational x) = hash "FromRational" `combine` hash x
+  hashWithSalt s (Sub x y)  = s `hashWithSalt` "Sub" `hashWithSalt` x `hashWithSalt` y
+  hashWithSalt s (Negate x)      = s `hashWithSalt` "Negate"      `hashWithSalt` x
+  hashWithSalt s (Abs x)         = s `hashWithSalt` "Abs"         `hashWithSalt` x
+  hashWithSalt s (Signum x)      = s `hashWithSalt` "Signum"      `hashWithSalt` x
+  hashWithSalt s (FromInteger x) = s `hashWithSalt` "FromInteger" `hashWithSalt` x
 
-instance Hashable a => Hashable (Floatings a) where
-  hash (Pow x y) = hash "Pow" `combine` hash x `combine` hash y
-  hash (LogBase x y) = hash "LogBase" `combine` hash x `combine` hash y
-  hash (Exp x)   = hash "Exp"   `combine` hash x
-  hash (Log x)   = hash "Log"   `combine` hash x
-  hash (Sin x)   = hash "Sin"   `combine` hash x
-  hash (Cos x)   = hash "Cos"   `combine` hash x
-  hash (ASin x)  = hash "ASin"  `combine` hash x
-  hash (ATan x)  = hash "ATan"  `combine` hash x
-  hash (ACos x)  = hash "ACos"  `combine` hash x
-  hash (Sinh x)  = hash "Sinh"  `combine` hash x
-  hash (Cosh x)  = hash "Cosh"  `combine` hash x
-  hash (Tanh x)  = hash "Tanh"  `combine` hash x
-  hash (ASinh x) = hash "ASinh" `combine` hash x
-  hash (ATanh x) = hash "ATanh" `combine` hash x
-  hash (ACosh x) = hash "ACosh" `combine` hash x
+instance Hashable a => Hashable (Fractionals a)
+instance Hashable a => Hashable (Floatings a)
 
 instance Hashable a => Hashable (Expr a) where
-  hash (ESym name)     = hash "ESym"        `combine` hash name
-  hash (EConst x)      = hash "EConst"      `combine` hash x
-  hash (ENum x)        = hash "ENum"        `combine` hash x
-  hash (EFractional x) = hash "EFractional" `combine` hash x
-  hash (EFloating x)   = hash "EFloating"   `combine` hash x
-
---deriving instance Enum a => Enum (Nums a)
---deriving instance Bounded a => Bounded (Nums a)
-
---deriving instance Enum a => Enum (Fractionals a)
---deriving instance Bounded a => Bounded (Fractionals a)
-
---deriving instance Enum a => Enum (Floatings a)
---deriving instance Bounded a => Bounded (Floatings a)
+  hashWithSalt s (ESym name)     = s `hashWithSalt` "ESym"        `hashWithSalt` name
+  hashWithSalt s (EConst x)      = s `hashWithSalt` "EConst"      `hashWithSalt` x
+  hashWithSalt s (ENum x)        = s `hashWithSalt` "ENum"        `hashWithSalt` x
+  hashWithSalt s (EFractional x) = s `hashWithSalt` "EFractional" `hashWithSalt` x
+  hashWithSalt s (EFloating x)   = s `hashWithSalt` "EFloating"   `hashWithSalt` x
 
-fromNeg :: (Num a, Ord a) => Expr a -> Maybe (Expr a)
+fromNeg :: Expr a -> Maybe (Expr a)
 fromNeg (ENum (Negate x)) = Just x
 fromNeg (ENum (FromInteger k))
   | k < 0 = Just (ENum (FromInteger (abs k)))
 fromNeg (EFractional (FromRational r))
   | r < 0 = Just (EFractional (FromRational (abs r)))
-fromNeg (EConst c)
-  | c < 0 = Just (EConst (abs c))
+--fromNeg (EConst c)
+--  | c < 0 = Just (EConst (abs c))
 fromNeg _ = Nothing
 
-
-instance (Num a, Ord a) => Num (Expr a) where
+instance (Eq a, Num a) => Num (Expr a) where
+  (*) x y
+    | isVal 0 x || isVal 0 y = 0
+    | isVal 1 x = y
+    | isVal 1 y = x
   (*) (EConst x) (EConst y) = EConst (x*y)
   (*) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx * ky)
   (*) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx * ry)
@@ -257,16 +225,16 @@
   (*) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx * y
   (*) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx * ry)
   (*) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx * fromInteger ky)
-  (*) x y
-    | isVal 0 x || isVal 0 y = 0
-    | isVal 1 x = y
-    | isVal 1 y = x
   (*) x y = case (fromNeg x, fromNeg y) of
               (Just x', Just y') -> x' * y'
               (Nothing, Just y') -> negate (x  * y')
               (Just x', Nothing) -> negate (x' * y )
               _ -> ENum $ Mul x y
 
+  (+) x y
+    | isVal 0 x = y
+    | isVal 0 y = x
+    | x == negate y = 0
   (+) (EConst x) (EConst y) = EConst (x+y)
   (+) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx + ky)
   (+) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx + ry)
@@ -276,16 +244,16 @@
   (+) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx + y
   (+) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx + ry)
   (+) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx + fromInteger ky)
-  (+) x y
-    | isVal 0 x = y
-    | isVal 0 y = x
-    | x == negate y = 0
   (+) x y = case (fromNeg x, fromNeg y) of
               (Just x', Just y') -> negate (x' + y')
               (Nothing, Just y') -> x  - y'
               (Just x', Nothing) -> y  - x'
               _ -> ENum $ Add x y
 
+  (-) x y
+    | isVal 0 x = negate y
+    | isVal 0 y = x
+    | x == y = 0
   (-) (EConst x) (EConst y) = EConst (x-y)
   (-) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx - ky)
   (-) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx - ry)
@@ -295,10 +263,6 @@
   (-) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx - y
   (-) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx - ry)
   (-) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx - fromInteger ky)
-  (-) x y
-    | isVal 0 x = negate y
-    | isVal 0 y = x
-    | x == y = 0
   (-) x y = case (fromNeg x, fromNeg y) of
               (Just x', Just y') -> y' - x' -- (-x) - (-y) == y - x
               (Nothing, Just y') -> x + y' -- (x) - (-y) == x + y
@@ -325,7 +289,12 @@
 
   fromInteger = ENum . FromInteger
 
-instance (Fractional a, Ord a) => Fractional (Expr a) where
+instance (Eq a, Fractional a) => Fractional (Expr a) where
+  (/) x y
+    | isVal 0 y = error "Fractional (Expr a) divide by zero"
+    | isVal 0 x = 0
+    | isVal 1 y = x
+    | x == y = 1
   (/) (EConst x) (EConst y) = EConst (x/y)
   (/) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = EFractional $ FromRational (kx % ky)
   (/) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx / ry)
@@ -335,10 +304,6 @@
   (/) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx / y
   (/) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx / ry)
   (/) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx / fromInteger ky)
-  (/) x y
-    | isVal 0 y = error "Fractional (Expr a) divide by zero"
-    | isVal 0 x = 0
-    | isVal 1 y = x
   (/) x y = case (fromNeg x, fromNeg y) of
               (Just x', Just y') -> x' / y'
               (Nothing, Just y') -> negate (x  / y')
@@ -347,14 +312,18 @@
 
   fromRational = EFractional . FromRational
 
-instance (Floating a, Ord a) => Floating (Expr a) where
+instance (Eq a, Floating a) => Floating (Expr a) where
   pi          = EConst pi
+  x ** 1      = x
+  0 ** 0      = error "Expr: 0 ** 0 indeterminate"
+  _ ** 0      = 1
   x ** y      = EFloating $ Pow x y
   logBase x y = EFloating $ LogBase x y
   exp         = applyFloatingUn (  exp,   Exp)
   log         = applyFloatingUn (  log,   Log)
   sin         = applyFloatingUn (  sin,   Sin)
   cos         = applyFloatingUn (  cos,   Cos)
+  tan         = applyFloatingUn (  tan,   Tan)
   asin        = applyFloatingUn ( asin,  ASin)
   atan        = applyFloatingUn ( atan,  ATan)
   acos        = applyFloatingUn ( acos,  ACos)
@@ -378,108 +347,43 @@
   GNum :: Num a => Nums b -> GExpr a b
   GFractional :: Fractional a => Fractionals b -> GExpr a b
   GFloating :: Floating a => Floatings b -> GExpr a b
-deriving instance (Ord a, Ord b) => Ord (GExpr a b)
 
--- you might use this to use Expr's nice Show instance
-gexprToExpr :: (b -> Expr a) -> GExpr a b -> Expr a
-gexprToExpr _ (GSym s@(Sym _)) = ESym s
-gexprToExpr _ (GSym sd@(SymDependent _ _ _)) = ESym sd
-gexprToExpr _ (GConst c) = EConst c
-gexprToExpr f (GNum (Mul x y))               = ENum (Mul (f x) (f y))
-gexprToExpr f (GNum (Add x y))               = ENum (Add (f x) (f y))
-gexprToExpr f (GNum (Sub x y))               = ENum (Sub (f x) (f y))
-gexprToExpr f (GNum (Negate x))              = ENum (Negate (f x))
-gexprToExpr f (GNum (Abs x))                 = ENum (Abs (f x))
-gexprToExpr f (GNum (Signum x))              = ENum (Signum (f x))
-gexprToExpr _ (GNum (FromInteger x))         = ENum (FromInteger x)
-gexprToExpr f (GFractional (Div x y))        = EFractional (Div (f x) (f y))
-gexprToExpr _ (GFractional (FromRational x)) = EFractional (FromRational x)
-gexprToExpr f (GFloating (Pow x y))          = EFloating (Pow (f x) (f y))
-gexprToExpr f (GFloating (LogBase x y))      = EFloating (LogBase (f x) (f y))
-gexprToExpr f (GFloating (Exp x))            = EFloating (Exp   (f x))
-gexprToExpr f (GFloating (Log x))            = EFloating (Log   (f x))
-gexprToExpr f (GFloating (Sin x))            = EFloating (Sin   (f x))
-gexprToExpr f (GFloating (Cos x))            = EFloating (Cos   (f x))
-gexprToExpr f (GFloating (ASin x))           = EFloating (ASin  (f x))
-gexprToExpr f (GFloating (ATan x))           = EFloating (ATan  (f x))
-gexprToExpr f (GFloating (ACos x))           = EFloating (ACos  (f x))
-gexprToExpr f (GFloating (Sinh x))           = EFloating (Sinh  (f x))
-gexprToExpr f (GFloating (Cosh x))           = EFloating (Cosh  (f x))
-gexprToExpr f (GFloating (Tanh x))           = EFloating (Tanh  (f x))
-gexprToExpr f (GFloating (ASinh x))          = EFloating (ASinh (f x))
-gexprToExpr f (GFloating (ATanh x))          = EFloating (ATanh (f x))
-gexprToExpr f (GFloating (ACosh x))          = EFloating (ACosh (f x))
-
-getParents :: GExpr a b -> [b]
-getParents (GSym _)                       = []
-getParents (GConst _)                     = []
-getParents (GNum (Mul x y))               = [x,y]
-getParents (GNum (Add x y))               = [x,y]
-getParents (GNum (Sub x y))               = [x,y]
-getParents (GNum (Negate x))              = [x]
-getParents (GNum (Abs x))                 = [x]
-getParents (GNum (Signum x))              = [x]
-getParents (GNum (FromInteger _))         = []
-getParents (GFractional (Div x y))        = [x,y]
-getParents (GFractional (FromRational _)) = []
-getParents (GFloating (Pow x y))          = [x,y]
-getParents (GFloating (LogBase x y))      = [x,y]
-getParents (GFloating (Exp x))            = [x]
-getParents (GFloating (Log x))            = [x]
-getParents (GFloating (Sin x))            = [x]
-getParents (GFloating (Cos x))            = [x]
-getParents (GFloating (ASin x))           = [x]
-getParents (GFloating (ATan x))           = [x]
-getParents (GFloating (ACos x))           = [x]
-getParents (GFloating (Sinh x))           = [x]
-getParents (GFloating (Cosh x))           = [x]
-getParents (GFloating (Tanh x))           = [x]
-getParents (GFloating (ASinh x))          = [x]
-getParents (GFloating (ATanh x))          = [x]
-getParents (GFloating (ACosh x))          = [x]
+deriving instance (Show a, Show b) => Show (GExpr a b)
+instance Functor (GExpr a) where
+  fmap _ (GSym s) = GSym s
+  fmap _ (GConst c) = GConst c
+  fmap f (GNum nums) = GNum (fmap f nums)
+  fmap f (GFractional fracs) = GFractional (fmap f fracs)
+  fmap f (GFloating floatings) = GFloating (fmap f floatings)
 
-instance (Show a, Show b) => Show (GExpr a b) where
-  show = show . (gexprToExpr (\x -> ESym (Sym ("{" ++ show x ++ "}"))))
-  
-deriving instance (Eq a, Eq b) => Eq (GExpr a b)
+instance F.Foldable (GExpr a) where
+  foldMap _ (GSym _) = mempty
+  foldMap _ (GConst _) = mempty
+  foldMap f (GNum nums) = F.foldMap f nums
+  foldMap f (GFractional fracs) = F.foldMap f fracs
+  foldMap f (GFloating floatings) = F.foldMap f floatings
 
-instance (Hashable a, Hashable b) => Hashable (GExpr a b) where
-  hash (GSym name)     = hash "GSym"        `combine` hash name
-  hash (GConst x)      = hash "GConst"      `combine` hash x
-  hash (GNum x)        = hash "GNum"        `combine` hash x
-  hash (GFractional x) = hash "GFractional" `combine` hash x
-  hash (GFloating x)   = hash "GFloating"   `combine` hash x
+  foldr _ z (GSym _) = z
+  foldr _ z (GConst _) = z
+  foldr f z (GNum nums) = F.foldr f z nums
+  foldr f z (GFractional fracs) = F.foldr f z fracs
+  foldr f z (GFloating floatings) = F.foldr f z floatings
 
-instance MuRef (Expr a) where
-  type DeRef (Expr a) = GExpr a
-  mapDeRef _ (ESym name) = pure (GSym name)
-  mapDeRef _ (EConst c)  = pure (GConst c)
-  mapDeRef f (ENum (Mul x y)) = GNum <$> (Mul <$> (f x) <*> (f y))
-  mapDeRef f (ENum (Add x y)) = GNum <$> (Add <$> (f x) <*> (f y))
-  mapDeRef f (ENum (Sub x y)) = GNum <$> (Sub <$> (f x) <*> (f y))
-  mapDeRef f (ENum (Negate x)) = GNum <$> (Negate <$> (f x))
-  mapDeRef f (ENum (Abs x)) = GNum <$> (Negate <$> (f x))
-  mapDeRef f (ENum (Signum x)) = GNum <$> (Signum <$> (f x))
-  mapDeRef _ (ENum (FromInteger k)) = pure $ GNum (FromInteger k)
+instance T.Traversable (GExpr a) where
+  traverse _ (GSym s) = pure (GSym s)
+  traverse _ (GConst c) = pure (GConst c)
+  traverse f (GNum nums) = GNum <$> T.traverse f nums
+  traverse f (GFractional fracs) = GFractional <$> T.traverse f fracs
+  traverse f (GFloating floatings) = GFloating <$> T.traverse f floatings
 
-  mapDeRef f (EFractional (Div x y)) = GFractional <$> (Div <$> (f x) <*> (f y))
-  mapDeRef _ (EFractional (FromRational x)) = pure $ GFractional (FromRational x)
+deriving instance (Eq a, Eq b) => Eq (GExpr a b)
 
-  mapDeRef f (EFloating (Pow x y))     = GFloating <$> (Pow <$> (f x) <*> (f y))
-  mapDeRef f (EFloating (LogBase x y)) = GFloating <$> (LogBase <$> (f x) <*> (f y))
-  mapDeRef f (EFloating (Exp   x))     = GFloating <$> (Exp   <$> (f x))
-  mapDeRef f (EFloating (Log   x))     = GFloating <$> (Log   <$> (f x))
-  mapDeRef f (EFloating (Sin   x))     = GFloating <$> (Sin   <$> (f x))
-  mapDeRef f (EFloating (Cos   x))     = GFloating <$> (Cos   <$> (f x))
-  mapDeRef f (EFloating (ASin  x))     = GFloating <$> (ASin  <$> (f x))
-  mapDeRef f (EFloating (ATan  x))     = GFloating <$> (ATan  <$> (f x))
-  mapDeRef f (EFloating (ACos  x))     = GFloating <$> (ACos  <$> (f x))
-  mapDeRef f (EFloating (Sinh  x))     = GFloating <$> (Sinh  <$> (f x))
-  mapDeRef f (EFloating (Cosh  x))     = GFloating <$> (Cosh  <$> (f x))
-  mapDeRef f (EFloating (Tanh  x))     = GFloating <$> (Tanh  <$> (f x))
-  mapDeRef f (EFloating (ASinh x))     = GFloating <$> (ASinh <$> (f x))
-  mapDeRef f (EFloating (ATanh x))     = GFloating <$> (ATanh <$> (f x))
-  mapDeRef f (EFloating (ACosh x))     = GFloating <$> (ACosh <$> (f x))
+instance (Hashable a, Hashable b) => Hashable (GExpr a b) where
+  hashWithSalt s (GSym name)     = s `hashWithSalt` "GSym"        `hashWithSalt` name
+  hashWithSalt s (GConst x)      = s `hashWithSalt` "GConst"      `hashWithSalt` x
+  hashWithSalt s (GNum x)        = s `hashWithSalt` "GNum"        `hashWithSalt` x
+  hashWithSalt s (GFractional x) = s `hashWithSalt` "GFractional" `hashWithSalt` x
+  hashWithSalt s (GFloating x)   = s `hashWithSalt` "GFloating"   `hashWithSalt` x
 
 substitute :: (Ord a, Hashable a, Show a) => Expr a -> [(Expr a, Expr a)] -> Expr a
 substitute expr subList
@@ -491,26 +395,27 @@
     nonSymInputs = filter (not . isSym . fst) subList
     lookup' e = let hm = HM.fromList subList in
       HM.lookupDefault e e hm
-    
+
     subs e@(ESym _) = lookup' e
     subs e@(EConst _) = e
     subs e@(ENum (FromInteger _)) = e
     subs e@(EFractional (FromRational _)) = e
-    subs (ENum (Mul x y)) = (subs x) * (subs y)
-    subs (ENum (Add x y)) = (subs x) + (subs y)
-    subs (ENum (Sub x y)) = (subs x) - (subs y)
+    subs (ENum (Mul x y)) = subs x * subs y
+    subs (ENum (Add x y)) = subs x + subs y
+    subs (ENum (Sub x y)) = subs x - subs y
     subs (ENum (Negate x)) = negate (subs x)
     subs (ENum (Abs x))    = abs (subs x)
     subs (ENum (Signum x)) = signum (subs x)
-    
-    subs (EFractional (Div x y)) = (subs x) / (subs y)
-    
-    subs (EFloating (Pow x y))     = (subs x) ** (subs y)
+
+    subs (EFractional (Div x y)) = subs x / subs y
+
+    subs (EFloating (Pow x y))     = subs x ** subs y
     subs (EFloating (LogBase x y)) = logBase (subs x) (subs y)
     subs (EFloating (Exp   x))     = exp   (subs x)
     subs (EFloating (Log   x))     = log   (subs x)
     subs (EFloating (Sin   x))     = sin   (subs x)
     subs (EFloating (Cos   x))     = cos   (subs x)
+    subs (EFloating (Tan   x))     = tan   (subs x)
     subs (EFloating (ASin  x))     = asin  (subs x)
     subs (EFloating (ATan  x))     = atan  (subs x)
     subs (EFloating (ACos  x))     = acos  (subs x)
@@ -521,7 +426,8 @@
     subs (EFloating (ATanh x))     = atanh (subs x)
     subs (EFloating (ACosh x))     = acosh (subs x)
 
--- | this substitute is sketchy because it doesn't perform simplifications that are often assumed to be done
+-- | this substitute is sketchy because it doesn't perform simplifications
+--   that are often assumed to be done
 sketchySubstitute :: (Eq a, Hashable a, Show a) => Expr a -> [(Expr a, Expr a)] -> Expr a
 sketchySubstitute expr subList
   | nonSymInputs /= [] = error $ "substitute got non-ESym input: " ++ show nonSymInputs
@@ -532,7 +438,7 @@
     nonSymInputs = filter (not . isSym . fst) subList
     lookup' e = let hm = HM.fromList subList in
       HM.lookupDefault e e hm
-    
+
     subs e@(ESym _) = lookup' e
     subs e@(EConst _)  = e
     subs e@(ENum (FromInteger _)) = e
@@ -543,15 +449,16 @@
     subs (ENum (Negate x)) = ENum (Negate (subs x))
     subs (ENum (Abs x)) = ENum (Negate (subs x))
     subs (ENum (Signum x)) = ENum (Signum (subs x))
-  
+
     subs (EFractional (Div x y)) = EFractional (Div (subs x) (subs y))
-  
+
     subs (EFloating (Pow x y))     = EFloating (Pow (subs x) (subs y))
     subs (EFloating (LogBase x y)) = EFloating (LogBase (subs x) (subs y))
     subs (EFloating (Exp   x))     = EFloating (Exp   (subs x))
     subs (EFloating (Log   x))     = EFloating (Log   (subs x))
     subs (EFloating (Sin   x))     = EFloating (Sin   (subs x))
     subs (EFloating (Cos   x))     = EFloating (Cos   (subs x))
+    subs (EFloating (Tan   x))     = EFloating (Tan   (subs x))
     subs (EFloating (ASin  x))     = EFloating (ASin  (subs x))
     subs (EFloating (ATan  x))     = EFloating (ATan  (subs x))
     subs (EFloating (ACos  x))     = EFloating (ACos  (subs x))
@@ -583,6 +490,7 @@
 foldExpr f acc (EFloating (Log x))       = foldExpr f acc x
 foldExpr f acc (EFloating (Sin x))       = foldExpr f acc x
 foldExpr f acc (EFloating (Cos x))       = foldExpr f acc x
+foldExpr f acc (EFloating (Tan x))       = foldExpr f acc x
 foldExpr f acc (EFloating (ASin x))      = foldExpr f acc x
 foldExpr f acc (EFloating (ATan x))      = foldExpr f acc x
 foldExpr f acc (EFloating (ACos x))      = foldExpr f acc x
@@ -618,6 +526,7 @@
 isVal v (ENum (FromInteger k)) = v == fromInteger k
 isVal v (EFractional (FromRational r)) = v == fromRational r
 isVal _ _ = False
+{-# INLINE isVal #-}
 
 -- | if the expression is a constant, a fromInteger, or a fromRational, return the constant part
 --   otherwise return nothing
@@ -629,7 +538,7 @@
 
 -- | Separate nonlinear and linear parts of an expression
 --   @extractLinearPart (fNonLin(x)+a*x) x == (fNonLin(x), a)
-extractLinearPart :: (Num a, Ord a, Show a) => Expr a -> Expr a -> (Expr a, a)
+extractLinearPart :: (Num a, Eq a, Show a) => Expr a -> Expr a -> (Expr a, a)
 extractLinearPart e@(EConst _) _ = (e,0)
 extractLinearPart e@(ENum (FromInteger _)) _ = (e,0)
 extractLinearPart e@(EFractional (FromRational _)) _ = (e,0)
@@ -649,13 +558,13 @@
     (xNonlin,xLin) = extractLinearPart x arg
 extractLinearPart e@(ENum (Mul x y)) arg = case (getConst x, getConst y) of
   (Nothing,Nothing) -> (e,0)
-  (Just cx, Nothing) -> let (yNl,yL) = extractLinearPart y arg in ((EConst cx)*yNl,cx*yL)
-  (Nothing, Just cy) -> let (xNl,xL) = extractLinearPart x arg in (xNl*(EConst cy),xL*cy)
+  (Just cx, Nothing) -> let (yNl,yL) = extractLinearPart y arg in (EConst cx * yNl,cx*yL)
+  (Nothing, Just cy) -> let (xNl,xL) = extractLinearPart x arg in (xNl * EConst cy,xL*cy)
   _ -> error $ "extractLinearPart got ENum (Mul x y) where x and y are both constants\n"++
        "x: " ++ show x ++ "\ny: " ++ show y
 extractLinearPart e@(EFractional (Div x y)) arg = case getConst y of
   Nothing -> (e,0)
-  Just cy -> let (xNl,xL) = extractLinearPart x arg in (xNl/(EConst cy),xL/cy)
+  Just cy -> let (xNl,xL) = extractLinearPart x arg in (xNl/EConst cy,xL/cy)
 extractLinearPart e@(ENum (Abs _))    _ = (e,0)
 extractLinearPart e@(ENum (Signum _)) _ = (e,0)
 extractLinearPart e@(EFloating (Pow _ _)) _ = (e,0)
@@ -664,6 +573,7 @@
 extractLinearPart e@(EFloating (Log _))   _ = (e,0)
 extractLinearPart e@(EFloating (Sin _))   _ = (e,0)
 extractLinearPart e@(EFloating (Cos _))   _ = (e,0)
+extractLinearPart e@(EFloating (Tan _))   _ = (e,0)
 extractLinearPart e@(EFloating (ASin _))  _ = (e,0)
 extractLinearPart e@(EFloating (ATan _))  _ = (e,0)
 extractLinearPart e@(EFloating (ACos _))  _ = (e,0)
@@ -685,7 +595,7 @@
 --            | Abs a
 --            | Signum a
 --            | FromInteger Integer
-  
+
 --instance Arbitrary a => Arbitrary (Expr a) where
 --   arbitrary = oneof [arbConst, arbUnary, arbBinary]
 --
diff --git a/Dvda/FunGraph.hs b/Dvda/FunGraph.hs
deleted file mode 100644
--- a/Dvda/FunGraph.hs
+++ /dev/null
@@ -1,164 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language TypeOperators #-}
-{-# Language TypeFamilies #-}
-{-# Language FlexibleInstances #-}
-
-module Dvda.FunGraph ( FunGraph
-                     , ToFunGraph
-                     , NumT
-                     , (:*)(..)
-                     , MVS(..)
-                     , toFunGraph
-                     , countNodes
-                     , fgInputs
-                     , fgOutputs
-                     , fgLookupGExpr
-                     , fgReified
-                     , topSort
---                     , fgGraph
-                     , nodelistToFunGraph
-                     , exprsToFunGraph
-                     ) where
-
-import Control.Applicative
-import Data.Foldable ( Foldable )
-import qualified Data.Foldable as F
-import qualified Data.Graph as Graph
-import Data.Hashable ( Hashable )
-import qualified Data.HashSet as HS
-import Data.Traversable ( Traversable )
-import qualified Data.Traversable as T
-
-import Dvda.Expr
-import Dvda.Reify ( ReifyGraph(..), reifyGraphs )
-
-data FunGraph a = FunGraph { fgGraph :: Graph.Graph
-                           , fgInputs :: [MVS (GExpr a Int)]
-                           , fgOutputs :: [MVS Int]
-                           , fgReified :: [(Int, GExpr a Int)]
-                           , fgLookupGExpr :: Int -> Maybe (GExpr a Int)
-                           , fgVertexFromKey :: Int -> Maybe Int
-                           , fgNodeFromVertex :: Int -> (GExpr a Int, Int, [Int])
-                           }
-
-instance Show a => Show (FunGraph a) where
-  show fg = "FunGraph\ninputs:\n" ++ show (fgInputs fg) ++ "\noutputs:\n" ++ show (fgOutputs fg) ++ "\ngraph:\n" ++ show (fgGraph fg)
-
----- | matrix or vector or scalar
-data MVS a = Mat [[a]] | Vec [a] | Sca a deriving Show
-
-instance Functor MVS where
-  fmap f (Sca x)  = Sca (f x)
-  fmap f (Vec xs) = Vec (map f xs)
-  fmap f (Mat xs) = Mat (map (map f) xs)
-
-instance Foldable MVS where
-  foldr f x0 (Sca x)  = foldr f x0 [x]
-  foldr f x0 (Vec xs) = foldr f x0 xs
-  foldr f x0 (Mat xs) = foldr f x0 (concat xs)
-
-instance Traversable MVS where
-  traverse f (Sca x)  = Sca <$> f x
-  traverse f (Vec xs) = Vec <$> T.traverse f xs
-  traverse f (Mat xs) = Mat <$> T.traverse (T.traverse f) xs
-
-class ToFunGraph a where
-  type NumT a
-  toMVSList :: a -> [MVS (Expr (NumT a))]
-instance ToFunGraph (Expr a) where
-  type NumT (Expr a) = a
-  toMVSList x = [Sca x]
-instance ToFunGraph [Expr a] where
-  type NumT [Expr a] = NumT (Expr a)
-  toMVSList x = [Vec x]
-instance ToFunGraph [[Expr a]] where
-  type NumT [[Expr a]] = NumT [Expr a]
-  toMVSList x = [Mat x]
-
-data a :* b = a :* b deriving Show
-infixr 6 :*
-instance (ToFunGraph a, ToFunGraph b, NumT a ~ NumT b) => ToFunGraph (a :* b) where
-  type NumT (a :* b) = NumT a
-  toMVSList (x :* y) = toMVSList x ++ toMVSList y
-
--- | find any symbols which are parents of outputs, but are not supplied by the user
-detectMissingInputs :: (Eq a, Hashable a, Show a) => [MVS (Expr a)] -> [(Int,GExpr a Int)] -> [GExpr a Int]
-detectMissingInputs exprs gr = HS.toList $ HS.difference allGraphInputs allUserInputs
-  where
-    allUserInputs = let f (ESym name) acc = (GSym name):acc
-                        f _ e = error $ "detectMissingInputs given non-ESym input \"" ++ show e ++ "\""
-                    in HS.fromList $ foldr f [] (concatMap F.toList exprs)
-
-    allGraphInputs = let f (_,(GSym name)) acc = (GSym name):acc
-                         f _ acc = acc
-                     in HS.fromList $ foldr f [] gr
-
--- | if the same input symbol (like ESym "x") is given at two different places throw an exception
-findConflictingInputs :: (Eq a, Hashable a, Show a) => [MVS (Expr a)] -> [Expr a]
-findConflictingInputs exprs = HS.toList redundant
-  where
-    redundant = snd $ foldl f (HS.empty, HS.empty) (concatMap F.toList exprs)
-      where
-        f (knownExprs, redundantExprs) expr@(ESym _)
-          | HS.member expr knownExprs = (knownExprs, HS.insert expr redundantExprs)
-          | otherwise = (HS.insert expr knownExprs, redundantExprs)
-        f _ e = error $ "findConflictingInputs saw non-ESym input \"" ++ show e ++ "\""
-
-
--- | Take inputs and outputs which are of classes ToFunGraph (heterogenous lists of @Expr a@)
---   and traverse the outputs reifying all expressions and creating a hashmap of StableNames (stable pointers).
---   Once the hashmap is created, lookup the provided inputs and return a FunGraph which contains an
---   expression graph, input/output indices, and other useful functions. StableNames is non-deterministic
---   so this function may return graphs with more or fewer CSE's eliminated.
---   If CSE is then performed on the graph, the result is deterministic.
-toFunGraph :: (Eq a, Hashable a, Show a, ToFunGraph b, ToFunGraph c, NumT b ~ a, NumT c ~ a)
-              => b -> c -> IO (FunGraph a)
-toFunGraph inputs outputs = mvsToFunGraph (toMVSList inputs) (toMVSList outputs)
-
-mvsToFunGraph :: (Eq a, Hashable a, Show a) => [MVS (Expr a)] -> [MVS (Expr a)] -> IO (FunGraph a)
-mvsToFunGraph inputMVSExprs outputMVSExprs = do
-  -- reify the outputs
-  (ReifyGraph rgr, outputMVSIndices) <- reifyGraphs outputMVSExprs
-  let fg = nodelistToFunGraph rgr inputMVSGExprs outputMVSIndices
-      inputMVSGExprs = map (fmap f) inputMVSExprs
-        where
-          f (ESym name) = (GSym name)
-          f x = error $ "ERROR: mvsToFunGraph given non-ESym input \"" ++ show x ++ "\""
-  return $ case (detectMissingInputs inputMVSExprs rgr, findConflictingInputs inputMVSExprs) of
-    ([],[]) -> fg
-    (xs,[]) -> error $ "mvsToFunGraph found inputs that were not provided by the user: " ++ show xs
-    ( _,xs) -> error $ "mvsToFunGraph found idential inputs set more than once: " ++ show xs
-
-nodelistToFunGraph :: [(Int,GExpr a Int)] -> [MVS (GExpr a Int)] -> [MVS Int] -> FunGraph a
-nodelistToFunGraph rgr inputMVSIndices outputMVSIndices =
-  FunGraph { fgGraph = gr
-           , fgInputs = inputMVSIndices
-           , fgOutputs = outputMVSIndices
-           , fgLookupGExpr = lookupG
-           , fgReified = rgr
-           , fgVertexFromKey = lookupKey
-           , fgNodeFromVertex = lookupVertex
-           }
-  where
-    -- make sure all the inputs are symbolic, and find their indices in the Expr graph
-    (gr, lookupVertex, lookupKey) = Graph.graphFromEdges $ map (\(k,gexpr) -> (gexpr, k, getParents gexpr)) rgr
-    lookupG k = (\(g,_,_) -> g) <$> lookupVertex <$> lookupKey k
-
-
----------------------------------- utilities -----------------------------
-countNodes :: FunGraph a -> Int
-countNodes = length . Graph.vertices . fgGraph
-
-topSort :: FunGraph a -> [Int]
-topSort fg = map ((\(_,k,_) -> k) . (fgNodeFromVertex fg)) $ Graph.topSort (fgGraph fg)
-
--- | make a FunGraph out of outputs, automatically detecting the proper inputs
-exprsToFunGraph :: (Eq a, Show a, Hashable a) => [Expr a] -> IO (FunGraph a)
-exprsToFunGraph outputs = do
-  let getSyms :: [Expr a] -> [Sym]
-      getSyms exprs = HS.toList $ foldr (\acc expr -> foldExpr f expr acc) HS.empty exprs
-        where
-          f (ESym s) hs = HS.insert s hs
-          f _ hs = hs
-      inputs = map ESym $ getSyms outputs
-  toFunGraph inputs outputs
diff --git a/Dvda/MultipleShooting/CoctaveTemplates.hs b/Dvda/MultipleShooting/CoctaveTemplates.hs
deleted file mode 100644
--- a/Dvda/MultipleShooting/CoctaveTemplates.hs
+++ /dev/null
@@ -1,122 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.MultipleShooting.CoctaveTemplates ( writeMexAll
-                                              , writeSetupSource
-                                              , writeUnstructConsts
-                                              , writeToStruct
-                                              , writeUnstruct
-                                              , writePlot
-                                              )where
-
-import Data.Maybe ( fromMaybe )
-import Data.Hashable ( Hashable )
-import Data.List ( elemIndex, transpose )
-
-import Dvda.Expr ( Expr(..), Sym(..) )
-import Dvda.HashMap ( HashMap )
-import qualified Dvda.HashMap as HM
-
-writeMexAll :: String -> String
-writeMexAll name = unlines $ map f ["time", "outputs", "sim", "cost", "constraints"]
-  where
-    f x = "tic\nfprintf('mexing " ++ file ++ "...  ')\n"++"mex " ++ file ++ "\nt1 = toc;\nfprintf('finished in %.2f seconds\\n', t1)"
-      where
-        file = name ++ "_" ++ x ++ ".c"
-
-
-writeSetupSource :: Show a => String -> [Expr a] -> [a] -> [a] -> String
-writeSetupSource name dvs lbs ubs =
-  unlines $
-  [ "function [x0, Aineq, bineq, Aeq, beq, lb, ub] = "++ name ++"_setup()"
-  , ""
-  , "x0 = zeros(" ++ show (length dvs) ++ ",1);"
-  , "Aineq = [];"
-  , "bineq = [];"
-  , "Aeq = [];"
-  , "beq = [];"
-  , "lb = " ++ show lbs ++ "';"
-  , "ub = " ++ show ubs ++ "';"
-  ]
-
-
--- take nice matlab structs and return vector of design constants
-writeUnstructConsts :: Eq a => String -> [Expr a] -> String
-writeUnstructConsts name constants =
-  unlines $
-  [ "function constants = " ++ name ++ "_unstructConstants(constStruct)\n"
-  , "constants = zeros(" ++ show (length constants) ++ ", 1);"
-  , ""
-  , concatMap fromConst constants
-  ]
-  where
-    readName e = case e of
-      ESym (Sym nm) -> nm
-      _ -> error "const not ESym Sym"
-    fromConst e = "constants(" ++ show (1 + (fromJustErr "fromConst error" $ e `elemIndex` constants)) ++ ") = constStruct." ++ readName e ++ ";\n"
-
-
----- take vector of design variables and vector of constants and return nice matlab struct
-writeToStruct :: (Eq a, Show a, Hashable a)
-                 => String -> [Expr a] -> [Expr a] -> [Expr a] -> HashMap String [Expr a] -> String
-writeToStruct name dvs params constants outputMap =
-  unlines $
-  ["function ret = " ++ name ++ "_struct(designVars,constants)"
-  , ""
-  , "ret.time = " ++ name ++ "_time(designVars, constants);"
-  , "outs = " ++ name ++ "_outputs(designVars, constants);"
-  , concat $ zipWith (\name' k -> "ret." ++name'++ " = outs("++show k++",:);\n") (HM.keys outputMap) [(1::Int)..]
-  ] ++
-  toStruct dvs "designVars" (map show params) (map (\x -> [x]) params) ++
-  toStruct constants "constants" (map show constants) (map (\x -> [x]) constants)
-    where
-      dvsToIdx dvs' = (fromJustErr "toStruct error") . (flip HM.lookup (HM.fromList (zip dvs' [(1::Int)..])))
-
-      toStruct dvs' nm = zipWith (\name' vars -> "ret." ++ name' ++ " = " ++ nm ++ "(" ++ show (map (dvsToIdx dvs') vars) ++ ");\n")
-
--- take nice matlab structs and return vector of design variables
-writeUnstruct :: (Eq a, Show a)
-                 => String
-                 -> [Expr a] -> [Expr a]
-                 -> [Expr a] -> [[Expr a]]
-                 -> [Expr a] -> [[Expr a]]
-                 -> String
-writeUnstruct name dvs params states allStates actions allActions =
-  unlines $
-  [ "function dvs = " ++ name ++ "_unstruct(dvStruct)\n"
-  , "dvs = zeros(" ++ show (length dvs) ++ ", 1);"
-  , ""
-  , concatMap fromParam params
-  , concat $ zipWith fromXU states  (transpose allStates)
-  , concat $ zipWith fromXU actions (transpose allActions)
-  ]
-  where
-    dvIdx e = fromMaybe (error $ "dvIdx error - " ++ show e ++ " is not a design variable")
-              (e `elemIndex` dvs)
-    fromParam e = "dvs(" ++ show (1 + dvIdx e) ++ ") = dvStruct." ++ show e ++ ";\n"
-    fromXU e es =
-      "dvs(" ++ show (map ((1 +) . dvIdx) es) ++ ") = dvStruct." ++ show e ++ ";\n"
-
-writePlot :: String -> HashMap String [Expr a] -> String
-writePlot name outputMap =
-  unlines $
-  [ "function " ++ name ++ "_plot(designVars, constants)\n"
-  , "x = " ++ name ++ "_struct(designVars, constants);\n"
-  , init $ unlines $ zipWith f (HM.keys outputMap) [(1::Int)..]
-  ]
-  where
-    rows = ceiling $ sqrt $ (fromIntegral ::Int -> Double) $ HM.size outputMap
-    cols = (HM.size outputMap `div` rows) + 1
-    f name' k = unlines $
-                [ "subplot(" ++ show rows ++ "," ++ show cols ++ ","++show k++");"
-                , "plot( x.time, x." ++ name' ++ " );"
-                , "xlabel('time');"
-                , "ylabel('" ++ name'' ++ "');"
-                , "title('"  ++ name'' ++ "');"
-                ]
-      where
-        name'' = foldl (\acc x -> if x == '_' then acc ++ "\\_" else acc ++ [x]) "" name'
-
-
-fromJustErr :: String -> Maybe a -> a
-fromJustErr _ (Just x) = x
-fromJustErr message Nothing = error $ "fromJustErr got Nothing, message: \"" ++ message ++ "\""
diff --git a/Dvda/MultipleShooting/MSCoctave.hs b/Dvda/MultipleShooting/MSCoctave.hs
deleted file mode 100644
--- a/Dvda/MultipleShooting/MSCoctave.hs
+++ /dev/null
@@ -1,283 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.MultipleShooting.MSCoctave ( msCoctave
-                                       , run
-                                       ) where
-
-import qualified Control.Monad.State as State
-import Data.Hashable ( Hashable )
-import qualified Data.HashSet as HS
-import Data.List ( zipWith6 )
-import Data.Maybe ( fromMaybe )
-
-import Dvda.AD ( rad )
-import Dvda.CGen  ( showMex )
-import Dvda.CSE ( cse )
-import Dvda.Codegen.WriteFile ( writeSourceFile )
-import Dvda.Expr ( Expr(..), sym, substitute )
-import Dvda.FunGraph ( (:*)(..), toFunGraph, countNodes )
-import Dvda.HashMap ( HashMap )
-import qualified Dvda.HashMap as HM
-import Dvda.MultipleShooting.CoctaveTemplates
-import Dvda.MultipleShooting.MSMonad
-import Dvda.MultipleShooting.Types
-
-{-
-    min f(x) st:
-    
-    c(x) <= 0
-    ceq(x) == 0
-    A*x <= b
-    Aeq*x == beq
-    lb <= x <= ub
--}
-type Integrator a = [Expr Double]
-                   -> [Expr Double]
-                   -> [Expr Double]
-                   -> [Expr Double]
-                   -> ([Expr Double]
-                       -> [Expr Double] -> [Expr Double])
-                   -> Expr Double
-                   -> [Expr Double]
-
--- take user provided bounds and make sure they're complete
--- return functions which will lookup bounds on given state/action @ timestep, and given param
-setupBounds :: (Eq a, Hashable a, Show a)
-               => [(Expr a, (a,a, BCTime))]
-               -> Int
-               -> (Expr a -> Int -> (a,a), Expr a -> (a,a))
-setupBounds userBounds nSteps = (lookupAll, lookupParam)
-  where
-    lookupAll x k
-      | k >= nSteps = error "don't ask for bounds at timestep >= number of total timesteps"
-      | otherwise = case HM.lookup (x,k) specificTimestepBounds of
-        Just bnd -> bnd
-        Nothing -> case HM.lookup x everyTimestepBounds of
-          Just bnd -> bnd
-          Nothing -> error $ "need to set bounds for \"" ++ show x ++ "\" at timestep " ++ show k
-
-    lookupParam x = case HM.lookup x everyTimestepBounds of
-        Just bnd -> bnd
-        Nothing -> error $ "need to set bounds for \"" ++ show x ++ "\""
-
-    -- bounds set at only one timestep
---    everyTimestepBounds :: HashMap (Expr a) (a,a)
-    everyTimestepBounds = let
-      everyTS (e,(lb,ub,ALWAYS)) = [(e,(lb,ub))]
-      everyTS _ = []
-      f (e,lbub) hm =
-        if HM.member e hm
-        then error $ "you set bounds twice for \"" ++ show e ++ "\""
-        else HM.insert e lbub hm
-      in foldr f HM.empty $ concatMap everyTS userBounds
-
-    -- bounds set at specific timestep
---    specificTimestepBounds :: HashMap (Expr a, Int) (a,a)
-    specificTimestepBounds = let
-      specificTS (e,(lb,ub,TIMESTEP k)) = [((e,k),(lb,ub))]
-      specificTS _ = []
-      f (e,lbub) hm =
-        if HM.member e hm
-        then error $ "you set bounds twice for \"" ++ show e ++ "\""
-        else HM.insert e lbub hm
-      in foldr f HM.empty $ concatMap specificTS userBounds
-
-vectorizeDvs :: [[a]] -> [[a]] -> [a] -> [a]
-vectorizeDvs allStates allActions params = concat allStates ++ concat allActions ++ params
-
-msCoctave ::
-  State (Step Double) b
-  -> Integrator Double
-  -> Int
-  -> String
-  -> FilePath
-  -> IO ()
-msCoctave userStep' odeError n funDir name = do
-  let step = State.execState userStep' $
-             Step { stepStates  = Nothing
-                  , stepActions = Nothing
-                  , stepDxdt = Nothing
-                  , stepDt = Nothing
-                  , stepLagrangeTerm = Nothing
-                  , stepMayerTerm = Nothing
-                  , stepBounds = []
-                  , stepConstraints = []
-                  , stepParams = HS.empty
-                  , stepConstants = HS.empty
-                  , stepOutputs = HM.empty
-                  , stepPeriodic = HS.empty
-                  }
-      getWithErr :: String -> (Step Double -> Maybe c) -> c
-      getWithErr fieldName f = case f step of
-        Nothing -> error $ "need to set " ++ fieldName
-        Just ret -> ret
-
-      actions = getWithErr "actions" stepActions
-      dt      = getWithErr "dt"      stepDt
-      (states,outputs,dxdt,lagrangeState) = let
-        states'  = getWithErr "states" stepStates
-        dxdt'    = getWithErr "dxdt"   stepDxdt
-        outputs' = stepOutputs step
-        in
-         case stepLagrangeTerm step of
-           Nothing -> (states',outputs',dxdt',Nothing)
-           Just (lagrangeTerm,(lb,ub)) ->
-             ( states' ++ [lagrangeState']
-             , HM.union outputs' $ HM.fromList
-               [(lagrangeStateName, lagrangeState'), (lagrangeTermName, lagrangeTerm)]
-             , dxdt'++[lagrangeTerm]
-             , Just (lagrangeState',(lb,ub)) )
-              where
-                lagrangeState' = sym lagrangeStateName
-        
-      params    = HS.toList (stepParams    step)
-      constants = HS.toList (stepConstants step)
-
-      allStates   = [[sym $ show x ++ "__" ++ show k | x <-  states] | k <- [0..(n-1)]]
-      allActions  = [[sym $ show u ++ "__" ++ show k | u <- actions] | k <- [0..(n-1)]]
-      dvs = vectorizeDvs allStates allActions params
-
-      outputMap :: HashMap String [Expr Double]
-      outputMap = HM.map f outputs
-        where
-          f output = zipWith (subStatesActions output) allStates allActions
-
-      subStatesActions f x u = substitute f (zip states x ++ zip actions u)
-
-      subAllTimesteps :: Expr Double -> [Expr Double]
-      subAllTimesteps something = zipWith (subStatesActions something) allStates allActions
-
-      (lbs,ubs) = unzip $ vectorizeDvs stateBounds actionBounds paramBounds
-        where
-          (getAllBounds,getParamBounds) = setupBounds bounds n
-          stateBounds  = [[getAllBounds x k | x <- states ] | k <- [0..(n-1)]]
-          actionBounds = [[getAllBounds u k | u <- actions] | k <- [0..(n-1)]]
-          paramBounds  = [getParamBounds p | p <- params]
-
-          bounds = stepBounds step ++ lagrangeBound
-            where
-              lagrangeBound = case lagrangeState of
-                Nothing -> []
-                Just (ls,(lb,ub)) -> [(ls,(0,0,TIMESTEP 0)),(ls, (lb, ub, ALWAYS))]
-
-      cost = subStatesActions finalCost (last allStates) (last allActions)
-        where
-          finalCost = case (stepMayerTerm step, lagrangeState) of
-            (Just mc, Nothing) -> mc
-            (Nothing, Just (ls,_)) -> ls
-            (Just mc, Just (ls,_)) -> mc + ls
-            (Nothing,Nothing) -> error "need to set cost function"
-
-      (ceq, cineq) = foldl f ([],[]) allConstraints
-        where
-          f (eqs,ineqs) (Constraint x EQ y) = (eqs ++ [x - y], ineqs)
-          f (eqs,ineqs) (Constraint x LT y) = (eqs, ineqs ++ [x - y])
-          f (eqs,ineqs) (Constraint x GT y) = (eqs, ineqs ++ [y - x])
-      
-          execDxdt x u = map (flip substitute (zip states x ++ zip actions u)) dxdt
-
-          dodeConstraints = map (Constraint 0 EQ) $ concat $
-                            zipWith6 odeError (init allStates) (init allActions) (tail allStates) (tail allActions)
-                            (repeat execDxdt) (repeat dt)
-
-          allConstraints = dodeConstraints ++ (concatMap (g . (fmap subAllTimesteps)) (stepConstraints step)) ++ periodicConstraints
-            where
-              g (Constraint [] _ _) = []
-              g (Constraint _ _ []) = []
-              g (Constraint (x:xs) ord (y:ys)) = Constraint x ord y : g (Constraint xs ord ys)
-            
-              periodicConstraints = map lookup' $ HS.toList (stepPeriodic step)
-                where
-                  lookup' x = fromMaybe (error $ "couldn't find periodic thing \"" ++ show x ++ "\" in hashmap")
-                              $ HM.lookup x xuMap
-                  xuMap = HM.fromList $ zip states  (zipWith setEqual (head  allStates) (last allStates )) ++
-                                        zip actions (zipWith setEqual (head allActions) (last allActions))
-                    where
-                      setEqual x y = Constraint x EQ y
-
-  (costSource,costFg0,costFg) <- do
-    let costGrad = rad cost dvs
-    fg0 <- toFunGraph (dvs :* constants) (cost :* costGrad)
-    let fg = cse fg0
-    return (showMex (name ++ "_cost") fg, fg0, fg)
-  
-  (constraintsSource,constraintsFg0,constraintsFg) <- do
-    let cineqJacob = map (flip rad dvs) cineq
-        ceqJacob   = map (flip rad dvs) ceq
-    fg0 <- toFunGraph (dvs :* constants) (cineq :* ceq :* cineqJacob :* ceqJacob)
-    let fg = cse fg0
-    return (showMex (name ++ "_constraints") fg, fg0, fg)
-
-  (timeSource,timeFg) <- do
-    fg <- toFunGraph (dvs :* constants) (take n $ scanl (+) 0 (repeat dt))
-    return (showMex (name ++ "_time") fg, fg)
-
-  (outputSource,outputFg) <- do
-    fg <- toFunGraph (dvs :* constants) (HM.elems outputMap)
-    return (showMex (name ++ "_outputs") fg, fg)
-
-  (simSource,simFg) <- do
-    fg <- toFunGraph (states :* actions :* params :* constants) dxdt
-    return (showMex (name ++ "_sim") fg, fg)
-      
-  let setupSource = writeSetupSource name dvs lbs ubs
-      mexAllSource = writeMexAll name
-      unstructConstsSource = writeUnstructConsts name constants
-      structSource = writeToStruct name dvs params constants outputMap
-      unstructSource = writeUnstruct name dvs params states allStates actions allActions
-      plotSource = writePlot name outputMap
-
-  _ <- writeSourceFile         mexAllSource funDir $ name ++ "_mex_all.m"
-  _ <- writeSourceFile          setupSource funDir $ name ++ "_setup.m"
-  _ <- writeSourceFile         structSource funDir $ name ++ "_struct.m"
-  _ <- writeSourceFile unstructConstsSource funDir $ name ++ "_unstructConstants.m"
-  _ <- writeSourceFile       unstructSource funDir $ name ++ "_unstruct.m"
-  _ <- writeSourceFile           plotSource funDir $ name ++ "_plot.m"
-
-  _ <- writeSourceFile           timeSource funDir $ name ++ "_time.c"
-  _ <- writeSourceFile         outputSource funDir $ name ++ "_outputs.c"
-  _ <- writeSourceFile            simSource funDir $ name ++ "_sim.c"
-  _ <- writeSourceFile           costSource funDir $ name ++ "_cost.c"
-  _ <- writeSourceFile    constraintsSource funDir $ name ++ "_constraints.c"
-
-  putStrLn $ "nodes in time:        " ++ show (countNodes timeFg)
-  putStrLn $ "nodes in output:      " ++ show (countNodes outputFg)
-  putStrLn $ "nodes in sim:         " ++ show (countNodes simFg)
-  putStrLn $ "nodes in cost:        " ++ show (countNodes costFg) ++
-    " (" ++ show (countNodes costFg0) ++ " before CSE)"
-  putStrLn $ "nodes in constraints: " ++ show (countNodes constraintsFg) ++
-    " (" ++ show (countNodes constraintsFg0) ++ " before CSE)"
-  
-
-spring :: State (Step Double) ()
-spring = do
-  [x, v] <- setStates ["x","v"]
-  [u]    <- setActions ["u"]
-  [k, b] <- addConstants ["k", "b"]
-  let cost = 2*x*x + 3*v*v + 10*u*u
-  setDxdt [v, -k*x - b*v + u]
-  setDt (tEnd/((fromIntegral n')-1))
-
-  setLagrangeTerm cost (-1,2000)
-
-  setBound x (5,5) (TIMESTEP 0)
-  setBound v (0,0) (TIMESTEP 0)
-  
-  setBound x (-5,5) ALWAYS
-  setBound v (-10,10) ALWAYS
-  setBound u (-200, 200) ALWAYS
-
-  setBound v (0,0) (TIMESTEP (n'-1))
-
-  setPeriodic x
-  setPeriodic u
-
-tEnd :: Expr Double
-tEnd = 1.5
-
-n' :: Int
-n' = 18
-
-run :: IO ()
-run = msCoctave spring simpsonsRuleError' n' "../Documents/MATLAB/" "spring"
---run = msCoctave spring eulerError' n' "../Documents/MATLAB/" "spring"
diff --git a/Dvda/MultipleShooting/MSMonad.hs b/Dvda/MultipleShooting/MSMonad.hs
deleted file mode 100644
--- a/Dvda/MultipleShooting/MSMonad.hs
+++ /dev/null
@@ -1,155 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.MultipleShooting.MSMonad ( State
-                                     , setStates
-                                     , setActions
-                                     , addParam
-                                     , addParams
-                                     , addConstant
-                                     , addConstants
-                                     , setDxdt
-                                     , setLagrangeTerm
-                                     , setMayerTerm
-                                     , setDt
-                                     , addOutput
-                                     , setPeriodic
-                                     , addConstraint
-                                     , setBound
-                                     , lagrangeStateName
-                                     , lagrangeTermName
-                                     ) where
-
-import Data.Hashable ( Hashable )
-import qualified Data.HashSet as HS
-import Data.List ( nub, sort )
-import Data.Maybe ( isJust, fromMaybe )
-import Data.Monoid ( mappend )
-import Control.Monad ( when, zipWithM_ )
-import Control.Monad.State ( State )
-import qualified Control.Monad.State as State
-
-import qualified Dvda.HashMap as HM
-
-import Dvda.Expr ( Expr(..), sym )
-import Dvda.MultipleShooting.Types
-
-lagrangeStateName,lagrangeTermName :: String
-lagrangeStateName = "lagrangeState"
-lagrangeTermName = "lagrangeTerm"
-
-failDuplicates :: [String] -> [String]
-failDuplicates names
-  | length names == length (nub names) = names
-  | otherwise = error $ "ERROR: saw duplicate names in: " ++ show (sort names)
-
-checkOctaveName :: String -> String
-checkOctaveName name
-  | any (`elem` badChars) name =
-    error $ "ERROR: saw illegal octave variable character in string: \"" ++ name ++
-    "\", illegal characters: " ++ badChars
-  | name == lagrangeStateName = error "don't call your variable \"" ++ lagrangeStateName ++ "\", it's reserved"
-  | name == lagrangeTermName = error "don't call your variable \"" ++ lagrangeTermName ++ "\", it's reserved"
-  | otherwise = name
-  where
-    badChars = "\"'~!@#$%^&*()+`-=[]{}\\|;:,.<>/?"
-
-setStates :: [String] -> State (Step a) [Expr a]
-setStates names' = do
-  step <- State.get
-  case stepStates step of Just _ -> error "states already set, don't call setStates twice"
-                          Nothing -> do
-                            let names = failDuplicates (map checkOctaveName names')
-                                syms = map sym (failDuplicates names)
-                            State.put $ step {stepStates = Just syms}
-                            zipWithM_ addOutput syms names
-                            return syms
-
-setActions :: [String] -> State (Step a) [Expr a]
-setActions names' = do
-  step <- State.get
-  case stepActions step of Just _ -> error "actions already set, don't call setActions twice"
-                           Nothing -> do
-                             let names = failDuplicates (map checkOctaveName names')
-                                 syms = map sym (failDuplicates names)
-                             State.put $ step {stepActions = Just syms}
-                             zipWithM_ addOutput syms names
-                             return syms
-
-addParam :: (Eq a, Hashable a) => String -> State (Step a) (Expr a)
-addParam name = do
-  [blah] <- addParams [name]
-  return blah
-
-addConstant :: (Eq a, Hashable a) => String -> State (Step a) (Expr a)
-addConstant name = do
-  [blah] <- addConstants [name]
-  return blah
-
-addParams :: (Eq a, Hashable a) => [String] -> State (Step a) [Expr a]
-addParams names = do
-  step  <- State.get
-  let syms = map (sym . checkOctaveName) names
-      params0 = stepParams step
-  State.put $ step {stepParams = HS.union params0 (HS.fromList syms)}
-  return syms
-
-addConstants :: (Eq a, Hashable a) => [String] -> State (Step a) [Expr a]
-addConstants names = do
-  step  <- State.get
-  let syms = map (sym . checkOctaveName) names
-      constants0 = stepConstants step
-  State.put $ step {stepConstants = HS.union constants0 (HS.fromList syms)}
-  return syms
-
-addOutput :: Expr a -> String -> State (Step a) ()
-addOutput var name = do
-  step <- State.get
-  let hm = stepOutputs step
-      err = error $ "ERROR: already have an output with name: \"" ++ name ++ "\""
-  State.put $ step {stepOutputs = HM.insertWith err (checkOctaveName name) var hm}
-
-setDt :: Expr a -> State (Step a) ()
-setDt expr = do
-  step  <- State.get
-  when (isJust (stepDt step)) $ error "dt already set, don't call setDt twice"
-  State.put $ step {stepDt = Just expr}
-
-setPeriodic :: (Eq a, Hashable a, Show a) => Expr a -> State (Step a) ()
-setPeriodic var = do
-  step <- State.get
-  let newPeriodic
-        | var `HS.member` (stepPeriodic step) = error $ "you called setPeriodic twice on \"" ++ show var ++ "\""
-        | not (var `elem` (fromMaybe [] (mappend (stepStates step) (stepActions step)))) =
-          error $ "you can only make states or actions periodic, you can't make \"" ++ show var ++ "\" periodic"
-        | otherwise = HS.insert var (stepPeriodic step)
-  State.put $ step {stepPeriodic = newPeriodic}
--------------------------------------------
-
-setDxdt :: [Expr a] -> State (Step a) ()
-setDxdt vars = do
-  step  <- State.get
-  when (isJust (stepDxdt step)) $ error "dxdt already set, don't call setDxdt twice"
-  State.put $ step {stepDxdt = Just vars}
-
-setLagrangeTerm :: Expr a -> (a,a) -> State (Step a) ()
-setLagrangeTerm var (lb,ub) = do
-  step  <- State.get
-  when (isJust (stepLagrangeTerm step)) $ error "Lagrange term already set, don't call setLagrangeTerm twice"
-  State.put $ step {stepLagrangeTerm = Just (var,(lb,ub))}
-
-setMayerTerm :: Expr a -> State (Step a) ()
-setMayerTerm var = do
-  step  <- State.get
-  when (isJust (stepMayerTerm step)) $ error "Mayer term already set, don't call setMayerTerm twice"
-  State.put $ step {stepMayerTerm = Just var}
-
-setBound :: (Show a, Eq a, Hashable a)
-            => Expr a -> (a, a) -> BCTime -> State (Step a) ()
-setBound var@(ESym _) (lb, ub) bctime = do
-  step <- State.get
-  State.put $ step {stepBounds = (var, (lb,ub,bctime)):(stepBounds step)}
-setBound _ _ _ = error "WARNING - setBound called on non-design variable, try addConstraint instead"
-
-addConstraint :: Expr a -> Ordering -> Expr a -> State (Step a) ()
-addConstraint x ordering y =
-  State.state (\step -> ((), step {stepConstraints = (stepConstraints step) ++ [Constraint x ordering y]}))
diff --git a/Dvda/MultipleShooting/Types.hs b/Dvda/MultipleShooting/Types.hs
deleted file mode 100644
--- a/Dvda/MultipleShooting/Types.hs
+++ /dev/null
@@ -1,90 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language FlexibleContexts #-}
-
-module Dvda.MultipleShooting.Types ( Step(..)
-                                   , Constraint(..)
-                                   , Ode(..)
-                                   , BCTime(..)
-                                   , eulerError
-                                   , simpsonsRuleError
-                                   , eulerError'
-                                   , simpsonsRuleError'
-                                   ) where
-
-import Data.HashSet ( HashSet )
-
-import Dvda.Expr ( Expr(..) )
-import Dvda.HashMap ( HashMap )
-import Dvda.SparseLA
-
-data BCTime = ALWAYS | TIMESTEP Int deriving (Show, Eq)
-
-data Constraint a = Constraint a Ordering a deriving Show
-instance Functor Constraint where
-  fmap f (Constraint x ordering y) = Constraint (f x) ordering (f y)
-  
-
-data Step a = Step { stepStates :: Maybe [Expr a]
-                   , stepActions :: Maybe [Expr a]
-                   , stepParams :: HashSet (Expr a)
-                   , stepConstants :: HashSet (Expr a)
-                   , stepDxdt :: Maybe [Expr a]
-                   , stepLagrangeTerm :: Maybe (Expr a, (a,a))
-                   , stepMayerTerm :: Maybe (Expr a)
-                   , stepDt :: Maybe (Expr a)
-                   , stepBounds :: [(Expr a, (a,a, BCTime))]
-                   , stepConstraints :: [Constraint (Expr a)]
-                   , stepOutputs :: HashMap String (Expr a)
-                   , stepPeriodic :: HashSet (Expr a)
-                   }
-
-data Ode a = Ode (SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a)) (Int,Int)
-
-wrapOdeError :: Fractional (Expr a)
-                => (SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> Ode a -> Expr a -> SparseVec (Expr a))
-                -> [Expr a] -> [Expr a] -> [Expr a] -> [Expr a]
-                -> ([Expr a] -> [Expr a] -> [Expr a])
-                -> Expr a
-                -> [Expr a]
-wrapOdeError odeError xk uk xkp1 ukp1 dxdt dt =
-  denseListFromSv $ odeError xk' uk' xkp1' ukp1' (Ode dxdt' (error "FUUUUCK")) dt
-  where
-    xk'   = svFromList xk
-    xkp1' = svFromList xkp1
-    uk'   = svFromList uk
-    ukp1' = svFromList ukp1
-    dxdt' x u = svFromList $ dxdt (denseListFromSv x) (denseListFromSv u)
-
-eulerError' :: Fractional (Expr a)
-               => [Expr a] -> [Expr a] -> [Expr a] -> [Expr a]
-               -> ([Expr a] -> [Expr a] -> [Expr a])
-               -> Expr a
-               -> [Expr a]
-eulerError' = wrapOdeError eulerError
-
-simpsonsRuleError' :: Fractional (Expr a)
-                      => [Expr a] -> [Expr a] -> [Expr a] -> [Expr a]
-                      -> ([Expr a] -> [Expr a] -> [Expr a])
-                      -> Expr a
-                      -> [Expr a]
-simpsonsRuleError' = wrapOdeError simpsonsRuleError
-
-eulerError :: Fractional (Expr a) => SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> Ode a -> Expr a -> SparseVec (Expr a)
-eulerError xk uk xkp1 _ (Ode ode _) dt = xkp1 - (xk + svScale dt f0)
-  where
-    f0 = ode xk uk
-
-simpsonsRuleError :: Fractional (Expr a) => SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> Ode a -> Expr a -> SparseVec (Expr a)
-simpsonsRuleError xk uk xkp1 ukp1 (Ode ode _) dt = xkp1 - xk - (svScale (dt/6.0) (f0 + fourFm + f1))
-  where
-    f0 = ode xk uk
-    f1 = ode xkp1 ukp1
-
-    um = svScale 0.5 (uk + ukp1)
-    xm = xm' - xm''
-      where
-        xm' = svScale 0.5 (xk + xkp1)
-        xm'' = svScale (0.125 * dt) (f1 - f0)
-
-    fm = ode xm um
-    fourFm = svScale 4 fm
diff --git a/Dvda/Reify.hs b/Dvda/Reify.hs
deleted file mode 100644
--- a/Dvda/Reify.hs
+++ /dev/null
@@ -1,88 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language RankNTypes #-}
-{-# Language TemplateHaskell #-}
-{-# Language TypeFamilies #-}
-
--- this file is a modified version from Andy Gill's data-reify package
-
-module Dvda.Reify ( MuRef(..)
-                  , ReifyGraph(..)
-                  , reifyGraphs
-                  ) where
-
-import Control.Concurrent.MVar ( newMVar, takeMVar, putMVar, MVar, readMVar )
-import Control.Applicative ( Applicative )
-import Data.Hashable ( Hashable, hash )
-import Data.Traversable ( Traversable )
-import qualified Data.Traversable as T
-import System.Mem.StableName ( StableName, makeStableName, hashStableName )
-import Unsafe.Coerce ( unsafeCoerce )
-
-import Dvda.ReifyGraph ( ReifyGraph(..) )
-
-import qualified Data.HashTable.IO as H
-type HashTable k v = H.CuckooHashTable k v
-
-class MuRef a where
-  type DeRef a :: * -> *
-  mapDeRef :: Applicative f
-              => (forall b . (MuRef b, DeRef a ~ DeRef b) => b -> f u)
-              -> a
-              -> f (DeRef a u)
-
--- | 'reifyGraph' takes a data structure that admits 'MuRef', and returns a 'ReifyGraph' that contains
--- the dereferenced nodes, with their children as 'Int' rather than recursive values.
-reifyGraphs :: (MuRef s, Traversable t) => [t s] -> IO (ReifyGraph (DeRef s), [t Int])
-reifyGraphs m = do
-  stableNameMap <- H.new >>= newMVar
-  graph <- newMVar []
-  uVar <- newMVar 0
-  roots <- mapM (T.mapM (findNodes stableNameMap graph uVar)) m
-  pairs <- readMVar graph
-  return (ReifyGraph pairs, roots)
-
-findNodes :: MuRef s
-          => MVar (HashTable DynStableName Int)
-          -> MVar [(Int,DeRef s Int)]
-          -> MVar Int
-          -> s
-          -> IO Int
-findNodes stableNameMap graph uVar j | j `seq` True = do
-  st <- makeDynStableName j
-  tab <- takeMVar stableNameMap
-  amIHere <- H.lookup tab st
-  case amIHere of
-    -- if the j's StableName is already in the table, return the element
-    Just var -> do putMVar stableNameMap tab
-                   return var
-    -- if j's StableName is not yet in the table, recursively call findNodes
-    Nothing -> do var <- newUnique uVar
-                  H.insert tab st var
-                  putMVar stableNameMap tab
-                  res <- mapDeRef (findNodes stableNameMap graph uVar) j
-                  tab' <- takeMVar graph
-                  putMVar graph $ (var,res) : tab'
-                  return var
-findNodes _ _ _ _ = error "findNodes: strictness seq function failed to return True"
-
-newUnique :: MVar Int -> IO Int
-newUnique var = do
-  v <- takeMVar var
-  let v' = succ v
-  putMVar var v'
-  return v'
-  
--- Stable names that not use phantom types.
--- As suggested by Ganesh Sittampalam.
-data DynStableName = DynStableName (StableName ())
-
-instance Hashable DynStableName where
-  hash (DynStableName sn) = hashStableName sn
-  
-instance Eq DynStableName where
-	(DynStableName sn1) == (DynStableName sn2) = sn1 == sn2
-
-makeDynStableName :: a -> IO DynStableName
-makeDynStableName a = do
-	st <- makeStableName a
-	return $ DynStableName (unsafeCoerce st)
diff --git a/Dvda/ReifyGraph.hs b/Dvda/ReifyGraph.hs
deleted file mode 100644
--- a/Dvda/ReifyGraph.hs
+++ /dev/null
@@ -1,16 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language FlexibleContexts #-}
-{-# Language UndecidableInstances #-}
-
--- this file is a modified version from Andy Gill's data-reify package
-
-module Dvda.ReifyGraph ( ReifyGraph(..)
-                       ) where
-
-data ReifyGraph e = ReifyGraph [(Unique,e Unique)]
-
-type Unique = Int
-
--- | If 'e' is s Functor, and 'e' is 'Show'-able, then we can 'Show' a 'Graph'.
-instance (Show (e Int)) => Show (ReifyGraph e) where
-  show (ReifyGraph netlist) = show [ (u,e) | (u,e) <- netlist]
diff --git a/Dvda/ShowExprTests.hs b/Dvda/ShowExprTests.hs
new file mode 100644
--- /dev/null
+++ b/Dvda/ShowExprTests.hs
@@ -0,0 +1,40 @@
+{-# OPTIONS_GHC -Wall #-}
+
+module Dvda.ShowExprTests ( runTests
+                          ) where
+
+import Data.Maybe ( mapMaybe )
+
+import Dvda.Expr
+
+someShows :: [(String, Expr Double)]
+someShows = [ ("x * y", x * y)
+            , ("x / y", x / y)
+            , ("(x * y) * z", x * y * z)
+            , ("(x * y) / z", x * y / z)
+            , ("x / (y * z)", x / (y * z))
+            , ("cos(x)", cos x)
+            , ("sin(cos(x))", sin (cos x))
+            , ("sin(x ** y)", sin (x ** y))
+            , ("sin(x + y)", sin (x + y))
+            , ("x ** sin(y)", x ** sin y)
+            , ("(x + y) * z", (x + y)*z)
+            , ("10 * x", 10*x)
+            ]
+  where
+    x = sym "x"
+    y = sym "y"
+    z = sym "z"
+
+testShows :: [(String, Expr Double)] -> IO ()
+testShows = putStrLn . unlines . map betterShow . mapMaybe testShow
+  where
+    betterShow (x,y) = x ++ " =/= " ++ y
+    testShow (str,expr)
+      | expr' == str = Nothing
+      | otherwise = Just (expr', str)
+      where
+        expr' = show expr
+
+runTests :: IO ()
+runTests = testShows someShows
diff --git a/Dvda/SparseLA.hs b/Dvda/SparseLA.hs
deleted file mode 100644
--- a/Dvda/SparseLA.hs
+++ /dev/null
@@ -1,245 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.SparseLA ( SparseVec
-                     , SparseMat
-                     , svFromList
-                     , smFromLists
-                     , svFromSparseList
-                     , smFromSparseList
-                     , denseListFromSv
-                     , sparseListFromSv
-                     , svZeros
-                     , smZeros
-                     , svSize
-                     , smSize
-                     , svMap
-                     , smMap
-                     , svBinary
-                     , smBinary
-                     , svAdd
-                     , svSub
-                     , svMul
-                     , smAdd
-                     , smSub
-                     , smMul
-                     , svScale
-                     , smScale
-                     , getRow
-                     , getCol
-                     , svCat
-                     , svCats
-                     , sVV
-                     , sMV
-                     ) where
-
-import Data.List ( foldl' )
-import Data.Maybe ( fromJust, fromMaybe ) --, isNothing )
---import qualified Data.Traversable as T
-import Data.IntMap ( IntMap )
-import qualified Data.IntMap as IM
-
--- map from row to (map from col to value)
-data SparseMat a = SparseMat (Int,Int) (IntMap (IntMap a))
-
-instance Show a => Show (SparseMat a) where
-  show (SparseMat rowsCols xs) = "SparseMat " ++ show vals ++ " " ++ show rowsCols
-    where
-      vals = concatMap f (IM.toList xs)
-      f (row,m) = map g (IM.toList m)
-        where
-          g (col, val) = ((row, col), val)
-    
-instance Num a => Num (SparseMat a) where
-  x + y = fromJust $ smAdd x y
-  x - y = fromJust $ smSub x y
-  x * y = fromJust $ smMul x y
-  abs = smMap abs
-  signum = smMap signum
-  fromInteger = error "fromInteger not declared for Num SparseMat"
-
--- puts zeroes where there aren't entries
-denseListFromSv :: Num a => SparseVec a -> [a]
-denseListFromSv v@(SparseVec _ im) = IM.elems $ IM.union im (IM.fromList $ zip [0..n-1] (repeat 0))
-  where
-    n = svSize v
-
-sparseListFromSv :: SparseVec a -> [a]
-sparseListFromSv (SparseVec _ im) = IM.elems im
-  
-svZeros :: Int -> SparseVec a
-svZeros n = SparseVec n IM.empty
-
-smZeros :: (Int, Int) -> SparseMat a
-smZeros rowsCols = SparseMat rowsCols IM.empty
-
-smSize :: SparseMat a -> (Int,Int)
-smSize (SparseMat rowsCols _) = rowsCols
-
-smMap :: (a -> b) -> SparseMat a -> SparseMat b
-smMap f (SparseMat sh maps) = SparseMat sh (IM.map (IM.map f) maps)
-
-smFromLists :: [[a]] -> SparseMat a
-smFromLists blah = smFromSparseList sparseList (rows, cols)
-  where
-    rows = length blah
-    cols = length (head blah)
-    sparseList = concat $ zipWith (\row xs -> zipWith (\col x -> ((row,col),x)) [0..] xs) [0..] blah
-
-smFromSparseList :: [((Int,Int),a)] -> (Int,Int) -> SparseMat a
-smFromSparseList xs' rowsCols = SparseMat rowsCols (foldr f IM.empty xs')
-  where
-    f ((row,col), val) = IM.insertWith g row (IM.singleton col val)
-      where
-        g = IM.union
---        g = IM.unionWith (error $ "smFromList got 2 values for entry: "++show (row,col))
-
----- more efficient using mergeWithKey, but needs containers 0.5 so wait till ghc 7.6 :(
--- smBinary :: (a -> b -> c) -> (IntMap a -> IntMap c) -> (IntMap b -> IntMap c)
---             -> SparseMat a -> SparseMat b -> Maybe (SparseMat c)
--- smBinary fBoth fLeft fRight (SparseMat shx xs) (SparseMat shy ys)
---   | shx /= shy = Nothing
---   | isNothing merged = Nothing
---   | otherwise = Just $ SparseMat shx (fromJust merged)
---   where
---     merged = T.sequence $ IM.mergeWithKey f (IM.map (Just . fLeft)) (IM.map (Just . fRight)) xs ys
---       where
---         cols = Repa.shapeOfList [head $ Repa.listOfShape shx]
---         f _ x y = case svBinary fBoth fLeft fRight (SparseVec cols x) (SparseVec cols y) of
---           Just (SparseVec _ im) -> Just (Just im)
---           Nothing -> Just Nothing
-
-smBinary :: (a -> a -> a) -> (IntMap a -> IntMap a) -> (IntMap a -> IntMap a)
-            -> SparseMat a -> SparseMat a -> Maybe (SparseMat a)
-smBinary fBoth fLeft fRight (SparseMat shx@(_,cols) xs) (SparseMat shy ys)
-  | shx /= shy = Nothing
-  | otherwise = Just $ SparseMat shx merged
-  where
-    merged = IM.unionWith f (IM.map fLeft xs) (IM.map fRight ys)
-      where
-        f x y = case svBinary fBoth fLeft fRight (SparseVec cols x) (SparseVec cols y) of
-          Just (SparseVec _ im) -> im
-          Nothing -> error "goons everywhere"
-
---------------------------------------------------------------------------------------
-data SparseVec a = SparseVec Int (IntMap a)
-
-svSize :: SparseVec a -> Int
-svSize (SparseVec sh _) = sh
-
-instance Show a => Show (SparseVec a) where
-  show sv@(SparseVec _ xs) = "SparseVec " ++ show vals ++ " " ++ show rows
-    where
-      rows = svSize sv
-      vals = IM.toList xs
-
-instance Num a => Num (SparseVec a) where
-  x + y = fromJust $ svAdd x y
-  x - y = fromJust $ svSub x y
-  x * y = fromJust $ svMul x y
-  abs = svMap abs
-  signum = svMap signum
-  fromInteger = error "fromInteger not declared for Num SparseVec"
-
-svFromList :: [a] -> SparseVec a
-svFromList xs = svFromSparseList (zip [0..] xs) (length xs)
-
-svFromSparseList :: [(Int,a)] -> Int -> SparseVec a
-svFromSparseList xs rows = SparseVec rows (IM.fromList xs)
-
-svMap :: (a -> b) -> SparseVec a -> SparseVec b
-svMap f (SparseVec sh maps) = SparseVec sh (IM.map f maps)
-
-svBinary :: (a -> b -> c) -> (IntMap a -> IntMap c) -> (IntMap b -> IntMap c)
-            -> SparseVec a -> SparseVec b -> Maybe (SparseVec c)
-svBinary fBoth fLeft fRight (SparseVec shx xs) (SparseVec shy ys)
-  | shx /= shy = Nothing
-  | otherwise = Just $ SparseVec shx merged
-  where
-    -- more efficient using mergeWithKey, but needs containers 0.5 so wait till ghc 7.6 :(
---    merged = IM.mergeWithKey (\_ x y -> Just (fBoth x y)) fLeft fRight xs ys
-    merged = IM.unionWithKey f (fLeft xs) (fRight ys)
-      where
-        f k _ _ = fBoth (fromJust $ IM.lookup k xs) (fromJust $ IM.lookup k ys)
-
----------------------------------------------------------------------------
-svAdd :: Num a => SparseVec a -> SparseVec a -> Maybe (SparseVec a)
-svAdd = svBinary (+) id id
-
-svSub :: Num a => SparseVec a -> SparseVec a -> Maybe (SparseVec a)
-svSub = svBinary (-) id (IM.map negate)
-
-svMul :: Num a => SparseVec a -> SparseVec a -> Maybe (SparseVec a)
-svMul = svBinary (*) (\_ -> IM.empty) (\_ -> IM.empty)
-
-
-smAdd :: Num a => SparseMat a -> SparseMat a -> Maybe (SparseMat a)
-smAdd = smBinary (+) id id
-
-smSub :: Num a => SparseMat a -> SparseMat a -> Maybe (SparseMat a)
-smSub = smBinary (-) id (IM.map negate)
-
-smMul :: Num a => SparseMat a -> SparseMat a -> Maybe (SparseMat a)
-smMul = smBinary (*) (\_ -> IM.empty) (\_ -> IM.empty)
-
---------------------------------------------------------------------------
-
-svScale :: Num a => a -> SparseVec a -> SparseVec a
-svScale x (SparseVec sh xs) = SparseVec sh (IM.map (x *) xs)
-
-smScale :: Num a => a -> SparseMat a -> SparseMat a
-smScale x (SparseMat sh xs) = SparseMat sh (IM.map (IM.map (x *)) xs)
-
-
---------------------------------------------------------------------------
-getRow :: Int -> SparseMat a -> SparseVec a
-getRow row sm@(SparseMat (_,cols) xs)
-  | row >= (\(rows,_) -> rows) (smSize sm) =
-    error $ "getRow saw out of bounds index " ++ show row ++ " for matrix size " ++ show (smSize sm)
-  | otherwise = SparseVec cols out
-  where
-    out = fromMaybe IM.empty (IM.lookup row xs)
-
-getCol :: Int -> SparseMat a -> SparseVec a
-getCol col sm@(SparseMat (rows,_) xs)
-  | col >= (\(_,cols) -> cols) (smSize sm) =
-    error $ "getCol saw out of bounds index " ++ show col ++ " for matrix size " ++ show (smSize sm)
-  | otherwise = SparseVec rows out
-  where
-    out = IM.mapMaybe (IM.lookup col) xs
-
----------------------------------------------------------------------------
-sVV :: Num a => SparseVec a -> SparseVec a -> Maybe a
-sVV x y = fmap (\(SparseVec _ xs) -> sum (IM.elems xs)) (svMul x y)
-
-sMV :: Num a => SparseMat a -> SparseVec a -> Maybe (SparseVec a)
-sMV (SparseMat (mrows,mcols) ms) vec@(SparseVec vsize _)
-  | mcols /= vsize = Nothing
-  | otherwise = Just $ SparseVec mrows out
-  where
-    out = IM.mapMaybe f ms
-      where
-        f im = sVV (SparseVec mcols im) vec
-
----------------------------------------------------------------------------
-svCat :: SparseVec a -> SparseVec a -> SparseVec a
-svCat svx@(SparseVec _ xs) svy@(SparseVec _ ys) = SparseVec (shx + shy) (IM.union xs newYs)
-  where
-    shx = svSize svx
-    shy = svSize svy
-    newYs = IM.fromList $ map (\(k,x) -> (k+shx, x)) $ IM.toList ys
-
-svCats :: [SparseVec a] -> SparseVec a
-svCats [] = SparseVec 0 IM.empty
-svCats (xs0:xs) = foldl' svCat xs0 xs
-
---mx' :: SparseMat Double
---mx' = smFromList [((0,0), 10), ((0,2), 20), ((1,0), 30)] (2,3)
---
---my' :: SparseMat Double
---my' = smFromList [((0,0), 1), ((0,1), 7)] (2,3)
---
---x' :: SparseVec Int
---x' = svFromList [(0,10), (1, 20)] 4
---
---y' :: SparseVec Int
---y' = svFromList [(0,7), (3, 30)] 4
diff --git a/Dvda/Vis.hs b/Dvda/Vis.hs
deleted file mode 100644
--- a/Dvda/Vis.hs
+++ /dev/null
@@ -1,81 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
-module Dvda.Vis ( previewGraph
-                , previewGraph'
-                ) where
-
-import Control.Concurrent ( threadDelay )
-import Data.GraphViz ( Labellable, toLabelValue, preview )
-import Data.GraphViz.Attributes.Complete ( Label )
-import qualified Data.Graph.Inductive as FGL
-
-import Dvda.Expr
-import Dvda.FunGraph
-
--- | show a nice Dot graph
-previewGraph :: (Ord a, Show a) => FunGraph a -> IO ()
-previewGraph fg = do
-  preview $ toFGLGraph fg
-  threadDelay 10000
-
--- | show a nice Dot graph with labeled edges
-previewGraph' :: (Ord a, Show a) => FunGraph a -> IO ()
-previewGraph' fg = do
-  preview $ FGL.emap (\(FGLEdge x) -> FGLEdge' x) $ toFGLGraph fg
-  threadDelay 10000
-
-toFGLGraph :: FunGraph a -> FGL.Gr (FGLNode a) (FGLEdge a)
-toFGLGraph fg = FGL.mkGraph fglNodes fglEdges
-  where
-    fglNodes = map (\(k,gexpr) -> (k, FGLNode (k, gexpr))) $ fgReified fg
-    fglEdges = concatMap nodeToEdges $ fgReified fg
-      where
-        nodeToEdges (k,gexpr) = map (\p -> (p,k,FGLEdge (p,k,gexpr))) (getParents gexpr)
-
-data FGLNode a = FGLNode (Int, GExpr a Int)
-data FGLEdge a = FGLEdge (Int, Int, GExpr a Int)
-data FGLEdge' a = FGLEdge' (Int, Int, GExpr a Int)
-instance Eq a => Eq (FGLEdge a) where
-  (==) (FGLEdge (p0,k0,g0)) (FGLEdge (p1,k1,g1)) = (==) (p0,k0,g0) (p1,k1,g1)
-instance Eq a => Eq (FGLEdge' a) where
-  (==) (FGLEdge' (p0,k0,g0)) (FGLEdge' (p1,k1,g1)) = (==) (p0,k0,g0) (p1,k1,g1)
-instance Ord a => Ord (FGLEdge a) where
-  compare (FGLEdge (p0,k0,g0)) (FGLEdge (p1,k1,g1)) = compare (p0,k0,g0) (p1,k1,g1)
-instance Ord a => Ord (FGLEdge' a) where
-  compare (FGLEdge' (p0,k0,g0)) (FGLEdge' (p1,k1,g1)) = compare (p0,k0,g0) (p1,k1,g1)
-
-instance Labellable (FGLEdge a) where
-  toLabelValue (FGLEdge (p,k,_)) = toLabelValue $ show p ++ " --> " ++ show k
-instance Show a => Labellable (FGLEdge' a) where
-  toLabelValue (FGLEdge' (_,_,gexpr)) = toLabelValue $ show gexpr
-
-tlv :: Int -> String -> Label
-tlv k s = toLabelValue $ show k ++ ": " ++ s
-
-instance Show a => Labellable (FGLNode a) where
-  toLabelValue (FGLNode (k, (GSym s)))                       = tlv k (show s)
-  toLabelValue (FGLNode (k, (GConst c)))                     = tlv k (show c)
-  toLabelValue (FGLNode (k, (GNum (Mul _ _))))               = tlv k "*"
-  toLabelValue (FGLNode (k, (GNum (Add _ _))))               = tlv k "+"
-  toLabelValue (FGLNode (k, (GNum (Sub _ _))))               = tlv k "-"
-  toLabelValue (FGLNode (k, (GNum (Negate _))))              = tlv k "-"
-  toLabelValue (FGLNode (k, (GNum (Abs _))))                 = tlv k "abs"
-  toLabelValue (FGLNode (k, (GNum (Signum _))))              = tlv k "signum"
-  toLabelValue (FGLNode (k, (GNum (FromInteger x))))         = tlv k (show x)
-  toLabelValue (FGLNode (k, (GFractional (Div _ _))))        = tlv k "/"
-  toLabelValue (FGLNode (k, (GFractional (FromRational x)))) = tlv k (show (fromRational x :: Double))
-  toLabelValue (FGLNode (k, (GFloating (Pow _ _))))          = tlv k "**"
-  toLabelValue (FGLNode (k, (GFloating (LogBase _ _))))      = tlv k "logBase"
-  toLabelValue (FGLNode (k, (GFloating (Exp _))))            = tlv k "exp"
-  toLabelValue (FGLNode (k, (GFloating (Log _))))            = tlv k "log"
-  toLabelValue (FGLNode (k, (GFloating (Sin _))))            = tlv k "sin"
-  toLabelValue (FGLNode (k, (GFloating (Cos _))))            = tlv k "cos"
-  toLabelValue (FGLNode (k, (GFloating (ASin _))))           = tlv k "asin"
-  toLabelValue (FGLNode (k, (GFloating (ATan _))))           = tlv k "atan"
-  toLabelValue (FGLNode (k, (GFloating (ACos _))))           = tlv k "acos"
-  toLabelValue (FGLNode (k, (GFloating (Sinh _))))           = tlv k "sinh"
-  toLabelValue (FGLNode (k, (GFloating (Cosh _))))           = tlv k "cosh"
-  toLabelValue (FGLNode (k, (GFloating (Tanh _))))           = tlv k "tanh"
-  toLabelValue (FGLNode (k, (GFloating (ASinh _))))          = tlv k "asinh"
-  toLabelValue (FGLNode (k, (GFloating (ATanh _))))          = tlv k "atanh"
-  toLabelValue (FGLNode (k, (GFloating (ACosh _))))          = tlv k "acosh"
diff --git a/dvda.cabal b/dvda.cabal
--- a/dvda.cabal
+++ b/dvda.cabal
@@ -1,5 +1,5 @@
 Name:                dvda
-Version:             0.3.2.1
+Version:             0.4
 License:             BSD3
 License-file:        LICENSE
 Author:              Greg Horn
@@ -49,49 +49,47 @@
 --                     Dvda.OctaveSyntax
 --                     Dvda.Tests.Function
 --                     Dvda.Tests.Unary
-                     Dvda.SparseLA
 
                      Dvda.AD
-                     Dvda.CGen
+                     Dvda.Algorithm
+                     Dvda.Algorithm.Construct
+                     Dvda.Algorithm.Eval
+                     Dvda.Algorithm.FunGraph
+                     Dvda.Algorithm.Reify
  --                    Dvda.Codegen.CPlugins
-                     Dvda.Codegen.Gcc
-                     Dvda.Codegen.WriteFile
-                     Dvda.CSE
+--                     Dvda.Codegen.CGen
+--                     Dvda.Codegen.Gcc
+--                     Dvda.Codegen.PythonGen
+--                     Dvda.Codegen.WriteFile
+--                     Dvda.CSE
                      Dvda.Expr
-                     Dvda.Examples
-                     Dvda.FunGraph
-                     Dvda.MultipleShooting.CoctaveTemplates
-                     Dvda.MultipleShooting.MSCoctave
-                     Dvda.MultipleShooting.MSMonad
-                     Dvda.MultipleShooting.Types
-                     Dvda.Reify
-                     Dvda.ReifyGraph
-                     Dvda.Vis
-
-  Other-modules:     Dvda.HashMap
+--                     Dvda.Examples
+                     Dvda.HashMap
+--                     Dvda.Vis
 
+  Other-modules:     Dvda.ShowExprTests
   Build-depends:     base       >= 4     && < 5,
-                     file-location >= 0.4.4 && < 0.5,
-                     hashable  >= 1.1 && < 1.2,
-                     containers >= 0.4 && < 0.5,
-                     unordered-containers  >= 0.2 && < 0.3,
-                     hashtables  >= 1.0.1.6 && < 1.1,
-                     graphviz >= 2999.12 && < 2999.13,
-                     fgl >= 5.4 && < 5.5,
-                     mtl >= 2.0 && < 2.1,
-                     directory >= 1.1 && < 1.2,
-                     QuickCheck == 2.4.*,
-                     test-framework-quickcheck2,
-                     test-framework,
-                     process >= 1.1 && < 1.2
---                     text >= 0.11 && < 0.12,
---                     plugins >= 1.5 && < 1.6,
+                     hashable  >= 1.2,
+                     vector >= 0.10,
+                     unordered-containers  >= 0.2,
+                     containers >= 0.5,
+                     hashtables  >= 1.1.0,
+                     mtl
+--                     file-location >= 0.4.5 && < 0.5
+--                     graphviz >= 2999.15 && < 2999.17
+--                     fgl >= 5.4 && < 5.5
+--                     directory >= 1.2 && < 1.3
+--                     QuickCheck == 2.5.*
+--                     test-framework-quickcheck2
+--                     test-framework
+--                     process >= 1.1 && < 1.2
+--                     text >= 0.11 && < 0.12
+--                     plugins >= 1.5 && < 1.6
 --                     unix
---                     text,
+--                     text
 
   Ghc-options:       -Wall -O2
   GHC-Prof-Options:  -Wall -O2 -prof -fprof-auto -fprof-cafs -rtsopts
-  GHC-Shared-Options: -fPIC
 
 
 flag test
@@ -99,23 +97,22 @@
   default:     False
 
 Test-suite test
-  type:		     exitcode-stdio-1.0
+  type:              exitcode-stdio-1.0
   hs-source-dirs:    .
   main-is:           TestMain.hs
   build-depends:     base,
                      dvda,
-                     file-location >= 0.4.4 && < 0.5,
-                     hashable  >= 1.1 && < 1.2,
-                     hashtables  >= 1.0.1.6 && < 1.1,
-                     containers >= 0.4 && < 0.5,
-                     unordered-containers  >= 0.2 && < 0.3,
-                     graphviz >= 2999.12 && < 2999.13,
-                     fgl >= 5.4 && < 5.5,
-                     mtl >= 2.0 && < 2.1,
-                     directory >= 1.1 && < 1.2,
-                     QuickCheck == 2.4.*,
-                     process >= 1.1 && < 1.2,
---                     directory >= 1.1 && < 1.2,
+                     file-location,
+                     hashable,
+                     hashtables,
+                     containers,
+                     unordered-containers,
+                     graphviz,
+                     fgl,
+                     mtl,
+                     directory,
+                     QuickCheck,
+                     process,
                      ad,
                      test-framework-quickcheck2,
                      test-framework