content: {-# OPTIONS_GHC -Wall #-}
  -----------------------------------------------------------------------------
  -- |
  -- Module      :  ToySolver.Combinatorial.HittingSet.FredmanKhachiyan1996
  -- Copyright   :  (c) Masahiro Sakai 2015
  -- License     :  BSD-style
  --
  -- Maintainer  :  masahiro.sakai@gmail.com
  -- Stability   :  provisional
  -- Portability :  portable
  --
  -- This modules provides functions to check whether two monotone boolean
  -- functions /f/ and /g/ given in DNFs are mutually dual (/i.e./ f(x1,…,xn) = ¬g(¬x1,…,¬xn)).
  --
  -- References:
  --
  -- * [FredmanKhachiyan1996] Michael L. Fredman and Leonid Khachiyan,
  --   On the Complexicy of Dualization of Monotone Disjunctifve Normal Forms,
  --   Journal of Algorithms, vol. 21, pp. 618-628, 1996.
  --   <http://www.sciencedirect.com/science/article/pii/S0196677496900620>
  --   <http://www.cs.tau.ac.il/~fiat/dmsem03/On%20the%20complexity%20of%20Dualization%20of%20monotone%20Disjunctive%20Normal%20Forms%20-%201996.pdf>
  --
  -----------------------------------------------------------------------------
  module ToySolver.Combinatorial.HittingSet.FredmanKhachiyan1996
    (
    -- * Main functions
      areDualDNFs
    , checkDuality
    , checkDualityA
    , checkDualityB
  
    -- * Redundancy
    -- $Redundancy
    , isRedundant
    , deleteRedundancy
  
    -- * Utilities for testing
    , isCounterExampleOf
    , occurFreq
  
    -- * Internal functions exported only for testing purpose
    , condition_1_1
    , condition_1_2
    , condition_1_3
    , condition_2_1
    , condition_1_1_solve
    , condition_1_2_solve
    , condition_1_3_solve
    , condition_2_1_solve
    ) where
  
  import Prelude hiding (all, any)
  import Control.Arrow ((***))
  import Control.Exception (assert)
  import Control.Monad
  import Data.Foldable (all, any)
  import Data.IntSet (IntSet)
  import qualified Data.IntSet as IntSet
  import Data.List hiding (all, any, intersect)
  import Data.Maybe
  import Data.Ratio
  import Data.Set (Set)
  import qualified Data.Set as Set
  
  -- | xhi n ** xhi n == n
  xhi :: Double -> Double
  xhi n = iterate f m !! 30
    where
      m = logn / logBase 2 logn
        where
          logn = logBase 2 n
      f x = m * ((logx + logBase 2 logx) / logx)
        where
          logx = logBase 2 x
  
  disjoint :: IntSet -> IntSet -> Bool
  disjoint xs ys = IntSet.null $ xs `IntSet.intersection` ys
  
  intersect :: IntSet -> IntSet -> Bool
  intersect xs ys = not $ disjoint xs ys
  
  isHittingSetOf :: IntSet -> Set IntSet -> Bool
  isHittingSetOf xs yss = all (xs `intersect`) yss
  
  isCounterExampleOf :: IntSet -> (Set IntSet, Set IntSet) -> Bool
  isCounterExampleOf xs (f,g) = lhs == rhs
    where
      lhs = or [is `IntSet.isSubsetOf` xs | is <- Set.toList f]
      rhs = or [xs `disjoint` js | js <- Set.toList g]
  
  _volume :: Set IntSet -> Set IntSet -> Int
  _volume f g = Set.size f * Set.size g
  
  condition_1_1 :: Set IntSet -> Set IntSet -> Bool
  condition_1_1 f g = all (\is -> all (\js -> is `intersect` js) g) f
  
  condition_1_1_solve :: Set IntSet -> Set IntSet -> Maybe IntSet
  condition_1_1_solve f g = listToMaybe $ do
    is <- Set.toList f
    js <- Set.toList g
    guard $ is `disjoint` js
    return is
  
  condition_1_2 :: Set IntSet -> Set IntSet -> Bool
  condition_1_2 f g = h f == h g
    where
      h = IntSet.unions . Set.toList
  
  -- | Find @xs@ such that @xs `isCounterExampleOf` (f,g)@ unless @'condition_1_2' f g@.
  condition_1_2_solve :: Set IntSet -> Set IntSet -> Maybe IntSet
  condition_1_2_solve f g
    | Just (v1,_) <- IntSet.minView d1 = Just $ head [IntSet.delete v1 is | is <- Set.toList f, v1 `IntSet.member` is]
    | Just (v2,_) <- IntSet.minView d2 = Just $ head [(vs `IntSet.difference` (IntSet.delete v2 js)) | js <- Set.toList g, v2 `IntSet.member` js]
    | otherwise = Nothing
    where
      f_vs = IntSet.unions $ Set.toList f
      g_vs = IntSet.unions $ Set.toList g
      vs = f_vs `IntSet.union` g_vs
      d1 = f_vs `IntSet.difference` g_vs
      d2 = g_vs `IntSet.difference` f_vs
  
  condition_1_3 :: Set IntSet -> Set IntSet -> Bool
  condition_1_3 f g = maxSize f <= Set.size g && maxSize g <= Set.size f
    where
      maxSize xs = foldl' max 0 [IntSet.size i | i <- Set.toList xs]
  
  condition_1_3_solve :: Set IntSet -> Set IntSet -> Maybe IntSet
  condition_1_3_solve f g = listToMaybe $
    [head [is' | i <- IntSet.toList is, let is' = IntSet.delete i is, is' `isHittingSetOf` g] | is <- Set.toList f, IntSet.size is > g_size] ++
    [xs `IntSet.difference` head [js' | j <- IntSet.toList js, let js' = IntSet.delete j js, js' `isHittingSetOf` f] | js <- Set.toList g, IntSet.size js > f_size]
    where
      xs = IntSet.unions $ Set.toList $ Set.union f g
      f_size = Set.size f
      g_size = Set.size g
  
  e :: Set IntSet -> Set IntSet -> Rational
  e f g = sum [1 % (2 ^ IntSet.size i) | i <- Set.toList f] +
          sum [1 % (2 ^ IntSet.size j) | j <- Set.toList g]
  
  condition_2_1 :: Set IntSet -> Set IntSet -> Bool
  condition_2_1 f g = e f g >= 1
  
  condition_2_1_solve :: Set IntSet -> Set IntSet -> Maybe IntSet
  condition_2_1_solve f g =
    if condition_2_1 f g
    then Nothing
    else Just $ loop (IntSet.toList vs) f g IntSet.empty
    where
      vs = IntSet.unions $ Set.toList $ Set.union f g
  
      loop :: [Int] -> Set IntSet -> Set IntSet -> IntSet -> IntSet
      loop [] _ _ ret = ret
      loop (v:vs) f g ret =
        if e f1 g1 <= e f2 g2
        then loop vs f1 g1 (IntSet.insert v ret)
        else loop vs f2 g2 ret
        where
          -- f = x f1 ∨ f2
          -- g = x g2 ∨ g1
          f1 = Set.map (IntSet.delete v) f
          g1 = Set.filter (v `IntSet.notMember`) g
          f2 = Set.filter (v `IntSet.notMember`) f
          g2 = Set.map (IntSet.delete v) g
  
  -- | @'isRedundant' F@ tests whether /F/ contains redundant implicants.
  isRedundant :: Set IntSet -> Bool
  isRedundant = loop . sortOn IntSet.size . Set.toList
    where
      loop :: [IntSet] -> Bool
      loop [] = False
      loop (xs:yss) = or [xs `IntSet.isSubsetOf` ys | ys <- yss] || loop yss
  
  isIrredundant :: Set IntSet -> Bool
  isIrredundant = not . isRedundant
  
  -- | Removes redundant implicants from a given DNF.
  deleteRedundancy :: Set IntSet -> Set IntSet
  deleteRedundancy = foldl' f Set.empty . sortOn IntSet.size . Set.toList
    where
      f :: Set IntSet -> IntSet -> Set IntSet
      f xss ys =
        if any (`IntSet.isSubsetOf` ys) xss
        then xss
        else Set.insert ys xss
  
  -- | @occurFreq  x F@ computes /|{I∈F | x∈I}| \/ |F|/
  occurFreq
    :: Fractional a
    => Int -- ^ a variable /x/
    -> Set IntSet -- ^ a DNF /F/
    -> a
  occurFreq x f = fromIntegral (sum [1 | is <- Set.toList f, x `IntSet.member` is] :: Int) / fromIntegral (Set.size f)
  
  -- | @areDualDNFs f g@ checks whether two monotone boolean functions /f/
  -- and /g/ are mutually dual (/i.e./ f(x1,…,xn) = ¬g(¬x1,…,xn)).
  --
  -- Note that this function does not care about redundancy of DNFs.
  --
  -- Complexity: /O(n^o(log n))/ where @n = 'IntSet.size' f + 'IntSet.size' g@.
  areDualDNFs
    :: Set IntSet -- ^ a monotone boolean function /f/ given in DNF
    -> Set IntSet -- ^ a monotone boolean function /g/ given in DNF
    -> Bool
  areDualDNFs f g = isNothing $ checkDualityB f g
  
  -- | Synonym for 'checkDualityB'.
  checkDuality :: Set IntSet -> Set IntSet -> Maybe IntSet
  checkDuality = checkDualityB
  
  -- | @checkDualityA f g@ checks whether two monotone boolean functions /f/
  -- and /g/ are mutually dual (/i.e./ f(x1,…,xn) = ¬g(¬x1,…,¬xn)) using
  -- “Algorithm A” of [FredmanKhachiyan1996].
  --
  -- If they are indeed mutually dual it returns @Nothing@, otherwise
  -- it returns @Just cs@ such that {xi ↦ (if xi∈cs then True else False) | i∈{1…n}}
  -- is a solution of f(x1,…,xn) = g(¬x1,…,¬xn)).
  --
  -- Note that this function does not care about redundancy of DNFs.
  --
  -- Complexity: /O(n^O(log^2 n))/ where @n = 'IntSet.size' f + 'IntSet.size' g@.
  checkDualityA
    :: Set IntSet -- ^ a monotone boolean function /f/ given in DNF
    -> Set IntSet -- ^ a monotone boolean function /g/ given in DNF
    -> Maybe IntSet
  checkDualityA f g
    | Just xs <- condition_1_1_solve f g = Just xs
    | otherwise = checkDualityA' (deleteRedundancy f) (deleteRedundancy g)
  
  checkDualityA' :: Set IntSet -> Set IntSet -> Maybe IntSet
  checkDualityA' f g
    | assert (isIrredundant f && isIrredundant g) False = undefined
    | Just s <- condition_1_2_solve f g = Just s
    | Just s <- condition_1_3_solve f g = Just s
    | Just s <- condition_2_1_solve f g = Just s
    | v <= 1 = solveSmall f g
    | otherwise = msum
        [ -- If x=0 then f(xs)=g(¬xs) ⇔ f1(xs) = g0(¬xs) ∨ g1(¬xs)
          checkDualityA' f1 (deleteRedundancy (g0 `Set.union` g1))
        , -- If x=1 then f(xs)=g(¬xs) ⇔ f0(xs) ∨ f1(xs) = g1(¬xs)
          liftM (IntSet.insert x) $ checkDualityA' (deleteRedundancy (f0 `Set.union` f1)) g1
        ]
    where
      size_f = Set.size f
      size_g = Set.size g
      n = size_f + size_g
      v = size_f * size_g
      xs = IntSet.unions $ Set.toList $ f `Set.union` g
      epsilon :: Double
      epsilon = 1 / logBase 2 (fromIntegral n)
      x = head [x | x <- IntSet.toList xs, occurFreq x f >= epsilon || occurFreq x g >= epsilon]
      -- f = x f0 ∨ f1
      (f0, f1) = Set.map (IntSet.delete x) *** id $ Set.partition (x `IntSet.member`) f
      -- g = x g0 ∨ g1
      (g0, g1) = Set.map (IntSet.delete x) *** id $ Set.partition (x `IntSet.member`) g
  
  solveSmall :: Set IntSet -> Set IntSet -> Maybe IntSet
  solveSmall f g
    | assert (isIrredundant f && isIrredundant g) False = undefined
    | Set.null f && Set.null g = Just IntSet.empty
    | Set.null f = assert (not (Set.null g)) $
        if IntSet.empty `Set.member` g
        then Nothing
        else Just (IntSet.unions (Set.toList g))
    | Set.null g = assert (not (Set.null f)) $
        if IntSet.empty `Set.member` f
        then Nothing
        else Just IntSet.empty
    | size_f == 1 && size_g == 1 =
        let
          is = Set.findMin f
          js = Set.findMin g
          is_size = IntSet.size is
          js_size = IntSet.size js
        in if is_size == 1 then
             if js_size == 1 then
               Nothing
             else
               Just (js `IntSet.difference` is)
           else
             Just (is `IntSet.difference` js)
    | otherwise = error "should not happen"
    where
      size_f = Set.size f
      size_g = Set.size g
  
  -- | @checkDualityB f g@ checks whether two monotone boolean functions /f/
  -- and /g/ are mutually dual (i.e. f(x1,…,xn) = ¬g(¬x1,…,xn)) using
  -- “Algorithm B” of [FredmanKhachiyan1996].
  --
  -- If they are indeed mutually dual it returns @Nothing@, otherwise
  -- it returns @Just cs@ such that {xi ↦ (if xi∈cs then True else False) | i∈{1…n}}
  -- is a solution of f(x1,…,xn) = g(¬x1,…,xn)).
  --
  -- Note that this function does not care about redundancy of DNFs.
  --
  -- Complexity: /O(n^o(log n))/ where @n = 'Set.size' f + 'Set.size' g@.
  checkDualityB
    :: Set IntSet -- ^ a monotone boolean function /f/ given in DNF
    -> Set IntSet -- ^ a monotone boolean function /g/ given in DNF
    -> Maybe IntSet
  checkDualityB f g
    | Just xs <- condition_1_1_solve f g = Just xs
    | otherwise = checkDualityB' (deleteRedundancy f) (deleteRedundancy g)
  
  checkDualityB' :: Set IntSet -> Set IntSet -> Maybe IntSet
  checkDualityB' f g
    | assert (isIrredundant f && isIrredundant g) False = undefined
    | Just s <- condition_1_2_solve f g = Just s
    | Just s <- condition_1_3_solve f g = Just s
    | Just s <- condition_2_1_solve f g = Just s
    | v <= 1 = solveSmall f g
  --  | min size_f size_g <= 2 = undefined
    | otherwise =
        let x = head [x | x <- IntSet.toList xs, occurFreq x f > (0::Rational), occurFreq x g > (0::Rational)]
            -- f = x f0 ∨ f1
            (f0, f1) = Set.map (IntSet.delete x) *** id $ Set.partition (x `IntSet.member`) f
            -- g = x g0 ∨ g1
            (g0, g1) = Set.map (IntSet.delete x) *** id $ Set.partition (x `IntSet.member`) g
            epsilon_x_f :: Double
            epsilon_x_f = occurFreq x f
            epsilon_x_g :: Double
            epsilon_x_g = occurFreq x g
        in if epsilon_x_f <= epsilon_v then
             msum
             [ {- f(xs ∪ {x↦0})=g(¬xs ∪ {x↦1}) ⇔ f1(xs) = g0(¬xs) ∨ g1(¬xs). -}
               checkDualityB' f1 (deleteRedundancy (g0 `Set.union` g1))
             , {- f(¬xs ∪ {x↦1}) = g(xs ∪ {x↦0})
                  ⇔ f0(¬xs) ∨ f1(¬xs) = g1(xs)
                  ⇔ f0(¬xs) ∨ (¬g0(xs)∧¬g1(xs)) = g1(xs) {- duality of f1 and g0∨g1 -}
                  ⇔ f0(¬xs) = g1(xs) and g0(xs) = 1
                     {- g0(xs)=0 is impossible, because
                        it implies f0(¬xs)∨¬g1(xs)=g1(xs) and then
                        f0(¬xs)=g1(xs)=1 which contradicts condition (1.1) -}
                -}
               msum $ do
                 js <- Set.toList g0
                 let f' = Set.filter (js `disjoint`) f0
                     g' = Set.map (`IntSet.difference` js) g1
                 return $ do
                   ys <- checkDualityB' (deleteRedundancy g') (deleteRedundancy f')
                   let ys2 = IntSet.insert x $ xs `IntSet.difference` (js `IntSet.union` ys)
                   assert (ys2 `isCounterExampleOf` (f,g)) $ return ()
                   return ys2
             ]
           else if epsilon_x_g <= epsilon_v then
             msum
             [ {- f(xs ∪ {x↦1}) = g(¬xs ∪ {x↦0})) ⇔ f0(xs) ∨ f1(xs) = g1(¬xs). -}
               liftM (IntSet.insert x) $ checkDualityB' (deleteRedundancy (f0 `Set.union` f1)) g1
             , {- f(xs ∪ {x↦0}) = g(¬xs ∪ {x↦1})
                  ⇔ f1(xs) = g0(¬xs) ∨ g1(¬xs)
                  ⇔ f1(xs) = g0(¬xs) ∨ (¬f0(xs)∧¬f1(xs))  {- duality of g1 and f0∨f1 -}
                  ⇔ f1(xs) = g0(¬xs) and f0(xs) = 1
                     {- f0(xs)=0 is impossible, because
                        it implies f1(xs)=g0(¬xs)∨¬f1(xs) and then
                        f1(xs)=g0(¬xs)=1 which contradicts condition (1.1) -}
                -}
               msum $ do
                 is <- Set.toList f0
                 let f' = Set.map (`IntSet.difference` is) f1
                     g' = Set.filter (is `disjoint`) g0
                 return $ do
                   ys <- checkDualityB' (deleteRedundancy f') (deleteRedundancy g')
                   let ys2 = is `IntSet.union` ys
                   assert (ys2 `isCounterExampleOf` (f,g)) $ return ()
                   return ys2
             ]
           else -- epsilon_v < min epsilon_x_f epsilon_x_g
             msum
             [ -- If x=0 then f(xs)=g(¬xs) ⇔ f1(xs) = g0(¬xs) ∨ g1(¬xs)
               checkDualityB' f1 (deleteRedundancy (g0 `Set.union` g1))
             , -- If x=1 then f(xs)=g(¬xs) ⇔ f0(xs) ∨ f1(xs) = g1(¬xs)
               liftM (IntSet.insert x) $ checkDualityB' (deleteRedundancy (f0 `Set.union` f1)) g1
             ]
    where
      size_f = Set.size f
      size_g = Set.size g
      v = size_f * size_g
      epsilon_v = 1 / xhi (fromIntegral v)
      xs = IntSet.unions $ Set.toList $ f `Set.union` g
  
  -- $Redundancy
  -- An implicant /I∈F/ of a DNF /F/ is redundant if /F/ contains proper subset of /I/.
  -- A DNF /F/ is called redundant if it contains redundant implicants.
  -- The main functions of this modules does not care about redundancy of DNFs.
  