Commit 31fc73e62686a393720a3b4c7a34ba8816827c2b

  • avatar
  • Ville Tirronen <ville.tirronen @j…u.fi> (Committer)
  • Mon Oct 31 17:40:53 EET 2011
  • avatar
  • Ville Tirronen <ville.tirronen @j…u.fi> (Author)
  • Mon Oct 31 17:40:53 EET 2011
All equations written. Not checked.
cmaes-test.hs
(63 / 26)
  
1{-# LANGUAGE NoMonomorphismRestriction, ViewPatterns, TypeFamilies, RecordWildCards #-}
1{-# LANGUAGE NoMonomorphismRestriction, ViewPatterns, TypeFamilies, RecordWildCards, ParallelListComp, ScopedTypeVariables #-}
22module Main where
33
44import Diagrams.Prelude hiding ((<.>),(<>),sample)
1212import Control.Monad
1313import Data.Monoid
1414import Data.Foldable (foldMap)
15import Data.Ord
1516import Control.Applicative
1617import Foreign.Storable(Storable)
1718import Control.Monad.Primitive
19import Data.List
20import Control.Arrow
1821import qualified Test.QuickCheck as QC
1922import Test.QuickCheck ((==>))
2023import qualified Data.Vector.Storable as V hiding (toList,fromList)
4949
5050
5151-- * CMA-ES algorithm
52data CMAESParameters = CMAES {λ, μ, μeff, c_σ, d_σ, c_c, c_1, c_μ :: Double
53 ,w :: Vector Double}
52data CMAESParameters = CMAES {λ, μ, μeff, c_σ, d_σ, c_c, c_1, c_μ, n :: Double
53 ,w :: [Double]}
5454
5555defaultParams :: Int -> CMAESParameters
5656defaultParams n' = CMAES{..}
6060 μ = fromIntegral (floor μ')
6161 μ' = λ / 2
6262 w = weights λ
63 μeff = 1 / V.sum (V.map (**2) w)
63 μeff = 1 / sum (map (**2) w)
6464 c_σ = (μeff + 2) / (n + μeff + 5) -- 46
6565 d_σ = 1 + 2 * max 0 (sqrt ((μeff-1)/(n+1)) + c_σ) --46
6666 c_c = (4 + μeff/n) / (n + 4 + 2*μeff/n) -- 47
6969 c_μ = min (1-c_1)
7070 (σ_μ*((μeff-2+1/μeff)
7171 /((n+2)**2+σ_μ*μeff/2)))
72{-
73cmaEs cov σ gen = do
74 -- Perform eigendecomposition of cov
75 (d,b) <- eigSH cov
72
73
74cmaEs :: CMAESParameters -> (Vector Double -> Double)
75 -> Matrix Double
76 -> Vector Double
77 -> Vector Double
78 -> Vector Double
79 -> Vector Double
80 -> Int
81 -> GenIO
82 -> IO (Matrix Double, Vector Double, Vector Double, Vector Double, Vector Double, Int)
83cmaEs CMAES{..} f cov p_σ p_c σ m g generator = do
84 -- Perform eigendecomposition of cov. TODO: This doesn't need to be done on every g
85 let (d,b) = eigSH cov
86
7687 -- Generate new population
77 zk <- replicateM λ (normalVector gen)
88 zk <- replicateM (floor λ) (normalVector dim generator)
7889 let yk = map (b <> diag d <>) zk
79 xk = m + σ * yk
90 xk = map ((m +) . (σ *)) yk
91
8092 -- Selection and recombination
81 let selected = take λ . map fst . sortBy (comparing fst) . map (f &&& id) $ yk
82 yw = sum [weight*y | y <- selected]
83 m' = m + σ*yw
84 -- Step size control
85 let sqrtc = b <> (negate d) <> trans b
86 p_σ' = (1-c_σ)*p_σ + sqrt (c_σ*(2-c_σ)*μeff) * sqrtc <> yw
87-}
93 let selected = take (floor λ) . map snd . sortBy (comparing fst) . map (f &&& id) $ yk
94 yw = sum [weight .* y | y <- selected | weight <- w ]
95 m' = m + σ * yw
8896
97 -- Step size control (TODO: Check if updated variables should be used)
98 let sqrtc = b <> (diag (negate d)) <> trans b
99 expectD = sqrt n * (1 - 1/(4*n) + 1/(21*n**2))
100 p_σ' = (1-c_σ) .* p_σ + sqrt (c_σ*(2-c_σ)*μeff) .* (sqrtc <> yw)
101 σ' = σ *. exp ((c_σ/d_σ)*(norm p_σ' / expectD - 1) )
102
103 -- Covariance Matrix Adaptation
104 let h_σ | (norm p_σ)/(sqrt (1-(1-c_σ)**(2*fromIntegral g+1))) < (1.4 + 2/(n+1)) * expectD = 1
105 | otherwise = 0
106 δ = (1 - h_σ) * c_c * (2 - c_c) -- Should be less than one
107
108 p_c' = (1-c_c) .* p_σ' + h_σ .* (sqrt (c_c*(2-c_c)*μeff) .* yw)
109 cov' = (1-c_1-c_μ) .* cov + c_1 .* (p_c `outer` p_c + (δ .* cov))
110 + c_μ .* sum [wi .* yi `outer` yi | wi <- w | yi <- selected]
111
112 return (cov', p_σ', p_c', σ', m', g+1)
113 where
114 norm = pnorm PNorm2
115 dim = V.length p_σ -- TODO: Dimensionality of the problem
116
117
89118-- * Random Numbers
90normalVector n gen = fromList <$> replicateM n gen normal
119normalVector :: Int -> GenIO -> IO (Vector Double)
120normalVector n gen = (fromList <$> replicateM n (normal gen))
91121
92122-- * CMA-ES specific auxiliaries
93123
94124lambda n = 4 + 3 * log n
95125
96prop_weightSum λ = λ > 3 ==> abs (sum (toList $ weights λ) - 1) < 0.00001
97prop_weightPos λ = λ > 3 ==> all (>0) (toList $ weights λ)
126prop_weightSum λ = λ > 3 ==> abs (sum (weights λ) - 1) < 0.00001
127prop_weightPos λ = λ > 3 ==> all (>0) (weights λ)
98128
99weights :: Double -> Vector Double
100weights λ = V.fromList . map (/sum w') $ w'
129weights :: Double -> [Double]
130weights λ = map (/sum w') $ w'
101131 where
102132 μ' = λ/2
103133 μ = floor μ'
135135
136136-- * Visualization
137137
138graph xs = system ((2><2) [1,0 , 0,1]) <§> (lc red . fromVertices . map P . zip [0..] . V.toList $ xs)
138graph xs = system ((2><2) [1,0 , 0,1]) <§> (lc red . fromVertices . map P . zip [0..] $ xs)
139139
140140points color = foldMap (\(pair -> t) -> translate t (circle 0.06 # fc color))
141141
142system :: (V a ~ (Double, Double), PathLike a, Transformable a) => Matrix Double -> a
143system m = mconcat [arr (pair v) | v <- toRows m]
142system :: (V a ~ (examplesDouble, Double), PathLike a, Transformable a) => Matrix Double -> a
143system m = mconcat [arrow (pair v) | v <- toRows m]
144144
145arr v = P (0,0) ~~ P v <§>
145arrow v = P (0,0) ~~ P v <§>
146146 (polygon with { polyType = PolyRegular 3 0.1
147147 , polyOrient = OrientTo (negateV v)
148148 } # translate v )
162162 return $ m + b <> diag d <> (s * fromList xs)
163163
164164-- * Numerical
165
166-- | short hand for scaling with scalar
167(*.) :: Container c e => e -> c e -> c e
168(*.) = flip N.scale
165169
166170empiricalCov :: [Vector Double] -> Matrix Double
167171empiricalCov xs = 1/(lambda-1) * sum [(xi - m) `outer` (xi - m) | xi <- xs]