module Math.ExpPairs
( optimize
, OptimizeResult
, optimalValue
, optimalPair
, optimalPath
, simulateOptimize
, simulateOptimize'
, LinearForm
, RationalForm (..)
, IneqType
, Constraint
, InitPair
, Path
, RatioInf (..)
, RationalInf
, pattern K
, pattern L
, pattern M
, (>.), (>=.), (<.), (<=.)
, scaleLF
) where
import Control.Arrow hiding ((<+>))
import Data.Function (on)
import Data.Ord (comparing)
import Data.List (minimumBy)
import Data.Monoid
import Data.Ratio
import Text.PrettyPrint.Leijen hiding ((<$>), (<>))
import qualified Text.PrettyPrint.Leijen as PP
import Text.Printf
import Math.ExpPairs.LinearForm
import Math.ExpPairs.Process
import Math.ExpPairs.Pair
import Math.ExpPairs.RatioInf
pattern K n = LinearForm n 0 0
pattern L n = LinearForm 0 n 0
pattern M n = LinearForm 0 0 n
(>.) :: Num t => LinearForm t -> LinearForm t -> Constraint t
lf1 >. lf2 = Constraint (lf1 lf2) Strict
infix 5 >.
(>=.) :: Num t => LinearForm t -> LinearForm t -> Constraint t
lf1 >=. lf2 = Constraint (lf1 lf2) NonStrict
infix 5 >=.
(<.) :: Num t => LinearForm t -> LinearForm t -> Constraint t
lf1 <. lf2 = Constraint (lf2 lf1) Strict
infix 5 <.
(<=.) :: Num t => LinearForm t -> LinearForm t -> Constraint t
lf1 <=. lf2 = Constraint (lf2 lf1) NonStrict
infix 5 <=.
evalFunctional :: [InitPair] -> [InitPair] -> [RationalForm Rational] -> [Constraint Rational] -> Path -> (RationalInf, InitPair)
evalFunctional corners interiors rfs cons path = case rs of
[] -> (InfPlus, error "evalFunctional: cannot find any exponential pair, which satisfies constraints")
_ -> minimumBy (comparing fst) rs
where
applyPath = map (evalPath path . initPairToProjValue &&& id)
corners' = applyPath corners
interiors' = applyPath interiors
predicate (p, _) = all (checkConstraint p) cons
qs
| all predicate corners' = corners'
| any predicate corners' = filter predicate interiors'
| otherwise = []
rs = map (first $ \p -> maximum (map (evalRF p) rfs)) qs
checkMConstraints :: Path -> [Constraint Rational] -> Bool
checkMConstraints path = all (\con -> any (\p -> checkConstraint (evalPath path p) con) triangleT) where
triangleT = [(0, 1, 1), (0, 1, 2), (1, 1, 2)]
data OptimizeResult = OptimizeResult {
optimalValue :: RationalInf,
optimalPair :: InitPair,
optimalPath :: Path
}
deriving (Show)
instance Pretty OptimizeResult where
pretty (OptimizeResult r' ip p) = pretty1 r' PP.<$>
(parens (pretty (k%m) PP.<> comma PP.<> pretty (l%m)) <+> equals <+> pretty p </> pretty ip)
where
pretty1 r@(Finite rr) = text (printf "%.6f" (fromRational rr :: Double)) <+> equals <+> pretty r
pretty1 r = pretty r
(k, l, m) = evalPath p $ initPairToProjValue ip
instance Eq OptimizeResult where
(==) = (==) `on` optimalValue
instance Ord OptimizeResult where
compare = compare `on` optimalValue
simulateOptimize :: Rational -> OptimizeResult
simulateOptimize r = OptimizeResult (Finite r) Corput01 mempty
simulateOptimize' :: RationalInf -> OptimizeResult
simulateOptimize' r = OptimizeResult r Corput01 mempty
optimize :: [RationalForm Rational] -> [Constraint Rational] -> OptimizeResult
optimize rfs cons = optimize' rfs cons (OptimizeResult r0 ip0 mempty) where
(r0, ip0) = evalFunctional [Corput01, Corput12] [Corput01, Corput12] rfs cons mempty
optimize' :: [RationalForm Rational] -> [Constraint Rational] -> OptimizeResult -> OptimizeResult
optimize' rfs cons ret@(OptimizeResult r _ path)
| lengthPath path > 100 = ret
| otherwise = retBA where
ret0@(OptimizeResult r0 ip0 _) = if r0' < r then OptimizeResult r0' ip0' path else ret where
(r0', ip0') = evalFunctional corners interiors rfs cons path
corners = [Mix 1 0, Mix 0 1, Mix 0 0]
interiors = initPairs
cons0 = if r0==InfPlus then cons else cons ++ map (consBuilder r0) rfs
retA@(OptimizeResult r1 ip1 _) = if checkMConstraints patha cons0 && r1' < r0 then branchA else ret0 where
patha = path <> aPath
branchA@(OptimizeResult r1' _ _) = optimize' rfs cons (OptimizeResult r0 ip0 patha)
cons1 = if r1==r0 then cons0 else cons ++ map (consBuilder r1) rfs
retBA = if checkMConstraints pathba cons1 && r2' < r1 then branchB else retA where
pathba = path <> baPath
branchB@(OptimizeResult r2' _ _) = optimize' rfs cons (OptimizeResult r1 ip1 pathba)
consBuilder rr (num :/: den) = (substituteLF (num, den, 1) (L (toRational rr) K 1)) >. 0