-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Matrix.Simplex
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Two-step simplex algorithm
--
-- I only guarantee that this module wastes inodes
--
-----------------------------------------------------------------------------

-- Originally based off the code in Sedgewick, but modified to match the
-- conventions from Papadimitriou and Steiglitz.

-- TODO: Is our column/row selection the same as Bland's anti-cycle
-- algorithm?

-- TODO: Add check for redundant rows in two-phase algorithm

-- TODO: Lots of testing

module Matrix.Simplex (Simplex(..), simplex, twophase) where

import Data.Array

eps :: Double
eps :: Double
eps = Double
1.0e-10

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

-- Pivot around a!(p,q)

pivot :: Int -> Int -> Array (Int,Int) Double -> Array (Int,Int) Double
pivot :: Int -> Int -> Array (Int, Int) Double -> Array (Int, Int) Double
pivot Int
p Int
q Array (Int, Int) Double
a0 = forall {e}. Num e => Array (Int, Int) e -> Array (Int, Int) e
step4 forall a b. (a -> b) -> a -> b
$ forall {e}.
Fractional e =>
Array (Int, Int) e -> Array (Int, Int) e
step3 forall a b. (a -> b) -> a -> b
$ forall {e}. Num e => Array (Int, Int) e -> Array (Int, Int) e
step2 forall a b. (a -> b) -> a -> b
$ forall {e}.
Fractional e =>
Array (Int, Int) e -> Array (Int, Int) e
step1 forall a b. (a -> b) -> a -> b
$ Array (Int, Int) Double
a0
    where step1 :: Array (Int, Int) e -> Array (Int, Int) e
step1 Array (Int, Int) e
a = Array (Int, Int) e
a forall i e. Ix i => Array i e -> [(i, e)] -> Array i e
// [ ((Int
j,Int
k), Array (Int, Int) e
aforall i e. Ix i => Array i e -> i -> e
!(Int
j,Int
k) forall a. Num a => a -> a -> a
- Array (Int, Int) e
aforall i e. Ix i => Array i e -> i -> e
!(Int
p,Int
k) forall a. Num a => a -> a -> a
* Array (Int, Int) e
aforall i e. Ix i => Array i e -> i -> e
!(Int
j,Int
q) forall a. Fractional a => a -> a -> a
/ Array (Int, Int) e
aforall i e. Ix i => Array i e -> i -> e
!(Int
p,Int
q)) | Int
k <- [Int
0..Int
m], Int
j <- [Int
ph..Int
n], Int
j forall a. Eq a => a -> a -> Bool
/= Int
p Bool -> Bool -> Bool
&& Int
k forall a. Eq a => a -> a -> Bool
/= Int
q ]
          step2 :: Array (Int, Int) e -> Array (Int, Int) e
step2 Array (Int, Int) e
a = Array (Int, Int) e
a forall i e. Ix i => Array i e -> [(i, e)] -> Array i e
// [ ((Int
j,Int
q),e
0) | Int
j <- [Int
ph..Int
n], Int
j forall a. Eq a => a -> a -> Bool
/= Int
p ]
          step3 :: Array (Int, Int) e -> Array (Int, Int) e
step3 Array (Int, Int) e
a = Array (Int, Int) e
a forall i e. Ix i => Array i e -> [(i, e)] -> Array i e
// [ ((Int
p,Int
k), Array (Int, Int) e
aforall i e. Ix i => Array i e -> i -> e
!(Int
p,Int
k) forall a. Fractional a => a -> a -> a
/ Array (Int, Int) e
aforall i e. Ix i => Array i e -> i -> e
!(Int
p,Int
q)) | Int
k <- [Int
0..Int
m], Int
k forall a. Eq a => a -> a -> Bool
/= Int
q ]
          step4 :: Array (Int, Int) e -> Array (Int, Int) e
step4 Array (Int, Int) e
a = Array (Int, Int) e
a forall i e. Ix i => Array i e -> [(i, e)] -> Array i e
// [ ((Int
p,Int
q),e
1) ]
          ((Int
ph,Int
_),(Int
n,Int
m)) = forall i e. Array i e -> (i, i)
bounds Array (Int, Int) Double
a0

-- chooseq picks the lowest numbered favorable column.  If there are no
-- favorable columns, then q==m is returned, and we have reached an
-- optimum.


chooseq :: (Ord b, Num b, Ix a, Ix b, Num a) =>
           Array (a, b) Double -> b
