work on Subtropical

src/Theory/NonLinearRealArithmatic/Expr.hs

@@ -32,6 +32,30 @@ varsIn expr = case expr of
   BinaryOp _ sub_expr1 sub_expr2 ->
     nubOrd (varsIn sub_expr1 <> varsIn sub_expr2)
 
+substitute :: Var -> Expr a -> Expr a -> Expr a
+substitute var expr_subst_with expr_subst_in = 
+  case expr_subst_in of
+    Const _ -> 
+      expr_subst_in
+    Var var' -> 
+      if var == var' then expr_subst_with else expr_subst_in
+    UnaryOp op sub_expr -> 
+      UnaryOp op (substitute var expr_subst_with sub_expr)
+    BinaryOp op sub_expr1 sub_expr2 -> 
+      BinaryOp op 
+        (substitute var expr_subst_with sub_expr1)
+        (substitute var expr_subst_with sub_expr2)
+
+instance Num a => Num (Expr a) where
+  expr1 + expr2 = BinaryOp Add expr1 expr2
+  expr1 * expr2 = BinaryOp Mul expr1 expr2
+  expr1 - expr2 = BinaryOp Sub expr1 expr2
+
+  abs expr = error "TODO: abs not defined for expr. What would make sense here?"
+  signum expr = error "TODO: signum not defined for expr. What would make sense here?"
+
+  fromInteger n = Const (fromInteger n)
+
 {-
 
 domainOf :: Bounded a => Var -> VarDomains a -> Interval a

src/Theory/NonLinearRealArithmatic/Polynomial.hs

@@ -13,7 +13,10 @@ module Theory.NonLinearRealArithmatic.Polynomial
   , isLinear
   , extractTerm
   , fromExpr
+  , toExpr
   , map
+  , degree
+  , termDegree
   , Assignable(..)
   , Assignment
   ) where
@@ -193,11 +196,11 @@ fromExpr expr =
     BinaryOp Div _ _ -> error "Division in user provided expressions not supported"
 
 toExpr :: forall a. (Ord a, Num a) => Polynomial a -> Expr a
-toExpr = foldr1 (BinaryOp Add) . fmap from_term . getTerms
+toExpr = sum . fmap from_term . getTerms
   where
     from_term :: Term a -> Expr a
     from_term (Term coeff monomial) = 
-      foldr (BinaryOp Mul) (Const coeff) (M.mapWithKey from_var monomial)
+      Const coeff * product (M.mapWithKey from_var monomial)
 
     from_var :: Var -> Int -> Expr a
     from_var var n 

src/Theory/NonLinearRealArithmatic/Subtropical.hs

@@ -1,14 +1,15 @@
 {-|
   Subtropical is a fast but incomplete method for solving non-linear real constraints.
 -}
-module Theory.NonLinearRealArithmatic.Subtropical 
-  ( subtropical 
-  ) where
+module Theory.NonLinearRealArithmatic.Subtropical (subtropical) where
 
 import qualified Data.List as L
 import qualified Data.IntMap as M
-import Data.IntMap ( IntMap )
+import qualified Data.IntMap.Merge.Lazy as MM
+import Theory.NonLinearRealArithmatic.Constraint (Constraint, ConstraintRelation (..))
+import Theory.NonLinearRealArithmatic.Expr (Expr (..), Var, BinaryOp (..), substitute)
 import Theory.NonLinearRealArithmatic.Polynomial
+import Control.Arrow ((&&&))
 
 -- |
 -- The frame of a polynomial is a set of points, obtained from the 
@@ -33,34 +34,34 @@ import Theory.NonLinearRealArithmatic.Polynomial
 frame :: (Ord a, Num a) => Polynomial a -> ([Monomial], [Monomial])
 frame polynomial = undefined -- TODO
   where
-    (pos_terms, neg_terms) 
-      = L.partition ((> 0) . getCoeff) 
+    (pos_terms, neg_terms)
+      = L.partition ((> 0) . getCoeff)
       $ L.filter ((/= 0) . getCoeff) (getTerms polynomial)
 
-findDominatingDirection :: (Num a, Ord a) => Polynomial a -> Maybe (IntMap Int)
+findDominatingDirection :: (Num a, Ord a) => Polynomial a -> Maybe (Assignment Int)
 findDominatingDirection terms = undefined
   where
     pos_terms = filter ((> 0) . getCoeff) (getTerms terms)
-    
+
 -- |
 -- Returns True iff the polynomial evaluates to something positive under 
 -- the given variable assignment.
-isPositiveSolution :: (Num a, Ord a, Assignable a) => Polynomial a -> IntMap a -> Bool
+isPositiveSolution :: (Num a, Ord a, Assignable a) => Polynomial a -> Assignment a -> Bool
 isPositiveSolution polynomial solution = eval solution polynomial > 0
 
 -- |
-positiveSolution :: (Num a, Ord a, Assignable a) => Polynomial a -> Maybe (IntMap a)
-positiveSolution polynomial = do 
+positiveSolution :: (Num a, Ord a, Assignable a) => Polynomial a -> Maybe (Assignment a)
+positiveSolution polynomial = do
   normal_vector <- findDominatingDirection polynomial
-  
+
   -- For a sufficiently large base the polynomial should evaluate 
   -- to something positive.
   let bases = [ 2^n | n <- [1..] ]
-  let candidates = [ M.map (b^) normal_vector | b <- bases ]      
+  let candidates = [ M.map (b^) normal_vector | b <- bases ]
 
   L.find (isPositiveSolution polynomial) candidates
