diff --git a/limp.cabal b/limp.cabal
--- a/limp.cabal
+++ b/limp.cabal
@@ -1,5 +1,5 @@
 name:                limp
-version:             0.3.2.0
+version:             0.3.2.1
 synopsis:            representation of Integer Linear Programs
 description:         so far, this package just provides two representations for linear programs: "Numeric.Limp.Program", which is what I expect end-users to use, and
                      "Numeric.Limp.Canon", which is simpler, but would be less nice for writing linear programs.
@@ -23,6 +23,9 @@
   hs-source-dirs: src
   exposed-modules:
         Numeric.Limp.Rep
+        Numeric.Limp.Rep.Rep
+        Numeric.Limp.Rep.Arbitrary
+        Numeric.Limp.Rep.IntDouble
         Numeric.Limp.Error
 
         Numeric.Limp.Program.Bounds
@@ -45,7 +48,10 @@
         Numeric.Limp.Canon.Simplify.Subst
         Numeric.Limp.Canon.Simplify
 
-  -- other-modules:       
+        Numeric.Limp.Solve.Simplex.StandardForm
+        Numeric.Limp.Solve.Simplex.Maps
+        Numeric.Limp.Solve.Branch.Simple
+
   build-depends:
         base        < 5,
         containers  == 0.5.*
@@ -59,6 +65,15 @@
   type: exitcode-stdio-1.0
   main-is: Main.hs
   hs-source-dirs: tests
+  other-modules:
+        Arbitrary.Assignment
+        Arbitrary.Program
+        Arbitrary.Var
+        BranchExample
+        Convert
+        SimplexExample
+        Simplexs
+        Simplify
   build-depends:
         base        < 5,
         containers  == 0.5.*,
diff --git a/src/Numeric/Limp/Rep.hs b/src/Numeric/Limp/Rep.hs
--- a/src/Numeric/Limp/Rep.hs
+++ b/src/Numeric/Limp/Rep.hs
@@ -5,85 +5,11 @@
 -- We bundle Z and R up into a single representation instead of abstracting over both,
 -- because we must be able to convert from Z to R without loss.
 --
-module Numeric.Limp.Rep where
-
-import Data.Map (Map)
-import qualified Data.Map as M
-
-import Data.Monoid
-
--- | The Representation class. Requires its members @Z c@ and @R c@ to be @Num@, @Ord@ and @Eq@.
---
--- For some reason, for type inference to work, the members must be @data@ instead of @type@.
--- This gives some minor annoyances when unpacking them. See 'unwrapR' below.
---
-class ( Num (Z c), Ord (Z c), Eq (Z c), Integral (Z c)
-      , Num (R c), Ord (R c), Eq (R c), RealFrac (R c)) => Rep c where
-
- -- | Integers
- data Z c
- -- | Real numbers
- data R c
-
- -- | Convert an integer to a real. This should not lose any precision.
- -- (whereas @fromIntegral 1000 :: Word8@ would lose precision)
- fromZ :: Z c -> R c
- fromZ = fromIntegral
-
-
--- | An assignment from variables to values.
--- Maps integer variables to integers, and real variables to reals.
-data Assignment z r c
- = Assignment (Map z (Z c)) (Map r (R c))
-
-deriving instance (Show (Z c), Show (R c), Show z, Show r) => Show (Assignment z r c)
-
-instance (Ord z, Ord r) => Monoid (Assignment z r c) where
- mempty = Assignment M.empty M.empty
- mappend (Assignment z1 r1) (Assignment z2 r2)
-  = Assignment (M.union z1 z2) (M.union r1 r2)
-
-
--- | Retrieve value of integer variable - or 0, if there is no value.
-zOf :: (Rep c, Ord z) => Assignment z r c -> z -> Z c
-zOf (Assignment zs _) z
- = maybe 0 id $ M.lookup z zs
-
--- | Retrieve value of real variable - or 0, if there is no value.
-rOf :: (Rep c, Ord r) => Assignment z r c -> r -> R c
-rOf (Assignment _ rs) r
- = maybe 0 id $ M.lookup r rs
-
--- | Retrieve value of an integer or real variable, with result cast to a real regardless.
-zrOf :: (Rep c, Ord z, Ord r) => Assignment z r c -> Either z r -> R c
-zrOf a = either (fromZ . zOf a) (rOf a)
-
-assSize :: Assignment z r c -> Int
-assSize (Assignment mz mr)
- = M.size mz + M.size mr
-
-
--- | A representation that uses native 64-bit ints and 64-bit doubles.
--- Really, this should be 32-bit ints.
-data IntDouble
-
-instance Rep IntDouble where
- -- | Automatically defer numeric operations to the native int.
- newtype Z IntDouble = Z Int
-    deriving (Ord,Eq,Integral,Real,Num,Enum)
- newtype R IntDouble = R Double
-    deriving (Ord,Eq,Num,Enum,Fractional,Real,RealFrac)
-
--- | Define show manually, so we can strip out the "Z" and "R" prefixes.
-instance Show (Z IntDouble) where
- show (Z i) = show i
-
-instance Show (R IntDouble) where
- show (R i) = show i
-
-
+module Numeric.Limp.Rep
+    ( module Numeric.Limp.Rep.Rep
+    , module Numeric.Limp.Rep.IntDouble )
+    where
 
--- | Convert a wrapped (R IntDouble) to an actual Double.
-unwrapR :: R IntDouble -> Double
-unwrapR (R d) = d
+import Numeric.Limp.Rep.Rep
+import Numeric.Limp.Rep.IntDouble
 
diff --git a/src/Numeric/Limp/Rep/Arbitrary.hs b/src/Numeric/Limp/Rep/Arbitrary.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Limp/Rep/Arbitrary.hs
@@ -0,0 +1,22 @@
+-- | Arbitrary precision number representation
+module Numeric.Limp.Rep.Arbitrary where
+import Numeric.Limp.Rep.Rep
+
+-- | A representation that uses arbitrary-sized Integers and Rationals
+data Arbitrary
+
+instance Rep Arbitrary where
+ -- | Automatically defer numeric operations to the native int.
+ newtype Z Arbitrary = Z Integer
+    deriving (Ord,Eq,Integral,Real,Num,Enum)
+ newtype R Arbitrary = R Rational
+    deriving (Ord,Eq,Num,Enum,Fractional,Real,RealFrac)
+
+-- | Define show manually, so we can strip out the "Z" and "R" prefixes.
+instance Show (Z Arbitrary) where
+ show (Z i) = show i
+
+instance Show (R Arbitrary) where
+ show (R i) = show i
+
+
diff --git a/src/Numeric/Limp/Rep/IntDouble.hs b/src/Numeric/Limp/Rep/IntDouble.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Limp/Rep/IntDouble.hs
@@ -0,0 +1,29 @@
+-- | Fixed/floating precision number representation
+module Numeric.Limp.Rep.IntDouble where
+import Numeric.Limp.Rep.Rep
+
+-- | A representation that uses native 64-bit ints and 64-bit doubles.
+-- Really, this should be 32-bit ints.
+data IntDouble
+
+instance Rep IntDouble where
+ -- | Automatically defer numeric operations to the native int.
+ newtype Z IntDouble = Z Int
+    deriving (Ord,Eq,Integral,Real,Num,Enum)
+ newtype R IntDouble = R Double
+    deriving (Ord,Eq,Num,Enum,Fractional,Real,RealFrac)
+
+-- | Define show manually, so we can strip out the "Z" and "R" prefixes.
+instance Show (Z IntDouble) where
+ show (Z i) = show i
+
+instance Show (R IntDouble) where
+ show (R i) = show i
+
+
+
+-- | Convert a wrapped (R IntDouble) to an actual Double.
+unwrapR :: R IntDouble -> Double
+unwrapR (R d) = d
+
+
diff --git a/src/Numeric/Limp/Rep/Rep.hs b/src/Numeric/Limp/Rep/Rep.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Limp/Rep/Rep.hs
@@ -0,0 +1,64 @@
+-- | Representation of integers (Z) and reals (R) of similar precision.
+-- Programs are abstracted over this, so that ideally in the future we could have a
+-- solver that produces Integers and Rationals, instead of just Ints and Doubles.
+--
+-- We bundle Z and R up into a single representation instead of abstracting over both,
+-- because we must be able to convert from Z to R without loss.
+--
+module Numeric.Limp.Rep.Rep where
+
+import Data.Map (Map)
+import qualified Data.Map as M
+
+import Data.Monoid
+
+-- | The Representation class. Requires its members @Z c@ and @R c@ to be @Num@, @Ord@ and @Eq@.
+--
+-- For some reason, for type inference to work, the members must be @data@ instead of @type@.
+-- This gives some minor annoyances when unpacking them. See 'unwrapR' below.
+--
+class ( Num (Z c), Ord (Z c), Eq (Z c), Integral (Z c)
+      , Num (R c), Ord (R c), Eq (R c), RealFrac (R c)) => Rep c where
+
+ -- | Integers
+ data Z c
+ -- | Real numbers
+ data R c
+
+ -- | Convert an integer to a real. This should not lose any precision.
+ -- (whereas @fromIntegral 1000 :: Word8@ would lose precision)
+ fromZ :: Z c -> R c
+ fromZ = fromIntegral
+
+
+-- | An assignment from variables to values.
+-- Maps integer variables to integers, and real variables to reals.
+data Assignment z r c
+ = Assignment (Map z (Z c)) (Map r (R c))
+
+deriving instance (Show (Z c), Show (R c), Show z, Show r) => Show (Assignment z r c)
+
+instance (Ord z, Ord r) => Monoid (Assignment z r c) where
+ mempty = Assignment M.empty M.empty
+ mappend (Assignment z1 r1) (Assignment z2 r2)
+  = Assignment (M.union z1 z2) (M.union r1 r2)
+
+
+-- | Retrieve value of integer variable - or 0, if there is no value.
+zOf :: (Rep c, Ord z) => Assignment z r c -> z -> Z c
+zOf (Assignment zs _) z
+ = maybe 0 id $ M.lookup z zs
+
+-- | Retrieve value of real variable - or 0, if there is no value.
+rOf :: (Rep c, Ord r) => Assignment z r c -> r -> R c
+rOf (Assignment _ rs) r
+ = maybe 0 id $ M.lookup r rs
+
+-- | Retrieve value of an integer or real variable, with result cast to a real regardless.
+zrOf :: (Rep c, Ord z, Ord r) => Assignment z r c -> Either z r -> R c
+zrOf a = either (fromZ . zOf a) (rOf a)
+
+assSize :: Assignment z r c -> Int
+assSize (Assignment mz mr)
+ = M.size mz + M.size mr
+
diff --git a/src/Numeric/Limp/Solve/Branch/Simple.hs b/src/Numeric/Limp/Solve/Branch/Simple.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Limp/Solve/Branch/Simple.hs
@@ -0,0 +1,84 @@
+-- | The simplest, stupidest possible branch and bound algorithm.
+--
+--
+module Numeric.Limp.Solve.Branch.Simple
+    (branch, makeIntegral)
+    where
+import Numeric.Limp.Canon.Program
+import Numeric.Limp.Canon.Simplify
+import Numeric.Limp.Rep
+
+import Control.Applicative
+import Control.Monad
+import qualified Data.Map as M
+import Data.Monoid
+
+branch
+    :: (Ord z, Ord r, Rep c)
+    => (Program z r c -> Maybe (Assignment () (Either z r) c, R c))
+    -> Program z r c
+    -> Maybe (Assignment z r c, R c)
+branch solver start_prog
+ = go mempty start_prog
+ where
+  go ass p
+   -- TODO:
+   -- simp can actually change the objective function
+   -- because Canon doesn't store a constant summand on the objective.
+   -- we really need to return the modified summand and take that into account when
+   -- choosing between two integer assignments.
+   | Right (ass', p') <- simplify' ass p
+   = do  (assRelax,co) <- solver p'
+         case makeIntegral assRelax of
+          Left (var, val)
+           -> branchon p' ass' (Left var) val
+          Right r
+           -> Just (ass' <> r, co)
+   | otherwise
+   = Nothing
+
+  branchon p ass var val
+   = let lo = addBound p var (Just (fromZ $ truncate val + 1), Nothing)
+         up = addBound p var (Nothing, Just (fromZ $ truncate val))
+         loB     = go ass lo
+         upB     = go ass up
+     in case (loB, upB) of
+        (Just (a1, o1), Just (a2, o2))
+         | o1 > o2
+         -> Just (a1, o1)
+         | otherwise
+         -> Just (a2, o2)
+        (Just r, Nothing)
+         -> Just r
+        (Nothing, Just r)
+         -> Just r
+        (Nothing, Nothing)
+         -> Nothing
+     
+
+  addBound p v b
+   = let bs = _bounds p
+         b' = maybe (Nothing,Nothing) id
+            $ M.lookup v bs
+     in  p { _bounds = M.insert v (mergeBounds b b') bs }
+
+makeIntegral
+    :: (Ord z, Ord r, Rep c)
+    => Assignment () (Either z r) c
+    -> Either (z, R c)
+              (Assignment z r c)
+makeIntegral (Assignment _ vs)
+ =   uncurry Assignment
+ <$> foldM go (M.empty, M.empty) (M.toList vs)
+ where
+  go (zs,rs) (var, val)
+   = case var of
+      Right r
+       -> return (zs, M.insert r val rs)
+      Left z
+       | val' <- truncate val
+       , val == fromZ val'
+       -> return (M.insert z val' zs, rs)
+       | otherwise
+       -> Left (z, val)
+
diff --git a/src/Numeric/Limp/Solve/Simplex/Maps.hs b/src/Numeric/Limp/Solve/Simplex/Maps.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Limp/Solve/Simplex/Maps.hs
@@ -0,0 +1,316 @@
+-- | The simplest, stupidest possible simplex algorithm.
+-- The idea here is to be slow, but "obviously correct" so other algorithms
+-- can be verified against it.
+--
+-- That's the plan, at least. For now this is just a first cut of trying to implement simplex.
+--
+module Numeric.Limp.Solve.Simplex.Maps
+    where
+import Numeric.Limp.Rep
+
+import Numeric.Limp.Solve.Simplex.StandardForm
+
+import Control.Arrow
+import qualified Data.Map as M
+import Data.Function (on)
+import Data.List (minimumBy, sortBy)
+
+
+-- | Result of a single pivot attempt
+data IterateResult z r c
+    -- | Maximum reached!
+    = Done
+    -- | Pivot was made
+    | Progress (Standard z r c)
+    -- | No progress can be made: unbounded along the objective
+    | Stuck
+
+deriving instance (Show z, Show r, Show (R c)) => Show (IterateResult z r c)
+
+
+-- | Try to find a pivot and then perform it.
+-- We're assuming, at this stage, that the existing solution is feasible.
+simplex1 :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> IterateResult z r c
+simplex1 s
+ -- Check if there are any positive columns in the objective:
+ = case pivotCols of
+    -- if there are none, we are already at the maximum
+    []
+     -> Done
+    -- there are some; try to find the first pivot row that works
+    _
+     -> go pivotCols
+ where
+
+  -- Check if there's any row worth pivoting on for this column.
+  -- We're trying to see if we can increase the value of this
+  -- column's variable from zero.
+  go ((pc,_):pcs)
+   = case pivotRowForCol s pc of
+       Nothing -> go pcs
+       Just pr
+        -> Progress
+        -- Perform the pivot.
+        -- This moves the variable pr out of the basis, and pc into the basis.
+         $ pivot s (pr,pc)
+
+  -- We've tried all the pivot columns and failed.
+  -- This means there's no edge we can take to increase our objective,
+  -- so it must be unbounded.
+  go []
+   = Stuck
+
+
+  -- We want to find some positive column from the objective.
+  -- In fact, find all of them and order descending.
+  pivotCols
+   = let ls  = M.toList $ fst $ _objective s
+         kvs = sortBy (compare `on` (negate . snd)) ls
+     in  filter ((>0) . snd) kvs
+
+
+-- | Find pivot row for given column.
+-- We're trying to find a way to increase the value of
+-- column from zero, and the returned row will be decreased to zero.
+-- Since all variables are >= 0, we cannot return a row that would set the column to negative.
+pivotRowForCol :: (Ord z, Ord r, Rep c)
+        => Standard z r c
+        -> StandardVar z r
+        -> Maybe (StandardVar z r)
+pivotRowForCol s col
+ = fmap   fst
+ $ minBy' (compare `on` snd)
+ $ concatMap (\(n,r)
+           -> let rv = lookupRow r col
+                  o  = objOfRow  r
+              in if    rv > 0
+                 then [(n, o / rv)]
+                 else [])
+ $ M.toList
+ $ _constraints s
+
+-- | Find minimum, or nothing if empty
+minBy' :: (a -> a -> Ordering) -> [a] -> Maybe a
+minBy' _ []
+ = Nothing
+minBy' f ls
+ = Just $ minimumBy f ls
+
+
+-- | Perform pivot for given row and column.
+-- We normalise row so that row.column = 1
+--
+-- > norm = row / row[column]
+--
+-- Then, for all other rows including the objective,
+-- we want to make sure its column entry is zero:
+--
+-- > row' = row - row[column]*norm
+--
+-- In the end, this means "column" will be an identity column, or a basis column.
+--
+pivot   :: (Ord z, Ord r, Rep c)
+        => Standard z r c
+        -> (StandardVar z r, StandardVar z r)
+        -> Standard z r c
+pivot s (pr,pc)
+ = let norm = normaliseRow
+       -- All other rows
+       rest = filter ((/=pr) . fst) $ M.toList $ _constraints s
+   in Standard
+    { _constraints = M.fromList ((pc, norm) : map (id *** fixup norm) rest)
+    , _objective   = fixup norm $ _objective s
+    , _substs      = _substs s }
+ where
+  -- norm = row / row[column]
+  normaliseRow
+   | Just row@(rm, ro) <- M.lookup pr $ _constraints s
+   = let c' = lookupRow row pc
+     in  (M.map (/c') rm, ro / c')
+
+   -- Pivot would not be chosen if row doesn't exist..
+   | otherwise
+   = (M.empty, 0)
+
+  -- row' = row - row[column]*norm
+  fixup (nm,no) row@(rm,ro)
+   = let co = lookupRow row pc
+     in  {- row' = row - co*norm -}
+         ( M.unionWith (+) rm (M.map ((-co)*) nm)
+         , ro - co * no )
+
+
+-- | Single phase of simplex.
+-- Keep repeating until no progress can be made.
+single_simplex :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> Maybe (Standard z r c)
+single_simplex s
+ = case simplex1 s of
+    Done        -> Just     s
+    Progress s' -> single_simplex  s'
+    Stuck       -> Nothing
+
+
+-- | Two phase:
+--  first, find a satisfying solution.
+--  then, solve simplex as normal.
+simplex
+        :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> Maybe (Standard z r c)
+simplex s
+ =   find_initial_sat s
+ >>= single_simplex
+
+-- | Find a satisfying solution.
+--   if there are any rows with negative values, this means their basic values are negative
+--   (which is not satisfying the x >= 0 constraint)
+--   these negative-valued rows must be pivoted around using modified pivot criteria
+find_initial_sat
+        :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> Maybe (Standard z r c)
+find_initial_sat s
+ = case negative_val_rows of
+    []      -> Just s
+    rs      -> go rs
+ where
+  -- Find all rows with negative values
+  -- because their current value is not feasible
+  negative_val_rows
+   = filter ((<0) . objOfRow . snd)
+   $ M.toList
+   $ _constraints s
+
+  -- Find largest negative (closest to zero) to pivot on:
+  -- pivoting on a negative will negate the value, setting it to positive
+  min_of_row (_,(rm,_))
+   = minBy' (compare `on` (negate . snd))
+   $ filter ((<0) . snd)
+   $ M.toList rm
+
+
+  -- There is no feasible solution
+  go []
+   = Nothing
+
+  -- Try pivoting on the rows 
+  go (r:rs)
+   | Just (pc,_) <- min_of_row r
+   , Just  pr    <- pivotRowForNegatives pc
+   = simplex
+   $ pivot s (pr, pc)
+
+   | otherwise
+   = go rs
+
+  -- opposite of pivotRowForCol...
+  pivotRowForNegatives col
+   = fmap   fst
+   $ minBy' (compare `on` (negate . snd))
+   $ concatMap (\(n,r)
+             -> let rv = lookupRow r col
+                    o  = objOfRow  r
+                in if    rv < 0
+                   then [(n, o / rv)]
+                   else [])
+   $ M.toList
+   $ _constraints s
+
+
+  
+
+-- Get map of each constraint's value
+assignmentAll :: (Rep c)
+        => Standard z r c
+        -> (M.Map (StandardVar z r) (R c), R c)
+assignmentAll s
+ = ( M.map    val (_constraints s)
+   , objOfRow     (_objective  s))
+ where
+  val (_, v)
+   = v
+
+-- Perform reverse substitution on constraint values
+-- to get original values (see StandardForm)
+assignment
+        :: (Ord z, Ord r, Rep c)
+        => Standard z r c
+        -> (Assignment () (Either z r) c, R c)
+assignment s
+ = ( Assignment M.empty $ M.union vs' rs'
+   , o )
+ where
+  (vs, o) = assignmentAll s
+
+  vs'     = M.fromList
+          $ concatMap only_svs
+          $ M.toList vs
+
+  rs'     = M.map eval $ _substs s
+
+  eval (lin,co)
+          = M.fold (+) co
+          $ M.mapWithKey (\k r -> r * (maybe 0 id $ M.lookup k vs))
+          $ lin
+
+  only_svs (SV v, val)
+   = [(v, val)]
+  only_svs _
+   = []
+
+
+
+-- Junk ---------------
+
+-- | Minimise whatever variables are 'basic' in given standard
+-- input must not already have an objective row "SVO",
+-- because the existing objective is added as a new row with that name
+minimise_basics
+        :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> Standard z r c
+minimise_basics s
+ = s
+ { _objective   = (M.map (const (1)) $ _constraints s, 0)
+ , _constraints = M.insert SVO (_objective s) (_constraints s)
+ }
+
+-- | Find the basic variables and "price them out" of the objective function,
+-- by subtracting multiples of the basic row from objective
+pricing_out 
+        :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> Standard z r c
+pricing_out s
+ = s
+ { _objective = M.foldWithKey  go
+                    (_objective   s)
+                    (_constraints s)
+ }
+ where
+  go v row@(rm,ro) obj@(om,oo)
+   | coeff <- lookupRow obj v
+   , coeff /= 0
+   , rowv  <- lookupRow row v
+   , mul   <- -(coeff / rowv)
+   = -- rowv = 1
+     -- obj' = obj - (coeff/rowv)*row
+     ( M.unionWith (+) om (M.map (mul*) rm)
+     , oo + mul*ro )
+   | otherwise
+   = obj
+
+-- | Pull the previously-hidden objective out of constraints, and use it
+drop_fake_objective
+        :: (Ord z, Ord r, Rep c)
+        => Standard z r c -> Standard z r c
+drop_fake_objective s
+ | cs     <- _constraints s
+ , Just o <- M.lookup SVO cs
+ = s
+ { _objective   = o
+ , _constraints = M.delete SVO cs }
+
+ | otherwise
+ = s
+
+
+
diff --git a/src/Numeric/Limp/Solve/Simplex/StandardForm.hs b/src/Numeric/Limp/Solve/Simplex/StandardForm.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Limp/Solve/Simplex/StandardForm.hs
@@ -0,0 +1,227 @@
+-- | Standard form for programs: only equalities and all variables >= 0
+-- To convert an arbitrary program to this form, we need to:
+--
+-- Convert unconstrained (-inf <= x <= +inf) variable into two separate parts, x+ and x-
+--  wherever x occurs, it will be replaced with "x+" - "x-".
+--
+-- Convert variables with non-zero lower bounds (c <= x) to a new variable x', so that
+--  x = x' + c
+--
+-- The opposite of these conversions must be performed when extracting a variable assignment
+-- from the solved program.
+--
+-- All constraints are converted into a less-than with a constant on the right, and then
+-- these less-than constraints (f <= c) have a slack variable s added such that
+--  f + s == c && s >= 0
+--
+module Numeric.Limp.Solve.Simplex.StandardForm
+    where
+import Numeric.Limp.Rep
+import Numeric.Limp.Canon.Constraint
+import Numeric.Limp.Canon.Linear
+import qualified Numeric.Limp.Canon.Program as C
+
+import qualified Data.Map as M
+import qualified Data.Set as S
+
+
+-- | A single linear function with a constant summand
+type StandardRow z r c
+    = (StandardLinear z r c, R c)
+
+-- | Entire program in standard form, as well as substitutions required to extract an assignment
+data Standard z r c
+    = Standard
+    { _objective   :: StandardRow z r c
+    , _constraints :: M.Map (StandardVar z r) (StandardRow z r c)
+    , _substs      :: StandardSubst z r c
+    }
+deriving instance (Show z, Show r, Show (R c)) => Show (Standard z r c)
+
+type StandardSubst  z r c
+    = M.Map (Either z r) (StandardRow z r c)
+
+type StandardLinear z r c
+    = M.Map (StandardVar z r) (R c)
+
+data StandardVar z r
+    -- | A normal variable
+    = SV (Either z r)
+
+    -- | A slack variable, introduced to make less-eq constraints into equalities
+    | SVS Int
+    -- | Magic objective, used when hiding an existing objective as a constraint
+    -- and creating a new objective
+    | SVO 
+
+    -- | When a variable has a lower bound other than 0, we replace all occurences with
+    -- with a new version minus the lower bound.
+    -- x >= 5
+    -- ==>
+    -- Lx - 5 >= 5
+    -- ==>
+    -- Lx >= 0
+    | SVLower (Either z r)
+
+    -- | When unconstrained variables are encountered, they are replaced with
+    -- x = SVPos x - SVNeg x
+    -- so both parts can be constrained to >= 0.
+    | SVPos (Either z r)
+    | SVNeg (Either z r)
+    deriving (Eq, Ord, Show)
+
+
+-- | Sum a list of linear functions together
+addLinears
+    :: (Ord z, Ord r, Rep c)
+    => [(StandardLinear z r c, R c)] -> (StandardLinear z r c, R c)
+addLinears []
+ = (M.empty, 0)
+addLinears ((lin,co):rs)
+ = let (lin',co') = addLinears rs
+   in  (M.unionWith (+) lin lin', co + co')
+
+
+-- | Perform substitution over a linear function/row
+substLinear
+    :: (Ord z, Ord r, Rep c)
+    => StandardSubst z r c -> (StandardLinear z r c, R c) -> (StandardLinear z r c, R c)
+substLinear sub (lin, co)
+ = let (lin', co') = addLinears 
+                   $ map subby 
+                   $ M.toList lin
+   in (lin', co + co')
+ where
+  subby (var, coeff)
+   = case var of
+      SV s
+       | Just (vs,cnst) <- M.lookup s sub
+       -> (M.map (*coeff) vs, -cnst * coeff)
+      _
+       -> (M.fromList [(var, coeff)], 0)
+
+
+-- | Convert canon program into standard form
+standard :: (Ord z, Ord r, Rep c)
+        => C.Program z r c
+        -> Standard z r c
+standard p
+ = Standard
+ { _objective   = objective
+ , _constraints = constraints
+ , _substs      = substs }
+ where
+  fv = C.varsOfProgram p
+  bs = C._bounds p
+
+  -- Objective is just negated
+  objective
+   = substLinear substs
+    ( M.map negate
+      $ standardOfLinear $ C._objective p
+    , 0)
+
+  -- Constraints are created for original program's bounds and constraints
+  -- and substitution is performed.
+  -- Each constraint/row receives its own slack variable.
+  constraints
+   = M.fromList
+   $ zipWith (\c s -> (s, substLinear substs $ c s))
+   ( constrs ++ bounds )
+   ( map SVS [1..] )
+
+  -- Union of all substitutions
+  substs
+   = M.fromList
+   $ concatMap substOf
+   $ S.toList fv
+
+  -- Substitution for "x" ==> "x+" - "x-"
+  negPos v
+   = [(v, (M.fromList [(SVPos v, 1), (SVNeg v, -1)], 0))]
+
+  -- Look at bounds of variables and decide
+  substOf v
+   = case M.lookup v bs of
+     -- Unconstrained, so it can be negative
+     Nothing
+      -> negPos v
+     Just (Nothing, Nothing)
+      -> negPos v
+     Just (Just 0, _)
+      -> []
+     -- Nonzero lower bound, so replace: v = v' + n
+     Just (Just n, _)
+      -> [(v, (M.fromList [(SVLower v, 1)], n)) ]
+     _
+      -> []
+
+  bounds
+   = concatMap linearOfBound
+   $ M.toList
+   $ C._bounds p
+
+  linearOfBound (v,binds)
+   = case binds of
+     (_, Just n)
+      -> [\s -> (M.fromList [(SV v, 1), (s, 1)], n)]
+     _
+      -> []
+
+  Constraint cs = C._constraints p
+  constrs
+   = concatMap linearOfConstraint cs
+  linearOfConstraint (C1 lo lin up)
+   = let lin' = standardOfLinear lin
+   in case (lo,up) of
+      (Nothing,Nothing)
+       -> []
+      (Just lo', Nothing)
+       -> [ lt lo' lin' ]
+      (Nothing,  Just up')
+       -> [ gt up' lin' ]
+      (Just lo', Just up')
+       -> [ lt lo' lin'
+          , gt up' lin' ]
+
+
+  lt lo lin s
+   = ( M.union (M.map negate lin) (M.fromList [(s,1)])
+     , negate lo )
+  gt up lin s
+   = ( M.union lin (M.fromList [(s, 1)])
+     , up )
+
+  standardOfLinear (Linear lin)
+   = M.mapKeysMonotonic SV lin
+
+
+--- 5 <= x1 <= 40
+-- ==>
+-- x1 subst Lx1+5
+-- Lx1 + 5 <= 40
+-- ==>
+-- Lx1 <= 35
+
+-- assignmentOfMap :: Standard z r c -> M.Map (StandardVar z r) (R c) -> Assignment z r c
+
+
+
+-- Simple helpers ----------
+
+-- | Get the coefficient of a variable in given row
+lookupRow :: (Ord z, Ord r, Rep c)
+    => StandardRow z r c
+    -> StandardVar z r
+    -> R c
+lookupRow (r,_) v
+ = case M.lookup v r of
+    Nothing -> 0
+    Just vv -> vv
+
+-- | Get objective or basis value of a row
+objOfRow
+    :: StandardRow z r c
+    -> R c
+objOfRow = snd
+
diff --git a/tests/Arbitrary/Assignment.hs b/tests/Arbitrary/Assignment.hs
new file mode 100644
--- /dev/null
+++ b/tests/Arbitrary/Assignment.hs
@@ -0,0 +1,33 @@
+module Arbitrary.Assignment where
+
+import Numeric.Limp.Rep
+
+import Arbitrary.Var
+
+import Test.QuickCheck
+import Control.Applicative
+import Data.Map (fromList)
+
+type Assignment' = Assignment ZVar RVar IntDouble
+
+instance Arbitrary (Z IntDouble) where
+ arbitrary = Z <$> arbitrary
+
+instance Arbitrary (R IntDouble) where
+ arbitrary = R <$> (fromIntegral <$> (arbitrary :: Gen Int))
+
+
+instance Arbitrary (Assignment ZVar RVar IntDouble) where
+ arbitrary = arbitrary >>= assignment
+
+
+assignment :: Vars -> Gen Assignment'
+assignment (Vars zs rs)
+ = do   zs' <- listOf (elements zs)
+        zvs <- infiniteListOf arbitrary
+
+        rs' <- listOf (elements rs)
+        rvs <- infiniteListOf arbitrary
+
+        return $ Assignment (fromList $ zs' `zip` zvs) (fromList $ rs' `zip` rvs)
+
diff --git a/tests/Arbitrary/Program.hs b/tests/Arbitrary/Program.hs
new file mode 100644
--- /dev/null
+++ b/tests/Arbitrary/Program.hs
@@ -0,0 +1,82 @@
+module Arbitrary.Program where
+
+import qualified Numeric.Limp.Program as P
+import Numeric.Limp.Rep
+
+import Arbitrary.Var
+import Arbitrary.Assignment
+
+import Test.QuickCheck
+import Control.Applicative
+
+type Program' = P.Program ZVar RVar IntDouble
+
+data ProgramAss = ProgramAss Program' Assignment'
+ deriving Show
+
+instance Arbitrary ProgramAss where
+ arbitrary
+  = do  a <- arbitrary
+        ProgramAss <$> program a <*> assignment a
+
+instance Arbitrary Program' where
+ arbitrary = arbitrary >>= program
+
+
+program :: Vars -> Gen Program'
+program vs
+ = do   dir  <- elements [P.Minimise, P.Maximise]
+        
+        obj  <- linearR        vs
+        cons <- constraints    vs
+        bnds <- listOf (bounds vs)
+
+        return $ P.program dir obj cons bnds
+
+
+linearR :: Vars -> Gen (P.Linear ZVar RVar IntDouble P.KR)
+linearR (Vars zs rs)
+ = do   let vs = map Left zs ++ map Right rs
+        vs' <- listOf1 (elements vs)
+        cs' <- infiniteListOf arbitrary
+        summand <- arbitrary
+        return $ P.LR (vs' `zip` cs') summand
+
+linearZ :: Vars -> Gen (P.Linear ZVar RVar IntDouble P.KZ)
+linearZ (Vars zs _rs)
+ = do   vs' <- listOf1 (elements zs)
+        cs' <- infiniteListOf arbitrary
+        summand <- arbitrary
+        return $ P.LZ (vs' `zip` cs') summand
+
+
+constraints :: Vars -> Gen (P.Constraint ZVar RVar IntDouble)
+constraints vs
+ = oneof
+ [ (P.:==)   <$> lR <*> lR
+ , (P.:<=)   <$> lR <*> lR
+ , (P.:<)    <$> lZ <*> lZ
+ , (P.:>=)   <$> lR <*> lR
+ , (P.:>)    <$> lZ <*> lZ
+ , P.Between <$> lR <*> lR <*> lR
+ , (P.:&&)   <$> constraints vs <*> constraints vs
+ , return P.CTrue ]
+ where
+  lR = linearR vs
+  lZ = linearZ vs
+
+
+bounds :: Vars -> Gen (P.Bounds ZVar RVar IntDouble)
+bounds (Vars zs rs)
+ = oneof [bZ, bR]
+ where
+  bZ = do   v <- elements zs
+            a <- arbitrary
+            b <- arbitrary
+            return $ P.BoundZ (a,v,b)
+
+  bR = do   v <- elements rs
+            a <- arbitrary
+            b <- arbitrary
+            return $ P.BoundR (a,v,b)
+
diff --git a/tests/Arbitrary/Var.hs b/tests/Arbitrary/Var.hs
new file mode 100644
--- /dev/null
+++ b/tests/Arbitrary/Var.hs
@@ -0,0 +1,38 @@
+module Arbitrary.Var where
+import Test.QuickCheck
+
+data ZVar = ZVar String
+ deriving (Eq,Ord)
+
+instance Show ZVar where
+ show (ZVar z) = "z$" ++ z
+
+instance Arbitrary ZVar where
+ arbitrary
+        -- 26 variables should be enough for anyone!
+  = do  c <- elements ['a'..'z']
+        return $ ZVar [c]
+
+
+data RVar = RVar String
+ deriving (Eq,Ord)
+
+instance Show RVar where
+ show (RVar r) = "r$" ++ r
+
+instance Arbitrary RVar where
+ arbitrary
+  = do  c <- elements ['a'..'z']
+        return $ RVar [c]
+
+
+data Vars = Vars [ZVar] [RVar]
+ deriving Show
+
+instance Arbitrary Vars where
+ arbitrary
+  = do  NonEmpty zs   <- arbitrary :: Gen (NonEmptyList ZVar)
+        NonEmpty rs   <- arbitrary :: Gen (NonEmptyList RVar)
+        return $ Vars zs rs
+
+
diff --git a/tests/BranchExample.hs b/tests/BranchExample.hs
new file mode 100644
--- /dev/null
+++ b/tests/BranchExample.hs
@@ -0,0 +1,94 @@
+module BranchExample where
+
+import Numeric.Limp.Rep.Rep     as R
+import Numeric.Limp.Rep.Arbitrary     as R
+import Numeric.Limp.Program as P
+import Numeric.Limp.Canon   as C
+import Numeric.Limp.Solve.Simplex.Maps   as SM
+import Numeric.Limp.Solve.Simplex.StandardForm   as ST
+import Numeric.Limp.Solve.Branch.Simple  as B
+
+import Numeric.Limp.Canon.Pretty
+import Debug.Trace
+
+import Control.Applicative
+
+-- Dead simple ones -------------------------
+-- x = 2
+prog1 :: P.Program String String R.Arbitrary
+prog1
+ = P.maximise
+    -- objective
+        (z "x" 1)
+    -- subject to
+     (   z "x"  2 :<= con 5
+     :&& z "x"  4 :>= con 7)
+    []
+
+-- x = 1, y = 2
+prog2 :: P.Program String String R.Arbitrary
+prog2
+ = P.minimise
+    -- objective
+        (z "x" 1 .+. z "y" 1)
+    -- subject to
+     (   z "x"  2 :<= con 5 -- z "y" 1 .+. con 1
+     :&& z "x"  1 :>= con 1 
+     :&& z "y"  1 :<= con 4
+     :&& z "y"  1 :>= con 1)
+    [ lowerZ 0 "x" 
+    , lowerZ 0 "y" ]
+
+
+xkcd :: Direction -> P.Program String String R.Arbitrary
+xkcd dir = P.program dir
+           ( z1 mf .+.
+             z1 ff .+.
+             z1 ss .+.
+             z1 hw .+.
+             z1 ms .+.
+             z1 sp )
+           ( z mf mfp .+.
+             z ff ffp .+.
+             z ss ssp .+.
+             z hw hwp .+.
+             z ms msp .+.
+             z sp spp :== con 1505 )
+           [ lowerZ 0 mf
+           , lowerZ 0 ff
+           , lowerZ 0 ss
+           , lowerZ 0 hw
+           , lowerZ 0 ms
+           , lowerZ 0 sp
+            ]
+  where
+    (mf, mfp) = ("mixed-fruit",       215)
+    (ff, ffp) = ("french-fries",      275)
+    (ss, ssp) = ("side-salad",        335)
+    (hw, hwp) = ("hot-wings",         355)
+    (ms, msp) = ("mozzarella-sticks", 420)
+    (sp, spp) = ("sampler-plate",     580)
+
+test :: (Show z, Show r, Ord z, Ord r)
+     => P.Program z r R.Arbitrary -> IO ()
+test prog
+ = let prog' = C.program prog
+       
+       simpl p = SM.simplex $ ST.standard p
+
+       solver p
+        | st <- ST.standard p
+        -- , trace (ppr show show p) True
+        , Just s' <- SM.simplex st
+        -- , trace ("SAT") True
+        , ass <- SM.assignment s'
+         = Just ass
+        | otherwise
+        -- , trace ("unsat") True
+        = Nothing
+       bb    = B.branch solver
+   in  do   
+            -- putStrLn (show (simpl prog'))
+            -- putStrLn (show (solver prog'))
+            putStrLn (show (bb prog'))
+
diff --git a/tests/Convert.hs b/tests/Convert.hs
new file mode 100644
--- /dev/null
+++ b/tests/Convert.hs
@@ -0,0 +1,19 @@
+module Convert where
+
+import Numeric.Limp.Program as P
+import Numeric.Limp.Canon   as C
+
+import Arbitrary.Program
+import Data.Monoid
+
+import Test.Tasty.QuickCheck
+import Test.Tasty.TH
+
+
+tests = $(testGroupGenerator)
+
+prop_constraints_converted :: ProgramAss -> Bool
+prop_constraints_converted (ProgramAss p a)
+ =  P.checkProgram a  p
+ == C.checkProgram a (C.program p)
+
diff --git a/tests/SimplexExample.hs b/tests/SimplexExample.hs
new file mode 100644
--- /dev/null
+++ b/tests/SimplexExample.hs
@@ -0,0 +1,55 @@
+module SimplexExample where
+
+import Numeric.Limp.Rep     as R
+import Numeric.Limp.Program as P
+import Numeric.Limp.Canon   as C
+import Numeric.Limp.Solve.Simplex.Maps      as SM
+import Numeric.Limp.Solve.Simplex.StandardForm  as ST
+
+import Control.Monad
+import qualified Data.Map as M
+
+
+data Xs = X1 | X2 | X3
+ deriving (Eq, Ord, Show)
+
+prog :: P.Program () Xs R.IntDouble
+prog
+ = P.maximise
+    -- objective
+        (r X1 60 .+. r X2  30 .+. r X3  20)
+    -- subject to
+     (   r X1  8 .+. r X2   6 .+. r X3   1 :<= con 48
+     :&& r X1  2 .+. r X2 1.5 .+. r X3 0.5 :<= con  8
+     :&& r X1  4 .+. r X2   2 .+. r X3 1.5 :<= con 20
+     :&&             r X2   1              :<= con  5)
+    -- bounds ommitted for now
+    [ lowerR 0 X1 , lowerR 0 X2 , lowerR 0 X3 ]
+    -- []
+
+test :: IO Bool
+test
+ = case SM.simplex $ ST.standard $ C.program prog of
+   Nothing
+    -> do   putStrLn "Error: simplex returned Nothing"
+            putStrLn (show $ ST.standard $ C.program prog)
+            putStrLn (show $ SM.simplex1 $ ST.standard $ C.program prog)
+            return False
+
+   Just s
+    -> do   let (Assignment _ vars,obj) = SM.assignment s
+            let vars'      = M.toList vars
+            let e_vars = [(Right X1, 2.0), (Right X3, 8.0)] :: [(Either () Xs, R IntDouble)]
+            let e_obj  = -280
+            putStrLn "Vars:"
+            putStrLn (show vars')
+            putStrLn "Obj:"
+            putStrLn (show obj)
+
+            when (obj /= e_obj) $
+                putStrLn ("Bad objective: should be " ++ show e_obj)
+            when (vars' /= e_vars) $
+                putStrLn ("Bad vars: should be "      ++ show e_vars)
+
+            return (obj == e_obj && vars' == e_vars)
+
diff --git a/tests/Simplexs.hs b/tests/Simplexs.hs
new file mode 100644
--- /dev/null
+++ b/tests/Simplexs.hs
@@ -0,0 +1,311 @@
+module Simplexs where
+
+import Numeric.Limp.Rep.Rep     as R
+import Numeric.Limp.Rep.Arbitrary     as R
+import Numeric.Limp.Program as P
+import Numeric.Limp.Canon   as C
+import Numeric.Limp.Solve.Simplex.Maps      as SM
+import Numeric.Limp.Solve.Simplex.StandardForm  as ST
+
+import qualified Data.Map as M
+
+
+data Xs = X1 | X2 | X3
+ deriving (Eq, Ord, Show)
+
+-- Dead simple ones -------------------------
+-- x1 = 10
+prog1 :: P.Program () Xs R.Arbitrary
+prog1
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 10)
+    -- bounds omitted for now
+    []
+
+-- x1 = 10
+prog2 :: P.Program () Xs R.Arbitrary
+prog2
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 10)
+    [ lowerR 0 X1 ]
+
+-- x1 = 0
+prog3 :: P.Program () Xs R.Arbitrary
+prog3
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 10)
+    [ lowerR 0 X1 ]
+
+-- Unbounded!
+prog4 :: P.Program () Xs R.Arbitrary
+prog4
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 10)
+    []
+
+
+-- Two constraints! --------------
+
+-- x = 10
+prog5 :: P.Program () Xs R.Arbitrary
+prog5
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 10
+     :&& r X1  1 :>= con (-10))
+    []
+
+-- x = -10
+prog6 :: P.Program () Xs R.Arbitrary
+prog6
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 10
+     :&& r X1  1 :>= con (-10))
+    []
+
+
+-- Now two variables -------------
+-- x1 = 20, x2 = 10
+prog7 :: P.Program () Xs R.Arbitrary
+prog7
+ = P.maximise
+    -- objective
+        (r X1 1 .+. r X2 1)
+    -- subject to
+     (   r X1  1 :<= r X2 2
+     :&& r X2  1 :<= con 10)
+    [lowerR 0 X1, lowerR 0 X2]
+
+-- x1 = 20, x2 = 10
+prog8 :: P.Program () Xs R.Arbitrary
+prog8
+ = P.maximise
+    -- objective
+        (r X1 1 .+. r X2 1)
+    -- subject to
+     (   r X1  1 :<= r X2 2
+     :&& r X2  1 :<= con 10)
+    [] -- [lowerR 0 X1, lowerR 0 X2]
+
+-- Something where vars=0 isn't sat ------
+-- x1 = 8
+prog9 :: P.Program () Xs R.Arbitrary
+prog9
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :>= con 8 
+     :&& r X1  1 :<= con 10)
+    [lowerR 0 X1]
+
+-- x1 = 10
+prog10 :: P.Program () Xs R.Arbitrary
+prog10
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :>= con 8 
+     :&& r X1  1 :<= con 10)
+    [lowerR 0 X1]
+
+
+
+-- An equality constraint ------------
+-- x1 = 10
+prog11 :: P.Program () Xs R.Arbitrary
+prog11
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :== con 10 )
+    [lowerR 0 X1]
+
+-- x1 = 10
+prog12 :: P.Program () Xs R.Arbitrary
+prog12
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :== con 10 )
+    [lowerR 0 X1]
+
+
+-- From wikipedia ----------------
+-- x1 = 2.142..., x3 = 3.571...
+prog13 :: P.Program () Xs R.Arbitrary
+prog13
+ = P.minimise
+    -- objective
+        (r X1 (-2) .+. r X2 (-3) .+. r X3 (-4))
+    -- subject to
+     (   r X1  3   .+. r X2 2    .+. r X3 1 :== con 10
+     :&& r X1  2   .+. r X2 5    .+. r X3 3 :== con 15)
+    [lowerR 0 X1
+    ,lowerR 0 X2
+    ,lowerR 0 X3]
+
+-- x1 = 1.818..., x2 = 2.272...
+prog14 :: P.Program () Xs R.Arbitrary
+prog14
+ = P.maximise
+    -- objective
+        (r X1 (-2) .+. r X2 (-3) .+. r X3 (-4))
+    -- subject to
+     (   r X1  3   .+. r X2 2    .+. r X3 1 :== con 10
+     :&& r X1  2   .+. r X2 5    .+. r X3 3 :== con 15)
+    [lowerR 0 X1
+    ,lowerR 0 X2
+    ,lowerR 0 X3]
+
+-- An equality constraint on unconstrained (+-) ------------
+-- x1 = 10
+prog15 :: P.Program () Xs R.Arbitrary
+prog15
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :== con 10 )
+    []
+
+-- A lower bound greater than zero ------------
+-- x1 = 5
+prog16 :: P.Program () Xs R.Arbitrary
+prog16
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 30 )
+    [lowerR 5 X1]
+
+-- Lower and upper bounds -------
+-- x1 = 5
+prog17 :: P.Program () Xs R.Arbitrary
+prog17
+ = P.minimise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 30 )
+    [lowerUpperR 5 X1 10]
+-- x1 = 10
+prog18 :: P.Program () Xs R.Arbitrary
+prog18
+ = P.maximise
+    -- objective
+        (r X1 1)
+    -- subject to
+     (   r X1  1 :<= con 30 )
+    [lowerUpperR 5 X1 10]
+
+-- x1 = 1, x2 = 2
+prog19 :: P.Program () Xs R.Arbitrary
+prog19
+ = P.minimise
+    (r X1 1 .+. r X2 1)
+    (    r X1 2 :<= r X2 1
+    :&&  r X1 1 :>= con 1)
+    [ lowerR 0 X1
+    , lowerR 0 X2]
+
+
+-- error uncovered by branch -------
+-- x1 = 1
+-- x2 = 1.870...
+prog20 :: P.Program () Xs R.Arbitrary
+prog20
+ = P.minimise
+    -- x1 = mozzarella
+    -- x2 = sampler plate
+    (r1 X1 .+. r1 X2)
+    (r X1 420 .+. r X2 580 :== con 1505)
+    [ lowerR 1 X1
+    , lowerUpperR 0 X2 2 ]
+
+{-
+Minimize
+	1.0 "french-fries" + 1.0 "hot-wings" + 1.0 "mixed-fruit" + 1.0 "mozzarella-sticks" + 1.0 "sampler-plate" + 1.0 "side-salad"
+Subject to
+	-275.0 "french-fries" - 355.0 "hot-wings" - 215.0 "mixed-fruit" - 420.0 "mozzarella-sticks" - 580.0 "sampler-plate" - 335.0 "side-salad" >= -1505.0
+	-275.0 "french-fries" - 355.0 "hot-wings" - 215.0 "mixed-fruit" - 420.0 "mozzarella-sticks" - 580.0 "sampler-plate" - 335.0 "side-salad" <= -1505.0
+
+Bounds
+	0.0 <= "french-fries"
+	0.0 <= "hot-wings"
+	0.0 <= "mixed-fruit"
+	1.0 <= "mozzarella-sticks"
+	0.0 <= "sampler-plate" <= 2.0
+	0.0 <= "side-salad"
+-}
+
+-- nonzero lower bound with non-1 coeff
+-- x1 = 2.5
+prog21 :: P.Program () Xs R.Arbitrary
+prog21
+ = P.minimise
+    (r1 X1)
+    (r X1 2 :>= con 5)
+    [ lowerR 1 X1 ]
+
+-- eq bound with non-1 coeff
+-- x1 = 1, x2 = 3
+prog22 :: P.Program () Xs R.Arbitrary
+prog22
+ = P.minimise
+    (r1 X1 .+. r1 X2)
+    (r X1 2 .+. r X2 1 :>= con 5)
+    [ lowerUpperR 1 X1 1
+    , lowerR 0 X2]
+
+
+std :: (Ord z, Ord r, Rep c) => P.Program z r c -> Standard z r c
+std = ST.standard . C.program
+
+
+
+
+test :: P.Program () Xs R.Arbitrary -> IO Bool
+test p
+ = case SM.simplex $ ST.standard $ C.program p of
+   Nothing
+    -> do   putStrLn "Error: simplex returned Nothing"
+            putStrLn (show $ ST.standard $ C.program p)
+            putStrLn (show $ SM.simplex1 $ ST.standard $ C.program p)
+            return False
+
+   Just s
+    -> do   let (Assignment _ vars,obj) = SM.assignment s
+            let vars'      = M.toList vars
+
+            putStrLn (show $ ST.standard $ C.program p)
+            putStrLn (show $ SM.simplex1 $ ST.standard $ C.program p)
+
+            putStrLn "Vars:"
+            putStrLn (show vars')
+            putStrLn "Obj:"
+            putStrLn (show obj)
+
+            return True
+
diff --git a/tests/Simplify.hs b/tests/Simplify.hs
new file mode 100644
--- /dev/null
+++ b/tests/Simplify.hs
@@ -0,0 +1,110 @@
+module Simplify where
+
+import Numeric.Limp.Program as P
+import Numeric.Limp.Canon   as C
+import Numeric.Limp.Canon.Simplify as CS
+import Numeric.Limp.Canon.Simplify.Subst as CS
+import Numeric.Limp.Canon.Simplify.Bounder as CS
+import Numeric.Limp.Canon.Simplify.Crunch as CS
+
+import Numeric.Limp.Canon.Pretty
+
+import Arbitrary.Assignment     as Arb
+import Arbitrary.Var            as Arb
+import Arbitrary.Program        as Arb
+import Data.Monoid
+
+import Test.Tasty.QuickCheck
+import Test.Tasty.TH
+
+import Debug.Trace
+
+tests = $(testGroupGenerator)
+
+prop_bounder :: ProgramAss -> Property
+prop_bounder (ProgramAss p a)
+ = let cp     = C.program p
+       cp'    = CS.bounderProgram cp
+       valcp  = C.checkProgram a cp
+       valcp'
+        | Right p' <- cp'
+        = C.checkProgram a p'
+        -- Infeasible, so assignment is false
+        | otherwise
+        = False
+   in  counterexample
+       (unlines
+         [ "CP: " ++ show cp
+         , "CP':" ++ show cp'
+         , "Val: " ++ show (valcp, valcp')])
+       $ valcp == valcp'
+
+
+prop_crunch :: ProgramAss -> Property
+prop_crunch (ProgramAss p a)
+ = let cp     = C.program p
+       cp'    = CS.crunchProgram cp
+       valcp  = C.checkProgram a cp
+       valcp' = C.checkProgram a cp'
+   in  counterexample
+       (unlines
+         [ "CP: " ++ show cp
+         , "CP':" ++ show cp'
+         , "Val: " ++ show (valcp, valcp')])
+       $ valcp == valcp'
+
+
+-- | I don't think this property is very interesting.
+-- The real property should be something like:
+--
+-- > solve cp == solve (simplify cp)
+--
+prop_simplify :: Program' -> Property
+prop_simplify p
+ = let cp = C.program p
+       simp = CS.simplify cp
+   in  case simp of
+       Left _
+        -> property True
+       Right (a', cp')
+        -> let valcp  = C.checkProgram a' cp
+               valcp' = C.checkProgram a' cp'
+           in  counterexample
+           (unlines
+             [ "CP: " ++ show cp
+             , "CP':" ++ show cp'
+             , "Ass:" ++ show a'
+             , "Val: " ++ show (valcp, valcp')])
+           $ valcp == valcp'
+
+
+prop_subst_linear :: Vars -> Property
+prop_subst_linear vs
+ = forAll (Arb.linearR    vs) $ \f ->
+   forAll (Arb.assignment vs) $ \a ->
+   forAll (Arb.assignment vs) $ \b ->
+     let (fc, _)   = C.linear f
+         (fc', c') = substLinear a fc
+     in  C.evalR (a <> b) fc == C.evalR b fc' + c'
+
+
+-- subst can actually make a failing program pass.
+-- so this test needs to be implication, not equivalence.
+prop_subst_program :: Vars -> Property
+prop_subst_program vs
+ = forAll (Arb.program    vs) $ \f ->
+   forAll (Arb.assignment vs) $ \a ->
+   forAll (Arb.assignment vs) $ \b ->
+     let fc    = C.program f
+         fc'   = substProgram a fc
+         both  = a <> b
+         valcp = C.checkProgram both fc 
+         valcp'= C.checkProgram b fc'
+     in counterexample 
+         (unlines
+         [ "CP: " ++ show fc
+         , "CP':" ++ show fc'
+         , "Ass:" ++ show both
+         , "Val: " ++ show (valcp, valcp')])
+       $ if valcp then valcp' else True
+