module Numeric.AD.Lagrangian.Internal where
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import qualified Data.Packed.Vector as V
import qualified Data.Packed.Matrix as M
import GHC.IO (unsafePerformIO)
import Numeric.AD
import Numeric.Optimization.Algorithms.HagerZhang05
import Numeric.LinearAlgebra.Algorithms
newtype Constraint = Constraint
{unConstraint :: forall a. (Floating a) => ([a] -> a, a)}
infixr 1 <=>
(<=>) :: (forall a. (Floating a) => [a] -> a)
-> (forall b. (Floating b) => b)
-> Constraint
f <=> c = Constraint (f, c)
minimize :: (forall a. (Floating a) => [a] -> a)
-> [Constraint]
-> Double
-> Int
-> Either (Result, Statistics) (V.Vector Double, V.Vector Double)
minimize f constraints tolerance argCount = result where
(sqGradLgn, gradSqGradLgn) = (fst . g, snd . g) where
g = grad' $ squaredGrad $ lagrangian f constraints argCount
guess = U.replicate (argCount + length constraints) 1
result = case unsafePerformIO $
optimize
(defaultParameters {printFinal = False})
tolerance
guess
(VFunction (sqGradLgn . U.toList))
(VGradient (U.fromList . gradSqGradLgn . U.toList))
Nothing of
(vs, ToleranceStatisfied, _) -> Right (S.take argCount vs,
S.drop argCount vs)
(_, x, y) -> Left (x, y)
maximize :: (forall a. (Floating a) => [a] -> a)
-> [Constraint]
-> Double
-> Int
-> Either (Result, Statistics) (V.Vector Double, V.Vector Double)
maximize f = minimize $ negate . f
lagrangian :: (Floating a)
=> (forall b. (Floating b) => [b] -> b)
-> [Constraint]
-> Int
-> [a]
-> a
lagrangian f constraints argCount argsAndLams = result where
args = take argCount argsAndLams
lams = drop argCount argsAndLams
appliedConstraints = fmap (\(Constraint (g, c)) -> g args c) constraints
result = (f args) + (sum . zipWith (*) lams $ appliedConstraints)
squaredGrad :: (Floating a)
=> (forall b. (Floating b) => [b] -> b)
-> [a]
-> a
squaredGrad f = sum . fmap square . grad f where
square x = x * x
feasible :: (Floating a, Field a, M.Element a)
=> (forall b. (Floating b) => [b] -> b)
->[Constraint]
-> [a]
-> Bool
feasible f constraints points = result where
sqGradLgn :: (Floating a) => [a] -> a
sqGradLgn = squaredGrad $ lagrangian f constraints $ length points
hessianMatrix = M.fromLists . hessian sqGradLgn $ points
result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix