{-# LANGUAGE LambdaCase #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}
module Moo.GeneticAlgorithm.Continuous
(
module Moo.GeneticAlgorithm.Types
, getRandomGenomes
, uniformGenomes
, rouletteSelect
, stochasticUniversalSampling
, tournamentSelect
, withPopulationTransform
, withScale
, rankScale
, withFitnessSharing
, distance1, distance2, distanceInf
, bestFirst
, blendCrossover
, unimodalCrossover
, unimodalCrossoverRP
, simulatedBinaryCrossover
, module Moo.GeneticAlgorithm.Crossover
, gaussianMutate
, module Moo.GeneticAlgorithm.Random
, module Moo.GeneticAlgorithm.Run
) where
import Control.Monad (liftM, replicateM)
import Data.List (genericLength, foldl')
import Moo.GeneticAlgorithm.Crossover
import Moo.GeneticAlgorithm.LinAlg
import Moo.GeneticAlgorithm.Random
import Moo.GeneticAlgorithm.Selection
import Moo.GeneticAlgorithm.Types
import Moo.GeneticAlgorithm.Run
import Moo.GeneticAlgorithm.Utilities (getRandomGenomes)
uniformGenomes :: Int -> [(Double,Double)] -> [Genome Double]
uniformGenomes popsize ranges =
let dims = map (uncurry subtract) ranges :: [Double]
ndims = length dims :: Int
vol = product dims
mdim = vol ** (1.0/fromIntegral ndims) :: Double
msamples = (fromIntegral popsize) ** (1.0/fromIntegral ndims) :: Double
ptsPerDim = map (\d -> round $ d*msamples/mdim) dims :: [Int]
ptsInLastDims = product $ drop 1 ptsPerDim :: Int
ptsInFirstDim = popsize `div` ptsInLastDims :: Int
ptsPerDim' = ptsInFirstDim : (drop 1 ptsPerDim) :: [Int]
linspaces = zipWith linspace ranges ptsPerDim' :: [[Double]]
in sproduct [[]] linspaces
where
linspace :: (Double, Double) -> Int -> [Double]
linspace (lo, hi) n = map (\i -> (fromIntegral i)*(hi-lo)/fromIntegral (n-1)) [0..n-1]
sproduct :: [[Double]] -> [[Double]] -> [[Double]]
sproduct gs [] = gs
sproduct gs (l:ls) =
let gs' = [x:g | g<-gs, x<-l]
in sproduct gs' ls
distance1 :: (Num a) => [a] -> [a] -> a
distance1 xs ys = sum . map abs $ zipWith (-) xs ys
distance2 :: (Floating a) => [a] -> [a] -> a
distance2 xs ys = sqrt . sum . map (^(2::Int)) $ zipWith (-) xs ys
distanceInf :: (Real a) => [a] -> [a] -> a
distanceInf xs ys = maximum . map abs $ zipWith (-) xs ys
blendCrossover :: Double
-> CrossoverOp Double
blendCrossover _ [] = return ([], [])
blendCrossover _ [celibate] = return ([],[celibate])
blendCrossover alpha (xs:ys:rest) = do
(xs',ys') <- unzip `liftM` mapM (blx alpha) (zip xs ys)
return ([xs',ys'], rest)
where
blx a (x,y) =
let l = min x y - a*d
u = max x y + a*d
d = abs (x - y)
in do
x' <- getRandomR (l, u)
y' <- getRandomR (l, u)
return (x', y')
unimodalCrossover :: Double
-> Double
-> CrossoverOp Double
unimodalCrossover sigma_xi sigma_eta (x1:x2:x3:rest) = do
let d = x2 `minus` x1
let x_mean = 0.5 `scale` (x1 `plus` x2)
let dist3 =
let v31 = x3 `minus` x1
v21 = x2 `minus` x1
base = norm2 v21
area = sqrt $ (dot v31 v31)*(dot v21 v21) - (dot v21 v31)^(2::Int)
h = area / base
in if isNaN h
then norm2 v31
else h
let n = length x1
(parCorr, orthCorrs) <-
if norm2 d > 1e-6
then do
let exs = drop 1 . mkBasis $ d
getNormals n >>= \case
(xi:etas) -> let
xi' = sigma_xi * xi
parCorr = xi' `scale` d
etas' = map (dist3 * sigma_eta *) etas
orthCorrs = zipWith scale etas' exs
in return (parCorr, orthCorrs)
_ -> error "Parameters too short"
else do
let exs = map (basisVector n) [0..n-1]
etas <- getNormals n
let etas' = map (dist3 * sigma_eta *) etas
let orthCorrs = zipWith scale etas' exs
let zeroCorr = replicate n 0.0
return (zeroCorr, orthCorrs)
let totalCorr = foldr plus parCorr orthCorrs
let child1 = x_mean `minus` totalCorr
let child2 = x_mean `plus` totalCorr
return ([child1, child2], x3:rest)
where
getNormals n = do
ps <- replicateM ((n + 1) `div` 2) getNormal2
return . take n $ concatMap (\(x,y) -> [x,y]) ps
basisVector n i = replicate (n-i-1) 0.0 ++ [1] ++ replicate i 0.0
mkBasis :: [Double] -> [[Double]]
mkBasis dir0 =
let n = length dir0
dims = [0..n-1]
ixs = map (basisVector n) dims
in map normalize . reverse $ foldr build [dir0] ixs
where
build ix exs =
let projs = map (proj ix) exs
rem = foldl' minus ix projs
in if norm2 rem <= 1e-6 * maximum (map norm2 exs)
then exs
else rem : exs
unimodalCrossover _ _ [] = return ([], [])
unimodalCrossover _ _ (x1:x2:[]) = return ([x1,x2], [])
unimodalCrossover _ _ [celibate] = return ([], [celibate])
unimodalCrossoverRP :: CrossoverOp Double
unimodalCrossoverRP [] = return ([], [])
unimodalCrossoverRP parents@(x1:_) =
let n = genericLength x1
sigma_xi = 0.5
sigma_eta = 0.35 / sqrt n
in unimodalCrossover sigma_xi sigma_eta parents
simulatedBinaryCrossover :: Double
-> CrossoverOp Double
simulatedBinaryCrossover n (x1:x2:rest) = do
let cdf beta | beta < 0 = 0.0
| beta <= 1.0 = 0.5*beta**(n+1)
| otherwise = 1.0-0.5/beta**(n+1)
u <- getDouble
let solve eps u = solve' 0.0 (upperB 2.0)
where
upperB b | cdf b < u = upperB (b*2)
| otherwise = b
solve' b1 b2 =
let b = 0.5*(b1+b2)
r = cdf b - u
in if abs r < eps
then b
else
if r >= 0
then solve' b1 b
else solve' b b2
let beta = solve 1e-6 u
let xmean = 0.5 `scale` (x1 `plus` x2)
let deltax = (0.5 * beta) `scale` (x2 `minus` x1)
let c1 = xmean `plus` deltax
let c2 = xmean `minus` deltax
return ([c1,c2], rest)
simulatedBinaryCrossover _ celibates = return ([], celibates)
gaussianMutate :: Double
-> Double
-> MutationOp Double
gaussianMutate p sigma vars = mapM mutate vars
where
mutate = withProbability p $ \v -> do
n <- getNormal
return (v + sigma*n)