moo 1.0 → 1.2

raw patch · 42 files changed

+4518/−4196 lines, 42 filesdep +MonadRandomdep +paralleldep +vectordep −monad-mersenne-randomsetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: MonadRandom, parallel, vector

Dependencies removed: monad-mersenne-random

API changes (from Hackage documentation)

Files

• LICENSE +32/−32
• Moo/GeneticAlgorithm.hs +146/−146
• Moo/GeneticAlgorithm/Binary.hs +264/−240
• Moo/GeneticAlgorithm/Constraints.hs +290/−290
• Moo/GeneticAlgorithm/Continuous.hs +276/−250
• Moo/GeneticAlgorithm/Crossover.hs +64/−64
• Moo/GeneticAlgorithm/LinAlg.hs +31/−31
• Moo/GeneticAlgorithm/Multiobjective.hs +21/−18
• Moo/GeneticAlgorithm/Multiobjective/Metrics.hs +143/−0
• Moo/GeneticAlgorithm/Multiobjective/NSGA2.hs +496/−495
• Moo/GeneticAlgorithm/Multiobjective/Types.hs +45/−45
• Moo/GeneticAlgorithm/Niching.hs +55/−55
• Moo/GeneticAlgorithm/Random.hs +143/−111
• Moo/GeneticAlgorithm/Run.hs +260/−252
• Moo/GeneticAlgorithm/Selection.hs +163/−158
• Moo/GeneticAlgorithm/Statistics.hs +75/−75
• Moo/GeneticAlgorithm/StopCondition.hs +30/−30
• Moo/GeneticAlgorithm/Types.hs +158/−157
• Moo/GeneticAlgorithm/Utilities.hs +77/−81
• README.md +176/−145
• Setup.hs +2/−2
• Tests/Common.hs +81/−87
• Tests/Internals/TestConstraints.hs +84/−84
• Tests/Internals/TestControl.hs +35/−35
• Tests/Internals/TestCrossover.hs +87/−83
• Tests/Internals/TestFundamentals.hs +45/−45
• Tests/Internals/TestMultiobjective.hs +147/−147
• Tests/Internals/TestSelection.hs +67/−66
• Tests/Problems/Rosenbrock.hs +91/−91
• examples/ExampleMain.hs +154/−154
• examples/README.md +72/−35
• examples/beale.hs +26/−26
• examples/cp_himmelblau.hs +64/−64
• examples/cp_sphere2.hs +46/−46
• examples/knapsack.hs +102/−102
• examples/mop_constr2.hs +45/−45
• examples/mop_kursawe.hs +48/−48
• examples/mop_minsum_maxprod.hs +57/−51
• examples/rosenbrock.hs +121/−121
• examples/schaffer2.hs +39/−39
• moo-tests.hs +26/−26
• moo.cabal +134/−124

@@ -1,32 +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.+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.

