-- SPDX-License-Identifier: MPL-2.0

{- | Seeded bounded CPU global search. Bounds are finite inclusive boxes.
Global optimality cannot be certified by these stochastic algorithms: the
result's convergence flag denotes population/step collapse only. Exhausted
generations return the best point with optConverged=False. Random state is
explicit and deterministic; no system entropy or global generator is used.
-}
module HMath.Optimization.Global (GlobalMethod (..), globalMinimize, basinHopping) where

import Control.Monad (foldM)
import Data.Complex (realPart)
import qualified Data.Foldable as Foldable (foldl')
import Data.List (sortOn)
import Data.Word (Word64)
import HMath.Error
import qualified HMath.Matrix as M
import qualified HMath.Matrix.LAPACK as LA
import HMath.Optimization.Common
import HMath.Optimization.Local (nelderMead)
import HMath.Random (Generator (PCG32), RNG, normal, sampleN, seedRNG, uniform)

{- | Seeded bounded stochastic search strategy. A successful finite run does not certify
global optimality.
-}
data GlobalMethod
    = SimulatedAnnealing
    | DifferentialEvolution
    | GeneticAlgorithm
    | ParticleSwarm
    | CMAES
    deriving (GlobalMethod -> GlobalMethod -> Bool
(GlobalMethod -> GlobalMethod -> Bool)
-> (GlobalMethod -> GlobalMethod -> Bool) -> Eq GlobalMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: GlobalMethod -> GlobalMethod -> Bool
== :: GlobalMethod -> GlobalMethod -> Bool
$c/= :: GlobalMethod -> GlobalMethod -> Bool
/= :: GlobalMethod -> GlobalMethod -> Bool
Eq, Int -> GlobalMethod -> ShowS
[GlobalMethod] -> ShowS
GlobalMethod -> String
(Int -> GlobalMethod -> ShowS)
-> (GlobalMethod -> String)
-> ([GlobalMethod] -> ShowS)
-> Show GlobalMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> GlobalMethod -> ShowS
showsPrec :: Int -> GlobalMethod -> ShowS
$cshow :: GlobalMethod -> String
show :: GlobalMethod -> String
$cshowList :: [GlobalMethod] -> ShowS
showList :: [GlobalMethod] -> ShowS
Show)

validateBounds :: [(Double, Double)] -> Either MathError ()
validateBounds :: [(Double, Double)] -> Either MathError ()
validateBounds [(Double, Double)]
bounds = do
    [Double] -> Either MathError ()
checkedPoint (((Double, Double) -> [Double]) -> [(Double, Double)] -> [Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\(Double
lo, Double
hi) -> [Double
lo, Double
hi, Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lo]) [(Double, Double)]
bounds)
    if [(Double, Double)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Double, Double)]
bounds Bool -> Bool -> Bool
|| ((Double, Double) -> Bool) -> [(Double, Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Double -> Double -> Bool) -> (Double, Double) -> Bool
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(>=)) [(Double, Double)]
bounds
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"global search requires finite nonempty lower < upper bounds")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
clip :: [(Double, Double)] -> [Double] -> [Double]
clip :: [(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds = ((Double, Double) -> Double -> Double)
-> [(Double, Double)] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\(Double
lo, Double
hi) Double
x -> Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
lo (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
hi Double
x)) [(Double, Double)]
bounds
uniformPoint :: [(Double, Double)] -> RNG -> ([Double], RNG)
uniformPoint :: [(Double, Double)] -> RNG -> ([Double], RNG)
uniformPoint [(Double, Double)]
bounds RNG
rng =
    (([Double], RNG) -> (Double, Double) -> ([Double], RNG))
-> ([Double], RNG) -> [(Double, Double)] -> ([Double], RNG)
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]
xs, RNG
s) (Double
lo, Double
hi) -> let (Double
u, RNG
t) = RNG -> (Double, RNG)
uniform RNG
s in ([Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lo) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u], RNG
t))
        ([], RNG
rng)
        [(Double, Double)]
bounds
normals :: Int -> RNG -> Either MathError ([Double], RNG)
normals :: Int -> RNG -> Either MathError ([Double], RNG)
normals Int
n = Int
-> (RNG -> Either MathError (Double, RNG))
-> RNG
-> Either MathError ([Double], RNG)
forall a.
Int
-> (RNG -> Either MathError (a, RNG))
-> RNG
-> Either MathError ([a], RNG)
sampleN Int
n (Double -> Double -> RNG -> Either MathError (Double, RNG)
normal Double
0 Double
1)
population :: Int -> [(Double, Double)] -> RNG -> ([[Double]], RNG)
population :: Int -> [(Double, Double)] -> RNG -> ([[Double]], RNG)
population Int
n [(Double, Double)]
bounds RNG
rng =
    (([[Double]], RNG) -> Int -> ([[Double]], RNG))
-> ([[Double]], RNG) -> [Int] -> ([[Double]], RNG)
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]]
xs, RNG
s) Int
_ -> let ([Double]
x, RNG
t) = [(Double, Double)] -> RNG -> ([Double], RNG)
uniformPoint [(Double, Double)]
bounds RNG
s in ([[Double]]
xs [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]
x], RNG
t))
        ([], RNG
rng)
        [Int
1 .. Int
n]
score :: Objective -> [[Double]] -> Either MathError [([Double], Double)]
score :: Objective -> [[Double]] -> Either MathError [([Double], Double)]
score Objective
f = ([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))
bestResult :: Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult :: Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost Bool
converged [([Double], Double)]
pop =
    let ([Double]
x, Double
y) = (([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)]
pop [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
        spread :: Double
spread = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [[Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs ([Double] -> [Double] -> [Double]
sub [Double]
z [Double]
x)) | ([Double]
z, Double
_) <- [([Double], Double)]
pop])
     in [Double]
-> Double -> Double -> Int -> Int -> Bool -> OptimizeResult
OptimizeResult [Double]
x Double
y Double
spread Int
k Int
cost Bool
converged

-- | Run a reproducible seeded stochastic algorithm in a finite box. size>=4; budget exhaustion returns the best finite incumbent with convergence=False.
globalMinimize ::
    GlobalMethod ->
    Int ->
    Word64 ->
    OptimizeOptions ->
    [(Double, Double)] ->
    Objective ->
    Either MathError OptimizeResult
globalMinimize :: GlobalMethod
-> Int
-> Word64
-> OptimizeOptions
-> [(Double, Double)]
-> Objective
-> Either MathError OptimizeResult
globalMinimize GlobalMethod
method Int
size Word64
seed OptimizeOptions
o [(Double, Double)]
bounds Objective
f = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    [(Double, Double)] -> Either MathError ()
validateBounds [(Double, Double)]
bounds
    if Int
size Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
4
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"global population must contain at least 4 points")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let rng :: RNG
rng = Generator -> Word64 -> RNG
seedRNG Generator
PCG32 Word64
seed
    case GlobalMethod
method of
        GlobalMethod
SimulatedAnnealing -> do
            let ([Double]
x, RNG
s) = [(Double, Double)] -> RNG -> ([Double], RNG)
uniformPoint [(Double, Double)]
bounds RNG
rng
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o Int
1
            y <- Objective -> Objective
checkedValue Objective
f [Double]
x
            anneal 0 1 s x y (x, y)
        GlobalMethod
CMAES ->
            let center :: [Double]
