{- |
Module : Solver.Internal
Description : Core functions to solve the distribution problem
Copyright : (c) Pablo Couto 2014
License : GPL-3
Maintainer : pablo@infty.in
Stability : experimental

This module contains the core functions used in solving the Generalized
Assignment Problem with the help of the GLPK toolkit.-}

{-# LANGUAGE CPP #-}

module Referees.Solver.Internal where

import Referees.Solver.Types.Internal
    ( Bounds(Bounds),
      ProfitFunction,
      ProfitMatrix,
      Col,
      Row,
      Index(Index, _idx),
      Copies,
      Capacity )

import Control.Exception ( SomeException, handle )
import Control.Monad ( guard )
import Control.Monad.LPMonad
    ( constrain,
      equalTo,
      execLPM,
      leqTo,
      setDirection,
      setObjective,
      setVarKind )
import Data.LinearProgram
    ( mipDefaults,
      ReturnCode,
      GLPOpts(msgLev),
      LP,
      VarKind(BinVar),
      Direction(Max),
      MsgLev(MsgAll, MsgOff),
      LinFunc,
      glpSolveVars,
      var,
      (*&),
      gsum )
import qualified Data.LinearProgram as LP ( Bounds(Bound) )
import Data.Map as Map ( Map, toList )
import Data.Matrix ( getElem, matrix, ncols, nrows )
import System.Exit ( exitFailure )

-- * Solving the Generalized Assignment Problem by approximation with GLPK

-- ** Formulation of the problem for GLPK

-- | 'lpGAP' is based on Martello & Toth 1990’s linear integer programming
-- formulation of the Generalized Assignment Problem (GAP). This function, in
-- addition, excludes all combinations whose profit is @0@. 'lpGAP' interfaces
-- with the <https://www.gnu.org/software/glpk/ GLPK> toolkit through @glpk-hs@.
--
-- In 'lpGAP' @profitM caps bounds@, @profitM@ is a profit matrix (which can be
-- computed by 'mkProfitMatrix', using a 'ProfitFunction'), @caps@ stands for a
-- list of items’ capacities (computable by 'toCap'), and @bounds@ encodes in a
-- 'Bounds' 'Copies' value (producible with 'mkBounds') both the lower and upper
-- bounds on the number of copies to distribute per each item.
--
-- * Cf. Martello, S. and Toth, P.: /Knapsack Problems: Algorithms and Computer/
-- /Implementations/, ch. 7, John Wiley & Sons, 1990.
--
lpGAP :: ProfitMatrix -> [Capacity] -> Bounds Copies -> LP String Double
lpGAP profitM cap (Bounds minCops maxCops) = execLPM $ do
  setDirection Max
  setObjective $ objFun profitM

  let (m, n) =
        (Index $ nrows profitM, Index $ ncols profitM) :: (Row, Col)
      setCap j =
        subFunCap profitM j `leqTo` fromIntegral (cap !! (_idx j - 1))
      setMult i =
        if minCops == maxCops
           then
             subFunMult profitM i `equalTo` fromIntegral maxCops
           else
             let bounds = LP.Bound (fromIntegral minCops) (fromIntegral maxCops)
             in subFunMult profitM i `constrain` bounds
      setVK (i, j) =
        setVarKind (x i j) BinVar
      setOff (i, j) =
        var (x i j) `equalTo` 0

  mapM_ setCap [1 .. n]
  mapM_ setMult [1 .. m]
  mapM_ setVK
    $ [(i, j) | i <- [1 .. m]
              , j <- [1 .. n]]
  mapM_ setOff
    $ [(i, j) | i <- [1 .. m]
              , j <- [1 .. n]
              , profitM `safeGetElem` (i, j) == 0]

-- | @'x' i j@, for @i = {1, …, m}@, @j = {1, …, n}@, where @m@ stands for the
-- number of items, and @n@ for the number of bins. These “variables” are to be
-- used only within the GAP formulation in 'lpGAP'.
--
x :: Row -> Col -> String
x i j = show (i, j)

-- | Objective function to maximize in 'lpGAP'. It stands for the overall profit
-- accrued by a given distribution of items among bins.
--
-- * Cf. Martello, S. and Toth, P.: op. cit.
--
objFun :: ProfitMatrix -> LinFunc String Double
objFun pM = gsum $ do
  let (m, n) = (Index $ nrows pM, Index $ ncols pM) :: (Row, Col)
  i <- [1 .. m]
  j <- [1 .. n]
  let p = pM `safeGetElem` (i, j)
  return $ p *& x i j

-- | Constraint function. It is used, in 'lpGAP', to specify the capacity of
-- each bin. Note that items are not divisible, and therefore their cost, as
-- detracted from a given bin’s capacity, is fixed at @1@.
--
-- * Cf. Martello, S. and Toth, P.: op. cit.
--
subFunCap :: ProfitMatrix -> Col -> LinFunc String Double
subFunCap pM j = gsum $ do
  let m = Index $ nrows pM :: Row
  i <- [1 .. m]
  return $ 1.0 *& x i j

-- | Constraint function. It is used, in 'lpGAP', to specify how many copies of
-- each item can be distributed.
--
-- * Cf. Martello, S. and Toth, P.: op. cit.
--
subFunMult :: ProfitMatrix -> Row -> LinFunc String Double
subFunMult pM i = gsum $ do
  let n = Index $ ncols pM :: Col
  j <- [1 .. n]
  return $ var $ x i j

-- ** Formulation accessories

-- | Given a 'ProfitFunction', 'mkProfitMatrix' computes a profit matrix between
-- @bins@ and @items@, optionally taking a @quality@ as being by default shared
-- between them. The values in a 'ProfitMatrix' serve to distribute items among
-- bins according to the capacities of the latter and the values of the former.
--
mkProfitMatrix :: ProfitFunction a b c
               -> [a] -- ^ Items
               -> [b] -- ^ Bins
               -> Maybe c -- ^ Optional quality
               -> ProfitMatrix
mkProfitMatrix f items bins quality =
  safeMatrix itemsNumber binsNumber profitFnMatrixWrapper
    where
      itemsNumber = Index $ length items :: Row
      binsNumber = Index $ length bins :: Col
      profitFnMatrixWrapper (i, b) = f i' b' quality
        where -- adapting Matrix to [] “indices”
          i' = items !! _idx (i - 1)
          b' = bins !! _idx (b - 1)

-- | Wrapper to use 'matrix' with extra type safety.
--
safeMatrix :: Row
           -> Col
           -> ((Row, Col) -> Double) -- ^ Specializes the @((Int, Int) -> a)@
                                     -- function used with 'matrix'
           -> ProfitMatrix
safeMatrix r c f = matrix (_idx r) (_idx c) f'
  where
    f' :: (Int, Int) -> Double
    f' (i, j) = f ((Index i, Index j) :: (Row, Col))

-- | Wrapper to use 'getElem' with extra type safety.
--
safeGetElem :: ProfitMatrix -> (Row, Col) -> Double
safeGetElem pM (r, c) = getElem (_idx r) (_idx c) pM

-- ** Runners and helpers

-- | Simple helper function to run the @glpk-hs@ interface to GLPK, as
-- constructed with 'lpGAP'. Takes the same arguments as 'lpGAP'.
--
run_lpGAP :: ProfitMatrix -> [Capacity] -> Bounds Copies
          -> IO (ReturnCode, Maybe (Double, Map String Double))
run_lpGAP profitM cap bnds =
  handle fuse $
#ifdef DEBUG
    glpSolveVars mipDefaults { msgLev = MsgAll }
#else
    glpSolveVars mipDefaults { msgLev = MsgOff }
#endif
      $ lpGAP profitM cap bnds
  where
    fuse e = flip const (e :: SomeException) $ do
      putStrLn $
        "Fatal error: unknown exception in solver. If this was"
        ++ "unexpected, please report the issue.\n"
        ++ "Error message: " ++ show e
      exitFailure

-- | 'fromGLPKtoList' turns the (unIOd) output of 'run_lpGAP' into a more usable
-- format.
--
fromGLPKtoList :: (ReturnCode, Maybe (Double, Map String Double))
               -> Maybe [(Col, Row)]
fromGLPKtoList (_, output) =
  case output of
    Just (_, output') ->
      Just $ do
        (s, v) <- Map.toList output'
        guard $ v == 1.0
        return $ fixIndices . read $ s :: [(Col, Row)]
    Nothing ->
      Just []
  where -- turn them around
    fixIndices (i, j) = (j - 1, i - 1)