-- SPDX-License-Identifier: MPL-2.0

{- | Bounded CPU reference local solvers. Dense quasi-Newton updates implement
their named secant formulas; L-BFGS uses a separate two-loop recursion.
Scalar Brent combines safeguarded parabolic interpolation and golden steps.
Newton-CG uses Hessian products and detects nonpositive curvature.
-}
module HMath.Optimization.Local (
    ScalarMethod (..),
    minimizeScalar,
    QuasiNewton (..),
    minimizeQuasiNewton,
    NewtonMethod (..),
    minimizeNewton,
    nonlinearCG,
    ResidualFunction,
    JacobianFunction,
    LeastSquaresLoss (..),
    leastSquares,
    nelderMead,
    powell,
) where

import Control.Monad (foldM)
import qualified Data.Foldable as Foldable (foldl')
import Data.List (sortOn)
import HMath.Error
import HMath.Optimization.Common

-- | Bounded one-dimensional minimization strategy (golden section,Brent,or ternary search).
data ScalarMethod = GoldenSection | TernarySearch | BrentSearch deriving (ScalarMethod -> ScalarMethod -> Bool
(ScalarMethod -> ScalarMethod -> Bool)
-> (ScalarMethod -> ScalarMethod -> Bool) -> Eq ScalarMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: ScalarMethod -> ScalarMethod -> Bool
== :: ScalarMethod -> ScalarMethod -> Bool
$c/= :: ScalarMethod -> ScalarMethod -> Bool
/= :: ScalarMethod -> ScalarMethod -> Bool
Eq, Int -> ScalarMethod -> ShowS
[ScalarMethod] -> ShowS
ScalarMethod -> String
(Int -> ScalarMethod -> ShowS)
-> (ScalarMethod -> String)
-> ([ScalarMethod] -> ShowS)
-> Show ScalarMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> ScalarMethod -> ShowS
showsPrec :: Int -> ScalarMethod -> ShowS
$cshow :: ScalarMethod -> String
show :: ScalarMethod -> String
$cshowList :: [ScalarMethod] -> ShowS
showList :: [ScalarMethod] -> ShowS
Show)

-- | Minimize on a finite bracket a<b. Result residual is interval width, not an unconstrained derivative.
minimizeScalar ::
    ScalarMethod ->
    OptimizeOptions ->
    (Double -> Either MathError Double) ->
    Double ->
    Double ->
    Either MathError OptimizeResult
minimizeScalar :: ScalarMethod
-> OptimizeOptions
-> (Double -> Either MathError Double)
-> Double
-> Double
-> Either MathError OptimizeResult
minimizeScalar ScalarMethod
method OptimizeOptions
o Double -> Either MathError Double
f Double
a0 Double
b0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    [Double] -> Either MathError ()
checkedPoint [Double
a0, Double
b0]
    if Double
a0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
b0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"scalar minimization requires a < b") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    case ScalarMethod
method of
        ScalarMethod
GoldenSection -> do
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o Int
2
            let c :: Double
c = Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ratio Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a0); d :: Double
d = Double
a0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ratio Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a0)
            fc <- Double -> Either MathError Double
eval Double
c
            fd <- eval d
            golden 0 2 a0 b0 c d fc fd
        ScalarMethod
TernarySearch -> Int -> Int -> Double -> Double -> Either MathError OptimizeResult
ternary Int
0 Int
0 Double
a0 Double
b0
        ScalarMethod
BrentSearch -> do
            let x :: Double
x = Double
a0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ratio) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a0)
            y <- Double -> Either MathError Double
eval Double
x
            brent 0 1 a0 b0 x x x y y y 0 0
  where
    ratio :: Double
ratio = (Double -> Double
forall a. Floating a => a -> a
sqrt Double
5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
    eval :: Double -> Either MathError Double
eval Double
x = do [Double] -> Either MathError ()
checkedPoint [Double
x]; y <- Double -> Either MathError Double
f Double
x; checkedPoint [y]; Right y
    done :: Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either a OptimizeResult
done Int
k Int
n Double
a Double
b Double
x Double
y = OptimizeResult -> Either a OptimizeResult
forall a b. b -> Either a b
Right ([Double]
-> Double -> Double -> Int -> Int -> Bool -> OptimizeResult
OptimizeResult [Double
x] Double
y (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Int
k Int
n Bool
True)
    golden :: Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Either MathError OptimizeResult
golden Int
k Int
n Double
a Double
b Double
c Double
d Double
fc Double
fd
        | Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) =
            if Double
fc Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
fd then Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either MathError OptimizeResult
forall {a}.
Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either a OptimizeResult
done Int
k Int
n Double
a Double
b Double
c Double
fc else Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either MathError OptimizeResult
forall {a}.
Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either a OptimizeResult
done Int
k Int
n Double
a Double
b Double
d Double
fd
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o = MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"golden section" Int
k)
        | Bool
otherwise = do
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
            if Double
fc Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
fd
                then do
                    let z :: Double
z = Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ratio Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)
                    fz <- Double -> Either MathError Double
eval Double
z
                    golden (k + 1) (n + 1) a d z c fz fc
                else do
                    let z :: Double
z = Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ratio Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c)
                    fz <- Double -> Either MathError Double
eval Double
z
                    golden (k + 1) (n + 1) c b d z fd fz
    ternary :: Int -> Int -> Double -> Double -> Either MathError OptimizeResult
ternary Int
k Int
n Double
a Double
b
        | Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) = do OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1); y <- Double -> Either MathError Double
eval ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2); done k (n + 1) a b ((a + b) / 2) y
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o = MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"ternary search" Int
k)
        | Bool
otherwise = do
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
            let c :: Double
c = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3; d :: Double
d = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
            fc <- Double -> Either MathError Double
eval Double
c
            fd <- eval d
            if fc < fd then ternary (k + 1) (n + 2) a d else ternary (k + 1) (n + 2) c b
    brent :: Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Either MathError OptimizeResult
brent Int
k Int
n Double
a Double
b Double
x Double
w Double
v Double
fx Double
fw Double
fv Double
step Double
old
        | Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs Double
x) = Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either MathError OptimizeResult
forall {a}.
Int
-> Int
-> Double
-> Double
-> Double
-> Double
-> Either a OptimizeResult
done Int
k Int
n Double
a Double
b Double
x Double
fx
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o = MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Brent minimum" Int
k)
        | Bool
otherwise = do
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
            let mid :: Double
mid = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                tol :: Double
tol = OptimizeOptions -> Double
optTolerance OptimizeOptions
o Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4
                r :: Double
r = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fx Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fv)
                q0 :: Double
q0 = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
fx Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fw)
                p0 :: Double
p0 = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
q0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r
                q1 :: Double
q1 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
q0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
r)
                p :: Double
p = if Double
q1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 then -Double
p0 else Double
p0
                q :: Double
q = Double -> Double
forall a. Num a => a -> a
abs Double
q1
                parabolic :: Bool
parabolic = Double -> Double
forall a. Num a => a -> a
abs Double
old Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
tol Bool -> Bool -> Bool
&& Double
q Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double -> Double
forall a. Num a => a -> a
abs (Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
old Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
                movement :: Double
movement = if Bool
parabolic then Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
q else (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ratio) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
mid then Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x else Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
                newOld :: Double
newOld = if Bool
parabolic then Double
step else if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
mid then Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x else Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x
                d :: Double
d = if Double -> Double
forall a. Num a => a -> a
abs Double
movement Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
tol then Double
movement else if Double
movement Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 then Double
tol else -Double
tol
                u :: Double
u = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d
            fu <- Double -> Either MathError Double
eval Double
u
            let (a', b') = if fu <= fx then if u < x then (a, x) else (x, b) else if u < x then (u, b) else (a, u)
            if fu <= fx
                then brent (k + 1) (n + 1) a' b' u x w fu fx fw d newOld
                else
                    if fu <= fw || w == x
                        then brent (k + 1) (n + 1) a' b' x u v fx fu fv d newOld
                        else
                            if fu <= fv || v == x || v == w
                                then brent (k + 1) (n + 1) a' b' x w u fx fw fu d newOld
                                else brent (k + 1) (n + 1) a' b' x w v fx fw fv d newOld

-- | Quasi-Newton update selector such as BFGS,DFP,SR1,or limited-memory BFGS.
data QuasiNewton = BFGS | LBFGS Int | DFP | SR1 deriving (QuasiNewton -> QuasiNewton -> Bool
(QuasiNewton -> QuasiNewton -> Bool)
-> (QuasiNewton -> QuasiNewton -> Bool) -> Eq QuasiNewton
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: QuasiNewton -> QuasiNewton -> Bool
== :: QuasiNewton -> QuasiNewton -> Bool
$c/= :: QuasiNewton -> QuasiNewton -> Bool
/= :: QuasiNewton -> QuasiNewton -> Bool
Eq, Int -> QuasiNewton -> ShowS
[QuasiNewton] -> ShowS
QuasiNewton -> String
(Int -> QuasiNewton -> ShowS)
-> (QuasiNewton -> String)
-> ([QuasiNewton] -> ShowS)
-> Show QuasiNewton
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> QuasiNewton -> ShowS
showsPrec :: Int -> QuasiNewton -> ShowS
$cshow :: QuasiNewton -> String
show :: QuasiNewton -> String
$cshowList :: [QuasiNewton] -> ShowS
showList :: [QuasiNewton] -> ShowS
Show)

-- | Run the named inverse-Hessian secant update with safeguarded curvature and a strong-Wolfe descent search.
minimizeQuasiNewton ::
    QuasiNewton ->
    OptimizeOptions ->
    Objective ->
    Gradient ->
    [Double] ->
    Either MathError OptimizeResult
minimizeQuasiNewton :: QuasiNewton
-> OptimizeOptions
-> Objective
-> Gradient
-> [Double]
-> Either MathError OptimizeResult
minimizeQuasiNewton QuasiNewton
method OptimizeOptions
o Objective
f Gradient
g [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    case QuasiNewton
method of
        LBFGS Int
m | Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 -> MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"positive L-BFGS memory required")
        QuasiNewton
_ -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o Int
2
    y <- Objective -> Objective
checkedValue Objective
f [Double]
x0
    gs <- checkedGradient g x0
    loop 0 2 x0 y gs (identity (length x0)) []
  where
    loop :: Int
-> Int
-> [Double]
-> Double
-> [Double]
-> [[Double]]
-> [([Double], [Double], Double)]
-> Either MathError OptimizeResult
loop Int
k Int
n [Double]
x Double
y [Double]
gs [[Double]]
h [([Double], [Double], Double)]
hist
        | [Double] -> Double
norm [Double]
gs Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o = OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right ([Double]
-> Double -> Double -> Int -> Int -> Bool -> OptimizeResult
OptimizeResult [Double]
x Double
y ([Double] -> Double
norm [Double]
gs) Int
k Int
n Bool
True)
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o = MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence (QuasiNewton -> String
forall a. Show a => a -> String
show QuasiNewton
method) Int
k)
        | Bool
otherwise = do
            let direction0 :: [Double]
