{-|
Module      : LPPaver.Decide.Linearisation
Description : Linearisations for conjunctions
Copyright   : (c) Junaid Rasheed, 2021-2022
License     : MPL
Maintainer  : jrasheed178@gmail.com
Stability   : experimental
Module defining linearisations for conjunctions of 'E.ESafe' terms.
-}
module LPPaver.Decide.Linearisation where

import MixedTypesNumPrelude
import qualified PropaFP.Expression as E
import PropaFP.VarMap
import AERN2.MP
import AERN2.MP.Precision
import AERN2.BoxFun.Type
import AERN2.BoxFun.Box
import qualified Data.PQueue.Prio.Max as Q
import Data.Maybe
import PropaFP.Translators.BoxFun
import Control.Parallel.Strategies
import Data.List
import qualified Data.Map as M
import Linear.Simplex.Simplex
import Linear.Simplex.Util
import qualified Linear.Simplex.Types as LT

import LPPaver.Decide.Util
import LPPaver.Constraint.Type
import LPPaver.Constraint.Util
import qualified AERN2.Linear.Vector.Type as V
import Data.Bifunctor

-- |Remove unsat areas from a conjunction arising from a DNF by weakening the conjunction using 'createConstraintsToRemoveConjunctionUnsatArea'.
-- The resulting linear system is solved and optimised by the two-phase simplex method.
-- If the linear system is infeasible, the entire conjunction was unsatisfiable.
removeConjunctionUnsatAreaWithSimplex 
  :: [(CN MPBall, CN MPBall, Box)]  -- ^ A list of values needed to linearise each term in the conjunction. 
                                    -- In each triple, the first item is the value of the term from the 'extreme' left corner of a 'VarMap', 
                                    -- the second item is the value of the term from the 'extreme' right corner of a 'VarMap', 
                                    -- and the third item are partial derivatives of the term over a 'VarMap'.
  -> VarMap                         -- ^ The VarMap over which we are examining the conjunction.
  -> (Maybe Bool, Maybe VarMap)     -- ^ The result of the simplex method on the resulting linear system.
                                    -- (Just False, Nothing) is returned if the system is infeasible.
                                    -- (Nothing, Just newArea) is returned if the system is feasible: newArea is an optimisation of the given 'VarMap'.
