msakai / toysolver / 513

Committed 24 Nov 2024 08:49AM UTC coverage: 68.675% (-1.0%) from 69.681%

Build # 513

Build Type

push

github

Committed by

web-flow

Commit Message

Merge 496ed4263 into 46e1509c4

Run Details

5 of 121 new or added lines in 7 files covered. (4.13%)

105 existing lines in 14 files now uncovered.

9745 of 14190 relevant lines covered (68.68%)

0.69 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

59.9

/src/ToySolver/Combinatorial/SubsetSum.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Combinatorial.SubsetSum
-- Copyright   :  (c) Masahiro Sakai 2015
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- References
--
-- * D. Pisinger, "An exact algorithm for large multiple knapsack problems,"
--   European Journal of Operational Research, vol. 114, no. 3, pp. 528-541,
--   May 1999. DOI:10.1016/s0377-2217(98)00120-9
--   <http://www.sciencedirect.com/science/article/pii/S0377221798001209>
--   <http://www.diku.dk/~pisinger/95-6.ps>
--
-----------------------------------------------------------------------------
module ToySolver.Combinatorial.SubsetSum
  ( Weight
  , subsetSum
  , maxSubsetSum
  , minSubsetSum
  ) where

import Control.Exception (assert)
import Control.Monad
import Control.Monad.ST
import Data.STRef
import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IntMap
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map
import Data.Vector.Generic ((!))
import qualified Data.Vector as V
import qualified Data.Vector.Generic as VG
import qualified Data.Vector.Generic.Mutable as VM
import qualified Data.Vector.Unboxed as VU

type Weight = Integer

-- | Maximize Σ_{i=1}^n wi xi subject to Σ_{i=1}^n wi xi ≤ c and xi ∈ {0,1}.
--
-- Note: 0 (resp. 1) is identified with False (resp. True) in the assignment.
maxSubsetSum
  :: VG.Vector v Weight
  => v Weight -- ^ weights @[w1, w2 .. wn]@
  -> Weight -- ^ capacity @c@
  -> Maybe (Weight, VU.Vector Bool)
  -- ^
  -- * the objective value Σ_{i=1}^n wi xi, and
  --
  -- * the assignment @[x1, x2 .. xn]@, identifying 0 (resp. 1) with @False@ (resp. @True@).
maxSubsetSum w c =
  case normalizeWeightsToPositive (w,c) of
    (w1, c1, trans1)
      | c1 < 0 -> Nothing
      | otherwise ->
          case normalize2 (w1, c1) of
            (w2, c2, trans2) ->
              case normalizeGCDLe (w2, c2) of
                (w3, c3, trans3) ->
                  Just $ trans1 $ trans2 $ trans3 $ maxSubsetSum' w3 c3