-  
--- newtype Solution a = Sol { getValues :: IntMap a }
+
+-- newtype Solution a = Sol { getValues :: Assignment a }
 
 -- instance Num a => Num (Solution a) where
 --   (Sol s1) + (Sol s2) = Sol $ M.unionWith (+) s1 s2
@@ -69,37 +70,93 @@ positiveSolution polynomial = do
 --   negate (Sol s1) = Sol $ M.map negate s1
 --   signum (Sol s1) = Sol $ M.map signum s1
 --   fromInteger i = error "TODO: define this"
+
+-- | Returns True if the polynomial evaluates to 0 under the given variable assignment.
+isRoot :: (Num a, Ord a, Assignable a) => Polynomial a -> Assignment a -> Bool
+isRoot polynomial solution = eval solution polynomial == 0
+
+{-|
+  Given a positive and a negative solution, we can find a root to any multivariate polynomial.
+  By the intermediate value theorem, a root lies somewhere on the line segment, connecting the 
+  two points. For example, consider the polynomial
+
+    x^2 + y^2 - 3
+
+  The points (1,1) and (2,2) give a negative- and positive solution respectively. Intuitively, 
+  the roots of the polynomial form a circle of radius 3. The negative solution correspoinds to 
+  a point in the interior of the circle and the positive solution to a point on the exterior.
+  So the line segment connecting the two points:
+
+    (1,1) + t*((2,2) - (1,1)) = (1+t, 1+t)               (where 0 <= t <= 1)
   
--- |
--- Returns True if the polynomial evaluates to 0 under the given 
--- variable assignment.
-isSolution :: (Num a, Ord a, Assignable a) => Polynomial a -> IntMap a -> Bool
-isSolution polynomial solution = eval solution polynomial == 0
-
+  must intersect the circle. We find the intersection points by substituting `x = 1+2t` / `y = 1+2t`
+  and solving the original polynomial, thereby reducing the problem to finding roots of a univariate 
+  polynomial in the interval (0 :..: 1).
+-}
+f :: forall a. (Num a, Ord a, Assignable a, Fractional a, Floating a) => Polynomial a -> Assignment a -> Assignment a -> Assignment a
+f polynomial neg_sol pos_sol =
+  let
+    -- An arbitrary ID for the variable t. The polynomial we construct is univariate in t, 
+    -- so there is no danger of variable ID collision.
+    t = 0
+
+    line_segment_component :: Var -> a -> a -> Expr a
+    line_segment_component _ neg_sol_component pos_sol_component =
+      Const neg_sol_component + Var t * Const (pos_sol_component - neg_sol_component)
+
+    line_segment_components :: Assignment (Expr a)
+    line_segment_components = MM.merge MM.dropMissing MM.dropMissing (MM.zipWithMatched line_segment_component) neg_sol pos_sol
+
+    substitute_all :: Assignment (Expr a) -> Expr a -> Expr a
+    substitute_all assignment expr = M.foldrWithKey substitute expr assignment
+
+    polynomial_in_t = fromExpr $ substitute_all line_segment_components $ toExpr polynomial
+
+    degree_coeff_pairs :: [(Int, a)]
+    degree_coeff_pairs = L.sortOn fst $ (termDegree &&& getCoeff) <$> getTerms polynomial_in_t
+
+    roots_in_t :: [a]
+    roots_in_t =
+      case degree_coeff_pairs of
+        -- linear polynomial ==> solve directly for t
+        [ (0, c), (1, b) ] -> [ - b/c ]
+        -- quadratic polynomial ==> apply quadratic equation
+        [ (0, c), (1, b), (2, a) ] -> [ (-b + sqrt (b^2 - 4*a*c)) / (2*a), (-b - sqrt (b^2 - 4*a*c)) / (2*a) ]
+        -- TODO:
+        _ -> error "TODO: what to do with higher degree polynomials?"
+
+    t_solution :: Assignment a
+    t_solution = 
+      case L.find (\t' -> 0 <= t' && t' <= 1) roots_in_t of
+        Just value -> M.singleton t value
+        Nothing -> error "No solution for t between 0 and 1"
+  in
+    M.map (eval t_solution . fromExpr) line_segment_components
+
 {-| 
+  TODO: deal with to multiple constraints
 -}
-subtropical :: forall a. (Ord a, Assignable a) => Polynomial a -> Maybe (Assignment a)
-subtropical polynomial = 
+subtropical :: forall a. (Ord a, Assignable a) => Constraint a -> Maybe (Assignment a)
+subtropical (Equals, polynomial) =
   let
-    -- solution that maps all variables to one
+    -- assignment that maps all variables to one
     one :: Assignment a
     one = M.fromList $ zip (varsIn polynomial) (repeat 1)
-    
+
     go :: Assignment a -> Polynomial a -> Maybe (Assignment a)
     go neg_sol polynomial = do
       pos_sol <- positiveSolution polynomial
-
-      -- TODO: solve for t element [0;1]
-      -- neg_sol + t * (pos_sol - neg_sol)
+
+      -- TODO: solve for t element [0;1] neg_sol + t * (pos_sol - neg_sol)
       let t = 1
-      
-      return 
-        $ M.unionWith (+) neg_sol 
-        $ M.map (* t) 
+
+      return
+        $ M.unionWith (+) neg_sol
+        $ M.map (* t)
         $ M.unionWith (-) pos_sol neg_sol
   in
     case eval one polynomial `compare` 0 of
       LT -> go one polynomial
       GT -> go one (negate polynomial)
       EQ -> Just one
-      
+subtropical _ = error "TODO: implement subtropical for inequality constraints"