Added equation pruning and tested most of it using QuickCheck

transformations/ArrayUsageCheck.hs

@@ -205,3 +205,48 @@ solveConstraints p ineq
 mygcd :: Integer -> Integer -> Integer
 mygcd 0 0 = 0
 mygcd x y = gcd x y
+
+
+-- | Prunes the inequalities.  It does what is described in section 2.3 of Pugh's ACM paper;
+-- it removes redundant inequalities, fails (evaluates to Nothing) if it finds a contradiction
+-- and turns any 2x + y <= 4, 2x + y >= 4 pairs into equalities.  The list of such equalities
+-- (which may well be an empty list) and the remaining inequalities is returned.
+pruneAndCheck :: InequalityProblem -> Maybe (EqualityProblem, InequalityProblem)
+pruneAndCheck ineq = let (opps,others) = splitEither $ groupOpposites $ map pruneGroup groupedIneq in
+                     mapM checkOpposite opps >>* (\x -> (x,others))
+  where
+    groupedIneq = groupBy (\x y -> EQ == coeffSort x y) $ sortBy coeffSort ineq
+
+    coeffSort :: InequalityConstraintEquation -> InequalityConstraintEquation -> Ordering
+    coeffSort x y = compare (tail $ elems x) (tail $ elems y)
+
+    -- | Takes in a group of inequalities with identical a_1 .. a_n coefficients
+    -- and returns the equation with the smallest unit coefficient.  Consider the standard equation:
+    -- a_1.x_1 + a_2.x_2 .. a_n.x_n >= -a_0.  We want one equation with the maximum value of -a_0
+    -- (this will be the strongest equation), which is therefore the minimum value of a_0.
+    -- This therefore automatically removes duplicate and redundant equations.
+    pruneGroup :: [InequalityConstraintEquation] -> InequalityConstraintEquation
+    pruneGroup = minimumBy (\x y -> compare (x ! 0) (y ! 0))
+
+    -- | Groups all equations with their opposites, if found.  Returns either a pair
+    -- or a singleton.  O(N^2), but there shouldn't be that many inequalities to process (<= 10, I expect).
+    -- Assumes equations have already been pruned, and that therefore for every unique a_1 .. a_n 
+    -- set, there is only one equation.
+    groupOpposites :: InequalityProblem -> [Either (InequalityConstraintEquation,InequalityConstraintEquation) InequalityConstraintEquation]
+    groupOpposites [] = []
+    groupOpposites (e:es) = case findOpposite e es of
+                              Just (opp,rest) -> (Left (e,opp)) : (groupOpposites rest)
+                              Nothing -> (Right e) : (groupOpposites es)
+
+    findOpposite :: InequalityConstraintEquation -> [InequalityConstraintEquation] -> Maybe (InequalityConstraintEquation,[InequalityConstraintEquation])
+    findOpposite _ [] = Nothing
+    findOpposite target (e:es) | negTarget == (tail $ elems e) = Just (e,es)
+                               | otherwise = case findOpposite target es of
+                                               Just (opp,rest) -> Just (opp,e:rest)
+                                               Nothing -> Nothing
+      where
+        negTarget = map negate $ tail $ elems target
+
+    checkOpposite :: (InequalityConstraintEquation,InequalityConstraintEquation) -> Maybe EqualityConstraintEquation
+    checkOpposite (x,y) | (x ! 0) == negate (y ! 0) = Just x -- doesn't matter which we pick to become the equality
+                        | otherwise = Nothing -- The inequalities can't hold

transformations/RainUsageCheckTest.hs

@@ -20,8 +20,10 @@ module RainUsageCheckTest (qcTests) where
 
 import Control.Monad.Identity
 import Data.Graph.Inductive
-import Data.Array
+import Data.Array.IArray
+import Data.List
 import qualified Data.Map as Map
+import Data.Maybe
 import qualified Data.Set as Set
 import Prelude hiding (fail)
 import Test.HUnit
