From fd68822fda2472f676d835d80d6b17668e52ee45 Mon Sep 17 00:00:00 2001 From: Justin Bedo Date: Tue, 7 Oct 2014 10:55:00 +1100 Subject: Added LPSolve as a solver --- Math/LinProg/Compile.hs | 40 ++++++++++++++++++++ Math/LinProg/LPSolve.hs | 67 ++++++++++++++++++++++++++++++++++ Math/LinProg/LPSolve/FFI.hs | 89 +++++++++++++++++++++++++++++++++++++++++++++ Math/LinProg/Types.hs | 8 ++++ 4 files changed, 204 insertions(+) create mode 100644 Math/LinProg/Compile.hs create mode 100644 Math/LinProg/LPSolve.hs create mode 100644 Math/LinProg/LPSolve/FFI.hs (limited to 'Math/LinProg') 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 >: -- cgit v1.2.3