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 --- LICENSE | 30 +++++++++++++++ Math/LinProg/Compile.hs | 40 ++++++++++++++++++++ Math/LinProg/LPSolve.hs | 67 ++++++++++++++++++++++++++++++++++ Math/LinProg/LPSolve/FFI.hs | 89 +++++++++++++++++++++++++++++++++++++++++++++ Math/LinProg/Types.hs | 8 ++++ default.nix | 7 +++- lp_solve.nix | 24 ++++++++++++ 7 files changed, 264 insertions(+), 1 deletion(-) create mode 100644 LICENSE create mode 100644 Math/LinProg/Compile.hs create mode 100644 Math/LinProg/LPSolve.hs create mode 100644 Math/LinProg/LPSolve/FFI.hs create mode 100644 lp_solve.nix diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d0e4fe3 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2014, Justin Bedo + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Justin Bedo nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 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 >: diff --git a/default.nix b/default.nix index 6cdf2b1..83c16f8 100644 --- a/default.nix +++ b/default.nix @@ -1,10 +1,14 @@ -{ haskellPackages ? (import {}).haskellPackages }: +{ stdenv ? (import {}).stdenv +, fetchurl ? (import {}).fetchurl +, haskellPackages ? (import {}).haskellPackages }: let inherit (haskellPackages) cabal recursionSchemes lens free; + lpSolve = import ./lp_solve.nix {inherit stdenv fetchurl;}; + in cabal.mkDerivation (self: { pname = "LinProg"; version = "0.0.0.1"; @@ -13,5 +17,6 @@ in cabal.mkDerivation (self: { recursionSchemes lens free + lpSolve ]; }) diff --git a/lp_solve.nix b/lp_solve.nix new file mode 100644 index 0000000..aad43fd --- /dev/null +++ b/lp_solve.nix @@ -0,0 +1,24 @@ +{ stdenv, fetchurl }: + +stdenv.mkDerivation rec { + name = "lp_solve-${version}"; + version = "5.5.2.0"; + + src = fetchurl { + url = "mirror://sourceforge/lpsolve/lpsolve/${version}/lp_solve_${version}_source.tar.gz"; + sha256 = "176c7f023mb6b8bfmv4rfqnrlw88lsg422ca74zjh19i2h5s69sq"; + }; + + buildPhase = '' + cd lpsolve55 + so=1 sh ccc + ''; + + installPhase = '' + mkdir -p $out/include + cp ../*.h $out/include + + mkdir -p $out/lib + cp bin/ux64/* $out/lib + ''; +} -- cgit v1.2.3