Added LPSolve as a solver

7 files changed, 264 insertions, 1 deletions

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 <nixpkgs> {}).haskellPackages }:
+{ stdenv ? (import <nixpkgs> {}).stdenv
+, fetchurl ? (import <nixpkgs> {}).fetchurl
+, haskellPackages ? (import <nixpkgs> {}).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
+  '';
+}