module Geometry.VertexEnum.VertexEnum
  ( vertexenum
  , checkConstraints
  , interiorPoint )
  where
import           Control.Monad                   ( unless, when, (<$!>) )
import           Data.Maybe                      ( isJust, fromJust )
import           Foreign.C.Types                 ( CDouble, CUInt )
import           Foreign.Marshal.Alloc           ( free, mallocBytes )
import           Foreign.Marshal.Array           ( peekArray, pokeArray )
import           Foreign.Storable                ( peek, sizeOf )
import           Geometry.VertexEnum.CVertexEnum ( c_intersections )
import           Geometry.VertexEnum.Constraint  ( Constraint, toRationalConstraint )
import           Geometry.VertexEnum.Internal    ( iPoint, normalizeConstraints, findSigns )

hsintersections :: [[Double]]     -- halfspaces

                 -> [Double]      -- interior point

                 -> Bool          -- print to stdout

                 -> IO [[Double]]
hsintersections :: [[Double]] -> [Double] -> Bool -> IO [[Double]]
hsintersections [[Double]]
halfspaces [Double]
ipoint Bool
stdout = do
  let n :: Int
n     = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
halfspaces
      dim :: Int
dim   = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ipoint
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
dimInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
halfspaces) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
    [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error [Char]
"the length of the point does not match the number of variables."
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
dim Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
    [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error [Char]
"dimension must be at least 2."
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
dim) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
    [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error [Char]
"insufficient number of halfspaces."
  Ptr CDouble
hsPtr <- Int -> IO (Ptr CDouble)
forall a. Int -> IO (Ptr a)
mallocBytes (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
dimInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* CDouble -> Int
forall a. Storable a => a -> Int
sizeOf (CDouble
forall a. HasCallStack => a
undefined :: CDouble))
  Ptr CDouble -> [CDouble] -> IO ()
forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CDouble
hsPtr (([Double] -> [CDouble]) -> [[Double]] -> [CDouble]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((Double -> CDouble) -> [Double] -> [CDouble]
forall a b. (a -> b) -> [a] -> [b]
map Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac) [[Double]]
halfspaces)
  Ptr CDouble
ipointPtr <- Int -> IO (Ptr CDouble)
forall a. Int -> IO (Ptr a)
mallocBytes (Int
dim Int -> Int -> Int
forall a. Num a => a -> a -> a
* CDouble -> Int
forall a. Storable a => a -> Int
sizeOf (CDouble
forall a. HasCallStack => a
undefined :: CDouble))
  Ptr CDouble -> [CDouble] -> IO ()
forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CDouble
ipointPtr ((Double -> CDouble) -> [Double] -> [CDouble]
forall a b. (a -> b) -> [a] -> [b]
map Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac [Double]
ipoint)
  Ptr CUInt
exitcodePtr <- Int -> IO (Ptr CUInt)
forall a. Int -> IO (Ptr a)
mallocBytes (CUInt -> Int
forall a. Storable a => a -> Int
sizeOf (CUInt
forall a. HasCallStack => a
undefined :: CUInt))
  Ptr CUInt
nintersectionsPtr <- Int -> IO (Ptr CUInt)
forall a. Int -> IO (Ptr a)
mallocBytes (CUInt -> Int
forall a. Storable a => a -> Int
sizeOf (CUInt
forall a. HasCallStack => a
undefined :: CUInt))
  Ptr (Ptr CDouble)
resultPtr <- Ptr CDouble
-> Ptr CDouble
-> CUInt
-> CUInt
-> Ptr CUInt
-> Ptr CUInt
-> CUInt
-> IO (Ptr (Ptr CDouble))
c_intersections Ptr CDouble
hsPtr Ptr CDouble
ipointPtr
               (Int -> CUInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dim) (Int -> CUInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
               Ptr CUInt
nintersectionsPtr Ptr CUInt
exitcodePtr (Int -> CUInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CUInt) -> Int -> CUInt
forall a b. (a -> b) -> a -> b
$ Bool -> Int
forall a. Enum a => a -> Int
fromEnum Bool
stdout)
  CUInt
exitcode <- Ptr CUInt -> IO CUInt
forall a. Storable a => Ptr a -> IO a
peek Ptr CUInt
exitcodePtr
  Ptr CUInt -> IO ()
forall a. Ptr a -> IO ()
free Ptr CUInt
exitcodePtr
  Ptr CDouble -> IO ()
forall a. Ptr a -> IO ()
free Ptr CDouble
hsPtr
  if CUInt
exitcode CUInt -> CUInt -> Bool
forall a. Eq a => a -> a -> Bool
/= CUInt
0
    then do
      Ptr (Ptr CDouble) -> IO ()
forall a. Ptr a -> IO ()
free Ptr (Ptr CDouble)
resultPtr
      Ptr CUInt -> IO ()
forall a. Ptr a -> IO ()
free Ptr CUInt
nintersectionsPtr
      [Char] -> IO [[Double]]
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO [[Double]]) -> [Char] -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ [Char]
"qhull returned an error (code " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ CUInt -> [Char]
forall a. Show a => a -> [Char]
show CUInt
exitcode [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
")."
    else do
      Int
nintersections <- (CUInt -> Int) -> IO CUInt -> IO Int
forall (m :: * -> *) a b. Monad m => (a -> b) -> m a -> m b
(<$!>) CUInt -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Ptr CUInt -> IO CUInt
forall a. Storable a => Ptr a -> IO a
peek Ptr CUInt
nintersectionsPtr)
      [[Double]]
result <- ([[CDouble]] -> [[Double]]) -> IO [[CDouble]] -> IO [[Double]]
forall (m :: * -> *) a b. Monad m => (a -> b) -> m a -> m b
(<$!>) (([CDouble] -> [Double]) -> [[CDouble]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((CDouble -> Double) -> [CDouble] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac))
                       (([Ptr CDouble] -> IO [[CDouble]])
-> IO [Ptr CDouble] -> IO [[CDouble]]
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
(=<<) ((Ptr CDouble -> IO [CDouble]) -> [Ptr CDouble] -> IO [[CDouble]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (Int -> Ptr CDouble -> IO [CDouble]
forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
dim))
                              (Int -> Ptr (Ptr CDouble) -> IO [Ptr CDouble]
forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
nintersections Ptr (Ptr CDouble)
resultPtr))
      Ptr (Ptr CDouble) -> IO ()
forall a. Ptr a -> IO ()
free Ptr (Ptr CDouble)
resultPtr
      Ptr CUInt -> IO ()
forall a. Ptr a -> IO ()
free Ptr CUInt
nintersectionsPtr
      [[Double]] -> IO [[Double]]
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return [[Double]]
result

-- | Vertex enumeration

vertexenum :: 
    Real a 
  => [Constraint a] -- ^ linear inequalities

  -> Maybe [Double] -- ^ point satisfying the inequalities, @Nothing@ for automatic point

  -> IO [[Double]]  -- ^ vertices of the polytope defined by the inequalities

vertexenum :: forall a.
Real a =>
[Constraint a] -> Maybe [Double] -> IO [[Double]]
vertexenum [Constraint a]
constraints Maybe [Double]
point = do
  let halfspacesMatrix :: [[Double]]
halfspacesMatrix = 
        ([a] -> [Double]) -> [[a]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> Double) -> [a] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map a -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac) ([Constraint a] -> [[a]]
forall a. Real a => [Constraint a] -> [[a]]
normalizeConstraints [Constraint a]
constraints)
  if Maybe [Double] -> Bool
forall a. Maybe a -> Bool
isJust Maybe [Double]
point
    then do
      let check :: [(Double, Bool)]
check = [Constraint a] -> [Double] -> [(Double, Bool)]
forall a. Real a => [Constraint a] -> [Double] -> [(Double, Bool)]
checkConstraints [Constraint a]
constraints (Maybe [Double] -> [Double]
forall a. HasCallStack => Maybe a -> a
fromJust Maybe [Double]
point)
      Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ ((Double, Bool) -> Bool) -> [(Double, Bool)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double, Bool) -> Bool
forall a b. (a, b) -> b
snd [(Double, Bool)]
check) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
        [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error [Char]
"vertexenum: the provided point does not fulfill the inequalities."
      [[Double]] -> [Double] -> Bool -> IO [[Double]]
hsintersections [[Double]]
halfspacesMatrix (Maybe [Double] -> [Double]
forall a. HasCallStack => Maybe a -> a
fromJust Maybe [Double]
point) Bool
False
    else do
      [Double]
ipoint <- [Constraint a] -> IO [Double]
forall a. Real a => [Constraint a] -> IO [Double]
interiorPoint [Constraint a]
constraints 
      [[Double]] -> [Double] -> Bool -> IO [[Double]]
hsintersections [[Double]]
halfspacesMatrix [Double]
ipoint Bool
False

-- | Checks whether a point fulfills some inequalities; returns the 

-- difference between the upper member and the lower member for each

-- inequality, which is positive in case if the inequality is fulfilled.

checkConstraints :: 
  Real a 
  => [Constraint a]   -- ^ linear inequalities

  -> [Double]         -- ^ point to be tested

  -> [(Double, Bool)] -- ^ difference and status for each constraint

checkConstraints :: forall a. Real a => [Constraint a] -> [Double] -> [(Double, Bool)]
checkConstraints [Constraint a]
constraints [Double]
point = 
  if Int
nvars Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
point Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
    then 
      [Double] -> [Bool] -> [(Double, Bool)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
differences ((Double -> Bool) -> [Double] -> [Bool]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) [Double]
differences)
    else 
      [Char] -> [(Double, Bool)]
forall a. HasCallStack => [Char] -> a
error [Char]
"checkConstraints: the length of the point does not match the number of variables."
  where
    halfspacesMatrix :: [[Double]]
halfspacesMatrix = 
      ([a] -> [Double]) -> [[a]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> Double) -> [a] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map a -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac) ([Constraint a] -> [[a]]
forall a. Real a => [Constraint a] -> [[a]]
normalizeConstraints [Constraint a]
constraints)
    nvars :: Int
nvars = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Double]]
halfspacesMatrix [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
    checkRow :: [a] -> [a] -> a
checkRow [a]
pt [a]
row = - [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
row ([a]
pt [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
1]))
    differences :: [Double]
differences = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
forall {a}. Num a => [a] -> [a] -> a
checkRow [Double]
point) [[Double]]
halfspacesMatrix

