module Numeric.AD.Lagrangian.Internal where
import Numeric.Optimization.Algorithms.HagerZhang05
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import Numeric.AD
import GHC.IO (unsafePerformIO)
import Numeric.AD.Types
import Numeric.AD.Internal.Classes
import Numeric.LinearAlgebra.Algorithms
import qualified Data.Packed.Vector as V
import qualified Data.Packed.Matrix as M
type Constraint a = ([a] -> a, a)
feasible :: (forall a. Floating a => ([a] -> a, [Constraint a], [a]))
-> Bool
feasible params = result where
obj :: Floating a => [a] -> a
obj argsAndLams = squaredGrad lang argsAndLams
lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a)
lang = lagrangian fAndGs (length point)
fAndGs :: (forall a. Floating a => ([a] -> a, [Constraint a]))
fAndGs = (\(x, y, _) -> (x, y)) params
point :: Floating a => [a]
point = (\(_, _, x) -> x) params
h :: [[Double]]
h = hessian obj point
hessianMatrix = M.fromLists h
result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix
solve :: (forall a. Floating a => ([a] -> a, [Constraint a]))
-> Int
-> Either (Result, Statistics) ([Double], [Double])
solve params argCount = result where
obj :: Floating a => [a] -> a
obj argsAndLams = squaredGrad lang argsAndLams
lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a)
lang = lagrangian params argCount
constraintCount = length (snd params)
guess = U.fromList $ replicate (argCount + constraintCount) (1.0 :: Double)
result = case unsafePerformIO (optimize (defaultParameters { printFinal = False })
0.00001 guess (toFunction obj) (toGradient obj)
Nothing) of
(vs, ToleranceStatisfied, _) -> Right (take argCount . S.toList $ vs,
drop argCount . S.toList $ vs)
(_, x, y) -> Left (x, y)
lagrangian :: Floating a
=> ([a] -> a, [Constraint a])
-> Int
-> [a]
-> a
lagrangian (f, constraints) argsLength argsAndLams = result where
result = f args + (sum $ zipWith (*) lams appliedConstraints)
appliedConstraints = map (\(f, c) -> f args c) constraints
args = take argsLength argsAndLams
lams = drop argsLength argsAndLams
sumMap f = sum . map f
squaredGrad :: Num a
=> (forall s. Mode s => [AD s a] -> AD s a) -> [a] -> a
squaredGrad f vs = sumMap (\x -> x*x) (grad f vs)
toFunction :: (forall a. Floating a => [a] -> a) -> Function Simple
toFunction f = VFunction (f . U.toList)
toGradient :: (forall a. Floating a => [a] -> a) -> Gradient Simple
toGradient f = VGradient (U.fromList . grad f . U.toList)