chooseq :: forall b a.
(Ord b, Num b, Ix a, Ix b, Num a) =>
Array (a, b) Double -> b
chooseq Array (a, b) Double
a0 = forall {a}. (Ix a, Num a) => b -> Array (a, b) Double -> b
chooseq' b
1 Array (a, b) Double
a0
    where chooseq' :: b -> Array (a, b) Double -> b
chooseq' b
q Array (a, b) Double
a | b
q forall a. Ord a => a -> a -> Bool
> b
m          = b
q
                       | Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
0,b
q) forall a. Ord a => a -> a -> Bool
< -Double
eps = b
q
                       | Bool
otherwise      = b -> Array (a, b) Double -> b
chooseq' (b
qforall a. Num a => a -> a -> a
+b
1) Array (a, b) Double
a
          ((a
_,b
_),(a
_,b
m)) = forall i e. Array i e -> (i, i)
bounds Array (a, b) Double
a0

-- choosep picks a row with a positive element in column q.  If no such
-- element exists, then the p==n is returned, and the problem is
-- unfeasible.

choosep :: (Ord a, Num a, Ix a, Ix b) =>
           b -> Array (a, b) Double -> a
choosep :: forall a b.
(Ord a, Num a, Ix a, Ix b) =>
b -> Array (a, b) Double -> a
choosep b
q Array (a, b) Double
a0 = a -> Array (a, b) Double -> a
choosep' a
1 Array (a, b) Double
a0
    where choosep' :: a -> Array (a, b) Double -> a
choosep' a
p Array (a, b) Double
a | a
p forall a. Ord a => a -> a -> Bool
> a
n         = a
p
                       | Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
p,b
q) forall a. Ord a => a -> a -> Bool
> Double
eps = a
p
                       | Bool
otherwise     = a -> Array (a, b) Double -> a
choosep' (a
pforall a. Num a => a -> a -> a
+a
1) Array (a, b) Double
a
          ((a
_,b
_),(a
n,b
_)) = forall i e. Array i e -> (i, i)
bounds Array (a, b) Double
a0

-- refinep picks the row using the ratio test.

refinep :: (Ord a, Num a, Num b, Ix a, Ix b) =>
           a -> b -> Array (a, b) Double -> a
refinep :: forall a b.
(Ord a, Num a, Num b, Ix a, Ix b) =>
a -> b -> Array (a, b) Double -> a
refinep a
p0 b
q Array (a, b) Double
a0 = a -> a -> Array (a, b) Double -> a
refinep' (a
p0forall a. Num a => a -> a -> a
+a
1) a
p0 Array (a, b) Double
a0
    where refinep' :: a -> a -> Array (a, b) Double -> a
refinep' a
i a
p Array (a, b) Double
a | a
i forall a. Ord a => a -> a -> Bool
> a
n = a
p
                         | Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
i,b
q) forall a. Ord a => a -> a -> Bool
> Double
eps Bool -> Bool -> Bool
&& Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
i,b
0) forall a. Fractional a => a -> a -> a
/ Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
i,b
q) forall a. Ord a => a -> a -> Bool
< Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
p,b
0) forall a. Fractional a => a -> a -> a
/ Array (a, b) Double
aforall i e. Ix i => Array i e -> i -> e
!(a
p,b
q) = a -> a -> Array (a, b) Double -> a
refinep' (a
iforall a. Num a => a -> a -> a
+a
1) a
i Array (a, b) Double
a
                         | Bool
otherwise = a -> a -> Array (a, b) Double -> a
refinep' (a
iforall a. Num a => a -> a -> a
+a
1) a
p Array (a, b) Double
a
          ((a
_,b
_),(a
n,b
_)) = forall i e. Array i e -> (i, i)
bounds Array (a, b) Double
a0

-- * Types

-- | Type for results of the simplex algorithm

data Simplex a = Unbounded | Infeasible | Optimal a deriving (ReadPrec [Simplex a]
ReadPrec (Simplex a)
ReadS [Simplex a]
forall a. Read a => ReadPrec [Simplex a]
forall a. Read a => ReadPrec (Simplex a)
forall a. Read a => Int -> ReadS (Simplex a)
forall a. Read a => ReadS [Simplex a]
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Simplex a]
$creadListPrec :: forall a. Read a => ReadPrec [Simplex a]
readPrec :: ReadPrec (Simplex a)
$creadPrec :: forall a. Read a => ReadPrec (Simplex a)
readList :: ReadS [Simplex a]
$creadList :: forall a. Read a => ReadS [Simplex a]
readsPrec :: Int -> ReadS (Simplex a)
$creadsPrec :: forall a. Read a => Int -> ReadS (Simplex a)
Read,Int -> Simplex a -> ShowS
forall a. Show a => Int -> Simplex a -> ShowS
forall a. Show a => [Simplex a] -> ShowS
forall a. Show a => Simplex a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Simplex a] -> ShowS
$cshowList :: forall a. Show a => [Simplex a] -> ShowS
show :: Simplex a -> String
$cshow :: forall a. Show a => Simplex a -> String
showsPrec :: Int -> Simplex a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Simplex a -> ShowS
Show)

-- * Functions

-- | The simplex algorithm for standard form:
--
-- min   c'x
--
-- where Ax = b, x >= 0
--
-- a!(0,0) = -z
--
-- a!(0,j) = c'
--
-- a!(i,0) = b
--
-- a!(i,j) = A_ij

simplex :: Array (Int,Int) Double -- ^ stating tableau
        -> Simplex (Array (Int,Int) Double) -- ^ solution

simplex :: Array (Int, Int) Double -> Simplex (Array (Int, Int) Double)
simplex Array (Int, Int) Double
a | Int
q forall a. Ord a => a -> a -> Bool
> Int
m      = forall a. a -> Simplex a
Optimal Array (Int, Int) Double
a
          | Int
p forall a. Ord a => a -> a -> Bool
> Int
n      = forall a. Simplex a
Unbounded
          | Bool
otherwise  = Array (Int, Int) Double -> Simplex (Array (Int, Int) Double)
simplex forall a b. (a -> b) -> a -> b
$ Int -> Int -> Array (Int, Int) Double -> Array (Int, Int) Double
pivot Int
p' Int
q forall a b. (a -> b) -> a -> b
$ Array (Int, Int) Double
a
    where q :: Int
q = forall b a.
(Ord b, Num b, Ix a, Ix b, Num a) =>
Array (a, b) Double -> b
chooseq Array (Int, Int) Double
a
          p :: Int
p = forall a b.
(Ord a, Num a, Ix a, Ix b) =>
b -> Array (a, b) Double -> a
choosep Int
q Array (Int, Int) Double
a
          p' :: Int
p' = forall a b.
(Ord a, Num a, Num b, Ix a, Ix b) =>
a -> b -> Array (a, b) Double -> a
refinep Int
p Int
q Array (Int, Int) Double
a
          ((Int
_,Int
_),(Int
n,Int
m)) = forall i e. Array i e -> (i, i)
bounds Array (Int, Int) Double
a

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

addart :: (Num e, Enum a, Ix a, Num a) =>
          Array (a, a) e -> Array (a, a) e
addart :: forall e a.
(Num e, Enum a, Ix a, Num a) =>
Array (a, a) e -> Array (a, a) e
addart Array (a, a) e
a = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array ((-a
1,a
0),(a
n,a
mforall a. Num a => a -> a -> a
+a
n)) forall a b. (a -> b) -> a -> b
$ [((a, a), e)]
z forall a. [a] -> [a] -> [a]
++ [((a, a), e)]
xsi forall a. [a] -> [a] -> [a]
++ [((a, a), e)]
b forall a. [a] -> [a] -> [a]
++ [((a, a), e)]
art forall a. [a] -> [a] -> [a]
++ [((a, a), e)]
x
    where z :: [((a, a), e)]
z = ((-a
1,a
0), Array (a, a) e
aforall i e. Ix i => Array i e -> i -> e
!(a
0,a
0)) forall a. a -> [a] -> [a]
: [ ((-a
1,a
j),e
0) | a
j <- [a
1..a
n] ] forall a. [a] -> [a] -> [a]
++ [ ((-a
1,a
jforall a. Num a => a -> a -> a
+a
n),Array (a, a) e
aforall i e. Ix i => Array i e -> i -> e
!(a
0,a
j)) | a
j <- [a
1..a
m] ]
          xsi :: [((a, a), e)]
