diff --git a/Algebra/AD.hs b/Algebra/AD.hs
deleted file mode 100644
--- a/Algebra/AD.hs
+++ /dev/null
@@ -1,112 +0,0 @@
-{-# LANGUAGE RankNTypes, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, RebindableSyntax #-}
-module Algebra.AD (D(..),E(..),subst,dVar,var,sqrtE) where
-
-import Algebra.Classes hiding ((:+))
-import Data.Map (Map)
-import qualified Data.Map.Strict as M
-import Prelude hiding (Num(..),(/),fromRational,recip)
--- import Data.Vector as V
-import Data.Function (on)
-
-data AST v c = V v
-  | AST v c :* AST v c
-  | AST v c :+ AST v c
-  | AST v c :- AST v c
-  | K c
-
-instance (Show v, Show c) => Show (AST v c) where
-  showsPrec p (V v) = shows v
-  showsPrec p (K c ) = shows c
-  showsPrec p (x :+ y) = parens (p>2) (showsPrec 2 x . showString " + " . showsPrec 2 y)
-  showsPrec p (x :* y) = parens (p>3) (showsPrec 3 x . showString " + " . showsPrec 3 y)
-
-parens True x = showString "(" . x . showString ")"
-parens False x = x
-
-data D v c = D {dValue :: !c
-               ,dDerivs :: !(Map v c)
-               }
-  deriving Show
-
-dVar :: forall v c. Ring c => v -> c -> D v c
-dVar v c = D c (M.singleton v 1)
-
-var :: (Multiplicative c, Additive c, Ord v) => v -> E v c
-var v = E $ \env -> env v
-
-instance (Ord v,Additive c) => Additive (D v c) where
-  zero = D zero zero
-  D v1 d1 + D v2 d2 = D (v1 + v2) (d1 + d2)
-
-instance (Ord v,Group c) => Group (D v c) where
-  negate (D x d) = D (negate x) (negate d)
-  D v1 d1 - D v2 d2 = D (v1 - v2) (d1 - d2)
-
-instance Ord c => Ord (D v c) where
-  compare = compare `on` dValue
-
-instance Eq c => Eq (D v c) where
-  (==) = (==) `on` dValue
-
-instance (Ord v,Ring c) => Multiplicative (D v c) where
-  one = D one zero
-  D v1 d1 * D v2 d2 = D (v1 * v2) (v2 *^ d1 + v1 *^ d2)
-
-instance (AbelianAdditive c,Ord v) => AbelianAdditive (D v c)
-instance (Ord v,Ring c) => Module c (D v c) where
-  k *^ D v d = D (k * v) (k *^ d)
-
-instance (Ord v,Ring c) => Module (D v c) (D v c) where
-  (*^) = (*)
-
-instance (Ord v, Ring c) => Ring (D v c) where
-  fromInteger k = D (fromInteger k) zero
-
-newtype E v c = E {fromE :: (v -> D v c) -> D v c}
-
-instance (Ord v,Additive c) => Additive (E v c) where
-  zero = E (const zero)
-  (+) = liftE2 (+)
-
-instance (Ord v,Group c) => Group (E v c) where
-  negate (E x) = E (negate . x)
-  (-) = liftE2 (-)
-
-instance (Ord v,Ring c) => Multiplicative (E v c) where
-  one = E (const one)
-  (*) = liftE2 (*)
-
-instance (Ord v,Ring c) => AbelianAdditive (E v c)
-instance (Ord v,Ring c) => Module c (E v c) where
-  k *^ E x = E ((k *^) . x)
-
-instance (Ord v,Ring c) => Module (E v c) (E v c) where
-  (*^) = (*)
-
-instance (Ord v, Ring c) => Ring (E v c) where
-  fromInteger k = E (\ _ -> fromInteger k)
-
-liftE2 :: forall t t1. (D t t1 -> D t t1 -> D t t1) -> E t t1 -> E t t1 -> E t t1
-liftE2 f (E x) (E y) = E (\e -> f (x e) (y e))
-
-liftE :: forall t t1. (D t t1 -> D t t1) -> E t t1 -> E t t1
-liftE f (E x) = E (\e -> f (x e))
-
-subst :: E v c -> (v -> E v c) -> E v c
-subst (E p) f = E $ \k -> p (\a -> fromE (f a) k)
-
-sqrtD :: (Ord v, Floating c, Field c) => D v c -> D v c
-sqrtD (D v d) = D (sqrtv) ((0.5/sqrtv) *^ d)
-  where sqrtv = sqrt v
-
-sqrtE :: forall t t1. (Floating t1, Ord t, Field t1) => E t t1 -> E t t1
-sqrtE = liftE sqrtD
-
-instance (Field c,Ord v) => Division (D v c) where
-  recip (D v d) = D (recip v) (negate (square iv) *^ d)
-    where square x = x*x
-          iv = recip v
-
-instance (Field c,Ord v) => Division (E v c) where
-  recip = liftE recip
-  (/) = liftE2 (/)
diff --git a/Algebra/Linear/GaussianElimination.hs b/Algebra/Linear/GaussianElimination.hs
deleted file mode 100644
--- a/Algebra/Linear/GaussianElimination.hs
+++ /dev/null
@@ -1,32 +0,0 @@
-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, RebindableSyntax, DeriveFunctor #-}
-module Algebra.Linear.GaussianElimination (linSolve) where
-
-import Prelude hiding (Num(..),(/))
-import Algebra.Classes
-import qualified Data.Map.Strict as M
-import Data.Maybe
-import Algebra.Linear
-
-gaussianElim :: (Ord v,Eq c,Field c) => [Constraint v c] -> [(v,LinFunc v c)]
-gaussianElim [] = []
-gaussianElim (co@(Func f k):cos) = case M.minViewWithKey f of
-  Nothing -> if k == 0 then gaussianElim cos
-             else error "gaussianElim: unsatisfiable"
-  Just ((_,0),_) -> error "gaussianElim: invariant not respected"
-  Just ((v,c),_) -> (v,clean $ var v - co') : gaussianElim (fmap t cos)
-    where co' = (1/c) *^ co
-          t co2@(Func m _) = case (M.lookup v m) of
-            Nothing -> co2
-            Just cv -> clean (co2 - cv *^ co')
-
-substitution :: (Ord v,Ring c, Eq c) => [(v,LinFunc v c)] -> [(v,LinFunc v c)]
-substitution [] = []
-substitution ((v,f):vfs) = (v,f) : substitution [(v',subst v f f') | (v',f') <- vfs]
-
-subst :: (Eq scalar, Ord v, Ring scalar) => v -> LinFunc v scalar -> LinFunc v scalar -> LinFunc v scalar
-subst v f f'@(Func m _) = clean $ f' - coef *^ var v + coef *^ f
-  where coef = fromMaybe 0 (M.lookup v m)
-
-
-linSolve :: (Eq c, Ord v, Field c) => [Constraint v c] -> [(v, LinFunc v c)]
-linSolve = substitution . reverse . gaussianElim . fmap clean
diff --git a/Graphics/Diagrams/Core.hs b/Graphics/Diagrams/Core.hs
--- a/Graphics/Diagrams/Core.hs
+++ b/Graphics/Diagrams/Core.hs
@@ -1,13 +1,10 @@
-{-# LANGUAGE TypeSynonymInstances, FlexibleContexts, FlexibleInstances, GeneralizedNewtypeDeriving, MultiParamTypeClasses, RecursiveDo, TypeFamilies, OverloadedStrings, RecordWildCards,UndecidableInstances, PackageImports, TemplateHaskell, RankNTypes, GADTs, ImpredicativeTypes, DeriveFunctor, ScopedTypeVariables, ConstraintKinds #-}
+{-# LANGUAGE TypeSynonymInstances, FlexibleContexts, FlexibleInstances, GeneralizedNewtypeDeriving, MultiParamTypeClasses, RecursiveDo, TypeFamilies, OverloadedStrings, RecordWildCards,UndecidableInstances, PackageImports, TemplateHaskell, RankNTypes, GADTs, ImpredicativeTypes, DeriveFunctor, ScopedTypeVariables, ConstraintKinds, OverloadedStrings #-}
 
 module Graphics.Diagrams.Core (
   module Graphics.Diagrams.Types,
-  Expr, constant, newVars,
+  Expr, constant, absE, newVars,
   minimize, maximize,
   (===), (>==), (<==), (=~=),
-  GExpr,
-  minimize', maximize',
-  fromLinear,
   Diagram(..), runDiagram,
   drawText,
   freeze, relax, tighten,
@@ -15,57 +12,106 @@
   ) where
 
 import Prelude hiding (sum,mapM_,mapM,concatMap,Num(..),(/),fromRational,recip,(/))
+import qualified Prelude
 import Control.Monad.RWS hiding (forM,forM_,mapM_,mapM)
 import Algebra.Classes as AC
-import Algebra.Linear as Linear
-import Algebra.Linear.GaussianElimination
 import Data.Map (Map)
 import qualified Data.Map.Strict as M
 import Control.Lens hiding (element)
 import Data.Traversable
 import Data.Foldable
 import System.IO.Unsafe
-import Data.Vector (Vector,(!))
-import qualified Data.Vector as V
-import qualified Data.Vector.Storable as S
-import Numeric.Optimization.Algorithms.HagerZhang05 hiding (optimize)
-import qualified Numeric.Optimization.Algorithms.HagerZhang05 as HagerZhang05
-import Algebra.AD as AD
 import Graphics.Diagrams.Types
-
+import Data.List (isPrefixOf,intercalate)
+import Data.String
+import System.Process
+import SMT.Model
 
 -- | Expressions are linear functions of the variables
-type Expr = LinFunc Var Constant
 
 newtype Var = Var Int
   deriving (Ord,Eq,Show,Enum)
 
--- | A non-linear expression.
-type GExpr = E Var Constant
+class IsDouble a where
+  fromDouble :: Double -> a
+  -- sqrtE :: a -> a
+  absE :: a -> a
 
-fromLinear :: Expr -> GExpr
-fromLinear m = fromLinear' one m AD.var
+instance IsDouble Constant where
+  fromDouble = id
+  -- sqrtE = sqrt
+  absE = Prelude.abs
 
-fromLinear' :: forall a scalar. (Module scalar a) => a -> LinFunc Var scalar -> (Var -> a) -> a
-fromLinear' unit (Func m c) f = c *^ unit + fromSum (M.foldMapWithKey (\v k -> AC.Sum (k *^ f v)) m)
+-- | S-Expression represented as strings
+newtype SExpr = S {fromS :: String}
 
-substLinear :: Expr -> (Var -> Expr) -> Expr
-substLinear = fromLinear' (Func M.empty 1)
+-- Generic environment
+newtype R env y = R {_fromR :: env -> y}
+  deriving Functor
+instance Applicative (R env) where
+  pure x = R (\_ -> x)
+  R f <*> R x = R (\rho -> (f rho) (x rho))
+liftA2 :: forall (f :: * -> *) a b a1.
+            Applicative f =>
+            (a1 -> a -> b) -> f a1 -> f a -> f b
+liftA2 f x y = f <$> x <*> y
+instance IsDouble x => IsDouble (R env x) where
+  fromDouble x = pure (fromDouble x)
+  -- sqrtE = fmap sqrtE
+  absE = fmap absE
+instance Additive x => Additive (R env x) where
+  zero = pure zero
+  (+) = liftA2 (+)
+instance Multiplicative x => Multiplicative (R env x) where
+  one = pure one
+  (*) = liftA2 (*)
+instance Division x => Division (R env x) where
+  recip = fmap recip
+  (/) = liftA2 (/)
+instance AbelianAdditive x => AbelianAdditive (R env x)
+instance Ring x => Ring (R env x) where
+  fromInteger x = pure (fromInteger x)
+instance Field x => Field (R env x) where
+  fromRational x = pure (fromRational x)
+instance Group x => Group (R env x) where
+  negate = fmap negate
+  (-) = liftA2 (-)
 
-rename :: (Var -> Var) -> Expr -> Expr
-rename f e = substLinear e (Linear.var . f)
+newtype Expr = E {_fromE :: forall x. (Field x,IsDouble x) => R (Var -> x) x}
 
+instance Additive Expr where
+  zero = E zero
+  E x + E y = E (x+y)
+instance Group Expr where
+  E x - E y = E (x-y)
+-- Not a decidable theory: avoid.
+-- instance Multiplicative Expr where
+--   one = E one
+--   E x * E y = E (x*y)
+-- instance Division Expr where
+--   E x / E y = E (x/y)
+-- instance Ring Expr
+-- instance Field Expr
+-- instance Module Expr Expr where
+--   (*^) = (*)
+instance AbelianAdditive Expr
+instance Module Rational Expr where
+  k *^ E x = E (fromRational k*x)
+instance Module Constant Expr where
+  k *^ E x = E (fromDouble k*x)
 
 -- | Some action to perform after a solution has been found.
 data Freeze m where
   Freeze :: forall t m. Functor t => (t Constant -> m ()) -> t Expr -> Freeze m
 
+type Constraint = SExpr
+
 data DiagramState = DiagramState
-  {_diaNextVar :: Var
-  ,_diaLinConstraints :: [Constraint Var Constant]
-  ,_diaObjective :: GExpr
+  {_diaNextVar :: Var  -- ^ next var to allocate
+  ,_diaConstraints :: [Constraint]
+  ,_diaObjective :: Expr -- ^ objective function
   ,_diaVarNames :: Map Var String
-  ,_diaNoOverlaps :: [Pair (Point' GExpr)]
+  ,_diaNoOverlaps :: [Pair (Point' Expr)]
   }
 
 $(makeLenses ''DiagramState)
@@ -82,9 +128,37 @@
 -- Diagrams
 
 
+runDiagram :: Monad m => Backend lab m -> Diagram lab m a -> m a
+runDiagram backend diag = do
+  let env = Env one defaultPathOptions backend
+  (a,finalState,ds) <- runRWST (fromDia $ do x<-diag;resolveNonOverlaps;return x) env $
+    DiagramState (Var 0) [] zero (M.empty) []
+  let maxVar =  finalState ^. diaNextVar
+      decls = [sexp ["declare-const", smtVar x, "Real"] | x <- [Var 0 .. maxVar]]
+      constrs = intercalate "\n" $ map fromS $
+        decls ++
+        (unop "assert" <$> (finalState ^. diaConstraints)) ++
+        [ unop "minimize" (renderExpr (finalState ^. diaObjective)),
+          sexp ["check-sat"],
+          sexp ["get-model"]
+        ]
+      solution = unsafePerformIO $ do
+        writeFile "problem.smt2" $ constrs
+        _exitCode <- system "z3 -smt2 problem.smt2 > result.smt2"
+        res <- readFile "result.smt2"
+        let modelText = unlines . dropWhile (not . ("(model" `isPrefixOf`)) . lines $ res
+        case readModel modelText of
+          Right model -> return $ M.fromList model
+          Left err -> do print err
+                         error "die."
+      lkMod m (Var v) = M.findWithDefault (error "variable not in model") ("x"++show v) m
+
+  forM_ ds (\(Freeze f x) -> f (fmap (\(E (R g)) -> g (lkMod solution)) x))
+  return a
+
 -- | Relax the optimisation functions by the given factor
 relax :: Monad m => Rational -> Diagram lab m a -> Diagram lab m a
-relax factor = tighten (one/factor)
+relax factor = tighten (recip factor)
 
 -- | Tighten the optimisation functions by the given factor
 tighten :: Monad m => Rational -> Diagram lab m a -> Diagram lab m a
@@ -93,12 +167,6 @@
 --------------
 -- Variables
 
-
-newVar :: Monad m => String -> Diagram lab m Expr
-newVar name = do
-  [v] <- newVars [name]
-  return v
-
 newVars :: Monad m => [String] -> Diagram lab m [Expr]
 newVars kinds = forM kinds $ \name -> do
   v <- rawNewVar name
@@ -117,60 +185,17 @@
 -- Expressions
 
 
-runDiagram :: Monad m => Backend lab m -> Diagram lab m a -> m a
-runDiagram backend diag = do
-  let env = Env one defaultPathOptions backend
-  (a,finalState,ds) <- runRWST (fromDia $ do x<-diag;resolveNonOverlaps;return x) env $
-    DiagramState (Var 0) [] zero (M.empty) []
-  let reducedConstraints = M.fromList $ linSolve (finalState ^. diaLinConstraints)
-      linSolvSubst v = case M.lookup v reducedConstraints of
-        Nothing -> Linear.var v
-        Just x -> x
-      solvSubst = fromLinear . linSolvSubst
-      maxVar =  finalState ^. diaNextVar
-      freeVars = filter (\v -> M.notMember v reducedConstraints) [Var 0..maxVar]
-      varsToFreeVars' :: Map Var Int
-      varsToFreeVars' = M.fromList (zip freeVars [0..])
-      freeVarsToVars :: V.Vector Var
-      freeVarsToVars = V.fromList freeVars
-
-      obj' :: GExpr
-      obj' = AD.subst (finalState ^. diaObjective) solvSubst
-      grad' rho = (y,fmap (\v -> M.findWithDefault zero v dy) freeVarsToVars)
-         where D y dy = fromE obj' (\v -> maybe (error "not found") id (M.lookup v env') )
-               env' = M.mapWithKey (\v i -> D (rho ! i) (M.singleton v 1)) varsToFreeVars'
-      (solution,_result,_statistics) = unsafePerformIO $ do
-       putStrLn $ "free vars: " ++ show freeVars
-       putStrLn $ "reducedConstraints: " ++ show reducedConstraints
-       putStrLn $ "optimizing ..."
-       opt@(s,_,_) <- HagerZhang05.optimize
-        defaultParameters {verbose = VeryVerbose}
-        0.01
-        (V.replicate (length freeVars) 0)
-        (VFunction (fst . grad')) (VGradient (snd . grad')) (Just (VCombined grad'))
-       putStrLn $ "done: "  ++ show (grad' (S.convert s))
-       return opt
-      solution' = fmap (solution S.!) varsToFreeVars'
-      fullSolution = M.union solution' (fmap (valueIn solution') reducedConstraints)
-
-  forM_ ds (\(Freeze f x) -> f (fmap (valueIn fullSolution) x))
-  return a
-
-
-valueIn :: Map Var Double -> Expr -> Double
-valueIn sol f = fromLinear' one f (\v -> M.findWithDefault 0 v sol)
-
 -- | Embed a variable in an expression
 variable :: Var -> Expr
-variable v = Func (M.singleton v 1) 0
+variable v = E (R $ \rho -> rho v)
 
 -- | Embed a constant in an expression
-constant :: Constant -> Expr
-constant c = Func M.empty c
+constant :: Double -> Expr
+constant c = E (R $ \_ -> fromDouble c)
 
 satAll :: Monad m => String -> (Expr -> a -> Diagram lab m b) -> [a] -> Diagram lab m Expr
 satAll name p xs = do
-  [m] <- newVars [(name)]
+  [m] <- newVars [name]
   mapM_ (p m) xs
   return m
 
@@ -181,52 +206,25 @@
 
 --------------
 -- Expression constraints
-(===), (>==), (<==) :: Expr -> Expr -> Monad m => Diagram lab m ()
-e1 <== e2 = do
-  let Func f c = e1 - e2
-      isFalse = M.null f && c < 0
-  when isFalse $ error "Diagrams.Core: inconsistent constraint!"
-  minimize' $ E $ \s ->
-    let [v1,v2] = map (($ s) . fromE . fromLinear) [e1,e2]
-    in if v1 <= v2 then zero else square (square (v2-v1))
 
-prettyExpr :: Monad m => Expr -> Diagram lab m String
-prettyExpr (Func f k) = do
-  vnames <- Dia (use diaVarNames)
-  let vname n = case M.lookup n vnames of
-        Nothing -> error ("prettyExpr: variable not found: " ++ show n)
-        Just nm -> nm
-  return $ prettySum ([prettyProd c (vname v) | (v,c) <- M.assocs f]  ++ [show k | k /= 0])
-  where prettySum [] = "0"
-        prettySum xs = foldr1 prettyPlus xs
-        prettyPlus a ('-':b) = a ++ ('-':b)
-        prettyPlus x y = x ++ "+" ++ y
-        prettyProd 1 v = show v
-        prettyProd (-1) v = '-' : show v
-        prettyProd c v = show c ++ show v
-
+(===), (>==), (<==) :: Expr -> Expr -> Monad m => Diagram lab m ()
+e1 <== e2 = assert (e1 .<= e2)
 (>==) = flip (<==)
-
-e1 === e2 = do
-  constrName <- (\x y -> x ++ " = " ++ y) <$> prettyExpr e1 <*> prettyExpr e2
-  diaLinConstraints %= (e1 - e2 :)
+e1 === e2 = assert (e1 .== e2)
 
 -- | minimize the distance between expressions
-(=~=) :: Monad m => GExpr -> GExpr -> Diagram lab m ()
-x =~= y = minimize' $ square (x-y)
+(=~=) :: Monad m => Expr -> Expr -> Diagram lab m ()
+x =~= y = minimize $ absE (x-y)
 
 -------------------------
 -- Expression objectives
 
+
 minimize,maximize :: Monad m => Expr -> Diagram lab m ()
-minimize = minimize' . fromLinear
 maximize = minimize . negate
-
-minimize',maximize' :: Monad m => GExpr -> Diagram lab m ()
-maximize' = minimize' . negate
-minimize' f = do
+minimize f = do
   tightness <- view diaTightness
-  diaObjective %= \o -> (fromRational tightness::Double) *^ f + o
+  diaObjective %= \o -> tightness *^ f + o
 
 
 drawText :: Monad m => Point' Expr -> lab -> Diagram lab m BoxSpec
@@ -241,22 +239,72 @@
 -- Non-overlapping things
 
 registerNonOverlap :: Monad m => Point' Expr -> Point' Expr -> Diagram lab m ()
-registerNonOverlap nw se = Dia $ diaNoOverlaps %= (Pair (fromLinear <$> nw) (fromLinear <$>  se):)
+registerNonOverlap nw se = Dia $ diaNoOverlaps %= (Pair nw se:)
 
-surface :: forall a. Multiplicative a => Point' a -> a
-surface (Point x y) = x*y
+allPairs :: forall a. [a] -> [Pair a]
+allPairs [] = []
+allPairs (x:xs) = [Pair x y | y <- xs] ++ allPairs xs
 
 resolveNonOverlaps :: Monad m => Diagram lab m ()
 resolveNonOverlaps = do
   noOvl <- Dia $ use diaNoOverlaps
-  minimize' $ E $ \s ->
-    add $ do
-      pair <- allPairs noOvl
-      let (Pair bx1 bx2) = fmap (fmap (fmap (($ s) . fromE))) pair
-          overlap = inters bx1 bx2
-      return $ if nonEmpty overlap then (square $ surface overlap) else zero
-    where
-      allPairs [] = []
-      allPairs (x:xs) = [Pair x y | y <- xs] ++ allPairs xs
-      inters (Pair p1 q1) (Pair p2 q2) = (min <$> q1 <*> q2) - (max <$> p1 <*> p2)
-      nonEmpty (Point a b) = a > zero && b > zero
+  forM_ (allPairs noOvl) $ \p -> do
+    assert (disj (ffmap xpart p) .|| disj (ffmap ypart p))
+ where disj (Pair (Pair p1 q1) (Pair p2 q2)) = (q1 .<= p2) .|| (q2 .<= p1)
+       ffmap f = fmap (fmap f)
+
+
+---------------------------------
+-- Constraint & SExpr utils
+
+assert :: Monad m => SExpr -> Diagram lab m ()
+assert x = diaConstraints %= (x:)
+
+(.<=),(.==) :: Expr -> Expr -> Constraint
+x .<= y = binop "<=" (renderExpr x) (renderExpr y)
+
+x .== y = binop "=" (renderExpr x) (renderExpr y)
+
+(.||) :: Constraint -> Constraint -> Constraint
+(.||) = binop "or"
+
+renderExpr :: Expr -> SExpr
+renderExpr (E (R x)) = x smtVar
+
+smtVar :: Var -> SExpr
+smtVar (Var x) = S ("x" ++ show x)
+
+parens :: forall a. IsString [a] => [a] -> [a]
+parens x = "(" ++ x ++ ")"
+sexp :: [SExpr] -> SExpr
+sexp xs = S $ parens $ intercalate " " $ map fromS xs
+binop :: String -> SExpr -> SExpr -> SExpr
+binop s x y = sexp [S s,x,y]
+unop :: String -> SExpr -> SExpr
+unop s x = sexp [S s,x]
+instance Multiplicative (SExpr) where
+  (*) = binop "*"
+  one = S "1"
+instance Division (SExpr) where
+  (/) = binop "/"
+instance Additive (SExpr) where
+  (+) = binop "+"
+  zero = S "0"
+instance AbelianAdditive (SExpr)
+instance Field (SExpr)
+instance Ring (SExpr)
+instance Group (SExpr) where
+  negate = unop "-"
+  (-) = binop "-"
+instance IsString SExpr where
+  fromString = S
+
+instance IsDouble Expr where
+  fromDouble d = E (R $ \_ -> fromDouble d)
+  -- sqrtE (E x) = E (sqrtE x)
+  absE (E x) = E (absE x)
+
+instance IsDouble SExpr where
+  fromDouble x = S $ show x
+  -- sqrtE x = binop "^" x "0.5"
+  absE = unop "abs"
diff --git a/Graphics/Diagrams/DerivationTrees.hs b/Graphics/Diagrams/DerivationTrees.hs
--- a/Graphics/Diagrams/DerivationTrees.hs
+++ b/Graphics/Diagrams/DerivationTrees.hs
@@ -100,7 +100,7 @@
     let pt = ptObj # S
     pt `eastOf` (concl # W)
     pt `westOf` (concl # E)
-    fromLinear (xpart pt) =~= fromLinear (xpart (concl # Center))
+    (xpart pt) =~= (xpart (concl # Center))
     let top = ypart (concl # S)
     ypart pt + (fromIntegral steps *- layerHeight) === top
     using linkStyle $ path $ polyline [ptObj # Base,Point (xpart pt) top]
@@ -169,9 +169,9 @@
   -- layout hints (not necessary for "correctness")
   let xd = xdiff (separ # W) (psGrp # W)
   xd   === xdiff (psGrp # E) (separ # E)
-  relax 2 $ fromLinear (2 *- xd) =~= fromLinear premisesDist
+  relax 2 $ (2 *- xd) =~= premisesDist
   -- centering of conclusion
-  relax 3 $ minimize' $ sqNorm $ fmap fromLinear $ (separ # Center) - (concl # Center)
+  relax 3 $ minimize $ orthonorm $ (separ # Center) - (concl # Center)
 
   -- draw the rule.
   using ruleStyle $ path $ polyline [separ # W,separ # E]
diff --git a/Graphics/Diagrams/Object.hs b/Graphics/Diagrams/Object.hs
--- a/Graphics/Diagrams/Object.hs
+++ b/Graphics/Diagrams/Object.hs
@@ -8,7 +8,7 @@
 import Control.Monad
 import Control.Lens (set,view)
 import Algebra.Classes hiding (normalize)
-import Prelude hiding (Num(..))
+import Prelude hiding (Num(..),(/))
 
 data Anchor = Center | N | NW | W | SW | S | SE | E | NE | BaseW | Base | BaseE
   deriving Show
@@ -17,7 +17,9 @@
 type Box = Object
 
 type Anchorage = Anchor -> Point
-data Object = Object {objectName :: String, objectOutline :: Path, anchors :: Anchorage}
+data Object = Object { objectName :: !String
+                     , objectOutline :: !Path
+                     , anchors :: !Anchorage}
 
 
 infix 8 #
@@ -164,6 +166,7 @@
       p = circlePath (bx # Center) radius
   pathObject $ Object name p (anchors bx)
 
+-- | Debug, by tracing the bounding box of the object in a certain color.
 traceBox :: (Monad m) => Color -> Object -> Diagram lab m ()
 traceBox c l = do
   stroke c $ path $ polygon (map (l #) [NW,NE,SE,SW])
@@ -229,14 +232,15 @@
 -- vector.
 autoLabelObj :: Monad m => Box -> OVector -> Diagram lab m ()
 autoLabelObj lab (OVector pt v) = do
-  let normalVector :: Point' GExpr
-      normalVector = normalize $ fromLinear <$> v
+  let normalVector :: Point' Expr
+      normalVector = v
   -- label must touch the point
   tighten 10 $ pt `insideBox` lab
+  minimize (orthonorm (pt+v- lab#Center))
   -- go as far as possible in the normal direction
-  maximize' $ dotProd (fromLinear <$> ((lab#Center) - pt)) normalVector
-  -- don't stray from the normal line
-  minimize' $ square $ dotProd (fromLinear <$> ((lab#Center) - pt)) (rotate90 normalVector)
+  -- maximize $ dotProd (((lab#Center) - pt)) normalVector
+  -- don't stray away from the normal line
+  -- minimize $ absE $ dotProd (((lab#Center) - pt)) (rotate90 normalVector)
   --
 
 -- | @autoLabel o i@ Layouts the label object @o@ at the given incidence
diff --git a/Graphics/Diagrams/Path.hs b/Graphics/Diagrams/Path.hs
--- a/Graphics/Diagrams/Path.hs
+++ b/Graphics/Diagrams/Path.hs
@@ -111,8 +111,8 @@
    CurveTo (pt (negate r) (negate k)) (pt (negate k) (negate r)) (pt zero (negate r)),
    CurveTo (pt k (negate r)) (pt r (negate k)) (pt r zero),
    Cycle]
- where k1 :: Constant
-       k1 = 4 * (sqrt 2 - 1) / 3
+ where k1 :: Double
+       k1 = fromInteger 4 * (sqrt (fromInteger 2) - (fromInteger 1)) / fromInteger 3
        k = k1 *^ r
        pt x y = center + (Point x y)
 
diff --git a/Graphics/Diagrams/Point.hs b/Graphics/Diagrams/Point.hs
--- a/Graphics/Diagrams/Point.hs
+++ b/Graphics/Diagrams/Point.hs
@@ -7,7 +7,6 @@
 import Data.List (transpose)
 import Prelude hiding (sum,mapM_,mapM,concatMap,maximum,minimum,Num(..),(/))
 import Algebra.Classes
-import Algebra.AD (sqrtE)
 
 infix 4 .=.
 ----------------
@@ -17,12 +16,18 @@
 
 type Point = Point' Expr
 
+orthonorm :: Point -> Expr
+orthonorm (Point x y) = absE x + absE y
+
+{-
+
+
 -- | Norm of a vector. Don't minimize this: the solver does not like functions
 -- with non-continuous derivatives (at zero in this case).
-norm :: Point' GExpr -> GExpr
+norm :: Point' Expr -> Expr
 norm p = sqrtE (sqNorm p)
 
-normalize :: Point' GExpr -> Point' GExpr
+normalize :: Point' Expr -> Point' Expr
 normalize x = (one/norm x) *^ x
 
 -- | Dot product
@@ -32,14 +37,14 @@
 -- | Squared norm of a vector
 sqNorm :: forall a. (Ring a) => Point' a -> a
 sqNorm p = dotProd p p
-
+-}
 -- | Rotate a vector 90 degres in the trigonometric direction.
 rotate90 :: forall a. Group a => Point' a -> Point' a
 rotate90 (Point x y) = Point (negate y) x
 
 -- | Rotate a vector 180 degres
 rotate180 :: forall a. Group a => Point' a -> Point' a
-rotate180 x = rotate90 . rotate90 $ x
+rotate180 = rotate90 . rotate90
 
 xdiff,ydiff :: Point -> Point -> Expr
 xdiff p q = xpart (q - p)
diff --git a/Graphics/Diagrams/Types.hs b/Graphics/Diagrams/Types.hs
--- a/Graphics/Diagrams/Types.hs
+++ b/Graphics/Diagrams/Types.hs
@@ -24,7 +24,7 @@
 
 -- | Average
 avg :: Module Constant a => [a] -> a
-avg xs = (1/fromIntegral (length xs)) *- add xs
+avg xs = (one/fromIntegral (length xs)::Constant) *^ add xs
 
 ------------
 -- Types
diff --git a/SMT/Model.hs b/SMT/Model.hs
new file mode 100644
--- /dev/null
+++ b/SMT/Model.hs
@@ -0,0 +1,49 @@
+module SMT.Model (readModel) where
+
+import Text.ParserCombinators.Parsek.Position
+import Data.Char (isSpace)
+
+type P a = Parser a
+
+tok :: String -> P ()
+tok s = spaces >> string s >> return ()
+
+many1 p = (:) <$> p <*> many p
+
+parseDouble :: P Double
+parseDouble = do
+  spaces
+  x <- many1 digit
+  string "."
+  y <- many1 digit
+  return $ read $ x ++ "." ++ y
+parseValue = parseDouble <|> parens (parseDiv <|> parseNeg)
+
+parseDiv = do
+  tok "/"
+  x <- parseValue
+  y <- parseValue
+  return (x/y)
+
+parseNeg = do
+  tok "-"
+  x <- parseValue
+  return (negate x)
+
+parseAssoc :: P (String,Double)
+parseAssoc = parens $ do
+  tok "define-fun"
+  spaces
+  v <- many1 (satisfy (not . isSpace))
+  parens (return ())
+  tok "Real"
+  x <- parseValue
+  return (v,x)
+
+parens = between (tok "(") (tok")")
+parseModel = parens $ do
+  tok "model"
+  many parseAssoc
+
+readModel :: String -> ParseResult SourcePos [(String, Double)]
+readModel = parse "<model>" parseModel longestResult
diff --git a/Setup.hs b/Setup.hs
--- a/Setup.hs
+++ b/Setup.hs
@@ -1,2 +1,4 @@
+#!/usr/bin/env runhaskell
 import Distribution.Simple
+main :: IO ()
 main = defaultMain
diff --git a/lp-diagrams.cabal b/lp-diagrams.cabal
--- a/lp-diagrams.cabal
+++ b/lp-diagrams.cabal
@@ -1,5 +1,5 @@
 name:                lp-diagrams
-version:             2.0
+version:             2.0.0
 synopsis:            An EDSL for diagrams based based on linear constraints
 license:             AGPL-3
 license-file:        LICENSE
@@ -10,11 +10,20 @@
 build-type:          Simple
 -- extra-source-files:  
 cabal-version:       >=1.18
+description:
+  A library to describe diagrams. The defining
+  feature of the package is the ability to use linear constraints to
+  specify layout, which are resolved using z3 (latest z3 must be installed).
+  Backends are provided either of the following packages lp-diagrams-svg (svg) or marxup (tikz).
 
 Flag graphviz
   Description: Enable graphviz support
   Default:     True
 
+source-repository head
+  type: git
+  location: https://github.com/jyp/lp-diagrams
+
 library
   if flag(graphviz)
     build-depends: graphviz
@@ -28,21 +37,20 @@
                        Graphics.Diagrams.Types,
                        Graphics.Diagrams.DerivationTrees,
                        Graphics.Diagrams
+  other-modules:       SMT.Model
 
-                       Algebra.AD
                        Algebra.Linear
-                       Algebra.Linear.GaussianElimination
   build-depends:       base >=4.8 && < 666,
-                       ad,
                        lens >=4.12,
                        text >=1.2 ,
                        typography-geometry >=1.0 ,
                        gasp,
-                       nonlinear-optimization,
                        reflection,
                        vector,
                        polynomials-bernstein,
                        mtl >=2.2 ,
                        containers >=0.5,
-                       labeled-tree
+                       labeled-tree,
+                       parsek,
+                       process
   default-language:    Haskell2010