module Moo.GeneticAlgorithm.Multiobjective.NSGA2 where
import Moo.GeneticAlgorithm.Types
import Moo.GeneticAlgorithm.Multiobjective.Types
import Moo.GeneticAlgorithm.Random
import Moo.GeneticAlgorithm.Utilities (doCrossovers)
import Moo.GeneticAlgorithm.Selection (tournamentSelect)
import Moo.GeneticAlgorithm.Constraints
import Moo.GeneticAlgorithm.Run (makeStoppable)
import Control.Monad (forM_, (<=<), when, liftM)
import Control.Monad.ST (ST)
import Data.Array (array, (!), elems, listArray)
import Data.Array.ST (STArray, runSTArray, newArray, readArray, writeArray, getElems, getBounds)
import Data.Function (on)
import Data.List (sortBy)
import Data.STRef
type DominationCmp a = MultiPhenotype a -> MultiPhenotype a -> Bool
domination :: [ProblemType]
-> DominationCmp a
domination ptypes p q =
let pvs = takeObjectiveValues p
qvs = takeObjectiveValues q
pqs = zip3 ptypes pvs qvs
qps = zip3 ptypes qvs pvs
in (any better1 pqs) && (all (not . better1) qps)
where
better1 :: (ProblemType, Objective, Objective) -> Bool
better1 (Minimizing, pv, qv) = pv < qv
better1 (Maximizing, pv, qv) = pv > qv
constrainedDomination :: (Real b, Real c)
=> [Constraint a b]
-> ([Constraint a b] -> Genome a -> c)
-> [ProblemType]
-> DominationCmp a
constrainedDomination constraints violation ptypes p q =
let pok = isFeasible constraints p
qok = isFeasible constraints q
in case (pok, qok) of
(True, True) -> domination ptypes p q
(False, True) -> False
(True, False) -> True
(False, False) ->
let pviolation = violation constraints (takeGenome p)
qviolation = violation constraints (takeGenome q)
in pviolation < qviolation
data RankedSolution a = RankedSolution {
rs'phenotype :: MultiPhenotype a
, rs'nondominationRank :: Int
, rs'localCrowdingDistnace :: Double
} deriving (Show, Eq)
nondominatedSort :: DominationCmp a -> [MultiPhenotype a] -> [[MultiPhenotype a]]
nondominatedSort dominates = nondominatedSortFast dominates
nondominatedSortFast :: DominationCmp a -> [MultiPhenotype a] -> [[MultiPhenotype a]]
nondominatedSortFast dominates gs =
let n = length gs
garray = listArray (0, n1) gs
fronts = runSTArray $ do
sp <- newArray ((0,0), (n1, (n+2)1)) 0 :: ST s (STArray s (Int,Int) Int)
fronts <- newArray ((0,0), (1,n1)) 0 :: ST s (STArray s (Int,Int) Int)
forM_ (zip gs [0..]) $ \(p, pi) -> do
forM_ (zip gs [0..]) $ \(q, qi) -> do
when ( p `dominates` q ) $
includeInSp sp pi qi
when ( q `dominates` p) $
incrementNp sp pi
np <- readArray sp (pi, 0)
when (np == 0) $
addToFront 0 fronts pi
buildFronts sp fronts 0
frontSizes = takeWhile (>0) . take n $ elems fronts
frontElems = map (\i -> garray ! i) . drop n $ elems fronts
in splitAll frontSizes frontElems
where
includeInSp sp pi qi = do
oldspsize <- readArray sp (pi, 1)
writeArray sp (pi, 2 + oldspsize) qi
writeArray sp (pi, 1) (oldspsize + 1)
incrementNp sp pi = do
oldnp <- readArray sp (pi, 0)
writeArray sp (pi, 0) (oldnp + 1)
frontSize fronts i =
readArray fronts (0, i)
frontStartIndex fronts frontno = do
startref <- newSTRef 0
forM_ [0..(frontno1)] $ \i -> do
oldstart <- readSTRef startref
l <- frontSize fronts i
writeSTRef startref (oldstart + l)
readSTRef startref
addToFront frontno fronts pi = do
start <- frontStartIndex fronts frontno
sz <- frontSize fronts frontno
writeArray fronts (1, start + sz) pi
writeArray fronts (0, frontno) (sz + 1)
frontElems fronts i = do
start <- frontStartIndex fronts i
sz <- frontSize fronts i
felems <- newArray (0, sz1) (1) :: ST s (STArray s Int Int)
forM_ [0..sz1] $ \elix ->
readArray fronts (1, start+elix) >>= writeArray felems elix
getElems felems
dominatedSet sp pi = do
sz <- readArray sp (pi, 1)
delems <- newArray (0, sz1) (1) :: ST s (STArray s Int Int)
forM_ [0..sz1] $ \elix ->
readArray sp (pi, 2+elix) >>= writeArray delems elix
getElems delems
buildFronts sp fronts i = do
maxI <- (snd . snd) `liftM` getBounds fronts
if (i >= maxI || i < 0)
then return fronts
else do
fsz <- frontSize fronts i
if fsz <= 0
then return fronts
else do
felems <- frontElems fronts i
forM_ felems $ \pi -> do
dominated <- dominatedSet sp pi
forM_ dominated $ \qi -> do
nq <- liftM (+ (1::Int)) $ readArray sp (qi, 0)
writeArray sp (qi, 0) nq
when (nq <= 0) $
addToFront (i+1) fronts qi
buildFronts sp fronts (i+1)
splitAll [] _ = []
splitAll _ [] = []
splitAll (sz:szs) els =
let (front, rest) = splitAt sz els
in front : (splitAll szs rest)
crowdingDistances :: [[Objective]] -> [Double]
crowdingDistances [] = []
crowdingDistances pop@(objvals:_) =
let m = length objvals
n = length pop
inf = 1.0/0.0 :: Double
ovTable = array ((0,0), (n1, m1))
[ ((i, objid), (pop !! i) !! objid)
| i <- [0..(n1)], objid <- [0..(m1)] ]
distances = runSTArray $ do
ss <- newArray (0, n1) 0.0
forM_ [0..(m1)] $ \objid -> do
let ixs = sortByObjective objid pop
forM_ (zip3 ixs (drop 1 ixs) (drop 2 ixs)) $ \(iprev, i, inext) -> do
sum_of_si <- readArray ss i
let si = (ovTable ! (inext, objid)) (ovTable ! (iprev, objid))
writeArray ss i (sum_of_si + si)
writeArray ss (head ixs) inf
writeArray ss (last ixs) inf
return ss
in elems distances
where
sortByObjective :: Int -> [[Objective]] -> [Int]
sortByObjective i pop = sortIndicesBy (compare `on` (!! i)) pop
crowdedCompare :: RankedSolution a -> RankedSolution a -> Ordering
crowdedCompare (RankedSolution _ ranki disti) (RankedSolution _ rankj distj) =
case (ranki < rankj, ranki == rankj, disti > distj) of
(True, _, _) -> LT
(_, True, True) -> LT
(_, True, False) -> if disti == distj
then EQ
else GT
_ -> GT
rankAllSolutions :: DominationCmp a -> [MultiPhenotype a] -> [[RankedSolution a]]
rankAllSolutions dominates genomes =
let
fronts = nondominatedSort dominates genomes
frontsDists = map (crowdingDistances . map snd) fronts
ranks = iterate (+1) 1
in map rankedSolutions1 (zip3 fronts ranks frontsDists)
where
rankedSolutions1 :: ([MultiPhenotype a], Int, [Double]) -> [RankedSolution a]
rankedSolutions1 (front, rank, dists) =
zipWith (\g d -> RankedSolution g rank d) front dists
nondominatedRanking
:: forall fn a . ObjectiveFunction fn a
=> DominationCmp a
-> MultiObjectiveProblem fn
-> [Genome a]
-> [(Genome a, Objective)]
nondominatedRanking dominates problems genomes =
let egs = evalAllObjectives problems genomes
fronts = nondominatedSort dominates egs
ranks = concatMap assignRanks (zip fronts (iterate (+1) 1))
in ranks
where
assignRanks :: ([MultiPhenotype a], Int) -> [(Genome a, Objective)]
assignRanks (gs, r) = map (\(eg, rank) -> (fst eg, fromIntegral rank)) $ zip gs (repeat r)
nsga2Ranking
:: forall fn a . ObjectiveFunction fn a
=> DominationCmp a
-> MultiObjectiveProblem fn
-> Int
-> [Genome a]
-> [(MultiPhenotype a, Double)]
nsga2Ranking dominates problems n genomes =
let evaledGenomes = evalAllObjectives problems genomes
fronts = rankAllSolutions dominates evaledGenomes
frontSizes = map length fronts
nFullFronts = length . takeWhile (< n) $ scanl1 (+) frontSizes
partialSize = n (sum (take nFullFronts frontSizes))
(frontsFull, frontsPartial) = splitAt nFullFronts fronts
fromFullFronts = concatMap (map assignRank) frontsFull
fromPartialFront = concatMap (map assignRank
. take partialSize
. sortBy crowdedCompare) $
take 1 frontsPartial
in fromFullFronts ++ fromPartialFront
where
assignRank eg =
let r = fromIntegral $ rs'nondominationRank eg
phenotype = rs'phenotype $ eg
in (phenotype, r)
sortIndicesBy :: (a -> a -> Ordering) -> [a] -> [Int]
sortIndicesBy cmp xs = map snd $ sortBy (cmp `on` fst) (zip xs (iterate (+1) 0))
stepNSGA2
:: forall fn a . ObjectiveFunction fn a
=> MultiObjectiveProblem fn
-> SelectionOp a
-> CrossoverOp a
-> MutationOp a
-> StepGA Rand a
stepNSGA2 problems select crossover mutate stop input = do
let dominates = domination (map fst problems)
case input of
(Left _) ->
stepNSGA2'firstGeneration dominates problems select crossover mutate stop input
(Right _) ->
stepNSGA2'nextGeneration dominates problems select crossover mutate stop input
stepNSGA2bt
:: forall fn a . ObjectiveFunction fn a
=> MultiObjectiveProblem fn
-> CrossoverOp a
-> MutationOp a
-> StepGA Rand a
stepNSGA2bt problems crossover mutate stop popstate =
let n = either length length popstate
select = tournamentSelect Minimizing 2 n
in stepNSGA2 problems select crossover mutate stop popstate
stepConstrainedNSGA2
:: forall fn a b c . (ObjectiveFunction fn a, Real b, Real c)
=> [Constraint a b]
-> ([Constraint a b] -> Genome a -> c)
-> MultiObjectiveProblem fn
-> SelectionOp a
-> CrossoverOp a
-> MutationOp a
-> StepGA Rand a
stepConstrainedNSGA2 constraints violation problems select crossover mutate stop input = do
let dominates = constrainedDomination constraints violation (map fst problems)
case input of
(Left _) ->
stepNSGA2'firstGeneration dominates problems select crossover mutate stop input
(Right _) ->
stepNSGA2'nextGeneration dominates problems select crossover mutate stop input
stepConstrainedNSGA2bt
:: forall fn a b c . (ObjectiveFunction fn a, Real b, Real c)
=> [Constraint a b]
-> ([Constraint a b] -> Genome a -> c)
-> MultiObjectiveProblem fn
-> CrossoverOp a
-> MutationOp a
-> StepGA Rand a
stepConstrainedNSGA2bt constraints violation problems crossover mutate stop popstate =
let n = either length length popstate
tournament = tournamentSelect Minimizing 2 n
in stepConstrainedNSGA2 constraints violation problems tournament crossover mutate stop popstate
stepNSGA2'firstGeneration
:: forall fn a . ObjectiveFunction fn a
=> DominationCmp a
-> MultiObjectiveProblem fn
-> SelectionOp a
-> CrossoverOp a
-> MutationOp a
-> StepGA Rand a
stepNSGA2'firstGeneration dominates problems select crossover mutate = do
let objective = nondominatedRanking dominates problems
makeStoppable objective $ \phenotypes -> do
let popsize = length phenotypes
let genomes = map takeGenome phenotypes
selected <- liftM (map takeGenome) $ (shuffle <=< select) phenotypes
newgenomes <- (mapM mutate) <=< (flip doCrossovers crossover) $ selected
let pool = newgenomes ++ genomes
return $ stepNSGA2'poolSelection dominates problems popsize pool
stepNSGA2'nextGeneration
:: forall fn a . ObjectiveFunction fn a
=> DominationCmp a
-> MultiObjectiveProblem fn
-> SelectionOp a
-> CrossoverOp a
-> MutationOp a
-> StepGA Rand a
stepNSGA2'nextGeneration dominates problems select crossover mutate = do
let dummyObjective = nondominatedRanking dominates problems
makeStoppable dummyObjective $ \rankedgenomes -> do
let popsize = length rankedgenomes
selected <- liftM (map takeGenome) $ select rankedgenomes
newgenomes <- (mapM mutate) <=< flip doCrossovers crossover <=< shuffle $ selected
let pool = (map takeGenome rankedgenomes) ++ newgenomes
return $ stepNSGA2'poolSelection dominates problems popsize pool
stepNSGA2'poolSelection
:: forall fn a . ObjectiveFunction fn a
=> DominationCmp a
-> MultiObjectiveProblem fn
-> Int
-> [Genome a]
-> [Phenotype a]
stepNSGA2'poolSelection dominates problems n pool =
let rankedgenomes = let grs = nsga2Ranking dominates problems n pool
in map (\(mp,r) -> (takeGenome mp, r)) grs
selected = take n rankedgenomes
in selected