removeConjunctionUnsatAreaWithSimplex :: [(CN MPBall, CN MPBall, Box)]
-> VarMap -> (Maybe Bool, Maybe VarMap)
removeConjunctionUnsatAreaWithSimplex [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives VarMap
varMap =
  case Maybe VarMap
mOptimizedVars of
    Just VarMap
optimizedVars -> (Maybe Bool
forall a. Maybe a
Nothing, VarMap -> Maybe VarMap
forall a. a -> Maybe a
Just VarMap
optimizedVars)
    Maybe VarMap
Nothing            -> (Bool -> Maybe Bool
forall a. a -> Maybe a
Just Bool
False, Maybe VarMap
forall a. Maybe a
Nothing)
  where
    ([PolyConstraint]
simplexSystem, Map String Integer
stringIntVarMap) = [Constraint] -> ([PolyConstraint], Map String Integer)
constraintsToSimplexConstraints ([Constraint] -> ([PolyConstraint], Map String Integer))
-> [Constraint] -> ([PolyConstraint], Map String Integer)
forall a b. (a -> b) -> a -> b
$ [(CN MPBall, CN MPBall, Box)] -> VarMap -> [Constraint]
createConstraintsToRemoveConjunctionUnsatArea [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives VarMap
varMap

    vars :: [String]
vars = ((String, (Rational, Rational)) -> String) -> VarMap -> [String]
forall a b. (a -> b) -> [a] -> [b]
map (String, (Rational, Rational)) -> String
forall a b. (a, b) -> a
fst VarMap
varMap

    mFeasibleSolution :: Maybe (DictionaryForm, [Integer], [Integer], Integer)
mFeasibleSolution = [PolyConstraint]
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
findFeasibleSolution [PolyConstraint]
simplexSystem


    -- Uses objective var to extract optimized values for each variable
    extractSimplexResult :: Maybe (Integer, [(Integer, Rational)]) -> Rational
    extractSimplexResult :: Maybe (Integer, [(Integer, Rational)]) -> Rational
extractSimplexResult Maybe (Integer, [(Integer, Rational)])
maybeResult =
      case Maybe (Integer, [(Integer, Rational)])
maybeResult of
        Just (Integer
optimizedIntVar, [(Integer, Rational)]
result) -> -- optimizedIntVar refers to the objective variable. We extract the value of the objective
                                          -- variable from the result
          case Integer -> [(Integer, Rational)] -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
optimizedIntVar [(Integer, Rational)]
result of
            Just Rational
optimizedVarResult -> Rational
optimizedVarResult
            Maybe Rational
Nothing -> String -> Rational
forall a. HasCallStack => String -> a
error String
"Extracting simplex result after finding feasible solution resulted in an infeasible result. This should not happen."
        Maybe (Integer, [(Integer, Rational)])
Nothing -> String -> Rational
forall a. HasCallStack => String -> a
error String
"Could not optimize feasible system. This should not happen."

    mOptimizedVars :: Maybe VarMap
mOptimizedVars =
      case Maybe (DictionaryForm, [Integer], [Integer], Integer)
mFeasibleSolution of
        Just (DictionaryForm
feasibleSystem, [Integer]
slackVars, [Integer]
artificialVars, Integer
objectiveVar) ->
          VarMap -> Maybe VarMap
forall a. a -> Maybe a
Just (VarMap -> Maybe VarMap) -> VarMap -> Maybe VarMap
forall a b. (a -> b) -> a -> b
$
          (String -> (String, (Rational, Rational))) -> [String] -> VarMap
forall a b. (a -> b) -> [a] -> [b]
map -- Optimize (minimize and maximize) all variables in the varMap
          (\String
var ->
            case String -> Map String Integer -> Maybe Integer
forall k a. Ord k => k -> Map k a -> Maybe a
M.lookup String
var Map String Integer
stringIntVarMap of
              Just Integer
intVar ->
                case String -> VarMap -> Maybe (Rational, Rational)
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup String
var VarMap
varMap of
                  Just (Rational
originalL, Rational
_) -> -- In the simplex system, the original lower bound of each var was shifted to 0. We undo this shift after optimization.
                    (
                      String
var,
                      (
                        Rational
originalL Rational -> Rational -> AddType Rational Rational
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Maybe (Integer, [(Integer, Rational)]) -> Rational
extractSimplexResult (ObjectiveFunction
-> DictionaryForm
-> [Integer]
-> [Integer]
-> Integer
-> Maybe (Integer, [(Integer, Rational)])
optimizeFeasibleSystem ([(Integer, Rational)] -> ObjectiveFunction
LT.Min [(Integer
intVar, Rational
1.0)]) DictionaryForm
feasibleSystem [Integer]
slackVars [Integer]
artificialVars Integer
objectiveVar),
                        Rational
originalL Rational -> Rational -> AddType Rational Rational
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Maybe (Integer, [(Integer, Rational)]) -> Rational
extractSimplexResult (ObjectiveFunction
-> DictionaryForm
-> [Integer]
-> [Integer]
-> Integer
-> Maybe (Integer, [(Integer, Rational)])
optimizeFeasibleSystem ([(Integer, Rational)] -> ObjectiveFunction
LT.Max [(Integer
intVar, Rational
1.0)]) DictionaryForm
feasibleSystem [Integer]
slackVars [Integer]
artificialVars Integer
objectiveVar)
                      )
                    )
                  Maybe (Rational, Rational)
Nothing -> String -> (String, (Rational, Rational))
forall a. HasCallStack => String -> a
error String
"Optimized var not found in original varMap. This should not happen."
              Maybe Integer
Nothing -> String -> (String, (Rational, Rational))
forall a. HasCallStack => String -> a
error String
"Integer version of var not found. This should not happen."
          )
          [String]
vars
        Maybe (DictionaryForm, [Integer], [Integer], Integer)
Nothing -> Maybe VarMap
forall a. Maybe a
Nothing

-- |Find a satisfiable point from a conjunction arising from a DNF by strengthening the conjunction using 'createConstraintsToFindSatSolution'.
-- The resulting linear system is solved by the first phase of the two-phase simplex method.
findConjunctionSatAreaWithSimplex 
  :: [(CN MPBall, CN MPBall, Box)]  -- ^ A list of values needed to linearise each term in the conjunction. 
                                    -- In each triple, the first item is the value of the term from the 'extreme' left corner of a 'VarMap', 
                                    -- the second item is the value of the term from the 'extreme' right corner of a 'VarMap', 
                                    -- and the third item are partial derivatives of the term over a 'VarMap'.
  -> VarMap                         -- ^ The VarMap over which we are examining the conjunction.
  -> Bool                           -- ^ A boolean used to determine which 'extreme' corner to strengthen the conjunction from.
                                    -- If true, linearise from the 'extreme' left corner and vice versa.
  -> Maybe VarMap                   -- ^ The result. If this is Nothing, no satisfiable point was found.
findConjunctionSatAreaWithSimplex :: [(CN MPBall, CN MPBall, Box)] -> VarMap -> Bool -> Maybe VarMap
findConjunctionSatAreaWithSimplex [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives VarMap
varMap Bool
isLeftCorner =
  case Maybe [(Integer, Rational)]
mFeasibleVars of
    Just [(Integer, Rational)]
newPoints ->
      VarMap -> Maybe VarMap
forall a. a -> Maybe a
Just
      (VarMap -> Maybe VarMap) -> VarMap -> Maybe VarMap
forall a b. (a -> b) -> a -> b
$
      (String -> (String, (Rational, Rational))) -> [String] -> VarMap
forall a b. (a -> b) -> [a] -> [b]
map
      (\String
var ->
        case String -> Map String Integer -> Maybe Integer
forall k a. Ord k => k -> Map k a -> Maybe a
M.lookup String
var Map String Integer
stringIntVarMap of
          Just Integer
intVar ->
            case String -> VarMap -> Maybe (Rational, Rational)
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup String
var VarMap
varMap of
              Just (Rational
originalL, Rational
_) -> -- In the simplex system, the original lower bound of each var was shifted to 0. We undo this shift after finding a feasible solution
                (
                  String
var,
                  let feasiblePoint :: AddType Rational Rational
feasiblePoint = Rational
originalL Rational -> Rational -> AddType Rational Rational
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Rational -> Maybe Rational -> Rational
forall a. a -> Maybe a -> a
fromMaybe Rational
0.0 (Integer -> [(Integer, Rational)] -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
intVar [(Integer, Rational)]
newPoints)
                  in (Rational
AddType Rational Rational
feasiblePoint, Rational
AddType Rational Rational
feasiblePoint)
                )
              Maybe (Rational, Rational)
Nothing -> String -> (String, (Rational, Rational))
forall a. HasCallStack => String -> a
error String
"Optimized var not found in original varMap. This should not happen."
          Maybe Integer
Nothing -> String -> (String, (Rational, Rational))
forall a. HasCallStack => String -> a
error String
"Integer version of var not found. This should not happen."
      )
      [String]
vars
    Maybe [(Integer, Rational)]
Nothing -> String -> Maybe VarMap -> Maybe VarMap
forall {p1} {p2}. p1 -> p2 -> p2
trace String
"no sat solution" Maybe VarMap
forall a. Maybe a
Nothing
  where
    ([PolyConstraint]
simplexSystem, Map String Integer
stringIntVarMap) = [Constraint] -> ([PolyConstraint], Map String Integer)
constraintsToSimplexConstraints ([Constraint] -> ([PolyConstraint], Map String Integer))
-> [Constraint] -> ([PolyConstraint], Map String Integer)
forall a b. (a -> b) -> a -> b
$ [(CN MPBall, CN MPBall, Box)] -> VarMap -> Bool -> [Constraint]
createConstraintsToFindSatSolution [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives VarMap
varMap Bool
isLeftCorner

    vars :: [String]
vars = ((String, (Rational, Rational)) -> String) -> VarMap -> [String]
forall a b. (a -> b) -> [a] -> [b]
map (String, (Rational, Rational)) -> String
forall a b. (a, b) -> a
fst VarMap
varMap

    mFeasibleSolution :: Maybe (DictionaryForm, [Integer], [Integer], Integer)
mFeasibleSolution = [PolyConstraint]
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
findFeasibleSolution [PolyConstraint]
simplexSystem

    mFeasibleVars :: Maybe [(Integer, Rational)]
mFeasibleVars =
      case Maybe (DictionaryForm, [Integer], [Integer], Integer)
mFeasibleSolution of
        Just (DictionaryForm
feasibleSystem, [Integer]
_slackVars, [Integer]
_artificialVars, Integer
_objectiveVar) -> [(Integer, Rational)] -> Maybe [(Integer, Rational)]
forall a. a -> Maybe a
Just ([(Integer, Rational)] -> Maybe [(Integer, Rational)])
-> [(Integer, Rational)] -> Maybe [(Integer, Rational)]
forall a b. (a -> b) -> a -> b
$ DictionaryForm -> [(Integer, Rational)]
displayDictionaryResults DictionaryForm
feasibleSystem
        Maybe (DictionaryForm, [Integer], [Integer], Integer)
Nothing -> Maybe [(Integer, Rational)]
forall a. Maybe a
Nothing

-- |Linearisations that weaken a conjunction of terms over some box.
createConstraintsToRemoveConjunctionUnsatArea 
  :: [(CN MPBall, CN MPBall, Box)]  -- ^ A list of values needed to linearise each term in the conjunction. 
                                    -- In each triple, the first item is the value of the term from the 'extreme' left corner of a 'VarMap', 
                                    -- the second item is the value of the term from the 'extreme' right corner of a 'VarMap', 
                                    -- and the third item are partial derivatives of the term over a 'VarMap'.
  -> VarMap                         -- ^ The VarMap over which we are examining the conjunction.
  -> [Constraint]                   -- ^ An implicit linear system that is a weakening of the conjunction.
createConstraintsToRemoveConjunctionUnsatArea :: [(CN MPBall, CN MPBall, Box)] -> VarMap -> [Constraint]
createConstraintsToRemoveConjunctionUnsatArea [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives VarMap
varMap =
  [Constraint]
domainConstraints [Constraint] -> [Constraint] -> [Constraint]
forall a. [a] -> [a] -> [a]
++ [Constraint]
functionConstraints
  where
    vars :: [String]
vars = ((String, (Rational, Rational)) -> String) -> VarMap -> [String]
forall a b. (a -> b) -> [a] -> [b]
map (String, (Rational, Rational)) -> String
forall a b. (a, b) -> a
fst VarMap
varMap
    varsNewUpperBounds :: [Rational]
varsNewUpperBounds = ((String, (Rational, Rational)) -> Rational)
-> VarMap -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (\(String
_, (Rational
l, Rational
r)) -> Rational
r Rational -> Rational -> SubType Rational Rational
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Rational
l) VarMap
varMap

    -- var >= varLower - varLower && var <= varUpper - varLower
    -- Since var >= 0 is assumed by the simplex method, var >= varLower - varLower is not needed
    domainConstraints :: [Constraint]
domainConstraints =
      ((String, (Rational, Rational)) -> Constraint)
-> VarMap -> [Constraint]
forall a b. (a -> b) -> [a] -> [b]
map
      (\(String
var, (Rational
varLower, Rational
varUpper)) ->
        [(String, Rational)] -> Rational -> Constraint
LEQ [(String
var, Rational
1.0)] (Rational -> Constraint) -> Rational -> Constraint
forall a b. (a -> b) -> a -> b
$ Rational
varUpper Rational -> Rational -> SubType Rational Rational
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Rational
varLower
      )
      VarMap
varMap

    -- The following constraints in this variable are...
    --  fn - (fnx1GradientR * x1) - .. - (fnxnGradientR * xn) <= fnLeftCorner + (fnx1GradientR * -x1L) + .. + (fnxnGradientR * -xnL)
    --  fn - (fnx1GradientL * x1) - .. - (fnxnGradientL * xn) <= fnRightCorner + (-fnx1GradientL * x1R) + .. + (-fnxnGradientL * xnR)
    -- and these are equivalent to...
    --  fn <= fnLeftCorner  + (fnx1GradientR  * (x1 - x1L)) + .. + (fnxnGradientR * (xn - xnL))
    --  fn <= fnRightCorner + (-fnx1GradientL * (x1R - x1)) + .. + (-fnxnGradientL * (xnR - xn))
    -- 
    functionConstraints :: [Constraint]
functionConstraints =
      ((Integer, (CN MPBall, CN MPBall, Box)) -> [Constraint])
-> [(Integer, (CN MPBall, CN MPBall, Box))] -> [Constraint]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap
      (\(Integer
fnInt, (CN MPBall
fLeftRange, CN MPBall
fRightRange, Box
fPartialDerivatives)) ->
        let
          fNegatedPartialDerivativesLowerBounds :: [Rational]
fNegatedPartialDerivativesLowerBounds = (CN MPBall -> Rational) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational
forall t. CanNeg t => t -> NegType t
negate (Rational -> Rational)
-> (CN MPBall -> Rational) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational, Rational) -> Rational
forall a b. (a, b) -> a
fst ((Rational, Rational) -> Rational)
-> (CN MPBall -> (Rational, Rational)) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CN MPBall -> (Rational, Rational)
mpBallToRational) ([CN MPBall] -> [Rational]) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> a -> b
$ Box -> [CN MPBall]
forall a. Vector a -> [a]
V.toList Box
fPartialDerivatives
          fNegatedPartialDerivativesUpperBounds :: [Rational]
fNegatedPartialDerivativesUpperBounds = (CN MPBall -> Rational) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational
forall t. CanNeg t => t -> NegType t
negate (Rational -> Rational)
-> (CN MPBall -> Rational) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational, Rational) -> Rational
forall a b. (a, b) -> b
snd ((Rational, Rational) -> Rational)
-> (CN MPBall -> (Rational, Rational)) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CN MPBall -> (Rational, Rational)
mpBallToRational) ([CN MPBall] -> [Rational]) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> a -> b
$ Box -> [CN MPBall]
forall a. Vector a -> [a]
V.toList Box
fPartialDerivatives
          fLeftUpperBound :: Rational
fLeftUpperBound = (Rational, Rational) -> Rational
forall a b. (a, b) -> b
snd ((Rational, Rational) -> Rational)
-> (Rational, Rational) -> Rational
forall a b. (a -> b) -> a -> b
$ CN MPBall -> (Rational, Rational)
mpBallToRational CN MPBall
fLeftRange
          fRightUpperBound :: Rational
fRightUpperBound = (Rational, Rational) -> Rational
forall a b. (a, b) -> b
snd ((Rational, Rational) -> Rational)
-> (Rational, Rational) -> Rational
forall a b. (a -> b) -> a -> b
$ CN MPBall -> (Rational, Rational)
mpBallToRational CN MPBall
fRightRange

          --FIXME: make this safe. Check if vars contain ^fSimplex[0-9]+$. If so, try ^fSimplex1[0-9]+$ or something
          fSimplexN :: String
fSimplexN = String
"fSimplex" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
fnInt
        in
          [
            -- f is definitely below this line from the left corner
            -- Multiplication with left corner and partial derivatives omitted, since left corner is zero
            [(String, Rational)] -> Rational -> Constraint
LEQ ((String
fSimplexN, Rational
1.0) (String, Rational) -> [(String, Rational)] -> [(String, Rational)]
forall a. a -> [a] -> [a]
: [String] -> [Rational] -> [(String, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [String]
vars [Rational]
fNegatedPartialDerivativesUpperBounds)
              Rational
fLeftUpperBound,
            -- f is definitely below this line from the right corner
            [(String, Rational)] -> Rational -> Constraint
LEQ ((String
fSimplexN, Rational
1.0) (String, Rational) -> [(String, Rational)] -> [(String, Rational)]
forall a. a -> [a] -> [a]
: [String] -> [Rational] -> [(String, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [String]
vars [Rational]
fNegatedPartialDerivativesLowerBounds)
              (Rational -> Constraint) -> Rational -> Constraint
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational -> Rational)
-> Rational -> [Rational] -> Rational
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Rational -> Rational -> Rational
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
add Rational
fRightUpperBound ([Rational] -> Rational) -> [Rational] -> Rational
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational -> Rational)
-> [Rational] -> [Rational] -> [Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Rational -> Rational -> Rational
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
mul [Rational]
varsNewUpperBounds [Rational]
fNegatedPartialDerivativesLowerBounds
          ]
      )
      ([(Integer, (CN MPBall, CN MPBall, Box))] -> [Constraint])
-> [(Integer, (CN MPBall, CN MPBall, Box))] -> [Constraint]
forall a b. (a -> b) -> a -> b
$
      [Integer]
-> [(CN MPBall, CN MPBall, Box)]
-> [(Integer, (CN MPBall, CN MPBall, Box))]
forall a b. [a] -> [b] -> [(a, b)]
zip
      [Integer
1..]
      [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives

    mpBallToRational :: CN MPBall -> (Rational, Rational)
    mpBallToRational :: CN MPBall -> (Rational, Rational)
mpBallToRational = (MPFloat -> Rational)
-> (MPFloat -> Rational)
-> (MPFloat, MPFloat)
-> (Rational, Rational)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap MPFloat -> Rational
forall t. CanBeRational t => t -> Rational
rational MPFloat -> Rational
forall t. CanBeRational t => t -> Rational
rational ((MPFloat, MPFloat) -> (Rational, Rational))
-> (CN MPBall -> (MPFloat, MPFloat))
-> CN MPBall
-> (Rational, Rational)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MPBall -> (MPFloat, MPFloat)
forall i.
IsInterval i =>
i -> (IntervalEndpoint i, IntervalEndpoint i)
endpoints (MPBall -> (MPFloat, MPFloat))
-> (CN MPBall -> MPBall) -> CN MPBall -> (MPFloat, MPFloat)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MPBall -> MPBall
reducePrecionIfInaccurate (MPBall -> MPBall) -> (CN MPBall -> MPBall) -> CN MPBall -> MPBall
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CN MPBall -> MPBall
forall p. CN p -> p
unCN

-- |Linearisations that strengthen a conjunction of terms over some box.
createConstraintsToFindSatSolution
  :: [(CN MPBall, CN MPBall, Box)]  -- ^ A list of values needed to linearise each term in the conjunction. 
                                    -- In each triple, the first item is the value of the term from the 'extreme' left corner of a 'VarMap', 
                                    -- the second item is the value of the term from the 'extreme' right corner of a 'VarMap', 
                                    -- and the third item are partial derivatives of the term over a 'VarMap'.
  -> VarMap                         -- ^ The VarMap over which we are examining the conjunction.
  -> Bool                           -- ^ A boolean used to determine which 'extreme' corner to strengthen the conjunction from.
                                    -- If true, linearise from the 'extreme' left corner and vice versa.
  -> [Constraint]                   -- ^ An implicit linear system that is a weakening of the conjunction.
createConstraintsToFindSatSolution :: [(CN MPBall, CN MPBall, Box)] -> VarMap -> Bool -> [Constraint]
createConstraintsToFindSatSolution [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives VarMap
varMap Bool
isLeftCorner =
  [Constraint]
domainConstraints [Constraint] -> [Constraint] -> [Constraint]
forall a. [a] -> [a] -> [a]
++ [Constraint]
functionConstraints
  where
    vars :: [String]
vars = ((String, (Rational, Rational)) -> String) -> VarMap -> [String]
forall a b. (a -> b) -> [a] -> [b]
map (String, (Rational, Rational)) -> String
forall a b. (a, b) -> a
fst VarMap
varMap
    varsNewUpperBounds :: [Rational]
varsNewUpperBounds = ((String, (Rational, Rational)) -> Rational)
-> VarMap -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (\(String
_, (Rational
l, Rational
r)) -> Rational
r Rational -> Rational -> SubType Rational Rational
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Rational
l) VarMap
varMap

    -- var >= varLower - varLower && var <= varUpper - varLower
    -- Since var >= 0 is assumed by the simplex method, var >= varLower - varLower is not needed
    domainConstraints :: [Constraint]
domainConstraints =
      ((String, (Rational, Rational)) -> Constraint)
-> VarMap -> [Constraint]
forall a b. (a -> b) -> [a] -> [b]
map
      (\(String
var, (Rational
varLower, Rational
varUpper)) ->
        [(String, Rational)] -> Rational -> Constraint
LEQ [(String
var, Rational
1.0)] (Rational -> Constraint) -> Rational -> Constraint
forall a b. (a -> b) -> a -> b
$ Rational
varUpper Rational -> Rational -> SubType Rational Rational
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Rational
varLower
      )
      VarMap
varMap

    -- The following constraints in this variable are...
    --  Left  corner: fn - (fnx1GradientL * x1) - .. - (fnxnGradientL * xn) <= fnLeftCorner  - (fnx1GradientL * x1L) - .. - (fnx1GradientL * x1L)
    --  Right corner: fn - (fnx1GradientR * x1) - .. - (fnxnGradientR * xn) <= fnRightCorner - (fnx1GradientR * x1R) - .. - (fnx1GradientR * x1R)
    -- and these are equivalent to...
    --  Left  corner: fn <= fnLeftCorner + (fnx1GradientL * (x1 - x1L)) + .. + (fnxnGradientL * (xn - xnL))
    --  Right corner: fn <= ynRightCorner + (-fnx1GradientR * (xR - x)) + .. + (-fnxnGradientR * (xnR - xn))
    functionConstraints :: [Constraint]
functionConstraints =
      (Integer -> (CN MPBall, CN MPBall, Box) -> Constraint)
-> [Integer] -> [(CN MPBall, CN MPBall, Box)] -> [Constraint]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
      (((Integer, (CN MPBall, CN MPBall, Box)) -> Constraint)
-> Integer -> (CN MPBall, CN MPBall, Box) -> Constraint
forall a b c. ((a, b) -> c) -> a -> b -> c
curry
        (\(Integer
fnInt, (CN MPBall
fLeftRange, CN MPBall
fRightRange, Box
fPartialDerivatives)) ->
          let
            fNegatedPartialDerivativesLowerBounds :: [Rational]
fNegatedPartialDerivativesLowerBounds = (CN MPBall -> Rational) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational
forall t. CanNeg t => t -> NegType t
negate (Rational -> Rational)
-> (CN MPBall -> Rational) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational, Rational) -> Rational
forall a b. (a, b) -> a
fst ((Rational, Rational) -> Rational)
-> (CN MPBall -> (Rational, Rational)) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CN MPBall -> (Rational, Rational)
mpBallToRational) ([CN MPBall] -> [Rational]) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> a -> b
$ Box -> [CN MPBall]
forall a. Vector a -> [a]
V.toList Box
fPartialDerivatives
            fNegatedPartialDerivativesUpperBounds :: [Rational]
fNegatedPartialDerivativesUpperBounds = (CN MPBall -> Rational) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational
forall t. CanNeg t => t -> NegType t
negate (Rational -> Rational)
-> (CN MPBall -> Rational) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational, Rational) -> Rational
forall a b. (a, b) -> b
snd ((Rational, Rational) -> Rational)
-> (CN MPBall -> (Rational, Rational)) -> CN MPBall -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CN MPBall -> (Rational, Rational)
mpBallToRational) ([CN MPBall] -> [Rational]) -> [CN MPBall] -> [Rational]
forall a b. (a -> b) -> a -> b
$ Box -> [CN MPBall]
forall a. Vector a -> [a]
V.toList Box
fPartialDerivatives
            fLeftLowerBound :: Rational
fLeftLowerBound = (Rational, Rational) -> Rational
forall a b. (a, b) -> a
fst ((Rational, Rational) -> Rational)
-> (Rational, Rational) -> Rational
forall a b. (a -> b) -> a -> b
$ CN MPBall -> (Rational, Rational)
mpBallToRational CN MPBall
fLeftRange
            fRightLowerBound :: Rational
fRightLowerBound = (Rational, Rational) -> Rational
forall a b. (a, b) -> a
fst ((Rational, Rational) -> Rational)
-> (Rational, Rational) -> Rational
forall a b. (a -> b) -> a -> b
$ CN MPBall -> (Rational, Rational)
mpBallToRational CN MPBall
fRightRange

            --FIXME: make this safe. Check if vars contain ^fSimplex[0-9]+$. If so, try ^fSimplex1[0-9]+$ or something
            fSimplexN :: String
fSimplexN = String
"fSimplex" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
fnInt
          in
            if Bool
isLeftCorner
              then
                -- f is definitely above this line from the left corner
                -- Multiplication with left corner and partial derivatives omitted, since left corner is zero
                [(String, Rational)] -> Rational -> Constraint
LEQ ((String
fSimplexN, Rational
1.0) (String, Rational) -> [(String, Rational)] -> [(String, Rational)]
forall a. a -> [a] -> [a]
: [String] -> [Rational] -> [(String, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [String]
vars [Rational]
fNegatedPartialDerivativesLowerBounds)
                  Rational
fLeftLowerBound
              else
                -- f is definitely above this line from the right corner
                [(String, Rational)] -> Rational -> Constraint
LEQ ((String
fSimplexN, Rational
1.0) (String, Rational) -> [(String, Rational)] -> [(String, Rational)]
forall a. a -> [a] -> [a]
: [String] -> [Rational] -> [(String, Rational)]
forall a b. [a] -> [b] -> [(a, b)]
zip [String]
vars [Rational]
fNegatedPartialDerivativesUpperBounds)
                  (Rational -> Constraint) -> Rational -> Constraint
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational -> Rational)
-> Rational -> [Rational] -> Rational
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Rational -> Rational -> Rational
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
add Rational
fRightLowerBound ([Rational] -> Rational) -> [Rational] -> Rational
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational -> Rational)
-> [Rational] -> [Rational] -> [Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Rational -> Rational -> Rational
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
mul [Rational]
varsNewUpperBounds [Rational]
fNegatedPartialDerivativesUpperBounds

        )
      )
      [Integer
1..]
      [(CN MPBall, CN MPBall, Box)]
cornerValuesWithDerivatives

    mpBallToRational :: CN MPBall -> (Rational, Rational)
    mpBallToRational :: CN MPBall -> (Rational, Rational)
mpBallToRational = (MPFloat -> Rational)
-> (MPFloat -> Rational)
-> (MPFloat, MPFloat)
-> (Rational, Rational)
forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap MPFloat -> Rational
forall t. CanBeRational t => t -> Rational
rational MPFloat -> Rational
forall t. CanBeRational t => t -> Rational
rational ((MPFloat, MPFloat) -> (Rational, Rational))
-> (CN MPBall -> (MPFloat, MPFloat))
-> CN MPBall
-> (Rational, Rational)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MPBall -> (MPFloat, MPFloat)
forall i.
IsInterval i =>
i -> (IntervalEndpoint i, IntervalEndpoint i)
endpoints (MPBall -> (MPFloat, MPFloat))
-> (CN MPBall -> MPBall) -> CN MPBall -> (MPFloat, MPFloat)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MPBall -> MPBall
reducePrecionIfInaccurate (MPBall -> MPBall) -> (CN MPBall -> MPBall) -> CN MPBall -> MPBall
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CN MPBall -> MPBall
forall p. CN p -> p
unCN