@@ -1,146 +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+{-# 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

@@ -1,240 +1,264 @@-{-# 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+{-# LANGUAGE CPP #-}
+{-# 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.
+#if MIN_VERSION_base(4, 7, 0)
+encodeGray :: (FiniteBits b, Bits b, Integral b) => (b, b) -> b -> [Bool]
+#else
+encodeGray :: (Bits b, Integral b) => (b, b) -> b -> [Bool]
+#endif
+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@.
+#if MIN_VERSION_base(4, 7, 0)
+decodeGray :: (FiniteBits b, Bits b, Integral b) => (b, b) -> [Bool] -> b
+#else
+decodeGray :: (Bits b, Integral b) => (b, b) -> [Bool] -> b
+#endif
+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.
+#if MIN_VERSION_base(4, 7, 0)
+encodeBinary :: (FiniteBits b, Bits b, Integral b) => (b, b) -> b -> [Bool]
+#else
+encodeBinary :: (Bits b, Integral b) => (b, b) -> b -> [Bool]
+#endif
+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@.
+#if MIN_VERSION_base(4, 7, 0)
+decodeBinary :: (FiniteBits b, Bits b, Integral b) => (b, b) -> [Bool] -> b
+#else
+decodeBinary :: (Bits b, Integral b) => (b, b) -> [Bool] -> b
+#endif
+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
+
+#if MIN_VERSION_base(4, 7, 0)
+encodeWithCode :: (FiniteBits b, Bits b, Integral b) => ([Bool] -> [Bool]) -> (b, b) -> b -> [Bool]
+#else
+encodeWithCode :: (Bits b, Integral b) => ([Bool] -> [Bool]) -> (b, b) -> b -> [Bool]
+#endif
+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') ++ (repeat False)
+
+#if MIN_VERSION_base(4, 7, 0)
+decodeWithCode :: (FiniteBits b, Bits b, Integral b) => ([Bool] -> [Bool]) -> (b, b) -> [Bool] -> b
+#else
+decodeWithCode :: (Bits b, Integral b) => ([Bool] -> [Bool]) -> (b, b) -> [Bool] -> b
+#endif
+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

@@ -1,290 +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)+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 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 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)

@@ -1,250 +1,276 @@-{-# 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)+{-# LANGUAGE LambdaCase #-}
+{-# 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
+  , uniformGenomes
+
+  -- * 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.Utilities (getRandomGenomes)
+
+
+-- | Generate at most @popsize@ genomes uniformly distributed in @ranges@.
+uniformGenomes :: Int -> [(Double,Double)] -> [Genome Double]
+uniformGenomes popsize ranges =
+    let dims = map (uncurry subtract) ranges :: [Double]
+        ndims = length dims :: Int
+        vol = product dims
+        mdim = vol ** (1.0/fromIntegral ndims) :: Double
+        msamples = (fromIntegral popsize) ** (1.0/fromIntegral ndims) :: Double
+        ptsPerDim = map (\d -> round $ d*msamples/mdim) dims :: [Int]
+        ptsInLastDims = product $ drop 1 ptsPerDim :: Int
+        ptsInFirstDim = popsize `div` ptsInLastDims :: Int
+        ptsPerDim' = ptsInFirstDim : (drop 1 ptsPerDim) :: [Int]
+        linspaces = zipWith linspace ranges ptsPerDim' :: [[Double]]
+    in  sproduct [[]] linspaces
+  where
+    linspace :: (Double, Double) -> Int -> [Double]
+    linspace (lo, hi) n = map (\i -> (fromIntegral i)*(hi-lo)/fromIntegral (n-1)) [0..n-1]
+    sproduct :: [[Double]] -> [[Double]] ->  [[Double]]
+    sproduct gs [] = gs
+    sproduct gs (l:ls) =
+           let gs' = [x:g | g<-gs, x<-l]
+           in  sproduct gs' ls
+
+
+-- | 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
+        getNormals n >>= \case
+          (xi:etas) -> let
+              xi' = sigma_xi * xi
+              parCorr = xi' `scale` d
+              etas' = map (dist3 * sigma_eta *) etas
+              orthCorrs = zipWith scale etas' exs
+            in return (parCorr, orthCorrs)
+          _ -> error "Parameters too short"
+      else do -- 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)

@@ -1,64 +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))+{- |
+
+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))

@@ -1,31 +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-+{- |
+
+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
+

@@ -1,18 +1,21 @@-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+module Moo.GeneticAlgorithm.Multiobjective
+    (
+    -- * Types
+      SingleObjectiveProblem
+    , MultiObjectiveProblem
+    , MultiPhenotype
+    -- * Evaluation
+    , evalAllObjectives
+    , takeObjectiveValues
+    -- * NSGA-II: A non-dominated sorting genetic algorithm
+    , stepNSGA2
+    , stepNSGA2bt
+    , stepConstrainedNSGA2
+    , stepConstrainedNSGA2bt
+    -- * Performance metrics
+    , hypervolume
+    ) where
+
+import Moo.GeneticAlgorithm.Multiobjective.Types
+import Moo.GeneticAlgorithm.Multiobjective.NSGA2
+import Moo.GeneticAlgorithm.Multiobjective.Metrics

@@ -0,0 +1,143 @@+{-# LANGUAGE RankNTypes #-}
+{- | Performance metrics for multiobjective problems.
+
+-}
+
+module Moo.GeneticAlgorithm.Multiobjective.Metrics where
+
+
+import Data.List (tails, sortBy)
+import Data.Function (on)
+
+
+import Moo.GeneticAlgorithm.Types
+import Moo.GeneticAlgorithm.Multiobjective.Types
+import Moo.GeneticAlgorithm.Multiobjective.NSGA2
+
+
+type Point = [Double]
+
+
+-- | Calculate the hypervolume indicator using WFG algorithm.
+--
+-- Reference:
+-- While, L., Bradstreet, L., & Barone, L. (2012). A fast way of
+-- calculating exact hypervolumes. Evolutionary Computation, IEEE
+-- Transactions on, 16(1), 86-95.
+--
+hypervolume :: forall fn a . ObjectiveFunction fn a
+            => MultiObjectiveProblem fn   -- ^ multiobjective problem @mop@
+            -> [Objective]                -- ^ reference point (the worst point)
+            -> [MultiPhenotype a]         -- ^ a set of solutions to evaluate
+            -> Double                     -- ^ hypervolume
+hypervolume mop refPoint solutions =
+    let ptypes = map fst mop :: [ProblemType]
+        points = map takeObjectiveValues solutions
+    in  wfgHypervolume_sort 0 ptypes refPoint points
+
+
+-- | Basic (non-optimized) WFG algorithm to calculate hypervolume.
+--
+-- Reference: While et al. (2012).
+wfgHypervolume :: [ProblemType]  -- ^ problem types
+               -> Point          -- ^ reference point (the @worst@ point)
+               -> [Point]        -- ^ a set of points
+               -> Double
+wfgHypervolume ptypes worst pts =
+    let ptsAndTails = zip pts (drop 1 (tails pts)) :: [(Point, [Point])]
+        exclusiveHvs = map
+                       (\(pt, rest) -> exclusiveHypervolume ptypes worst pt rest)
+                       ptsAndTails
+    in  sum exclusiveHvs
+
+
+-- | WFG algorithm to calculate hypervolume with sorting optimization.
+wfgHypervolume_sort :: Int            -- ^ index of the objective to sort
+                    -> [ProblemType]  -- ^ problem types
+                    -> Point          -- ^ reference point (the @worst@ point)
+                    -> [Point]        -- ^ a set of points
+                    -> Double
+wfgHypervolume_sort k ptypes worst pts
+    | null ptypes || length ptypes <= k || k < 0 =
+        wfgHypervolume_sort 0 ptypes worst pts  -- bad input, sort the first objective
+    | otherwise =
+        let ptype = ptypes !! k
+            pts' = sortBy (flip compare `on` get ptype k) pts
+        in  wfgHypervolume ptypes worst pts'
+    where
+      get :: ProblemType -> Int -> [Double] -> Double
+      get Minimizing k objvals
+          | length objvals > k = objvals !! k
+          | otherwise          = inf
+      get Maximizing k objvals
+          | length objvals > k = objvals !! k
+          | otherwise          = - inf
+      inf :: Double
+      inf = 1/0
+
+
+-- | Construct a limited set (a step of the WFG algorithm).
+--
+-- @
+--     limitSet(S, p) = { limit(x, p) | x \in S }
+--     where limit(<s1, ..., sn>, <p1, ..., pn>) = < worse(s1,p1), ..., worse(sn, pn)>.
+-- @
+limitSet :: [ProblemType] -- ^ problem types
+         -> Point         -- ^ reference point
+         -> [Point]       -- ^ original set
+         -> [Point]       -- ^ limited set
+limitSet ptypes refPoint =
+    map (zipWith3 worst ptypes refPoint)
+  where
+    worst :: ProblemType -> Double -> Double -> Double
+    worst Minimizing x y | x > y     = x
+                         | otherwise = y
+    worst Maximizing x y | x < y     = x
+                         | otherwise = y
+
+
+-- | Construct a non-dominated subset (a step of the WFG algorithm).
+nondominatedSet :: [ProblemType]  -- ^ problem types
+                -> [Point]        -- ^ original set
+                -> [Point]        -- ^ a non-dominated subset
+nondominatedSet ptypes points =
+    let dominates = domination ptypes
+        dummySolutions = map (\objvals -> ([], objvals)) points :: [MultiPhenotype Double]
+        fronts = nondominatedSort dominates dummySolutions :: [[MultiPhenotype Double]]
+    in  case fronts of
+          (nds:_) -> map takeObjectiveValues nds
+          _       -> []
+
+
+-- | Calculate inclusive hypervolume of a point @p@ (the size of the
+-- part of the objective space dominated by @p@ alone).
+inclusiveHypervolume :: [ProblemType]  -- ^ problem types
+                     -> Point          -- ^ reference point (the @worst@ point)
+                     -> Point          -- ^ a point @p@ to evaluate
+                     -> Double         -- ^ inclusive hypervolume
+inclusiveHypervolume ptypes worst p =
+    product $ zipWith3 hyperside ptypes worst p
+ where
+    hyperside :: ProblemType -> Double -> Double -> Double
+    hyperside Minimizing upper x = pos $ upper - x
+    hyperside Maximizing lower x = pos $ x - lower
+    -- Positive part: to truncate the hypervolume if an unsuitable
+    -- reference point is given (not the worst one possible)
+    pos :: Double -> Double
+    pos x = 0.5 * (x + abs x)
+
+
+-- | Calculate exclusive hypervolume of a point @p@ relative to the
+-- @underlying@ set (the size of the part of the objective space that
+-- is dominated by @p@, but is not dominated by any member of the
+-- @underlying@ set).
+exclusiveHypervolume :: [ProblemType]  -- ^ problem types
+                     -> Point          -- ^ reference point (the @worst@ point)
+                     -> Point          -- ^ a point @p@ to evaluate
+                     -> [Point]        -- ^ an @underlying@ set of points
+                     -> Double         -- ^ exclusive hypervolume
+exclusiveHypervolume ptypes worst p underlying =
+    let inclusiveHv = inclusiveHypervolume ptypes worst p
+        nds = nondominatedSet ptypes $ limitSet ptypes p underlying
+        underlyingHv = wfgHypervolume ptypes worst nds
+    in  inclusiveHv - underlyingHv

@@ -1,495 +1,496 @@-{-# 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+{-# LANGUAGE Rank2Types, ConstraintKinds #-}
+{-# LANGUAGE FlexibleContexts #-}
+{- |
+
+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

@@ -1,45 +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+{-# 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

@@ -1,55 +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+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

@@ -1,111 +1,143 @@-{-# 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+{- | 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
+    , randomSampleIndices
+    , shuffle
+    -- * Building blocks
+    , withProbability
+    -- * Re-exports from random number generator packages
+    , getBool, getInt, getWord, getInt64, getWord64, getDouble
+    , runRand, evalRand, newPureMT, liftRand
+    , Rand, Random, PureMT
+    ) where
+
+import Control.Monad (liftM)
+import qualified Control.Monad.Random.Strict as MonadRandom
+import Control.Monad.Random.Strict (liftRand, runRand, evalRand)
+import Data.Complex (Complex (..))
+import Data.Int (Int64)
+import Data.Word (Word64)
+import System.Random (RandomGen, Random(..))
+import System.Random.Mersenne.Pure64
+import qualified System.Random.Shuffle as S
+import qualified Data.Set as Set
+
+type Rand = MonadRandom.Rand PureMT
+
+-- | 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 = liftRand $ \s -> randomR range s
+
+-- | Yield a new randomly selected value of type @a@.
+-- See 'System.Random.random' for details.
+getRandom :: Random a => Rand a
+getRandom = liftRand random
+
+getBool :: Rand Bool
+getBool = getRandom
+getDouble :: Rand Double
+getDouble = getRandom
+getWord :: Rand Word
+getWord = getRandom
+getInt :: Rand Int
+getInt = getRandom
+getInt64 :: Rand Int64
+getInt64 = getRandom
+getWord64 :: Rand Word64
+getWord64 = getRandom
+
+-- | 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 =
+  liftRand $ \g -> select g n (length xs) xs []
+  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)
+
+-- | Select @sampleSize@ numbers in the range from @0@ to @(populationSize-1)@.
+-- The function works best when @sampleSize@ is much smaller than @populationSize@.
+randomSampleIndices :: Int -> Int -> Rand [Int]
+randomSampleIndices sampleSize populationSize =
+    liftRand $ \g ->
+        let (sampleSet, g') = buildSampleSet g sampleSize Set.empty
+        in  (Set.toList sampleSet, g')
+  where
+    buildSampleSet g n s
+        | n <= 0 = (s, g)
+        | otherwise =
+            let (i, g') = randomR (0, populationSize-1) g
+            in  if (i `Set.member` s)
+                then buildSampleSet g' n s
+                else buildSampleSet g' (n-1) (Set.insert i s)
+
+-- | Randomly reorder the list.
+shuffle :: [a] -> Rand [a]
+shuffle xs = liftRand $ \g -> randomShuffle xs (length 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

@@ -1,252 +1,260 @@-{-# 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+{-# LANGUAGE BangPatterns, Rank2Types #-}
+{-# LANGUAGE GADTs #-}
+{- |
+
+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') = runRand initialize rng
+  return $ evalRand (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') = runRand 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.  Usually only the first two
+-- arguments are given, and the result is passed to 'runGA'.
+{-# 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.
+-- Usually only the first three arguments are given, and the result is
+-- passed to 'runGA'.
+{-# 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 via an IORef. Usually only the first three
+-- arguments are given, and the result is passed to 'runIO'.
+{-# 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') = runRand (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 LogHook a m w where
+    WriteEvery :: (Monad m, Monoid w)
+               => Int
+               -> (Int -> Population a -> w)
+               -> LogHook a m 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

@@ -1,158 +1,163 @@-{- |--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+{- |
+
+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)
+import qualified Data.Vector as V
+
+
+
+-- | 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 = do
+    let popvec = V.fromList xs
+    let popsize = V.length popvec
+    groups <- replicateM n $ randomSampleIndices size popsize
+    return $ map (tournament1 problem popvec) groups
+  where
+    tournament1 problem popvec group =
+      let contestants = map (popvec V.!) group
+          best = head $ bestFirst problem contestants
+      in  best
+
+-- | 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

@@ -1,76 +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]+{-# 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

@@ -1,30 +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+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

@@ -1,157 +1,158 @@-{-# 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)+{-# 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
+import Control.Parallel.Strategies (parMap, rseq)
+
+-- | 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 in parallel.
+instance (a1 ~ a2) =>
+    ObjectiveFunction (Genome a1 -> Objective) a2 where
+        evalObjective f gs = parMap rseq (\g -> (g, f g)) gs
+
+-- | 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)

@@ -1,81 +1,77 @@-{-# 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)+{- |
+
+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 (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 runRand 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 =
+    liftRand $ \rng -> randomGenomes rng n ranges
+
+
+-- | 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 parents'
+     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)

@@ -1,145 +1,176 @@-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+Moo
+===
+
+     ------------------------------------------------
+    < Moo. Breeding Genetic Algorithms with Haskell. >
+     ------------------------------------------------
+            \   ^__^
+             \  (oo)\_______
+                (__)\       )\/\
+                    ||----w |
+                    ||     ||
+
+
+Installation
+------------
+
+### Installation from Hackage
+
+Hackage is a Haskell community's package archive. This is where the
+latest versions of packages are published first.
+To install Moo from Hackage use Cabal-Install:
+
+  * install Haskell Platform or GHC and Cabal-Install
+  * run `cabal update`
+  * run `cabal install moo`
+
+
+### Installation with Stack
+
+Stackage is a stable package archive. Stackage builds are supposed to
+be reproducible. Stackage also provides Long Term Support releases.
+To build Moo with Stackage dependencies, use the `stack` tool:
+
+  * install [`stack`](https://docs.haskellstack.org/)
+  * if necessary, install GHC: run `stack setup`
+  * run: `stack update`
+  * in the project source directory run: `stack build`
+  * to run tests: `stack test`
+
+
+### Build Status
+
+[![Build Status](https://travis-ci.org/astanin/moo.svg?branch=master)](https://travis-ci.org/astanin/moo)
+
+
+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

@@ -1,2 +1,2 @@-import Distribution.Simple-main = defaultMain+import Distribution.Simple
+main = defaultMain

@@ -1,87 +1,81 @@-{-# 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)+{-# 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
+
+
+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 = return $ uniformGenomes (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)

@@ -1,84 +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]-    ]+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 evalRand (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 evalRand (pureMT 1) $ ctournament badvsugly
+        assertEqual "lesser degree of violation is preferred"
+                    (replicate n (-1.0)) $ (map (head . takeGenome) result)
+        let result = flip evalRand (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]
+    ]

@@ -1,35 +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]])-    ]+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 evalRand (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 evalRand (pureMT 1) $
+                             (nextGeneration Minimizing objective select 0 noCrossover noMutation)
+                             (Generations 1) (Left [[1],[2,2]])
+    ]

@@ -1,83 +1,87 @@-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-    ]+module Tests.Internals.TestCrossover where
+
+
+import Test.HUnit
+import System.Random.Mersenne.Pure64 (pureMT)
+import Data.List (group, transpose)
+
+
+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 evalRand (pureMT 1) $
+                      doNCrossovers 4 genomes (onePointCrossover 0.5)
+        let expected4 = [[1,0,0,0],[0,1,1,1],[1,1,0,0],[0,0,1,1]]
+        assertEqual "4 crossovers" expected4 result4
+        let genesums4 = map sum . transpose $ result4
+        assertEqual "gene-sums (4 genomes)" [2,2,2,2] genesums4
+        let result3 = flip evalRand (pureMT 1) $
+                      doNCrossovers 3 genomes (onePointCrossover 0.5)
+        let expected3 = [[1,0,0,0],[0,1,1,1],[1,1,0,0]]
+        assertEqual "3 crossovers" expected3 result3
+    , "do all crossovers" ~: do
+        let genomes = [[1,1,1,1],[0,0,0,0]] :: [[Int]]
+        let result = flip evalRand (pureMT 1) $
+                     doCrossovers genomes (onePointCrossover 0.5)
+        let expected = [[1,1,0,0],[0,0,1,1]]
+        assertEqual "all crossovers (2 genomes)" expected result
+        let genesums2 = map sum . transpose $ result
+        assertEqual "gene-sums (2 genomes)" [1,1,1,1] genesums2
+        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 evalRand (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 evalRand (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 evalRand (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 evalRand (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 evalRand (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 evalRand (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 evalRand (pureMT 1) mkChildren :: [Genome Double]
+        assertEqual "change approximately half" True $
+                    all (\s -> (s >= mu - 4*sigma && s <= mu + 4*sigma)) . map sum $
+                        children
+    ]

@@ -1,45 +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-    ]+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 evalRand (pureMT 1) $
+                    withProbability 0 (return . (+1)) 42
+        assertEqual "probability 1" 43 $
+                    flip evalRand (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 evalRand (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 evalRand (pureMT 1) $
+                         return . length . filter id =<< asymmetricMutate 1 0 g
+        assertEqual "flip all ones" 0 $
+                    flip evalRand (pureMT 1) $
+                         return . length . filter id =<< asymmetricMutate 0 1 g
+        assertEqual "flip all" 8 $
+                    flip evalRand (pureMT 1) $
+                         return . length . filter id =<< asymmetricMutate 1 1 g
+    ]

@@ -1,147 +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)-    ]+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 evalRand (pureMT 1) $
+                     loop (Generations 1)
+                     (stepNSGA2bt mp noCrossover noMutation) gs
+        assertEqual "solutions and ranking" (Set.fromList expected) (Set.fromList result)
+    ]

@@ -1,66 +1,67 @@-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-    ]+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 evalRand (pureMT 1) $
+                        tournamentSelect Minimizing 3 4 $
+                        map dummyGenome [3,2,4]
+        let resultMax = flip evalRand (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 evalRand (pureMT 1) $
+                         tournamentSelect Maximizing tsize times $
+                         genomes
+        let objVals = map takeObjectiveValue resultMany
+        -- take the same samples again with the same see
+        let samples = map (map (genomes !!)) $
+                      flip evalRand (pureMT 1) $
+                           replicateM times (randomSampleIndices tsize (length 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 evalRand (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 evalRand (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
+    ]

@@ -1,91 +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))-  ]+{- 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 400 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 dist ++ " >= " ++ 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))
+  ]

@@ -1,154 +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+-- | 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

@@ -1,35 +1,72 @@-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
+========
+
+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, logging hypervolume evolution in IO)
+
+  * [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)
+
+  * [fourbittrap.hs](fourbittrap.hs) concatenation of N-bit trap
+    functions is a difficult problem for genetic algorithms, and
+    requires to use large populations
+    (A basic GA with convergence check)
+
+Examples of integer-coded GAs:
+
+  * [ilp.hs](ilp.hs) an integer programming problem (a constrained GA
+    with genomes as lists of integers and a custom mutation operator)
+
+
+How to build examples within Cabal sandbox
+------------------------------------------
+
+For sandboxed builds, if you initialized the sandbox in the top-level
+directory of the moo source distribution as
+
+    cabal sandbox init
+
+and compiled the library with
+
+    cabal install
+
+then
+
+  * In the `examples/` directory run
+
+        cabal sandbox init --sandbox=../.cabal-sandbox
+
+  * Build examples like
+
+        cabal exec ghc -- --make example_file.hs
+
+    instead of
+
+        ghc --make example_file.hs

@@ -1,27 +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)+{- 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)

@@ -1,64 +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+{- 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

@@ -1,46 +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+{- 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

@@ -1,102 +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"+{-
+  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 evalRand 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"