@@ -391,11 +393,11 @@ newtype OmegaTestInput = OMI (EqualityProblem, InequalityProblem) deriving (Show
 
 -- | Generates an Omega test problem with between 1 and 10 variables (incl), where the solutions
 -- are numbers between -20 and + 20 (incl).
-generateProblem :: Gen OmegaTestInput
-generateProblem = (choose (1,10) >>= (\n -> replicateM n $ choose (-20,20)) >>= generateEqualities) >>* OMI
+generateProblem :: Gen (EqualityProblem, InequalityProblem)
+generateProblem = choose (1,10) >>= (\n -> replicateM n $ choose (-20,20)) >>= generateEqualities
 
 instance Arbitrary OmegaTestInput where
-  arbitrary = generateProblem
+  arbitrary = generateProblem >>* OMI
 
 qcOmegaEquality :: [QuickCheckTest]
 qcOmegaEquality = [scaleQC (40,200,2000,10000) prop]
@@ -406,6 +408,69 @@ qcOmegaEquality = [scaleQC (40,200,2000,10000) prop]
         omegaCheck (Just ineqs) = all (all (== 0) . elems) ineqs
         omegaCheck Nothing = False
 
+type MutatedEquation =
+  (InequalityProblem
+  ,[(EqualityConstraintEquation,EqualityConstraintEquation)]
+  ,InequalityProblem)
+
+-- | Given a distinct equation list, mutates each one at random using one of these mutations:
+-- * Unchanged
+-- * Generates similar but redundant equations
+-- * Generates its dual (to be transformed into an equality equation)
+-- * Generates an inconsistent partner (rare - 20% chance of existing in the returned problem).
+-- The equations passed in do not have to be consistent, merely unique and normalised.
+-- Returns the input, and the expected output.
+mutateEquations :: InequalityProblem -> Gen MutatedEquation
+mutateEquations ineq = do ineq' <- mapM mutate ineq >>*
+                            foldl (\(a,b,c) (x,y,z) -> (a++x,b++y,c++z)) ([],[],[])
+                          return ineq'
+  --TODO add the inconsistent option in as described in the documentation (and test for it)
+  where
+    mutate :: InequalityConstraintEquation -> Gen MutatedEquation
+    mutate ineq = oneof
+                    [
+                      return ([ineq],[],[ineq])
+                     ,addRedundant ineq
+                     ,return $ addDual ineq
+                    ]
+
+    addDual :: InequalityConstraintEquation -> MutatedEquation
+    addDual eq = ([eq,neg],[(eq,neg)],[]) where neg = amap negate eq
+
+    addRedundant :: InequalityConstraintEquation -> Gen MutatedEquation
+    addRedundant ineq = do i <- choose (1,5) -- number of redundant equations to add
+                           newIneqs <- replicateM i addRedundant'
+                           return (ineq : newIneqs, [], [ineq])
+                             where
+                               -- A redundant equation is one with a bigger unit coefficient:
+                               addRedundant' = do n <- choose (1,100)
+                                                  return $ ineq // [(0,n + (ineq ! 0))]
+
+newtype OmegaPruneInput = OPI MutatedEquation deriving (Show)
+
+instance Arbitrary OmegaPruneInput where
+  arbitrary = (generateProblem  >>= (return . snd) >>= mutateEquations) >>* OPI
+
+qcOmegaPrune :: [QuickCheckTest]
+qcOmegaPrune = [scaleQC (10,100,1000,10000) prop]
+  where
+    --We perform the map assocs because we can't compare arrays using *==*
+    -- (toConstr fails in the pretty-printing!).
+    prop (OPI (inp,outEq,outIneq)) =
+      (sort (map assocs (snd result)) *==* sort (map assocs outIneq))
+      *&&* (checkEq (fst result) outEq)
+      where
+        Just result = pruneAndCheck inp
+
+    checkEq :: [EqualityConstraintEquation] ->
+      [(EqualityConstraintEquation,EqualityConstraintEquation)] -> Result
+    checkEq [] [] = mkPassResult
+    checkEq eqs [] = mkFailResult $ "checkEq: " ++ show eqs
+    checkEq eqs ((x,y):xys)
+      = case findAndRemove (\z -> z == x || z == y) eqs of
+          (Just _, eqs') -> checkEq eqs' xys
+          _ -> mkFailResult $ "checkEq: " ++ show eqs ++ " could not match: " ++ show (x,y)
+
 qcTests :: (Test, [QuickCheckTest])
 qcTests = (TestList
  [
@@ -415,7 +480,7 @@ qcTests = (TestList
   ,testReachDef
   ,testArrayCheck
  ]
- ,qcOmegaEquality)
+ ,qcOmegaEquality ++ qcOmegaPrune)