-- | Returns a point fulfilling a list of inequalities

interiorPoint :: 
  Real a 
  => [Constraint a] -- ^ linear inequalities

  -> IO [Double]    -- ^ point fulfilling the inequaities

interiorPoint :: forall a. Real a => [Constraint a] -> IO [Double]
interiorPoint [Constraint a]
constraints = do 
  let
    constraints' :: [Constraint Rational]
constraints' = (Constraint a -> Constraint Rational)
-> [Constraint a] -> [Constraint Rational]
forall a b. (a -> b) -> [a] -> [b]
map Constraint a -> Constraint Rational
forall a. Real a => Constraint a -> Constraint Rational
toRationalConstraint [Constraint a]
constraints
    halfspacesMatrix :: [[Rational]]
halfspacesMatrix = [Constraint Rational] -> [[Rational]]
forall a. Real a => [Constraint a] -> [[a]]
normalizeConstraints [Constraint Rational]
constraints'
  [Bool]
signs <- [[Rational]] -> IO [Bool]
findSigns [[Rational]]
halfspacesMatrix 
  Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when ([Bool] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Bool]
signs) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
    [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error [Char]
"interiorPoint: no feasible point."
  [[Rational]] -> [Bool] -> IO [Double]
iPoint [[Rational]]
halfspacesMatrix [Bool]
signs