direction0 = case QuasiNewton
method of LBFGS Int
_ -> Double -> [Double] -> [Double]
scale (-Double
1) ([([Double], [Double], Double)] -> [Double] -> [Double]
twoLoop [([Double], [Double], Double)]
hist [Double]
gs); QuasiNewton
_ -> Double -> [Double] -> [Double]
scale (-Double
1) ([[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
gs)
                direction :: [Double]
direction = if [Double] -> [Double] -> Double
dot [Double]
direction0 [Double]
gs Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then [Double]
direction0 else Double -> [Double] -> [Double]
scale (-Double
1) [Double]
gs
            (alpha, yn, gn, cost) <-
                LineSearch
-> Int
-> Objective
-> Gradient
-> [Double]
-> [Double]
-> Double
-> [Double]
-> Either MathError (Double, Double, [Double], Int)
lineSearch LineSearch
StrongWolfe (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
80 (OptimizeOptions -> Int
optMaxEvaluations OptimizeOptions
o Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n)) Objective
f Gradient
g [Double]
x [Double]
direction Double
y [Double]
gs
            checkBudget o (n + cost)
            let xn = [Double] -> [Double] -> [Double]
add [Double]
x (Double -> [Double] -> [Double]
scale Double
alpha [Double]
direction)
                s = [Double] -> [Double] -> [Double]
sub [Double]
xn [Double]
x
                z = [Double] -> [Double] -> [Double]
sub [Double]
gn [Double]
gs
                sz = [Double] -> [Double] -> Double
dot [Double]
s [Double]
z
                hz = [[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
z
                zhz = [Double] -> [Double] -> Double
dot [Double]
z [Double]
hz
                r = [Double] -> [Double] -> [Double]
sub [Double]
s [Double]
hz
                update [[c]]
a [[c]]
b = ([c] -> [c] -> [c]) -> [[c]] -> [[c]] -> [[c]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((c -> c -> c) -> [c] -> [c] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith c -> c -> c
forall a. Num a => a -> a -> a
(+)) [[c]]
a [[c]]
b
                minus [[c]]
a [[c]]
b = ([c] -> [c] -> [c]) -> [[c]] -> [[c]] -> [[c]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((c -> c -> c) -> [c] -> [c] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-)) [[c]]
a [[c]]
b
                h' = case QuasiNewton
method of
                    QuasiNewton
BFGS
                        | Double
sz Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
z ->
                            let v :: [[Double]]
v = [[Double]] -> [[Double]] -> [[Double]]
forall {c}. Num c => [[c]] -> [[c]] -> [[c]]
minus (Int -> [[Double]]
identity ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x)) (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sz)) ([Double] -> [Double] -> [[Double]]
outer [Double]
s [Double]
z))
                             in [[Double]] -> [[Double]] -> [[Double]]
forall {c}. Num c => [[c]] -> [[c]] -> [[c]]
update ([[Double]] -> [[Double]] -> [[Double]]
matMul ([[Double]] -> [[Double]] -> [[Double]]
matMul [[Double]]
v [[Double]]
h) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
v)) (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sz)) ([Double] -> [Double] -> [[Double]]
outer [Double]
s [Double]
s))
                    QuasiNewton
DFP
                        | Double
sz Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
z Bool -> Bool -> Bool
&& Double
zhz Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 ->
                            [[Double]] -> [[Double]] -> [[Double]]
forall {c}. Num c => [[c]] -> [[c]] -> [[c]]
update ([[Double]] -> [[Double]] -> [[Double]]
forall {c}. Num c => [[c]] -> [[c]] -> [[c]]
minus [[Double]]
h (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
zhz)) ([Double] -> [Double] -> [[Double]]
outer [Double]
hz [Double]
hz))) (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sz)) ([Double] -> [Double] -> [[Double]]
outer [Double]
s [Double]
s))
                    QuasiNewton
SR1 | Double -> Double
forall a. Num a => a -> a
abs ([Double] -> [Double] -> Double
dot [Double]
r [Double]
z) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-8 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
z -> [[Double]] -> [[Double]] -> [[Double]]
forall {c}. Num c => [[c]] -> [[c]] -> [[c]]
update [[Double]]
h (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> [Double] -> Double
dot [Double]
r [Double]
z)) ([Double] -> [Double] -> [[Double]]
outer [Double]
r [Double]
r))
                    QuasiNewton
_ -> [[Double]]
h
                hist' = case QuasiNewton
method of LBFGS Int
m | Double
sz Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
norm [Double]
z -> Int
-> [([Double], [Double], Double)] -> [([Double], [Double], Double)]
forall a. Int -> [a] -> [a]
take Int
m (([Double]
s, [Double]
z, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sz) ([Double], [Double], Double)
-> [([Double], [Double], Double)] -> [([Double], [Double], Double)]
forall a. a -> [a] -> [a]
: [([Double], [Double], Double)]
hist); QuasiNewton
_ -> [([Double], [Double], Double)]
hist
            loop (k + 1) (n + cost) xn yn gn h' hist'
    twoLoop :: [([Double], [Double], Double)] -> [Double] -> [Double]
twoLoop [] [Double]
gs = [Double]
gs
    twoLoop hs :: [([Double], [Double], Double)]
hs@(([Double]
s0, [Double]
y0, Double
_) : [([Double], [Double], Double)]
_) [Double]
gs =
        let ([Double]
q, [Double]
alphas) =
                (([Double], [Double])
 -> ([Double], [Double], Double) -> ([Double], [Double]))
-> ([Double], [Double])
-> [([Double], [Double], Double)]
-> ([Double], [Double])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl'
                    (\([Double]
v, [Double]
as) ([Double]
s, [Double]
y, Double
rho) -> let a :: Double
a = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dot [Double]
s [Double]
v in ([Double] -> [Double] -> [Double]
sub [Double]
v (Double -> [Double] -> [Double]
scale Double
a [Double]
y), [Double]
as [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
a]))
                    ([Double]
gs, [])
                    [([Double], [Double], Double)]
