Added LPSolve as a solver

4 files changed, 204 insertions, 0 deletions

diff --git a/Math/LinProg/Compile.hs b/Math/LinProg/Compile.hs
new file mode 100644
index 0000000..07398e7
--- /dev/null
+++ b/Math/LinProg/Compile.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE TemplateHaskell, FlexibleInstances, ScopedTypeVariables #-}
+
+module Math.LinProg.Compile (
+  compile
+  ,Equation
+  ,CompilerS(..)
+  ,objective
+  ,equals
+  ,leqs
+) where
+
+import Data.List
+import Math.LinProg.Types
+import Control.Lens
+import Control.Monad.State
+import Control.Monad.Free
+import Data.Maybe
+
+type Equation t v = (LinExpr t v, t) -- LHS and RHS
+
+data CompilerS t v = CompilerS {
+  _objective :: LinExpr t v
+  ,_equals :: [Equation t v]
+  ,_leqs :: [Equation t v]
+} deriving (Eq)
+
+$(makeLenses ''CompilerS)
+
+compile :: (Num t, Show t, Ord t, Eq v) => LinProg t v () -> CompilerS t v
+compile ast = compile' ast initCompilerS where
+  compile' (Free (Objective a c)) state = compile' c $ state & objective +~ a
+  compile' (Free (EqConstraint a b c)) state = compile' c $ state & equals %~ (split (a-b):)
+  compile' (Free (LeqConstraint a b c)) state = compile' c $ state & leqs %~ (split (a-b):)
+  compile' _ state = state
+
+  initCompilerS = CompilerS
+    0
+    []
+    []
+
diff --git a/Math/LinProg/LPSolve.hs b/Math/LinProg/LPSolve.hs
new file mode 100644
index 0000000..23c3949
--- /dev/null
+++ b/Math/LinProg/LPSolve.hs
@@ -0,0 +1,67 @@
+{-# LANGUAGE ViewPatterns #-}
+
+module Math.LinProg.LPSolve where
+
+import Control.Applicative
+import Control.Monad
+import Data.List
+import Control.Lens
+import Math.LinProg.LPSolve.FFI hiding (solve)
+import qualified Math.LinProg.LPSolve.FFI as F
+import Math.LinProg.Compile
+import Math.LinProg.Types
+import qualified Data.Map as M
+import Prelude hiding (EQ)
+
+solve :: (Eq v, Ord v) => LinProg Double v () -> IO (Either (Maybe ResultCode) [(v, Double)])
+solve (compile -> lp) = do
+    model <- makeLP nconstr nvars
+    case model of
+      Nothing -> return (Left Nothing)
+      Just m' -> with m' $ \m -> do
+
+        -- Eqs
+        forM_ (zip [1..] $ lp ^. equals) $ \(i, eq) ->
+          forM_ (M.keys varLUT) $ \v -> do
+            let w = getVar v $ fst eq
+                c = negate $ snd eq
+            when (w /= 0) $ do
+              setMat m i (varLUT M.! v) w
+              setConstrType m i EQ
+              setRHS m i c
+              return ()
+
+        -- Leqs
+        forM_ (zip [1+nequals..] $ lp ^. leqs) $ \(i, eq) ->
+          forM_ (M.keys varLUT) $ \v -> do
+            let w = getVar v $ fst eq
+                c = negate $ snd eq
+            when (w /= 0) $ do
+              setMat m i (varLUT M.! v) w
+              setConstrType m i LE
+              setRHS m i c
+              return ()
+
+        -- Objective
+        forM_ (M.keys varLUT) $ \v -> do
+          let w = getVar v $ lp ^. objective
+          when (w /= 0) $ void $ setMat m 0 (varLUT M.! v) w
+
+        res <- F.solve m
+        case res of
+          Optimal -> do
+            sol <- snd <$> getSol nvars m
+            return $ Right (zip (M.keys varLUT) sol)
+          _ -> return $ Left (Just res)
+  where
+    nconstr = length allConstr
+    nvars = M.size varLUT
+    nequals = length (lp ^. equals)
+
+    allConstr = (lp ^. equals) ++ (lp ^. leqs)
+    varLUT = M.fromList $ zip (sort $ nub $ concatMap (vars . fst) allConstr) [1..]
+
+    with m f = do
+      r <- f m
+      freeLP m
+      return r
diff --git a/Math/LinProg/LPSolve/FFI.hs b/Math/LinProg/LPSolve/FFI.hs
new file mode 100644
index 0000000..eb818df
--- /dev/null
+++ b/Math/LinProg/LPSolve/FFI.hs
@@ -0,0 +1,89 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+module Math.LinProg.LPSolve.FFI where
+
+import Foreign
+import Foreign.C
+import Foreign.Marshal.Array
+import Control.Applicative ((<$>))
+import qualified Data.Map as M
+
+type LPRec = Ptr ()
+
+data ResultCode =
+  NoMemory
+  |Optimal
+  |SubOptimal
+  |Infeasible
+  |Unbounded
+  |Degenerate
+  |NumFailure
+  |UserAbort
+  |Timeout
+  |PreSolved
+  |ProcFail
+  |ProcBreak
+  |FeasFound
+  |NoFeasFound
+  deriving (Show, Eq, Ord)
+
+data ConstraintType =
+  Fr
+  |LE
+  |GE
+  |EQ
+  deriving (Show, Eq, Ord, Enum)
+
+foreign import ccall "make_lp" c_make_lp :: CInt -> CInt -> IO LPRec
+foreign import ccall "free_lp" c_free_lp :: Ptr LPRec -> IO ()
+foreign import ccall "set_mat" c_set_mat :: LPRec -> CInt -> CInt -> CDouble -> IO CChar
+foreign import ccall "set_rh" c_set_rh :: LPRec -> CInt -> CDouble -> IO CChar
+foreign import ccall "solve" c_solve :: LPRec -> IO CInt
+foreign import ccall "get_variables" c_get_variables :: LPRec -> Ptr CDouble -> IO CChar
+foreign import ccall "set_constr_type" c_set_constr_type :: LPRec -> CInt -> CInt -> IO CChar
+
+setConstrType :: LPRec -> Int -> ConstraintType -> IO Word8
+setConstrType lp i t = fromIntegral <$> c_set_constr_type lp (fromIntegral i) (fromIntegral $ fromEnum t)
+
+makeLP :: Int -> Int -> IO (Maybe LPRec)
+makeLP n m = do
+  m <- c_make_lp (fromIntegral n) (fromIntegral m)
+  return $ if m == nullPtr then
+    Nothing
+  else
+    Just m
+
+freeLP :: LPRec -> IO ()
+freeLP m = with m $ \m' -> c_free_lp m'
+
+setMat :: LPRec -> Int -> Int -> Double -> IO Word8
+setMat a b c d = fromIntegral <$> c_set_mat a (fromIntegral b) (fromIntegral c) (realToFrac d)
+
+setRHS :: LPRec -> Int -> Double -> IO Word8
+setRHS a b c = fromIntegral <$> c_set_rh a (fromIntegral b) (realToFrac c)
+
+solve :: LPRec -> IO ResultCode
+solve lp = (lut M.!) . fromIntegral <$> c_solve lp
+  where
+    lut = M.fromList [
+        (-2, NoMemory)
+        ,(0, Optimal)
+        ,(1, SubOptimal)
+        ,(2, Infeasible)
+        ,(3, Unbounded)
+        ,(4, Degenerate)
+        ,(5, NumFailure)
+        ,(6, UserAbort)
+        ,(7, Timeout)
+        ,(9, PreSolved)
+        ,(10, ProcFail)
+        ,(11, ProcBreak)
+        ,(12, FeasFound)
+        ,(13, NoFeasFound)
+      ]
+
+getSol :: Int -> LPRec -> IO (Word8, [Double])
+getSol n lp = allocaArray (1+n) $ \arr -> do
+    res <- c_get_variables lp arr
+    vals <- peekArray (1+n) arr
+    let vs = map realToFrac vals
+    return (fromIntegral res, vs)
diff --git a/Math/LinProg/Types.hs b/Math/LinProg/Types.hs
index d96fcb5..8b95ef0 100644
--- a/Math/LinProg/Types.hs
+++ b/Math/LinProg/Types.hs
@@ -34,6 +34,10 @@ type LinExpr t v = Fix (LinExpr' t v)
 
 var = Fix . Var
 
+instance Fractional t => Fractional (LinExpr t v) where
+  a / b = Fix (Mul a (1/b))
+  fromRational a = Fix (Lit (fromRational a))
+
 instance Num t => Num (LinExpr t v) where
   a * b = Fix (Mul a b)
   a + b = Fix (Add a b)
@@ -104,3 +108,7 @@ geq b a = liftF (LeqConstraint a b ())
 a =: b = eq a b
 a <: b = leq a b
 a >: b = geq a b
+
+infix 4 =:
+infix 4 <:
+infix 4 >: