Added inequality elimination to the Omega Test, but need to also add the code for checking for integer solutions

transformations/ArrayUsageCheck.hs

@@ -307,3 +307,73 @@ numVariables :: InequalityProblem -> Int
 numVariables ineq = length (nub $ concatMap findVars ineq)
   where
     findVars = map fst . filter ((/= 0) . snd) . tail . assocs
+
+-- | Eliminating the inequalities works as follows:
+--
+-- Rearrange (and normalise) equations for a particular variable x_k to eliminate so that 
+-- a_k is always positive and you have:
+--  A: a_Ak . x_k <= sum (i is 0 to n, without k) a_Ai . x_i
+--  B: a_Bk . x_k >= sum (i is 0 to n, without k) a_Bi . x_i
+--  C: equations where a_k is zero.
+--
+-- Determine if there is an integer solution for x_k:
+--
+-- TODO
+--
+-- Form lots of new equations:
+--  Given  a_Ak . x_k <= RHS(A)
+--         a_Bk . x_k >= RHS(B)
+--  We can get (since a_Ak and a_bk are positive):
+--         a_Ak . A_Bk . x_k <= A_Bk . RHS(A)
+--         a_Ak . A_Bk . x_k >= A_Ak . RHS(B)
+--  For every combination of the RHS(A) and RHS(B) generate an inequality: a_Bk . RHS(A) - a_Ak . RHS(B) >=0
+-- Add these new equations to the set C, and iterate
+
+fmElimination :: InequalityProblem -> InequalityProblem
+fmElimination ineq = fmElimination' (presentItems ineq) (map normaliseIneq ineq)
+  where
+    -- Finds all variables that have at least one non-zero coefficient in the equation set.
+    -- a_0 is ignored; 0 will never be in the returned list
+    presentItems :: InequalityProblem -> [Int]
+    presentItems = nub . map fst . filter ((/= 0) . snd) . concatMap (tail . assocs)
+
+    fmElimination' :: [Int] -> InequalityProblem -> InequalityProblem
+    fmElimination' [] ineq = ineq
+    fmElimination' (k:ks) ineq = fmElimination' ks (map normaliseIneq $ eliminate k ineq)
+    
+    --TODO should we still be checking for redundant equations in the new ones we generate?
+    
+    eliminate :: Int -> InequalityProblem -> InequalityProblem
+    eliminate k ineq = eqC ++ map (uncurry pairIneqs) (product2 (eqA,eqB))
+      where
+        (eqA,eqB,eqC) = partition ineq
+      
+        -- We need to partition the related equations into sets A,B and C.
+        -- C is straightforward (a_k is zero).
+        -- In set B, a_k > 0, and "RHS(B)" (as described above) is the negation of the other
+        -- coefficients.  Therefore "-RHS(B)" is the other coefficients as-is.
+        -- In set A, a_k < 0,  and "RHS(A)" (as described above) is the other coefficients, untouched
+        -- So we simply zero out a_k, and return the rest, associated with the absolute value of a_k.
+        partition :: InequalityProblem -> ([(Integer, InequalityConstraintEquation)], [(Integer,InequalityConstraintEquation)], [InequalityConstraintEquation])
+        partition = (\(x,y,z) -> (concat x, concat y, concat z)) . unzip3 . map partition'
+          where
+            partition' e | a_k == 0 = ([],[],[e])
+                         | a_k <  0 = ([(abs a_k, e // [(k,0)])],[],[])
+                         | a_k >  0 = ([],[(abs a_k, e // [(k,0)])],[])
+              where
+                a_k = e ! k
+        
+        pairIneqs :: (Integer, InequalityConstraintEquation) -> (Integer, InequalityConstraintEquation) -> InequalityConstraintEquation
+        pairIneqs (x,ex) (y,ey) = arrayZipWith (+) (amap (* y) ex) (amap (* x) ey)
+
+    normaliseIneq :: InequalityConstraintEquation -> InequalityConstraintEquation
+    normaliseIneq ineq | g > 1     = arrayMapWithIndex norm ineq
+                       | otherwise = ineq
+      where
+        norm ind val | ind == 0  = normaliseUnits val
+                     | otherwise = val `div` g
+      
+        g = mygcdList $ tail $ elems ineq
+        -- I think div would do here, because g will always be positive, but
+        -- I feel safer using the mathematical description:
+        normaliseUnits a_0 = floor $ (fromInteger a_0 :: Double) / (fromInteger g)

transformations/RainUsageCheckTest.hs

@@ -367,16 +367,17 @@ testIndexes = TestList
    -- TODO need to run the below exampls on a better test (they are not "easily" solved):
    
    -- Can (j + 1 % 10 == i + 1 % 10)?
-   ,notSolveable $ withKIsMod (i ++ con 1) 10 $ withNIsMod (j ++ con 1) 10 $ (4, [k === n], i_j_constraint 0 9)
+   ,neverSolveable $ withKIsMod (i ++ con 1) 10 $ withNIsMod (j ++ con 1) 10 $ (4, [k === n], i_j_constraint 0 9)
    -- Off by one (i + 1 % 9)
-   ,easilySolved $ withKIsMod (i ++ con 1) 9 $ withNIsMod (j ++ con 1) 9 $ (5, [k === n], i_j_constraint 0 9)
+   ,hardSolved $ withKIsMod (i ++ con 1) 9 $ withNIsMod (j ++ con 1) 9 $ (5, [k === n], i_j_constraint 0 9)
    
+   -- The "nightmare" example from the Omega Test paper:
+   ,neverSolveable (6,[],leq [con 27, 11 ** i ++ 13 ** j, con 45] &&& leq [con (-10), 7 ** i ++ (-9) ** j, con 4])
    
    ,safeParTest 100 True (0,10) [i]
    ,safeParTest 120 False (0,10) [i,i ++ con 1]
    ,safeParTest 140 True (0,10) [2 ** i, 2 ** i ++ con 1]
    
-   --TODO deal with modulo in future
   ]
   where
     -- Given some indexes using "i", this function checks whether these can
@@ -406,6 +407,8 @@ testIndexes = TestList
     i_j_constraint :: Integer -> Integer -> [HandyIneq]
     i_j_constraint low high = [con low <== i, i ++ con 1 <== j, j <== con high]
     
+    --TODO clear up this mess of easilySolved/hardSolved helper functions
+    
     findSolveable :: [(Int, [HandyEq], [HandyIneq])] -> [(Int, [HandyEq], [HandyIneq])]
     findSolveable = filter isSolveable
     
@@ -424,10 +427,21 @@ testIndexes = TestList
             -- we can't give a useful test failure here:
             else assertFailure $ "testIndexes " ++ show ind ++ " more than one variable left after solving"
 
+    hardSolved :: (Int, [HandyEq], [HandyIneq]) -> Test
+    hardSolved (ind, eq, ineq) = TestCase $
+      assertBool ("testIndexes " ++ show ind) $ isJust $ 
+        (uncurry solveAndPrune) (makeConsistent eq ineq) >>= (pruneAndCheck . fmElimination)
+
     notSolveable :: (Int, [HandyEq], [HandyIneq]) -> Test
     notSolveable (ind, eq, ineq) = TestCase $ assertEqual ("testIndexes " ++ show ind) Nothing $
       (uncurry solveAndPrune) (makeConsistent eq ineq) >>* ((<= 1) . numVariables)
 
+
+    neverSolveable :: (Int, [HandyEq], [HandyIneq]) -> Test
+    neverSolveable (ind, eq, ineq) = TestCase $ assertEqual ("testIndexes " ++ show ind) Nothing $
+      (uncurry solveAndPrune) (makeConsistent eq ineq) >>= (pruneAndCheck . fmElimination)
+
+
     -- The easy way of writing equations is built on the following Haskell magic.
     -- Essentially, everything is a list of (index, coefficient).  You can scale
     -- with the ** operator, and you can form equalities and inequalities with