normalizeWeightsToPositive
  :: VG.Vector v Weight
  => (v Weight, Weight)
  -> (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
normalizeWeightsToPositive (w,c)
  | VG.all (>=0) w = (VG.convert w, c, id)
  | otherwise = runST $ do
      w2 <- VM.new (VG.length w)
      let loop !i !offset
            | i >= VG.length w = return offset
            | otherwise = do
                let wi = w ! i
                if wi < 0 then do
                  VM.write w2 i (- wi)
                  loop (i+1) (offset + wi)
                else do
                  VM.write w2 i wi
                  loop (i+1) offset
      offset <- loop 0 (0::Integer)
      w2' <- VG.unsafeFreeze w2
      let trans (obj, bs) = (obj + offset, bs2)
            where
              bs2 = VU.imap (\i bi -> if w ! i < 0 then not bi else bi) bs
      return (w2', c - offset, trans)

normalize2
  :: (V.Vector Weight, Weight)
  -> (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
normalize2 (w,c)
  | VG.all (\wi -> 0<wi && wi<=c) w = (w, c, id)
  | otherwise = (VG.filter (\wi -> 0<wi && wi<=c) w, c, trans)
  where
    trans (obj, bs) = (obj, bs2)
      where
        bs2 = VU.create $ do
          v <- VM.new (VG.length w)
          let loop !i !j =
                when (i < VG.length w) $ do
                  let wi = w ! i
                  if 0 < wi && wi <= c then do
                    VM.write v i (bs ! j)
                    loop (i+1) (j+1)
                  else do
                    VM.write v i False
                    loop (i+1) j
          loop 0 0
          return v

normalizeGCDLe
  :: (V.Vector Weight, Weight)
  -> (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
normalizeGCDLe (w,c)
  | VG.null w || d == 1 = (w, c, id)
  | otherwise = (VG.map (`div` d) w, c `div` d, trans)
  where
    d = VG.foldl1' gcd w
    trans (obj, bs) = (obj * d, bs)

normalizeGCDEq
  :: (V.Vector Weight, Weight)
  -> Maybe (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
normalizeGCDEq (w,c)
  | VG.null w || d == 1 = Just (w, c, id)
  | c `mod` d == 0 = Just (VG.map (`div` d) w, c `div` d, trans)
  | otherwise = Nothing
  where
    d = VG.foldl1' gcd w
    trans (obj, bs) = (obj * d, bs)

maxSubsetSum' :: V.Vector Weight -> Weight -> (Weight, VU.Vector Bool)
maxSubsetSum' !w !c
  | wsum <= c = (wsum, VG.replicate (VG.length w) True)
  | c <= fromIntegral (maxBound :: Int) =
      maxSubsetSumInt' (VG.generate (VG.length w) (\i -> fromIntegral (w VG.! i))) (fromIntegral c) wsum
  | otherwise =
      maxSubsetSumInteger' w c wsum
  where
    wsum = VG.sum w

maxSubsetSumInteger' :: V.Vector Weight -> Weight -> Weight -> (Weight, VU.Vector Bool)
maxSubsetSumInteger' w !c wsum = assert (wbar <= c) $ assert (wbar + (w ! b) > c) $ runST $ do
  objRef <- newSTRef (wbar, [], [])
  let updateObj gs ft = do
        let loop [] _ = return ()
            loop _ [] = return ()
            loop xxs@((gobj,gsol):xs) yys@((fobj,fsol):ys)
              | c < gobj + fobj = loop xs yys
              | otherwise = do
                  (curr, _, _) <- readSTRef objRef
                  when (curr < gobj + fobj) $ writeSTRef objRef (gobj + fobj, gsol, fsol)
                  loop xxs ys
        loop (Map.toDescList gs) (Map.toAscList ft)

  let loop !s !t !gs !ft !flag = do
        (obj, gsol, fsol) <- readSTRef objRef
        if obj == c || (s == 0 && t == n-1) then do
          let sol = VG.create $ do
                bs <- VM.new n
                forM_ [0..b-1] $ \i -> VM.write bs i True
                forM_ [b..n-1] $ \i -> VM.write bs i False
                forM_ fsol $ \i -> VM.write bs i True
                forM_ gsol $ \i -> VM.write bs i False
                return bs
          return (obj, sol)
        else do
          let updateF = do
                -- Compute f_{t+1} from f_t
                let t' = t + 1
                    wt' = w ! t'
                    m = Map.mapKeysMonotonic (+ wt') $ Map.map (t' :) $ splitLE (c - wt') ft
                    ft' = ft `Map.union` m
                updateObj gs m
                loop s t' gs ft' (not flag)
              updateG = do
                -- Compute g_{s-1} from g_s
                let s' = s - 1
                    ws = w ! s'
                    m = Map.map (s' :) $ g_drop $ Map.mapKeysMonotonic (subtract ws) $ gs
                    gs' = gs `Map.union` m
                updateObj m ft
                loop s' t gs' ft (not flag)
          if s == 0 then
            updateF
          else if t == n-1 then
            updateG
          else
            if flag then updateG else updateF

  let -- f_{b-1}
      fb' :: Map Integer [Int]
      fb' = Map.singleton 0 []
      -- g_{b}
      gb :: Map Integer [Int]
      gb = Map.singleton wbar []
  loop b (b-1) gb fb' True

  where
    n = VG.length w

    b :: Int
    b = loop (-1) 0
      where
        loop :: Int -> Integer -> Int
        loop !i !s
          | s > c = i
          | otherwise = loop (i+1) (s + (w ! (i+1)))

    wbar :: Weight
    wbar = VG.sum $ VG.slice 0 b w

    max_f :: Weight
    max_f = wsum - fromIntegral wbar

    min_g :: Weight
    min_g = 0 `max` (c - max_f)

    g_drop :: Map Integer [Int] -> Map Integer [Int]
    g_drop g =
      case Map.splitLookup min_g g of
        (lo, _, _) | Map.null lo -> g
        (_, Just v, hi) -> Map.insert min_g v hi
        (lo, Nothing, hi) ->
          case Map.findMax lo of
            (k,v) -> Map.insert k v hi

    splitLE :: Ord k => k -> Map k v -> Map k v
    splitLE k m =
      case Map.splitLookup k m of
        (lo, Nothing, _) -> lo
        (lo, Just v, _) -> Map.insert k v lo

maxSubsetSumInt' :: VU.Vector Int -> Int -> Weight -> (Weight, VU.Vector Bool)
maxSubsetSumInt' w !c wsum = assert (wbar <= c) $ assert (wbar + (w ! b) > c) $ runST $ do
  objRef <- newSTRef (wbar, [], [])
  let updateObj gs ft = do
        let loop [] _ = return ()
            loop _ [] = return ()
            loop xxs@((gobj,gsol):xs) yys@((fobj,fsol):ys)
              | c < gobj + fobj = loop xs yys
              | otherwise = do
                  (curr, _, _) <- readSTRef objRef
                  when (curr < gobj + fobj) $ writeSTRef objRef (gobj + fobj, gsol, fsol)
                  loop xxs ys
        loop (IntMap.toDescList gs) (IntMap.toAscList ft)

  let loop !s !t !gs !ft !flag = do
        (obj, gsol, fsol) <- readSTRef objRef
        if obj == c || (s == 0 && t == n-1) then do
          let sol = VG.create $ do
                bs <- VM.new n
                forM_ [0..b-1] $ \i -> VM.write bs i True
                forM_ [b..n-1] $ \i -> VM.write bs i False
                forM_ fsol $ \i -> VM.write bs i True
                forM_ gsol $ \i -> VM.write bs i False
                return bs
          return (fromIntegral obj, sol)
        else do
          let updateF = do
                -- Compute f_{t+1} from f_t
                let t' = t + 1
                    wt' = w ! t'
                    m = IntMap.mapKeysMonotonic (+ wt') $ IntMap.map (t' :) $ splitLE (c - wt') ft
                    ft' = ft `IntMap.union` m
                updateObj gs m
                loop s t' gs ft' (not flag)
              updateG = do
                -- Compute g_{s-1} from g_s
                let s' = s - 1
                    ws = w ! s'
                    m = IntMap.map (s' :) $ g_drop $ IntMap.mapKeysMonotonic (subtract ws) $ gs
                    gs' = gs `IntMap.union` m
                updateObj m ft
                loop s' t gs' ft (not flag)
          if s == 0 then
            updateF
          else if t == n-1 then
            updateG
          else
            if flag then updateG else updateF

  let -- f_{b-1}
      fb' :: IntMap [Int]
      fb' = IntMap.singleton 0 []
      -- g_{b}
      gb :: IntMap [Int]
      gb = IntMap.singleton wbar []
  loop b (b-1) gb fb' True

  where
    n = VG.length w

    b :: Int
    b = loop (-1) 0
      where
        loop :: Int -> Integer -> Int
        loop !i !s
          | s > fromIntegral c = i
          | otherwise = loop (i+1) (s + fromIntegral (w ! (i+1)))

    wbar :: Int
    wbar = VG.sum $ VG.slice 0 b w

    max_f :: Integer
    max_f = wsum - fromIntegral wbar

    min_g :: Int
    min_g = if max_f < fromIntegral c then c - fromIntegral max_f else 0

    g_drop :: IntMap [Int] -> IntMap [Int]
    g_drop g =
      case IntMap.splitLookup min_g g of
        (lo, _, _) | IntMap.null lo -> g
        (_, Just v, hi) -> IntMap.insert min_g v hi
        (lo, Nothing, hi) ->
          case IntMap.findMax lo of
            (k,v) -> IntMap.insert k v hi

    splitLE :: Int -> IntMap v -> IntMap v
    splitLE k m =
      case IntMap.splitLookup k m of
        (lo, Nothing, _) -> lo
        (lo, Just v, _) -> IntMap.insert k v lo

-- | Minimize Σ_{i=1}^n wi xi subject to Σ_{i=1}^n wi xi ≥ l and xi ∈ {0,1}.
--
-- Note: 0 (resp. 1) is identified with False (resp. True) in the assignment.
minSubsetSum
  :: VG.Vector v Weight
  => v Weight -- ^ weights @[w1, w2 .. wn]@
  -> Weight -- ^ @l@
  -> Maybe (Weight, VU.Vector Bool)
  -- ^
  -- * the objective value Σ_{i=1}^n wi xi, and
  --
  -- * the assignment @[x1, x2 .. xn]@, identifying 0 (resp. 1) with @False@ (resp. @True@).
minSubsetSum w l =
  case maxSubsetSum w (wsum - l) of
    Nothing -> Nothing
    Just (obj, bs) -> Just (wsum - obj, VG.map not bs)
  where
    wsum = VG.sum w

{-
minimize Σ wi xi = Σ wi (1 - ¬xi) = Σ wi - (Σ wi ¬xi)
subject to Σ wi xi ≥ n

maximize Σ wi ¬xi
subject to Σ wi ¬xi ≤ (Σ wi) - n

Σ wi xi ≥ n
Σ wi (1 - ¬xi) ≥ n
(Σ wi) - (Σ wi ¬xi) ≥ n
(Σ wi ¬xi) ≤ (Σ wi) - n
-}

-- | Solve Σ_{i=1}^n wi xi = c and xi ∈ {0,1}.
--
-- Note that this is different from usual definition of the subset sum problem,
-- as this definition allows all xi to be zero.
--
-- Note: 0 (resp. 1) is identified with False (resp. True) in the assignment.
subsetSum
  :: VG.Vector v Weight
  => v Weight -- ^ weights @[w1, w2 .. wn]@
  -> Weight -- ^ @l@
  -> Maybe (VU.Vector Bool)
  -- ^
  -- the assignment @[x1, x2 .. xn]@, identifying 0 (resp. 1) with @False@ (resp. @True@).
subsetSum w c =
  case normalizeWeightsToPositive (w,c) of
    (w1, c1, trans1)
      | c1 < 0 -> Nothing
      | otherwise ->
          case normalize2 (w1, c1) of
            (w2, c2, trans2) -> do
              (w3, c3, trans3) <- normalizeGCDEq (w2,c2)
              let (obj, sol) = maxSubsetSum' w3 c3
              guard $ obj == c3
              return $ snd $ trans1 $ trans2 $ trans3 (obj, sol)

1	{-# LANGUAGE BangPatterns #-}
2	{-# LANGUAGE FlexibleContexts #-}
3	{-# LANGUAGE ScopedTypeVariables #-}
4	{-# OPTIONS_GHC -Wall #-}
5	{-# OPTIONS_HADDOCK show-extensions #-}
6	-----------------------------------------------------------------------------
7	-- \|
8	-- Module : ToySolver.Combinatorial.SubsetSum
9	-- Copyright : (c) Masahiro Sakai 2015
10	-- License : BSD-style
11	--
12	-- Maintainer : masahiro.sakai@gmail.com
13	-- Stability : provisional
14	-- Portability : non-portable
15	--
16	-- References
17	--
18	-- * D. Pisinger, "An exact algorithm for large multiple knapsack problems,"
19	-- European Journal of Operational Research, vol. 114, no. 3, pp. 528-541,
20	-- May 1999. DOI:10.1016/s0377-2217(98)00120-9
21	-- <http://www.sciencedirect.com/science/article/pii/S0377221798001209>
22	-- <http://www.diku.dk/~pisinger/95-6.ps>
23	--
24	-----------------------------------------------------------------------------
25	module ToySolver.Combinatorial.SubsetSum
26	( Weight
27	, subsetSum
28	, maxSubsetSum
29	, minSubsetSum
30	) where
31
32	import Control.Exception (assert)
33	import Control.Monad
34	import Control.Monad.ST
35	import Data.STRef
36	import Data.IntMap.Strict (IntMap)
37	import qualified Data.IntMap.Strict as IntMap
38	import Data.Map.Strict (Map)
39	import qualified Data.Map.Strict as Map
40	import Data.Vector.Generic ((!))
41	import qualified Data.Vector as V
42	import qualified Data.Vector.Generic as VG
43	import qualified Data.Vector.Generic.Mutable as VM
44	import qualified Data.Vector.Unboxed as VU
45
46	type Weight = Integer
47
48	-- \| Maximize Σ_{i=1}^n wi xi subject to Σ_{i=1}^n wi xi ≤ c and xi ∈ {0,1}.
49	--
50	-- Note: 0 (resp. 1) is identified with False (resp. True) in the assignment.
51	maxSubsetSum
52	:: VG.Vector v Weight
53	=> v Weight -- ^ weights @[w1, w2 .. wn]@
54	-> Weight -- ^ capacity @c@
55	-> Maybe (Weight, VU.Vector Bool)
56	-- ^
57	-- * the objective value Σ_{i=1}^n wi xi, and
58	--
59	-- * the assignment @[x1, x2 .. xn]@, identifying 0 (resp. 1) with @False@ (resp. @True@).
60	maxSubsetSum w c =	1✔
61	case normalizeWeightsToPositive (w,c) of	1✔
62	(w1, c1, trans1)
63	\| c1 < 0 -> Nothing	1✔
64	\| otherwise ->	×
65	case normalize2 (w1, c1) of	1✔
66	(w2, c2, trans2) ->
67	case normalizeGCDLe (w2, c2) of	1✔
68	(w3, c3, trans3) ->
69	Just $ trans1 $ trans2 $ trans3 $ maxSubsetSum' w3 c3	1✔
70
71	normalizeWeightsToPositive
72	:: VG.Vector v Weight
73	=> (v Weight, Weight)
74	-> (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
75	normalizeWeightsToPositive (w,c)	1✔
76	\| VG.all (>=0) w = (VG.convert w, c, id)	1✔
77	\| otherwise = runST $ do	×
78	w2 <- VM.new (VG.length w)	1✔
79	let loop !i !offset	1✔
80	\| i >= VG.length w = return offset	1✔
81	\| otherwise = do	×
82	let wi = w ! i	1✔
83	if wi < 0 then do	1✔
84	VM.write w2 i (- wi)	1✔
85	loop (i+1) (offset + wi)	1✔
86	else do	1✔
87	VM.write w2 i wi	1✔
88	loop (i+1) offset	1✔
89	offset <- loop 0 (0::Integer)	1✔
90	w2' <- VG.unsafeFreeze w2	1✔
91	let trans (obj, bs) = (obj + offset, bs2)	1✔
92	where
93	bs2 = VU.imap (\i bi -> if w ! i < 0 then not bi else bi) bs	1✔
94	return (w2', c - offset, trans)	1✔
95
96	normalize2
97	:: (V.Vector Weight, Weight)
98	-> (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
99	normalize2 (w,c)	1✔
100	\| VG.all (\wi -> 0<wi && wi<=c) w = (w, c, id)	1✔
101	\| otherwise = (VG.filter (\wi -> 0<wi && wi<=c) w, c, trans)	×
102	where
103	trans (obj, bs) = (obj, bs2)	1✔
104	where
105	bs2 = VU.create $ do	1✔
106	v <- VM.new (VG.length w)	1✔
107	let loop !i !j =	1✔
108	when (i < VG.length w) $ do	1✔
109	let wi = w ! i	1✔
110	if 0 < wi && wi <= c then do	1✔
111	VM.write v i (bs ! j)	1✔
112	loop (i+1) (j+1)	1✔
113	else do	1✔
114	VM.write v i False	1✔
115	loop (i+1) j	1✔
116	loop 0 0	1✔
117	return v	1✔
118
119	normalizeGCDLe
120	:: (V.Vector Weight, Weight)
121	-> (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
122	normalizeGCDLe (w,c)	1✔
123	\| VG.null w \|\| d == 1 = (w, c, id)	1✔
124	\| otherwise = (VG.map (`div` d) w, c `div` d, trans)	×
125	where
126	d = VG.foldl1' gcd w	1✔
127	trans (obj, bs) = (obj * d, bs)	1✔
128
129	normalizeGCDEq
130	:: (V.Vector Weight, Weight)
131	-> Maybe (V.Vector Weight, Weight, (Weight, VU.Vector Bool) -> (Weight, VU.Vector Bool))
132	normalizeGCDEq (w,c)	1✔
133	\| VG.null w \|\| d == 1 = Just (w, c, id)	1✔
UNCOV 134	\| c `mod` d == 0 = Just (VG.map (`div` d) w, c `div` d, trans)	×
135	\| otherwise = Nothing	×
136	where
137	d = VG.foldl1' gcd w	1✔
138	trans (obj, bs) = (obj * d, bs)	×
139
140	maxSubsetSum' :: V.Vector Weight -> Weight -> (Weight, VU.Vector Bool)
141	maxSubsetSum' !w !c	1✔
142	\| wsum <= c = (wsum, VG.replicate (VG.length w) True)	1✔
143	\| c <= fromIntegral (maxBound :: Int) =	×
144	maxSubsetSumInt' (VG.generate (VG.length w) (\i -> fromIntegral (w VG.! i))) (fromIntegral c) wsum	1✔
145	\| otherwise =	×
146	maxSubsetSumInteger' w c wsum	×
147	where
148	wsum = VG.sum w	1✔
149
150	maxSubsetSumInteger' :: V.Vector Weight -> Weight -> Weight -> (Weight, VU.Vector Bool)
151	maxSubsetSumInteger' w !c wsum = assert (wbar <= c) $ assert (wbar + (w ! b) > c) $ runST $ do	×
152	objRef <- newSTRef (wbar, [], [])	×
153	let updateObj gs ft = do	×
154	let loop [] _ = return ()	×
155	loop _ [] = return ()	×
156	loop xxs@((gobj,gsol):xs) yys@((fobj,fsol):ys)
157	\| c < gobj + fobj = loop xs yys	×
158	\| otherwise = do	×
159	(curr, _, _) <- readSTRef objRef	×
160	when (curr < gobj + fobj) $ writeSTRef objRef (gobj + fobj, gsol, fsol)	×
161	loop xxs ys	×
162	loop (Map.toDescList gs) (Map.toAscList ft)	×
163
164	let loop !s !t !gs !ft !flag = do	×
165	(obj, gsol, fsol) <- readSTRef objRef	×
166	if obj == c \|\| (s == 0 && t == n-1) then do	×
167	let sol = VG.create $ do	×
168	bs <- VM.new n	×
169	forM_ [0..b-1] $ \i -> VM.write bs i True	×
170	forM_ [b..n-1] $ \i -> VM.write bs i False	×
171	forM_ fsol $ \i -> VM.write bs i True	×
172	forM_ gsol $ \i -> VM.write bs i False	×
173	return bs	×
174	return (obj, sol)	×
175	else do	×
176	let updateF = do	×
177	-- Compute f_{t+1} from f_t
178	let t' = t + 1	×
179	wt' = w ! t'	×
180	m = Map.mapKeysMonotonic (+ wt') $ Map.map (t' :) $ splitLE (c - wt') ft	×
181	ft' = ft `Map.union` m	×
182	updateObj gs m	×
183	loop s t' gs ft' (not flag)	×
184	updateG = do	×
185	-- Compute g_{s-1} from g_s
186	let s' = s - 1	×
187	ws = w ! s'	×
188	m = Map.map (s' :) $ g_drop $ Map.mapKeysMonotonic (subtract ws) $ gs	×
189	gs' = gs `Map.union` m	×
190	updateObj m ft	×
191	loop s' t gs' ft (not flag)	×
192	if s == 0 then	×
193	updateF	×
194	else if t == n-1 then	×
195	updateG	×
196	else
197	if flag then updateG else updateF	×
198
199	let -- f_{b-1}
200	fb' :: Map Integer [Int]
201	fb' = Map.singleton 0 []	×
202	-- g_{b}
203	gb :: Map Integer [Int]
204	gb = Map.singleton wbar []	×
205	loop b (b-1) gb fb' True	×
206
207	where
208	n = VG.length w	×
209
210	b :: Int
211	b = loop (-1) 0	×
212	where
213	loop :: Int -> Integer -> Int
214	loop !i !s	×
215	\| s > c = i	×
216	\| otherwise = loop (i+1) (s + (w ! (i+1)))	×
217
218	wbar :: Weight
219	wbar = VG.sum $ VG.slice 0 b w	×
220
221	max_f :: Weight
222	max_f = wsum - fromIntegral wbar	×
223
224	min_g :: Weight
225	min_g = 0 `max` (c - max_f)	×
226
227	g_drop :: Map Integer [Int] -> Map Integer [Int]
228	g_drop g =	×
229	case Map.splitLookup min_g g of	×
230	(lo, _, _) \| Map.null lo -> g	×
231	(_, Just v, hi) -> Map.insert min_g v hi	×
232	(lo, Nothing, hi) ->
233	case Map.findMax lo of	×
234	(k,v) -> Map.insert k v hi	×
235
236	splitLE :: Ord k => k -> Map k v -> Map k v
237	splitLE k m =	×
238	case Map.splitLookup k m of	×
239	(lo, Nothing, _) -> lo	×
240	(lo, Just v, _) -> Map.insert k v lo	×
241
242	maxSubsetSumInt' :: VU.Vector Int -> Int -> Weight -> (Weight, VU.Vector Bool)
243	maxSubsetSumInt' w !c wsum = assert (wbar <= c) $ assert (wbar + (w ! b) > c) $ runST $ do	×
244	objRef <- newSTRef (wbar, [], [])	1✔
245	let updateObj gs ft = do	1✔
246	let loop [] _ = return ()	×
247	loop _ [] = return ()	×
248	loop xxs@((gobj,gsol):xs) yys@((fobj,fsol):ys)
249	\| c < gobj + fobj = loop xs yys	1✔
250	\| otherwise = do	×
251	(curr, _, _) <- readSTRef objRef	1✔
252	when (curr < gobj + fobj) $ writeSTRef objRef (gobj + fobj, gsol, fsol)	1✔
253	loop xxs ys	1✔
254	loop (IntMap.toDescList gs) (IntMap.toAscList ft)	1✔
255
256	let loop !s !t !gs !ft !flag = do	1✔
257	(obj, gsol, fsol) <- readSTRef objRef	1✔
258	if obj == c \|\| (s == 0 && t == n-1) then do	1✔
259	let sol = VG.create $ do	1✔
260	bs <- VM.new n	1✔
261	forM_ [0..b-1] $ \i -> VM.write bs i True	1✔
262	forM_ [b..n-1] $ \i -> VM.write bs i False	1✔
263	forM_ fsol $ \i -> VM.write bs i True	1✔
264	forM_ gsol $ \i -> VM.write bs i False	1✔
265	return bs	1✔
266	return (fromIntegral obj, sol)	1✔
267	else do	1✔
268	let updateF = do	1✔
269	-- Compute f_{t+1} from f_t
270	let t' = t + 1	1✔
271	wt' = w ! t'	1✔
272	m = IntMap.mapKeysMonotonic (+ wt') $ IntMap.map (t' :) $ splitLE (c - wt') ft	1✔
273	ft' = ft `IntMap.union` m	1✔
274	updateObj gs m	1✔
275	loop s t' gs ft' (not flag)	1✔
276	updateG = do	1✔
277	-- Compute g_{s-1} from g_s
278	let s' = s - 1	1✔
279	ws = w ! s'	1✔
280	m = IntMap.map (s' :) $ g_drop $ IntMap.mapKeysMonotonic (subtract ws) $ gs	1✔
281	gs' = gs `IntMap.union` m	1✔
282	updateObj m ft	1✔
283	loop s' t gs' ft (not flag)	1✔
284	if s == 0 then	1✔
285	updateF	1✔
286	else if t == n-1 then	1✔
287	updateG	1✔
288	else
289	if flag then updateG else updateF	1✔
290
291	let -- f_{b-1}
292	fb' :: IntMap [Int]
293	fb' = IntMap.singleton 0 []	1✔
294	-- g_{b}
295	gb :: IntMap [Int]
296	gb = IntMap.singleton wbar []	1✔
297	loop b (b-1) gb fb' True	1✔
298
299	where
300	n = VG.length w	1✔
301
302	b :: Int
303	b = loop (-1) 0	1✔
304	where
305	loop :: Int -> Integer -> Int
306	loop !i !s	1✔
307	\| s > fromIntegral c = i	1✔
308	\| otherwise = loop (i+1) (s + fromIntegral (w ! (i+1)))	×
309
310	wbar :: Int
311	wbar = VG.sum $ VG.slice 0 b w	1✔
312
313	max_f :: Integer
314	max_f = wsum - fromIntegral wbar	1✔
315
316	min_g :: Int
317	min_g = if max_f < fromIntegral c then c - fromIntegral max_f else 0	1✔
318
319	g_drop :: IntMap [Int] -> IntMap [Int]
320	g_drop g =	1✔
321	case IntMap.splitLookup min_g g of	1✔
322	(lo, _, _) \| IntMap.null lo -> g	1✔
323	(_, Just v, hi) -> IntMap.insert min_g v hi	1✔
324	(lo, Nothing, hi) ->
325	case IntMap.findMax lo of	1✔
326	(k,v) -> IntMap.insert k v hi	1✔
327
328	splitLE :: Int -> IntMap v -> IntMap v
329	splitLE k m =	1✔
330	case IntMap.splitLookup k m of	1✔
331	(lo, Nothing, _) -> lo	1✔
332	(lo, Just v, _) -> IntMap.insert k v lo	1✔
333
334	-- \| Minimize Σ_{i=1}^n wi xi subject to Σ_{i=1}^n wi xi ≥ l and xi ∈ {0,1}.
335	--
336	-- Note: 0 (resp. 1) is identified with False (resp. True) in the assignment.
337	minSubsetSum
338	:: VG.Vector v Weight
339	=> v Weight -- ^ weights @[w1, w2 .. wn]@
340	-> Weight -- ^ @l@
341	-> Maybe (Weight, VU.Vector Bool)
342	-- ^
343	-- * the objective value Σ_{i=1}^n wi xi, and
344	--
345	-- * the assignment @[x1, x2 .. xn]@, identifying 0 (resp. 1) with @False@ (resp. @True@).
346	minSubsetSum w l =	1✔
347	case maxSubsetSum w (wsum - l) of	1✔
348	Nothing -> Nothing	1✔
349	Just (obj, bs) -> Just (wsum - obj, VG.map not bs)	1✔
350	where
351	wsum = VG.sum w	1✔
352
353	{-
354	minimize Σ wi xi = Σ wi (1 - ¬xi) = Σ wi - (Σ wi ¬xi)
355	subject to Σ wi xi ≥ n
356
357	maximize Σ wi ¬xi
358	subject to Σ wi ¬xi ≤ (Σ wi) - n
359
360	Σ wi xi ≥ n
361	Σ wi (1 - ¬xi) ≥ n
362	(Σ wi) - (Σ wi ¬xi) ≥ n
363	(Σ wi ¬xi) ≤ (Σ wi) - n
364	-}
365
366	-- \| Solve Σ_{i=1}^n wi xi = c and xi ∈ {0,1}.
367	--
368	-- Note that this is different from usual definition of the subset sum problem,
369	-- as this definition allows all xi to be zero.
370	--
371	-- Note: 0 (resp. 1) is identified with False (resp. True) in the assignment.
372	subsetSum
373	:: VG.Vector v Weight
374	=> v Weight -- ^ weights @[w1, w2 .. wn]@
375	-> Weight -- ^ @l@
376	-> Maybe (VU.Vector Bool)
377	-- ^
378	-- the assignment @[x1, x2 .. xn]@, identifying 0 (resp. 1) with @False@ (resp. @True@).
379	subsetSum w c =	1✔
380	case normalizeWeightsToPositive (w,c) of	1✔
381	(w1, c1, trans1)
382	\| c1 < 0 -> Nothing	1✔
383	\| otherwise ->	×
384	case normalize2 (w1, c1) of	1✔
385	(w2, c2, trans2) -> do	1✔
386	(w3, c3, trans3) <- normalizeGCDEq (w2,c2)	1✔
387	let (obj, sol) = maxSubsetSum' w3 c3	1✔
388	guard $ obj == c3	1✔
389	return $ snd $ trans1 $ trans2 $ trans3 (obj, sol)	×

msakai / toysolver / 513

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous