Generator.hs

{-# LANGUAGE TypeFamilies #-}
module Test.Numeric.GLPK.Generator where

import qualified Numeric.GLPK as LP
import Numeric.GLPK ((<=.), (>=.))

import qualified Test.QuickCheck as QC
import System.Random (Random)

import qualified Data.Array.Comfort.Boxed as BoxedArray
import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape
import qualified Data.NonEmpty as NonEmpty
import qualified Data.List.HT as ListHT
import qualified Data.Ix as Ix
import Data.Array.Comfort.Storable (Array, (!))
import Data.Traversable (sequenceA, for)
import Data.Tuple.HT (mapSnd)
import Data.Maybe (fromMaybe)
import Data.Int (Int64)

import Control.Applicative (liftA2)

import Text.Printf (PrintfArg, printf)

import Foreign.Storable (Storable)


{- |
Generate constraints in the form of a polyhedron
which contains warrantedly the zero vector.
That is, there is an admissible solution.
In order to assert that the polyhedron is closed,
we bound all variables by a hypercube.
-}
genProblem ::
   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
   Array sh a -> QC.Gen (LP.Bounds ix, LP.Constraints ix)
genProblem origin =
   liftA2 (,)
      (for (Array.toAssociations origin) $ \(ix,x) ->
         LP.Inequality ix <$>
         liftA2 LP.Between
            (doubleFromElement . (x+) <$> QC.choose (-100,-50))
            (doubleFromElement . (x+) <$> QC.choose (50,100)))
      (do
         numConstraints <- QC.choose (1,20)
         QC.vectorOf numConstraints $ do
            ixs <- QC.sublistOf $ Shape.indices $ Array.shape origin
            terms <- for ixs $ \ix -> do
               coeff <- QC.choose (-10,10)
               return (coeff, ix)
            let offset = scalarProduct terms origin
            let deviation = 25
            LP.Inequality
               (map (uncurry (LP.Term . doubleFromElement)) terms)
               <$>
               QC.oneof (
                  (do bound <- QC.choose (offset-deviation, offset+deviation)
                      return $
                         if bound > offset
                            then LP.LessEqual    $ doubleFromElement bound
                            else LP.GreaterEqual $ doubleFromElement bound) :
                  (liftA2 LP.Between
                     (doubleFromElement <$>
                        QC.choose (offset-deviation, offset))
                     (doubleFromElement <$>
                        QC.choose (offset, offset+deviation))) :
                  []))

scalarProduct ::
   (Shape.Indexed sh, Shape.Index sh ~ ix, Storable a, Num a) =>
   [(a,ix)] -> Array sh a -> a
scalarProduct terms origin =
   sum $ map (\(coeff, ix) -> coeff * origin!ix) terms

genVarShape :: QC.Gen (Shape.Range Char)
genVarShape = Shape.Range 'a' <$> QC.choose ('a','j')

genOrigin :: QC.Gen (Array (Shape.Range Char) Int64)
genOrigin = genVector =<< genVarShape

_genOrigin :: QC.Gen (Array (Shape.Range Char) Double)
_genOrigin = genVector =<< genVarShape


shrinkVarShape :: Shape.Range Char -> [Shape.Range Char]
shrinkVarShape (Shape.Range from to) =
   if from<to then [Shape.Range from (pred to)] else []

shrinkOrigin ::
   (Storable a) => Array (Shape.Range Char) a -> [Array (Shape.Range Char) a]
shrinkOrigin vec =
   case Array.shape vec of
      Shape.Range from to ->
         if from<to
            then [Array.sample (Shape.Range from (pred to)) (vec!)]
            else []


class (Storable a, Random a, Num a, Ord a) => Element a where
   doubleFromElement :: a -> Double

instance Element Double where
   doubleFromElement = id

instance Element Int64 where
   doubleFromElement = fromIntegral

genObjective ::
   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
   Array sh a -> QC.Gen (LP.Direction, LP.Objective sh)
genObjective origin =
   liftA2 (,) QC.arbitraryBoundedEnum
      (fmap (Array.map doubleFromElement . flip asTypeOf origin) $
       genVector $ Array.shape origin)

genVector :: (Shape.Indexed sh, Element a) => sh -> QC.Gen (Array sh a)
genVector shape =
   fmap Array.fromBoxed $ sequenceA $
   BoxedArray.fromAssociations (QC.choose (-10,10)) shape []

shrinkProblem ::
   (LP.Bounds ix, LP.Constraints ix) ->
   [(LP.Bounds ix, LP.Constraints ix)]
shrinkProblem (bounds, constraints) =
   map (\shrinked -> (bounds, shrinked)) $
   filter (not . null) $ QC.shrinkList (const []) constraints

genObjectives ::
   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>
   Array sh a -> QC.Gen (NonEmpty.T [] (LP.Direction, [LP.Term ix]))
genObjectives origin = do
   let shape = Array.shape origin
   let stageRange :: (Int,Int)
       stageRange = (0,3)
   stages <- for (Shape.indices shape) $ \ix -> (,) ix <$> QC.choose stageRange
   let varSets =
         fromMaybe (error "there should be at least one stage") $
         NonEmpty.fetch $
         filter (not . null) $
         map (\k -> map fst $ filter ((k==) . snd) stages) $
         Ix.range stageRange
   let asTypeOfElement :: a -> f a -> a
       asTypeOfElement = const
   for varSets $ \varSet ->
      liftA2 (,)
         QC.arbitraryBoundedEnum
         (for varSet $ \ix ->
            flip LP.Term ix . doubleFromElement
               <$> QC.choose (-10, 10 `asTypeOfElement` origin))

shrinkObjectives ::
   NonEmpty.T [] (LP.Direction, [LP.Term ix]) ->
   [NonEmpty.T [] (LP.Direction, [LP.Term ix])]
shrinkObjectives (NonEmpty.Cons obj objs) =
   map (NonEmpty.Cons obj) $
   QC.shrinkList
      (\(dir,terms) ->
         map ((,) dir) $ filter (not . null) $
         QC.shrinkList (const []) terms)
      objs

successiveObjectives ::
   (Shape.Indexed sh, Shape.Index sh ~ ix) =>
   Array sh a -> Double ->
   NonEmpty.T [] (LP.Direction, [LP.Term ix]) ->
   ((LP.Direction, LP.Objective sh),
    [((LP.SolutionType, (Double, Array sh Double)) -> LP.Constraints ix,
      (LP.Direction, LP.Objective sh))])
successiveObjectives origin tol xs =
   let shape = Array.shape origin in
   (mapSnd (LP.objectiveFromTerms shape) $ NonEmpty.head xs,
    NonEmpty.mapAdjacent
      (\(dir0,obj0) y1 ->
         (\(_sol,(opt,_vec)) ->
            case dir0 of
               LP.Minimize -> [obj0 <=. opt + tol]
               LP.Maximize -> [obj0 >=. opt - tol],
          mapSnd (LP.objectiveFromTerms shape) y1))
      xs)


approxReal :: (Ord a, Num a) => a -> a -> a -> Bool
approxReal tol x y = abs (x-y) <= tol

approx :: (PrintfArg a, Ord a, Num a) => String -> a -> a -> a -> QC.Property
approx name tol x y =
   QC.counterexample (printf "%s: %f - %f" name x y) (approxReal tol x y)

approxSuccession ::
   (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>
   a ->
   Either LP.NoSolutionType
      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->
   Either LP.NoSolutionType
      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->
   QC.Property
approxSuccession tol x y =
   QC.counterexample (show x) $
   QC.counterexample (show y) $
   case (x,y) of
      (Left sx, Left sy) -> sx==sy
      (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->
         let equalSol (solX, (optX, _)) (solY, (optY, _)) =
               solX == solY && approxReal tol optX optY
         in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys
      _ -> False