msakai / toysolver / 496

Committed 10 Nov 2024 11:05AM UTC coverage: 69.994% (-1.1%) from 71.113%

Build # 496

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #117 from msakai/update-coveralls-and-haddock

GitHub Actions: Update coveralls and haddock configuration

Run Details

9872 of 14104 relevant lines covered (69.99%)

0.7 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

65.22

/src/ToySolver/Arith/Simplex/Textbook.hs

{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.Simplex.Textbook
-- Copyright   :  (c) Masahiro Sakai 2011
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Naïve implementation of Simplex method
--
-- Reference:
--
-- * <http://www.math.cuhk.edu.hk/~wei/lpch3.pdf>
--
-----------------------------------------------------------------------------
module ToySolver.Arith.Simplex.Textbook
  (
  -- * The @Tableau@ type
    Tableau
  , RowIndex
  , ColIndex
  , Row
  , emptyTableau
  , objRowIndex
  , pivot
  , lookupRow
  , addRow
  , setObjFun

  -- * Optimization direction
  , module Data.OptDir

  -- * Reading status
  , currentValue
  , currentObjValue
  , isFeasible
  , isOptimal

  -- * High-level solving functions
  , simplex
  , dualSimplex
  , phaseI
  , primalDualSimplex

  -- * For debugging
  , isValidTableau
  , toCSV
  ) where

import Data.Ord
import Data.List
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.OptDir
import Data.VectorSpace
import Control.Exception

import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar

-- ---------------------------------------------------------------------------

type Tableau r = VarMap (Row r)
{-
tbl ! v == (m, val)
==>
var v .+. m .==. constant val
-}

-- | Basic variables
type RowIndex = Int

-- | Non-basic variables
type ColIndex = Int

type Row r = (VarMap r, r)

data PivotResult r = PivotUnbounded | PivotFinished | PivotSuccess (Tableau r)
  deriving (Show, Eq, Ord)

emptyTableau :: Tableau r
emptyTableau = IM.empty

objRowIndex :: RowIndex
objRowIndex = -1

pivot :: (Fractional r, Eq r) => RowIndex -> ColIndex -> Tableau r -> Tableau r
{-# INLINE pivot #-}
{-# SPECIALIZE pivot :: RowIndex -> ColIndex -> Tableau Rational -> Tableau Rational #-}
{-# SPECIALIZE pivot :: RowIndex -> ColIndex -> Tableau Double -> Tableau Double #-}
pivot r s tbl =
    assert (isValidTableau tbl) $  -- precondition
    assert (isValidTableau tbl') $ -- postcondition
      tbl'
  where
    tbl' = IM.insert s row_s $ IM.map f $ IM.delete r $ tbl
    f orig@(row_i, row_i_val) =
      case IM.lookup s row_i of
        Nothing -> orig
        Just c ->
          ( IM.filter (0/=) $ IM.unionWith (+) (IM.delete s row_i) (IM.map (negate c *) row_r')
          , row_i_val - c*row_r_val'
          )
    (row_r, row_r_val) = lookupRow r tbl
    a_rs = row_r IM.! s
    row_r' = IM.map (/ a_rs) $ IM.insert r 1 $ IM.delete s row_r
    row_r_val' = row_r_val / a_rs
    row_s = (row_r', row_r_val')

-- | Lookup a row by basic variable
lookupRow :: RowIndex -> Tableau r -> Row r
lookupRow r m = m IM.! r

-- 行の基底変数の列が0になるように変形
normalizeRow :: (Num r, Eq r) => Tableau r -> Row r -> Row r
normalizeRow a (row0,val0) = obj'
  where
    obj' = g $ foldl' f (IM.empty, val0) $
           [ case IM.lookup j a of
               Nothing -> (IM.singleton j x, 0)
               Just (row,val) -> (IM.map ((-x)*) (IM.delete j row), -x*val)
           | (j,x) <- IM.toList row0 ]
    f (m1,v1) (m2,v2) = (IM.unionWith (+) m1 m2, v1+v2)
    g (m,v) = (IM.filter (0/=) m, v)

setRow :: (Num r, Eq r) => Tableau r -> RowIndex -> Row r -> Tableau r
setRow tbl i row = assert (isValidTableau tbl) $ assert (isValidTableau tbl') $ tbl'
  where
    tbl' = IM.insert i (normalizeRow tbl row) tbl

addRow :: (Num r, Eq r) => Tableau r -> RowIndex -> Row r -> Tableau r
addRow tbl i row = assert (i `IM.notMember` tbl) $ setRow tbl i row

setObjFun :: (Num r, Eq r) => Tableau r -> LA.Expr r -> Tableau r
setObjFun tbl e = addRow tbl objRowIndex row
  where
    row =
      case LA.extract LA.unitVar e of
        (c, e') -> (LA.coeffMap (negateV e'), c)

copyObjRow :: (Num r, Eq r) => Tableau r -> Tableau r -> Tableau r
copyObjRow from to =
  case IM.lookup objRowIndex from of
    Nothing -> IM.delete objRowIndex to
    Just row -> addRow to objRowIndex row

currentValue :: Num r => Tableau r -> Var -> r
currentValue tbl v =
  case IM.lookup v tbl of
    Nothing -> 0
    Just (_, val) -> val

currentObjValue :: Tableau r -> r
currentObjValue = snd . lookupRow objRowIndex

isValidTableau :: Tableau r -> Bool
isValidTableau tbl =
  and [v `IM.notMember` m | (v, (m,_)) <- IM.toList tbl, v /= objRowIndex] &&
  and [IS.null (IM.keysSet m `IS.intersection` vs) | (m,_) <- IM.elems tbl']
  where
    tbl' = IM.delete objRowIndex tbl
    vs = IM.keysSet tbl'

isFeasible :: Real r => Tableau r -> Bool
isFeasible tbl =
  and [val >= 0 | (v, (_,val)) <- IM.toList tbl, v /= objRowIndex]

isOptimal :: Real r => OptDir -> Tableau r -> Bool
isOptimal optdir tbl =
  and [not (cmp cj) | cj <- IM.elems (fst (lookupRow objRowIndex tbl))]
  where
    cmp = case optdir of
            OptMin -> (0<)
            OptMax -> (0>)

isImproving :: Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptMin from to = currentObjValue to <= currentObjValue from
isImproving OptMax from to = currentObjValue to >= currentObjValue from

-- ---------------------------------------------------------------------------
-- primal simplex

simplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
{-# SPECIALIZE simplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE simplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
simplex optdir = go
  where
    go tbl = assert (isFeasible tbl) $
      case primalPivot optdir tbl of
        PivotFinished  -> assert (isOptimal optdir tbl) (True, tbl)
        PivotUnbounded -> (False, tbl)
        PivotSuccess tbl' -> assert (isImproving optdir tbl tbl') $ go tbl'

primalPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> PivotResult r
{-# INLINE primalPivot #-}
primalPivot optdir tbl
  | null cs   = PivotFinished
  | null rs   = PivotUnbounded
  | otherwise = PivotSuccess (pivot r s tbl)
  where
    cmp = case optdir of
            OptMin -> (0<)
            OptMax -> (0>)
    cs = [(j,cj) | (j,cj) <- IM.toList (fst (lookupRow objRowIndex tbl)), cmp cj]
    -- smallest subscript rule
    s = fst $ head cs
    -- classical rule
    --s = fst $ (if optdir==OptMin then maximumBy else minimumBy) (compare `on` snd) cs
    rs = [ (i, y_i0 / y_is)
         | (i, (row_i, y_i0)) <- IM.toList tbl, i /= objRowIndex
         , let y_is = IM.findWithDefault 0 s row_i, y_is > 0
         ]
    r = fst $ minimumBy (comparing snd) rs

-- ---------------------------------------------------------------------------
-- dual simplex

dualSimplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
{-# SPECIALIZE dualSimplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE dualSimplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
dualSimplex optdir = go
  where
    go tbl = assert (isOptimal optdir tbl) $
      case dualPivot optdir tbl of
        PivotFinished  -> assert (isFeasible tbl) $ (True, tbl)
        PivotUnbounded -> (False, tbl)
        PivotSuccess tbl' -> assert (isImproving optdir tbl' tbl) $ go tbl'

dualPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> PivotResult r
{-# INLINE dualPivot #-}
dualPivot optdir tbl
  | null rs   = PivotFinished
  | null cs   = PivotUnbounded
  | otherwise = PivotSuccess (pivot r s tbl)
  where
    rs = [(i, row_i) | (i, (row_i, y_i0)) <- IM.toList tbl, i /= objRowIndex, 0 > y_i0]
    (r, row_r) = head rs
    cs = [ (j, if optdir==OptMin then y_0j / y_rj else - y_0j / y_rj)
         | (j, y_rj) <- IM.toList row_r
         , y_rj < 0
         , let y_0j = IM.findWithDefault 0 j obj
         ]
    s = fst $ minimumBy (comparing snd) cs
    (obj,_) = lookupRow objRowIndex tbl

-- ---------------------------------------------------------------------------
-- phase I of the two-phased method

phaseI :: (Real r, Fractional r) => Tableau r -> VarSet -> (Bool, Tableau r)
{-# SPECIALIZE phaseI :: Tableau Rational -> VarSet -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE phaseI :: Tableau Double -> VarSet -> (Bool, Tableau Double) #-}
phaseI tbl avs
  | currentObjValue tbl1' /= 0 = (False, tbl1')
  | otherwise = (True, copyObjRow tbl $ removeArtificialVariables avs $ tbl1')
  where
    optdir = OptMax
    tbl1 = setObjFun tbl $ negateV $ sumV [LA.var v | v <- IS.toList avs]
    tbl1' = go tbl1
    go tbl2
      | currentObjValue tbl2 == 0 = tbl2
      | otherwise =
        case primalPivot optdir tbl2 of
          PivotSuccess tbl2' -> assert (isImproving optdir tbl2 tbl2') $ go tbl2'
          PivotFinished -> assert (isOptimal optdir tbl2) tbl2
          PivotUnbounded -> error "phaseI: should not happen"

-- post-processing of phaseI
-- pivotを使ってartificial variablesを基底から除いて、削除
removeArtificialVariables :: (Real r, Fractional r) => VarSet -> Tableau r -> Tableau r
removeArtificialVariables avs tbl0 = tbl2
  where
    tbl1 = foldl' f (IM.delete objRowIndex tbl0) (IS.toList avs)
    tbl2 = IM.map (\(row,val) ->  (IM.filterWithKey (\j _ -> j `IS.notMember` avs) row, val)) tbl1
    f tbl i =
      case IM.lookup i tbl of
        Nothing -> tbl
        Just row ->
          case [j | (j,c) <- IM.toList (fst row), c /= 0, j `IS.notMember` avs] of
            [] -> IM.delete i tbl
            j:_ -> pivot i j tbl

-- ---------------------------------------------------------------------------
-- primal-dual simplex

data PDResult = PDUnsat | PDOptimal | PDUnbounded

primalDualSimplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
{-# SPECIALIZE primalDualSimplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE primalDualSimplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
primalDualSimplex optdir = go
  where
    go tbl =
      case pdPivot optdir tbl of
        Left PDOptimal   -> assert (isFeasible tbl) $ assert (isOptimal optdir tbl) $ (True, tbl)
        Left PDUnsat     -> assert (not (isFeasible tbl)) $ (False, tbl)
        Left PDUnbounded -> assert (not (isOptimal optdir tbl)) $ (False, tbl)
        Right tbl' -> go tbl'

pdPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> Either PDResult (Tableau r)
pdPivot optdir tbl
  | null ps && null qs = Left PDOptimal -- optimal
  | otherwise =
      case ret of
        Left p -> -- primal update
          let rs = [ (i, (bi - t) / y_ip)
                   | (i, (row_i, bi)) <- IM.toList tbl, i /= objRowIndex
                   , let y_ip = IM.findWithDefault 0 p row_i, y_ip > 0
                   ]
              q = fst $ minimumBy (comparing snd) rs
          in if null rs
             then Left PDUnsat
             else Right (pivot q p tbl)
        Right q -> -- dual update
          let (row_q, _bq) = (tbl IM.! q)
              cs = [ (j, (cj'-t) / (-y_qj))
                   | (j, y_qj) <- IM.toList row_q
                   , y_qj < 0
                   , let cj  = IM.findWithDefault 0 j obj
                   , let cj' = if optdir==OptMax then cj else -cj
                   ]
              p = fst $ maximumBy (comparing snd) cs
              (obj,_) = lookupRow objRowIndex tbl
          in if null cs
             then Left PDUnbounded -- dual infeasible
             else Right (pivot q p tbl)
  where
    qs = [ (Right i, bi) | (i, (_row_i, bi)) <- IM.toList tbl, i /= objRowIndex, 0 > bi ]
    ps = [ (Left j, cj')
         | (j,cj) <- IM.toList (fst (lookupRow objRowIndex tbl))
         , let cj' = if optdir==OptMax then cj else -cj
         , 0 > cj' ]
    (ret, t) = minimumBy (comparing snd) (qs ++ ps)

-- ---------------------------------------------------------------------------

toCSV :: (Num r) => (r -> String) -> Tableau r -> String
toCSV showCell tbl = unlines . map (concat . intersperse ",") $ header : body
  where
    header :: [String]
    header = "" : map colName cols ++ [""]

    body :: [[String]]
    body = [showRow i (lookupRow i tbl) | i <- rows]

    rows :: [RowIndex]
    rows = IM.keys (IM.delete objRowIndex tbl) ++ [objRowIndex]

    cols :: [ColIndex]
    cols = [0..colMax]
      where
        colMax = maximum (-1 : [c | (row, _) <- IM.elems tbl, c <- IM.keys row])

    rowName :: RowIndex -> String
    rowName i = if i==objRowIndex then "obj" else "x" ++ show i

    colName :: ColIndex -> String
    colName j = "x" ++ show j

    showRow i (row, row_val) = rowName i : [showCell (IM.findWithDefault 0 j row') | j <- cols] ++ [showCell row_val]
      where row' = IM.insert i 1 row

-- ---------------------------------------------------------------------------

1	{-# LANGUAGE ScopedTypeVariables #-}
2	{-# OPTIONS_GHC -Wall #-}
3	{-# OPTIONS_HADDOCK show-extensions #-}
4	-----------------------------------------------------------------------------
5	-- \|
6	-- Module : ToySolver.Arith.Simplex.Textbook
7	-- Copyright : (c) Masahiro Sakai 2011
8	-- License : BSD-style
9	--
10	-- Maintainer : masahiro.sakai@gmail.com
11	-- Stability : provisional
12	-- Portability : non-portable
13	--
14	-- Naïve implementation of Simplex method
15	--
16	-- Reference:
17	--
18	-- * <http://www.math.cuhk.edu.hk/~wei/lpch3.pdf>
19	--
20	-----------------------------------------------------------------------------
21	module ToySolver.Arith.Simplex.Textbook
22	(
23	-- * The @Tableau@ type
24	Tableau
25	, RowIndex
26	, ColIndex
27	, Row
28	, emptyTableau
29	, objRowIndex
30	, pivot
31	, lookupRow
32	, addRow
33	, setObjFun
34
35	-- * Optimization direction
36	, module Data.OptDir
37
38	-- * Reading status
39	, currentValue
40	, currentObjValue
41	, isFeasible
42	, isOptimal
43
44	-- * High-level solving functions
45	, simplex
46	, dualSimplex
47	, phaseI
48	, primalDualSimplex
49
50	-- * For debugging
51	, isValidTableau
52	, toCSV
53	) where
54
55	import Data.Ord
56	import Data.List
57	import qualified Data.IntMap as IM
58	import qualified Data.IntSet as IS
59	import Data.OptDir
60	import Data.VectorSpace
61	import Control.Exception
62
63	import qualified ToySolver.Data.LA as LA
64	import ToySolver.Data.IntVar
65
66	-- ---------------------------------------------------------------------------
67
68	type Tableau r = VarMap (Row r)
69	{-
70	tbl ! v == (m, val)
71	==>
72	var v .+. m .==. constant val
73	-}
74
75	-- \| Basic variables
76	type RowIndex = Int
77
78	-- \| Non-basic variables
79	type ColIndex = Int
80
81	type Row r = (VarMap r, r)
82
83	data PivotResult r = PivotUnbounded \| PivotFinished \| PivotSuccess (Tableau r)
84	deriving (Show, Eq, Ord)	×
85
86	emptyTableau :: Tableau r
87	emptyTableau = IM.empty	1✔
88
89	objRowIndex :: RowIndex
90	objRowIndex = -1	1✔
91
92	pivot :: (Fractional r, Eq r) => RowIndex -> ColIndex -> Tableau r -> Tableau r
93	{-# INLINE pivot #-}
94	{-# SPECIALIZE pivot :: RowIndex -> ColIndex -> Tableau Rational -> Tableau Rational #-}
95	{-# SPECIALIZE pivot :: RowIndex -> ColIndex -> Tableau Double -> Tableau Double #-}
96	pivot r s tbl =	1✔
97	assert (isValidTableau tbl) $ -- precondition	×
98	assert (isValidTableau tbl') $ -- postcondition	×
99	tbl'	1✔
100	where
101	tbl' = IM.insert s row_s $ IM.map f $ IM.delete r $ tbl	1✔
102	f orig@(row_i, row_i_val) =	1✔
103	case IM.lookup s row_i of	1✔
104	Nothing -> orig	1✔
105	Just c ->
106	( IM.filter (0/=) $ IM.unionWith (+) (IM.delete s row_i) (IM.map (negate c *) row_r')	1✔
107	, row_i_val - c*row_r_val'	1✔
108	)
109	(row_r, row_r_val) = lookupRow r tbl	1✔
110	a_rs = row_r IM.! s	1✔
111	row_r' = IM.map (/ a_rs) $ IM.insert r 1 $ IM.delete s row_r	1✔
112	row_r_val' = row_r_val / a_rs	1✔
113	row_s = (row_r', row_r_val')	1✔
114
115	-- \| Lookup a row by basic variable
116	lookupRow :: RowIndex -> Tableau r -> Row r
117	lookupRow r m = m IM.! r	1✔
118
119	-- 行の基底変数の列が0になるように変形
120	normalizeRow :: (Num r, Eq r) => Tableau r -> Row r -> Row r
121	normalizeRow a (row0,val0) = obj'	1✔
122	where
123	obj' = g $ foldl' f (IM.empty, val0) $	1✔
124	[ case IM.lookup j a of	1✔
125	Nothing -> (IM.singleton j x, 0)	1✔
126	Just (row,val) -> (IM.map ((-x)) (IM.delete j row), -xval)	1✔
127	\| (j,x) <- IM.toList row0 ]	1✔
128	f (m1,v1) (m2,v2) = (IM.unionWith (+) m1 m2, v1+v2)	1✔
129	g (m,v) = (IM.filter (0/=) m, v)	1✔
130
131	setRow :: (Num r, Eq r) => Tableau r -> RowIndex -> Row r -> Tableau r
132	setRow tbl i row = assert (isValidTableau tbl) $ assert (isValidTableau tbl') $ tbl'	×
133	where
134	tbl' = IM.insert i (normalizeRow tbl row) tbl	1✔
135
136	addRow :: (Num r, Eq r) => Tableau r -> RowIndex -> Row r -> Tableau r
137	addRow tbl i row = assert (i `IM.notMember` tbl) $ setRow tbl i row	×
138
139	setObjFun :: (Num r, Eq r) => Tableau r -> LA.Expr r -> Tableau r
140	setObjFun tbl e = addRow tbl objRowIndex row	1✔
141	where
142	row =	1✔
143	case LA.extract LA.unitVar e of	1✔
144	(c, e') -> (LA.coeffMap (negateV e'), c)	1✔
145
146	copyObjRow :: (Num r, Eq r) => Tableau r -> Tableau r -> Tableau r
147	copyObjRow from to =	1✔
148	case IM.lookup objRowIndex from of	1✔
149	Nothing -> IM.delete objRowIndex to	1✔
150	Just row -> addRow to objRowIndex row	×
151
152	currentValue :: Num r => Tableau r -> Var -> r
153	currentValue tbl v =	1✔
154	case IM.lookup v tbl of	1✔
155	Nothing -> 0	1✔
156	Just (_, val) -> val	1✔
157
158	currentObjValue :: Tableau r -> r
159	currentObjValue = snd . lookupRow objRowIndex	1✔
160
161	isValidTableau :: Tableau r -> Bool
162	isValidTableau tbl =	1✔
163	and [v `IM.notMember` m \| (v, (m,_)) <- IM.toList tbl, v /= objRowIndex] &&	1✔
164	and [IS.null (IM.keysSet m `IS.intersection` vs) \| (m,_) <- IM.elems tbl']	1✔
165	where
166	tbl' = IM.delete objRowIndex tbl	1✔
167	vs = IM.keysSet tbl'	1✔
168
169	isFeasible :: Real r => Tableau r -> Bool
170	isFeasible tbl =	1✔
171	and [val >= 0 \| (v, (_,val)) <- IM.toList tbl, v /= objRowIndex]	1✔
172
173	isOptimal :: Real r => OptDir -> Tableau r -> Bool
174	isOptimal optdir tbl =	1✔
175	and [not (cmp cj) \| cj <- IM.elems (fst (lookupRow objRowIndex tbl))]	1✔
176	where
177	cmp = case optdir of	1✔
178	OptMin -> (0<)	×
179	OptMax -> (0>)	1✔
180
181	isImproving :: Real r => OptDir -> Tableau r -> Tableau r -> Bool
182	isImproving OptMin from to = currentObjValue to <= currentObjValue from	×
183	isImproving OptMax from to = currentObjValue to >= currentObjValue from	×
184
185	-- ---------------------------------------------------------------------------
186	-- primal simplex
187
188	simplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
189	{-# SPECIALIZE simplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
190	{-# SPECIALIZE simplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
191	simplex optdir = go	1✔
192	where
193	go tbl = assert (isFeasible tbl) $	×
194	case primalPivot optdir tbl of	1✔
195	PivotFinished -> assert (isOptimal optdir tbl) (True, tbl)	×
196	PivotUnbounded -> (False, tbl)	×
197	PivotSuccess tbl' -> assert (isImproving optdir tbl tbl') $ go tbl'	×
198
199	primalPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> PivotResult r
200	{-# INLINE primalPivot #-}
201	primalPivot optdir tbl	1✔
202	\| null cs = PivotFinished	1✔
203	\| null rs = PivotUnbounded	×
204	\| otherwise = PivotSuccess (pivot r s tbl)	×
205	where
206	cmp = case optdir of	1✔
207	OptMin -> (0<)	1✔
208	OptMax -> (0>)	1✔
209	cs = [(j,cj) \| (j,cj) <- IM.toList (fst (lookupRow objRowIndex tbl)), cmp cj]	×
210	-- smallest subscript rule
211	s = fst $ head cs	1✔
212	-- classical rule
213	--s = fst $ (if optdir==OptMin then maximumBy else minimumBy) (compare `on` snd) cs
214	rs = [ (i, y_i0 / y_is)	1✔
215	\| (i, (row_i, y_i0)) <- IM.toList tbl, i /= objRowIndex	1✔
216	, let y_is = IM.findWithDefault 0 s row_i, y_is > 0	×
217	]
218	r = fst $ minimumBy (comparing snd) rs	1✔
219
220	-- ---------------------------------------------------------------------------
221	-- dual simplex
222
223	dualSimplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
224	{-# SPECIALIZE dualSimplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
225	{-# SPECIALIZE dualSimplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
226	dualSimplex optdir = go	1✔
227	where
228	go tbl = assert (isOptimal optdir tbl) $	×
229	case dualPivot optdir tbl of	1✔
230	PivotFinished -> assert (isFeasible tbl) $ (True, tbl)	×
231	PivotUnbounded -> (False, tbl)	×
232	PivotSuccess tbl' -> assert (isImproving optdir tbl' tbl) $ go tbl'	×
233
234	dualPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> PivotResult r
235	{-# INLINE dualPivot #-}
236	dualPivot optdir tbl	1✔
237	\| null rs = PivotFinished	1✔
238	\| null cs = PivotUnbounded	×
239	\| otherwise = PivotSuccess (pivot r s tbl)	×
240	where
241	rs = [(i, row_i) \| (i, (row_i, y_i0)) <- IM.toList tbl, i /= objRowIndex, 0 > y_i0]	1✔
242	(r, row_r) = head rs	1✔
243	cs = [ (j, if optdir==OptMin then y_0j / y_rj else - y_0j / y_rj)	×
244	\| (j, y_rj) <- IM.toList row_r	1✔
245	, y_rj < 0	1✔
246	, let y_0j = IM.findWithDefault 0 j obj	×
247	]
248	s = fst $ minimumBy (comparing snd) cs	1✔
249	(obj,_) = lookupRow objRowIndex tbl	1✔
250
251	-- ---------------------------------------------------------------------------
252	-- phase I of the two-phased method
253
254	phaseI :: (Real r, Fractional r) => Tableau r -> VarSet -> (Bool, Tableau r)
255	{-# SPECIALIZE phaseI :: Tableau Rational -> VarSet -> (Bool, Tableau Rational) #-}
256	{-# SPECIALIZE phaseI :: Tableau Double -> VarSet -> (Bool, Tableau Double) #-}
257	phaseI tbl avs	1✔
258	\| currentObjValue tbl1' /= 0 = (False, tbl1')	×
259	\| otherwise = (True, copyObjRow tbl $ removeArtificialVariables avs $ tbl1')	×
260	where
261	optdir = OptMax	1✔
262	tbl1 = setObjFun tbl $ negateV $ sumV [LA.var v \| v <- IS.toList avs]	1✔
263	tbl1' = go tbl1	1✔
264	go tbl2	1✔
265	\| currentObjValue tbl2 == 0 = tbl2	1✔
266	\| otherwise =	×
267	case primalPivot optdir tbl2 of	1✔
268	PivotSuccess tbl2' -> assert (isImproving optdir tbl2 tbl2') $ go tbl2'	×
269	PivotFinished -> assert (isOptimal optdir tbl2) tbl2	×
270	PivotUnbounded -> error "phaseI: should not happen"	×
271
272	-- post-processing of phaseI
273	-- pivotを使ってartificial variablesを基底から除いて、削除
274	removeArtificialVariables :: (Real r, Fractional r) => VarSet -> Tableau r -> Tableau r
275	removeArtificialVariables avs tbl0 = tbl2	1✔
276	where
277	tbl1 = foldl' f (IM.delete objRowIndex tbl0) (IS.toList avs)	1✔
278	tbl2 = IM.map (\(row,val) -> (IM.filterWithKey (\j _ -> j `IS.notMember` avs) row, val)) tbl1	1✔
279	f tbl i =	1✔
280	case IM.lookup i tbl of	1✔
281	Nothing -> tbl	1✔
282	Just row ->
283	case [j \| (j,c) <- IM.toList (fst row), c /= 0, j `IS.notMember` avs] of	×
284	[] -> IM.delete i tbl	×
285	j:_ -> pivot i j tbl	×
286
287	-- ---------------------------------------------------------------------------
288	-- primal-dual simplex
289
290	data PDResult = PDUnsat \| PDOptimal \| PDUnbounded
291
292	primalDualSimplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
293	{-# SPECIALIZE primalDualSimplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
294	{-# SPECIALIZE primalDualSimplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
295	primalDualSimplex optdir = go	1✔
296	where
297	go tbl =	1✔
298	case pdPivot optdir tbl of	1✔
299	Left PDOptimal -> assert (isFeasible tbl) $ assert (isOptimal optdir tbl) $ (True, tbl)	×
300	Left PDUnsat -> assert (not (isFeasible tbl)) $ (False, tbl)	×
301	Left PDUnbounded -> assert (not (isOptimal optdir tbl)) $ (False, tbl)	×
302	Right tbl' -> go tbl'	1✔
303
304	pdPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> Either PDResult (Tableau r)
305	pdPivot optdir tbl	1✔
306	\| null ps && null qs = Left PDOptimal -- optimal	1✔
307	\| otherwise =	×
308	case ret of	1✔
309	Left p -> -- primal update
310	let rs = [ (i, (bi - t) / y_ip)	1✔
311	\| (i, (row_i, bi)) <- IM.toList tbl, i /= objRowIndex	1✔
312	, let y_ip = IM.findWithDefault 0 p row_i, y_ip > 0	1✔
313	]
314	q = fst $ minimumBy (comparing snd) rs	1✔
315	in if null rs	×
316	then Left PDUnsat	×
317	else Right (pivot q p tbl)	1✔
318	Right q -> -- dual update
319	let (row_q, _bq) = (tbl IM.! q)	1✔
320	cs = [ (j, (cj'-t) / (-y_qj))	1✔
321	\| (j, y_qj) <- IM.toList row_q	1✔
322	, y_qj < 0	×
323	, let cj = IM.findWithDefault 0 j obj	×
324	, let cj' = if optdir==OptMax then cj else -cj	×
325	]
326	p = fst $ maximumBy (comparing snd) cs	1✔
327	(obj,_) = lookupRow objRowIndex tbl	1✔
328	in if null cs	×
329	then Left PDUnbounded -- dual infeasible	×
330	else Right (pivot q p tbl)	1✔
331	where
332	qs = [ (Right i, bi) \| (i, (_row_i, bi)) <- IM.toList tbl, i /= objRowIndex, 0 > bi ]	1✔
333	ps = [ (Left j, cj')	1✔
334	\| (j,cj) <- IM.toList (fst (lookupRow objRowIndex tbl))	1✔
335	, let cj' = if optdir==OptMax then cj else -cj	×
336	, 0 > cj' ]	1✔
337	(ret, t) = minimumBy (comparing snd) (qs ++ ps)	1✔
338
339	-- ---------------------------------------------------------------------------
340
341	toCSV :: (Num r) => (r -> String) -> Tableau r -> String
342	toCSV showCell tbl = unlines . map (concat . intersperse ",") $ header : body	×
343	where
344	header :: [String]
345	header = "" : map colName cols ++ [""]	×
346
347	body :: [[String]]
348	body = [showRow i (lookupRow i tbl) \| i <- rows]	×
349
350	rows :: [RowIndex]
351	rows = IM.keys (IM.delete objRowIndex tbl) ++ [objRowIndex]	×
352
353	cols :: [ColIndex]
354	cols = [0..colMax]	×
355	where
356	colMax = maximum (-1 : [c \| (row, _) <- IM.elems tbl, c <- IM.keys row])	×
357
358	rowName :: RowIndex -> String
359	rowName i = if i==objRowIndex then "obj" else "x" ++ show i	×
360
361	colName :: ColIndex -> String
362	colName j = "x" ++ show j	×
363
364	showRow i (row, row_val) = rowName i : [showCell (IM.findWithDefault 0 j row') \| j <- cols] ++ [showCell row_val]	×
365	where row' = IM.insert i 1 row	×
366
367	-- ---------------------------------------------------------------------------

msakai / toysolver / 496

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous