{-- 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 <  || c2 <  = return $ -1
| otherwise = return 
knap_dp ((,_),_)= return 
knap_dp ((_,),_)= return 
knap_dp ((c1,c2),i:is)
| c1 <  || c2 < = 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 <  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 >  =LT
| a*d - b*c <  =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 >  =LT
|(v1*w2-v2*w1)*c1*u1*u2+(v1*u2-v2*u1)*c2*w1*w2 <  =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) >=  = 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'(_,[])=
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,),is)
is2 = itemSort' ((,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' _ (_,[])=
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 <  || c2 <  = s
| otherwise = max v s
knap_bb (c@(c1,c2),i:is) v s
| c1 <  || c2 <  = 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')  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"