center = ((Double, Double) -> Double) -> [(Double, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\(Double
a, Double
b) -> (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) [(Double, Double)]
bounds; sigma :: Double
sigma = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (((Double, Double) -> Double) -> [(Double, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\(Double
a, Double
b) -> (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
4) [(Double, Double)]
bounds)
             in Int
-> Int
-> RNG
-> [Double]
-> Double
-> [[Double]]
-> [Double]
-> [Double]
-> Maybe ([Double], Double)
-> Either MathError OptimizeResult
cma Int
0 Int
0 RNG
rng [Double]
center Double
sigma (Int -> [[Double]]
identity Int
dim) (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
dim Double
0) (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
dim Double
0) Maybe ([Double], Double)
forall a. Maybe a
Nothing
        GlobalMethod
_ -> do
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o Int
size
            let ([[Double]]
xs, RNG
s) = Int -> [(Double, Double)] -> RNG -> ([[Double]], RNG)
population Int
size [(Double, Double)]
bounds RNG
rng
            pop <- Objective -> [[Double]] -> Either MathError [([Double], Double)]
score Objective
f [[Double]]
xs
            case method of
                GlobalMethod
DifferentialEvolution -> Int
-> Int
-> RNG
-> [([Double], Double)]
-> Either MathError OptimizeResult
de Int
0 Int
size RNG
s [([Double], Double)]
pop
                GlobalMethod
GeneticAlgorithm -> Int
-> Int
-> RNG
-> [([Double], Double)]
-> Either MathError OptimizeResult
ga Int
0 Int
size RNG
s [([Double], Double)]
pop
                GlobalMethod
ParticleSwarm -> Int
-> Int
-> RNG
-> [([Double], Double)]
-> [[Double]]
-> [([Double], Double)]
-> Either MathError OptimizeResult
pso Int
0 Int
size RNG
s [([Double], Double)]
pop (Int -> [Double] -> [[Double]]
forall a. Int -> a -> [a]
replicate Int
size (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
dim Double
0)) [([Double], Double)]
pop
  where
    dim :: Int
dim = [(Double, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Double)]
bounds
    width :: [Double]
width = ((Double, Double) -> Double) -> [(Double, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double -> Double) -> (Double, Double) -> Double
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((Double -> Double -> Double) -> Double -> Double -> Double
forall a b c. (a -> b -> c) -> b -> a -> c
flip (-))) [(Double, Double)]
bounds
    spread :: [([Double], b)] -> Double
spread [([Double], b)]
pop = let ([Double]
best, b
_) = (([Double], b) -> b) -> [([Double], b)] -> [([Double], b)]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn ([Double], b) -> b
forall a b. (a, b) -> b
snd [([Double], b)]
pop [([Double], b)] -> Int -> ([Double], b)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 in [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (([Double], b) -> Double) -> [([Double], b)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Double
norm ([Double] -> Double)
-> (([Double], b) -> [Double]) -> ([Double], b) -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> [Double] -> [Double]
`sub` [Double]
best) ([Double] -> [Double])
-> (([Double], b) -> [Double]) -> ([Double], b) -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double], b) -> [Double]
forall a b. (a, b) -> a
fst) [([Double], b)]
pop)
    finished :: Int -> Int -> [([Double], b)] -> Bool
finished Int
k Int
cost [([Double], b)]
pop = Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o Bool -> Bool -> Bool
|| Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
size Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> OptimizeOptions -> Int
optMaxEvaluations OptimizeOptions
o Bool -> Bool -> Bool
|| [([Double], b)] -> Double
forall {b}. Ord b => [([Double], b)] -> Double
spread [([Double], b)]
pop Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o
    anneal :: Int
-> Int
-> RNG
-> [Double]
-> Double
-> ([Double], Double)
-> Either MathError OptimizeResult
anneal Int
k Int
cost RNG
rng [Double]
x Double
y ([Double], Double)
best
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o Bool -> Bool -> Bool
|| Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> OptimizeOptions -> Int
optMaxEvaluations OptimizeOptions
o =
            OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right
                ( (Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost Bool
False [([Double], Double)
best])
                    { optGradientNorm = 0.1 * sqrt (max 1e-12 (0.95 ^ k)) * maximum width
                    }
                )
        | Bool
otherwise = do
            (zs, s) <- Int -> RNG -> Either MathError ([Double], RNG)
normals Int
dim RNG
rng
            let temperature = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-12 (Double
0.95 Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
k)
                trial = [(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds ([Double] -> [Double] -> [Double]
add [Double]
x ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
w Double
z -> Double
0.1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
temperature Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) [Double]
width [Double]
zs))
                (u, t) = uniform s
            yn <- checkedValue f trial
            let accepted = Double
yn Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
y Bool -> Bool -> Bool
|| Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-300 Double
u) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (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
temperature
                best' = if Double
yn Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
best then ([Double]
trial, Double
yn) else ([Double], Double)
best
            anneal (k + 1) (cost + 1) t (if accepted then trial else x) (if accepted then yn else y) best'
    de :: Int
-> Int
-> RNG
-> [([Double], Double)]
-> Either MathError OptimizeResult
de Int
k Int
cost RNG
rng [([Double], Double)]
pop
        | Int -> Int -> [([Double], Double)] -> Bool
forall {b}. Ord b => Int -> Int -> [([Double], b)] -> Bool
finished Int
k Int
cost [([Double], Double)]
pop = OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right (Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost ([([Double], Double)] -> Double
forall {b}. Ord b => [([Double], b)] -> Double
spread [([Double], Double)]
pop Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o) [([Double], Double)]
pop)
        | Bool
otherwise = do
            (trials, s) <-
                (([[Double]], RNG)
 -> (Int, ([Double], Double)) -> Either MathError ([[Double]], RNG))
-> ([[Double]], RNG)
-> [(Int, ([Double], Double))]
-> Either MathError ([[Double]], RNG)
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM
                    ( \([[Double]]
acc, RNG
state) (Int
i, ([Double]
x, Double
_)) -> do
                        let others :: [Int]
others = [Int
j | Int
j <- [Int
0 .. Int
size Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
                            (Double
u1, RNG
s1) = RNG -> (Double, RNG)
uniform RNG
state
                            p1 :: Int
p1 = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
others Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
u1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
others)))
                            a :: Int
a = [Int]
others [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
p1
                            rest :: [Int]
rest = (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
a) [Int]
others
                            (Double
u2, RNG
s2) = RNG -> (Double, RNG)
uniform RNG
s1
                            b :: Int
b = [Int]
rest [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
rest Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
u2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
rest)))
                            rest2 :: [Int]