xsi = ((a
0,a
0), -forall e a b.
(Num e, Num a, Enum a, Ix a, Ix b) =>
Array (a, b) e -> b -> e
colsum Array (a, a) e
a a
0) forall a. a -> [a] -> [a]
: [ ((a
0,a
j),e
0) | a
j <- [a
1..a
n] ] forall a. [a] -> [a] -> [a]
++ [ ((a
0,a
jforall a. Num a => a -> a -> a
+a
n), -forall e a b.
(Num e, Num a, Enum a, Ix a, Ix b) =>
Array (a, b) e -> b -> e
colsum Array (a, a) e
a a
j) | a
j <- [a
1..a
m] ]
          b :: [((a, a), e)]
b = [ ((a
i,a
0), Array (a, a) e
aforall i e. Ix i => Array i e -> i -> e
!(a
i,a
0)) | a
i <- [a
1..a
n] ]
          art :: [((a, a), e)]
art = [ ((a
i,a
j), if a
i forall a. Eq a => a -> a -> Bool
== a
j then e
1 else e
0) | a
i <- [a
1..a
n], a
j <- [a
1..a
n] ]
          x :: [((a, a), e)]
x = [ ((a
i,a
jforall a. Num a => a -> a -> a
+a
n), Array (a, a) e
aforall i e. Ix i => Array i e -> i -> e
!(a
i,a
j)) | a
i <- [a
1..a
n], a
j <- [a
1..a
m] ]
          ((a
_,a
_),(a
n,a
m)) = forall i e. Array i e -> (i, i)
bounds Array (a, a) e
a

colsum :: (Num e, Num a, Enum a, Ix a, Ix b) =>
          Array (a, b) e -> b -> e
colsum :: forall e a b.
(Num e, Num a, Enum a, Ix a, Ix b) =>
Array (a, b) e -> b -> e
colsum Array (a, b) e
a b
j = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array (a, b) e
aforall i e. Ix i => Array i e -> i -> e
!(a
i,b
j) | a
i <- [a
1..a
n] ]
    where ((a
_,b
_),(a
n,b
_)) = forall i e. Array i e -> (i, i)
bounds Array (a, b) e
a

delart :: (Enum a, Ix a, Num a) =>
          Array (a, a) e -> Array (a, a) e -> Array (a, a) e
delart :: forall a e.
(Enum a, Ix a, Num a) =>
Array (a, a) e -> Array (a, a) e -> Array (a, a) e
delart Array (a, a) e
a Array (a, a) e
a' = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array ((a
0,a
0),(a
n,a
m)) forall a b. (a -> b) -> a -> b
$ [((a, a), e)]
z forall a. [a] -> [a] -> [a]
++ [((a, a), e)]
b forall a. [a] -> [a] -> [a]
++ [((a, a), e)]
x
    where z :: [((a, a), e)]
z = ((a
0,a
0), Array (a, a) e
a'forall i e. Ix i => Array i e -> i -> e
!(-a
1,a
0)) forall a. a -> [a] -> [a]
: [ ((a
0,a
j), Array (a, a) e
aforall i e. Ix i => Array i e -> i -> e
!(a
0,a
j)) | a
j <- [a
1..a
m] ]
          b :: [((a, a), e)]
b = [ ((a
i,a
0), Array (a, a) e
a'forall i e. Ix i => Array i e -> i -> e
!(a
i,a
0)) | a
i <- [a
1..a
n] ]
          x :: [((a, a), e)]
x = [ ((a
i,a
j), Array (a, a) e
a'forall i e. Ix i => Array i e -> i -> e
!(a
i,a
jforall a. Num a => a -> a -> a
+a
n)) | a
i <- [a
1..a
n], a
j <- [a
1..a
m] ]
          ((a
_,a
_),(a
n,a
m)) = forall i e. Array i e -> (i, i)
bounds Array (a, a) e
a

-- | The two-phase simplex algorithm

twophase :: Array (Int,Int) Double -- ^ stating tableau
         -> Simplex (Array (Int,Int) Double) -- ^ solution

twophase :: Array (Int, Int) Double -> Simplex (Array (Int, Int) Double)
twophase Array (Int, Int) Double
a | forall e a b.
(Num e, Ix a, Ix b, Num a, Num b) =>
Simplex (Array (a, b) e) -> e
cost Simplex (Array (Int, Int) Double)
a' forall a. Ord a => a -> a -> Bool
> Double
eps = forall a. Simplex a
Infeasible
           | Bool
otherwise     = Array (Int, Int) Double -> Simplex (Array (Int, Int) Double)
simplex forall a b. (a -> b) -> a -> b
$ forall a e.
(Enum a, Ix a, Num a) =>
Array (a, a) e -> Array (a, a) e -> Array (a, a) e
delart Array (Int, Int) Double
a (forall a. Simplex a -> a
gettab Simplex (Array (Int, Int) Double)
a')
    where a' :: Simplex (Array (Int, Int) Double)
a' = Array (Int, Int) Double -> Simplex (Array (Int, Int) Double)
simplex forall a b. (a -> b) -> a -> b
$ forall e a.
(Num e, Enum a, Ix a, Num a) =>
Array (a, a) e -> Array (a, a) e
addart forall a b. (a -> b) -> a -> b
$ Array (Int, Int) Double
a


-- How to handle cases where 'simplex' does not return Optimal?
gettab :: Simplex a -> a
gettab :: forall a. Simplex a -> a
gettab (Optimal a
a) = a
a

cost :: (Num e, Ix a, Ix b, Num a, Num b) =>
        Simplex (Array (a, b) e) -> e
cost :: forall e a b.
(Num e, Ix a, Ix b, Num a, Num b) =>
Simplex (Array (a, b) e) -> e
cost (Optimal Array (a, b) e
a) = forall a. Num a => a -> a
negate forall a b. (a -> b) -> a -> b
$ Array (a, b) e
aforall i e. Ix i => Array i e -> i -> e
!(a
0,b
0)

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

{-

Test vectors

This is from Sedgewick

> x1 = listArray ((0,0),(5,8)) [  0, -1, -1, -1, 0, 0, 0, 0, 0,
>                                 5, -1,  1,  0, 1, 0, 0, 0, 0,
>                                45,  1,  4,  0, 0, 1, 0, 0, 0,
>                                27,  2,  1,  0, 0, 0, 1, 0, 0,
>                                24,  3, -4,  0, 0, 0, 0, 1, 0,
>                                 4,  0,  0,  1, 0, 0, 0, 0, 1 ] :: Array (Int,Int) Double

P&S, Example 2.6

> x2 = listArray ((0,0),(3,5)) [ 0, 1, 1, 1, 1, 1,
>                                1, 3, 2, 1, 0, 0,
>                                3, 5, 1, 1, 1, 0,
>                                4, 2, 5, 1, 0, 1 ] :: Array (Int,Int) Double

P&S, Example 2.6 (after BFS selection)

> x2' = listArray ((0,0),(3,5)) [ -6, -3, -3,  0,  0,  0,
>                                 1,  3,  2,  1,  0,  0,
>                                 2,  2, -1,  0,  1,  0,
>                                 3, -1,  3,  0,  0,  1 ] :: Array (Int,Int) Double

P&S, Example 2.2 / Section 2.9

> x3 = listArray ((0,0),(4,7)) [ -34, -1, -14, -6, 0, 0, 0, 0,
>                                  4,  1,   1,  1, 1, 0, 0, 0,
>                                  2,  1,   0,  0, 0, 1, 0, 0,
>                                  3,  0,   0,  1, 0, 0, 1, 0,
>                                  6,  0,   3,  1, 0, 0, 0, 1 ] :: Array (Int,Int) Double

P&S, Example 2.7

> x4 = listArray ((0,0),(3,7)) [ 3, -3/4,  20, -1/2, 6, 0, 0, 0,
>                                0,  1/4,  -8,   -1, 9, 1, 0, 0,
>                                0,  1/2, -12, -1/2, 3, 0, 1, 0,
>                                1,    0,   0,    1, 0, 0, 0, 1 ] :: Array (Int,Int) Double

These come in handy for testing

> row j a = listArray (0,m) [ a!(j,k) | k <- [0..m] ]
>    where ((_,_),(n,m)) = bounds a

> column k a = listArray (0,n) [ a!(j,k) | j <- [0..n] ]
>    where ((_,_),(n,m)) = bounds a

> solution (Optimal a) = listArray (1,m) $ [ find a j | j <- [1..m] ]
>    where ((_,_),(n,m)) = bounds a

> find a j = findone' a 1 j
>     where findone' a i j | i > n          = 0
>                          | a!(i,j) == 1.0 = b!i
>                          | otherwise      = findone' a (i+1) j
>           b = listArray (1,n) [ a!(i,0) | i <- [1..n] ]
>           ((_,_),(n,m)) = bounds a

-}