hs
            gamma :: Double
gamma = [Double] -> [Double] -> Double
dot [Double]
s0 [Double]
y0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> [Double] -> Double
dot [Double]
y0 [Double]
y0
         in ([Double] -> (([Double], [Double], Double), Double) -> [Double])
-> [Double] -> [(([Double], [Double], Double), Double)] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl'
                (\[Double]
v (([Double]
s, [Double]
y, Double
rho), Double
a) -> [Double] -> [Double] -> [Double]
add [Double]
v (Double -> [Double] -> [Double]
scale (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dot [Double]
y [Double]
v) [Double]
s))
                (Double -> [Double] -> [Double]
scale Double
gamma [Double]
q)
                ([(([Double], [Double], Double), Double)]
-> [(([Double], [Double], Double), Double)]
forall a. [a] -> [a]
reverse ([([Double], [Double], Double)]
-> [Double] -> [(([Double], [Double], Double), Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [([Double], [Double], Double)]
hs [Double]
alphas))

-- | Newton/trust-region strategy selected for checked local minimization.
data NewtonMethod = Newton | NewtonCG | TrustRegion Double | Dogleg Double deriving (NewtonMethod -> NewtonMethod -> Bool
(NewtonMethod -> NewtonMethod -> Bool)
-> (NewtonMethod -> NewtonMethod -> Bool) -> Eq NewtonMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: NewtonMethod -> NewtonMethod -> Bool
== :: NewtonMethod -> NewtonMethod -> Bool
$c/= :: NewtonMethod -> NewtonMethod -> Bool
/= :: NewtonMethod -> NewtonMethod -> Bool
Eq, Int -> NewtonMethod -> ShowS
[NewtonMethod] -> ShowS
NewtonMethod -> String
(Int -> NewtonMethod -> ShowS)
-> (NewtonMethod -> String)
-> ([NewtonMethod] -> ShowS)
-> Show NewtonMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> NewtonMethod -> ShowS
showsPrec :: Int -> NewtonMethod -> ShowS
$cshow :: NewtonMethod -> String
show :: NewtonMethod -> String
$cshowList :: [NewtonMethod] -> ShowS
showList :: [NewtonMethod] -> ShowS
Show)

-- | Use an analytic square Hessian with Newton/CG or trust-region acceptance; singular systems and collapsed regions fail.
minimizeNewton ::
    NewtonMethod ->
    OptimizeOptions ->
    Objective ->
    Gradient ->
    Hessian ->
    [Double] ->
    Either MathError OptimizeResult
minimizeNewton :: NewtonMethod
-> OptimizeOptions
-> Objective
-> Gradient
-> Hessian
-> [Double]
-> Either MathError OptimizeResult
minimizeNewton NewtonMethod
method OptimizeOptions
o Objective
f Gradient
g Hessian
hess [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    case NewtonMethod
method of { TrustRegion Double
r -> Double -> Either MathError ()
positive Double
r; Dogleg Double
r -> Double -> Either MathError ()
positive Double
r; NewtonMethod
_ -> () -> Either MathError ()
forall a b. b -> Either a b
Right () }
    Int -> Int -> Double -> [Double] -> Either MathError OptimizeResult
loop Int
0 Int
0 (case NewtonMethod
method of TrustRegion Double
r -> Double
r; Dogleg Double
r -> Double
r; NewtonMethod
_ -> Double
1) [Double]
x0
  where
    positive :: Double -> Either MathError ()
positive Double
r = do
        [Double] -> Either MathError ()
checkedPoint [Double
r]
        if Double
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"positive trust radius required") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    loop :: Int -> Int -> Double -> [Double] -> Either MathError OptimizeResult
loop Int
k Int
n Double
radius [Double]
x = do
        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3)
        y <- Objective -> Objective
checkedValue Objective
f [Double]
x
        gs <- checkedGradient g x
        if norm gs <= optTolerance o
            then Right (OptimizeResult x y (norm gs) k (n + 2) True)
            else
                if k >= optMaxIterations o
                    then Left (NonConvergence (show method) k)
                    else do
                        h <- hess x
                        let dim = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x
                        if length h /= dim || any ((/= dim) . length) h
                            then Left (InvalidArgument "Hessian shape")
                            else checkedPoint (concat h)
                        let cauchy = let hg :: [Double]
hg = [[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
gs; d :: Double
d = [Double] -> [Double] -> Double
dot [Double]
gs [Double]
hg in Double -> [Double] -> [Double]
scale (if Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 then -Double
1 else -[Double] -> [Double] -> Double
dot [Double]
gs [Double]
gs Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d) [Double]
gs
                            cg = [[Double]] -> [Double] -> Double -> Int -> [Double]
truncatedCG [[Double]]
h [Double]
gs (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-12 (OptimizeOptions -> Double
optTolerance OptimizeOptions
o)) Int
dim
                        raw <- case method of NewtonMethod
NewtonCG -> Gradient
forall a b. b -> Either a b
Right [Double]
cg; TrustRegion Double
_ -> Gradient
forall a b. b -> Either a b
Right [Double]
cg; NewtonMethod
_ -> [[Double]] -> Gradient
solveDense [[Double]]
h (Double -> [Double] -> [Double]
scale (-Double
1) [Double]
gs)
                        let descending = if [Double] -> [Double] -> Double
dot [Double]
raw [Double]
gs Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then [Double]
raw else Double -> [Double] -> [Double]
scale (-Double
1) [Double]
gs
                            p = case NewtonMethod
method of
                                TrustRegion Double
_ -> Double -> [Double] -> [Double]
cap Double
radius [Double]
descending
                                Dogleg Double
_ -> Double -> [Double] -> [Double] -> [Double]
dogleg Double
radius [Double]
cauchy [Double]
descending
                                NewtonMethod
_ -> [Double]
descending
                        case method of
                            TrustRegion Double
_ -> Int
-> Int
-> Double
-> [Double]
-> Double
-> [Double]
-> [[Double]]
-> [Double]
-> Either MathError OptimizeResult
accept Int
k (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3) Double
radius [Double]
x Double
y [Double]
gs [[Double]]
h [Double]
p
                            Dogleg Double
_ -> Int
-> Int
-> Double
-> [Double]
-> Double
-> [Double]
-> [[Double]]
-> [Double]
-> Either MathError OptimizeResult
accept Int
k (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3) Double
radius [Double]
x Double
y [Double]
gs [[Double]]
h [Double]
p
                            NewtonMethod
_ -> do
                                (a, _, _, cost) <- LineSearch
-> Int
-> Objective
-> Gradient
-> [Double]
-> [Double]
-> Double
-> [Double]
-> Either MathError (Double, Double, [Double], Int)
lineSearch LineSearch
Armijo (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
80 (OptimizeOptions -> Int
optMaxEvaluations OptimizeOptions
o Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
3)) Objective
f Gradient
g [Double]
x [Double]
p Double
y [Double]
gs
                                loop (k + 1) (n + 3 + cost) radius (add x (scale a p))
    accept :: Int
-> Int
-> Double
-> [Double]
-> Double
-> [Double]
-> [[Double]]
-> [Double]
-> Either MathError OptimizeResult
accept Int
k Int
n Double
radius [Double]
x Double
y [Double]
gs [[Double]]
h [Double]
p = do
        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
        yn <- Objective -> Objective
checkedValue Objective
f ([Double] -> [Double] -> [Double]
add [Double]
x [Double]
p)
        let predicted = -[Double] -> [Double] -> Double
dot [Double]
gs [Double]
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dot [Double]
p ([[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
p)
            rho = if Double
predicted Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 then -Double
1 else (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
yn) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
predicted
            radius' =
                if Double
rho Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.25
                    then Double
radius Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4
                    else if Double
rho Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.75 Bool -> Bool -> Bool
&& [Double] -> Double
norm [Double]
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
radius then Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
1e6 (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
radius) else Double
radius
        if radius' < 1e-16
            then Left (NonConvergence "trust region collapsed" k)
            else loop (k + 1) (n + 1) radius' (if rho > 0.1 then add x p else x)
    cap :: Double -> [Double] -> [Double]
cap Double
r [Double]
p = if [Double] -> Double
norm [Double]
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
r then Double -> [Double] -> [Double]
scale (Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
norm [Double]
p) [Double]
p else [Double]
p
    dogleg :: Double -> [Double] -> [Double] -> [Double]
dogleg Double
r [Double]
pu [Double]
pb
        | [Double] -> Double
norm [Double]
pb Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
r = [Double]
pb
        | [Double] -> Double
norm [Double]
pu Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
r = Double -> [Double] -> [Double]
cap Double
r [Double]
pu
        | Bool
otherwise =
            let d :: [Double]
d = [Double] -> [Double] -> [Double]
sub [Double]
pb [Double]
pu
                a :: Double
a = [Double] -> [Double] -> Double
dot [Double]
d [Double]
d
                b :: Double
b = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dot [Double]
pu [Double]
d
                c :: Double
c = [Double] -> [Double] -> Double
dot [Double]
pu [Double]
pu Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r
                t :: Double
t = (-Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a)
             in [Double] -> [Double] -> [Double]
add [Double]
pu (Double -> [Double] -> [Double]
scale Double
t [Double]
d)

truncatedCG :: [[Double]] -> [Double] -> Double -> Int -> [Double]
truncatedCG :: [[Double]] -> [Double] -> Double -> Int -> [Double]
truncatedCG [[Double]]
h [Double]
g Double
tol Int
limit = Int -> [Double] -> [Double] -> [Double] -> [Double]
go Int
0 (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
g) Double
0) (Double -> [Double] -> [Double]
scale (-Double
1) [Double]
g) (Double -> [Double] -> [Double]
scale (-Double
1) [Double]
g)
  where
    go :: Int -> [Double] -> [Double] -> [Double] -> [Double]
go Int
k [Double]
x [Double]
r [Double]
p
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
limit Bool -> Bool -> Bool
|| [Double] -> Double
norm [Double]
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tol = [Double]
x
        | Double
curvature Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 = if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double -> [Double] -> [Double]
scale (-Double
1) [Double]
g else [Double]
x
        | Bool
otherwise =
            let alpha :: Double
alpha = [Double] -> [Double] -> Double
dot [Double]
r [Double]
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
curvature
                x' :: [Double]
x' = [Double] -> [Double] -> [Double]
add [Double]
x (Double -> [Double] -> [Double]
scale Double
alpha [Double]
p)
                r' :: [Double]
r' = [Double] -> [Double] -> [Double]
sub [Double]
r (Double -> [Double] -> [Double]
scale Double
alpha [Double]
hp)
                beta :: Double
beta = [Double] -> [Double] -> Double
dot [Double]
r' [Double]
r' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> [Double] -> Double
dot [Double]
r [Double]
r
             in Int -> [Double] -> [Double] -> [Double] -> [Double]
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
x' [Double]
r' ([Double] -> [Double] -> [Double]
add [Double]
r' (Double -> [Double] -> [Double]
scale Double
beta [Double]
p))
      where
        hp :: [Double]
hp = [[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
p; curvature :: Double
curvature = [Double] -> [Double] -> Double
dot [Double]
p [Double]
hp

-- | Polak–Ribiere+ nonlinear CG with descent restart and strong Wolfe search.
nonlinearCG ::
    OptimizeOptions -> Objective -> Gradient -> [Double] -> Either MathError OptimizeResult
nonlinearCG :: OptimizeOptions
-> Objective
-> Gradient
-> [Double]
-> Either MathError OptimizeResult
nonlinearCG OptimizeOptions
o Objective
f Gradient
g [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o Int
2
    y <- Objective -> Objective
checkedValue Objective
f [Double]
x0
    gs <- checkedGradient g x0
    loop 0 2 x0 y gs (scale (-1) gs)
  where
    loop :: Int
-> Int
-> [Double]
-> Double
-> [Double]
-> [Double]
-> Either MathError OptimizeResult
loop Int
k Int
n [Double]
x Double
y [Double]
gs [Double]
p
        | [Double] -> Double
norm [Double]
gs Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o = OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right ([Double]
-> Double -> Double -> Int -> Int -> Bool -> OptimizeResult
OptimizeResult [Double]
x Double
y ([Double] -> Double
norm [Double]
gs) Int
k Int
n Bool
True)
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o = MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"nonlinear CG" Int
k)
        | Bool
otherwise = do
            let d :: [Double]
d = if [Double] -> [Double] -> Double
dot [Double]
p [Double]
gs Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then [Double]
p else Double -> [Double] -> [Double]
scale (-Double
1) [Double]
gs
            (a, yn, gn, cost) <- LineSearch
-> Int
-> Objective
-> Gradient
-> [Double]
-> [Double]
-> Double
-> [Double]
-> Either MathError (Double, Double, [Double], Int)
lineSearch LineSearch
StrongWolfe (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
80 (OptimizeOptions -> Int
optMaxEvaluations OptimizeOptions
o Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n)) Objective
f Gradient
g [Double]
x [Double]
d Double
y [Double]
gs
            let beta = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 ([Double] -> [Double] -> Double
dot [Double]
gn ([Double] -> [Double] -> [Double]
sub [Double]
gn [Double]
gs) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> [Double] -> Double
dot [Double]
gs [Double]
gs)
            loop (k + 1) (n + cost) (add x (scale a d)) yn gn (add (scale (-1) gn) (scale beta d))

-- | Caller-supplied residual vector for nonlinear least-squares fitting.
type ResidualFunction = [Double] -> Either MathError [Double]

{- | Caller-supplied residual Jacobian with rows corresponding to residual components and
columns to parameters.
-}
type JacobianFunction = [Double] -> Either MathError [[Double]]

{- | Robust residual penalty selector. Parameterized Huber/Cauchy/SoftL1 choices carry the
selected scale; checked solvers validate it.
-}
data LeastSquaresLoss = SquaredLoss | HuberLoss Double | CauchyLoss Double | SoftL1 Double
    deriving (LeastSquaresLoss -> LeastSquaresLoss -> Bool
(LeastSquaresLoss -> LeastSquaresLoss -> Bool)
-> (LeastSquaresLoss -> LeastSquaresLoss -> Bool)
-> Eq LeastSquaresLoss
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LeastSquaresLoss -> LeastSquaresLoss -> Bool
== :: LeastSquaresLoss -> LeastSquaresLoss -> Bool
$c/= :: LeastSquaresLoss -> LeastSquaresLoss -> Bool
/= :: LeastSquaresLoss -> LeastSquaresLoss -> Bool
Eq, Int -> LeastSquaresLoss -> ShowS
[LeastSquaresLoss] -> ShowS
LeastSquaresLoss -> String
(Int -> LeastSquaresLoss -> ShowS)
-> (LeastSquaresLoss -> String)
-> ([LeastSquaresLoss] -> ShowS)
-> Show LeastSquaresLoss
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LeastSquaresLoss -> ShowS
showsPrec :: Int -> LeastSquaresLoss -> ShowS
$cshow :: LeastSquaresLoss -> String
show :: LeastSquaresLoss -> String
$cshowList :: [LeastSquaresLoss] -> ShowS
showList :: [LeastSquaresLoss] -> ShowS
Show)

{- | lambda=0 selects Gauss–Newton, positive lambda selects LM. Robust losses
use iteratively reweighted least squares with the actual robust objective.
-}
leastSquares ::
    LeastSquaresLoss ->
    Double ->
    OptimizeOptions ->
    ResidualFunction ->
    JacobianFunction ->
    [Double] ->
    Either MathError OptimizeResult
leastSquares :: LeastSquaresLoss
-> Double
-> OptimizeOptions
-> Gradient
-> Hessian
-> [Double]
-> Either MathError OptimizeResult
leastSquares LeastSquaresLoss
loss Double
lambda0 OptimizeOptions
o Gradient
residual Hessian
jac [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    [Double] -> Either MathError ()
checkedPoint [Double
lambda0]
    if Double
lambda0 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative LM damping") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    case LeastSquaresLoss
loss of
        LeastSquaresLoss
SquaredLoss -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()
        HuberLoss Double
c -> Double -> Either MathError ()
positive Double
c
        CauchyLoss Double
c -> Double -> Either MathError ()
positive Double
c
        SoftL1 Double
c -> Double -> Either MathError ()
positive Double
c
    Int -> Int -> Double -> [Double] -> Either MathError OptimizeResult
loop Int
0 Int
0 Double
lambda0 [Double]
x0
  where
    positive :: Double -> Either MathError ()
positive Double
c = do [Double] -> Either MathError ()
checkedPoint [Double
c]; if Double
c Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"robust loss scale") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    rho :: Double -> Double
rho Double
r = case LeastSquaresLoss
loss of
        LeastSquaresLoss
SquaredLoss -> Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        HuberLoss Double
c -> if Double -> Double
forall a. Num a => a -> a
abs Double
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
c then Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 else Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Num a => a -> a
abs Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
        CauchyLoss Double
c -> Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
        SoftL1 Double
c -> Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
    weight :: Double -> Double
weight Double
r = case LeastSquaresLoss
loss of
        LeastSquaresLoss
SquaredLoss -> Double
1
        HuberLoss Double
c -> if Double -> Double
forall a. Num a => a -> a
abs Double
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
c then Double
1 else Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Num a => a -> a
abs Double
r
        CauchyLoss Double
c -> Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
        SoftL1 Double
c -> Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
    loop :: Int -> Int -> Double -> [Double] -> Either MathError OptimizeResult
loop Int
k Int
n Double
lambda [Double]
x = do
        [Double] -> Either MathError ()
checkedPoint [Double]
x
        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
        rs <- Gradient
residual [Double]
x
        checkedPoint rs
        j <- jac x
        if length j /= length rs || any ((/= length x) . length) j
            then Left (InvalidArgument "least squares Jacobian shape")
            else checkedPoint (concat j)
        let ws = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
weight [Double]
rs
            wj = (Double -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> [Double] -> [Double]
scale [Double]
ws [[Double]]
j
            jt = [[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
j
            grad = [[Double]] -> [Double] -> [Double]
matVec [[Double]]
jt ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
ws [Double]
rs)
            normal = [[Double]] -> [[Double]] -> [[Double]]
matMul [[Double]]
jt [[Double]]
wj
            damped = ([Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) [[Double]]
normal (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale Double
lambda) (Int -> [[Double]]
identity ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x)))
            y = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
rho [Double]
rs)
        checkedPoint (y : norm grad : ws ++ grad ++ concat normal ++ concat damped)
        if norm grad <= optTolerance o
            then Right (OptimizeResult x y (norm grad) k (n + 2) True)
            else
                if k >= optMaxIterations o
                    then Left (NonConvergence "least squares" k)
                    else do
                        step <- solveDense damped (scale (-1) grad)
                        checkBudget o (n + 3)
                        checkedPoint (add x step)
                        rn <- residual (add x step)
                        checkedPoint rn
                        if length rn /= length rs then Left (InvalidArgument "residual length changed") else Right ()
                        let yn = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
rho [Double]
rn)
                        checkedPoint [yn]
                        if yn < y
                            then loop (k + 1) (n + 3) (if lambda0 == 0 then 0 else max 1e-15 (lambda / 3)) (add x step)
                            else
                                if lambda0 == 0
                                    then Left (NonConvergence "Gauss Newton step failed to decrease objective" k)
                                    else loop (k + 1) (n + 3) (min 1e20 (lambda * 10)) x

-- | Optimize an axis-initialized simplex with reflection, expansion, contraction, and shrink steps; report simplex spread.
nelderMead :: OptimizeOptions -> Objective -> [Double] -> Either MathError OptimizeResult
nelderMead :: OptimizeOptions
-> Objective -> [Double] -> Either MathError OptimizeResult
nelderMead OptimizeOptions
o Objective
f [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    if [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
x0 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"empty simplex") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let points :: [[Double]]
points =
            [Double]
x0
                [Double] -> [[Double]] -> [[Double]]
forall a. a -> [a] -> [a]
: [ [Double] -> [Double] -> [Double]
add [Double]
x0 [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then OptimizeOptions -> Double
optLearningRate OptimizeOptions
o else Double
0 | Int
j <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                  | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                  ]
    OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
points)
    simplex <- ([Double] -> Either MathError ([Double], Double))
-> [[Double]] -> Either MathError [([Double], Double)]
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 (\[Double]
x -> do y <- Objective -> Objective
checkedValue Objective
f [Double]
x; Right (x, y)) [[Double]]
points
    loop 0 (length simplex) simplex
  where
    loop :: Int
-> Int -> [([Double], Double)] -> Either MathError OptimizeResult
loop Int
k Int
n [([Double], Double)]
simplex = do
        let sorted :: [([Double], Double)]
sorted = (([Double], Double) -> Double)
-> [([Double], Double)] -> [([Double], Double)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn ([Double], Double) -> Double
forall a b. (a, b) -> b
snd [([Double], Double)]
simplex
            ([Double]
best, Double
fb) = [([Double], Double)]
sorted [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
            ([Double]
worst, Double
fw) = [([Double], Double)]
sorted [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0
            second :: Double
second = ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([([Double], Double)]
sorted [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
            rest :: [([Double], Double)]
rest = Int -> [([Double], Double)] -> [([Double], Double)]
forall a. Int -> [a] -> [a]
take ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0) [([Double], Double)]
sorted
            centroid :: [Double]
centroid =
                Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0)) (([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' [Double] -> [Double] -> [Double]
add (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0) Double
0) ((([Double], Double) -> [Double])
-> [([Double], Double)] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst [([Double], Double)]
rest))
            spread :: Double
spread = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((([Double], Double) -> Double) -> [([Double], Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\([Double]
x, Double
y) -> Double -> Double -> Double
forall a. Ord a => a -> a -> a
max ([Double] -> Double
norm ([Double] -> [Double] -> [Double]
sub [Double]
x [Double]
best)) (Double -> Double
forall a. Num a => a -> a
abs (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fb))) [([Double], Double)]
sorted)
        if Double
spread Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o
            then OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right ([Double]
-> Double -> Double -> Int -> Int -> Bool -> OptimizeResult
OptimizeResult [Double]
best Double
fb Double
spread Int
k Int
n Bool
True)
            else
                if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o
                    then MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Nelder Mead" Int
k)
                    else do
                        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0)
                        let reflected :: [Double]
reflected = [Double] -> [Double] -> [Double]
add [Double]
centroid ([Double] -> [Double] -> [Double]
sub [Double]
centroid [Double]
worst)
                        fr <- Objective -> Objective
checkedValue Objective
f [Double]
reflected
                        if fr < fb
                            then do
                                let expanded = [Double] -> [Double] -> [Double]
add [Double]
centroid (Double -> [Double] -> [Double]
scale Double
2 ([Double] -> [Double] -> [Double]
sub [Double]
reflected [Double]
centroid))
                                fe <- checkedValue f expanded
                                loop (k + 1) (n + 2) (rest ++ [if fe < fr then (expanded, fe) else (reflected, fr)])
                            else
                                if fr < second
                                    then loop (k + 1) (n + 1) (rest ++ [(reflected, fr)])
                                    else do
                                        let contracted =
                                                if Double
fr Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
fw
                                                    then [Double] -> [Double] -> [Double]
add [Double]
centroid (Double -> [Double] -> [Double]
scale Double
0.5 ([Double] -> [Double] -> [Double]
sub [Double]
reflected [Double]
centroid))
                                                    else [Double] -> [Double] -> [Double]
add [Double]
centroid (Double -> [Double] -> [Double]
scale Double
0.5 ([Double] -> [Double] -> [Double]
sub [Double]
worst [Double]
centroid))
                                        fc <- checkedValue f contracted
                                        if fc < min fr fw
                                            then loop (k + 1) (n + 2) (rest ++ [(contracted, fc)])
                                            else do
                                                shrunk <-
                                                    mapM
                                                        (\([Double]
x, Double
_) -> let y :: [Double]
y = [Double] -> [Double] -> [Double]
add [Double]
best (Double -> [Double] -> [Double]
scale Double
0.5 ([Double] -> [Double] -> [Double]
sub [Double]
x [Double]
best)) in do fy <- Objective -> Objective
checkedValue Objective
f [Double]
y; Right (y, fy))
                                                        (drop 1 sorted)
                                                loop (k + 1) (n + 2 + length shrunk) ((best, fb) : shrunk)

{- | Powell direction-set method with bounded line minimization on [-radius,radius].
Each outer iteration replaces the direction with greatest decrease by its
net displacement, provided that displacement improves the objective.
-}
powell :: Double -> OptimizeOptions -> Objective -> [Double] -> Either MathError OptimizeResult
powell :: Double
-> OptimizeOptions
-> Objective
-> [Double]
-> Either MathError OptimizeResult
powell Double
radius OptimizeOptions
o Objective
f [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    [Double] -> Either MathError ()
checkedPoint [Double
radius]
    if Double
radius Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
x0
        then MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"Powell radius or dimension")
        else Int
-> Int -> [Double] -> [[Double]] -> Either MathError OptimizeResult
loop Int
0 Int
0 [Double]
x0 (Int -> [[Double]]
identity ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0))
  where
    loop :: Int
-> Int -> [Double] -> [[Double]] -> Either MathError OptimizeResult
loop Int
k Int
n [Double]
x [[Double]]
directions = do
        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
        y <- Objective -> Objective
checkedValue Objective
f [Double]
x
        (xn, yn, cost, changes) <-
            foldM
                ( \([Double]
z, Double
old, Int
count, [Double]
cs) [Double]
d -> do
                    let local :: OptimizeOptions
local = OptimizeOptions
o{optMaxEvaluations = optMaxEvaluations o - count}
                    r <- ScalarMethod
-> OptimizeOptions
-> (Double -> Either MathError Double)
-> Double
-> Double
-> Either MathError OptimizeResult
minimizeScalar ScalarMethod
BrentSearch OptimizeOptions
local (\Double
a -> Objective -> Objective
checkedValue Objective
f ([Double] -> [Double] -> [Double]
add [Double]
z (Double -> [Double] -> [Double]
scale Double
a [Double]
d))) (-Double
radius) Double
radius
                    Right
                        (add z (scale (optPoint r !! 0) d), optValue r, count + optEvaluations r, cs ++ [old - optValue r])
                )
                (x, y, n + 1, [])
                directions
        let displacement = [Double] -> [Double] -> [Double]
sub [Double]
xn [Double]
x
        if norm displacement <= optTolerance o * (1 + norm x) && abs (y - yn) <= optTolerance o
            then Right (OptimizeResult xn yn (norm displacement) k cost True)
            else
                if k >= optMaxIterations o
                    then Left (NonConvergence "Powell" k)
                    else do
                        let index = (Double, Int) -> Double
forall a b. (a, b) -> a
fst ([(Double, Int)] -> (Double, Int)
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [(Double
c, Int
i) | (Int
i, Double
c) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
changes])
                            -- Index is selected below from (decrease,index) to avoid a
                            -- numeric conversion between objective and coordinate types.
                            chosen = (Double, Int) -> Int
forall a b. (a, b) -> b
snd ([(Double, Int)] -> (Double, Int)
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [(Double
c, Int
i) | (Int
i, Double
c) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
changes])
                            dirs =
                                if Double
index Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& [Double] -> Double
norm [Double]
displacement Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
                                    then
                                        [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
chosen then Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
norm [Double]
displacement) [Double]
displacement else [Double]
d
                                        | (Int
i, [Double]
d) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
directions
                                        ]
                                    else [[Double]]
directions
                        loop (k + 1) cost xn dirs