rest2 = (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
b) [Int]
rest
                            (Double
u3, RNG
s3) = RNG -> (Double, RNG)
uniform RNG
s2
                            c :: Int
c = [Int]
rest2 [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
rest2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
u3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
rest2)))
                            mutant :: [Double]
mutant = [Double] -> [Double] -> [Double]
add (([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ([([Double], Double)]
pop [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
a)) (Double -> [Double] -> [Double]
scale Double
0.8 ([Double] -> [Double] -> [Double]
sub (([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ([([Double], Double)]
pop [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
b)) (([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ([([Double], Double)]
pop [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
c))))
                            (Double
forcedU, RNG
s4) = RNG -> (Double, RNG)
uniform RNG
s3
                            forced :: Int
forced = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
dim Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
forcedU Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dim))
                            ([Double]
crossed, RNG
s5) =
                                (([Double], RNG) -> (Int, (Double, Double)) -> ([Double], RNG))
-> ([Double], RNG) -> [(Int, (Double, Double))] -> ([Double], RNG)
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]
vs, RNG
st) (Int
j, (Double
orig, Double
new)) -> let (Double
u, RNG
t) = RNG -> (Double, RNG)
uniform RNG
st in ([Double]
vs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
forced Bool -> Bool -> Bool
|| Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.9 then Double
new else Double
orig], RNG
t)
                                    )
                                    ([], RNG
s4)
                                    ([Int] -> [(Double, Double)] -> [(Int, (Double, Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
x [Double]
mutant))
                        ([[Double]], RNG) -> Either MathError ([[Double]], RNG)
forall a b. b -> Either a b
Right ([[Double]]
acc [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds [Double]
crossed], RNG
s5)
                    )
                    ([], RNG
rng)
                    ([Int] -> [([Double], Double)] -> [(Int, ([Double], Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [([Double], Double)]
pop)
            scored <- score f trials
            de (k + 1) (cost + size) s (zipWith (\([Double], Double)
old ([Double], Double)
new -> if ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
new Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
old then ([Double], Double)
new else ([Double], Double)
old) pop scored)
    ga :: Int
-> Int
-> RNG
-> [([Double], Double)]
-> Either MathError OptimizeResult
ga Int
k Int
cost RNG
rng [([Double], Double)]
pop
        | Int -> Int -> [([Double], Double)] -> Bool
forall {b}. Ord b => Int -> Int -> [([Double], b)] -> Bool
finished Int
k Int
cost [([Double], Double)]
pop = OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right (Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost ([([Double], Double)] -> Double
forall {b}. Ord b => [([Double], b)] -> Double
spread [([Double], Double)]
pop Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o) [([Double], Double)]
pop)
        | Bool
otherwise = 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)]
pop; elite :: ([Double], Double)
elite = [([Double], Double)]
sorted [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0; pool :: [([Double], Double)]
pool = Int -> [([Double], Double)] -> [([Double], Double)]
forall a. Int -> [a] -> [a]
take (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
2 (Int
size Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)) [([Double], Double)]
sorted
            (children, s) <-
                (([[Double]], RNG) -> Int -> Either MathError ([[Double]], RNG))
-> ([[Double]], RNG) -> [Int] -> Either MathError ([[Double]], RNG)
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM
                    ( \([[Double]]
acc, RNG
state) Int
_ -> do
                        let (Double
ua, RNG
s1) = RNG -> (Double, RNG)
uniform RNG
state
                            (Double
ub, RNG
s2) = RNG -> (Double, RNG)
uniform RNG
s1
                            a :: [Double]
a = ([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ([([Double], Double)]
pool [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([([Double], Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [([Double], Double)]
pool Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
ua Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([([Double], Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [([Double], Double)]
pool))))
                            b :: [Double]
b = ([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ([([Double], Double)]
pool [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([([Double], Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [([Double], Double)]
pool Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
ub Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([([Double], Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [([Double], Double)]
pool))))
                            (Double
alpha, RNG
s3) = RNG -> (Double, RNG)
uniform RNG
s2
                        (zs, s4) <- Int -> RNG -> Either MathError ([Double], RNG)
normals Int
dim RNG
s3
                        let child =
                                [(Double, Double)] -> [Double] -> [Double]
clip
                                    [(Double, Double)]
bounds
                                    ((Double -> Double -> Double -> Double)
-> [Double] -> [Double] -> [Double] -> [Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\Double
w Double
z Double
base -> Double
base Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.05 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) [Double]
width [Double]
zs ([Double] -> [Double] -> [Double]
add (Double -> [Double] -> [Double]
scale Double
alpha [Double]
a) (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha) [Double]
b)))
                        Right (acc ++ [child], s4)
                    )
                    ([], RNG
rng)
                    [Int
1 .. Int
size Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ys <- score f children
            ga (k + 1) (cost + size - 1) s (elite : ys)
    pso :: Int
-> Int
-> RNG
-> [([Double], Double)]
-> [[Double]]
-> [([Double], Double)]
-> Either MathError OptimizeResult
pso Int
k Int
cost RNG
rng [([Double], Double)]
pop [[Double]]
velocities [([Double], Double)]
personal
        | Int -> Int -> [([Double], Double)] -> Bool
forall {b}. Ord b => Int -> Int -> [([Double], b)] -> Bool
finished Int
k Int
cost [([Double], Double)]
pop = OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right (Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost ([([Double], Double)] -> Double
forall {b}. Ord b => [([Double], b)] -> Double
spread [([Double], Double)]
pop Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o) [([Double], Double)]
personal)
        | Bool
otherwise = do
            let global :: [Double]
global = ([Double], Double) -> [Double]
forall a b. (a, b) -> a
fst ((([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)]
personal [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
                ([([Double], [Double])]
pairs, RNG
s) =
                    (([([Double], [Double])], RNG)
 -> (([Double], Double), ([Double], ([Double], Double)))
 -> ([([Double], [Double])], RNG))
-> ([([Double], [Double])], RNG)
-> [(([Double], Double), ([Double], ([Double], Double)))]
-> ([([Double], [Double])], RNG)
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])]
acc, RNG
state) (([Double]
x, Double
_), ([Double]
vel, ([Double]
best, Double
_))) ->
                            let (Double
u, RNG
s1) = RNG -> (Double, RNG)
uniform RNG
state
                                (Double
v, RNG
s2) = RNG -> (Double, RNG)
uniform RNG
s1
                                velocity :: [Double]
velocity =
                                    [Double] -> [Double] -> [Double]
add (Double -> [Double] -> [Double]
scale Double
0.729 [Double]
vel) ([Double] -> [Double] -> [Double]
add (Double -> [Double] -> [Double]
scale (Double
1.49445 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u) ([Double] -> [Double] -> [Double]
sub [Double]
best [Double]
x)) (Double -> [Double] -> [Double]
scale (Double
1.49445 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v) ([Double] -> [Double] -> [Double]
sub [Double]
global [Double]
x)))
                             in ([([Double], [Double])]
acc [([Double], [Double])]
-> [([Double], [Double])] -> [([Double], [Double])]
forall a. [a] -> [a] -> [a]
++ [([(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds ([Double] -> [Double] -> [Double]
add [Double]
x [Double]
velocity), [Double]
velocity)], RNG
s2)
                        )
                        ([], RNG
rng)
                        ([([Double], Double)]
-> [([Double], ([Double], Double))]
-> [(([Double], Double), ([Double], ([Double], Double)))]
forall a b. [a] -> [b] -> [(a, b)]
zip [([Double], Double)]
pop ([[Double]]
-> [([Double], Double)] -> [([Double], ([Double], Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
velocities [([Double], Double)]
personal))
            ys <- Objective -> [[Double]] -> Either MathError [([Double], Double)]
score Objective
f ((([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])]
pairs)
            let pb = (([Double], Double) -> ([Double], Double) -> ([Double], Double))
-> [([Double], Double)]
-> [([Double], Double)]
-> [([Double], Double)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\([Double], Double)
old ([Double], Double)
new -> if ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
new Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
old then ([Double], Double)
new else ([Double], Double)
old) [([Double], Double)]
personal [([Double], Double)]
ys
            pso (k + 1) (cost + size) s ys (map snd pairs) pb
    -- Rank-mu + rank-one covariance adaptation with cumulative step-size
    -- control and h-sigma gating, following Hansen, arXiv:1604.00772.
    -- Symmetric inverse square roots make the cumulated ps path independent
    -- of the orientation chosen by the eigensolver.
    cma :: Int
-> Int
-> RNG
-> [Double]
-> Double
-> [[Double]]
-> [Double]
-> [Double]
-> Maybe ([Double], Double)
-> Either MathError OptimizeResult
cma Int
k Int
cost RNG
rng [Double]
mean Double
sigma [[Double]]
covariance [Double]
path [Double]
pathSigma Maybe ([Double], Double)
best
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o Bool -> Bool -> Bool
|| Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
size Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> OptimizeOptions -> Int
optMaxEvaluations OptimizeOptions
o = case Maybe ([Double], Double)
best of
            Maybe ([Double], Double)
Nothing -> MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"CMA-ES budget" Int
k)
            Just ([Double], Double)
b ->
                OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right
                    ( (Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost Bool
False [([Double], Double)
b])
                        { optGradientNorm = sigma * sqrt (maximum (0 : [row !! i | (i, row) <- zip [0 ..] covariance]))
                        }
                    )
        | Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [[Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | (Int
i, [Double]
row) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
covariance])) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= OptimizeOptions -> Double
optTolerance OptimizeOptions
o = case Maybe ([Double], Double)
best of
            Maybe ([Double], Double)
Nothing -> MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"CMA-ES initial collapse" Int
k)
            Just ([Double], Double)
b ->
                OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right
                    ( (Int -> Int -> Bool -> [([Double], Double)] -> OptimizeResult
bestResult Int
k Int
cost Bool
True [([Double], Double)
b])
                        { optGradientNorm = sigma * sqrt (maximum (0 : [row !! i | (i, row) <- zip [0 ..] covariance]))
                        }
                    )
        | Bool
otherwise = do
            (l, inverseRoot) <- [[Double]] -> Either MathError ([[Double]], [[Double]])
covarianceRoots [[Double]]
covariance
            (xs, s) <-
                foldM
                    ( \([[Double]]
acc, RNG
state) Int
_ -> do
                        (z, t) <- Int -> RNG -> Either MathError ([Double], RNG)
normals Int
dim RNG
state; Right (acc ++ [clip bounds (add mean (scale sigma (matVec l z)))], t)
                    )
                    ([], rng)
                    [1 .. size]
            scored <- score f xs
            let selected = Int -> [([Double], Double)] -> [([Double], Double)]
forall a. Int -> [a] -> [a]
take (Int
size Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) ((([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)]
scored)
                mu = [([Double], Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [([Double], Double)]
selected
                rawWeights = [Double -> Double
forall a. Floating a => a -> a
log (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
mu 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
forall a. Floating a => a -> a
log (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) | Int
i <- [Int
1 .. Int
mu]]
                weights = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
rawWeights) [Double]
rawWeights
                muEff = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> [Double] -> Double
dot [Double]
weights [Double]
weights
                mean' = ([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 Int
dim Double
0) ((Double -> ([Double], Double) -> [Double])
-> [Double] -> [([Double], Double)] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
w ([Double]
x, Double
_) -> Double -> [Double] -> [Double]
scale Double
w [Double]
x) [Double]
weights [([Double], Double)]
selected)
                y = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sigma) ([Double] -> [Double] -> [Double]
sub [Double]
mean' [Double]
mean)
                dimension = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dim
                cc = (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
muEff Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dimension) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
dimension Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
muEff Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dimension)
                cs = (Double
muEff Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
dimension Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
muEff Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
5)
                damping = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
muEff Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
dimension Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cs
                chi = Double -> Double
forall a. Floating a => a -> a
sqrt Double
dimension Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dimension) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
21 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dimension Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)))
                pathSigma' = [Double] -> [Double] -> [Double]
add (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cs) [Double]
pathSigma) (Double -> [Double] -> [Double]
scale (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
cs Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cs) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
muEff)) ([[Double]] -> [Double] -> [Double]
matVec [[Double]]
inverseRoot [Double]
y))
                hSigma = [Double] -> Double
norm [Double]
pathSigma' 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
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cs) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
chi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1.4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
dimension Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
                path' = [Double] -> [Double] -> [Double]
add (Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cc) [Double]
path) (Double -> [Double] -> [Double]
scale (if Bool
hSigma then Double -> Double
forall a. Floating a => a -> a
sqrt (Double
cc Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cc) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
muEff) else Double
0) [Double]
y)
                c1 = Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dim Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1.3) 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
muEff)
                cmu = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c1) (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
muEff Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
muEff) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dim Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2) 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
muEff))
                rankMu =
                    ([[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])
-> [[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
(+)))
                        (Int -> [Double] -> [[Double]]
forall a. Int -> a -> [a]
replicate Int
dim (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
dim Double
0))
                        [ ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scale Double
w) ([Double] -> [Double] -> [[Double]]
outer [Double]
d [Double]
d)
                        | (Double
w, ([Double]
x, Double
_)) <- [Double] -> [([Double], Double)] -> [(Double, ([Double], Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
weights [([Double], Double)]
selected
                        , let d :: [Double]
d = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sigma) ([Double] -> [Double] -> [Double]
sub [Double]
x [Double]
mean)
                        ]
                covariance' =
                    ([Double] -> [Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]] -> [[Double]]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3
                        ( (Double -> Double -> Double -> Double)
-> [Double] -> [Double] -> [Double] -> [Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3
                            (\Double
a Double
b Double
c -> (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cmu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if Bool
hSigma then Double
0 else Double
c1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cc Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cc))) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cmu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c)
                        )
                        [[Double]]
covariance
                        ([Double] -> [Double] -> [[Double]]
outer [Double]
path' [Double]
path')
                        [[Double]]
rankMu
            let sigma' = Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (Double
cs Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
damping Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double] -> Double
norm [Double]
pathSigma' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
chi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
                generationBest = (([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)]
scored [([Double], Double)] -> Int -> ([Double], Double)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
                best' = case Maybe ([Double], Double)
best of Just ([Double], Double)
old | ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
old Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< ([Double], Double) -> Double
forall a b. (a, b) -> b
snd ([Double], Double)
generationBest -> ([Double], Double)
old; Maybe ([Double], Double)
_ -> ([Double], Double)
generationBest
            checkedPoint (sigma' : concat covariance')
            cma (k + 1) (cost + size) s mean' sigma' covariance' path' pathSigma' (Just best')

covarianceRoots :: [[Double]] -> Either MathError ([[Double]], [[Double]])
covarianceRoots :: [[Double]] -> Either MathError ([[Double]], [[Double]])
covarianceRoots [[Double]]
covariance = do
    let n :: Int
n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
covariance
    matrix <- Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
n Int
n ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
covariance)
    decomposition <- LA.symmetricEigen matrix
    let eigenvalues = (Complex Double -> Double) -> [Complex Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Double
forall a. Complex a -> a
realPart (Eigen -> [Complex Double]
LA.eigenvalues Eigen
decomposition)
        vectors = (Complex Double -> Double) -> [Complex Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Double
forall a. Complex a -> a
realPart (Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
M.toList (Eigen -> Matrix (Complex Double)
LA.eigenvectors Eigen
decomposition))
        basis = [[[Double]
vectors [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    if any (<= 0) eigenvalues then Left (NotPositiveDefinite 0) else Right ()
    let square =
            [[Double]] -> [[Double]] -> [[Double]]
matMul
                [[([[Double]]
basis [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ([Double]
eigenvalues [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
basis)
        inverse =
            [[Double]] -> [[Double]] -> [[Double]]
matMul
                [[([[Double]]
basis [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt ([Double]
eigenvalues [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
basis)
    checkedPoint (concat square ++ concat inverse)
    Right (square, inverse)

{- | Basin hopping uses genuine local Nelder–Mead descents followed by seeded
Metropolis transitions. Bounds are enforced by rejecting out-of-box trials.
-}
basinHopping ::
    Word64 ->
    Double ->
    Int ->
    OptimizeOptions ->
    [(Double, Double)] ->
    Objective ->
    [Double] ->
    Either MathError OptimizeResult
basinHopping :: Word64
-> Double
-> Int
-> OptimizeOptions
-> [(Double, Double)]
-> Objective
-> [Double]
-> Either MathError OptimizeResult
basinHopping Word64
seed Double
temperature Int
hops OptimizeOptions
o [(Double, Double)]
bounds Objective
f [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    [(Double, Double)] -> Either MathError ()
validateBounds [(Double, Double)]
bounds
    [Double] -> Either MathError ()
checkedPoint (Double
temperature Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
x0)
    if Double
temperature Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Int
hops Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [(Double, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Double)]
bounds Bool -> Bool -> Bool
|| [(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds [Double]
x0 [Double] -> [Double] -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double]
x0
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"basin hopping parameters")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let guarded :: Objective
guarded [Double]
x = if [(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds [Double]
x [Double] -> [Double] -> Bool
forall a. Eq a => a -> a -> Bool
== [Double]
x then Objective
f [Double]
x else Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
1e100
    first <- OptimizeOptions
-> Objective -> [Double] -> Either MathError OptimizeResult
nelderMead OptimizeOptions
o Objective
guarded [Double]
x0
    go guarded 0 (optEvaluations first) (seedRNG PCG32 seed) first first
  where
    go :: Objective
-> Int
-> Int
-> RNG
-> OptimizeResult
-> OptimizeResult
-> Either MathError OptimizeResult
go Objective
guarded Int
k Int
cost RNG
rng OptimizeResult
current OptimizeResult
best
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
hops = OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right OptimizeResult
best{optIterations = k, optEvaluations = cost, optConverged = False}
        | Bool
otherwise = do
            (z, s) <- Int -> RNG -> Either MathError ([Double], RNG)
normals ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x0) RNG
rng
            let trial = [(Double, Double)] -> [Double] -> [Double]
clip [(Double, Double)]
bounds ([Double] -> [Double] -> [Double]
add (OptimizeResult -> [Double]
optPoint OptimizeResult
current) (Double -> [Double] -> [Double]
scale (OptimizeOptions -> Double
optLearningRate OptimizeOptions
o) [Double]
z)); (u, t) = uniform s
            r <- nelderMead (o{optMaxEvaluations = optMaxEvaluations o - cost}) guarded trial
            let accept =
                    OptimizeResult -> Double
optValue OptimizeResult
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< OptimizeResult -> Double
optValue OptimizeResult
current Bool -> Bool -> Bool
|| Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-300 Double
u) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (OptimizeResult -> Double
optValue OptimizeResult
current Double -> Double -> Double
forall a. Num a => a -> a -> a
- OptimizeResult -> Double
optValue OptimizeResult
r) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
temperature
                next = if Bool
accept then OptimizeResult
r else OptimizeResult
current
                best' = if OptimizeResult -> Double
optValue OptimizeResult
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< OptimizeResult -> Double
optValue OptimizeResult
best then OptimizeResult
r else OptimizeResult
best
            go guarded (k + 1) (cost + optEvaluations r) t next best'