-- maxElim :: (Ord v, Eq c, IsNat c, SemiRing c) => Constraint (MaxPoly v c) -> [ConditionalConstraint (P.Polynomial v c)]
-- maxElim ieq = walk [([],ieq)] where
--   walk [] = []
--   walk ((ps,ieq):ceqs) =
--     case splitIEQ ieq of
--       Nothing -> CC { premises = ps, constraint = fmap toPoly ieq } : walk ceqs
--       Just (p,q,ctx) -> walk (cp:cq:ceqs)
--         where
--           p' = toPoly p
--           q' = toPoly q
--           cp = ((p' :>=: q')          : ps,ctx p)
--           cq = ((q' :>=: (p' .+ one)) : ps,ctx q)
                      
--   toPoly (Var v)    = P.variable v
--   toPoly (Const c)  = P.coefficient c
--   toPoly (Plus p q) = toPoly p .+ toPoly q
--   toPoly (Mult p q) = toPoly p .* toPoly q
--   toPoly (Max {})     = error "maxElim: polynomial still contains max"
  
--   splitIEQ (p1 :>=: p2) =
--     case (splitMaxPoly p1, splitMaxPoly p2) of
--       (Right (p,q,ctx), _              ) -> Just (p,q, \pi -> ctx pi :>=: p2)
--       (_              , Right (p,q,ctx)) -> Just (p,q, \pi -> p1 :>=: ctx pi)    
--       (_              , _              ) -> Nothing
--   splitBinOp c p1 p2 = 
--     case (splitMaxPoly p1, splitMaxPoly p2) of
--       (Right (p1',p2',ctx), _                  ) -> Right (p1',p2', \ pi' -> c (ctx pi') p2)
--       (_                  , Right (p1',p2',ctx)) -> Right (p1',p2', \ pi' -> c p1 (ctx pi'))
--       (_                  , _                  ) -> Left (c p1 p2)
--   splitMaxPoly p@(Var _) = Left p
--   splitMaxPoly p@(Const _) = Left p
--   splitMaxPoly (Max p1 p2)  =
--     case splitBinOp Max p1 p2 of {Left _ -> Right (p1, p2, id); s -> s}
--   splitMaxPoly (Plus p1 p2) = splitBinOp Plus p1 p2
--   splitMaxPoly (Mult p1 p2) = splitBinOp Mult p1 p2  
      

simp :: (Ord c, Ord v, IsNat c, SemiRing c) => MaxPoly v c -> MaxPoly v c
simp = fromPolyList . filterSubsumed . nub . splitMax where -- 
  fromPolyList [] = zero
  fromPolyList ps = maximumA (map fromPoly ps)

  fromPoly = P.fromPolynomial variable constant
  -- TODO  max(x0 + x1 + x2,1 + x0)
  filterSubsumed ps = foldr (\ p -> filter (not . subsumes p)) ps ps
  p1 `subsumes` p2 =
    and [ c1 >= c2 | (c1 :>=: c2) <- P.absolutePositive (p1 `P.minus` p2)]
    && p1 /= p2
simp :: (Ord c, Ord v, IsNat c, SemiRing c) => MaxPoly v c -> MaxPoly v c
simp = fromPolyList . filterSubsumed . nub . splitMax where -- 
  fromPolyList [] = zero
  fromPolyList ps = maximumA (map fromPoly ps)

  fromPoly = P.fromPolynomial variable constant
  -- TODO  max(x0 + x1 + x2,1 + x0)
  filterSubsumed ps = foldr (\ p -> filter (not . subsumes p)) ps ps
  p1 `subsumes` p2 =
    and [ c1 >= c2 | (c1 :>=: c2) <- P.absolutePositive (p1 `P.minus` p2)]
    && p1 /= p2
simp = fromPolyList . filterSubsumed . nub . splitMax where -- 
  fromPolyList [] = zero
  fromPolyList ps = maximumA (map fromPoly ps)

  fromPoly = P.fromPolynomial variable constant
  -- TODO  max(x0 + x1 + x2,1 + x0)
  filterSubsumed ps = foldr (\ p -> filter (not . subsumes p)) ps ps
  p1 `subsumes` p2 =
    and [ c1 >= c2 | (c1 :>=: c2) <- P.absolutePositive (p1 `P.minus` p2)]
    && p1 /= p2