diff --git a/cplex-hs.cabal b/cplex-hs.cabal --- a/cplex-hs.cabal +++ b/cplex-hs.cabal @@ -1,5 +1,5 @@ name: cplex-hs -version: 0.3.0.0 +version: 0.4.0.4 synopsis: high-level CPLEX interface -- description: License: BSD3 @@ -21,12 +21,12 @@ exposed-modules: CPLEX.Core, CPLEX.Param CPLEX.Bindings - LSolver.Bindings - LSolver.Dummy - LSolver.Backend.Cplex - LSolver.Problems.MinCostMulticom + Data.LP + Data.LP.Backend.Cplex -- other-modules: Build-depends: base < 5.0 + , unordered-containers == 0.2.7.1 + , hashable , containers < 6.0 , mtl < 2.3 , primitive < 0.7.0.0 diff --git a/src/CPLEX/Bindings.hs b/src/CPLEX/Bindings.hs --- a/src/CPLEX/Bindings.hs +++ b/src/CPLEX/Bindings.hs @@ -29,6 +29,7 @@ , c_CPXgetnumrows , c_CPXgetobjval , c_CPXgetx + , c_CPXgetbase , c_CPXgetstat , c_CPXgetslack , c_CPXcopylp @@ -52,6 +53,8 @@ , c_CPXsetincumbentcallbackfunc , c_CPXsetcutcallbackfunc , c_CPXsetlazyconstraintcallbackfunc + -- More MIP + , c_CPXgetmiprelgap -- MIP cuts , c_CPXaddmipstarts , c_CPXcutcallbackadd @@ -59,6 +62,7 @@ , c_CPXaddlazyconstraints , c_CPXgetcallbacknodex , c_CPXgetcallbacknodelp + , c_CPXgetcallbackinfo ) where import Foreign.C (CChar (..), CDouble (..), CInt (..)) @@ -225,6 +229,11 @@ foreign import ccall safe "cplex.h CPXgetcallbacknodelp" c_CPXgetcallbacknodelp :: Ptr CpxEnv' -> Ptr () -> CInt -> Ptr (Ptr CpxLp') -> IO CInt +--int CPXXgetcallbackinfo( CPXCENVptr env, void * cbdata, int wherefrom, int +--whichinfo, void * result_p ) +foreign import ccall safe "cplex.h CPXgetcallbackinfo" c_CPXgetcallbackinfo :: + Ptr CpxEnv' -> Ptr () -> CInt -> CInt -> Ptr () -> IO CInt + --int CPXgetcallbacknodex(CPXCENVptr env, void * cbdata, int wherefrom, double * x, int begin, int end) foreign import ccall safe "cplex.h CPXgetcallbacknodex" c_CPXgetcallbacknodex :: Ptr CpxEnv' -> Ptr () -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CInt @@ -300,4 +309,8 @@ c_createCutCallbackPtr :: CCutCallback -> IO (FunPtr (CCutCallback)) --int CPXsetincumbentcallbackfunc(CPXENVptr env, int(*)(CALLBACK_INCUMBENT_ARGS) incumbentcallback, void * cbhandle) +foreign import ccall safe "cplex.h CPXgetmiprelgap" c_CPXgetmiprelgap :: + Ptr CpxEnv' -> Ptr CpxLp' -> Ptr CDouble -> IO CInt +foreign import ccall safe "cplex.h CPXgetbase" c_CPXgetbase :: + Ptr CpxEnv' -> Ptr CpxLp' -> Ptr CInt -> Ptr CInt -> IO CInt diff --git a/src/CPLEX/Core.hs b/src/CPLEX/Core.hs --- a/src/CPLEX/Core.hs +++ b/src/CPLEX/Core.hs @@ -1,5 +1,6 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE ScopedTypeVariables #-} module CPLEX.Core ( CpxEnv(..) , CpxLp @@ -29,8 +30,6 @@ , dualopt , siftopt , hybnetopt - , getSolution - , getMIPSolution , writeprob -- * change things , changeCoefList @@ -45,6 +44,7 @@ , getNumRows , getErrorString , getStatString + , getBaseVars -- MIP , setIncumbentCallback , setCutCallback @@ -56,6 +56,10 @@ , getCallbackNodeX , addCutFromCallback , addSingleMIPStart + , getSolution + , getMIPSolution + , getMipRelGap + , getMipBestInteger -- * convenience wrappers , withEnv , withLp @@ -79,6 +83,7 @@ import CPLEX.Bindings import CPLEX.Param import Data.Char(ord) +import Unsafe.Coerce newtype CpxEnv = CpxEnv (Ptr CpxEnv') newtype CpxLp = CpxLp (Ptr CpxLp') @@ -333,11 +338,40 @@ k -> fmap Left (getErrorString env (CpxRet k)) k -> fmap Left (getErrorString env (CpxRet k)) + +getMipRelGap :: CpxEnv -> Ptr () -> CInt -> IO Double +getMipRelGap (CpxEnv env') cbdata wherefrom = do + gap_p :: Ptr CDouble <- malloc + status <- c_CPXgetcallbackinfo env' cbdata wherefrom 125 (unsafeCoerce gap_p) + objVal :: CDouble <- peek gap_p + free gap_p + return $ if status == 0 then realToFrac objVal else fromIntegral status + +getMipBestInteger :: CpxEnv -> Ptr () -> CInt -> IO Double +getMipBestInteger (CpxEnv env') cbdata wherefrom = do + gap_p :: Ptr CDouble <- malloc + status <- c_CPXgetcallbackinfo env' cbdata wherefrom 101 (unsafeCoerce gap_p) + objVal :: CDouble <- peek gap_p + free gap_p + return $ if status == 0 then realToFrac objVal else fromIntegral status + writeprob :: CpxEnv -> CpxLp -> String -> IO (Maybe String) writeprob env@(CpxEnv env') lp@(CpxLp lp') filename = do fn <- newCAString filename status <- c_CPXwriteprob env' lp' fn nullPtr getErrorStatus env status + +getBaseVars :: CpxEnv -> CpxLp -> IO (Maybe (Vector Int)) +getBaseVars env@(CpxEnv env') lp@(CpxLp lp') = do + numcols <- getNumCols env lp + x <- VSM.new numcols + VSM.unsafeWith x $ \x' -> do + status <- c_CPXgetbase env' lp' x' nullPtr + case status of + 0 -> do vec <- VS.freeze x + let vec' = VS.map fromIntegral vec + return $ Just $ vec' + _ -> return Nothing toCpxError :: CpxEnv -> CInt -> IO (Maybe String) toCpxError env 0 = return Nothing diff --git a/src/CPLEX/Param.hs b/src/CPLEX/Param.hs --- a/src/CPLEX/Param.hs +++ b/src/CPLEX/Param.hs @@ -85,6 +85,11 @@ | CPX_PARAM_SOLUTIONTARGET | CPX_PARAM_CLONELOG | CPX_PARAM_MIPCBREDLP + | CPX_PARAM_MIPDISPLAY + | CPX_PARAM_EPGAP + | CPX_PARAM_EPAGAP + | CPX_PARAM_INTSOLLIM + | CPX_PARAM_NODELIM deriving Show paramToInt :: Num a => CPX_PARAM -> a @@ -155,6 +160,11 @@ paramToInt CPX_PARAM_APIENCODING = 1130 paramToInt CPX_PARAM_SOLUTIONTARGET = 1131 paramToInt CPX_PARAM_CLONELOG = 1132 +paramToInt CPX_PARAM_EPAGAP = 2008 +paramToInt CPX_PARAM_EPGAP = 2009 +paramToInt CPX_PARAM_MIPDISPLAY = 2012 +paramToInt CPX_PARAM_INTSOLLIM = 2015 +paramToInt CPX_PARAM_NODELIM = 2017 paramToInt CPX_PARAM_MIPCBREDLP = 2055 data CPX_PROB_TYPE = CPX_PROB_LP diff --git a/src/Data/LP.hs b/src/Data/LP.hs new file mode 100644 --- /dev/null +++ b/src/Data/LP.hs @@ -0,0 +1,104 @@ +{-# LANGUAGE FlexibleInstances #-} + +module Data.LP(Variable(..) + ,Bound(..) + ,Constraints(..) + ,Optimization(..) + ,(<+>) + ,Bounds(..) + ,Type(..) + ,MixedIntegerProblem(..) + ,LinearProblem(..) + ,MIPSolution(..) + ,LPSolution(..) + ,simplifyConstraints + ,removeEmptyConstraints + ) where + +import Data.List (intercalate) +import qualified Data.Vector as V +import qualified Data.HashMap.Strict as M +import Data.Hashable +import Data.Monoid +import qualified Data.HashSet as S + +type Map k v = M.HashMap k v + +data Variable a = Double :# a + +data Bound x = x :< Double + | x :> Double + | x := Double + deriving Show + +newtype Constraints a = Constraints [ Bound [Variable a] ] + +simplifyVars :: (Eq a, Hashable a) => [Variable a] -> [Variable a] +simplifyVars vars = map (\(v,c) -> c :# v) $ M.toList $ + foldr (\(c :# v) m -> M.insertWith (+) v c m) M.empty vars +simplifyBounds (xs :< b) = (simplifyVars xs) :< b +simplifyBounds (xs := b) = (simplifyVars xs) := b +simplifyBounds (xs :> b) = (simplifyVars xs) :> b + +simplifyConstraints :: (Eq a, Hashable a) => Constraints a -> Constraints a +simplifyConstraints (Constraints cs) = Constraints $ map simplifyBounds cs + +removeEmptyConstraints :: (Eq a, Hashable a) => Constraints a -> Constraints a +removeEmptyConstraints (Constraints cs) = Constraints $ filter isNonEmpty cs + where + isNonEmpty ([] :< b) = False + isNonEmpty ([] := b) = False + isNonEmpty ([] :> b) = False + isNonEmpty _ = True + +instance Monoid a => Monoid (Constraints a) where + (Constraints xs) `mappend` (Constraints ys) = Constraints $ xs <> ys + mempty = Constraints [] + +Constraints v1 <+> Constraints v2 = Constraints $ v1 ++ v2 + +data Optimization a = Maximize [Variable a] + | Minimize [Variable a] + +data Type = TContinuous | TInteger | TBinary + +instance (Show a) => Show (Variable a) where + show (d :# v) + | d == (-1) = "-" ++ (show v) + | d == 1 = (show v) + | otherwise = (show d) ++ "x" ++ (show v) + +instance Show a => Show (Optimization a) where + show (Minimize xs) = "Minimize\n\t" ++ (intercalate "+" $ map show xs) + show (Maximize xs) = "Maximize\n\t" ++ (intercalate "+" $ map show xs) + +showVars xs = intercalate " + " $ map show $ zipWith (:#) xs [0..] + +instance (Show a) => Show (Constraints a) where + show (Constraints bounds) = "\nSubject to\n" ++ (unlines $ map (\a -> "\t" ++ a) $ + map getVarSigns bounds) + +printVars xs = intercalate " + " $ map show xs +getVarSigns (x :< v) = (printVars x) ++ " <= " ++ (show v) +getVarSigns (x :> v) = (printVars x) ++ " >= " ++ (show v) +getVarSigns (x := v) = (printVars x) ++ " == " ++ (show v) + +instance Show Type where + show TContinuous = "Continous" + show TInteger = "Integer" + show TBinary = "Binary" + + +type Bounds = [Bound Int] + +data LinearProblem a = LP (Optimization a) (Constraints a) [(a, Maybe Double, Maybe Double)] + deriving Show + +data MixedIntegerProblem a = MILP (Optimization a) (Constraints a) [(a, Maybe Double, Maybe Double)] + [(a,Type)] + deriving Show + +data MIPSolution a = MIPSolution { mipOptimalSol :: Bool, mipObjVal :: Double, mipVars :: Map a Double } deriving (Show) + +data LPSolution a = LPSolution { lpOptimalSol :: Bool, lpObjVal :: Double, lpVars :: Map a Double, lpDualVars :: V.Vector Double, lpBasisVars :: Maybe (S.HashSet a)} deriving (Show) + diff --git a/src/Data/LP/Backend/Cplex.hs b/src/Data/LP/Backend/Cplex.hs new file mode 100644 --- /dev/null +++ b/src/Data/LP/Backend/Cplex.hs @@ -0,0 +1,340 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE RecordWildCards #-} + +module Data.LP.Backend.Cplex(solLP, standardBounds, defaultCallBacks, getCallBackLp, getIncCallBackXs, getCallBackXs + ,addCallBackCut + ,getCallBackGap + ,getCallBackBestObjI + ,solMIP + ,solMIP' + ,UserCutCallBack, CutCallBackM, UserIncumbentCallBack, IncumbentCallBackM, CallBacks(..)) where + +import qualified Data.Vector as V +import CPLEX.Bindings +import CPLEX.Param +import CPLEX.Core hiding (Bound) +--import Foreign.C (CInt) +import Data.LP +import qualified Data.Vector.Storable as VS +import Foreign.Ptr +import Foreign.ForeignPtr(newForeignPtr_) +import Foreign.Storable +import Foreign.C +import Control.Monad +import Control.Monad.Reader +import qualified Data.HashMap.Strict as M +import Data.Hashable +import Data.List (sortBy) +import Data.Ord (comparing) +import qualified Data.HashSet as S + + +type Map k v = M.HashMap k v + +data CallBacks a = ActiveCallBacks {cutcb :: Maybe (UserCutCallBack a), inccb :: Maybe (UserIncumbentCallBack), + lazycb :: Maybe (UserCutCallBack a) } +defaultCallBacks :: CallBacks a +defaultCallBacks = ActiveCallBacks {cutcb = Nothing, inccb = Nothing, lazycb = Nothing} + +type ParamValues = [(CPX_PARAM, Int)] + +type VarDic a = Map a Int +type RevDic a = Map Int a +data CutCallBackArgs a = CutCallBackArgs {env :: CpxEnv, cbdata :: Ptr (), wherefrom :: CInt, cbhandle :: Ptr (), userdata :: Ptr Int, vardic :: VarDic a, revdic :: RevDic a} +type CutCallBackM b a = (ReaderT (CutCallBackArgs b) IO a) +type UserCutCallBack a = CutCallBackM a Int + +data IncumbentCallBackArgs = IncumbentCallBackArgs {envi :: CpxEnv, cbdatai :: Ptr (), wherefromi :: CInt, cbhandlei :: Ptr (), + objVal :: CDouble, xs :: Ptr CDouble, isfeas :: Ptr Int , useraction :: Ptr Int} +type IncumbentCallBackM a = (ReaderT IncumbentCallBackArgs IO a) +type UserIncumbentCallBack = Double -> VS.Vector Double -> IncumbentCallBackM Bool + + + +incumbentcallback :: UserIncumbentCallBack -> CIncumbentCallback +incumbentcallback usercb env' cbdata wherefrom cbhandle objVal xs isfeas useraction = do + let env = CpxEnv env' + let oval = realToFrac objVal + foreignPtr <- newForeignPtr_ xs + lp <- getCallbackLP env cbdata (fromIntegral wherefrom) + xs' <- case lp of + Right lp' -> do + colCount <- getNumCols env lp' + return $ VS.map realToFrac $ VS.unsafeFromForeignPtr0 foreignPtr colCount + Left _ -> return VS.empty + isFeas <- runReaderT (usercb oval xs') $ IncumbentCallBackArgs env cbdata wherefrom cbhandle objVal xs isfeas useraction + poke isfeas (if isFeas then 1 else 0) + return 0 + + +cutcallback :: (Eq a, Hashable a) => VarDic a -> RevDic a -> UserCutCallBack a -> CCutCallback +cutcallback vardic revdic usercb env' cbdata wherefrom cbhandle ptrUser = do + let env = CpxEnv env' + runReaderT usercb $ CutCallBackArgs env cbdata wherefrom cbhandle ptrUser vardic revdic + +lazycallback :: (Eq a, Hashable a) => VarDic a -> RevDic a -> UserCutCallBack a -> CCutCallback +lazycallback vardic revdic usercb env' cbdata wherefrom cbhandle ptrUser = do + let env = CpxEnv env' + runReaderT usercb $ CutCallBackArgs env cbdata wherefrom cbhandle ptrUser vardic revdic + + +getCallBackLp :: (Eq a, Hashable a) => CutCallBackM a CpxLp +getCallBackLp = do + CutCallBackArgs{..} <- ask + res <- liftIO $ getCallbackLP env cbdata (fromIntegral wherefrom) + case res of + Right lp -> return lp + _ -> error "I cant get the LP for some reason" + +getIncCallBackXs :: IncumbentCallBackM (V.Vector Double) +getIncCallBackXs = do + IncumbentCallBackArgs{..} <- ask + + lp <- liftIO $ getCallbackLP envi cbdatai (fromIntegral wherefromi) + xs' <- case lp of + Right lp' -> do + colCount <- liftIO $ getNumCols envi lp' + (stat, xsVS) <- liftIO $ getCallbackNodeX envi cbdatai (fromIntegral wherefromi) 0 (colCount-1) + return $ VS.convert xsVS + Left msg -> return $ V.empty + return xs' + +getCallBackXs :: (Eq a, Hashable a) => CutCallBackM a (Map a Double) +getCallBackXs = do + CutCallBackArgs{..} <- ask + + lp <- liftIO $ getCallbackLP env cbdata (fromIntegral wherefrom) + xs' <- case lp of + Right lp' -> do + colCount <- liftIO $ getNumCols env lp' + (stat, xsVS) <- liftIO $ getCallbackNodeX env cbdata (fromIntegral wherefrom) 0 (colCount-1) + return $ VS.convert xsVS + Left msg -> return $ V.empty + let vars = V.toList xs' + let m = M.fromList $ zip (map (revdic M.!) [0..length vars - 1]) vars + return m + +getCallBackGap :: (Eq a, Hashable a) => CutCallBackM a Double +getCallBackGap = do + CutCallBackArgs{..} <- ask + gap <- liftIO $ getMipRelGap env cbdata wherefrom + return gap + +getCallBackBestObjI :: (Eq a, Hashable a) => CutCallBackM a Double +getCallBackBestObjI = do + CutCallBackArgs{..} <- ask + gap <- liftIO $ getMipBestInteger env cbdata wherefrom + return gap + +addCallBackCut :: (Eq a, Hashable a) => Bound [Variable a] -> CutCallBackM a (Maybe String) +addCallBackCut st_ = do + CutCallBackArgs{..} <- ask + let st = tokenizeVars st_ vardic + let (cnstrs,rhs) = toForm st + liftIO $ addCutFromCallback env cbdata wherefrom (fromIntegral $ length cnstrs) rhs cnstrs CPX_USECUT_FORCE + where + toForm (vars :< b) = (map (\(v :# i) -> (Col i, v)) vars, L b) + toForm (vars :> b) = (map (\(v :# i) -> (Col i, v)) vars, G b) + toForm (vars := b) = (map (\(v :# i) -> (Col i, v)) vars, E b) + +standardBounds :: (Int, Int) -> [(Int, Maybe Double, Maybe Double)] +standardBounds (i,j) = [(i',Just 0, Nothing) | i' <- [i..j] ] + +toBounds :: [(Int, Maybe Double, Maybe Double)] -> (Int, Int) -> V.Vector (Maybe Double, Maybe Double) +toBounds bounds vr@(a,b) = def V.// (map (\(q,w,e) -> (q,(w,e))) bounds) + where + def = V.fromList [k | _ <- [a..b], let k = (Just 0, Nothing)] + +toConstraints :: Constraints Int -> (Int, Int) -> ([(Row, Col, Double)], V.Vector Sense) +toConstraints constraints varRange = let (st, rhs) = toStandard constraints 0 [] [] varRange + in (st, V.fromList rhs) + + +toStandard :: Constraints Int -> Int -> [(Row, Col, Double)] -> [Sense] -> (Int, Int) + -> ([(Row, Col, Double)], [Sense]) +toStandard (Constraints []) _ accSt accRhs _ = (reverse $ accSt, reverse $ accRhs) +toStandard (Constraints (b:bs)) rowI accSt accRhs varRange = case b of + vars :< boundVal -> addRow vars L boundVal + vars := boundVal -> addRow vars E boundVal + vars :> boundVal -> addRow vars G boundVal + where addRow vars s boundVal = toStandard (Constraints bs) (rowI+1) (generateRow vars ++ accSt) + ((s boundVal) : accRhs) varRange + generateRow [] = [] + generateRow ((v :# i):vs) = (Row rowI, Col i, v):generateRow vs + + +toObj :: Optimization Int -> (ObjSense, V.Vector Double) +toObj (Maximize dd) = (CPX_MAX, varsToVector dd) +toObj (Minimize dd) = (CPX_MIN, varsToVector dd) + +-- This assumes that all elements are in !, zero pad boys +varsToVector :: [Variable Int] -> V.Vector Double +varsToVector vs = V.fromList $ map snd $ sortBy (comparing fst) $ map (\(c :# i) -> (i,c)) vs + + +solLP :: (Eq a, Hashable a) => LinearProblem a -> ParamValues -> IO (LPSolution a) +solLP (LP objective_ constraints_ bounds_) params = withEnv $ \env -> do + --setIntParam env CPX_PARAM_SCRIND cpx_ON + --setIntParam env CPX_PARAM_DATACHECK cpx_ON + mapM_ (\(p,v) -> setIntParam env p (fromIntegral v)) params + withLp env "testprob" $ \lp -> do + let + dic = generateVarDic constraints_ objective_ bounds_ + revDic = M.fromList $ map (\(a,b) -> (b,a)) $ M.toList dic + objective = tokenizeObj objective_ dic + constraints = tokenizeConstraints constraints_ dic + bounds = tokenizeBounds bounds_ dic + varRange = (0,M.size dic - 1) + + (varA, varB) = varRange + varCount = varB - varA + 1 + (objsen, obj) = toObj objective + (cnstrs,rhs) = toConstraints constraints varRange + xbnds = toBounds bounds varRange + + statusLp <- copyLp env lp objsen obj rhs cnstrs xbnds + + case statusLp of + Nothing -> return () + Just msg -> error $ "CPXcopylp error: " ++ msg + + statusOpt <- qpopt env lp + case statusOpt of + Nothing -> return () + Just msg -> error $ "CPXqpopt error: " ++ msg + + statusSol <- getSolution env lp + case statusSol of + Left msg -> error $ "CPXsolution error: " ++ msg + Right sol -> do + let vars = V.toList $ VS.convert $ solX sol + let m = M.fromList $ zip (map (revDic M.!) [0..length vars - 1]) vars + basism <- getBaseVars env lp + let basis' = basism >>= \basis -> Just $ S.fromList $ map (\(i,c) -> revDic M.! i) $ filter(\(i,c) -> c == 1) $ zip [0..] $ V.toList $ VS.convert basis + return $ LPSolution (solStat sol == CPX_STAT_OPTIMAL) (solObj sol) ( m ) (VS.convert $ solPi sol) basis' + +solMIP :: (Eq a, Hashable a) => MixedIntegerProblem a -> ParamValues -> CallBacks a -> IO (MIPSolution a) +solMIP = solMIP' M.empty +solMIP' :: (Eq a, Hashable a) => Map a Double -> MixedIntegerProblem a -> ParamValues -> CallBacks a -> IO (MIPSolution a) +solMIP' warmStart (MILP objective_ constraints_ bounds_ types_ ) params (ActiveCallBacks {..}) = withEnv $ \env -> do +-- setIntParam env CPX_PARAM_SCRIND 1 + -- setIntParam env CPX_PARAM_DATACHECK 1 + mapM_ (\(p,v) -> setIntParam env p (fromIntegral v)) params + withLp env "clu" $ \lp -> do + let + dic = generateVarDic constraints_ objective_ bounds_ + revDic = M.fromList $ map (\(a,b) -> (b,a)) $ M.toList dic + objective = tokenizeObj objective_ dic + constraints = tokenizeConstraints constraints_ dic + bounds = tokenizeBounds bounds_ dic + types = tokenizeTypes types_ dic + varRange = (0,M.size dic - 1) + + (varA, varB) = varRange + varCount = varB - varA + 1 + (objsen, obj) = toObj objective + (cnstrs,rhs) = toConstraints constraints varRange + xbnds = toBounds bounds varRange + + types' = V.fromList (replicate varCount CPX_CONTINUOUS) V.// (map (\(a,t) -> (a, typeToCPX t)) types) + + statusLp <- copyMip env lp objsen obj rhs cnstrs xbnds types' + + case statusLp of + Nothing -> return () + Just msg -> error $ "CPXcopyMIP error: " ++ msg + + case inccb of + Just cb -> do statusIncCB <- setIncumbentCallback env (incumbentcallback cb) + case statusIncCB of + Nothing -> return () + Just msg -> error $ "CPXIncumbentCallBackSet Error: " ++ msg + Nothing -> return () + + case lazycb of + Just cb -> do statusLazyCB <- setLazyConstraintCallback env (lazycallback dic revDic cb) + case statusLazyCB of + Nothing -> return () + Just msg -> error $ "CPXLazyConstraintCallBackSet Error: " ++ msg + Nothing -> return () + + case cutcb of + Just cb -> do statusCutCB <- setCutCallback env (cutcallback dic revDic cb) + case statusCutCB of + Nothing -> return () + Just msg -> error $ "CPXCutCallBackSet Error: " ++ msg + Nothing -> return () + + let warmVars = M.toList warmStart + let warmInd = map (dic M.!) $ map fst warmVars + let warmVals = map snd warmVars + when (M.size warmStart > 0) $ addSingleMIPStart env lp warmInd warmVals CPX_MIPSTART_AUTO >> return () + + statusOpt <- mipopt env lp + case statusOpt of + Nothing -> return () + Just msg -> error $ "CPXmipopt error: " ++ msg + + statusSol <- getMIPSolution env lp + case statusSol of + Left msg -> error $ "CPXsolution error: " ++ msg + Right sol -> do -- I call this imperative do notation + let vars = V.toList $ VS.convert $ solX sol + let m = M.fromList $ zip (map (revDic M.!) [0..length vars - 1]) vars + return $ MIPSolution (solStat sol == CPXMIP_OPTIMAL) (solObj sol) m + +typeToCPX :: Data.LP.Type -> CPLEX.Core.Type +typeToCPX (TInteger) = CPX_INTEGER +typeToCPX (TContinuous) = CPX_CONTINUOUS +typeToCPX (TBinary) = CPX_BINARY + + + + +generateVarDic :: (Eq a, Hashable a) => Constraints a -> Optimization a -> [(a, Maybe Double, Maybe Double)] + -> VarDic a +generateVarDic (Constraints bounds) opt bs = let allvars = (concat $ map getBounds bounds) ++ (getOptim opt) ++ (map fst' bs) + in go allvars 0 M.empty + where + getVar (_ :# v) = v + getBounds (vs :< _ ) = map getVar vs + getBounds (vs :> _ ) = map getVar vs + getBounds (vs := _ ) = map getVar vs + getOptim (Maximize a) = map getVar a + getOptim (Minimize a) = map getVar a + fst' (a,b,c) = a + go :: (Eq a, Hashable a) => [a] -> Int -> Map a Int -> Map a Int + go [] i m = m + go (v:vs) i m + | v `M.member` m = go vs i m + | otherwise = go vs (i+1) $ M.insertWith (\new old -> old) v i m + +-- Change variables in bound to have ID +tokenizeConstraints :: (Eq a, Hashable a) => Constraints a -> Map a Int -> Constraints Int +tokenizeConstraints (Constraints bounds) dic = Constraints $ map b2b bounds + where + v2v (d :# v) = d :# (dic M.! v) + b2b (vs :< d) = map v2v vs :< d + b2b (vs :> d) = map v2v vs :> d + b2b (vs := d) = map v2v vs := d + +tokenizeVars :: (Eq a, Hashable a) => Bound [Variable a] -> Map a Int -> Bound [Variable Int] +tokenizeVars bounds dic = b2b bounds + where + v2v (d :# v) = d :# (dic M.! v) + b2b (vs :< d) = map v2v vs :< d + b2b (vs :> d) = map v2v vs :> d + b2b (vs := d) = map v2v vs := d + +tokenizeObj :: (Eq a, Hashable a) => Optimization a -> Map a Int -> Optimization Int +tokenizeObj obj dic = o2o obj + where + v2v (d :# v) = d :# (dic M.! v) + o2o (Minimize vs) = Minimize $ map v2v vs + o2o (Maximize vs) = Maximize $ map v2v vs + +tokenizeBounds :: (Eq a, Hashable a) => [(a,b,b)] -> Map a Int -> [(Int,b,b)] +tokenizeBounds xs dic = map (\(a,b,c) -> (dic M.! a, b, c) ) xs + +tokenizeTypes :: (Eq a, Hashable a) => [(a,b)] -> Map a Int -> [(Int, b)] +tokenizeTypes xs dic = map (\(a,b) -> (dic M.! a, b) ) xs diff --git a/src/LSolver/Backend/Cplex.hs b/src/LSolver/Backend/Cplex.hs deleted file mode 100644 --- a/src/LSolver/Backend/Cplex.hs +++ /dev/null @@ -1,287 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} -{-# LANGUAGE RecordWildCards #-} - -module LSolver.Backend.Cplex(solLP, standardBounds, solMIP, defaultCallBacks, getCallBackLp, getIncCallBackXs, getCallBackXs, - addCallBackCut, UserCutCallBack, CutCallBackM, UserIncumbentCallBack, IncumbentCallBackM, CallBacks(..)) where - -import Data.Ix as I -import qualified Data.Vector as V -import qualified Data.Sequence as S -import Data.Foldable as F -import CPLEX.Bindings -import CPLEX.Param -import CPLEX.Core hiding (Bound) ---import Foreign.C (CInt) -import LSolver.Bindings -import qualified Data.Vector.Storable as VS -import Foreign.Ptr -import Foreign.ForeignPtr(newForeignPtr_) -import Foreign.Storable -import Foreign.C -import Control.Monad -import Control.Monad.Reader -import qualified Data.Map as M -import Data.List (sortBy) -import Data.Ord (comparing) - -data CallBacks = ActiveCallBacks {cutcb :: Maybe UserCutCallBack, inccb :: Maybe UserIncumbentCallBack, - lazycb :: Maybe UserCutCallBack} -defaultCallBacks = ActiveCallBacks {cutcb = Nothing, inccb = Nothing, lazycb = Nothing} - -type ParamValues = [(CPX_PARAM, Int)] - -data CutCallBackArgs = CutCallBackArgs {env :: CpxEnv, cbdata :: Ptr (), wherefrom :: CInt, cbhandle :: Ptr (), userdata :: Ptr Int} -type CutCallBackM a = (ReaderT CutCallBackArgs IO a) -type UserCutCallBack = CutCallBackM Int - -data IncumbentCallBackArgs = IncumbentCallBackArgs {envi :: CpxEnv, cbdatai :: Ptr (), wherefromi :: CInt, cbhandlei :: Ptr (), - objVal :: CDouble, xs :: Ptr CDouble, isfeas :: Ptr Int , useraction :: Ptr Int} -type IncumbentCallBackM a = (ReaderT IncumbentCallBackArgs IO a) -type UserIncumbentCallBack = Double -> VS.Vector Double -> IncumbentCallBackM Bool - - - -incumbentcallback :: UserIncumbentCallBack -> CIncumbentCallback -incumbentcallback usercb env' cbdata wherefrom cbhandle objVal xs isfeas useraction = do - let env = CpxEnv env' - let oval = realToFrac objVal - foreignPtr <- newForeignPtr_ xs - lp <- getCallbackLP env cbdata (fromIntegral wherefrom) - xs' <- case lp of - Right lp' -> do - colCount <- getNumCols env lp' - return $ VS.map realToFrac $ VS.unsafeFromForeignPtr0 foreignPtr colCount - Left _ -> return VS.empty - isFeas <- runReaderT (usercb oval xs') $ IncumbentCallBackArgs env cbdata wherefrom cbhandle objVal xs isfeas useraction - poke isfeas (if isFeas then 1 else 0) - return 0 - - -cutcallback :: UserCutCallBack -> CCutCallback -cutcallback usercb env' cbdata wherefrom cbhandle ptrUser = do - let env = CpxEnv env' - runReaderT usercb $ CutCallBackArgs env cbdata wherefrom cbhandle ptrUser - -lazycallback :: UserCutCallBack -> CCutCallback -lazycallback usercb env' cbdata wherefrom cbhandle ptrUser = do - let env = CpxEnv env' - runReaderT usercb $ CutCallBackArgs env cbdata wherefrom cbhandle ptrUser - - -getCallBackLp :: CutCallBackM CpxLp -getCallBackLp = do - CutCallBackArgs{..} <- ask - res <- liftIO $ getCallbackLP env cbdata (fromIntegral wherefrom) - case res of - Right lp -> return lp - _ -> error "I cant get the LP for some reason" - -getIncCallBackXs :: IncumbentCallBackM (V.Vector Double) -getIncCallBackXs = do - IncumbentCallBackArgs{..} <- ask - - lp <- liftIO $ getCallbackLP envi cbdatai (fromIntegral wherefromi) - xs' <- case lp of - Right lp' -> do - colCount <- liftIO $ getNumCols envi lp' - (stat, xsVS) <- liftIO $ getCallbackNodeX envi cbdatai (fromIntegral wherefromi) 0 (colCount-1) - return $ VS.convert xsVS - Left msg -> return $ V.empty - return xs' - -getCallBackXs :: CutCallBackM (V.Vector Double) -getCallBackXs = do - CutCallBackArgs{..} <- ask - - lp <- liftIO $ getCallbackLP env cbdata (fromIntegral wherefrom) - xs' <- case lp of - Right lp' -> do - colCount <- liftIO $ getNumCols env lp' - (stat, xsVS) <- liftIO $ getCallbackNodeX env cbdata (fromIntegral wherefrom) 0 (colCount-1) - return $ VS.convert xsVS - Left msg -> return $ V.empty - return xs' - -addCallBackCut :: Bound [Variable Int] -> CutCallBackM (Maybe String) -addCallBackCut st = do - CutCallBackArgs{..} <- ask - let (cnstrs,rhs) = toForm st - liftIO $ addCutFromCallback env cbdata wherefrom (fromIntegral $ length cnstrs) rhs cnstrs CPX_USECUT_FORCE - where - toForm (vars :< b) = (map (\(v :# i) -> (Col i, v)) vars, L b) - toForm (vars :> b) = (map (\(v :# i) -> (Col i, v)) vars, G b) - toForm (vars := b) = (map (\(v :# i) -> (Col i, v)) vars, E b) - -standardBounds :: (Enum t, Num a) => (t, t) -> [(t, Maybe a, Maybe a1)] -standardBounds (i,j) = map (\i' -> (i', Just 0, Nothing)) [i..j] - -toBounds :: (Num a, Ix t) => [(t, Maybe a, Maybe a1)] -> (t, t) -> [(Maybe a, Maybe a1)] -toBounds bounds varRange = F.toList $ aux bounds def - where - def = S.fromList [k | _ <- I.range varRange, let k = (Just 0, Nothing)] - aux [] s = s - aux ((b,lb,ub):bs) s = aux bs (S.update (I.index varRange b) (lb,ub) s) - -toConstraints :: Ix a => Constraints a -> (a, a) -> ([(Row, Col, Double)], V.Vector Sense) -toConstraints constraints varRange = let (st, rhs) = toStandard constraints 0 [] [] varRange - in (st, V.fromList rhs) - - -toStandard :: Ix a => Constraints a -> Int -> [(Row, Col, Double)] -> [Sense] -> (a, a) - -> ([(Row, Col, Double)], [Sense]) -toStandard (Constraints []) _ accSt accRhs _ = (reverse $ accSt, reverse $ accRhs) -toStandard (Constraints (b:bs)) rowI accSt accRhs varRange = case b of - vars :< boundVal -> addRow vars L boundVal - vars := boundVal -> addRow vars E boundVal - vars :> boundVal -> addRow vars G boundVal - where addRow vars s boundVal = toStandard (Constraints bs) (rowI+1) (generateRow vars ++ accSt) - ((s boundVal) : accRhs) varRange - generateRow [] = [] - generateRow ((v :# i):vs) = (Row rowI, Col $ I.index varRange i, v):generateRow vs - - -toObj :: Optimization Int -> (ObjSense, V.Vector Double) -toObj (Maximize dd) = (CPX_MAX, varsToVector dd) -toObj (Minimize dd) = (CPX_MIN, varsToVector dd) - --- This assumes that all elements are in !, zero pad boys -varsToVector :: [Variable Int] -> V.Vector Double -varsToVector vs = V.fromList $ map snd $ sortBy (comparing fst) $ map (\(c :# i) -> (i,c)) vs - - -solLP :: (Eq a, Ord a) => LinearProblem a -> ParamValues -> IO (LPSolution a) -solLP (LP objective_ constraints_ bounds_) params = withEnv $ \env -> do - --setIntParam env CPX_PARAM_SCRIND cpx_ON - --setIntParam env CPX_PARAM_DATACHECK cpx_ON - mapM_ (\(p,v) -> setIntParam env p (fromIntegral v)) params - withLp env "testprob" $ \lp -> do - let - dic = generateVarDic constraints_ - revDic = M.fromList $ map (\(a,b) -> (b,a)) $ M.toList dic - objective = tokenizeObj objective_ dic - constraints = tokenizeConstraints constraints_ dic - bounds = tokenizeBounds bounds_ dic - varRange = (0,M.size dic) - (objsen, obj) = toObj objective - (cnstrs,rhs) = toConstraints constraints varRange - xbnds = toBounds bounds varRange - - statusLp <- copyLp env lp objsen obj rhs cnstrs (V.fromList xbnds) - - case statusLp of - Nothing -> return () - Just msg -> error $ "CPXcopylp error: " ++ msg - - statusOpt <- qpopt env lp - case statusOpt of - Nothing -> return () - Just msg -> error $ "CPXqpopt error: " ++ msg - - statusSol <- getSolution env lp - case statusSol of - Left msg -> error $ "CPXsolution error: " ++ msg - Right sol -> do - let vars = V.toList $ VS.convert $ solX sol - let m = M.fromList $ zip (map (revDic M.!) [0..length vars - 1]) vars - return $ LPSolution (solStat sol == CPX_STAT_OPTIMAL) (solObj sol) ( m ) (VS.convert $ solPi sol) - -solMIP :: (Ord a, Eq a) => MixedIntegerProblem a -> ParamValues -> CallBacks -> IO (MIPSolution a) -solMIP (MILP objective_ constraints_ bounds_ types_ ) params (ActiveCallBacks {..}) = withEnv $ \env -> do --- setIntParam env CPX_PARAM_SCRIND 1 - -- setIntParam env CPX_PARAM_DATACHECK 1 - mapM_ (\(p,v) -> setIntParam env p (fromIntegral v)) params - withLp env "clu" $ \lp -> do - let - dic = generateVarDic constraints_ - revDic = M.fromList $ map (\(a,b) -> (b,a)) $ M.toList dic - objective = tokenizeObj objective_ dic - constraints = tokenizeConstraints constraints_ dic - bounds = tokenizeBounds bounds_ dic - types = tokenizeTypes types_ dic - varRange = (0,M.size dic) - - (varA, varB) = varRange - varCount = varB - varA + 1 - (objsen, obj) = toObj objective - (cnstrs,rhs) = toConstraints constraints varRange - xbnds = toBounds bounds varRange - - types' = V.fromList (replicate varCount CPX_CONTINUOUS) V.// (map (\(a,t) -> (a, typeToCPX t)) types) - - statusLp <- copyMip env lp objsen obj rhs cnstrs (V.fromList xbnds) types' - - case statusLp of - Nothing -> return () - Just msg -> error $ "CPXcopyMIP error: " ++ msg - - case inccb of - Just cb -> do statusIncCB <- setIncumbentCallback env (incumbentcallback cb) - case statusIncCB of - Nothing -> return () - Just msg -> error $ "CPXIncumbentCallBackSet Error: " ++ msg - Nothing -> return () - - case lazycb of - Just cb -> do statusLazyCB <- setLazyConstraintCallback env (lazycallback cb) - case statusLazyCB of - Nothing -> return () - Just msg -> error $ "CPXLazyConstraintCallBackSet Error: " ++ msg - Nothing -> return () - - case cutcb of - Just cb -> do statusCutCB <- setCutCallback env (cutcallback cb) - case statusCutCB of - Nothing -> return () - Just msg -> error $ "CPXCutCallBackSet Error: " ++ msg - Nothing -> return () - - statusOpt <- mipopt env lp - case statusOpt of - Nothing -> return () - Just msg -> error $ "CPXmipopt error: " ++ msg - - statusSol <- getMIPSolution env lp - case statusSol of - Left msg -> error $ "CPXsolution error: " ++ msg - Right sol -> do -- I call this imperative do notation - let vars = V.toList $ VS.convert $ solX sol - let m = M.fromList $ zip (map (revDic M.!) [0..length vars - 1]) vars - return $ MIPSolution (solStat sol == CPXMIP_OPTIMAL) (solObj sol) m - -typeToCPX :: LSolver.Bindings.Type -> CPLEX.Core.Type -typeToCPX (TInteger) = CPX_INTEGER -typeToCPX (TContinuous) = CPX_CONTINUOUS -typeToCPX (TBinary) = CPX_BINARY - - - - -generateVarDic :: (Eq a, Ord a) => Constraints a -> M.Map a Int -generateVarDic (Constraints bounds) = foldr addBoundToDic M.empty bounds - where - addToDic (d :# v) m = M.insertWith (\new old -> old) v (M.size m) m - addBoundToDic (vs :< _ ) m = foldr addToDic m vs - addBoundToDic (vs :> _ ) m = foldr addToDic m vs - addBoundToDic (vs := _ ) m = foldr addToDic m vs - --- Change variables in bound to have ID -tokenizeConstraints :: (Ord a, Eq a) => Constraints a -> M.Map a Int -> Constraints Int -tokenizeConstraints (Constraints bounds) dic = Constraints $ map b2b bounds - where - v2v (d :# v) = d :# (dic M.! v) - b2b (vs :< d) = map v2v vs :< d - b2b (vs :> d) = map v2v vs :> d - b2b (vs := d) = map v2v vs := d - -tokenizeObj :: (Ord a, Eq a) => Optimization a -> M.Map a Int -> Optimization Int -tokenizeObj obj dic = o2o obj - where - v2v (d :# v) = d :# (dic M.! v) - o2o (Minimize vs) = Minimize $ map v2v vs - o2o (Maximize vs) = Maximize $ map v2v vs - -tokenizeBounds :: (Ord a, Eq a) => [(a,b,b)] -> M.Map a Int -> [(Int,b,b)] -tokenizeBounds xs dic = map (\(a,b,c) -> (dic M.! a, b, c) ) xs - -tokenizeTypes :: (Ord a, Eq a) => [(a,b)] -> M.Map a Int -> [(Int, b)] -tokenizeTypes xs dic = map (\(a,b) -> (dic M.! a, b) ) xs diff --git a/src/LSolver/Bindings.hs b/src/LSolver/Bindings.hs deleted file mode 100644 --- a/src/LSolver/Bindings.hs +++ /dev/null @@ -1,70 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} - -module LSolver.Bindings(Variable(..), Bound(..), Constraints(..), Optimization(..), (<+>), - Bounds(..), Type(..), MixedIntegerProblem(..), LinearProblem(..), MIPSolution(..), LPSolution(..)) where - -import Data.List (intercalate) -import qualified Data.Vector as V -import qualified Data.Map as M -import Data.Monoid - -data Variable a = Double :# a - -data Bound x = x :< Double - | x :> Double - | x := Double - deriving Show - -newtype Constraints a = Constraints [ Bound [Variable a] ] - -instance Monoid a => Monoid (Constraints a) where - (Constraints xs) `mappend` (Constraints ys) = Constraints $ xs <> ys - mempty = Constraints [] - -Constraints v1 <+> Constraints v2 = Constraints $ v1 ++ v2 - -data Optimization a = Maximize [Variable a] - | Minimize [Variable a] - -data Type = TContinuous | TInteger | TBinary - -instance (Show a) => Show (Variable a) where - show (d :# v) - | d == (-1) = "-" ++ (show v) - | d == 1 = (show v) - | otherwise = (show d) ++ "x" ++ (show v) - -instance Show a => Show (Optimization a) where - show (Minimize xs) = "Minimize\n\t" ++ (intercalate "+" $ map show xs) - show (Maximize xs) = "Maximize\n\t" ++ (intercalate "+" $ map show xs) - -showVars xs = intercalate " + " $ map show $ zipWith (:#) xs [0..] - -instance (Show a) => Show (Constraints a) where - show (Constraints bounds) = "\nSubject to\n" ++ (unlines $ map (\a -> "\t" ++ a) $ - map getVarSigns bounds) - -printVars xs = intercalate " + " $ map show xs -getVarSigns (x :< v) = (printVars x) ++ " <= " ++ (show v) -getVarSigns (x :> v) = (printVars x) ++ " >= " ++ (show v) -getVarSigns (x := v) = (printVars x) ++ " == " ++ (show v) - -instance Show Type where - show TContinuous = "Continous" - show TInteger = "Integer" - show TBinary = "Binary" - - -type Bounds = [Bound Int] - -data LinearProblem a = LP (Optimization a) (Constraints a) [(a, Maybe Double, Maybe Double)] - deriving Show - -data MixedIntegerProblem a = MILP (Optimization a) (Constraints a) [(a, Maybe Double, Maybe Double)] - [(a,Type)] - deriving Show - -data MIPSolution a = MIPSolution { mipOptimalSol :: Bool, mipObjVal :: Double, mipVars :: M.Map a Double } deriving (Show) - -data LPSolution a = LPSolution { lpOptimalSol :: Bool, lpObjVal :: Double, lpVars :: M.Map a Double, lpDualVars :: V.Vector Double} deriving (Show) - diff --git a/src/LSolver/Dummy.hs b/src/LSolver/Dummy.hs deleted file mode 100644 --- a/src/LSolver/Dummy.hs +++ /dev/null @@ -1,49 +0,0 @@ -module LSolver.Dummy(sol', standardBounds) where - -import Data.Ix as I -import qualified Data.Vector as V -import qualified Data.Sequence as S -import Data.Foldable as F -import LSolver.Bindings - - -data Equality = L Double | E Double | G Double -data Mat = Row Int | Col Int - -data Var = X | Y | Z - -objF :: Optimization Int -objF = Maximize [1 :# 0 , 2 :# 1, 3 :# 2] - -st :: Constraints Int -st = Constraints [ - [(-1):#1, 1:#2, 1:#3] :< 20, - [1:#1, (-3):#2, 1:#3] :< 30 - ] - -bnds = [ (1, Just 0, Just 40), - (2,Just 0,Nothing), - (3, Just 0, Nothing)] - -standardBounds (i,j) = map (\i' -> (i', Just 0, Nothing)) [i..j] - -toBounds bounds varRange = F.toList $ aux bounds def - where - def = S.fromList [k | i <- I.range varRange, let k = (Just 0, Nothing)] - aux [] s = s - aux ((b,lb,ub):bs) s = aux bs (S.update (I.index varRange b) (lb,ub) s) - -toConstraints constraints varRange = let (st, rhs) = toStandard constraints 0 [] [] varRange - in (st, V.fromList rhs) -toStandard (Constraints []) _ accSt accRhs varRange = (reverse $ accSt, reverse $ accRhs) -toStandard (Constraints (b:bs)) rowI accSt accRhs varRange = case b of - vars :< boundVal -> addRow vars L boundVal - vars := boundVal -> addRow vars E boundVal - vars :> boundVal -> addRow vars G boundVal - where addRow vars s boundVal = toStandard (Constraints bs) (rowI+1) (generateRow vars ++ accSt) - ((s boundVal) : accRhs) varRange - generateRow [] = [] - generateRow ((v :# i):vs) = (Row rowI, Col $ I.index varRange i, v):generateRow vs - -sol' = do - putStrLn "Test" diff --git a/src/LSolver/Problems/MinCostMulticom.hs b/src/LSolver/Problems/MinCostMulticom.hs deleted file mode 100644 --- a/src/LSolver/Problems/MinCostMulticom.hs +++ /dev/null @@ -1,88 +0,0 @@ -module LSolver.Problems.MinCostMulticom(generateLinearProblem, getAllPaths, buildAdjacency) where - -import qualified Data.Vector as V -import qualified Data.Map as M -import Data.Vector((!),(//)) -import Control.Monad -import LSolver.Bindings -import LSolver.Backend.Cplex - -generateLinearProblem :: [(Int, Int)] -> [Double] -> [Double] -> [(Int, Int)] -> [Double] - -> LinearProblem Int -generateLinearProblem edges costs capacity commodities demand = - let - n = length costs - varRange = (0,n-1) - adj = buildAdjacency edges - (pathsAll, pathsCom) = genPaths adj commodities - edgeCosts = genEdgeCosts edges costs - edgeCapacity = genEdgeCapacity edges capacity - objFunc = genObjectiveFunction pathsAll edgeCosts - constraints = genConstraints pathsAll pathsCom edges edgeCapacity demand - in LP objFunc constraints (standardBounds varRange) - -loadFromFile :: String -> IO (LinearProblem Int) -loadFromFile fileName = do - contents <- readFile fileName - let [n, e, cost, cap, k, d] = lines $ contents :: [String] - return $ generateLinearProblem (read e) (read cost) (read cap) (read k) (read d) - -genObjectiveFunction pathsAll edgeCosts = Minimize $ zipWith (:#) (V.toList $ (\p -> sum $ - (\a -> edgeCosts M.! a) <$> p) <$> pathsAll) [0..] - -genEdgeCosts edges costs = M.fromList $ zip edges costs -genEdgeCapacity edges capacity = M.fromList $ zip edges capacity -genPaths adj commodities = - let pathsCom = V.fromList $ map V.fromList $ map (\(s,t) -> getAllPaths s t adj) commodities - pathsAll = V.foldr (V.++) V.empty pathsCom - in (pathsAll, pathsCom) - -genConstraints pathsAll pathsCom edges edgeCapacity demand = - Constraints $ (genConstraints1 pathsCom pathsAll) ++ (genConstraints2 pathsAll edges) - where - genConstraints1 pathsCom pathsAll = zipWith (:=) (createOneArrs pathsCom pathsAll) demand - - genConstraints2 pathsAll edges = let pathsEdges = genPathEdges pathsAll edges in - map (\a -> - (map (\p -> 1:#p) $ - pathsEdges M.! a) :< (edgeCapacity M.! a) ) edges --- helper, generates ones and zeroes vector to indicate active path -createOneArrs pathsCom pathsAll = map (\xs -> map (\i -> 1:#i) xs) $ createOneArrs' (V.toList pathsCom) 0 (length pathsAll) -createOneArrs' [] _ _ = [] -createOneArrs' (p:ps) i n = let pn = V.length p in - [i..i+pn-1]:createOneArrs' ps (pn+i) n - - --- Hilfe functions -genPathEdges paths edges = foldr (\p m -> genPathEdges' p edges m) M.empty [0..V.length paths-1] - where - genPathEdges' pid [] m = m - genPathEdges' pid (e:edges) m - | e `elem` (paths ! pid) = genPathEdges' pid edges $ M.insertWith (\new old -> if old /= new then old ++ new else new) e [pid] m - | otherwise = genPathEdges' pid edges m - - -getBounds :: Ord a => [(a,a)] -> (a,a) -getBounds edges = let lst = [fst,snd] <*> edges - in (minimum lst, maximum lst) - - -getAllPaths s e adj = let paths = filter (\xs -> e == last xs) $ filter (\(x:xs) -> x == s) $ getAllPaths' s e adj [] - in (\p -> tail $ zip (0:p) p) <$> paths -getAllPaths' s e _ paths - | s == e = map (\a -> reverse $ e:a) paths -getAllPaths' s e adj paths = - let startedPaths = liftM (s:) ([]:paths) - ps = filter (/= startedPaths) $ (\n -> getAllPaths' n e adj startedPaths) <$> neighbors - in concat $ ps - where neighbors = adj ! s - - ---buildAdjacency :: [(Int, Int)] -> Array (Int, Int) Int -buildAdjacency edges = let n = length edges - in buildAdjacency' edges $ V.replicate n [] -buildAdjacency' [] arr = arr -buildAdjacency' ((i,j):cs) arr = - let cur = (arr ! i) - added = if j `elem` cur then cur else j : cur - in buildAdjacency' cs (arr // [(i,added)])