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.2.0.2 +version: 0.3.0.0 synopsis: high-level CPLEX interface -- description: License: BSD3 @@ -15,18 +15,23 @@ High level interface to CPLEX build-type: Simple -cabal-version: >=1.8 +cabal-version: >=1.24 Library + default-language: Haskell2010 exposed-modules: CPLEX.Core, CPLEX.Param CPLEX.Bindings + LSolver.Bindings + LSolver.Dummy + LSolver.Backend.Cplex + LSolver.Problems.MinCostMulticom -- other-modules: - Build-depends: base >= 4 && < 5 - , containers - , mtl - , primitive >= 0.6.1.0 - , transformers >= 0.4.2.0 - , vector + Build-depends: base < 5.0 + , containers < 6.0 + , mtl < 2.3 + , primitive < 0.7.0.0 + , transformers < 0.5.3.0 + , vector < 0.12.0.0 hs-source-dirs: src if os(linux) diff --git a/src/CPLEX/Core.hs b/src/CPLEX/Core.hs --- a/src/CPLEX/Core.hs +++ b/src/CPLEX/Core.hs @@ -31,6 +31,7 @@ , hybnetopt , getSolution , getMIPSolution + , writeprob -- * change things , changeCoefList , changeObj @@ -332,6 +333,12 @@ k -> fmap Left (getErrorString env (CpxRet k)) k -> fmap Left (getErrorString env (CpxRet k)) +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 + toCpxError :: CpxEnv -> CInt -> IO (Maybe String) toCpxError env 0 = return Nothing toCpxError env k = do @@ -417,6 +424,44 @@ preorder = map (\(row,col,val) -> (col, [(row, val)])) amat +toRowForm :: Int -> [(Row,Col,Double)] -> (Vector CInt, Vector CInt, Vector CInt, Vector CDouble) +toRowForm numrows amat = (matbeg, matcnt, matind, matval) + where + matbeg = VS.fromList $ map fromIntegral begs + matcnt = VS.fromList $ map fromIntegral cnts + matind = VS.fromList $ map (fromIntegral . unCol) inds + matval = VS.fromList $ map realToFrac vals + + -- sort colMap into the from CPLEX wants + inds :: [Col] + vals :: [Double] + (inds,vals) = unzip $ concat cols + + begs :: [Int] + cnts :: [Int] + cols :: [[(Col,Double)]] + (begs,cnts,cols) = unzip3 $ rowMapInfo' 0 $ M.elems rowMap + + rowMapInfo' :: Int -> [[(Col,Double)]] -> [(Int,Int,[(Col,Double)])] + rowMapInfo' beg (col:xs) = (beg,cnt,col) : rowMapInfo' (beg + cnt) xs + where + cnt = length col + rowMapInfo' _ [] = [] + + -- add Rows with no entries in case some are missing + rowMap = M.union rowMap' emptyRowMap + + emptyRowMap :: M.Map Row [(Col,Double)] + emptyRowMap = M.fromList $ take numrows $ zip (map Row [0..]) (repeat []) + + -- a map from Row to all (Col,Double) pairs + rowMap' :: M.Map Row [(Col,Double)] + rowMap' = M.fromListWith (++) preorder + + -- reorganize the (Row,Col,Double) into (Row, [(Col,Double)]) with only 1 (Row,Double) + preorder :: [(Row,[(Col,Double)])] + preorder = map (\(row,col,val) -> (row, [(col, val)])) amat + copyLp :: CpxEnv -> CpxLp -> ObjSense -> V.Vector Double -> V.Vector Sense -> [(Row,Col,Double)] -> V.Vector (Maybe Double, Maybe Double) -> IO (Maybe String) copyLp = copyLpWithFun' c_CPXcopylp @@ -778,7 +823,7 @@ sense = VS.fromList $ V.toList sense' rhs = VS.fromList $ V.toList rhs' - (matbeg, matcnt, matind, matval) = toColForm nzcnt aMat + (matbeg, matcnt, matind, matval) = toRowForm rcnt aMat getErrorStatus env status = case status of 0 -> return Nothing diff --git a/src/LSolver/Backend/Cplex.hs b/src/LSolver/Backend/Cplex.hs new file mode 100644 --- /dev/null +++ b/src/LSolver/Backend/Cplex.hs @@ -0,0 +1,287 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/src/LSolver/Bindings.hs @@ -0,0 +1,70 @@ +{-# 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 new file mode 100644 --- /dev/null +++ b/src/LSolver/Dummy.hs @@ -0,0 +1,49 @@ +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 new file mode 100644 --- /dev/null +++ b/src/LSolver/Problems/MinCostMulticom.hs @@ -0,0 +1,88 @@ +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)])