Haskellで2制約ナップサック

{-- Time-stamp: <2008-10-11 15:24:47>

ナップサック問題の動的計画法と分枝限定法
Usage: ./knapsack2.exe bb < knapsack_problem/problem
       ./knapsack2.exe dp < knapsack_problem/problem

入力データ使用
(アイテム数) (ナップサックの容量)(ナップサックの容量2)
(1) (v_1) (w_1)   (u_1)
(2) (v_2) (w_2)   (u_2)
...
(n) (v_n) (w_n)   (u_n)
--}
import Control.Monad.State (State,get,put,evalState)
import Control.Monad (liftM,replicateM)
import Prelude hiding(lookup)
import Data.Map (Map,lookup,insert,empty)
import Data.List (sortBy)
import System (getArgs)
import System.CPUTime (getCPUTime)
-- ====動的計画法====
type Table k v = Map k v
type Memo a b = State (Table a b) b
type Capacity = (Int,Int)
type Item = (Int,(Int,Int))
-- メモ化のテクニック
memoise :: Ord a => (a -> Memo a b) -> a -> Memo a b
memoise f x = do
  table <- get
  case (lookup x table) of
    Just y -> return y
    Nothing -> do fx <- f x
                  table' <- get
                  put(insert x fx table')
                  return fx
runM :: (a -> Memo a b) -> a -> b
runM m v = evalState (m v) empty
-- 動的計画法本体
knap_dp ::(Capacity,[Item]) -> Memo (Capacity,[Item]) Int
knap_dp ((c1,c2),[])
    | c1 < 0 || c2 < 0 = return $ -1
    | otherwise = return 0
knap_dp ((0,_),_)= return 0
knap_dp ((_,0),_)= return 0
knap_dp ((c1,c2),i:is)
    | c1 < 0 || c2 <0 = return $ -1
    | otherwise = memoise knap_dp' ((c1,c2),i:is)
    where 
      knap_dp' ((c1,c2),i:is) = 
          do x <- knap_dp((c1,c2),is)
             y <- knap_dp((c1-  (fst.snd $ i),c2 - (snd.snd $ i)),is)
             return $ if y < 0 then x  else max x $ y+fst i
run_dp = runM knap_dp
-- ==== 分枝限定法 =====
itemSort = sortBy g
    where g (a,b) (c,d)
              | a*d - b*c > 0 =LT
              | a*d - b*c < 0 =GT
              | otherwise = EQ
itemSort' ((c1,c2),is) = sortBy h is
    where h (v1,(w1,u1)) (v2,(w2,u2))
              |(v1*w2-v2*w1)*c1*u1*u2+(v1*u2-v2*u1)*c2*w1*w2 > 0 =LT
              |(v1*w2-v2*w1)*c1*u1*u2+(v1*u2-v2*u1)*c2*w1*w2 < 0 =GT
              | otherwise = EQ
itemSort'' _ [] = []
itemSort'' f (i:is) =w:(itemSort'' (not f) is')
    where wMax x ys [] = x:ys
          wMax x ys (z:zs)
              | fst x * (g . snd $ z) - fst z * (g . snd $ x) >= 0 = wMax x (z:ys) zs
              | otherwise = wMax z (x:ys) zs
          g = if f then fst else snd
          (w:is') = wMax i [] is
-- 連続緩和問題を解いて、上界を求める
relax ::(Capacity,[Item]) -> Int
relax ((c1,c2),is) = let is1 = itemSort $ map (\(x,y)-> (x,fst y)) is
                         is2 = itemSort $ map (\(x,y)-> (x,snd y)) is
                         is3 = itemSort $ map (\(x,(y,z)) -> (x,y+z)) is
                     in minimum [relax' (c1,is1),relax' (c2,is2),relax' (c1+c2,is3)]
relax' ::(Int,[(Int,Int)])->Int
relax'(_,[])=0
relax'(c,i:is)
    | c < snd i = fst i * c `div` snd i
    | otherwise = relax'(c-snd i,is) + fst i
-- 貪欲法を用いて、下界を求める
greedy ::(Capacity,[Item]) -> Int
greedy (c@(c1,c2),is) =let is1 = itemSort' ((1,0),is)
                           is2 = itemSort' ((0,1),is)
                       in maximum [greedy' snd (c,is1),greedy' snd (c,is2), greedy' itemSort' (c,is),maximum.map fst $ is]
greedy' ::((Capacity,[Item])->[Item]) ->(Capacity,[Item])-> Int 
greedy' _ (_,[])=0
greedy' sort ((c1,c2),i:is)
    | c1 < w || c2 < u = greedy' sort ((c1,c2),is')
    | otherwise = greedy' sort ((c1-w,c2-u),is') + v
    where ((v,(w,u)):is') = sort ((c1,c2),(i:is))
-- 分枝限定法本体
knap_bb ::(Capacity,[Item])->Int->Int->Int
knap_bb ((c1,c2),[]) v s 
    | c1 < 0 || c2 < 0 = s
    | otherwise = max v s
knap_bb (c@(c1,c2),i:is) v s
    | c1 < 0 || c2 < 0 = s
    | ((+v) $! relax (c,i:is)) <= s = s
    | otherwise = let !s1 =knap_bb((c1-(fst.snd$i),c2-(snd.snd$i)),is) (v+ fst i) s
                      !s2 =knap_bb(c,is) v s1
                  in max s1 s2
run_bb ((c1,c2),is) = knap_bb ((c1,c2), iSort is') 0 v
    where is' = filter (\i-> (fst.snd $ i) <= c1 && (snd.snd$i) <=c2) is
          v = greedy ((c1,c2),is')
          iSort | c1 < c2 = itemSort'' True
                | otherwise = itemSort'' False
-- テストデータ
items ::[Item]
items =[(10,(2,5)),(6,(19,2)),(29,(13,3)),(12,(7,9))]
cap :: Capacity
cap =(19,20)
ins =(cap,items)
-- 入出力
getInts = liftM (map read .words) getLine :: IO [Int]
getItem = liftM (toTriple . tail . map read . words) getLine :: IO Item
    where toTriple (x:y:z:ws) = (x,(y,z))
main :: IO ()
main = do
  start <- getCPUTime
  arg <- getArgs
  [num,c1,c2] <- getInts
  is <- replicateM num getItem
  case head arg of
    "bb" -> putStr $ "opt:"++(show.run_bb) ((c1,c2),is)
    "dp" -> putStr $ "opt:"++ (show.run_dp) ((c1,c2),is)
  end <- getCPUTime
  putStr $ ",low:" ++ (show.greedy) ((c1,c2),is) ++",up:"++(show.relax) ((c1,c2),is)
  putStrLn $ ",time:" ++ show  ((end -start) `div` 10^9) ++ "ms"