Changeset 21

Show
Ignore:
Timestamp:
07/20/08 20:02:07 (2 years ago)
Author:
thesz
Message:

I wrote a Fourier-Motzkin implementation. It seem that there is a bug out there (prob1 doesn't fail), but I'm too tired to find it out right now.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • fort/IConstr.hs

    r20 r21  
    99import Control.Monad.State 
    1010 
    11 import Data.List (sortBy) 
     11import Data.List (sortBy,groupBy,nub) 
    1212import qualified Data.Map as Map 
    1313import Data.Maybe 
    1414import qualified Data.Set as Set 
    1515 
     16import Debug.Trace 
     17 
    1618------------------------------------------------------------------------------- 
    1719-- Helpers. 
    …… 
    2224sortBySnd list = sortBy (\(_,a) (_,b) -> compare a b) list 
    2325 
    24 type Coef = Rational 
     26type Coef = Integer 
    2527type MAP k a = Map.Map k a 
    2628 
    …… 
    3032data Vec v = 
    3133        Vec (MAP v Coef) Coef   -- a set of v*coef and constant. 
    32         deriving Show 
     34        deriving (Eq,Show) 
     35 
     36vec vcs co = foldr vecPlus (vecC co) $ map (uncurry vecCV) vcs 
    3337 
    3438vecC = Vec Map.empty 
    …… 
    6569vecDelVar v (Vec vcs c) = Vec (Map.delete v vcs) c 
    6670 
     71vecDivBy (Vec vcs co) gcd = Vec (Map.map (`div` gcd) vcs) (div co gcd) 
     72 
     73vecCoefsGCD (Vec vcs co) = case coefs of 
     74        [] -> abs co 
     75        xs -> foldr gcd (abs co) xs 
     76        where 
     77                coefs = Map.elems vcs 
     78 
    6779vecDef vec@(Vec vcs c) = case Map.toList vcs of 
    6880        [] -> Nothing 
    69         ((v,c):_) -> Just (v,vecScale (negate (1/c)) $ vecDelVar v vec) 
    70  
    71 vecSubstVar v def vec@(Vec vcs c) = vecMinus vec $ vecScale vc def 
     81        ((v,c):_) -> Just ((abs c,v),vecScale (negate $ signum c) $ vecReduceGCD $ vecDelVar v vec) 
     82 
     83vecDefVar v vec@(Vec vcs c) = case Map.lookup v vcs of 
     84        Nothing -> Nothing 
     85        Just c -> Just ((abs c,v),vecScale (negate $ signum c) $ vecReduceGCD $ vecDelVar v vec) 
     86 
     87vecReduceGCD vec 
     88        | gcd /= 0 = vecDivBy vec gcd 
     89        | otherwise = vec 
     90        where 
     91                gcd = vecCoefsGCD vec 
     92 
     93vecVarCoef v (Vec vcs _) = maybe 0 id (Map.lookup v vcs) 
     94 
     95vecSubstVar c v def vec@(Vec vcs _) = vecReduceGCD $ vecPlus (vecScale c $ vecDelVar v vec) (vecScale vc def) 
    7296        where 
    7397                vc = Map.findWithDefault 0 v vcs 
    …… 
    76100        ((v,c):_) -> (v,c) 
    77101        _ -> error "vecHeadVar: constant vector." 
     102 
     103vecVars (Vec vcs _) = Map.keys vcs 
     104 
     105vecHasVar v vec@(Vec vcs _) = fmap (\c -> (c,v,vecReduceGCD $ vecScale (negate $ signum c) $ vecDelVar v vec)) $ Map.lookup v vcs 
    78106 
    79107-------------------------------------------------------------------------------- 
    …… 
    112140                        HiBound x -> HiBound (x*c) 
    113141                        PInf -> PInf 
    114                  
     142 
     143-------------------------------------------------------------------------------- 
     144-- Inequalities. 
     145 
     146data Rel = LE | GE deriving (Eq, Ord, Show) 
    115147 
    116148-------------------------------------------------------------------------------- 
    …… 
    118150 
    119151data SolverS v = SolverS { 
    120                  ssDefs         :: [(v,Vec v)] 
     152                 ssDefs         :: [((Coef,v),Vec v)] 
     153                -- result of FM elim. 
     154                -- ((A,x),(u,v)) means that max u<=Ax<=min v 
     155                ,ssRanges       :: [((Coef,v),([Vec v],[Vec v]))] 
    121156                ,ssEqs          :: [Vec v] 
    122                 ,ssIneqs        :: [(Vec v,Range Coef)] 
     157                ,ssIneqs        :: [(Vec v,Rel)]       -- vec >= 0 or vec <= 0. 
    123158        } 
    124159        deriving Show 
    125160 
    126 emptySolverS = SolverS [] [] [] 
     161emptySolverS = SolverS [] [] [] [] 
    127162 
    128163type SolverM v a = StateT (SolverS v) [] a 
    129164 
    130 runSolver :: Ord v => SolverM v a -> [((),SolverS v)] 
     165runSolver :: (Show v, Ord v) => SolverM v a -> [((),SolverS v)] 
    131166runSolver problem = runStateT (problem >> solve) emptySolverS 
    132167 
    …… 
    143178-- Solution process. 
    144179 
    145 solve :: Ord v => SolverM v () 
     180solve :: (Show v, Ord v) => SolverM v () 
    146181solve = do 
    147182        checkEqs        -- we have to checkEqs once at start. 
    …… 
    149184        removeEqs 
    150185        -- then we ought to check and combine inequalities. 
    151         checkCombineInequalities 
     186        removeIneqVars 
    152187        -- and restore ranges. 
    153188        return () 
    154189 
    155 removeEqs :: Ord v => SolverM v () 
     190-------------------------------------------------------------------------------- 
     191-- Dealing with equalities. 
     192 
     193removeEqs :: (Show v, Ord v) => SolverM v () 
    156194removeEqs = do 
    157195        ss <- getSS 
    …… 
    159197                [] -> return () 
    160198                (eq:_) -> let 
    161                                 def = createDef eq 
     199                                def = trace ("creating def from "++show eq) $ createDef eq 
    162200                        in do 
    163                                 addDef def 
     201                                trace ("adding def "++show def) $ addDef def 
     202                                checkEqs 
    164203                                removeEqs 
    165204        where 
    …… 
    184223addDef def = do 
    185224        ss <- getSS 
    186         let ss' = ss { 
     225        let ss' = trace ("ss before removal: "++show ss) $ ss { 
    187226                         ssDefs = (def:) $ removeDef def snd (\(v,_) vec' -> (v,vec')) $ ssDefs ss 
    188227                        ,ssEqs = removeDef def id (flip const) $ ssEqs ss 
    189228                        ,ssIneqs = removeDef def fst (\(_,r) vec -> (vec,r)) $ ssIneqs ss 
    190229                } 
    191         putSS ss' 
    192  
    193 removeDef (v,defVec) get chg list = map (\x -> chg x $ vecSubstVar v defVec (get x)) list 
    194  
    195 checkCombineInequalities :: Ord v => SolverM v () 
    196 checkCombineInequalities = do 
    197         ss <- getSS 
    198         failMultiVar $ ssIneqs ss 
    199         let nIneqs = map normIneq $ ssIneqs ss 
    200         error "We should group normalized inequalities and tighten them and check contradiction." 
    201         putSS $ ss { ssIneqs = nIneqs } 
    202         where 
    203                 normIneq (vec,r) = (vecScale m vecMinusCo, rangeScale m rMinusCo) 
     230        trace ("putting ss after def removal: "++show ss') $ putSS ss' 
     231 
     232removeDef ((c,v),defVec) get chg list = map (\x -> chg x $ vecSubstVar c v defVec (get x)) list 
     233 
     234eq = vecPlus (vecMinus (vecV "i") (vecV "i1")) (vecV "n") 
     235tt = removeDef ((1,"i"),vecMinus (vecV "i1") (vecV "n")) id (flip const) [eq] 
     236 
     237-------------------------------------------------------------------------------- 
     238-- Dealing with inequalities (Fourier-Motzkin). 
     239 
     240removeIneqVars :: Ord v => SolverM v () 
     241removeIneqVars = do 
     242        checkIneqs 
     243        ss <- getSS 
     244        case selectVar $ ssIneqs ss of 
     245                Just ((c,v),(lows,highs),ineqs) -> do 
     246                        let add = [ (vecMinus l h,LE) | l <- lows, h <- highs] 
     247                        let rng = ((c,v),(lows,highs)) 
     248                        putSS $ ss { 
     249                                         ssRanges = rng : ssRanges ss 
     250                                        ,ssIneqs = ineqs ++ add 
     251                                } 
     252                        removeIneqVars 
     253                Nothing -> return () 
     254 
     255checkIneqs = do 
     256        ss <- getSS 
     257        let ies = ssIneqs ss 
     258        ies' <- mapM checkIE ies 
     259        putSS $ ss { ssIneqs = concat ies' } 
     260        where 
     261                checkIE ie@(vec,rel) 
     262                        | Just c <- vecIsConst vec = 
     263                                if contradict c rel then failSS else return [] 
     264                        | otherwise = return [ie] 
     265                contradict c LE = c > 0 
     266                contradict c GE = c < 0 
     267 
     268selectVar ineqs = case allGroups of 
     269        [] -> Nothing 
     270        (r:_) -> Just r 
     271        where 
     272                allVars = nub $ concatMap (vecVars . fst) ineqs 
     273                allGroups = filter nonEmpty $ map group' allVars 
     274                nonEmpty (_,(h,l),_) = not (null h) && not (null l) 
     275                group' var = group [] [] [] [] var ineqs 
     276                group coefs lows highs ineqs var [] = ((m,var),(map scl lows,map scl highs),ineqs) 
    204277                        where 
    205                                 co = vecConst vec 
    206                                 vecMinusCo = vecMinus vec (vecC co) 
    207                                 rMinusCo = rangeSubConst r co 
    208                                 (v,c) = vecHeadVar vec 
    209                                 m = signum c 
    210                 failMultiVar [] = return () 
    211                 failMultiVar ((v,r):vrs) 
    212                         | vecNumVars v > 1 = failSS 
    213                         | otherwise = failMultiVar vrs 
     278                                m = foldr lcm 1 coefs 
     279                                scl (c,def) = vecScale (div m c) def 
     280                group coefs lows highs ineqs var (ie@(vec,rel):ies) 
     281                        | Just (c,v,def) <- vecHasVar var vec = let 
     282                                        (c',lows',highs') = selectRel c def 
     283                                in group (c':coefs) lows' highs' ineqs var ies 
     284                        | otherwise = group coefs lows highs (ie:ineqs) var ies 
     285                        where 
     286                                negateRel c rel 
     287                                        | c < 0 = (abs c,invRel) 
     288                                        | otherwise = (c,rel) 
     289                                        where 
     290                                                invRel = case rel of 
     291                                                        LE -> GE 
     292                                                        GE -> LE 
     293                                selectRel c def = case rel' of 
     294                                        -- cx<=0+def. Adding def to highs. 
     295                                        LE -> (c',lows,(c,def):highs) 
     296                                        -- vice versa. 
     297                                        GE -> (c',(c,def):lows,highs) 
     298                                        where 
     299                                                (c',rel') = negateRel c rel 
     300 
     301tv1 = Vec (Map.fromList [("i",1),("i1",2)]) 4 
     302tv2 = Vec (Map.fromList [("i",-1),("i1",2)]) 4 
     303 
     304ties1 = [ 
     305                 (vec [(1,"i")] (-1),GE) 
     306                ,(vecMinus (vecV "i") (vecV "n"), LE) 
     307        ] 
     308 
     309-- SolverS {ssDefs = [((1,"i"),Vec (fromList [("i1",1),("n",-1)]) 0)], ssRanges = [], ssEqs = [Vec (fromList []) 0], ssIneqs = [(Vec (fromList [("n",-1)]) 1,LE),(Vec (fromList [("i1",1),("n",-1)]) 0,LE),(Vec (fromList [("i1",-1)]) 1,LE),(Vec (fromList [("i1",1),("n",-2)]) 0,LE),(Vec (fromList [("i1",-1),("n",1)]) 1,LE)]} 
     310ties2 = [(Vec (Map.fromList [("i1",1),("n",-1)]) 0,LE),(Vec (Map.fromList [("i1",-1)]) 1,LE),(Vec (Map.fromList [("i1",1),("n",-2)]) 0,LE),(Vec (Map.fromList [("i1",-1),("n",1)]) 1,LE)] 
     311-------------------------------------------------------------------------------- 
     312-- Building a problem. 
     313 
     314insertIE a rel b = do 
     315        ss <- getSS 
     316        putSS $ ss { ssIneqs = (vecMinus a b,rel):ssIneqs ss } 
     317 
     318insertEQ a b = do 
     319        ss <- getSS 
     320        putSS $ ss { ssEqs = vecMinus a b:ssEqs ss } 
     321 
     322between x lo hi = do 
     323        insertIE lo LE x 
     324        insertIE x  LE hi 
     325 
     326betweens xs lo hi = do 
     327        mapM (\x -> between x lo hi) xs 
     328        return () 
     329 
     330-- check solutions for (do i=1,n { a[i+n] = a[i]}) problem. 
     331tprob1 = do 
     332        let [i,i1,n] = map vecV ["i","i1","n"] 
     333        let one = vecC 1 
     334        betweens [i,i1] one n 
     335--      insertIE (vecPlus one i) LE i1 
     336        insertEQ (vecPlus i n) i1 
    214337 
    215338--------------------------------------------------------------------------------