packages feed

moo (empty) → 1.0

raw patch · 41 files changed

+4201/−0 lines, 41 filesdep +HUnitdep +arraydep +basesetup-changed

Dependencies added: HUnit, array, base, containers, gray-code, mersenne-random-pure64, monad-mersenne-random, moo, mtl, random, random-shuffle, time

Files

+ LICENSE view
@@ -0,0 +1,32 @@+Copyright (c)2011-2013, Sergey Astanin+Copyright (c)2011, Erlend Hamberg++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Erlend Hamberg, nor the name of Sergey+      Astanin, nor the names of other contributors may be used to+      endorse or promote products derived from this software without+      specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Moo/GeneticAlgorithm.hs view
@@ -0,0 +1,146 @@+{-# OPTIONS_GHC -fno-warn-unused-imports #-}++{- |+Copyright    : 2010-2011 Erlend Hamberg, 2011-2013 Sergey Astanin+License      : BSD3+Stability    : experimental+Portability  : portable++A library for custom genetic algorithms.++@+-----------+Quick Start+-----------+@++Import++  * either "Moo.GeneticAlgorithm.Binary"++  * or "Moo.GeneticAlgorithm.Continuous"++Genetic algorithms are used to find good solutions to optimization+and search problems. They mimic the process of natural evolution+and selection.++A genetic algorithm deals with a /population/ of candidate solutions.+Each candidate solution is represented with a 'Genome'. On every+iteration the best genomes are /selected/ ('SelectionOp'). The next+generation is produced through /crossover/ (recombination of the+parents, 'CrossoverOp') and /mutation/ (a random change in the genome,+'MutationOp') of the selected genomes. This process of selection --+crossover -- mutation is repeated until a good solution appears or all+hope is lost.++Genetic algorithms are often defined in terms of minimizing a cost+function or maximizing fitness. This library refers to observed+performance of a genome as 'Objective', which can be minimized as well+as maximized.+++@+--------------------------------+How to write a genetic algorithm+--------------------------------+@++  1. Provide an encoding and decoding functions to convert from model+     variables to genomes and back. See /How to choose encoding/ below.++  2. Write a custom objective function. Its type should be an instance+     of 'ObjectiveFunction' @a@. Functions of type @Genome a -> Objective@+     are commonly used.++  3. Optionally write custom selection ('SelectionOp'), crossover+     ('CrossoverOp') and mutation ('MutationOp') operators or just use+     some standard operators provided by this library. Operators specific+     to binary or continuous algorithms are provided by+     "Moo.GeneticAlgorithm.Binary" and "Moo.GeneticAlgorithm.Continuous"+     modules respectively.++  4. Use 'nextGeneration' or 'nextSteadyState' to create a single step+     of the algorithm, control the iterative process with 'loop',+     'loopWithLog', or 'loopIO'.++  5. Write a function to generate an initial population; for random+     uniform initialization use 'getRandomGenomes'+     or 'getRandomBinaryGenomes'.++Library functions which need access to random number generator work in+'Rand' monad.  You may use a high-level wrapper 'runGA' (or+'runIO' if you used 'loopIO'), which takes care of creating a new random+number generator and running the entire algorithm.++To solve constrained optimization problems, modify initialization and+selection operators (see "Moo.GeneticAlgorithm.Constraints").++To solve multi-objective optimization problems, use NSGA-II algorithm+(see "Moo.GeneticAlgorithm.Multiobjective").++@+----------------------+How to choose encoding+----------------------+@++ * For problems with discrete search space, binary (or Gray)+   encoding of the bit-string is usually used.+   A bit-string is represented as a list of @Bool@ values (@[Bool]@).+   To build a binary genetic algorithm, import "Moo.GeneticAlgorithm.Binary".++ * For problems with continuous search space, it is possible to use a+   vector of real variables as a genome.+   Such a genome is represented as a list of @Double@ or @Float@ values.+   Special crossover and mutation operators should be used.+   To build a continuous genetic algorithm, import+   "Moo.GeneticAlgorithm.Continuous".+++@+--------+Examples+--------+@++Minimizing Beale's function:++@+import Moo.GeneticAlgorithm.Continuous+++beale :: [Double] -> Double+beale [x, y] = (1.5 - x + x*y)**2 + (2.25 - x + x*y*y)**2 + (2.625 - x + x*y*y*y)**2+++popsize = 101+elitesize = 1+tolerance = 1e-6+++selection = tournamentSelect Minimizing 2 (popsize - elitesize)+crossover = unimodalCrossoverRP+mutation = gaussianMutate 0.25 0.1+step = nextGeneration Minimizing beale selection elitesize crossover mutation+stop = IfObjective (\\values -> (minimum values) < tolerance)+initialize = getRandomGenomes popsize [(-4.5, 4.5), (-4.5, 4.5)]+++main = do+  population <- runGA initialize (loop stop step)+  print (head . bestFirst Minimizing $ population)+@++See @examples/@ folder of the source distribution for more examples.++-}++module Moo.GeneticAlgorithm (+) where++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Run+import Moo.GeneticAlgorithm.Utilities+import Moo.GeneticAlgorithm.Binary+import Moo.GeneticAlgorithm.Continuous
+ Moo/GeneticAlgorithm/Binary.hs view
@@ -0,0 +1,240 @@+{-# LANGUAGE BangPatterns #-}+{-# OPTIONS_GHC -fno-warn-unused-imports #-}+{- |++Binary genetic algorithms. Candidates solutions are represented as bit-strings.++Choose Gray code if sudden changes to the variable value after a point+mutation are undesirable, choose binary code otherwise.  In Gray code+two successive variable values differ in only one bit, it may help to+prevent premature convergence.++To apply binary genetic algorithms to real-valued problems, the real+variable may be discretized ('encodeGrayReal' and+'decodeGrayReal'). Another approach is to use continuous genetic+algorithms, see "Moo.GeneticAlgorithm.Continuous".++To encode more than one variable, just concatenate their codes.+++-}++module Moo.GeneticAlgorithm.Binary (+  -- * Types+    module Moo.GeneticAlgorithm.Types++  -- * Encoding+  , encodeGray+  , decodeGray+  , encodeBinary+  , decodeBinary+  , encodeGrayReal+  , decodeGrayReal+  , bitsNeeded+  , splitEvery++  -- * Initialization+  , getRandomBinaryGenomes++  -- * Selection+  , rouletteSelect+  , stochasticUniversalSampling+  , tournamentSelect+  -- ** Scaling and niching+  , withPopulationTransform+  , withScale+  , rankScale+  , withFitnessSharing+  , hammingDistance+  -- ** Sorting+  , bestFirst+++  -- * Crossover+  , module Moo.GeneticAlgorithm.Crossover++  -- * Mutation+  , pointMutate+  , asymmetricMutate+  , constFrequencyMutate++  -- * Control+  , module Moo.GeneticAlgorithm.Random+  , module Moo.GeneticAlgorithm.Run+) where++import Codec.Binary.Gray.List+import Data.Bits+import Data.List (genericLength)++import Moo.GeneticAlgorithm.Crossover+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Selection+import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Run+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Utilities (getRandomGenomes)++-- | How many bits are needed to represent a range of integer numbers+-- @(from, to)@ (inclusive).+bitsNeeded :: (Integral a, Integral b) => (a, a) -> b+bitsNeeded (from, to) =+    let from' = min from to+        to'= max from to+    in  ceiling . logBase (2::Double) . fromIntegral $ (to' - from' + 1)++-- | Encode an integer number in the range @(from, to)@ (inclusive) as+-- binary sequence of minimal length. Use of Gray code means that a+-- single point mutation leads to incremental change of the encoded+-- value.+encodeGray :: (Bits b, Integral b) => (b, b) -> b -> [Bool]+encodeGray = encodeWithCode gray++-- | Decode a binary sequence using Gray code to an integer in the+-- range @(from, to)@ (inclusive). This is an inverse of 'encodeGray'.+-- Actual value returned may be greater than @to@.+decodeGray :: (Bits b, Integral b) => (b, b) -> [Bool] -> b+decodeGray = decodeWithCode binary++-- | Encode an integer number in the range @(from, to)@ (inclusive)+-- as a binary sequence of minimal length. Use of binary encoding+-- means that a single point mutation may lead to sudden big change+-- of the encoded value.+encodeBinary :: (Bits b, Integral b) => (b, b) -> b -> [Bool]+encodeBinary = encodeWithCode id++-- | Decode a binary sequence to an integer in the range @(from, to)@+-- (inclusive). This is an inverse of 'encodeBinary'.  Actual value+-- returned may be greater than @to@.+decodeBinary :: (Bits b, Integral b) => (b, b) -> [Bool] -> b+decodeBinary = decodeWithCode id++-- | Encode a real number in the range @(from, to)@ (inclusive)+-- with @n@ equally spaced discrete values in binary Gray code.+encodeGrayReal :: (RealFrac a) => (a, a) -> Int -> a -> [Bool]+encodeGrayReal range n = encodeGray (0, n-1) . toDiscreteR range n++-- | Decode a binary sequence using Gray code to a real value in the+-- range @(from, to)@, assuming it was discretized with @n@ equally+-- spaced values (see 'encodeGrayReal').+decodeGrayReal :: (RealFrac a) => (a, a) -> Int -> [Bool] -> a+decodeGrayReal range n = fromDiscreteR range n . decodeGray (0, n-1)++-- | Represent a range @(from, to)@ of real numbers with @n@ equally+-- spaced values.  Use it to discretize a real number @val@.+toDiscreteR :: (RealFrac a)+         => (a, a) -- ^ @(from, to)@, the range to be encoded+         -> Int    -- ^ @n@, how many discrete numbers from the range to consider+         -> a      -- ^ a real number in the range @(from, to)@  to discretize+         -> Int    -- ^ a discrete value (normally in the range @(0, n-1)@)+toDiscreteR range n val =+    let from = uncurry min range+        to = uncurry max range+        dx = (to - from) / (fromIntegral (n - 1))+    in  round $ (val - from) / dx++-- | Take a range @(from, to)@ of real numbers with @n@ equally spaced values.+-- Convert @i@-th value to a real number. This is an inverse of 'toDiscreteR'.+fromDiscreteR :: (RealFrac a)+       => (a, a)  -- ^ @(from, to)@, the encoded range+       -> Int     -- ^ @n@, how many discrete numbers from the range to consider+       -> Int     -- ^ a discrete value in the range @(0, n-1)@+       -> a       -- ^ a real number from the range+fromDiscreteR range n i =+    let from = uncurry min range+        to = uncurry max range+        dx = (to - from) / (fromIntegral (n - 1))+    in  from + (fromIntegral i) * dx++-- | Split a list into pieces of size @n@. This may be useful to split+-- the genome into distinct equally sized “genes” which encode+-- distinct properties of the solution.+splitEvery :: Int -> [a] -> [[a]]+splitEvery _ [] = []+splitEvery n xs = let (nxs,rest) = splitAt n xs in nxs : splitEvery n rest++encodeWithCode :: (Bits b, Integral b) => ([Bool] -> [Bool]) -> (b, b) -> b -> [Bool]+encodeWithCode code (from, to) n =+    let from' = min from to+        to' = max from to+        nbits = bitsNeeded (from', to')+    in  code . take nbits . toList' $ n - from'++decodeWithCode :: (Bits b, Integral b) => ([Bool] -> [Bool]) -> (b, b) -> [Bool] -> b+decodeWithCode decode (from, to) bits =+    let from' = min from to+    in  (from' +) . fromList . decode $ bits+++-- | Generate @n@ random binary genomes of length @len@.+-- Return a list of genomes.+getRandomBinaryGenomes :: Int -- ^ how many genomes to generate+                       -> Int -- ^ genome length+                       -> Rand ([Genome Bool])+getRandomBinaryGenomes n len = getRandomGenomes n (replicate len (False,True))+++-- |Flips a random bit along the length of the genome with probability @p@.+-- With probability @(1 - p)@ the genome remains unaffected.+pointMutate :: Double -> MutationOp Bool+pointMutate p = withProbability p $ \bits -> do+       r <- getRandomR (0, length bits - 1)+       let (before, (bit:after)) = splitAt r bits+       return (before ++ (not bit:after))+++-- |Flip @1@s and @0@s with different probabilities. This may help to control+-- the relative frequencies of @1@s and @0@s in the genome.+asymmetricMutate :: Double   -- ^ probability of a @False@ bit to become @True@+                 -> Double   -- ^ probability of a @True@ bit to become @False@+                 -> MutationOp Bool+asymmetricMutate prob0to1 prob1to0 = mapM flipbit+    where+      flipbit False = withProbability prob0to1 (return . not) False+      flipbit True  = withProbability prob1to0 (return . not) True+++-- Preserving the relative frequencies of ones and zeros:+--+-- ones' = p0*(n-ones) + (1-p1)*ones+-- ones + p0*ones + (p1 - 1)*ones = p0*n+-- p0 + p1 = p0 * n / ones+--+-- zeros' = (1-p0)*zeros + p1*(n-zeros)+-- zeros + (p0 - 1)*zeros + p1*zeros = n*p1+-- p0 + p1 = p1 * n / zeros+--+-- => p0 * zeros = p1 * ones+--+-- Average number of changed bits:+--+-- m = p0*zeros + p1*ones+--+-- => p0 = m / (2*zeros)+--    p1 = m / (2*ones)+--+-- Probability of changing a bit:+--+-- p = m / n+--++-- |Flip @m@ bits on average, keeping the relative frequency of @0@s+-- and @1@s in the genome constant.+constFrequencyMutate :: Real a+                     => a                -- ^ average number of bits to change+                     -> MutationOp Bool+constFrequencyMutate m bits =+    let (ones, zeros) = foldr (\b (o,z) -> if b then (o+1,z) else (o,z+1)) (0,0) bits+        p0to1 = fromRational $ 0.5 * (toRational m) / zeros+        p1to0 = fromRational $ 0.5 * (toRational m) / ones+    in  asymmetricMutate p0to1 p1to0 bits+++-- | Hamming distance between @x@ and @y@ is the number of coordinates+-- for which @x_i@ and @y_i@ are different.+--+-- Reference: Hamming, Richard W. (1950), “Error detecting and error+-- correcting codes”, Bell System Technical Journal 29 (2): 147–160,+-- MR 0035935.+hammingDistance :: (Eq a, Num i) => [a] -> [a] -> i+hammingDistance xs ys = genericLength . filter id $ zipWith (/=) xs ys
+ Moo/GeneticAlgorithm/Constraints.hs view
@@ -0,0 +1,290 @@+module Moo.GeneticAlgorithm.Constraints+    (+      ConstraintFunction+    , Constraint()+    , isFeasible+    -- *** Simple equalities and inequalities+    , (.<.), (.<=.), (.>.), (.>=.), (.==.)+    -- *** Double inequalities+    , LeftHandSideInequality()+    , (.<), (.<=), (<.), (<=.)+    -- ** Constrained initalization+    , getConstrainedGenomes+    , getConstrainedBinaryGenomes+    -- ** Constrained selection+    , withDeathPenalty+    , withFinalDeathPenalty+    , withConstraints+    , numberOfViolations+    , degreeOfViolation+    ) where+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Utilities (getRandomGenomes)+import Moo.GeneticAlgorithm.Selection (withPopulationTransform, bestFirst)+++type ConstraintFunction a b = Genome a -> b+++-- Defining a constraint as a pair of function and its boundary value+-- (vs just a boolean valued function) allows for estimating the+-- degree of constraint violation when necessary.++-- | Define constraints using '.<.', '.<=.', '.>.', '.>=.', and '.==.'+-- operators, with a 'ConstraintFunction' on the left hand side.+--+-- For double inequality constraints use pairs of '.<', '<.' and+-- '.<=', '<=.' respectively, with a 'ConstraintFunction' in the middle.+--+-- Examples:+--+-- @+-- function .>=. lowerBound+-- lowerBound .<= function <=. upperBound+-- @+data (Real b) => Constraint a b+    = LessThan (ConstraintFunction a b) b+    -- ^ strict inequality constraint,+    -- function value is less than the constraint value+    | LessThanOrEqual (ConstraintFunction a b) b+    -- ^ non-strict inequality constraint,+    -- function value is less than or equal to the constraint value+    | Equal (ConstraintFunction a b) b+    -- ^ equality constraint,+    -- function value is equal to the constraint value+    | InInterval (ConstraintFunction a b) (Bool, b) (Bool, b)+    -- ^ double inequality, boolean flags indicate if the+    -- bound is inclusive.+++(.<.) :: (Real b) => ConstraintFunction a b -> b -> Constraint a b+(.<.) = LessThan++(.<=.) :: (Real b) => ConstraintFunction a b -> b -> Constraint a b+(.<=.) = LessThanOrEqual++(.>.) :: (Real b) => ConstraintFunction a b -> b -> Constraint a b+(.>.) f v = LessThan (negate . f) (negate v)++(.>=.) :: (Real b) => ConstraintFunction a b -> b -> Constraint a b+(.>=.) f v = LessThanOrEqual (negate . f) (negate v)++(.==.) :: (Real b) => ConstraintFunction a b -> b -> Constraint a b+(.==.) = Equal+++-- Left hand side of the double inequality defined in the form:+-- @lowerBound .<= function <=. upperBound@.+data (Real b) => LeftHandSideInequality a b+    = LeftHandSideInequality (ConstraintFunction a b) (Bool, b)+    -- ^ boolean flag indicates if the bound is inclusive++(.<=) :: (Real b) => b -> ConstraintFunction a b -> LeftHandSideInequality a b+lval .<= f = LeftHandSideInequality f (True, lval)++(.<) :: (Real b) => b -> ConstraintFunction a b -> LeftHandSideInequality a b+lval .< f  = LeftHandSideInequality f (False, lval)++(<.) :: (Real b) => LeftHandSideInequality a b -> b -> Constraint a b+(LeftHandSideInequality f l) <. rval  = InInterval f l (False, rval)++(<=.) :: (Real b) => LeftHandSideInequality a b -> b -> Constraint a b+(LeftHandSideInequality f l) <=. rval = InInterval f l (True,  rval)++++-- | Returns @True@ if a @genome@ represents a feasible solution+-- with respect to the @constraint@.+satisfiesConstraint :: (Real b)+          => Genome a        -- ^ @genome@+          -> Constraint a b  -- ^ @constraint@+          -> Bool+satisfiesConstraint g (LessThan f v)  = f g < v+satisfiesConstraint g (LessThanOrEqual f v) = f g <= v+satisfiesConstraint g (Equal f v) = f g == v+satisfiesConstraint g (InInterval f (inclusive1,v1) (inclusive2,v2)) =+    let v' = f g+        c1 = if inclusive1 then v1 <= v' else v1 < v'+        c2 = if inclusive2 then v' <= v2 else v' < v2+    in  c1 && c2++++-- | Returns @True@ if a @genome@ represents a feasible solution,+-- i.e. satisfies all @constraints@.+isFeasible :: (GenomeState gt a, Real b)+           => [Constraint a b]  -- ^ constraints+           -> gt                -- ^ genome+           -> Bool+isFeasible constraints genome = all ((takeGenome genome) `satisfiesConstraint`) constraints+++-- | Generate @n@ feasible random genomes with individual genome elements+-- bounded by @ranges@.+getConstrainedGenomes :: (Random a, Ord a, Real b)+    => [Constraint a b]   -- ^ constraints+    -> Int                -- ^ @n@, how many genomes to generate+    -> [(a, a)]           -- ^ ranges for individual genome elements+    -> Rand ([Genome a])  -- ^ random feasible genomes+getConstrainedGenomes constraints n ranges+  | n <= 0            = return []+  | otherwise         = do+  candidates <- getRandomGenomes n ranges+  let feasible = filter (isFeasible constraints) candidates+  let found = length feasible+  more <- getConstrainedGenomes constraints (n - found) ranges+  return $ feasible ++ more+++-- | Generate @n@ feasible random binary genomes.+getConstrainedBinaryGenomes :: (Real b)+    => [Constraint Bool b]  -- ^ constraints+    -> Int                  -- ^ @n@, how many genomes to generate+    -> Int                  -- ^ @L@, genome length+    -> Rand [Genome Bool]   -- ^ random feasible genomes+getConstrainedBinaryGenomes constraints n len =+    getConstrainedGenomes constraints n (replicate len (False,True))+++-- | A simple estimate of the degree of (in)feasibility.+--+-- Count the number of constraint violations. Return @0@ if the solution is feasible.+numberOfViolations :: (Real b)+                   => [Constraint a b]  -- ^ constraints+                   -> Genome a  -- ^ genome+                   -> Int  -- ^ the number of violated constraints+numberOfViolations constraints genome =+    let satisfied = map (genome `satisfiesConstraint`) constraints+    in  length $ filter not satisfied+++-- | An estimate of the degree of (in)feasibility.+--+-- Given @f_j@ is the excess of @j@-th constraint function value,+-- return @sum |f_j|^beta@.  For strict inequality constraints, return+-- @sum (|f_j|^beta + eta)@.  Return @0.0@ if the solution is+-- feasible.+--+degreeOfViolation :: Double  -- ^ beta, single violation exponent+                  -> Double  -- ^ eta, equality penalty in strict inequalities+                  -> [Constraint a Double] -- ^ constrains+                  -> Genome a  -- ^ genome+                  -> Double    -- ^ total degree of violation+degreeOfViolation beta eta constraints genome =+    sum $ map violation constraints+  where+    violation (LessThan f v) =+        let v' = f genome+        in  if v' < v+            then 0.0+            else (abs $ v' - v) ** beta + eta+    violation (LessThanOrEqual f v) =+        let v' = f genome+        in  if v' <= v+            then 0.0+            else (abs $ v' - v) ** beta+    violation (Equal f v) =+        let v' = f genome+        in  if v' == v+            then 0.0+            else (abs $ v' - v) ** beta+    violation (InInterval f (incleft, l) (incright, r)) =+        let v' = f genome+            leftok = if incleft+                     then l <= v'+                     else l < v'+            rightok = if incright+                      then r >= v'+                      else r > v'+        in  case (leftok, rightok) of+            (True, True) -> 0.0+            (False, _)   -> (abs $ l - v') ** beta+                            + (fromIntegral . fromEnum . not $ incleft) * eta+            (_, False)   -> (abs $ v' - r) ** beta+                            + (fromIntegral . fromEnum . not $ incright) * eta+++-- | Modify objective function in such a way that 1) any feasible+-- solution is preferred to any infeasible solution, 2) among two+-- feasible solutions the one having better objective function value+-- is preferred, 3) among two infeasible solution the one having+-- smaller constraint violation is preferred.+--+-- Reference: Deb, K. (2000). An efficient constraint handling method+-- for genetic algorithms. Computer methods in applied mechanics and+-- engineering, 186(2), 311-338.+withConstraints :: (Real b, Real c)+    => [Constraint a b]                      -- ^ constraints+    -> ([Constraint a b] -> Genome a -> c)   -- ^ non-negative degree of violation,+                                             -- see 'numberOfViolations' and 'degreeOfViolation'+    -> ProblemType+    -> SelectionOp a+    -> SelectionOp a+withConstraints constraints violation ptype =+    withPopulationTransform (penalizeInfeasible constraints violation ptype)+++penalizeInfeasible :: (Real b, Real c)+    => [Constraint a b]+    -> ([Constraint a b] -> Genome a -> c)+    -> ProblemType+    -> Population a+    -> Population a+penalizeInfeasible constraints violation ptype phenotypes =+        let worst = takeObjectiveValue . head . worstFirst ptype $ phenotypes+            penalize p = let g = takeGenome p+                             v = fromRational . toRational . violation constraints $ g+                         in  if (v > 0)+                             then (g, worst `worsen` v)+                             else p+        in  map penalize phenotypes+   where+    worstFirst Minimizing = bestFirst Maximizing+    worstFirst Maximizing = bestFirst Minimizing++    worsen x delta = if ptype == Minimizing+                     then x + delta+                     else x - delta+++-- | Kill all infeasible solutions after every step of the genetic algorithm.+--+-- “Death penalty is very popular within the evolution strategies community,+-- but it is limited to problems in which the feasible search space is convex+-- and constitutes a reasonably large portion of the whole search space,” --+-- (Coello 1999).+--+-- Coello, C. A. C., & Carlos, A. (1999). A survey of constraint+-- handling techniques used with evolutionary algorithms.+-- Lania-RI-99-04, Laboratorio Nacional de Informática Avanzada.+withDeathPenalty :: (Monad m, Real b)+                 => [Constraint a b]  -- ^ constraints+                 -> StepGA m a        -- ^ unconstrained step+                 -> StepGA m a        -- ^ constrained step+withDeathPenalty cs step =+    \stop popstate -> do+      stepresult <- step stop popstate+      case stepresult of+        StopGA pop -> return (StopGA (filterFeasible cs pop))+        ContinueGA pop -> return (ContinueGA (filterFeasible cs pop))+++-- | Kill all infeasible solutions once after the last step of the+-- genetic algorithm. See also 'withDeathPenalty'.+withFinalDeathPenalty :: (Monad m, Real b)+                      => [Constraint a b]  -- ^ constriants+                      -> StepGA m a        -- ^ unconstrained step+                      -> StepGA m a        -- ^ constrained step+withFinalDeathPenalty cs step =+    \stop popstate -> do+      result <- step stop popstate+      case result of+        (ContinueGA _) -> return result+        (StopGA pop) -> return (StopGA (filterFeasible cs pop))+++filterFeasible :: (Real b) => [Constraint a b] -> Population a -> Population a+filterFeasible cs = filter (isFeasible cs . takeGenome)
+ Moo/GeneticAlgorithm/Continuous.hs view
@@ -0,0 +1,250 @@+{-# LANGUAGE BangPatterns #-}+{-# OPTIONS_GHC -fno-warn-unused-imports #-}+{- |++Continuous (real-coded) genetic algorithms. Candidate solutions are+represented as lists of real variables.++-}+++module Moo.GeneticAlgorithm.Continuous+  (+  -- * Types+    module Moo.GeneticAlgorithm.Types++  -- * Initialization+  , getRandomGenomes++  -- * Selection+  , rouletteSelect+  , stochasticUniversalSampling+  , tournamentSelect+  -- ** Scaling and niching+  , withPopulationTransform+  , withScale+  , rankScale+  , withFitnessSharing+  , distance1, distance2, distanceInf+  -- ** Sorting+  , bestFirst++  -- * Crossover+  -- ** Neighborhood-based operators+  , blendCrossover+  , unimodalCrossover+  , unimodalCrossoverRP+  , simulatedBinaryCrossover+  , module Moo.GeneticAlgorithm.Crossover++  -- * Mutation+  , gaussianMutate++  -- * Control+  , 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.Random+import Moo.GeneticAlgorithm.Utilities (getRandomGenomes)+++-- | 1-norm distance: @sum |x_i - y-i|@.+distance1 :: (Num a) => [a] -> [a] -> a+distance1 xs ys = sum . map abs $ zipWith (-) xs ys+++-- | 2-norm distance: @(sum (x_i - y_i)^2)^(1/2)@.+distance2 :: (Floating a) => [a] -> [a] -> a+distance2 xs ys = sqrt . sum . map (^(2::Int)) $ zipWith (-) xs ys+++-- | Infinity norm distance: @max |x_i - y_i|@.+distanceInf :: (Real a) => [a] -> [a] -> a+distanceInf xs ys = maximum . map abs $ zipWith (-) xs ys+++-- | Blend crossover (BLX-alpha) for continuous genetic algorithms.  For+-- each component let @x@ and @y@ be its values in the first and the+-- second parent respectively. Choose corresponding component values+-- of the children independently from the uniform distribution in the+-- range (L,U), where @L = min (x,y) - alpha * d@, @U = max+-- (x,y) + alpha * d@, and @d = abs (x - y)@. @alpha@ is usually+-- 0.5. Takahashi in [10.1109/CEC.2001.934452] suggests 0.366.+blendCrossover :: Double -- ^ @alpha@, range expansion parameter+               -> 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')++-- | Unimodal normal distributed crossover (UNDX) for continuous+-- genetic algorithms. Recommended parameters according to [ISBN+-- 978-3-540-43330-9] are @sigma_xi = 0.5@, @sigma_eta =+-- 0.35/sqrt(n)@, where @n@ is the number of variables (dimensionality+-- of the search space). UNDX mixes three parents, producing normally+-- distributed children along the line between first two parents, and using+-- the third parent to build a supplementary orthogonal correction+-- component.+--+-- UNDX preserves the mean of the offspring, and also the+-- covariance matrix of the offspring if @sigma_xi^2 = 0.25@.  By+-- preserving distribution of the offspring, /the UNDX can efficiently+-- search in along the valleys where parents are distributed in+-- functions with strong epistasis among parameters/ (idem).+unimodalCrossover :: Double  -- ^ @sigma_xi@, the standard deviation of+                            -- the mix between two principal parents+                  -> Double  -- ^ @sigma_eta@, the standard deviation+                            -- of the single orthogonal component+                  -> CrossoverOp Double+unimodalCrossover sigma_xi sigma_eta (x1:x2:x3:rest) = do+  let d = x2 `minus` x1  -- vector between parents+  let x_mean = 0.5 `scale` (x1 `plus` x2)  -- parents' average+   -- distance to the 3rd parent in the orthogonal subspace+  let dist3 =+          let v31 = x3 `minus` x1+              v21 = x2 `minus` x1+              base = norm2 v21+              -- twice the triangle area+              area = sqrt $ (dot v31 v31)*(dot v21 v21) - (dot v21 v31)^(2::Int)+              h = area / base+          in  if isNaN h    -- if x1 and x2 coincide+                then norm2 v31+                else h+  let n = length x1+  (parCorr, orthCorrs) <-+      if norm2 d > 1e-6+      then do -- distinct parents+        let exs = drop 1 . mkBasis $ d+        (xi:etas) <- getNormals n+        let xi' = sigma_xi * xi+        let parCorr = xi' `scale` d+        let etas' = map (dist3 * sigma_eta *) etas+        let orthCorrs = zipWith scale etas' exs+        return (parCorr, orthCorrs)+      else do -- identical parents, direction d is undefined+        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+  -- drop only two parents of the three, to keep the number of children the same+  return ([child1, child2], x3:rest)+  where+    -- generate a list of n normally distributed random vars+    getNormals n = do+      ps <- replicateM ((n + 1) `div` 2) getNormal2+      return . take n $ concatMap (\(x,y) -> [x,y]) ps+    -- i-th basis vector in n-dimensional space+    basisVector n i = replicate (n-i-1) 0.0 ++ [1] ++ replicate i 0.0+    -- generate orthonormal bases starting from direction dir0+    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   -- skip this vector, as linear depenent with dir0+                else rem : exs  -- add to the list of orthogonalized vectors+unimodalCrossover _ _ [] = return ([], [])+unimodalCrossover _ _ (x1:x2:[]) = return ([x1,x2], [])  -- FIXME the last two+unimodalCrossover _ _ [celibate]  = return ([], [celibate])++-- | Run 'unimodalCrossover' with default recommended parameters.+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++-- | Simulated binary crossover (SBX) operator for continuous genetic+-- algorithms. SBX preserves the average of the parents and has a+-- spread factor distribution similar to single-point crossover of the+-- binary genetic algorithms. If @n > 0@, then the heighest+-- probability density is assigned to the same distance between+-- children as that of the parents.+--+-- The performance of real-coded genetic algorithm with SBX is similar+-- to that of binary GA with a single-point crossover. For details see+-- Simulated Binary Crossover for Continuous Search Space (1995) Agrawal etal.+simulatedBinaryCrossover :: Double  -- ^ non-negative distribution+                                   -- parameter @n@, usually in the+                                   -- range from 2 to 5; for small+                                   -- values of @n@ children far away+                                   -- from the parents are more likely+                                   -- to be chosen.+                         -> CrossoverOp Double+simulatedBinaryCrossover n (x1:x2:rest) = do+  -- let pdf beta | beta >  1.0 = 0.5*(n+1)/beta**(n+2)+  --              | beta >= 0.0 = 0.5*(n+1)*beta**n+  --              | otherwise   = 0.0   -- beta < 0+  let cdf beta | beta < 0    = 0.0+               | beta <= 1.0 = 0.5*beta**(n+1)+               | otherwise   = 1.0-0.5/beta**(n+1)  -- beta > 1.0+  u <- getDouble  -- uniform random variable in [0,1]+  -- solve cdf(beta) = u with absolute residual less than eps > 0+  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)+++-- |For every variable in the genome with probability @p@ replace its+-- value @v@ with @v + sigma*N(0,1)@, where @N(0,1)@ is a normally+-- distributed random variable with mean equal 0 and variance equal 1.+-- With probability @(1 - p)^n@, where @n@ is the number+-- of variables, the genome remains unaffected.+gaussianMutate :: Double  -- ^ probability @p@+               -> Double  -- ^ @sigma@+               -> MutationOp Double+gaussianMutate p sigma vars = mapM mutate vars+  where+    mutate = withProbability p $ \v -> do+               n <- getNormal+               return (v + sigma*n)
+ Moo/GeneticAlgorithm/Crossover.hs view
@@ -0,0 +1,64 @@+{- |++Common crossover operators for genetic algorithms.++-}++module Moo.GeneticAlgorithm.Crossover+  (+  -- ** Discrete operators+    onePointCrossover+  , twoPointCrossover+  , uniformCrossover+  , noCrossover+  -- ** Application+  , doCrossovers+  , doNCrossovers+) where++import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Utilities++import Control.Monad (liftM)++-- | Crossover two lists in exactly @n@ random points.+nPointCrossover :: Int -> ([a], [a]) -> Rand ([a], [a])+nPointCrossover n (xs,ys)+    | n <= 0 = return (xs,ys)+    | otherwise =+  let len = min (length xs) (length ys)+  in  do+    pos <- getRandomR (0, len-n)+    let (hxs, txs) = splitAt pos xs+    let (hys, tys) = splitAt pos ys+    (rxs, rys) <- nPointCrossover (n-1) (tys, txs) -- FIXME: not tail recursive+    return (hxs ++ rxs, hys ++ rys)++-- |Select a random point in two genomes, and swap them beyond this point.+-- Apply with probability @p@.+onePointCrossover :: Double -> CrossoverOp a+onePointCrossover _ []  = return ([],[])+onePointCrossover _ [celibate] = return ([],[celibate])+onePointCrossover p (g1:g2:rest) = do+  (h1,h2) <- withProbability p (nPointCrossover 1) (g1, g2)+  return ([h1,h2], rest)++-- |Select two random points in two genomes, and swap everything in between.+-- Apply with probability @p@.+twoPointCrossover :: Double -> CrossoverOp a+twoPointCrossover _ []  = return ([], [])+twoPointCrossover _ [celibate] = return ([],[celibate])+twoPointCrossover p (g1:g2:rest) = do+  (h1,h2) <- withProbability p (nPointCrossover 2) (g1,g2)+  return ([h1,h2], rest)++-- |Swap individual bits of two genomes with probability @p@.+uniformCrossover :: Double -> CrossoverOp a+uniformCrossover _ []  = return ([], [])+uniformCrossover _ [celibate] = return ([],[celibate])+uniformCrossover p (g1:g2:rest) = do+  (h1, h2) <- unzip `liftM` mapM swap (zip g1 g2)+  return ([h1,h2], rest)+  where+    swap = withProbability p (\(a,b) -> return (b,a))
+ Moo/GeneticAlgorithm/LinAlg.hs view
@@ -0,0 +1,31 @@+{- |++Ersatz linear algebra.++-}++module Moo.GeneticAlgorithm.LinAlg+  ( minus+  , plus+  , scale+  , dot+  , norm2+  , proj+  , normalize+  ) where++minus :: Num a => [a] -> [a] -> [a]+minus xs ys  = zipWith (-) xs ys+plus :: Num a => [a] -> [a] -> [a]+plus xs ys   = zipWith (+) xs ys+scale :: Num a => a -> [a] -> [a]+scale a xs   = map (a*) xs+dot :: Num a => [a] -> [a] -> a+dot xs ys    = sum $ zipWith (*) xs ys+norm2 :: (Num a, Floating a) => [a] -> a+norm2 xs     = sqrt $ dot xs xs+proj :: (Num a, Fractional a) => [a] -> [a] -> [a]+proj xs dir  = ( dot xs dir / dot dir dir ) `scale` dir+normalize :: (Num a, Floating a, Fractional a) => [a] -> [a]+normalize xs = let a = norm2 xs in (1.0/a) `scale` xs+
+ Moo/GeneticAlgorithm/Multiobjective.hs view
@@ -0,0 +1,18 @@+module Moo.GeneticAlgorithm.Multiobjective+    (+    -- * Types+      SingleObjectiveProblem+    , MultiObjectiveProblem+    , MultiPhenotype+    -- * Evaluation+    , evalAllObjectives+    , takeObjectiveValues+    -- * NSGA-II: A non-dominated sorting genetic algorithm+    , stepNSGA2+    , stepNSGA2bt+    , stepConstrainedNSGA2+    , stepConstrainedNSGA2bt+    ) where++import Moo.GeneticAlgorithm.Multiobjective.Types+import Moo.GeneticAlgorithm.Multiobjective.NSGA2
+ Moo/GeneticAlgorithm/Multiobjective/NSGA2.hs view
@@ -0,0 +1,495 @@+{-# LANGUAGE Rank2Types, ConstraintKinds #-}+{- |++NSGA-II. A Fast Elitist Non-Dominated Sorting Genetic+Algorithm for Multi-Objective Optimization.++Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002). A+fast and elitist multiobjective genetic algorithm:+NSGA-II. Evolutionary Computation, IEEE Transactions on, 6(2),+182-197.++Functions to be used:++  'stepNSGA2', 'stepNSGA2bt',+  'stepConstrainedNSGA2', 'stepConstrainedNSGA2bt'++The other functions are exported for testing only.++-}++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+++-- | Returns @True@ if the first solution dominates the second one in+-- some sense.+type DominationCmp a = MultiPhenotype a -> MultiPhenotype a -> Bool+++-- | A solution @p@ dominates another solution @q@ if at least one 'Objective'+-- values of @p@ is better than the respective value of @q@, and the other+-- are not worse.+domination :: [ProblemType] -- ^ problem types per every objective+           -> 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+++-- | A solution p is said to constrain-dominate a solution q, if any of the+-- following is true: 1) Solution p is feasible and q is not. 2) Solutions+-- p and q are both infeasible but solution p has a smaller overall constraint+-- violation. 3) Solutions p and q are feasible, and solution p dominates solution q.+--+-- Reference: (Deb, 2002).+constrainedDomination :: (Real b, Real c)+                      => [Constraint a b]  -- ^ constraints+                      -> ([Constraint a b] -> Genome a -> c)  -- ^ non-negative degree of violation+                      -> [ProblemType]     -- ^ problem types per every objective+                      -> 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+++-- | Solution and its non-dominated rank and local crowding distance.+data RankedSolution a = RankedSolution {+      rs'phenotype :: MultiPhenotype a+    , rs'nondominationRank :: Int  -- ^ @0@ is the best+    , rs'localCrowdingDistnace :: Double  -- ^ @Infinity@ for less-crowded boundary points+    } deriving (Show, Eq)+++-- | Fast non-dominated sort from (Deb et al. 2002).+-- It is should be O(m N^2), with storage requirements of O(N^2).+nondominatedSort :: DominationCmp a -> [MultiPhenotype a] -> [[MultiPhenotype a]]+nondominatedSort dominates = nondominatedSortFast dominates+++-- | This is a direct translation of the pseudocode from (Deb et al. 2002).+nondominatedSortFast :: DominationCmp a -> [MultiPhenotype a] -> [[MultiPhenotype a]]+nondominatedSortFast dominates gs =+    let n = length gs   -- number of genomes+        garray = listArray (0, n-1) gs+        fronts = runSTArray $ do+                     -- structure of sp array:+                     -- sp [pi][0]    -- n_p, number of genomes dominating pi-th genome+                     -- sp [pi][1]    -- size of S_p, how many genomes pi-th genome dominates+                     -- sp [pi][2..]  -- indices of the genomes dominated by pi-th genome+                     --               -- where pi in [0..n-1]+                     --+                     -- structure of the fronts array:+                     -- fronts [0][i]        -- size of the i-th front+                     -- fronts [1][start..start+fsizes[i]-1] -- indices of the elements of the i-th front+                     --                                      -- where start = sum (take (i-1) fsizes)+                     --+                     -- domination table+                     sp <- newArray ((0,0), (n-1, (n+2)-1)) 0 :: ST s (STArray s (Int,Int) Int)+                     -- at most n fronts with 1 element each+                     fronts <- newArray ((0,0), (1,n-1)) 0 :: ST s (STArray s (Int,Int) Int)+                     forM_ (zip gs [0..]) $ \(p, pi) -> do  -- for each p in P+                       forM_ (zip gs [0..]) $ \(q, qi) -> do  -- for each q in P+                         when ( p `dominates` q ) $+                              -- if p dominates q, include q in S_p+                              includeInSp sp pi qi+                         when ( q `dominates` p) $+                              -- if q dominates p, increment n_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)++    -- size of the i-th front+    frontSize fronts i =+        readArray fronts (0, i)++    frontStartIndex fronts frontno = do+      -- start = sum (take (frontno-1) fsizes)+      startref <- newSTRef 0+      forM_ [0..(frontno-1)] $ \i -> do+          oldstart <- readSTRef startref+          l <- frontSize fronts i+          writeSTRef startref (oldstart + l)+      readSTRef startref++    -- adjust fronts array by updating frontno-th front size and appending+    -- pi to its elements; frontno should be the last front!+    addToFront frontno fronts pi = do+      -- update i-th front size and write an index in the correct position+      start <- frontStartIndex fronts frontno+      sz <- frontSize fronts frontno+      writeArray fronts (1, start + sz) pi+      writeArray fronts (0, frontno) (sz + 1)++    -- elements of the i-th front+    frontElems fronts i = do+      start <- frontStartIndex fronts i+      sz <- frontSize fronts i+      felems <- newArray (0, sz-1) (-1) :: ST s (STArray s Int Int)+      forM_ [0..sz-1] $ \elix ->+          readArray fronts (1, start+elix) >>= writeArray felems elix+      getElems felems++    -- elements which are dominated by the element pi+    dominatedSet sp pi = do+      sz <- readArray sp (pi, 1)+      delems <- newArray (0, sz-1) (-1) :: ST s (STArray s Int Int)+      forM_ [0..sz-1] $ \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) -- all fronts are singletons and the last is already built+         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   -- for each member p in F_i+          dominated <- dominatedSet sp pi+          forM_ dominated $ \qi -> do  -- modify each member from the set S_p+               nq <- liftM (+ (-1::Int)) $ readArray sp (qi, 0)  -- decrement n_q by one+               writeArray sp (qi, 0) nq+               when (nq <= 0) $  -- if n_q is zero, q is a member of the next front+                    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)+++-- | Crowding distance of a point @p@, as defined by Deb et+-- al. (2002), is an estimate (the sum of dimensions in their+-- pseudocode) of the largest cuboid enclosing the point without+-- including any other point in the population.+crowdingDistances :: [[Objective]] -> [Double]+crowdingDistances [] = []+crowdingDistances pop@(objvals:_) =+    let m = length objvals  -- number of objectives+        n = length pop      -- number of genomes+        inf = 1.0/0.0 :: Double+        -- (genome-idx, objective-idx) -> objective value+        ovTable = array ((0,0), (n-1, m-1))+                  [ ((i, objid), (pop !! i) !! objid)+                  | i <- [0..(n-1)], objid <- [0..(m-1)] ]+        -- calculate crowding distances+        distances = runSTArray $ do+          ss <- newArray (0, n-1) 0.0  -- initialize distances+          forM_ [0..(m-1)] $ \objid -> do    -- for every objective+            let ixs = sortByObjective objid pop+              -- for all inner points+            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   -- boundary points have infinite cuboids+            writeArray ss (last ixs) inf+          return ss+    in elems distances+  where+    sortByObjective :: Int -> [[Objective]] -> [Int]+    sortByObjective i pop = sortIndicesBy (compare `on` (!! i)) pop++-- | Given there is non-domination rank @rank_i@, and local crowding+-- distance @distance_i@ assigned to every individual @i@, the partial+-- order between individuals @i@ and @q@ is defined by relation+--+-- @i ~ j@ if @rank_i < rank_j@ or (@rank_i = rank_j@ and @distance_i@+-- @>@ @distance_j@).+--+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+++-- | Assign non-domination rank and crowding distances to all solutions.+-- Return a list of non-domination fronts.+rankAllSolutions :: DominationCmp a -> [MultiPhenotype a] -> [[RankedSolution a]]+rankAllSolutions dominates genomes =+    let -- non-dominated fronts+        fronts = nondominatedSort dominates genomes+        -- for every non-dominated front+        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+++-- | To every genome in the population, assign a single objective+-- value according to its non-domination rank. This ranking is+-- supposed to be used once in the beginning of the NSGA-II algorithm.+--+-- Note: 'nondominatedRanking' reorders the genomes.+nondominatedRanking+    :: forall fn a . ObjectiveFunction fn a+    => DominationCmp a+    -> MultiObjectiveProblem fn     -- ^ list of @problems@+    -> [Genome a]                   -- ^ a population of raw @genomes@+    -> [(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)+++-- | To every genome in the population, assign a single objective value+-- equal to its non-domination rank, and sort genomes by the decreasing+-- local crowding distance within every rank+-- (i.e. sort the population with NSGA-II crowded comparision+-- operator)+nsga2Ranking+    :: forall fn a . ObjectiveFunction fn a+    => DominationCmp a+    -> MultiObjectiveProblem fn    -- ^ a list of @objective@ functions+    -> Int                          -- ^ @n@, number of top-ranked genomes to select+    -> [Genome a]                   -- ^ a population of raw @genomes@+    -> [(MultiPhenotype a, Double)] -- ^ selected genomes with their non-domination ranks+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))++-- | A single step of the NSGA-II algorithm (Non-Dominated Sorting+-- Genetic Algorithm for Multi-Objective Optimization).+--+-- The next population is selected from a common pool of parents and+-- their children minimizing the non-domination rank and maximizing+-- the crowding distance within the same rank.+-- The first generation of children is produced without taking+-- crowding into account.+-- Every solution is assigned a single objective value which is its+-- sequence number after sorting with the crowded comparison operator.+-- The smaller value corresponds to solutions which are not worse+-- the one with the bigger value. Use 'evalAllObjectives' to restore+-- individual objective values.+--+-- Reference:+-- Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002). A+-- fast and elitist multiobjective genetic algorithm:+-- NSGA-II. Evolutionary Computation, IEEE Transactions on, 6(2),+-- 182-197.+--+-- Deb et al. used a binary tournament selection, base on crowded+-- comparison operator. To achieve the same effect, use+-- 'stepNSGA2bt' (or 'stepNSGA2' with 'tournamentSelect'+-- @Minimizing 2 n@, where @n@ is the size of the population).+--+stepNSGA2+    :: forall fn a . ObjectiveFunction fn a+    => MultiObjectiveProblem fn    -- ^ a list of @objective@ functions+    -> 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 _) ->  -- raw genomes => it's the first generation+        stepNSGA2'firstGeneration dominates problems select crossover mutate stop input+    (Right _) ->  -- ranked genomes => it's the second or later generation+        stepNSGA2'nextGeneration dominates problems select crossover mutate stop input+++-- | A single step of NSGA-II algorithm with binary tournament selection.+-- See also 'stepNSGA2'.+stepNSGA2bt+    :: forall fn a . ObjectiveFunction fn a+    => MultiObjectiveProblem fn    -- ^ a list of @objective@ functions+    -> 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+++-- | A single step of the constrained NSGA-II algorithm, which uses a+-- constraint-domination rule:+--+-- “A solution @i@ is said to constrain-dominate a solution @j@, if any of the+-- following is true: 1) Solution @i@ is feasible and @j@ is not. 2) Solutions+-- @i@ and @j@ are both infeasible but solution @i@ has a smaller overall constraint+-- violation. 3) Solutions @i@ and @j@ are feasible, and solution @i@ dominates solution @j@.”+--+-- Reference: (Deb, 2002).+--+stepConstrainedNSGA2+    :: forall fn a b c . (ObjectiveFunction fn a, Real b, Real c)+    => [Constraint a b]                     -- ^ constraints+    -> ([Constraint a b] -> Genome a -> c)  -- ^ non-negative degree of violation+    -> MultiObjectiveProblem fn             -- ^ a list of @objective@ functions+    -> 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+++-- | A single step of the constrained NSGA-II algorithm with binary tournament+-- selection. See also 'stepConstrainedNSGA2'.+stepConstrainedNSGA2bt+    :: forall fn a b c . (ObjectiveFunction fn a, Real b, Real c)+    => [Constraint a b]                     -- ^ constraints+    -> ([Constraint a b] -> Genome a -> c)  -- ^ non-negative degree of violation+    -> MultiObjectiveProblem fn             -- ^ a list of @objective@ functions+    -> 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    -- ^ a list of @objective@ functions+    -> 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+++-- | Use normal selection, crossover, mutation to produce new+-- children.  Select from a common pool of parents and children the+-- best according to the least non-domination rank and crowding.+stepNSGA2'nextGeneration+     :: forall fn a . ObjectiveFunction fn a+     => DominationCmp a+     -> MultiObjectiveProblem fn   -- ^ a list of objective functions+     -> SelectionOp a+     -> CrossoverOp a+     -> MutationOp a+     -> StepGA Rand a+stepNSGA2'nextGeneration dominates problems select crossover mutate = do+  -- nextGeneration is never called with raw genomes,+  -- => dummyObjective is never evaluated;+  -- nondominatedRanking is required to type-check+  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+++-- | Take a pool of phenotypes of size 2N, ordered by the crowded+-- comparison operator, and select N best.+stepNSGA2'poolSelection+    :: forall fn a . ObjectiveFunction fn a+    => DominationCmp a+    -> MultiObjectiveProblem fn    -- ^ a list of @objective@ functions+    -> Int                         -- ^ @n@, the number of solutions to select+    -> [Genome a]                  -- ^ @pool@ of genomes to select from+    -> [Phenotype a]               -- ^ @n@ best phenotypes+stepNSGA2'poolSelection dominates problems n pool =+    -- nsga2Ranking returns genomes properly sorted already+    let rankedgenomes = let grs = nsga2Ranking dominates problems n pool+                        in  map (\(mp,r) -> (takeGenome mp, r)) grs+        selected = take n rankedgenomes  -- :: [Phenotype a]+    in  selected
+ Moo/GeneticAlgorithm/Multiobjective/Types.hs view
@@ -0,0 +1,45 @@+{-# LANGUAGE MultiParamTypeClasses, Rank2Types, GADTs, FlexibleInstances #-}++module Moo.GeneticAlgorithm.Multiobjective.Types+    ( SingleObjectiveProblem+    , MultiObjectiveProblem+    , MultiPhenotype+    , evalAllObjectives+    , takeObjectiveValues+    ) where+++import Moo.GeneticAlgorithm.Types+++import Data.List (transpose)+++type SingleObjectiveProblem fn = ( ProblemType , fn )+type MultiObjectiveProblem fn = [SingleObjectiveProblem fn]+++-- | An individual with all objective functions evaluated.+type MultiPhenotype a = (Genome a, [Objective])+++instance a1 ~ a2 => GenomeState (MultiPhenotype a1) a2 where+    takeGenome = fst+++takeObjectiveValues :: MultiPhenotype a -> [Objective]+takeObjectiveValues = snd+++-- | Calculate multiple objective per every genome in the population.+evalAllObjectives+    :: forall fn gt a . (ObjectiveFunction fn a, GenomeState gt a)+    => MultiObjectiveProblem fn    -- ^ a list of @problems@+    -> [gt]                        -- ^ a population of @genomes@+    -> [MultiPhenotype a]+evalAllObjectives problems genomes =+    let rawgenomes = map takeGenome genomes+        pops_per_objective = map (\(_, f) -> evalObjective f rawgenomes) problems+        ovs_per_objective = map (map takeObjectiveValue) pops_per_objective+        ovs_per_genome = transpose ovs_per_objective+    in  zip rawgenomes ovs_per_genome
+ Moo/GeneticAlgorithm/Niching.hs view
@@ -0,0 +1,55 @@+module Moo.GeneticAlgorithm.Niching+    ( fitnessSharing+    ) where+++import Moo.GeneticAlgorithm.Types+++-- | A popular niching method proposed by D. Goldberg and+-- J. Richardson in 1987. The shared fitness of the individual is inversely+-- protoptional to its niche count.+-- The method expects the objective function to be non-negative.+--+-- An extension for minimization problems is implemented by+-- making the fitnes proportional to its niche count (rather than+-- inversely proportional).+--+-- Reference: Chen, J. H., Goldberg, D. E., Ho, S. Y., & Sastry,+-- K. (2002, July). Fitness inheritance in multiobjective+-- optimization. In Proceedings of the Genetic and Evolutionary+-- Computation Conference (pp. 319-326). Morgan Kaufmann Publishers+-- Inc..+fitnessSharing ::+    (Phenotype a -> Phenotype a -> Double)  -- ^ distance function+    -> Double                        -- ^ niche radius+    -> Double                        -- ^ niche compression exponent @alpha@ (usually 1)+    -> ProblemType                   -- ^ type of the optimization problem+    -> Population a+    -> Population a+fitnessSharing dist r alpha Maximizing phenotypes =+    let ms = map (nicheCount dist r alpha phenotypes) phenotypes+    in  zipWith (\(genome, value) m -> (genome, value/m)) phenotypes ms+fitnessSharing dist r alpha Minimizing phenotypes =+    let ms = map (nicheCount dist r alpha phenotypes) phenotypes+    in  zipWith (\(genome, value) m -> (genome, value*m)) phenotypes ms+++type DistanceFunction a = Phenotype a -> Phenotype a -> Double+++nicheCount :: DistanceFunction a+           -> Double -> Double+           -> Population a -> Phenotype a -> Double+nicheCount dist r alpha population phenotype =+    sum $ map (sharing dist r alpha phenotype) population+++sharing :: DistanceFunction a+        -> Double -> Double+        -> DistanceFunction a+sharing dist r alpha pi pj =+    let dij = dist pi pj+    in  if dij < r+        then 1.0 - (dij/r)**alpha+        else 0.0
+ Moo/GeneticAlgorithm/Random.hs view
@@ -0,0 +1,111 @@+{-# LANGUAGE BangPatterns #-}++{- | Some extra facilities to work with 'Rand' monad and 'PureMT'+     random number generator.+-}++module Moo.GeneticAlgorithm.Random+    (+    -- * Random numbers from given range+      getRandomR+    , getRandom+    -- * Probability distributions+    , getNormal2+    , getNormal+    -- * Random samples and shuffles+    , randomSample+    , shuffle+    -- * Building blocks+    , withProbability+    -- * Re-exports from random number generator packages+    , getBool, getInt, getWord, getInt64, getWord64, getDouble+    , runRandom, evalRandom, newPureMT+    , Rand, Random, PureMT+    ) where++import Control.Monad (liftM)+import Control.Monad.Mersenne.Random+import Data.Complex (Complex (..))+import System.Random (RandomGen, Random(..))+import System.Random.Mersenne.Pure64+import qualified System.Random.Shuffle as S++-- | Yield a new randomly selected value of type @a@ in the range @(lo, hi)@.+-- See 'System.Random.randomR' for details.+getRandomR :: Random a => (a, a) -> Rand a+getRandomR range = Rand $ \s -> let (r, s') = randomR range s in R r s'++-- | Yield a new randomly selected value of type @a@.+-- See 'System.Random.random' for details.+getRandom :: Random a => Rand a+getRandom = Rand $ \g -> let (r, g') = random g in R r g'++-- | Yield two randomly selected values which follow standard+-- normal distribution.+getNormal2 :: Rand (Double, Double)+getNormal2 = do+  -- Box-Muller method+  u <- getDouble+  v <- getDouble+  let (c :+ s) = exp (0 :+ (2*pi*v))+  let r = sqrt $ (-2) * log u+  return (r*c, r*s)++-- | Yield one randomly selected value from standard normal distribution.+getNormal :: Rand Double+getNormal = fst `liftM` getNormal2++-- | Take at most n random elements from the list. Preserve order.+randomSample :: Int -> [a] -> Rand [a]+randomSample n xs =+  Rand $ \g -> case select g n (length xs) xs [] of (xs', g') -> R xs' g'+  where+    select rng _ _ [] acc = (reverse acc, rng)+    select rng n m xs acc+        | n <= 0     = (reverse acc, rng)+        | otherwise  =+            let (k, rng') = randomR (0, m - n) rng+                (x:rest) = drop k xs+            in  select rng' (n-1) (m-k-1) rest (x:acc)+++-- | Randomly reorder the list.+shuffle :: [a] -> Rand [a]+shuffle xs = Rand $ \g ->+             let (xs', g') = randomShuffle xs (length xs) g in  R xs' g'++-- | Given a sequence (e1,...en) to shuffle, its length, and a random+-- generator, compute the corresponding permutation of the input+-- sequence, return the permutation and the new state of the+-- random generator.+randomShuffle :: RandomGen gen => [a] -> Int -> gen -> ([a], gen)+randomShuffle elements len g =+    let (rs, g') = rseq len g+    in  (S.shuffle elements rs, g')+  where+  -- | The sequence (r1,...r[n-1]) of numbers such that r[i] is an+  -- independent sample from a uniform random distribution+  -- [0..n-i]+  rseq :: RandomGen gen => Int -> gen -> ([Int], gen)+  rseq n g = second lastGen . unzip $ rseq' (n - 1) g+      where+        rseq' :: RandomGen gen => Int -> gen -> [(Int, gen)]+        rseq' i gen+          | i <= 0    = []+          | otherwise = let (j, gen') = randomR (0, i) gen+                        in  (j, gen') : rseq' (i - 1) gen'+        -- apply a function on the second element of a pair+        second :: (b -> c) -> (a, b) -> (a, c)+        second f (x,y) = (x, f y)+        -- the last returned random number generator+        lastGen [] = g   -- didn't use the generator yet+        lastGen (lst:[]) = lst+        lastGen gens = lastGen (drop 1 gens)++-- |Modify value with probability @p@. Return the unchanged value with probability @1-p@.+withProbability :: Double -> (a -> Rand a) -> (a -> Rand a)+withProbability p modify x = do+  t <- getDouble+  if t < p+     then modify x+     else return x
+ Moo/GeneticAlgorithm/Run.hs view
@@ -0,0 +1,252 @@+{-# LANGUAGE BangPatterns, Rank2Types #-}+{- |++Helper functions to run genetic algorithms and control iterations.++-}++module Moo.GeneticAlgorithm.Run (+  -- * Running algorithm+    runGA+  , runIO+  , nextGeneration+  , nextSteadyState+  , makeStoppable+  -- * Iteration control+  , loop, loopWithLog, loopIO+  , Cond(..), LogHook(..), IOHook(..)+) where++import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Selection (bestFirst)+import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.StopCondition+import Moo.GeneticAlgorithm.Utilities (doCrossovers, doNCrossovers)++import Data.Monoid (Monoid, mempty, mappend)+import Data.Time.Clock.POSIX (getPOSIXTime)+import Data.IORef (IORef, newIORef, readIORef, writeIORef)+import Control.Monad (liftM, when)++-- | Helper function to run the entire algorithm in the 'Rand' monad.+-- It takes care of generating a new random number generator.+runGA :: Rand [Genome a]             -- ^ function to create initial population+      -> ([Genome a] -> Rand b)       -- ^ genetic algorithm, see also 'loop' and 'loopWithLog'+      -> IO b                        -- ^ final population+runGA initialize ga = do+  rng <- newPureMT+  let (genomes0, rng') = runRandom initialize rng+  return $ evalRandom (ga genomes0) rng'++-- | Helper function to run the entire algorithm in the 'IO' monad.+runIO :: Rand [Genome a]                  -- ^ function to create initial population+      -> (IORef PureMT -> [Genome a] -> IO (Population a))+                                          -- ^ genetic algorithm, see also 'loopIO'+      -> IO (Population a)                -- ^ final population+runIO initialize gaIO = do+  rng <- newPureMT+  let (genomes0, rng') = runRandom initialize rng+  rngref <- newIORef rng'+  gaIO rngref genomes0++-- | Construct a single step of the genetic algorithm.+--+-- See "Moo.GeneticAlgorithm.Binary" and "Moo.GeneticAlgorithm.Continuous"+-- for the building blocks of the algorithm.+--+nextGeneration+    :: (ObjectiveFunction objectivefn a)+    => ProblemType          -- ^ a type of the optimization @problem@+    -> objectivefn          -- ^ objective function+    -> SelectionOp a        -- ^ selection operator+    -> Int                  -- ^ @elite@, the number of genomes to keep intact+    -> CrossoverOp a        -- ^ crossover operator+    -> MutationOp a         -- ^ mutation operator+    -> StepGA Rand a+nextGeneration problem objective selectOp elite xoverOp mutationOp =+  makeStoppable objective $ \pop -> do+    genomes' <- liftM (map takeGenome) $ withElite problem elite selectOp pop+    let top = take elite genomes'+    let rest = drop elite genomes'+    genomes' <- shuffle rest         -- just in case if @selectOp@ preserves order+    genomes' <- doCrossovers genomes' xoverOp+    genomes' <- mapM mutationOp genomes'+    return $ evalObjective objective (top ++ genomes')+++-- | Construct a single step of the incremental (steady-steate) genetic algorithm.+-- Exactly @n@ worst solutions are replaced with newly born children.+--+-- See "Moo.GeneticAlgorithm.Binary" and "Moo.GeneticAlgorithm.Continuous"+-- for the building blocks of the algorithm.+--+nextSteadyState+    :: (ObjectiveFunction objectivefn a)+    => Int                  -- ^ @n@, number of worst solutions to replace+    -> ProblemType          -- ^ a type of the optimization @problem@+    -> objectivefn          -- ^ objective function+    -> SelectionOp a        -- ^ selection operator+    -> CrossoverOp a        -- ^ crossover operator+    -> MutationOp a         -- ^ mutation operator+    -> StepGA Rand a+nextSteadyState n problem objective selectOp crossoverOp mutationOp =+    makeStoppable objective $ \pop -> do+      let popsize = length pop+      parents <- liftM (map takeGenome) (selectOp pop)+      children <- mapM mutationOp =<< doNCrossovers n parents crossoverOp+      let sortedPop = bestFirst problem pop+      let cpop = evalObjective objective children+      return . take popsize $ cpop ++ sortedPop+++-- | Wrap a population transformation with pre- and post-conditions+-- to indicate the end of simulation.+--+-- Use this function to define custom replacement strategies+-- in addition to 'nextGeneration' and 'nextSteadyState'.+makeStoppable+    :: (ObjectiveFunction objectivefn a, Monad m)+    => objectivefn+    -> (Population a -> m (Population a))  -- single step+    -> StepGA m a+makeStoppable objective onestep stop input = do+  let pop = either (evalObjective objective) id input+  if isGenomes input && evalCond stop pop+     then return $ StopGA pop   -- stop before the first iteration+     else do+       newpop <- onestep pop+       return $ if evalCond stop newpop+                then StopGA newpop+                else ContinueGA newpop+  where+    isGenomes (Left _) = True+    isGenomes (Right _) = False+++-- | Select @n@ best genomes, then select more genomes from the+-- /entire/ population (elite genomes inclusive). Elite genomes will+-- be the first in the list.+withElite :: ProblemType -> Int -> SelectionOp a -> SelectionOp a+withElite problem n select = \population -> do+  let elite = take n . eliteGenomes $ population+  selected <- select population+  return (elite ++ selected)+  where+    eliteGenomes = bestFirst problem++-- | Run strict iterations of the genetic algorithm defined by @step@.+-- Return the result of the last step.+{-# INLINE loop #-}+loop :: (Monad m)+     => Cond a+     -- ^ termination condition @cond@+     -> StepGA m a+     -- ^ @step@ function to produce the next generation+     -> [Genome a]+     -- ^ initial population+     -> m (Population a)+      -- ^ final population+loop cond step genomes0 = go cond (Left genomes0)+  where+    go cond !x = do+       x' <- step cond x+       case x' of+         (StopGA pop) -> return pop+         (ContinueGA pop) -> go (updateCond pop cond) (Right pop)++-- | GA iteration interleaved with the same-monad logging hooks.+{-# INLINE loopWithLog #-}+loopWithLog :: (Monad m, Monoid w)+     => LogHook a m w+     -- ^ periodic logging action+     -> Cond a+     -- ^ termination condition @cond@+     -> StepGA m a+     -- ^ @step@ function to produce the next generation+     -> [Genome a]+     -- ^ initial population+     -> m (Population a, w)+     -- ^ final population+loopWithLog hook cond step genomes0 = go cond 0 mempty (Left genomes0)+  where+    go cond !i !w !x = do+      x' <- step cond x+      case x' of+        (StopGA pop) -> return (pop, w)+        (ContinueGA pop) -> do+                         let w' = mappend w (runHook i pop hook)+                         let cond' = updateCond pop cond+                         go cond' (i+1) w' (Right pop)++    runHook !i !x (WriteEvery n write)+        | (rem i n) == 0 = write i x+        | otherwise      = mempty+++-- | GA iteration interleaved with IO (for logging or saving the+-- intermediate results); it takes and returns the updated random+-- number generator explicitly.+{-# INLINE loopIO #-}+loopIO+     :: [IOHook a]+     -- ^ input-output actions, special and time-dependent stop conditions+     -> Cond a+     -- ^ termination condition @cond@+     -> StepGA Rand a+     -- ^ @step@ function to produce the next generation+     -> IORef PureMT+     -- ^ reference to the random number generator+     -> [Genome a]+     -- ^ initial population @pop0@+     -> IO (Population a)+     -- ^ final population+loopIO hooks cond step rngref genomes0 = do+  rng <- readIORef rngref+  start <- realToFrac `liftM` getPOSIXTime+  (pop, rng') <- go start cond 0 rng (Left genomes0)+  writeIORef rngref rng'+  return pop+  where+    go start cond !i !rng !x = do+      stop <- (any id) `liftM` (mapM (runhook start i x) hooks)+      if (stop || either (const False) (evalCond cond) x)+         then return (asPopulation x, rng)+         else do+           let (x', rng') = runRandom (step cond x) rng+           case x' of+             (StopGA pop) -> return (pop, rng')+             (ContinueGA pop) ->+                 do+                   let i' = i + 1+                   let cond' = updateCond pop cond+                   go start cond' i' rng' (Right pop)++    -- runhook returns True to terminate the loop+    runhook _ i x (DoEvery n io) = do+             when ((rem i n) == 0) (io i (asPopulation x))+             return False+    runhook _ _ _ (StopWhen iotest)  = iotest+    runhook start _ _ (TimeLimit limit)  = do+             now <- realToFrac `liftM` getPOSIXTime+             return (now >= start + limit)++    -- assign dummy objective value to a genome+    dummyObjective :: Genome a -> Phenotype a+    dummyObjective g = (g, 0.0)++    asPopulation = either (map dummyObjective) id++-- | Logging to run every @n@th iteration starting from 0 (the first parameter).+-- The logging function takes the current generation count and population.+data (Monad m, Monoid w) => LogHook a m w =+    WriteEvery Int (Int -> Population a -> w)++-- | Input-output actions, interactive and time-dependent stop conditions.+data IOHook a+    = DoEvery { io'n :: Int, io'action :: (Int -> Population a -> IO ()) }+    -- ^ action to run every @n@th iteration, starting from 0;+    -- initially (at iteration 0) the objective value is zero.+    | StopWhen (IO Bool)+    -- ^ custom or interactive stop condition+    | TimeLimit { io't :: Double }+    -- ^ terminate iteration after @t@ seconds
+ Moo/GeneticAlgorithm/Selection.hs view
@@ -0,0 +1,158 @@+{- |++Selection operators for genetic algorithms.++-}++module Moo.GeneticAlgorithm.Selection+  (+    rouletteSelect+  , stochasticUniversalSampling+  , tournamentSelect+  -- ** Scaling and niching+  , withPopulationTransform+  , withScale+  , rankScale+  , withFitnessSharing+  -- ** Sorting+  , bestFirst+  ) where+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Niching (fitnessSharing)+++import Control.Monad (liftM, replicateM)+import Control.Arrow (second)+import Data.List (sortBy)+import Data.Function (on)++++-- | Apply given scaling or other transform to population before selection.+withPopulationTransform :: (Population a -> Population a) -> SelectionOp a -> SelectionOp a+withPopulationTransform transform select = \pop -> select (transform pop)+++-- | Transform objective function values before seletion.+withScale :: (Objective -> Objective) -> SelectionOp a -> SelectionOp a+withScale f select =+    let scale = map (second f)+    in  withPopulationTransform scale select++-- | Replace objective function values in the population with their+-- ranks.  For a population of size @n@, a genome with the best value+-- of objective function has rank @n' <= n@, and a genome with the+-- worst value of objective function gets rank @1@.+--+-- 'rankScale' may be useful to avoid domination of few super-genomes+-- in 'rouletteSelect' or to apply 'rouletteSelect' when an objective+-- function is not necessarily positive.+rankScale :: ProblemType -> Population a -> Population a+rankScale problem pop =+    let sorted = bestFirst (opposite problem) pop  -- worst first+        worst = takeObjectiveValue . head $ sorted+    in  ranks 1 worst sorted+    where+      ranks _ _ [] = []+      ranks rank worst ((genome,objective):rest)+          | worst == objective  = (genome,rank)   : ranks rank worst rest+          | otherwise           = (genome,rank+1) : ranks (rank+1) objective rest+      opposite Minimizing = Maximizing+      opposite Maximizing = Minimizing+++-- | A popular niching method proposed by D. Goldberg and+-- J. Richardson in 1987. The shared fitness of the individual is inversely+-- protoptional to its niche count.+-- The method expects the objective function to be non-negative.+--+-- An extension for minimization problems is implemented by+-- making the fitnes proportional to its niche count (rather than+-- inversely proportional).+--+-- Reference: Chen, J. H., Goldberg, D. E., Ho, S. Y., & Sastry,+-- K. (2002, July). Fitness inheritance in multiobjective+-- optimization. In Proceedings of the Genetic and Evolutionary+-- Computation Conference (pp. 319-326). Morgan Kaufmann Publishers+-- Inc..+withFitnessSharing ::+    (Phenotype a -> Phenotype a -> Double)  -- ^ distance function+    -> Double  -- ^ niche radius+    -> Double  -- ^ niche compression exponent @alpha@ (usually 1)+    -> ProblemType   -- ^ type of the optimization problem+    -> (SelectionOp a -> SelectionOp a)+withFitnessSharing dist r alpha ptype =+    withPopulationTransform (fitnessSharing dist r alpha ptype)+++-- |Objective-proportionate (roulette wheel) selection: select @n@+-- random items with each item's chance of being selected is+-- proportional to its objective function (fitness).+-- Objective function should be non-negative.+rouletteSelect :: Int -> SelectionOp a+rouletteSelect n xs = replicateM n roulette1+  where+  fs = map takeObjectiveValue xs+  xs' = zip xs (scanl1 (+) fs)+  sumScores = (snd . last) xs'+  roulette1 = do+    rand <- (sumScores*) `liftM` getDouble+    return $ (fst . head . dropWhile ((rand >) . snd)) xs'++-- |Performs tournament selection among @size@ individuals and+-- returns the winner. Repeat @n@ times.+tournamentSelect :: ProblemType  -- ^ type of the optimization problem+                 -> Int -- ^ size of the tournament group+                 -> Int -- ^ how many tournaments to run+                 -> SelectionOp a+tournamentSelect problem size n xs = replicateM n tournament1+  where+  tournament1 = do+    contestants <- randomSample size xs+    let winner = head $ bestFirst problem contestants+    return winner++-- | Stochastic universal sampling (SUS) is a selection technique+-- similar to roulette wheel selection. It gives weaker members a fair+-- chance to be selected, which is proportinal to their+-- fitness. Objective function should be non-negative.+stochasticUniversalSampling :: Int  -- ^ how many genomes to select+                            -> SelectionOp a+stochasticUniversalSampling n phenotypes = do+    let total = sum . map takeObjectiveValue $ phenotypes+    let step = total / (fromIntegral n)+    start <- getRandomR (0, step)+    let stops = [start + (fromIntegral i)*step | i <- [0..(n-1)]]+    let cumsums = scanl1 (+) (map takeObjectiveValue phenotypes)+    let ranges = zip (0:cumsums) cumsums+    -- for every stop select a phenotype with left cumsum <= stop < right cumsum+    return $ selectAtStops [] phenotypes stops ranges+  where+    selectAtStops selected _ [] _ = selected  -- no more stop points+    selectAtStops selected [] _ _ = selected  -- no more phenotypes+    selectAtStops selected phenotypes@(x:xs) stops@(s:ss) ranges@((l,r):lrs)+       | (l <= s && s < r) = selectAtStops (x:selected) phenotypes ss ranges  -- select a phenotype+       | s >= r = selectAtStops selected xs stops lrs  -- skip a phenotype AND the range+       | s < l  = error "stochasticUniformSampling: stop < leftSum"  -- should never happen+    selectAtStops _ _ _ _ = error "stochasticUniversalSampling: unbalanced ranges?"  -- should never happen++-- | Sort population by decreasing objective function (also known as+-- fitness for maximization problems). The genomes with the highest+-- fitness are put in the head of the list.+sortByFitnessDesc :: Population a -> Population a+sortByFitnessDesc = sortBy (flip compare `on` snd)++-- | Sort population by increasing objective function (also known as+-- cost for minimization problems). The genomes with the smallest+-- cost are put in the head of the list.+sortByCostAsc :: Population a -> Population a+sortByCostAsc = sortBy (compare `on` snd)++-- | Reorders a list of individual solutions,+-- by putting the best in the head of the list.+bestFirst :: ProblemType -> Population a -> Population a+bestFirst Maximizing = sortByFitnessDesc+bestFirst Minimizing = sortByCostAsc
+ Moo/GeneticAlgorithm/Statistics.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE BangPatterns #-}+{- |++Basic statistics for lists.++-}++module Moo.GeneticAlgorithm.Statistics+  ( average+  , variance+  , quantiles+  , median+  , iqr+  ) where++import Data.List (sort, foldl')++-- |Average+average :: (Num a, Fractional a) => [a] -> a+average = uncurry (/) . foldl' (\(!s, !c) x -> (s+x, c+1)) (0, 0)++-- |Population variance (divided by n).+variance :: (Floating a) => [a] -> a+variance xs = let (n, _, q) = foldr go (0, 0, 0) xs+              in  q / fromIntegral n+    where+    -- Algorithm by Chan et al.+    -- ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf+    go :: Floating a => a -> (Int, a, a) -> (Int, a, a)+    go x (n, sa, qa)+        | n == 0 = (1, x, 0)+        | otherwise =+            let na = fromIntegral n+                delta = x - sa/na+                sa' = sa + x+                qa' = qa + delta*delta*na/(na+1)+            in  (n + 1, sa', qa')+++-- | Compute empirical qunatiles (using R type 7 continuous sample quantile).+quantiles :: (Real a, RealFrac a)+             => [a]  -- ^ samples+             -> [a]  -- ^ probabilities in the range (0, 1)+             -> [a]  -- ^ estimated quantiles' values+quantiles xs probs =+    let xs' = sort xs+        n = length xs'+    in  map (quantile7 n xs') probs++-- | Estimate continuous quantile (like R's default type 7, SciPy (1,1), Excel).+quantile7 :: (Real a, RealFrac a)+             => Int       -- ^ @n@ the number of samples+             -> [a]       -- ^ @xs@ samples+             -> a         -- ^ @prob@ numeric probability (0, 1)+             -> a         -- ^ estimated quantile value+quantile7 n xs prob =+    let h = fromIntegral (n-1) * prob + 1+        i = floor h+        x1 = xs !! (i-1)+        x2 = xs !! (i)+    in  case (i >= n, i < 1) of+          (True, _) -> xs !! (i-1) -- prob >= 1+          (_, True) -> xs !! 0     -- prob < 0+          _         -> x1 + (h - fromIntegral i)*(x2 -x1)+++-- | Median+median :: (Real a, RealFrac a) => [a] -> a+median xs = head $ quantiles xs [0.5]+++-- | Interquartile range.+iqr :: (Real a, RealFrac a) => [a] -> a+iqr xs =+    let [q1,q2] = quantiles xs [0.25, 0.75]+    in  q2 - q1
+ Moo/GeneticAlgorithm/StopCondition.hs view
@@ -0,0 +1,30 @@+module Moo.GeneticAlgorithm.StopCondition where+++import Moo.GeneticAlgorithm.Types+++evalCond :: (Cond a) -> Population a -> Bool+evalCond (Generations n) _  = n <= 0+evalCond (IfObjective cond) p = cond . map takeObjectiveValue $ p+evalCond (GensNoChange n _ Nothing) _ = n <= 1+evalCond (GensNoChange n f (Just (prev, count))) p =+    let new = f . map takeObjectiveValue $ p+    in  (new == prev) && (count + 1 > n)+evalCond (Or c1 c2) x = evalCond c1 x || evalCond c2 x+evalCond (And c1 c2) x = evalCond c1 x && evalCond c2 x+++updateCond :: Population a -> Cond a -> Cond a+updateCond _ (Generations n) = Generations (n-1)+updateCond p (GensNoChange n f Nothing) =+     -- called 1st time _after_ the 1st iteration+    let counter = (Just (f (map takeObjectiveValue p), 1))+    in GensNoChange n f counter+updateCond p (GensNoChange n f (Just (v, c))) =+    let v' = f (map takeObjectiveValue p) in if v' == v+       then GensNoChange n f (Just (v, c+1))+       else GensNoChange n f (Just (v', 1))+updateCond p (And c1 c2) = And (updateCond p c1) (updateCond p c2)+updateCond p (Or c1 c2) = Or (updateCond p c1) (updateCond p c2)+updateCond _ c = c
+ Moo/GeneticAlgorithm/Types.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, GADTs, ExistentialQuantification #-}++module Moo.GeneticAlgorithm.Types+    (+    -- * Data structures+      Genome+    , Objective+    , Phenotype+    , Population+    , GenomeState(..)+    , takeObjectiveValue+    -- * GA operators+    , ProblemType (..)+    , ObjectiveFunction(..)+    , SelectionOp+    , CrossoverOp+    , MutationOp+    -- * Dummy operators+    , noMutation+    , noCrossover+    -- * Life cycle+    , StepGA+    , Cond(..)+    , PopulationState+    , StepResult(..)+    ) where++import Moo.GeneticAlgorithm.Random++-- | A genetic representation of an individual solution.+type Genome a = [a]++-- | A measure of the observed performance. It may be called cost+-- for minimization problems, or fitness for maximization problems.+type Objective = Double++-- | A genome associated with its observed 'Objective' value.+type Phenotype a = (Genome a, Objective)++-- | An entire population of observed 'Phenotype's.+type Population a = [Phenotype a]+++-- | 'takeGenome' extracts a raw genome from any type which embeds it.+class GenomeState gt a where+    takeGenome :: gt -> Genome a+++instance (a1 ~ a2) => GenomeState (Genome a1) a2 where+    takeGenome = id+++instance (a1 ~ a2) => GenomeState (Phenotype a1) a2 where+    takeGenome = fst+++takeObjectiveValue :: Phenotype a -> Objective+takeObjectiveValue = snd++-- | A type of optimization problem: whether the objective function+-- has to be miminized, or maximized.+data ProblemType = Minimizing | Maximizing deriving (Show, Eq)++-- | A function to evaluate a genome should be an instance of+-- 'ObjectiveFunction' class. It may be called a cost function for+-- minimization problems, or a fitness function for maximization+-- problems.+--+-- Some genetic algorithm operators, like 'rouletteSelect', require+-- the 'Objective' to be non-negative.+class ObjectiveFunction f a where+    evalObjective :: f -> [Genome a] -> Population a++-- | Evaluate fitness (cost) values genome per genome.+instance (a1 ~ a2) =>+    ObjectiveFunction (Genome a1 -> Objective) a2 where+        evalObjective f = map (\g -> (g, f g))++-- | Evaluate all fitness (cost) values at once.+instance (a1 ~ a2) =>+    ObjectiveFunction ([Genome a1] -> [Objective]) a2 where+        evalObjective f gs = zip gs (f gs)++-- | Evaluate fitness (cost) of all genomes, possibly changing their+-- order.+instance (a1 ~ a2) =>+    ObjectiveFunction ([Genome a1] -> [(Genome a1, Objective)]) a2 where+        evalObjective f gs = f gs++-- | A selection operator selects a subset (probably with repetition)+-- of genomes for reproduction via crossover and mutation.+type SelectionOp a = Population a -> Rand (Population a)++-- | A crossover operator takes some /parent/ genomes and returns some+-- /children/ along with the remaining parents. Many crossover+-- operators use only two parents, but some require three (like UNDX)+-- or more. Crossover operator should consume as many parents as+-- necessary and stop when the list of parents is empty.+type CrossoverOp a = [Genome a] -> Rand ([Genome a], [Genome a])++-- | A mutation operator takes a genome and returns an altered copy of it.+type MutationOp a = Genome a -> Rand (Genome a)++-- | Don't crossover.+noCrossover :: CrossoverOp a+noCrossover genomes = return (genomes, [])++-- | Don't mutate.+noMutation :: MutationOp a+noMutation = return+++-- | A single step of the genetic algorithm. See also 'nextGeneration'.+type StepGA m a = Cond a              -- ^ stop condition+                -> PopulationState a  -- ^ population of the current generation+                -> m (StepResult (Population a))  -- ^ population of the next generation+++-- | Iterations stop when the condition evaluates as @True@.+data Cond a =+      Generations Int                   -- ^ stop after @n@ generations+    | IfObjective ([Objective] -> Bool) -- ^ stop when objective values satisfy the @predicate@+    | forall b . Eq b => GensNoChange+      { c'maxgens ::  Int                 -- ^ max number of generations for an indicator to be the same+      , c'indicator ::  [Objective] -> b  -- ^ stall indicator function+      , c'counter :: Maybe (b, Int)       -- ^ a counter (initially @Nothing@)+      }                                 -- ^ terminate when evolution stalls+    | Or (Cond a) (Cond a)              -- ^ stop when at least one of two conditions holds+    | And (Cond a) (Cond a)             -- ^ stop when both conditions hold+++{-| On life cycle of the genetic algorithm:++>+>   [ start ]+>       |+>       v+>   (genomes) --> [calculate objective] --> (evaluated genomes) --> [ stop ]+>       ^  ^                                       |+>       |  |                                       |+>       |  `-----------.                           |+>       |               \                          v+>   [ mutate ]        (elite) <-------------- [ select ]+>       ^                                          |+>       |                                          |+>       |                                          |+>       |                                          v+>   (genomes) <----- [ crossover ] <-------- (evaluted genomes)+>++PopulationState can represent either @genomes@ or @evaluated genomed@.+-}+type PopulationState a = Either [Genome a] [Phenotype a]+++-- | A data type to distinguish the last and intermediate steps results.+data StepResult a = StopGA a | ContinueGA a deriving (Show)
+ Moo/GeneticAlgorithm/Utilities.hs view
@@ -0,0 +1,81 @@+{-# LANGUAGE BangPatterns #-}+{- |++Common utility functions.++-}++module Moo.GeneticAlgorithm.Utilities+  (+  -- * Non-deterministic functions+    getRandomGenomes+  , doCrossovers+  , doNCrossovers+) where++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Random+++import Control.Monad.Mersenne.Random+import Control.Monad (replicateM)+++-- | Generate @n@ random genomes made of elements in the+-- hyperrectangle ranges @[(from_i,to_i)]@. Return a list of genomes+-- and a new state of random number generator.+randomGenomes :: (Random a, Ord a)+              => PureMT  -- ^ random number generator+              -> Int     -- ^ n, number of genomes to generate+              -> [(a, a)]  -- ^ ranges for individual genome elements+              ->  ([Genome a], PureMT)+randomGenomes rng n ranges =+    let sortRange (r1,r2) = (min r1 r2, max r1 r2)+        ranges' = map sortRange ranges+    in  flip runRandom rng $+        replicateM n $ mapM getRandomR ranges'+++-- | Generate @n@ uniform random genomes with individual genome+-- elements bounded by @ranges@. This corresponds to random uniform+-- sampling of points (genomes) from a hyperrectangle with a bounding+-- box @ranges@.+getRandomGenomes :: (Random a, Ord a)+                         => Int  -- ^ @n@, how many genomes to generate+                         -> [(a, a)]  -- ^ ranges for individual genome elements+                         -> Rand ([Genome a])  -- ^ random genomes+getRandomGenomes n ranges =+    Rand $ \rng ->+        let (gs, rng') = randomGenomes rng n ranges+        in  R gs rng'+++-- | Crossover all available parents. Parents are not repeated.+doCrossovers :: [Genome a] -> CrossoverOp a -> Rand [Genome a]+doCrossovers []      _     = return []+doCrossovers parents xover = do+  (children', parents') <- xover parents+  if null children'+     then return []+     else do+       rest <- doCrossovers parents' xover+       return $ children' ++ rest+++-- | Produce exactly @n@ offsprings by repeatedly running the @crossover@+-- operator between randomly selected parents (possibly repeated).+doNCrossovers :: Int   -- ^ @n@, number of offsprings to generate+              -> [Genome a]  -- ^ @parents@' genomes+              -> CrossoverOp a  -- ^ @crossover@ operator+              -> Rand [Genome a]+doNCrossovers _ [] _ = return []+doNCrossovers n parents xover =+    doAnotherNCrossovers n []+  where+    doAnotherNCrossovers i children+        | i <= 0     = return . take n . concat $ children+        | otherwise  = do+      (children', _) <- xover =<< shuffle parents+      if (null children')+          then doAnotherNCrossovers 0 children  -- no more children+          else doAnotherNCrossovers (i - length children') (children':children)
+ README.md view
@@ -0,0 +1,145 @@+Moo+===++     ------------------------------------------------+    < Moo. Breeding Genetic Algorithms with Haskell. >+     ------------------------------------------------+            \   ^__^+             \  (oo)\_______+                (__)\       )\/\+                    ||----w |+                    ||     ||++++Features+--------++    |                       | Binary GA            | Continuous GA            |+    |-----------------------+----------------------+--------------------------|+    |Encoding               | binary bit-string    | sequence of real values  |+    |                       | Gray bit-string      |                          |+    |-----------------------+----------------------+--------------------------|+    |Initialization         |            random uniform                       |+    |                       |            constrained random uniform           |+    |                       |            arbitrary custom                     |+    |-----------------------+-------------------------------------------------|+    |Objective              |            minimization and maximiation         |+    |                       |            optional scaling                     |+    |                       |            optional ranking                     |+    |                       |            optional niching (fitness sharing)   |+    |-----------------------+-------------------------------------------------|+    |Selection              |            roulette                             |+    |                       |            stochastic universal sampling        |+    |                       |            tournament                           |+    |                       |            optional elitism                     |+    |                       |            optionally constrained               |+    |                       |            custom non-adaptive ^                |+    |-----------------------+-------------------------------------------------|+    |Crossover              |            one-point                            |+    |                       |            two-point                            |+    |                       |            uniform                              |+    |                       |            custom non-adaptive ^                |+    |                       +----------------------+--------------------------|+    |                       |                      | BLX-α (blend)            |+    |                       |                      | SBX (simulated binary)   |+    |                       |                      | UNDX (unimodal normally  |+    |                       |                      | distributed)             |+    |-----------------------+----------------------+--------------------------|+    |Mutation               | point                | Gaussian                 |+    |                       | asymmetric           |                          |+    |                       | constant frequency   |                          |+    |                       +----------------------+--------------------------|+    |                       |            custom non-adaptive ^                |+    |-----------------------+-------------------------------------------------|+    |Replacement            |            generational with elitism            |+    |                       |            steady state                         |+    |-----------------------+-------------------------------------------------|+    |Stop                   |            number of generations                |+    |condition              |            values of objective function         |+    |                       |            stall of objective function          |+    |                       |            custom or interactive (`loopIO`)     |+    |                       |            time limit (`loopIO`)                |+    |                       |            compound conditions (`And`, `Or`)    |+    |-----------------------+-------------------------------------------------|+    |Logging                |            pure periodic (any monoid)           |+    |                       |            periodic with `IO`                   |+    |-----------------------+-------------------------------------------------|+    |Constrainted           |            constrained initialization           |+    |optimization           |            constrained selection                |+    |                       |            death penalty                        |+    |-----------------------+-------------------------------------------------|+    |Multiobjective         |            NSGA-II                              |+    |optimization           |            constrained NSGA-II                  |+++`^` non-adaptive: any function which doesn't depend on generation number++There are other possible encodings which are possible to represent+with list-like genomes (`type Genome a = [a]`):++  * permutation encodings (`a` being an integer, or other `Enum` type)+  * tree encodings (`a` being a subtree type)+  * hybrid encodings (`a` being a sum type)+++Contributing+------------++There are many ways you can help developing the library:++  * I'm not a native speaker of English. If you are, please proof-read+    and correct the comments and the documentation.++  * Moo is designed with possibility of implementing custom genetic+    operators in mind. If you write new operators (`SelectionOp`,+    `CrossoverOp`, `MutationOp`) or replacement strategies+    (`StepGA`), consider contributing them to the library.+    In the comments please give a reference to an academic+    work which introduces or studies the method. Explain when or why+    it should be used. Provide tests and examples if possible.++  * Implementing some methods (like adaptive genetic algorithms) will+    require to change some library types. Please discuss your approach+    first.++  * Contribute examples. Solutions of known problems with known optima+    and interesting properties. Try to avoid examples which are too+    contrived.++++An example+----------++Minimizing [Beale's function][test-functions] (optimal value f(3, 0.5) = 0):++```haskell+import Moo.GeneticAlgorithm.Continuous+++beale :: [Double] -> Double+beale [x, y] = (1.5 - x + x*y)**2 + (2.25 - x + x*y*y)**2 + (2.625 - x + x*y*y*y)**2+++popsize = 101+elitesize = 1+tolerance = 1e-6+++selection = tournamentSelect Minimizing 2 (popsize - elitesize)+crossover = unimodalCrossoverRP+mutation = gaussianMutate 0.25 0.1+step = nextGeneration Minimizing beale selection elitesize crossover mutation+stop = IfObjective (\values -> (minimum values) < tolerance)+initialize = getRandomGenomes popsize [(-4.5, 4.5), (-4.5, 4.5)]+++main = do+  population <- runGA initialize (loop stop step)+  print (head . bestFirst Minimizing $ population)+```++For more examples, see [examples/](examples/) folder.++[test-functions]: http://en.wikipedia.org/wiki/Test_functions_for_optimization
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ Tests/Common.hs view
@@ -0,0 +1,87 @@+{-# LANGUAGE BangPatterns #-}+module Tests.Common where++import Moo.GeneticAlgorithm.Run+import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Random++import Data.List (foldl')+import Control.Monad (replicateM)+++type RealFunctionND = [Double] -> Double++data RealProblem = RealMinimize {+      minimizeFunction :: RealFunctionND      -- ^ function to minimize+    , minimizeVarRange :: [(Double, Double)]  -- ^ search space+    , minimizeSolution :: [Double]            -- ^ problem solution+    }+++-- Unit Gaussian mutation, 1/2 per genome+gauss sigma nvars =+    let p = 0.5/fromIntegral nvars+    in  gaussianMutate p sigma+++-- BLX-0.5 crossover+blxa = blendCrossover 0.5+++-- UNDX crossover+undx = unimodalCrossoverRP+++-- SBX crossover+sbx = simulatedBinaryCrossover 2+++randomGenomesReal :: Int -> [(Double,Double)] -> Rand [Genome Double]+randomGenomesReal popsize ranges = replicateM popsize randomGenome+    where+      randomGenome = mapM (\varRange -> getRandomR varRange) ranges+++data (ObjectiveFunction objectivefn a) => Solver objectivefn a = Solver {+      s'popsize :: Int+    , s'elitesize :: Int+    , s'objective :: objectivefn+    , s'select :: SelectionOp a+    , s'crossover :: CrossoverOp a+    , s'mutate :: MutationOp a+    , s'stopcond :: Cond a+    }+++-- default solver for real-valued problems+solverReal :: RealProblem -> Int -> Int -> CrossoverOp Double -> Cond Double+           -> Solver RealFunctionND Double+solverReal (RealMinimize f vranges sol) popsize elitesize crossover stopcond =+    let nvars = length vranges+        s = 0.1 * average (map (uncurry subtract) vranges)+        mutate = gauss s nvars+        select = tournamentSelect Minimizing 3 (popsize - elitesize)+    in  Solver popsize elitesize f select crossover mutate stopcond+++runSolverReal :: RealProblem+              -> Solver RealFunctionND Double+              -> IO (Population Double, Double)+              -- ^ final population and euclidean distance from the known solution+runSolverReal problem solver = do+    let ptype = Minimizing+    let init = randomGenomesReal (s'popsize solver) (minimizeVarRange problem)+    let step = nextGeneration  ptype (s'objective solver)+               (s'select solver) (s'elitesize solver)+               (s'crossover solver) (s'mutate solver)+    let ga   = loop (s'stopcond solver) step+    pop <- runGA init ga+    let best = takeGenome . head $ bestFirst ptype pop+    let dist = sqrt . sum . map (^2) $ zipWith (-) best (minimizeSolution problem)+    return (pop, dist)+++-- |Average+average :: (Num a, Fractional a) => [a] -> a+average = uncurry (/) . foldl' (\(!s, !c) x -> (s+x, c+1)) (0, 0)
+ Tests/Internals/TestConstraints.hs view
@@ -0,0 +1,84 @@+module Tests.Internals.TestConstraints where+++import Control.Monad (replicateM)+import Test.HUnit+import System.Random.Mersenne.Pure64 (pureMT)+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Selection+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Constraints+import Moo.GeneticAlgorithm.Binary+++testConstraints =+    TestList+    [ "constraint satisfaction" ~: do+        let gs =  [[-1],[0],[1],[2],[3::Int]]+        assertEqual ".<." [True, True, False, False, False] $+                    map (isFeasible [head .<. 1]) gs+        assertEqual ".<=." [True, True, True, False, False] $+                    map (isFeasible [head .<=. 1]) gs+        assertEqual ".>." [False, False, False, True, True] $+                    map (isFeasible [head .>. 1]) gs+        assertEqual ".>=." [False, False, True, True, True] $+                    map (isFeasible [head .>=. 1]) gs+        assertEqual ".==." [False, False, True, False, False] $+                    map (isFeasible [head .==. 1]) gs+        assertEqual "non-strict double inequality" [False, True, True, True, False] $+                    map (isFeasible [0 .<= head <=. 2]) gs+        assertEqual "strict double inequality" [False, False, True, False, False] $+                    map (isFeasible [0 .< head <. 2]) gs+    , "constrained initialization" ~: do+        let fI = fromIntegral :: Int -> Double+        let constraints = [ 1 .<= (fI . decodeBinary (0,255)) <=. 42 ]+        let n = 200+        let genomes = flip evalRandom (pureMT 1) $+                      getConstrainedBinaryGenomes constraints n 8+        assertEqual "exactly n genomes" n $+                    length genomes+        assertEqual "first constraint (<= .. <=)" True $+                    flip all genomes $ \bits ->+                        let x = fI $ decodeBinary (0,255) bits+                        in (x >= 0) && (x <= (42::Double))+    , "constrained selection (minimizing)" ~: do+        let n = 10+        let tournament2 = tournamentSelect Minimizing 2 n+        let constraints = [head .>=. 0, head .>=. (-1)]+        let ctournament = withConstraints constraints numberOfViolations Minimizing $+                          tournament2+        -- out of two solutions, one violates both constraints, another one only one+        let badvsugly = map (\x -> ([x], x)) [-1, -2]+        -- out of two solutions, one is feasible, the other is not+        let goodvsbad = map (\x -> ([x], x)) [0, -1]+        let result = flip evalRandom (pureMT 1) $ ctournament badvsugly+        assertEqual "lesser degree of violation is preferred"+                    (replicate n (-1.0)) $ (map (head . takeGenome) result)+        let result = flip evalRandom (pureMT 1) $ ctournament goodvsbad+        assertEqual "feasible solution is preferred"+                    (replicate n (0.0)) $ (map (head . takeGenome) result)+    , "numberOfViolations" ~: do+        let constraints = [head .>=. 0, head .>=. (-1)]+        assertEqual "1 violation" 1 $+                    numberOfViolations constraints [-1]+        assertEqual "2 violations" [2, 2] $+                    map (numberOfViolations constraints) [ [-2], [-3] ]+        assertEqual "no violations" 0 $+                    numberOfViolations constraints [0]+    , "degreeOfViolation" ~: do+        let constraints = [head .>=. 0, (negate . head) .<. (1)]+        assertEqual "no violation" 0 $+                    degreeOfViolation 2.0 0.5 constraints [0]+        assertEqual "1 non-strict violation" 0.25 $+                    degreeOfViolation 2.0 0.5 constraints [-0.5]+        assertEqual "1 non-strict and 1 strict violations" 1.5 $+                    degreeOfViolation 2.0 0.5 constraints [-1.0]+        assertEqual "non-strict double inequality"+                    [3.0,2.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,2.0,3.0] $+                    map (degreeOfViolation 1 0.5 [0 .<= head <=. 6]) $ map (:[]) [-3..9]+        assertEqual "strict double inequality"+                    [3.5,2.5,1.5,0.5,0.0,0.0,0.0,0.0,0.0,0.5,1.5,2.5,3.5] $+                    map (degreeOfViolation 1 0.5 [0 .< head <. 6]) $ map (:[]) [-3..9]+    ]
+ Tests/Internals/TestControl.hs view
@@ -0,0 +1,35 @@+module Tests.Internals.TestControl where+++import Test.HUnit+import System.Random.Mersenne.Pure64 (pureMT)+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Binary+import Moo.GeneticAlgorithm.Random+++instance (Eq a) => Eq (StepResult a) where+    (==) (StopGA xs) (StopGA ys) = xs == ys+    (==) (ContinueGA xs) (ContinueGA ys) = xs == ys+    (==) _ _ = False+++testControl =+    TestList+    [ "nextGeneration" ~: do+        let select = tournamentSelect Minimizing 2 8+        let objective = (fromIntegral . length) :: [Int] -> Double+        assertEqual "stop at initial population"  -- initial population is not changed+                    (StopGA [([1],1.0),([2,2],2.0)]) $+                    flip evalRandom (pureMT 1) $+                             (nextGeneration Minimizing objective select 0 noCrossover noMutation)+                             (Generations 0) (Left [[1],[2,2]])+        assertEqual "do at least one step"+                    (ContinueGA [([1],1.0),([1],1.0),([1],1.0),([1],1.0)+                                ,([1],1.0),([1],1.0),([1],1.0),([1],1.0)]) $+                    flip evalRandom (pureMT 1) $+                             (nextGeneration Minimizing objective select 0 noCrossover noMutation)+                             (Generations 1) (Left [[1],[2,2]])+    ]
+ Tests/Internals/TestCrossover.hs view
@@ -0,0 +1,83 @@+module Tests.Internals.TestCrossover where+++import Test.HUnit+import System.Random.Mersenne.Pure64 (pureMT)+import Data.List (group)+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Crossover+import Moo.GeneticAlgorithm.Random++++testCrossover =+    TestList+    [ "do N crossovers" ~: do+        let genomes = [[1,1,1,1],[0,0,0,0]] :: [[Int]]+        let result4 = flip evalRandom (pureMT 1) $+                      doNCrossovers 4 genomes (onePointCrossover 0.5)+        let expected4 = [[0,0,1,1],[1,1,0,0],[0,0,0,1],[1,1,1,0]]+        assertEqual "4 crossovers" expected4 result4+        let result3 = flip evalRandom (pureMT 1) $+                      doNCrossovers 3 genomes (onePointCrossover 0.5)+        let expected3 = [[0,0,1,1],[1,1,0,0],[0,0,0,1]]+        assertEqual "3 crossovers" expected3 result3+    , "do all crossovers" ~: do+        let genomes = [[1,1,1,1],[0,0,0,0]] :: [[Int]]+        let result = flip evalRandom (pureMT 1) $+                     doCrossovers genomes (onePointCrossover 0.5)+        let expected = [[1,1,1,0],[0,0,0,1]]+        assertEqual "all crossovers (2 genomes)" expected result+        let genomes3 = [[1,1,1,1],[0,0,0,0],[2,2,2,2]] :: [[Int]]+        -- genes from the last "celibate" genome are lost+        let result3 = filter (==2) . concat . map concat . flip map [0..100] $+                      \i -> flip evalRandom (pureMT i) $+                      doCrossovers genomes (onePointCrossover 1.0)+        assertEqual "discard last genomes without a pair" [] result3+    , "simple crossover" ~: do+        let ones = replicate 8 1+        let zeros = replicate 8 0+        let genomes = [ones, zeros]+        let n = 1000+        assertEqual "exactly one crossover point" True $+                    all (<=2) . map (length . group) $+                        flip evalRandom (pureMT 1) (doNCrossovers n genomes (onePointCrossover 1))+    , "simple crossover" ~: do+        let ones = replicate 8 1+        let zeros = replicate 8 0+        let genomes = [ones, zeros]+        let n = 1000+        assertEqual "exactly one crossover point" True $+                    all (<=2) . map (length . group) $+                        flip evalRandom (pureMT 1) (doNCrossovers n genomes (onePointCrossover 1))+    , "two-point crossover" ~: do+        let ones = replicate 8 1+        let zeros = replicate 8 0+        let genomes = [ones, zeros]+        let n = 1000+        assertEqual "exactly two crossover point" True $+                    all (<=3) . map (length . group) $+                        flip evalRandom (pureMT 1) (doNCrossovers n genomes (twoPointCrossover 1))+    , "uniform crossover" ~: do+        assertEqual "change all points"+                    ([[0,0,0,0,0,0,0,0,0,0],[1,1,1,1,1,1,1,1,1,1]],[]) $+                    flip evalRandom (pureMT 1) $+                             (uniformCrossover 1) [replicate 10 1,replicate 10 (0::Int)]+        assertEqual "change nothing"+                    ([[1,1,1,1,1,1,1,1,1,1],[0,0,0,0,0,0,0,0,0,0]],[]) $+                    flip evalRandom (pureMT 1) $+                             (uniformCrossover 0) [replicate 10 1,replicate 10 (0::Int)]+        let n = 1000+        let mu = 0.5*n+        let sigma = sqrt(n*0.5*(1-0.5))  -- normal approx to binomial distribution+        let genomes = [ replicate (round n) 1+                      , replicate (round n) 0]+        let xover = uniformCrossover 0.5 :: CrossoverOp Double+        let mkChildren = doNCrossovers 1000 genomes xover :: Rand [Genome Double]+        let children = flip evalRandom (pureMT 1) mkChildren :: [Genome Double]+        assertEqual "change approximately half" True $+                    all (\s -> (s >= mu - 4*sigma && s <= mu + 4*sigma)) . map sum $+                        children+    ]
+ Tests/Internals/TestFundamentals.hs view
@@ -0,0 +1,45 @@+module Tests.Internals.TestFundamentals where+++import Control.Monad (replicateM)+import Test.HUnit+import System.Random.Mersenne.Pure64 (pureMT)+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Multiobjective.Types+import Moo.GeneticAlgorithm.Random+import Moo.GeneticAlgorithm.Binary+++testFundamentals =+    TestList+    [ "takeGenome" ~: do+        assertEqual "raw genome" [True] $ takeGenome [True]+        assertEqual "phenotype" [True,True] $ takeGenome ([True,True], 42.0::Double)+        assertEqual "multiobjective phenotype" [False] $ takeGenome ([False], [42.0::Double])+    , "withProbability" ~: do+        assertEqual "probability 0" 42 $+                    flip evalRandom (pureMT 1) $+                    withProbability 0 (return . (+1)) 42+        assertEqual "probability 1" 43 $+                    flip evalRandom (pureMT 1) $+                    withProbability 1 (return . (+1)) 42+    , "pointMutate" ~: do+        let zeros = map (=='1') (replicate 16 '0')+        assertEqual "just 1 bit is changed" (replicate 10 1) $+                    flip evalRandom (pureMT 1) $+                         replicateM 10 $+                         return . length . filter id =<< pointMutate 1 zeros+    , "asymmetricMutate" ~: do+        let g = map (=='1') "0000000011111111"  -- 8 bits set+        assertEqual "flip all zeros" 16 $+                    flip evalRandom (pureMT 1) $+                         return . length . filter id =<< asymmetricMutate 1 0 g+        assertEqual "flip all ones" 0 $+                    flip evalRandom (pureMT 1) $+                         return . length . filter id =<< asymmetricMutate 0 1 g+        assertEqual "flip all" 8 $+                    flip evalRandom (pureMT 1) $+                         return . length . filter id =<< asymmetricMutate 1 1 g+    ]
+ Tests/Internals/TestMultiobjective.hs view
@@ -0,0 +1,147 @@+module Tests.Internals.TestMultiobjective where+++import Test.HUnit+import Control.Monad (forM_)+import Data.Function (on)+import Data.List (sortBy)+import qualified Data.Set as Set+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Multiobjective.Types+import Moo.GeneticAlgorithm.Multiobjective.NSGA2+import Moo.GeneticAlgorithm.Constraints+++import System.Random.Mersenne.Pure64 (pureMT)+++dummyGenome :: [Objective] -> MultiPhenotype Double+dummyGenome ovs = (ovs, ovs)+++testMultiobjective =+    TestList+    [ "domination predicate" ~: do+        let problems = [Minimizing, Maximizing, Minimizing]+        let worst = dummyGenome [100, 0, 100]+        let good1 = dummyGenome [0, 50, 50]+        let good23 = dummyGenome [50, 100, 0]+        let best = dummyGenome [0, 100, 0]+        assertEqual "good dominates worst"+                    True (domination problems good1 worst)+        assertEqual "good23 doesn't dominate good1"+                    False (domination problems good23 good1)+        assertEqual "good1 doesn't dominate good23"+                    False (domination problems good1 good23)+        assertEqual "best dominates good23"+                    True (domination problems best good23)+        assertEqual "worst doesn't dominate best"+                    False (domination problems worst best)+    , "constraint-domination predicate" ~: do+        let problems = [Minimizing]+        let constraints = [head .>=. 2, head .>=. 4]+        let feasible = dummyGenome [4]+        let worse = dummyGenome [5]  -- also feasible+        let infeasible = dummyGenome [3]+        let infeasible2 = dummyGenome [1]+        let dominates = constrainedDomination constraints numberOfViolations problems+        assertEqual "feasible dominates infeasible" [True, True, False] $+                    [ feasible `dominates` infeasible+                    , feasible `dominates` infeasible2+                    , infeasible `dominates` feasible ]+        assertEqual "less-infeasible dominates more-infeasible" [True,False] $+                    [ infeasible `dominates` infeasible2+                    , infeasible2 `dominates` infeasible ]+        assertEqual "better dominates worse" [True, False] $+                    [ feasible `dominates` worse+                    , worse `dominates` feasible ]+    , "non-dominated sort" ~: do+        let dominatesFn = domination [Minimizing, Minimizing]+        let genomes = [ ([1], [2, 2]), ([2], [3, 2]), ([2,2], [2,3])+                      , ([3], [1,1.5]), ([3,3], [1.5, 0.5]), ([4], [0,0::Double])]+        assertEqual "non-dominated fronts"+                    [[[4]],[[3],[3,3]],[[1]],[[2],[2,2]]]+                    (map (map fst) $ nondominatedSort dominatesFn genomes)+    , "non-dominated sort (singleton fronts)" ~: do+        let dominates1 = domination [Maximizing]+        let genomes1 = map (\x -> ([x],[x])) [3,1,2]+        assertEqual "singleton fronts"+                    [[3],[2],[1]]+                    (map (map (head . fst)) $ nondominatedSort dominates1 genomes1)+    , "calculate crowding distance" ~: do+        let inf = 1.0/0.0 :: Double+        assertEqual "two points" [inf, inf] $ crowdingDistances [[1],[2]]+        assertEqual "4 points" [inf, 2.5, inf, 2.0] $ crowdingDistances [[1.0], [2.0], [4.0], [3.5]]+        assertEqual "4 points 2D" [inf, 2.0, inf, 0.75, 2.0] $+                    crowdingDistances [[3,1], [1.75,1.75], [1,3], [2,2], [2.125,2.125]]+    , "rank with crowding" ~: do+        let dominatesFn = domination [Minimizing, Minimizing]+        let gs = map (\x -> ([], x)) [[2,1],[1,2],[3,1],[1.9,1.9],[1,3]]+        let rs = concat $ rankAllSolutions dominatesFn gs+        let inf = 1.0/0.0 :: Double+        assertEqual "non-dom ranks" [1,1,1,2,2]+                    (map rs'nondominationRank rs)+        assertEqual "in-front crowding distance" [inf, inf, 2.0, inf, inf]+                    (map rs'localCrowdingDistnace rs)+    , "calculate all objectives for all genomes" ~: do+        let genomes = [[8, 2], [2.0, 1.0], [1.0, 2.0], [4,4]]+        let objectives = [(Minimizing, sum), (Maximizing, product)]+                       :: [(ProblemType, [Double] -> Double)]+        let correct = [([8.0,2.0],[10.0,16.0]),([2.0,1.0],[3.0,2.0])+                      ,([1.0,2.0],[3.0,2.0]),([4.0,4.0],[8.0,16.0])]+        assertEqual "two objective functions" correct $+                    evalAllObjectives objectives genomes+    , "NSGA-II ranking with crowding" ~: do+        let dominatesFn = domination [Minimizing, Minimizing]+        let mp = [ (Minimizing, (!!0))+                 , (Minimizing, (!!1))+                 ] :: [(ProblemType, [Double] -> Double)]+        let gs = [ [5,1], [1,5], [2,4], [3,3]  -- first front+                 , [6,6]                       -- third front+                 , [6,2], [5,3], [4,4], [2,6]  -- second front+                 ] :: [[Double]]+        let expected7 = [(([5.0,1.0],[5.0,1.0]),1.0)+                        ,(([1.0,5.0],[1.0,5.0]),1.0) -- order is preserved in the first front:+                        ,(([2.0,4.0],[2.0,4.0]),1.0) -- [2,4] is more crowded than [3,3]+                        ,(([3.0,3.0],[3.0,3.0]),1.0) -- but it doesn't matter for full fronts+                        ,(([6.0,2.0],[6.0,2.0]),2.0)+                        ,(([2.0,6.0],[2.0,6.0]),2.0) -- is front boundary point, and goes before [4,4]+                        ,(([4.0,4.0],[4.0,4.0]),2.0) -- is less crowded than [5,3]+                        -- [5,3] is more crowded and is truncated+                        -- [6,6] is in the third front and is truncated+                        ]+        let result7 = nsga2Ranking dominatesFn mp 7 gs+        assertEqual "7 solutions" expected7 result7+    , "NSGA-II ranking (output length)" ~: do+        let dominatesFn = domination [Minimizing, Minimizing]+        let mp = [ (Minimizing, (!!0))+                 , (Minimizing, (!!1))+                 ] :: [(ProblemType, [Double] -> Double)]+        let gs = [ [5,1], [1,5], [2,4], [3,3]  -- first front+                 , [6,6]                       -- third front+                 , [6,2], [5,3], [4,4], [2,6]  -- second front+                 ] :: [[Double]]+        forM_ [0..(length gs)] $ \n -> do+          assertEqual (show n ++ " solutions") n $+                      length (nsga2Ranking dominatesFn mp n gs)+        assertEqual "max # of solutions" (length gs) $+                    length (nsga2Ranking dominatesFn mp maxBound gs)+    , "two NSGA-II steps" ~: do+        let mp = [ (Minimizing, (!!0))+                 , (Minimizing, (!!1))+                 ] :: [(ProblemType, [Double] -> Double)]+        let gs = [ [5,1], [1,5], [2,4], [3,3]  -- first front+                 , [6,6]                       -- third front+                 , [6,2], [5,3], [4,4], [2,6]  -- second front+                 ] :: [[Double]]+        let expected = [([1.0,5.0],1.0),([5.0,1.0],1.0),([1.0,5.0],1.0)+                       ,([5.0,1.0],1.0),([3.0,3.0],1.0),([3.0,3.0],1.0)+                       ,([2.0,4.0],1.0),([2.0,4.0],1.0),([1.0,5.0],1.0)]+        let result = flip evalRandom (pureMT 1) $+                     loop (Generations 1)+                     (stepNSGA2bt mp noCrossover noMutation) gs+        assertEqual "solutions and ranking" (Set.fromList expected) (Set.fromList result)+    ]
+ Tests/Internals/TestSelection.hs view
@@ -0,0 +1,66 @@+module Tests.Internals.TestSelection where+++import Test.HUnit+import System.Random.Mersenne.Pure64 (pureMT)+import Control.Monad (replicateM)+++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Selection+import Moo.GeneticAlgorithm.Random++++dummyGenome :: Objective -> Phenotype ()+dummyGenome objval = ([], objval)+++testSelection =+    TestList+    [ "tournamentSelect" ~: do+        let resultMin = flip evalRandom (pureMT 1) $+                        tournamentSelect Minimizing 3 4 $+                        map dummyGenome [3,2,4]+        let resultMax = flip evalRandom (pureMT 1) $+                        tournamentSelect Maximizing 2 3 $+                        map dummyGenome [2,3]+        assertEqual "4 times best of 3" [2,2,2,2] $+                    map takeObjectiveValue resultMin+        assertEqual "3 times best of 2" [3,3,3] $+                    map takeObjectiveValue resultMax+    , "tournamentSelect (10 times best of 4, seed 1)" ~: do+        let times = 10+        let tsize = 4+        let genomes = map dummyGenome [1..10]+        let resultMany = flip evalRandom (pureMT 1) $+                         tournamentSelect Maximizing tsize times $+                         genomes+        let objVals = map takeObjectiveValue resultMany+        -- take the same samples again with the same see+        let samples = flip evalRandom (pureMT 1) $+                           replicateM times (randomSample tsize genomes)+        assertEqual "maximum is selected every time" (replicate times True)  $+                    zipWith (\selected xs -> selected == (maximum . map takeObjectiveValue $ xs))+                            objVals samples+    , "rouletteSelect" ~: do+       let gs = map dummyGenome [1, 9]+       let tries = 100 * 1000 :: Int+       let numOfNines = length . filter (==9.0) . map takeObjectiveValue+                        . flip evalRandom (pureMT 1) $ rouletteSelect tries $ gs+       assertEqual "9 is selected from [1,9] 90% of time" 90 (numOfNines `div` 1000)+    , "stochasticUniversalSampling" ~: do+        let gs = map dummyGenome [2,1]+        let selected = flip evalRandom (pureMT 1) $+                       stochasticUniversalSampling 12 gs+        assertEqual "counts are fitness proportional" [4, 8] $+             map length [ (filter ((==1) . takeObjectiveValue) selected)+                        , (filter ((==2) . takeObjectiveValue) selected) ]+    , "rankScale" ~: do+        let expected = [([30.0],1.0),([10.0],2.0),([2.0],3.0),([0.0],4.0)]+        let expectedMax = [([0.0],1.0),([2.0],2.0),([10.0],3.0),([30.0],4.0)]+        let result = rankScale Minimizing (map (\x -> ([x],x)) [2,10,0,30])+        let resultMax = rankScale Maximizing (map (\x -> ([x],x)) [2,10,0,30])+        assertEqual "min problem" expected result+        assertEqual "max problem" expectedMax resultMax+    ]
+ Tests/Problems/Rosenbrock.hs view
@@ -0,0 +1,91 @@+{- Minimize Rosenbrock function using real-valued genetic algorithm.+   Optimal value x* = (1,...,1). F(x*) = 0.+-}++module Tests.Problems.Rosenbrock where++import Test.HUnit++import Text.Printf+import Data.List (intercalate)+import System.IO (hPutStrLn, stderr)+import Control.Monad (replicateM)++import Tests.Common++import Moo.GeneticAlgorithm.Types+import Moo.GeneticAlgorithm.Selection+import Moo.GeneticAlgorithm.Run+import Moo.GeneticAlgorithm.Random++pr _ = return ()+-- pr = hPutStrLn stderr+++rosenbrock :: [Double] -> Double+rosenbrock xs = sum . map f $ zip xs (drop 1 xs)+  where+   f (x1, x2) = 100 * (x2 - x1^2)^2 + (x1 - 1)^2+++testRosenbrock = TestList+  [ "Rosenbrock 2D GM/UNDX/500 gens" ~: do+      let tolerance = 1e-3  -- solution error+      let maxiters = 500+      let problem = RealMinimize rosenbrock [(-10,10),(-20,20)]  [1,1]+      let solver = solverReal problem 101 11 undx (Generations maxiters)+      (pop, dist) <- runSolverReal problem solver+      let best = takeGenome . head $ bestFirst Minimizing pop+      pr ""+      pr $ "best:    " ++ (intercalate " " (map (printf "%.5f") best))+      pr $ "error:   " ++ (printf "%.5g" dist)+      assertBool ("error >= " ++ show tolerance) (dist < tolerance)+  , "Rosenbrock 2D GM/SBX/min residual, max 500 gens" ~: do+      let tolerance = 1e-6  -- objective residual+      let maxiters = 500+      let problem = RealMinimize rosenbrock [(-20,20),(-20,20)] [1,1]+      let stop = Generations maxiters `Or` IfObjective ((>= -tolerance) . maximum)+      let solver = solverReal problem 101 11 sbx stop+      (pop, dist) <- runSolverReal problem solver+      let best = head $ bestFirst Minimizing pop+      let bestG = takeGenome best+      let bestF = takeObjectiveValue  best+      pr ""+      pr $ "best:    " ++ (intercalate " " (map (printf "%.5f") bestG))+      pr $ "residual: " ++ (printf "%.5g" bestF)+      assertBool ("residual < " ++ show (negate tolerance)) (bestF >= -tolerance)+  , "Rosenbrock 2D GM/BLX-0.5/min residual, max 500 gens" ~: do+      let tolerance = 1e-3  -- solution error+      let maxiters = 500+      let problem = RealMinimize rosenbrock [(-20,20),(-20,20)] [1,1]+      let stop = Generations maxiters+      let solver = solverReal problem 101 11 blxa stop+      (pop, dist) <- runSolverReal problem solver+      let bestG = takeGenome . head $ bestFirst Minimizing pop+      pr ""+      pr $ "best:    " ++ (intercalate " " (map (printf "%.5f") bestG))+      pr $ "error:   " ++ (printf "%.5g" dist)+      assertBool ("error >= " ++ show tolerance) (dist < tolerance)+  , "Rosenbrock 2D GM/UNDX/GensNoChange 10" ~: do+      let maxiters = 5000+      let popsize = 101+      let elite = 11+      let nochange = 10+      let select = tournamentSelect Minimizing 3 (popsize - elite)+      let stop = (GensNoChange nochange (round.(*1e3).maximum) Nothing) `Or` (Generations maxiters)+      let step = nextGeneration Minimizing rosenbrock select elite undx (gauss 1.0 2)+      let log = WriteEvery 1 (\_ p -> [minimum . map takeObjectiveValue $ p])+      let ga = loopWithLog log stop step+      let init = replicateM popsize . replicateM 2 $ getRandomR (-10,10)++      (pop, hist) <- runGA init ga++      let best = takeGenome . head $ bestFirst Minimizing pop+      pr ""+      pr $ "best:    " ++ (intercalate " " (map (printf "%.5f") best))+      let lastbest = take nochange (reverse hist)+      pr $ "last best: "+      mapM_ pr (map show $ reverse lastbest)+      assertBool "false positive on GensNoChange"+                     (all id $ zipWith (==) lastbest (drop 1 lastbest))+  ]
+ examples/ExampleMain.hs view
@@ -0,0 +1,154 @@+-- | The boring part common to many examples: command line options+-- and pretty-printing the results.+module ExampleMain+    ( exampleMain+    , ExampleDefaults(..)+    , exampleDefaults+    ) where+++import Moo.GeneticAlgorithm.Binary+import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Multiobjective+import Moo.GeneticAlgorithm.Statistics+++import Control.Monad (liftM, when)+import Data.List (intercalate)+import System.Console.GetOpt+import System.Environment (getArgs, getProgName)+import System.Exit (exitSuccess)+import Text.Printf+++data Flag = RunGenerations Int+          | PrintBest Bool+          | PrintStats Bool+          | DumpAll Bool+          | ShowHelp+            deriving (Show, Eq)+++data ExampleDefaults = ExampleDefaults+    { numGenerations :: Int+    , printBest :: Bool+    , printStats :: Bool+    , dumpAll :: Bool+    } deriving (Show, Eq)+++exampleDefaults = ExampleDefaults {+                    numGenerations = 100+                  , printBest = True+                  , printStats = False+                  , dumpAll = False+                  }+++exampleOptions :: ExampleDefaults -> [OptDescr Flag]+exampleOptions c =+    [ Option "gn" ["generations"]+                 (ReqArg (RunGenerations . read) "N")+                 ("number of generations (default: " ++ show (numGenerations c) ++ ")")+    , Option "b" ["best"]+                 (NoArg $ PrintBest True)+                 ("print the best solution" ++ (isDefault (printBest c)))+    , Option ""  ["no-best"]+                 (NoArg $ PrintBest False)+                 ("don't print the best solution" ++ (isDefault (not . printBest $ c)))+    , Option "d" ["dump"]+                 (NoArg $ DumpAll True)+                 ("dump the entire population and its objective values" ++ isDefault (dumpAll c))+    , Option ""  ["no-dump"]+                 (NoArg $ DumpAll False)+                 ("don't dump the entire population" ++ isDefault (not . dumpAll $ c))+    , Option "s" ["stats"]+                 (NoArg $ PrintStats True)+                 ("print population statistics" ++ isDefault (printStats c))+    , Option ""  ["no-stats"]+                 (NoArg $ PrintStats False)+                 ("don't print population statistics" ++ isDefault (not . printStats $ c))+    , Option "h" ["help"]+                 (NoArg ShowHelp)+                 "show help"+    ]+   where+   isDefault :: Bool -> String+   isDefault True = " (default)"+   isDefault False = ""+++updateDefaults :: ExampleDefaults -> [Flag] -> ExampleDefaults+updateDefaults d (RunGenerations n:opts) = updateDefaults (d { numGenerations = n }) opts+updateDefaults d (PrintBest b:opts) = updateDefaults (d { printBest = b }) opts+-- --stats overrid --dump, and vice versa+updateDefaults d (DumpAll b:opts) =+    let ps = printStats d+    in  flip updateDefaults opts (d { dumpAll = b, printStats = ps && (not b)})+updateDefaults d (PrintStats b:opts) =+    let da = dumpAll d+    in  flip updateDefaults opts (d { printStats = b, dumpAll = da && (not b)})+updateDefaults d [] = d++++printHeader conf = do+  when (printStats conf) $ putStrLn "# best, median"+  when (dumpAll conf) $ putStrLn "# x1, x2, ..., objective1, objective2, ..."+++printSnapshot conf sorted = do+  when (printBest conf) $+    if null sorted+       then putStrLn "# no solutions"+       else putStrLn $ "# best found: " ++ fmtPt (head sorted)++  when (printStats conf) $ do+    printHeader conf+    let ovs = map takeObjectiveValue sorted+    let obest = head ovs+    let omedian = median ovs+    putStrLn $ fmtXs " " [obest, omedian]++  when (dumpAll conf) $ do+    printHeader conf+    -- print the best solution last;+    -- (for scatter-plotting it above the others)+    flip mapM_ (reverse sorted) $ \p -> putStrLn $ fmtPtOneline p+    putStrLn ""++  where++    fmtPt :: (Show a, Real a, PrintfArg a) => Phenotype a -> String+    fmtPt (xs, v) = (printf "%.3g @ [" v) ++ fmtXs ", " xs ++ "]"++    fmtPtOneline :: (Show a, Real a, PrintfArg a) => Phenotype a -> String+    fmtPtOneline p = let xs = map (fromRational.toRational) . takeGenome $ p+                         vs = [takeObjectiveValue p]+                     in  fmtXs " " $ xs ++ vs++    fmtXs :: (Show a, Real a, PrintfArg a) => String -> [a] -> String+    fmtXs sep xs =  intercalate sep $ map (printf "%.3g") xs++++-- | Run a genetic algorithm defined by @problemtype@, and @step@.+-- Process command line options to change the number of iterations+-- and logging behaviour.+exampleMain :: (Show a, Real a, PrintfArg a)+            => ExampleDefaults -> ProblemType -> Rand [Genome a] -> StepGA Rand a -> IO ()+exampleMain defaults problemtype initialize step = do++  let options = exampleOptions defaults+  (opts, args, msgs) <- liftM (getOpt Permute options) getArgs+  when (ShowHelp `elem` opts) $ do+    progname <- getProgName+    let header = "usage: " ++ progname ++ " [options]\n\nOptions:\n"+    putStrLn (usageInfo header options)+    exitSuccess++  let conf = updateDefaults defaults opts+  let gens = numGenerations conf+  result <- runGA initialize (loop (Generations gens) step)+  let sorted = bestFirst problemtype $ result+  printSnapshot conf sorted
+ examples/README.md view
@@ -0,0 +1,35 @@+Examples+========++Examples of real-coded GAs:++  * [beale.hs](beale.hs) Beale function+    (basic GA)++  * [rosenbrock.hs](rosenbrock.hs) Rosenbrock function+    (basic GA with pure logging)++  * [schaffer2.hs](schaffer2.hs) Schaffer function #2+    (steady-state GA with niching)++  * [cp_sphere2.hs](cp_sphere2.hs) constrained 2D sphere function over a convex set+    (GA with a death penalty)++  * [cp_himmelblau.hs](cp_himmelblau.hs) constrained Himmelblau function over a non-convex set+    (GA with niching and constrained tournament selection)++  * [mop_minsum_maxprod.hs](mop_minsum_maxprod.hs) a simple multiobjective problem+    (basic NSGA-II)++  * [mop_kursawe.hs](mop_kursawe.hs) Kursawe function, a multiobjective problem+    with a discontinuous and non-convex Pareto set+    (constrained NSGA-II)++  * [mop_constr2.hs](mop_constr2.hs) a constrained multiobjective problem from (Deb, 2002),+    a part of the unconstrained Pareto-optimal region is not feasible+    (constrained NSGA-II with niching)++Examples of binary GAs:++  * [knapsack.hs](knapsack.hs) 0-1 knapsack problem.+    (A basic GA with logging in IO and time limit)
+ examples/beale.hs view
@@ -0,0 +1,27 @@+{- Minimize Beale function using real-valued genetic algorithm.+   Optimal value x* = [3, 0.5]. F(x*) = 0.+-}++import Moo.GeneticAlgorithm.Continuous+++beale :: [Double] -> Double+beale [x, y] = (1.5 - x + x*y)**2 + (2.25 - x + x*y*y)**2 + (2.625 - x + x*y*y*y)**2+++popsize = 101+elitesize = 1+tolerance = 1e-6+++selection = tournamentSelect Minimizing 2 (popsize - elitesize)+crossover = unimodalCrossoverRP+mutation = gaussianMutate 0.25 0.1+step = nextGeneration Minimizing beale selection elitesize crossover mutation+stop = IfObjective (\values -> (minimum values) < tolerance)+initialize = getRandomGenomes popsize [(-4.5, 4.5), (-4.5, 4.5)]+++main = do+  population <- runGA initialize (loop stop step)+  print (head . bestFirst Minimizing $ population)
+ examples/cp_himmelblau.hs view
@@ -0,0 +1,64 @@+{- Constrained Himmelblau function over a non-convex set.+++Test problem #1 from Deb, K. (2000). An efficient constraint+handling method for genetic algorithms. Computer methods in applied+mechanics and engineering, 186(2), 311-338.++Unconstrained optimum: (3,2)+Constrained optimum: (2.246826, 2.381865)++Running and visualizing in bash/zsh:++N=100 ; ghc --make cp_himmelblau && ./cp_himmelblau -b -d -g $N > output.txt && ( gnuplot -persist <<< "set view map; unset key ; set isosamples 100 ; set logscale cb ; splot [0:6][0:6] (x**2 + y - 11)**2 + (x + y*y - 7)**2 w pm3d, 'output.txt' u 1:2:(0) w p lc 2 pt 4; set xlabel 'x' ; set ylabel 'y' ; set title 'generation $N' ; replot " ; head -1 output.txt)+++-}+++import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Constraints+++import ExampleMain+++import Data.Function (on)+++f :: [Double] -> Double+f [x, y] = (x**2 + y - 11)**2 + (x + y**2 - 7)**2+xvar [x,_] = x+yvar [_,y] = y+g1 [x,y] = 4.84 - (x-0.05)**2 - (y-2.5)**2+g2 [x,y] = x**2 + (y-2.5)**2 - 4.84+++constraints = [ 0 .<= xvar <=. 6+              , 0 .<= yvar <=. 6+              , g1 .>=. 0+              , g2 .>=. 0 ]+++popsize = 100+initialize = getRandomGenomes popsize [(0,6),(0,6)]+select = withFitnessSharing (distance2 `on` takeGenome) 0.025 1 Minimizing $+         withConstraints constraints (degreeOfViolation 1.0 0.0) Minimizing $+         tournamentSelect Minimizing 2 popsize+step = withFinalDeathPenalty constraints $+       nextGeneration Minimizing f select 0+       (simulatedBinaryCrossover 0.5)+       (gaussianMutate 0.05 0.025)+++{-+-- exampleMain takes care of command line options and pretty printing.+-- If you don't need that, a bare bones main function looks like this:++main = do+  results <- runGA initialize (loop (Generations 100) step)+  print . head . bestFirst Minimizing $ results++-}+main = exampleMain (exampleDefaults { numGenerations = 100 } )+       Minimizing initialize step
+ examples/cp_sphere2.hs view
@@ -0,0 +1,46 @@+{- Constrained problem++   min (x^2 + y^2)++   with x + y >= 1.++-}++import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Constraints+++import ExampleMain+++f :: [Double] -> Double+f [x, y] = x*x + y*y+++constraints = [ sum .>=. 1 ]+++popsize = 100+++initialize = getRandomGenomes popsize [(-10,10),(-5,5)]+select = tournamentSelect Minimizing 2 popsize+crossover = unimodalCrossoverRP+mutation = noMutation+++step = withDeathPenalty constraints $+       nextGeneration Minimizing f select 2 crossover mutation+++{-+-- exampleMain takes care of command line options and pretty printing.+-- If you don't need that, a bare bones main function looks like this:++main = do+  results <- runGA initialize (loop (Generations 25) step)+  print . head . bestFirst Minimizing $ results++-}+main = exampleMain (exampleDefaults { numGenerations = 25 })+       Minimizing initialize step
+ examples/knapsack.hs view
@@ -0,0 +1,102 @@+{-+  The 0-1 knapsack problem. Given a set of items with given weight and value,+  choose which items to put into collection to maximize collection value+  with given maximum weight constraint.++  It is a binary genetic algorithm. This example interleaves computation+  with logging in IO monad, and terminates by reaching a time limit.++  To run:++      ghc --make knapsack.hs+      ./knapsack > output.txt++  To visualize the output in gnuplot:++      % gnuplot+      > plot 'output.txt' u 1:2 w l t 'median value', '' u 1:3 w l t 'best value' lt 3+-}++import Moo.GeneticAlgorithm.Binary++import Control.Monad+import Data.List (intercalate)++type Weight = Int+type Value = Int+type Problem = [(Weight, Value)]++items = 42+itemWeight = (1,9 :: Weight)+itemValue = (0,9 :: Value)+maxTotalWeight = items*2 :: Weight++popsize = 11+elitesize = 1++-- fitness function to maximize+totalValue :: Problem -> [Bool] -> Objective+totalValue things taken = fromIntegral . snd $ totalWeithtAndValue things taken++totalWeithtAndValue :: Problem -> Genome Bool -> (Weight, Value)+totalWeithtAndValue things taken = sumVals (0,0) $ zip taken things+  where+    sumVals (totalW, totalV) ((True, (w,v)):rest)  -- item is taken+        | totalW + w > maxTotalWeight  = (totalW, totalV)  -- weight limit exceeded+        | otherwise                    = sumVals (totalW+w,totalV+v) rest+    sumVals acc ((False, _):rest)      = sumVals acc rest+    sumVals (totalW, totalV) []        = (totalW, totalV)  -- all items in the knapsack+++select = tournamentSelect Maximizing 2 (popsize-elitesize)++-- generate items to choose from: [(weight, value)]+randomProblem ::  IO Problem+randomProblem = do+  rng <- newPureMT+  return . flip evalRandom rng $ do+                      weights <- replicateM items $ getRandomR itemWeight+                      values <- replicateM items $ getRandomR itemValue+                      return $ zip weights values++geneticAlgorithm :: Problem -> IO (Population Bool)+geneticAlgorithm things = do+  let initialize = replicateM popsize $ replicateM items getRandom+  let fitness = totalValue things+  let nextGen = nextGeneration Maximizing fitness select elitesize+                          (onePointCrossover 0.5) (pointMutate 0.5)+  runIO initialize $ loopIO+         [DoEvery 10 logStats, TimeLimit 0.1]  -- stop after 100 ms+         (Generations maxBound)  -- effectively, forever; unless an IOHook condition triggers+         nextGen++  where++    logStats :: Int -> Population Bool -> IO ()+    logStats iterno pop = do+      when (iterno == 0) $+           putStrLn "# generation medianValue bestValue"+      let gs = map takeGenome . bestFirst Maximizing $ pop  -- genomes+      let best = head gs+      let median = gs !! (length gs `div` 2)+      let bvalue = snd $ totalWeithtAndValue things best+      let mvalue = snd $ totalWeithtAndValue things median+      putStrLn $ intercalate " " (map show [iterno, mvalue, bvalue])+++main = do+  things <- randomProblem+  pop <- geneticAlgorithm things+  putStrLn "# final population:"+  let best = takeGenome . head . bestFirst Maximizing $ pop+  let bestthings = zip best things+  let taken = intercalate ", " . map (showItem . snd) $ filter fst bestthings+  let left = intercalate ", " . map (showItem . snd) $ filter (not . fst) bestthings+  putStrLn $ showPop pop+  putStrLn $ "# taken: " ++ taken+  putStrLn $ "# left: " ++ left++  where+    showPop = intercalate "\n" . map showG+    showG (bs,v) = "# " ++ (concatMap (show . fromEnum) bs) ++ " " ++ show v+    showItem (w, v) = "$" ++ show v ++ "/" ++ show w ++ "oz"
+ examples/mop_constr2.hs view
@@ -0,0 +1,46 @@+{- CONSTR2 problem from (Deb. 2002).+   A part of the unconstrained Pareto-optimal region is not feasible.+-}+++import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Constraints+import Moo.GeneticAlgorithm.Multiobjective+++popsize = 100+generations = 100+++mop :: MultiObjectiveProblem ([Double] -> Double)+mop = [ (Minimizing, \[x1,_] -> x1)+      , (Minimizing, \[x1,x2] -> (1+x2)/x1) ]+++constraints = [ 0.1 .<= x1 <=. 1.0+              , 0.0 .<= x2 <=. 5.0+              , g1 .>=. 6.0+              , g2 .>=. 1.0 ]+  where+    x1 [x,_] = x+    x2 [_,y] = y+    g1 [x1,x2] = 9*x1 + x2+    g2 [x1,x2] = 9*x1 - x2++++initialize = getConstrainedGenomes constraints popsize [(0.1,1.0),(0.0,5.0)]+tournament = tournamentSelect Minimizing 2 popsize+++step :: StepGA Rand Double+step = stepConstrainedNSGA2 constraints (degreeOfViolation 1 0)+       mop tournament (blendCrossover 0.1) noMutation -- (gaussianMutate 0.5 0.5)+++main = do+  result <- runGA initialize $ loop (Generations generations) step+  let solutions = map takeGenome $ takeWhile ((<= 10.0) . takeObjectiveValue) result+  let ovs = map takeObjectiveValues $ evalAllObjectives mop solutions+  flip mapM_ ovs $ \[x1,x2] ->+      putStrLn $ show x1 ++ "\t" ++ show x2
+ examples/mop_kursawe.hs view
@@ -0,0 +1,49 @@+{- Kursawe function++A multiobjective optimization problem with a discontinuous and+non-convex Pareto front.++Kursawe, F. (1991). A variant of evolution strategies for vector+optimization. In Parallel Problem Solving from Nature+(pp. 193-197). Springer Berlin Heidelberg.++-}+++import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Constraints+import Moo.GeneticAlgorithm.Multiobjective+++n = 3+popsize = 100+generations = 100+++mop :: MultiObjectiveProblem ([Double] -> Double)+mop = [ (Minimizing,+         \xs -> sum (map (\i -> -10*exp(-0.2*sqrt(((xs!!i)**2 + (xs!!(i+1))**2)))) [0..(n-2)]))+      , (Minimizing,+         \xs -> sum (map (\x -> abs(x)**0.8 + 5*sin(x**3)) xs)) ]+++constraints :: [Constraint Double Double]+constraints = [ (-5.0) .<= (!!0) <=. 5.0+              , (-5.0) .<= (!!1) <=. 5.0+              , (-5.0) .<= (!!2) <=. 5.0 ]+++initialize = getRandomGenomes popsize (replicate 3 (-5.0, 5.0))+++step :: StepGA Rand Double+step = stepConstrainedNSGA2bt constraints (degreeOfViolation 1 0)+       mop unimodalCrossoverRP (gaussianMutate 0.01 0.5)+++main = do+  result <- runGA initialize $ loop (Generations generations) step+  let solutions = map takeGenome $ takeWhile ((<= 10.0) . takeObjectiveValue) result+  let ovs = map takeObjectiveValues $ evalAllObjectives mop solutions+  flip mapM_ ovs $ \[x1,x2] ->+      putStrLn $ show x1 ++ "\t" ++ show x2
+ examples/mop_minsum_maxprod.hs view
@@ -0,0 +1,52 @@+{- A simple multiobjective problem:++  minimize f_1 = x + y+  maximize f_2 = x * y++  s.t. x >= 0, y >=0. -}+++import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Constraints+import Moo.GeneticAlgorithm.Multiobjective+++import Text.Printf (printf)+++mop :: MultiObjectiveProblem ([Double] -> Double)+mop = [ (Minimizing, sum :: [Double] -> Double)+      , (Maximizing, product)]+++constraints = [ xvar .>=. 0+              , yvar .>=. 0 ]+xvar [x,_] = x+yvar [_,y] = y+++genomes :: [[Double]]+genomes = [[3,3], [9,1], [1,4], [2,2], [1,9], [4,1], [1,1], [4,2]]+++popsize :: Int+popsize = 50+step :: StepGA Rand Double+step = withDeathPenalty constraints $+       stepNSGA2bt mop noCrossover (gaussianMutate 0.1 0.5)+++main = do+  putStrLn $ "# population size: " ++ show popsize+  result <- runGA+            (return . take popsize . cycle $ genomes) $+            (loop (Generations 100) step)+  putStrLn $ "# best:"+  printPareto result+++printPareto result = do+  let paretoGenomes = map takeGenome . takeWhile ((== 1.0) . takeObjectiveValue) $ result+  let paretoObjectives = map takeObjectiveValues $ evalAllObjectives mop paretoGenomes+  putStr $ unlines $+       map (\[x,y] -> printf "%12.3f\t%12.3f" x y ) paretoObjectives
+ examples/rosenbrock.hs view
@@ -0,0 +1,121 @@+{- Minimize Rosenbrock function using real-valued genetic algorithm.+   Optimal value x* = (1,...,1). F(x*) = 0.++   It is a real-values genetic algorithm. The user may choose a+   mutation and crossover operators.  This example uses hooks to save+   evolution history.++   To run:++       ghc --make rosenbrock.hs+       ./rosenbrock gm undx > output.txt++   To visualize the output in gnuplot:++       % gnuplot+       > set logscale y ; set xlabel 'generation' ;+       > plot 'output.txt' u 1:2 w l t 'median', '' u 1:3 w l t 'best' lt 3+++-}++import Moo.GeneticAlgorithm.Continuous++import Control.Monad+import Data.List+import System.Environment (getArgs)+import System.Exit (exitWith, ExitCode(..))+import Text.Printf (printf)++rosenbrock :: [Double] -> Double+rosenbrock xs = sum . map f $ zip xs (drop 1 xs)+  where+   f (x1, x2) = 100.0 * (x2 - x1^(2::Int))^(2::Int) + (x1 - 1)^(2::Int)++nvariables = 3+xrange = (-30.0, 30.0)+popsize = 100+precision = 1e-5+maxiters = 4000 :: Int+elitesize = 10++-- Rosenbrock function is minimized+objective :: [Double] -> Objective+objective xs = rosenbrock xs++-- selection: tournament selection+select = tournamentSelect Minimizing 3 (popsize-elitesize)++-- Gaussian mutation, mutate fraction @genomeschanged@ of the population+gm genomeschanged =+    let p = 1.0 - (1.0 - genomeschanged)**(1.0 / fromIntegral nvariables)+        s = 0.01*(snd xrange - fst xrange)+    in  gaussianMutate p s++mutationOps = [ ("gm", gm 0.33) ]++-- BLX-0.5 crossover+blxa = blendCrossover 0.5+-- UNDX crossover+undx = unimodalCrossoverRP+-- SBX crossover+sbx = simulatedBinaryCrossover 2++crossoverOps = [ ("blxa", blxa), ("undx", undx), ("sbx", sbx) ]++printUsage = do+  putStrLn usage+  exitWith (ExitFailure 1)+  where+  usage = intercalate " " [ "rosenbrock", mops, xops ]+  mops = intercalate "|" (map fst mutationOps)+  xops = intercalate "|" (map fst crossoverOps)++logStats = WriteEvery 10 $ \iterno pop ->+             let pop' =  bestFirst Minimizing pop+                 bestobjval = takeObjectiveValue $ head pop'+                 medianobjval = takeObjectiveValue $ pop' !! (length pop' `div` 2)+             in  [(iterno, medianobjval, bestobjval)]++printStats :: [(Int, Objective, Objective)] -> IO ()+printStats stats = do+  printf "# %-10s %15s %15s\n" "generation" "median" "best"+  flip mapM_ stats $ \(iterno, median, best) ->+      printf "%12d %15.3g %15.3g\n" iterno median best++geneticAlgorithm mutate crossover = do+  -- initial population+  let initialize = replicateM popsize $ replicateM nvariables (getRandomR xrange)+  let stop = IfObjective ((<= precision) . minimum) `Or` Generations maxiters+  let step = nextGeneration Minimizing objective select elitesize crossover mutate+  --+  let ga = loopWithLog logStats stop step+  runGA initialize ga+++printBest :: Population Double -> IO ()+printBest pop = do+  let bestGenome = takeGenome . head $ bestFirst Minimizing pop+  let vals = map (\x -> printf "%.5f" x) bestGenome+  putStrLn $ "# best solution: " ++ (intercalate ", " vals)++-- usage: rosenbrock mutationOperator crossoverOperator+main = do+  args <- getArgs+  conf <- case args of+           []       -> return (lookup "gm" mutationOps, lookup "undx" crossoverOps)+           (m:x:[]) -> return (lookup m mutationOps, lookup x crossoverOps)+           _        -> printUsage+  case conf of+    (Just mutate, Just crossover) -> do+       (pop, stats) <- geneticAlgorithm mutate crossover+       printStats stats+       printBest pop+       -- exit status depends on convergence+       let bestF = takeObjectiveValue . head $ bestFirst Minimizing pop+       if (abs bestF <= precision)+          then exitWith ExitSuccess+          else do+            printf "# failed to converge: best residual=%.5g, target=%g\n" bestF precision+            exitWith (ExitFailure 2)  -- failed to find a solution+    _ -> printUsage
+ examples/schaffer2.hs view
@@ -0,0 +1,39 @@+{- Schaffer function #2. Minimium at (0,0), equal to 0. -}++import Moo.GeneticAlgorithm.Continuous+import Moo.GeneticAlgorithm.Constraints+++import Data.Function (on)+++import ExampleMain+++schafferN2 :: [Double] -> Double+schafferN2 [x,y] = 0.5 + (sin(x*x-y*y)**2 - 0.5)/(1+0.001*(x*x+y*y))**2+xvar [x,_] = x+yvar [_,y] = y+++popsize = 100+initialize = getRandomGenomes popsize (replicate 2 (-100,100))+select = withFitnessSharing (distance2 `on` takeGenome) 1.0 1 Minimizing $+         tournamentSelect Minimizing 2 popsize+crossover = unimodalCrossoverRP+mutate = gaussianMutate 0.05 0.1+step = nextSteadyState (popsize `div` 100) Minimizing schafferN2+       select crossover mutate+++{-+-- exampleMain takes care of command line options and pretty printing.+-- If you don't need that, a bare bones main function looks like this:++main = do+  results <- runGA initialize (loop (Generations 1000) step)+  print . head . bestFirst Minimizing $ results++-}+main = exampleMain (exampleDefaults { numGenerations = 1000 })+       Minimizing initialize step
+ moo-tests.hs view
@@ -0,0 +1,26 @@+import System.Exit+import Test.HUnit++import Tests.Internals.TestFundamentals (testFundamentals)+import Tests.Internals.TestControl (testControl)+import Tests.Internals.TestSelection (testSelection)+import Tests.Internals.TestCrossover (testCrossover)+import Tests.Internals.TestConstraints (testConstraints)+import Tests.Internals.TestMultiobjective (testMultiobjective)+import Tests.Problems.Rosenbrock (testRosenbrock)++allTests = TestList+  [ testFundamentals+  , testControl+  , testSelection+  , testCrossover+  , testConstraints+  , testRosenbrock+  , testMultiobjective+  ]++main = do+  result <- runTestTT allTests+  if (errors result + failures result) > 0+    then exitFailure+    else exitSuccess
+ moo.cabal view
@@ -0,0 +1,124 @@+name:               moo+category:           AI, Algorithms, Optimisation, Optimization+build-type:         Simple+version:            1.0+synopsis:           Genetic algorithm library+description:        Moo library provides building blocks to build custom+                    genetic algorithms in Haskell. They can be used to+                    find solutions to optimization and search problems.+                    .+                    Variants supported out of the box: binary (using+                    bit-strings) and continuous (real-coded).+                    Potentially supported variants: permutation,+                    tree, hybrid encodings (require customizations).+                    .+                    Binary GAs: binary and Gray encoding; point mutation;+                    one-point, two-point, and uniform crossover.+                    Continuous GAs: Gaussian mutation; BLX-α, UNDX, and+                    SBX crossover.+                    Selection operators: roulette, and tournament;+                    with optional niching and scaling.+                    Replacement strategies: generational with elitism+                    and steady state.+                    Constrained optimization: random constrained+                    initialization, death penalty, constrained+                    selection without a penalty function.+                    Multi-objective optimization: NSGA-II+                    and constrained NSGA-II.++license:            BSD3+License-file:       LICENSE+maintainer:         Sergey Astanin <s.astanin@gmail.com>+author:             Sergey Astanin <s.astanin@gmail.com>+stability:          experimental+homepage:           http://www.github.com/astanin/moo/+cabal-version:       >=1.8+extra-source-files: README.md+                  , examples/README.md+                  , examples/ExampleMain.hs+                  , examples/beale.hs+                  , examples/cp_himmelblau.hs+                  , examples/cp_sphere2.hs+                  , examples/knapsack.hs+                  , examples/mop_constr2.hs+                  , examples/mop_kursawe.hs+                  , examples/mop_minsum_maxprod.hs+                  , examples/rosenbrock.hs+                  , examples/schaffer2.hs+++Library+    build-depends:      base >=4 && < 5+                      , monad-mersenne-random+                      , mersenne-random-pure64+                      , gray-code >= 0.2.1+                      , random >= 0.1+                      , random-shuffle >= 0.0.2+                      , mtl >= 2+                      , time+                      , array+    ghc-options:        -Wall -fno-warn-name-shadowing -fno-warn-orphans+    exposed-modules:    Moo.GeneticAlgorithm+                      , Moo.GeneticAlgorithm.Binary+                      , Moo.GeneticAlgorithm.Constraints+                      , Moo.GeneticAlgorithm.Continuous+                      , Moo.GeneticAlgorithm.Multiobjective+                      , Moo.GeneticAlgorithm.Random+                      , Moo.GeneticAlgorithm.Run+                      , Moo.GeneticAlgorithm.Statistics+                      , Moo.GeneticAlgorithm.Types+    other-modules:      Moo.GeneticAlgorithm.Crossover+                      , Moo.GeneticAlgorithm.LinAlg+                      , Moo.GeneticAlgorithm.Multiobjective.NSGA2+                      , Moo.GeneticAlgorithm.Multiobjective.Types+                      , Moo.GeneticAlgorithm.Selection+                      , Moo.GeneticAlgorithm.StopCondition+                      , Moo.GeneticAlgorithm.Utilities+                      , Moo.GeneticAlgorithm.Crossover+                      , Moo.GeneticAlgorithm.Niching++Test-Suite moo-tests+  Type:                 exitcode-stdio-1.0+  Main-Is:              moo-tests.hs+  Other-Modules:        Tests.Common+                      , Tests.Internals.TestControl+                      , Tests.Internals.TestCrossover+                      , Tests.Internals.TestFundamentals+                      , Tests.Internals.TestMultiobjective+                      , Tests.Internals.TestSelection+                      , Tests.Internals.TestConstraints+                      , Tests.Problems.Rosenbrock+                      , Moo.GeneticAlgorithm+                      , Moo.GeneticAlgorithm.Binary+                      , Moo.GeneticAlgorithm.Constraints+                      , Moo.GeneticAlgorithm.Continuous+                      , Moo.GeneticAlgorithm.Crossover+                      , Moo.GeneticAlgorithm.Niching+                      , Moo.GeneticAlgorithm.Run+                      , Moo.GeneticAlgorithm.Random+                      , Moo.GeneticAlgorithm.Utilities+                      , Moo.GeneticAlgorithm.LinAlg+                      , Moo.GeneticAlgorithm.Multiobjective+                      , Moo.GeneticAlgorithm.Multiobjective.NSGA2+                      , Moo.GeneticAlgorithm.Multiobjective.Types+                      , Moo.GeneticAlgorithm.Selection+                      , Moo.GeneticAlgorithm.Statistics+                      , Moo.GeneticAlgorithm.StopCondition+                      , Moo.GeneticAlgorithm.Types+  Build-Depends:+      moo+    , base < 5+    , HUnit+    , random >= 0.1+    , random-shuffle >= 0.0.2+    , monad-mersenne-random+    , mersenne-random-pure64+    , gray-code >= 0.2.1+    , mtl+    , time+    , array+    , containers++source-repository head+  type:     git+  location: git://github.com/astanin/moo.git