packages feed

genprog 0.1 → 0.1.0.1

raw patch · 7 files changed

+957/−951 lines, 7 filesdep ~MonadRandomdep ~basedep ~syb

Dependency ranges changed: MonadRandom, base, syb, syz

Files

− GenProg.hs
@@ -1,721 +0,0 @@--- |--- Module      :  GenProg--- Copyright   :  (c) 2010 Jan Snajder--- License     :  BSD-3 (see the LICENSE file)------ Maintainer  :  Jan Snajder <jan.snajder@fer.hr>--- Stability   :  experimental--- Portability :  non-portable------ The Genetic Programming Library.------ /Genetic programming/ is an evolutionary optimization technique--- inspired by biological evolution. It is similar to /genetic algorithms/--- except that the individual solutions are programs (or, more generally, --- /expressions/) representing a solution to a given problem. A genetic --- program is represented as an /abstract syntax tree/ and associated --- with a custom-defined /fitness/ value indicating the quality of the --- solution. Starting from a randomly generated initial population of --- genetic programs, the genetic operators of /selection/, /crossover/, --- and (occasionally) /mutation/ are used to evolve programs of --- increasingly better quality.------ Standard reference is: John Koza. /Genetic programming:/--- /On the Programming of Computers by Means of Natural Selection/.--- MIT Press, 1992.------ In GenProg, a genetic program is represented by a value of an--- algebraic datatype. To use a datatype as a genetic program, it--- suffices to define it as an instance of the 'GenProg' typeclass.--- A custom datatype can be made an instance of the 'GenProg'--- typeclass, provided it is an instance of the 'Data' typeclass (see--- "GenProg.GenExpr.Data").------ An example of how to use this library is given below.-----------------------------------------------------------------------------------{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies,-    NoMonomorphismRestriction #-}--module GenProg (-  -- * Genetic programs-  GenProg (..),-  -- * Expressions-  generateFullExpr,-  generateGrownExpr,-  depth,-  nodes,-  -- * Individuals-  Ind,-  unInd,-  mkInd,-  aFitness,-  sFitness,-  -- * Population-  Pop,-  unPop,-  mkPop,-  generatePop,-  replenishPop,-  mergePop,-  best,-  avgFitness,-  avgDepth,-  avgNodes,-  -- * Genetic operators-  -- | The following functions are not meant to be used directly.-  -- They are exposed for debugging purposes.-  crossoverInd,-  mutateInd,-  crossoverPop,-  mutatePop,-  -- * Evolution state-  EvolState (..),-  -- * Control parameters-  Fitness,-  Mutate,-  defaultMutation,-  Terminate,-  tSuccess,-  tFitness,-  tGeneration,-  EvolParams (..),-  defaultEvolParams,-  -- * Evolution-  evolve,-  evolveFrom,-  evolveTrace,-  evolveTraceFrom-  -- * Example-  -- $Example-  ) where--import Data.List-import Data.Ord-import Data.Maybe-import Control.Monad-import Control.Monad.Random-import GenProg.GenExpr.Data---- | A typeclass defining a genetic program interface.  Datatypes @e@--- that are to be used as genetic programs must be instances of the--- 'GenExpr' typeclass and must implement this interface. -class (Eq e, GenExpr e, MonadRandom m) => GenProg m e | e -> m where-  -- | Generates a random terminal @T@.-  terminal :: m e-  -- | Generates a random nonterminal (functional) node @F(T,...,T)@ whose-  -- arguments are again terminals (this condition is not verified).-  nonterminal :: m e---------------------------------------------------------------------------------- Expressions---- | Generates a random expression of a given maximum depth.-generateExpr :: (GenProg m e) => m e -> Int -> m e-generateExpr g d-  | d < 1     = error "GenProg.generateExpr: Invalid expression depth"-  | otherwise = nonterminal >>= step (d - 1)-  where step 0 _ = terminal-        step d e = nodeMapM (const g >=> step (d - 1)) e---- | Generates a random expression fully expanded to the specified depth.-generateFullExpr :: (GenProg m e) => Int -> m e-generateFullExpr = generateExpr nonterminal---- | Generates a random expression of limited depth. The maximum depth of--- the resulting expression may be less than the specified depth--- limit, and paths may be of different length.-generateGrownExpr :: (GenProg m e) => Int -> m e-generateGrownExpr d = do-  t <- getRandom-  generateExpr (if t then terminal else nonterminal) d---------------------------------------------------------------------------------- Individuals---- | A genetically programmed individual, representing a basic unit--- of evolution. (Basically a wrapper around a genetically programmable--- expression.)-data Ind e = Ind {-  -- | Returns the expression wrapped by an individual.-  unInd :: e,-  -- | Adjusted fitness of an individual. Adjusted fitness equals-  -- @1/(1+s)@, where @s@ is the standardized fitness as computed by-  -- 'fitness'. To reduce computational costs, this value is computed-  -- only once and then cached.-  aFitness :: Double,-  -- The indices of inner (functional) nodes of an individual's expression.-  iNodes :: [Int],-  -- The indices of external (terminal) nodes of an individual's expression.-  eNodes :: [Int] }-  deriving (Show)--instance (Eq e) => Eq (Ind e) where-  i1 == i2 = unInd i1 == unInd i2--instance (Eq e) => Ord (Ind e) where-  compare = comparing aFitness---- | Wraps an expression into an individual.-mkInd :: (GenProg m e) => Fitness e -> e -> Ind e-mkInd f e = Ind e (adjust $ f e) fs ts-  where (fs,ts) = nodeIndices e---- Adjusts fitness.-adjust :: Double -> Double-adjust f = 1 / (1 + max 0 f)---- Unadjusts fitness (the inverse of adjustFitness).-unadjust :: Double -> Double-unadjust f = 1 / f - 1---- | Standardized fitness of an individual as computed by 'fitness'-sFitness :: Ind e -> Double-sFitness = unadjust . aFitness---------------------------------------------------------------------------------- Population---- | A population of individuals. (Basically a wrapper around a list of--- individuals.)-data Pop e = Pop-  { unPop  :: [Ind e]   -- ^ Unwraps a population.-  , dist_  :: [Double]  -- ^ Fitness distribution.-  } deriving (Show, Eq)---- | Wraps a list of individuals into a population.-mkPop :: [Ind e] -> Pop e-mkPop is = Pop is ds-  where ds = map snd . distribution $-             map (\i -> (unInd i, aFitness i)) is---- | Generate population of given size and given depth limit using--- /ramped half-and-half/ method (Koza, 1992): for each depth value from 0 to--- the initial depth limit 'iDepth', 50% of individuals are generated using--- 'generateFullExpr' and 50% are generated using--- 'generateGrownExpr'. Afterwards, duplicates are removed, thus the--- size of the resulting population may actually be less than the--- specified size.-generatePop :: (GenProg m e) => EvolParams m e -> m (Pop e)-generatePop p-  | s < 2 || n==0 = error "GenProg.generatePop: Invalid population size"-  | otherwise = do-    iss <- forM [2..di] $ \i -> do-      is1 <- replicateM n (mkInd (fitness p) `liftM` generateFullExpr di)-      is2 <- replicateM n (mkInd (fitness p) `liftM` generateGrownExpr di)-      return $ is1 ++ is2-    return . mkPop . nub $ concat iss-  where n  = s `div` (2 * (di - 1))-        s  = popSize p-        di = iDepth p---- | Replenishes a population up to 'popSize' by randomly--- generating new individuals.-replenishPop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Pop e)-replenishPop p pop1 = do-  pop2 <- generatePop p-  return . mkPop $ unPop pop1 ++ drop s (unPop pop2)-  where s = length $ unPop pop1---- | Merges two populations by taking 'popSize' best-fitted individuals--- from the union of the two populations.-mergePop :: (GenProg m e) => EvolParams m e -> Pop e -> Pop e -> Pop e-mergePop p pop1 pop2 = mkPop $ take (popSize p) is-  where is = sortBy (flip $ comparing aFitness) $ unPop pop1 ++ unPop pop2---- | Population's best-fitted individual.-best :: Pop e -> Ind e-best = maximumBy (comparing aFitness) . unPop--avg :: (Fractional a) => [a] -> a-avg xs = sum xs / realToFrac n-  where n = length xs---- | Population's average standardized fitness.-avgFitness :: Pop e -> Double-avgFitness = avg . map (unadjust . aFitness) . unPop---- | Average depth of expressions in the population.-avgDepth :: (GenProg m e) => Pop e -> Double-avgDepth = avg . map (realToFrac . depth . unInd) . unPop---- | Average number of expression nodes in the population.-avgNodes :: (GenProg m e) => Pop e -> Double-avgNodes = avg . map (realToFrac . nodes . unInd) . unPop---------------------------------------------------------------------------------- Genetic operators---- Selects at random an index of an expression node. Functional--- (internal) nodes are selected with probability 'pci', whereas--- terminal (external) nodes are selecred with probability '1-pi'.-selectNode :: (GenProg m e, MonadRandom m) => Double -> Ind e -> m Int-selectNode pi i-  | null $ iNodes i = oneof $ eNodes i-  | otherwise       = choice pi (oneof $ iNodes i) (oneof $ eNodes i)---- | Crossover operation of two individuals, resulting in two--- offsprings. Crossover is performed by choosing at random two nodes--- in each expressions, and then by exchanging the subexpressions--- rooted at these nodes between the two individuals. The probability--- that an internal (functional) node is chosen as crossover point is--- set by the 'ciProb' parameter in 'EvolParams', whereas the--- probability that an external (terminal) node is chosen equals--- @1-ciProb@. Among internal and external nodes, nodes are chosen--- uniformly at random. If the depth of a created offspring exceeds--- the depth limit 'cDepth' specified by evolution parameters--- 'EvolParams', that offspring is discarded and a parent is--- reproduced (i.e., copied as-is).-crossoverInd :: (GenProg m e) =>-  EvolParams m e -> Ind e -> Ind e -> m (Ind e, Ind e)-crossoverInd p i1 i2 = do-  n1 <- selectNode (ciProb p) i1-  n2 <- selectNode (ciProb p) i2-  let (r1,r2) = exchange (unInd i1) n1 (unInd i2) n2-  return (if depth r1 <= cDepth p then mkInd (fitness p) r1 else i1,-          if depth r2 <= cDepth p then mkInd (fitness p) r2 else i2)---- | Mutates an individual by applying the mutation function @mutate@--- to a randomly selected node. The probability that an internal--- (functional) node is chosen for muration is set by the 'miProb'--- parameter in 'EvolParams', whereas the probability that an external--- (terminal) node is chosen equals @1-miProb@. Among internal and--- external nodes, nodes are chosen uniformly at random. If the depth--- of the mutated expression exceeds the depth limit 'cDepth'--- specified by evolution parameters 'EvolParams', the individual is--- left unaltered.-mutateInd :: (GenProg m e) => EvolParams m e -> Ind e -> m (Ind e)-mutateInd p i = do-  n  <- selectNode (miProb p) i-  e2 <- adjustM (mutate p) e1 n-  return . mkInd (fitness p) $ if depth e2 <= cDepth p then e2 else e1-  where e1 = unInd i---- Discrete distribution.-type Distribution a = [(a, Double)]---- Computes distribution from a weighted list.--- The weights need not sum to 1.-distribution :: [(a, Double)] -> Distribution a-distribution xs = [(x,f i) | ((x,_),i) <- zip xs [1..]]-  where f i = sum . map snd $ take i ys-        s   = sum $ map snd xs-        ys  = map (\(x, w) -> (x, w/s)) xs---- Samples a value from a discrete distribution.-choose :: (MonadRandom m) => Distribution a -> m a-choose xs = do-  p <- getRandomR (0,1)-  return . fst . fromJust $ find ((>= p) . snd) xs---- Chose first action with probability 'p' and second with probability--- 1-p.-choice :: (MonadRandom m) => Double -> m a -> m a -> m a-choice p a1 a2 = do-  r <- getRandomR (0,1)-  if r <= p then a1 else a2--oneof :: (MonadRandom m) => [a] -> m a-oneof xs = (xs!!) `liftM` getRandomR (0,length xs-1)---- Fitness-proportionate selection of an individual from a population.-selectInd :: (MonadRandom m) => Pop e -> m (Ind e)-selectInd pop = choose (zip (unPop pop) (dist_ pop))--reproducePop :: (MonadRandom m) => Pop e -> m (Ind e)-reproducePop = selectInd---- | Applies crossover to two randomly chosen individuals from a--- population. The probability of an individual being chosen as parent--- is fitness-proportionate (individuals with better fitness have--- better chanches of being chosen for crossover).-crossoverPop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Ind e,Ind e)-crossoverPop p pop = do-  i1 <- selectInd pop-  i2 <- selectInd pop-  crossoverInd p i1 i2---- | Applies mutation operation to individuals from a population. The--- probability of mutating each individual is determined by 'mProb' parameter--- from 'EvalParams'.-mutatePop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Pop e)-mutatePop p pop-  | mProb p == 0 = return pop-  | otherwise    = liftM mkPop . forM (unPop pop) $ \i ->-                     choice (mProb p) (mutateInd p i) (return i)---------------------------------------------------------------------------------- Evolution state---- | The state of the evolution.-data EvolState e = EvolState-  { pop        :: Pop e    -- ^ Current population.-  , iter       :: Int      -- ^ Iteration (current generation number).-  , cachedBest :: Ind e    -- ^ Best individual evolved so far.-  } deriving (Show,Eq)--initState :: Pop e -> EvolState e-initState pop =-  EvolState { pop = pop, iter = 0, cachedBest = best pop }---- | Advances to next evolution state.-nextState :: (GenProg m e ) =>-  EvolParams m e -> EvolState e -> m (EvolState e)-nextState p es1 = do-  pop2 <- evolvePop p pop1-  return $ es1 { pop = pop2, iter = iter es1 + 1,-                 cachedBest = max (cachedBest es1) (best pop1) }-  where pop1 = pop es1---------------------------------------------------------------------------------- Control parameters---- | Standardized fitness. It takes on values from 0 (best fitness) to--- +infinity (worst fitness).-type Fitness e = e -> Double---- | A function to mutate a chosen expression node.-type Mutate m e = e -> m e---- | Default mutation. Replaces a node, irrespective of its value,--- with a randomly generated subexpression whose depth is limited to--- 'iDepth'.-defaultMutation :: (GenProg m e) => EvolParams m e -> Mutate m e-defaultMutation p = const $ generateGrownExpr (iDepth p)---- | Termination predicate.-type Terminate e = EvolState e -> Bool---- | Termination predicate: terminate if any individual satisfies the--- specified predicate.-tSuccess :: (e -> Bool) -> Terminate e-tSuccess c = any (c . unInd) . unPop . pop---- | Termination predicate: terminate if best individual's--- standardized fitness is greater than or equal to the specified value.-tFitness :: (GenProg m e) => Double -> Terminate e-tFitness f = (>= f) . unadjust . aFitness . cachedBest---- | Termination predicate: terminate after running for the specified--- number of iterations.-tGeneration :: Int -> Terminate e-tGeneration n = (>=n) . iter---- | Parameters governing the evolution.------ Default evolution parameters,--- as used in (Koza, 1992), are defined by 'defaultEvolParams'--- and indicated below. At least the fitness function 'fitness' should--- be overriden.-data EvolParams m e = EvolParams {-  -- | Population size (number of individuals). Default is @500@.-  popSize   :: Int,-  -- | Depth of expressions in initial population. Default is @6@.-  iDepth    :: Int,-  -- | Maximum depth of expressions created during the evolution.-  -- Default is @17@.-  cDepth    :: Int,-  -- | Probability of crossover. Default is @0.9@. If crossover is not-  -- chosen, an individual is simply reproduced (copied as-is) into-  -- the next generation.-  cProb     :: Double,-  -- | Probability that an internal (functional) node is chosen as a-  -- crossover point. Default is @0.9@. If an internal node is not-  -- chosen, an external (terminal) node is-  -- chosen.-  ciProb    :: Double,-  -- | Probability that an individual gets mutated. Default is @0@-  -- (no mutation).-  mProb     :: Double,-  -- | Probability that an internal (functional) node is chosen for-  -- mutation. Default is @0.1@.-  miProb    :: Double,-  -- | Standardized fitness function. Default value is @undefined@-  -- (must be overriden).-  fitness   :: Fitness e,-  -- | Mutation function. Defines how to change a randomly chosen-  -- node. Default is @defaultMutation defaultEvolParams@-  -- (replacement of a chosen node with a randomly generated subexpression).-  mutate    :: Mutate m e,-  -- | Elitist factor: number of best-fitted individuals that are preserved-  -- from each generation (reproduced as-is into next evolution state).-  -- Default is @0@.-  elitists  :: Int,-  -- | Termination predicate. Default is @50@ (terminate after 50 generations).-  terminate :: Terminate e }--defaultEvolParams = EvolParams-  { popSize   = 500-  , iDepth    = 6-  , cDepth    = 17-  , cProb     = 0.9-  , ciProb    = 0.9-  , mProb     = 0.0-  , miProb    = 0.1-  , terminate = tGeneration 50-  , fitness   = error "GenProg.defaultEvolParams: fitness function is undefined"-  , mutate    = const $ generateGrownExpr (iDepth defaultEvolParams)-  , elitists  = 0 }---------------------------------------------------------------------------------- Evolution--untilM :: (Monad m) => (a -> Bool) -> (a -> m a) -> a -> m a-untilM p f x | p x       = return x-             | otherwise = f x >>= untilM p f--iterateUntilM :: (Monad m) => (a -> Bool) -> (a -> m a) -> a -> m [a]-iterateUntilM p f x-  | p x       = return []-  | otherwise = do y  <- f x-                   ys <- iterateUntilM p f y-                   return (y:ys)---- | Evolves one population from another one by performing a single--- evolution step.-evolvePop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Pop e)-evolvePop p pop1 = do-     pop2 <- mkPop `liftM` untilM ((>= s) . length) step []-     pop3 <- mutatePop p pop2-     return $ mkPop (elite ++ unPop pop3)-  where s = popSize p - length elite-        elite = take (elitists p) topRanked-        topRanked = sortBy (flip $ comparing aFitness) $ unPop pop1-        step is | length is == s - 1 = (:is) `liftM` reproducePop pop1-                | otherwise = choice (cProb p)-                    (do (i1,i2) <- crossoverPop p pop1; return (i1:i2:is))-                    ((:is) `liftM` reproducePop pop1)---- | Creates an initial population and evolves it until termination--- predicate is satisfied, returning the last evolution state.-evolve :: (GenProg m e) => EvolParams m e -> m (EvolState e)-evolve p = -- generatePop p >>= evolveFrom p-  last `liftM` evolveTrace p---- | Evolves a given initial population until termination--- predicate is satisfied, returning the last evolution state.--- If the size of the initial population is less than--- 'popSize', the population will be replenished (see 'replenishPop').-evolveFrom :: (GenProg m e) => EvolParams m e -> Pop e -> m (EvolState e)-evolveFrom p pop = -- untilM (terminate p) (nextState p) . initState-  last `liftM` evolveTraceFrom p pop---- | Runs evolution on a given initial population until termination--- predicate is satisfied and returns a list of successive evolution--- states. If the size of the initial population is less than--- 'popSize', the population will be replenished (see 'replenishPop').-evolveTraceFrom :: (GenProg m e) => EvolParams m e -> Pop e -> m [EvolState e]-evolveTraceFrom p pop1 =-  iterateUntilM (terminate p) (nextState p) . initState =<< replenishPop p pop1---- | Creates an initial population and runs evolution until--- termination predicate is satisfied. Returns a list of successive--- evolution states.-evolveTrace :: (GenProg m e) => EvolParams m e -> m [EvolState e]-evolveTrace p = generatePop p >>= evolveTraceFrom p---------------------------------------------------------------------------------- Example--{- $Example--This is a simple, worked through example of how to use the GenProg-library. Given a target number @n@, out aim is to evolve an arithmetic-expression that evaluates to @n@. For example, given @13@ as the-target number, one possible solution is @(3 * 5) - 2@. The constants-allowed to appear in the expression are restricted to integers from 1-to 9. The allowed operations are @+@, @-@, @*@, and integer division-without remainder.--We begin by defining the datatype for the genetically programed-expression:--@--- The following language extensions need to be enabled:--- DeriveDataTypeable, FlexibleInstances, MultiParamTypeClasses--import GenProg-import Data.Generics-import Control.Monad-import Control.Monad.Random--data E = Plus E E-       | Minus E E-       | Times E E-       | Div E E-       | Const Int-       deriving (Typeable,Data,Eq,Show)-@--In order to evolve arithmetic expressions, we need to be able to-compute their values. To this end we define--@-eval :: E -> Maybe Int-eval (Const c)     = Just c-eval (Plus e1 e2)  = liftM2 (+) (eval e1) (eval e2)-eval (Minus e1 e2) = liftM2 (-) (eval e1) (eval e2)-eval (Times e1 e2) = liftM2 (*) (eval e1) (eval e2)-eval (Div e1 e2) | ok        = liftM2 div x1 x2-                 | otherwise = Nothing-  where (x1,x2) = (eval e1,eval e2)-        ok = x2 /= Just 0 && liftM2 mod x1 x2 == Just 0-@--Dividing by zero and dividing with a remainder are not allowed and in-such cases we return @Nothing@.--Because we have made @E@ an instance of the 'Data' typeclass, it can-be readily used as a genetically programmable expression. Next step is-to make 'E' an instance of the 'GenProg' typeclass:--@-instance GenProg (Rand StdGen) E where-  terminal    = Const `liftM` getRandomR (1,9)-  nonterminal = do-    r <- getRandomR (0,3)-    [liftM2 Plus terminal terminal,-     liftM2 Minus terminal terminal,-     liftM2 Times terminal terminal,-     liftM2 Div terminal terminal] !! r-@--Thus, a random terminal node contains one of the constants from 1 to-9. A nonterminal node can be one of the four arithmetic operations,-each with terminal nodes as arguments.  Note that computations are run-within the standard random generator monad (@Rand StdGen@).--The fitness function evaluates the accurateness of the arithmetic-expression with respect to the target number. If the value of the-expression is far off from the target number @n@, the standardized-fitness should be high. Moreover, we would like to keep the expression-as simple as possible. To this end, we include a /parsimony factor/-that is proportional to the number of nodes an expression has. We-define the overall standardized fitness as--@-myFitness :: Double => Int -> E -> Double-myFitness n e = error + size-  where error = realToFrac $ maybe maxBound (abs . (n-)) (eval e)-        size  = (realToFrac $ nodes e) / 100-@--The number of nodes is divided by a factor of 100 to make it less-important than the numeric accuracy of the expression.--We now have everything in place to get the evolution going. We will use-default evolution parameters and choose @12345@ as the target number:-->>> let params = defaultEvolParams { fitness = myFitness 12345 }--Let us first create a random number generator: -->>> let g = mkStdGen 0--We are doing this because we want our results to be reproducible, and-because we want to be able to compare the results of different-evolution runs. Normally, you would use @getStdGen@ to get a random-generator with random seed.--To run the evolution and get the best evolved individual, we type-->>> let i = cachedBest $ evalRand (evolve params) g--To check out its standardized fitness, we type-->>> sFitness i-39.61--Let us see how the actual expression looks like:-->>> unInd i-Times (Minus (Minus (Minus (Plus (Const 4) (Const 4)) (Plus (Const 6) -(Const 7))) (Minus (Minus (Const 5) (Const 9)) (Plus (Minus (Const 5) -(Const 9)) (Minus (Const 4) (Const 4))))) (Plus (Times (Plus (Const 5) -(Const 1)) (Const 6)) (Times (Plus (Const 9) (Const 3)) (Minus (Const 1) -(Const 8))))) (Div (Times (Plus (Plus (Const 3) (Const 5)) (Times (Const 4) -(Const 7))) (Plus (Const 4) (Const 4))) (Minus (Minus (Plus (Const 2) -(Const 8)) (Plus (Const 6) (Const 7))) (Plus (Minus (Const 5) (Const 9)) -(Minus (Const 4) (Const 4)))))--The number of nodes is-->>> nodes $ unInd i-61--Let us see to what number the expression evaluates:-->>> eval $ unInd i-Just 12384--So in this run we didn't get a perfect match, but we were close. Let-us see if we can do better.--When doing genetic programming, it is always a good idea to experiment-a bit with the parameters. There are no parameters that work best for-any given problem. You can learn a lot about how parameters influence-the evolution by analysing how the evolution progresses in time. This-can be accomplised by evolving an evolution trace:-->>> let trace = evalRand (evolveTrace params) g--We can now analyse how the standardized fitness of the-best individual improves during the evolution:-->>> map (sFitness . best . pop) trace-[9591.35,2343.59,1935.59,2343.59,903.51,903.45,585.59,585.59,327.45,225.41,-225.41,135.43,57.49,39.61,39.61,39.61,39.61,39.61,57.43,57.47,57.43,57.45,-57.33,57.43,57.43,57.45,57.43,57.43,57.35,57.35,57.43,57.27,57.33,57.33,57.43,-57.29,57.33,57.41,57.29,57.43,57.33,57.35,57.35,57.33,57.39,57.39,57.39,57.33,-57.37,57.37]--We see that at some point the fitness decreases and then increases-again. This indicates that the best fitted individual was lost by-evolving from one generation to the other. We can prevent this by-employing the /elitist strategy/. Let us see what happens if we-preserve a best fitted individual in each generation:-->>> let trace = evalRand (evolveTrace params {elitists = 1}) g ->>> map (sFitness . best . pop) trace-[9591.35,2343.59,711.61,711.61,711.61,711.61,57.55,57.53,57.39,57.39,57.39,-57.39,57.37,57.37,57.37,57.37,57.37,57.37,57.37,57.37,57.35,57.35,57.35,-57.35,57.35,57.35,57.35,57.35,57.35,57.35,57.33,57.33,57.33,57.33,57.33,-57.33,57.33,57.33,57.33,25.31,25.31,25.31,25.31,25.31,25.31,25.296,25.296,-25.296,25.296,25.296]--This gives us better fitness, but still not an exact match:-->>> let i = cachedBest $ last trace->>> eval $ unInd i-Just 12320--In the previous evolution run fitness converged relatively fast, but then-remained stuck. To stir up things a little, let us allow for some-mutation. Setting mutation probability to 5%, while retaining the-elitist strategy, we get-->>> let trace = evalRand (evolveTrace params {elitists = 1, mProb = 0.05}) g->>> map (sFitness . best . pop) trace-[9591.35,9591.35,9591.35,9591.35,9591.35,9591.35,9159.35,8403.23,7239.11,-6087.15,6087.15,1479.13,819.21,60.13,51.19,5.19,5.19,5.19,5.19,5.19,1.23,-1.23,1.23,1.23,1.23,1.23,1.21,1.21,1.21,1.21,0.23998,0.23998,0.23998,0.23998,-0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,-0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998]--This time we've got a perfect match:-->>> let i = cachedBest $ last trace->>> eval $ unInd i-Just 12345--while at the same time the expression is rather compact:-->>> unInd i-Plus (Times (Const 4) (Plus (Const 9) (Const 4))) (Plus (Plus (Times -(Plus (Const 4) (Const 3)) (Times (Times (Const 3) (Const 9)) (Times -(Const 5) (Plus (Const 9) (Const 4))))) (Const 3)) (Const 5))->>> nodes $ unInd i-23---}
− GenProg/GenExpr.hs
@@ -1,59 +0,0 @@--- |--- Module      :  GenProg.GenExpr--- Copyright   :  (c) 2010 Jan Snajder--- License     :  BSD-3 (see the LICENSE file)------ Maintainer  :  Jan Snajder <jan.snajder@fer.hr>--- Stability   :  experimental--- Portability :  non-portable------ An interface to genetically programmable expressions.-----------------------------------------------------------------------------------module GenProg.GenExpr (-  GenExpr (..)) where--import Control.Monad---- | This typeclass defines an interface to expressions--- that can be genetically programmed.  The operations that must be--- provided by instances of this class are used for the generation--- of random individuals as well as crossover and mutation operations.--- (An instance for members of the @Data@ typeclass is provided in--- "GenProg.GenExpr.Data".)------ Minimal complete definition: 'exchange', 'nodeMapM', 'nodeMapQ',--- and 'nodeIndices'.-class GenExpr e where-  -- | Exchanges subtrees of two expressions:-  -- @exchange e1 n1 e2 n2@ replaces the subexpression of @e1@ rooted in node-  -- @n1@ with the subexpression of @e2@ rooted in @n2@, and vice versa.-  exchange :: e -> Int -> e -> Int -> (e, e)-  -- | Maps a monadic transformation function over the immediate-  -- children of the given node.-  nodeMapM :: Monad m => (e -> m e) -> e -> m e-  -- | Maps a query function over the immediate children of the given-  -- node and returns a list of results.-  nodeMapQ :: (e -> a) -> e -> [a]-  -- | A list of indices of internal (functional) and external-  -- (terminal) nodes of an expression.-  nodeIndices :: e -> ([Int], [Int])-  -- | Adjusts a subexpression rooted at the given node by applying a-  -- monadic transformation function.-  adjustM :: (Monad m) => (e -> m e) -> e -> Int -> m e-  -- | Number of nodes an expression has.-  nodes :: e -> Int-  -- | The depth of an expression. Equals 1 for single-node expressions.-  depth :: e -> Int---  -- | Default method (expensive because it calls exchange twice).-  adjustM f e n = replace e n `liftM` f (get e n)-    where get e n = fst $ exchange e 0 e n-          replace e1 n1 e2 = fst $ exchange e1 n1 e2 0--  nodes = (+1) . foldr (+) 0 . nodeMapQ nodes --  depth = (+1) . foldr max 0 . nodeMapQ depth- 
− GenProg/GenExpr/Data.hs
@@ -1,153 +0,0 @@--- |--- Module      :  GenProg.GenExpr.Data--- Copyright   :  (c) 2010 Jan Snajder--- License     :  BSD-3 (see the LICENSE file)------ Maintainer  :  Jan Snajder <jan.snajder@fer.hr>--- Stability   :  experimental--- Portability :  non-portable------ Implementation of the @GenProg.GenExpr@ interface for members of--- the 'Data' typeclass. The implementation is based on SYB and SYZ--- generic programming frameworks (see--- <http://hackage.haskell.org/package/syb> and--- <http://hackage.haskell.org/package/syz> for details).------ NB: Subexpressions that are candidates for crossover points or--- mutation must be of the same type as the expression itself, and--- must be reachable from the root node by type-preserving traversal.--- See below for an example.-----------------------------------------------------------------------------------{-# LANGUAGE ScopedTypeVariables, FlexibleInstances, Rank2Types,-    UndecidableInstances, DeriveDataTypeable #-}--module GenProg.GenExpr.Data (-  -- | This module re-exports @GenExpr@ typeclass.-  GenExpr (..)-  -- * Example-  -- $Example-  ) where--import Data.Generics-import Data.Generics.Zipper-import Data.Maybe-import Control.Monad-import GenProg.GenExpr--moduleName = "GenProg.GenExpr.Data"--instance (Data a) => GenExpr a where--  -- | Exchanges two expression nodes. Works by using two generic-  -- zippers and exchanging their holes.-  exchange e1 n1 e2 n2 = (fromZipper y1, fromZipper y2)-    where z1 = typeMoveForUnsafe n1 $ toZipper e1-          z2 = typeMoveForUnsafe n2 $ toZipper e2-          (y1,y2) = exchangeHoles z1 z2--  -- | Adjust an expression node. Works by applying a monadic-  -- tranformation on a zipper hole.-  adjustM f e n = fromZipper `liftM` transM (mkM f) z-    where z = typeMoveForUnsafe n (toZipper e)--  nodeMapM f = gmapM (mkM f)--  nodeMapQ q (x::a) = concat $ gmapQ ([] `mkQ` (\(y::a) -> [q y])) x--  nodeIndices = index 0 [] [] . toZipper---- Zipper moves--type Move a = Zipper a -> Maybe (Zipper a)--backtrack :: (Typeable a) => Move a-backtrack z = do-  z2 <- up z-  right z2 `mplus` backtrack z2--repeatM :: (Monad m) => Int -> (a -> m a) -> a -> m a-repeatM 0 _ x = return x-repeatM n f x = f x >>= repeatM (n - 1) f---- Moves zipper to next node in DFS order, but does not move down the--- zipper if node satisfies query 'q'.-nextDfsQ :: Typeable a => GenericQ Bool -> Move a-nextDfsQ q z = (if query q z then Nothing else down' z)-  `mplus` right z `mplus` backtrack z---- Moves the zipper to node 'n' from current position in DFS order,--- skipping nodes not satisfying query 'q2' and descending only down--- the nodes satisfying query 'q1'.-moveForQ :: (Typeable a) => GenericQ Bool -> GenericQ Bool -> Int -> Move a-moveForQ _  _  0 z = Just z-moveForQ q1 q2 n z = do-  z2 <- nextDfsQ q1 z-  moveForQ q1 q2 (if query q2 z2 then n - 1 else n) z2---- Moves the zipper to node 'n' from current position in DFS order,--- counting only nodes of type 'a', and not descending down the nodes--- of other type.-typeMoveFor :: (Typeable a) => Int -> Move a-typeMoveFor n (z::Zipper a) =-  moveForQ (True `mkQ` (\(_::a) -> False)) (False `mkQ` (\(_::a) -> True)) n z---- | Same as typeMoveFor, but throws an error if node index is out of--- bound.-typeMoveForUnsafe :: (Typeable a) => Int -> Zipper a -> Zipper a-typeMoveForUnsafe n z = fromMaybe-  (error $ moduleName ++ ".typeMoveForUnsafe: Nonexisting node.")-  (typeMoveFor n z)---- | Exchanges two zipper holes.-exchangeHoles :: (Data a) => Zipper a -> Zipper a -> (Zipper a, Zipper a)-exchangeHoles (z1::Zipper a) (z2::Zipper a) = (y1,y2)-  where Just h1 = getHole z1 :: Maybe a-        Just h2 = getHole z2 :: Maybe a-        y1 = setHole h2 z1-        y2 = setHole h1 z2--index :: (Data a) => Int -> [Int] -> [Int] -> Zipper a -> ([Int], [Int])-index i is es (z :: Zipper a) =-  maybe (is2,es2) (index (i + 1) is2 es2) (typeMoveFor 1 z)-  where Just h = getHole z :: Maybe a-        (is2,es2) = if terminalQ h then (is,i:es) else (i:is,es)--terminalQ :: (Data a) => a -> Bool-terminalQ = null . nodeMapQ id--{- $Example--Suppose you have a datatype defined as--@-data E = A E E-       | B String [E]-       | C- deriving (Eq,Show,Typeable,Data)-@--and an expression defined as--@-e = A (A C C) (B \"abc\" [C,C])-@--The subexpressions of a @e@ are considered to be only the subvalues of-@e@ that are of the same type as @e@.  Thus, the number of nodes of-expression @e@ is-->>> nodes e-5- -because subvalues of node @B@ are of different type than expression-@e@ and therefore not considered as subexpressions. --Consequently, during a genetic programming run, subexpressions that-are of a different type than the expression itself, or subexpression-that cannot be reached from the root node by a type-preserving-traversal, cannot be chosen as crossover points nor can they be-mutated.---}
genprog.cabal view
@@ -1,16 +1,17 @@-Name:                genprog-Version:             0.1-Synopsis:            Genetic programming library-License:             BSD3-License-file:        LICENSE-Author:              Jan Snajder-Maintainer:          jan.snajder@fer.hr-Copyright:           (c) 2010 Jan Snajder-Category:            AI, Algorithms, Optimisation-Stability:           Experimental-Build-type:          Simple-Cabal-version:       >= 1.6-Description:+name:                genprog+version:             0.1.0.1+synopsis:            Genetic programming library+license:             BSD3+license-file:        LICENSE+author:              Jan Snajder+maintainer:          jan.snajder@fer.hr+homepage:            http://github.com/jsnajder/genprog+copyright:           (c) 2010 Jan Snajder+category:            AI, Algorithms, Optimisation+stability:           Experimental+build-type:          Simple+cabal-version:       >= 1.8+description:   This package provides a /genetic programming/ framework. Genetic   programming is an evolutionary technique, inspired by biological   evolution, to evolve programs for solving specific problems. A genetic@@ -21,12 +22,17 @@   and /mutation/ are used to evolve programs of increasingly better   quality. -Library-  Exposed-modules: +library+  exposed-modules:      GenProg, GenProg.GenExpr, GenProg.GenExpr.Data-  Build-depends: -    base == 4.*, syb >= 0.1.0.2, syz >= 0.2, MonadRandom >= 0.1.5-  Extensions:+  build-depends: +    base == 4.5.*, syb == 0.4.*, syz == 0.2.*, MonadRandom == 0.1.*+  hs-source-dirs: src+  extensions:     MultiParamTypeClasses, FunctionalDependencies,      NoMonomorphismRestriction, ScopedTypeVariables, FlexibleInstances,      Rank2Types, UndecidableInstances, DeriveDataTypeable++source-repository head+  type:     git+  location: https://github.com/jsnajder/genprog
+ src/GenProg.hs view
@@ -0,0 +1,721 @@+-- |+-- Module      :  GenProg+-- Copyright   :  (c) 2010 Jan Snajder+-- License     :  BSD-3 (see the LICENSE file)+--+-- Maintainer  :  Jan Snajder <jan.snajder@fer.hr>+-- Stability   :  experimental+-- Portability :  non-portable+--+-- The Genetic Programming Library.+--+-- /Genetic programming/ is an evolutionary optimization technique+-- inspired by biological evolution. It is similar to /genetic algorithms/+-- except that the individual solutions are programs (or, more generally, +-- /expressions/) representing a solution to a given problem. A genetic +-- program is represented as an /abstract syntax tree/ and associated +-- with a custom-defined /fitness/ value indicating the quality of the +-- solution. Starting from a randomly generated initial population of +-- genetic programs, the genetic operators of /selection/, /crossover/, +-- and (occasionally) /mutation/ are used to evolve programs of +-- increasingly better quality.+--+-- Standard reference is: John Koza. /Genetic programming:/+-- /On the Programming of Computers by Means of Natural Selection/.+-- MIT Press, 1992.+--+-- In GenProg, a genetic program is represented by a value of an+-- algebraic datatype. To use a datatype as a genetic program, it+-- suffices to define it as an instance of the 'GenProg' typeclass.+-- A custom datatype can be made an instance of the 'GenProg'+-- typeclass, provided it is an instance of the 'Data' typeclass (see+-- "GenProg.GenExpr.Data").+--+-- An example of how to use this library is given below.+--+-----------------------------------------------------------------------------++{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies,+    NoMonomorphismRestriction #-}++module GenProg (+  -- * Genetic programs+  GenProg (..),+  -- * Expressions+  generateFullExpr,+  generateGrownExpr,+  depth,+  nodes,+  -- * Individuals+  Ind,+  unInd,+  mkInd,+  aFitness,+  sFitness,+  -- * Population+  Pop,+  unPop,+  mkPop,+  generatePop,+  replenishPop,+  mergePop,+  best,+  avgFitness,+  avgDepth,+  avgNodes,+  -- * Genetic operators+  -- | The following functions are not meant to be used directly.+  -- They are exposed for debugging purposes.+  crossoverInd,+  mutateInd,+  crossoverPop,+  mutatePop,+  -- * Evolution state+  EvolState (..),+  -- * Control parameters+  Fitness,+  Mutate,+  defaultMutation,+  Terminate,+  tSuccess,+  tFitness,+  tGeneration,+  EvolParams (..),+  defaultEvolParams,+  -- * Evolution+  evolve,+  evolveFrom,+  evolveTrace,+  evolveTraceFrom+  -- * Example+  -- $Example+  ) where++import Data.List+import Data.Ord+import Data.Maybe+import Control.Monad+import Control.Monad.Random+import GenProg.GenExpr.Data++-- | A typeclass defining a genetic program interface.  Datatypes @e@+-- that are to be used as genetic programs must be instances of the+-- 'GenExpr' typeclass and must implement this interface. +class (Eq e, GenExpr e, MonadRandom m) => GenProg m e | e -> m where+  -- | Generates a random terminal @T@.+  terminal :: m e+  -- | Generates a random nonterminal (functional) node @F(T,...,T)@ whose+  -- arguments are again terminals (this condition is not verified).+  nonterminal :: m e++-----------------------------------------------------------------------------+-- Expressions++-- | Generates a random expression of a given maximum depth.+generateExpr :: (GenProg m e) => m e -> Int -> m e+generateExpr g d+  | d < 1     = error "GenProg.generateExpr: Invalid expression depth"+  | otherwise = nonterminal >>= step (d - 1)+  where step 0 _ = terminal+        step d e = nodeMapM (const g >=> step (d - 1)) e++-- | Generates a random expression fully expanded to the specified depth.+generateFullExpr :: (GenProg m e) => Int -> m e+generateFullExpr = generateExpr nonterminal++-- | Generates a random expression of limited depth. The maximum depth of+-- the resulting expression may be less than the specified depth+-- limit, and paths may be of different length.+generateGrownExpr :: (GenProg m e) => Int -> m e+generateGrownExpr d = do+  t <- getRandom+  generateExpr (if t then terminal else nonterminal) d++-----------------------------------------------------------------------------+-- Individuals++-- | A genetically programmed individual, representing a basic unit+-- of evolution. (Basically a wrapper around a genetically programmable+-- expression.)+data Ind e = Ind {+  -- | Returns the expression wrapped by an individual.+  unInd :: e,+  -- | Adjusted fitness of an individual. Adjusted fitness equals+  -- @1/(1+s)@, where @s@ is the standardized fitness as computed by+  -- 'fitness'. To reduce computational costs, this value is computed+  -- only once and then cached.+  aFitness :: Double,+  -- The indices of inner (functional) nodes of an individual's expression.+  iNodes :: [Int],+  -- The indices of external (terminal) nodes of an individual's expression.+  eNodes :: [Int] }+  deriving (Show)++instance (Eq e) => Eq (Ind e) where+  i1 == i2 = unInd i1 == unInd i2++instance (Eq e) => Ord (Ind e) where+  compare = comparing aFitness++-- | Wraps an expression into an individual.+mkInd :: (GenProg m e) => Fitness e -> e -> Ind e+mkInd f e = Ind e (adjust $ f e) fs ts+  where (fs,ts) = nodeIndices e++-- Adjusts fitness.+adjust :: Double -> Double+adjust f = 1 / (1 + max 0 f)++-- Unadjusts fitness (the inverse of adjustFitness).+unadjust :: Double -> Double+unadjust f = 1 / f - 1++-- | Standardized fitness of an individual as computed by 'fitness'+sFitness :: Ind e -> Double+sFitness = unadjust . aFitness++-----------------------------------------------------------------------------+-- Population++-- | A population of individuals. (Basically a wrapper around a list of+-- individuals.)+data Pop e = Pop+  { unPop  :: [Ind e]   -- ^ Unwraps a population.+  , dist_  :: [Double]  -- ^ Fitness distribution.+  } deriving (Show, Eq)++-- | Wraps a list of individuals into a population.+mkPop :: [Ind e] -> Pop e+mkPop is = Pop is ds+  where ds = map snd . distribution $+             map (\i -> (unInd i, aFitness i)) is++-- | Generate population of given size and given depth limit using+-- /ramped half-and-half/ method (Koza, 1992): for each depth value from 0 to+-- the initial depth limit 'iDepth', 50% of individuals are generated using+-- 'generateFullExpr' and 50% are generated using+-- 'generateGrownExpr'. Afterwards, duplicates are removed, thus the+-- size of the resulting population may actually be less than the+-- specified size.+generatePop :: (GenProg m e) => EvolParams m e -> m (Pop e)+generatePop p+  | s < 2 || n==0 = error "GenProg.generatePop: Invalid population size"+  | otherwise = do+    iss <- forM [2..di] $ \i -> do+      is1 <- replicateM n (mkInd (fitness p) `liftM` generateFullExpr di)+      is2 <- replicateM n (mkInd (fitness p) `liftM` generateGrownExpr di)+      return $ is1 ++ is2+    return . mkPop . nub $ concat iss+  where n  = s `div` (2 * (di - 1))+        s  = popSize p+        di = iDepth p++-- | Replenishes a population up to 'popSize' by randomly+-- generating new individuals.+replenishPop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Pop e)+replenishPop p pop1 = do+  pop2 <- generatePop p+  return . mkPop $ unPop pop1 ++ drop s (unPop pop2)+  where s = length $ unPop pop1++-- | Merges two populations by taking 'popSize' best-fitted individuals+-- from the union of the two populations.+mergePop :: (GenProg m e) => EvolParams m e -> Pop e -> Pop e -> Pop e+mergePop p pop1 pop2 = mkPop $ take (popSize p) is+  where is = sortBy (flip $ comparing aFitness) $ unPop pop1 ++ unPop pop2++-- | Population's best-fitted individual.+best :: Pop e -> Ind e+best = maximumBy (comparing aFitness) . unPop++avg :: (Fractional a) => [a] -> a+avg xs = sum xs / realToFrac n+  where n = length xs++-- | Population's average standardized fitness.+avgFitness :: Pop e -> Double+avgFitness = avg . map (unadjust . aFitness) . unPop++-- | Average depth of expressions in the population.+avgDepth :: (GenProg m e) => Pop e -> Double+avgDepth = avg . map (realToFrac . depth . unInd) . unPop++-- | Average number of expression nodes in the population.+avgNodes :: (GenProg m e) => Pop e -> Double+avgNodes = avg . map (realToFrac . nodes . unInd) . unPop++-----------------------------------------------------------------------------+-- Genetic operators++-- Selects at random an index of an expression node. Functional+-- (internal) nodes are selected with probability 'pci', whereas+-- terminal (external) nodes are selecred with probability '1-pi'.+selectNode :: (GenProg m e, MonadRandom m) => Double -> Ind e -> m Int+selectNode pi i+  | null $ iNodes i = oneof $ eNodes i+  | otherwise       = choice pi (oneof $ iNodes i) (oneof $ eNodes i)++-- | Crossover operation of two individuals, resulting in two+-- offsprings. Crossover is performed by choosing at random two nodes+-- in each expressions, and then by exchanging the subexpressions+-- rooted at these nodes between the two individuals. The probability+-- that an internal (functional) node is chosen as crossover point is+-- set by the 'ciProb' parameter in 'EvolParams', whereas the+-- probability that an external (terminal) node is chosen equals+-- @1-ciProb@. Among internal and external nodes, nodes are chosen+-- uniformly at random. If the depth of a created offspring exceeds+-- the depth limit 'cDepth' specified by evolution parameters+-- 'EvolParams', that offspring is discarded and a parent is+-- reproduced (i.e., copied as-is).+crossoverInd :: (GenProg m e) =>+  EvolParams m e -> Ind e -> Ind e -> m (Ind e, Ind e)+crossoverInd p i1 i2 = do+  n1 <- selectNode (ciProb p) i1+  n2 <- selectNode (ciProb p) i2+  let (r1,r2) = exchange (unInd i1) n1 (unInd i2) n2+  return (if depth r1 <= cDepth p then mkInd (fitness p) r1 else i1,+          if depth r2 <= cDepth p then mkInd (fitness p) r2 else i2)++-- | Mutates an individual by applying the mutation function @mutate@+-- to a randomly selected node. The probability that an internal+-- (functional) node is chosen for muration is set by the 'miProb'+-- parameter in 'EvolParams', whereas the probability that an external+-- (terminal) node is chosen equals @1-miProb@. Among internal and+-- external nodes, nodes are chosen uniformly at random. If the depth+-- of the mutated expression exceeds the depth limit 'cDepth'+-- specified by evolution parameters 'EvolParams', the individual is+-- left unaltered.+mutateInd :: (GenProg m e) => EvolParams m e -> Ind e -> m (Ind e)+mutateInd p i = do+  n  <- selectNode (miProb p) i+  e2 <- adjustM (mutate p) e1 n+  return . mkInd (fitness p) $ if depth e2 <= cDepth p then e2 else e1+  where e1 = unInd i++-- Discrete distribution.+type Distribution a = [(a, Double)]++-- Computes distribution from a weighted list.+-- The weights need not sum to 1.+distribution :: [(a, Double)] -> Distribution a+distribution xs = [(x,f i) | ((x,_),i) <- zip xs [1..]]+  where f i = sum . map snd $ take i ys+        s   = sum $ map snd xs+        ys  = map (\(x, w) -> (x, w/s)) xs++-- Samples a value from a discrete distribution.+choose :: (MonadRandom m) => Distribution a -> m a+choose xs = do+  p <- getRandomR (0,1)+  return . fst . fromJust $ find ((>= p) . snd) xs++-- Chose first action with probability 'p' and second with probability+-- 1-p.+choice :: (MonadRandom m) => Double -> m a -> m a -> m a+choice p a1 a2 = do+  r <- getRandomR (0,1)+  if r <= p then a1 else a2++oneof :: (MonadRandom m) => [a] -> m a+oneof xs = (xs!!) `liftM` getRandomR (0,length xs-1)++-- Fitness-proportionate selection of an individual from a population.+selectInd :: (MonadRandom m) => Pop e -> m (Ind e)+selectInd pop = choose (zip (unPop pop) (dist_ pop))++reproducePop :: (MonadRandom m) => Pop e -> m (Ind e)+reproducePop = selectInd++-- | Applies crossover to two randomly chosen individuals from a+-- population. The probability of an individual being chosen as parent+-- is fitness-proportionate (individuals with better fitness have+-- better chanches of being chosen for crossover).+crossoverPop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Ind e,Ind e)+crossoverPop p pop = do+  i1 <- selectInd pop+  i2 <- selectInd pop+  crossoverInd p i1 i2++-- | Applies mutation operation to individuals from a population. The+-- probability of mutating each individual is determined by 'mProb' parameter+-- from 'EvalParams'.+mutatePop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Pop e)+mutatePop p pop+  | mProb p == 0 = return pop+  | otherwise    = liftM mkPop . forM (unPop pop) $ \i ->+                     choice (mProb p) (mutateInd p i) (return i)++-----------------------------------------------------------------------------+-- Evolution state++-- | The state of the evolution.+data EvolState e = EvolState+  { pop        :: Pop e    -- ^ Current population.+  , iter       :: Int      -- ^ Iteration (current generation number).+  , cachedBest :: Ind e    -- ^ Best individual evolved so far.+  } deriving (Show,Eq)++initState :: Pop e -> EvolState e+initState pop =+  EvolState { pop = pop, iter = 0, cachedBest = best pop }++-- | Advances to next evolution state.+nextState :: (GenProg m e ) =>+  EvolParams m e -> EvolState e -> m (EvolState e)+nextState p es1 = do+  pop2 <- evolvePop p pop1+  return $ es1 { pop = pop2, iter = iter es1 + 1,+                 cachedBest = max (cachedBest es1) (best pop1) }+  where pop1 = pop es1++-----------------------------------------------------------------------------+-- Control parameters++-- | Standardized fitness. It takes on values from 0 (best fitness) to+-- +infinity (worst fitness).+type Fitness e = e -> Double++-- | A function to mutate a chosen expression node.+type Mutate m e = e -> m e++-- | Default mutation. Replaces a node, irrespective of its value,+-- with a randomly generated subexpression whose depth is limited to+-- 'iDepth'.+defaultMutation :: (GenProg m e) => EvolParams m e -> Mutate m e+defaultMutation p = const $ generateGrownExpr (iDepth p)++-- | Termination predicate.+type Terminate e = EvolState e -> Bool++-- | Termination predicate: terminate if any individual satisfies the+-- specified predicate.+tSuccess :: (e -> Bool) -> Terminate e+tSuccess c = any (c . unInd) . unPop . pop++-- | Termination predicate: terminate if best individual's+-- standardized fitness is greater than or equal to the specified value.+tFitness :: (GenProg m e) => Double -> Terminate e+tFitness f = (>= f) . unadjust . aFitness . cachedBest++-- | Termination predicate: terminate after running for the specified+-- number of iterations.+tGeneration :: Int -> Terminate e+tGeneration n = (>=n) . iter++-- | Parameters governing the evolution.+--+-- Default evolution parameters,+-- as used in (Koza, 1992), are defined by 'defaultEvolParams'+-- and indicated below. At least the fitness function 'fitness' should+-- be overriden.+data EvolParams m e = EvolParams {+  -- | Population size (number of individuals). Default is @500@.+  popSize   :: Int,+  -- | Depth of expressions in initial population. Default is @6@.+  iDepth    :: Int,+  -- | Maximum depth of expressions created during the evolution.+  -- Default is @17@.+  cDepth    :: Int,+  -- | Probability of crossover. Default is @0.9@. If crossover is not+  -- chosen, an individual is simply reproduced (copied as-is) into+  -- the next generation.+  cProb     :: Double,+  -- | Probability that an internal (functional) node is chosen as a+  -- crossover point. Default is @0.9@. If an internal node is not+  -- chosen, an external (terminal) node is+  -- chosen.+  ciProb    :: Double,+  -- | Probability that an individual gets mutated. Default is @0@+  -- (no mutation).+  mProb     :: Double,+  -- | Probability that an internal (functional) node is chosen for+  -- mutation. Default is @0.1@.+  miProb    :: Double,+  -- | Standardized fitness function. Default value is @undefined@+  -- (must be overriden).+  fitness   :: Fitness e,+  -- | Mutation function. Defines how to change a randomly chosen+  -- node. Default is @defaultMutation defaultEvolParams@+  -- (replacement of a chosen node with a randomly generated subexpression).+  mutate    :: Mutate m e,+  -- | Elitist factor: number of best-fitted individuals that are preserved+  -- from each generation (reproduced as-is into next evolution state).+  -- Default is @0@.+  elitists  :: Int,+  -- | Termination predicate. Default is @50@ (terminate after 50 generations).+  terminate :: Terminate e }++defaultEvolParams = EvolParams+  { popSize   = 500+  , iDepth    = 6+  , cDepth    = 17+  , cProb     = 0.9+  , ciProb    = 0.9+  , mProb     = 0.0+  , miProb    = 0.1+  , terminate = tGeneration 50+  , fitness   = error "GenProg.defaultEvolParams: fitness function is undefined"+  , mutate    = const $ generateGrownExpr (iDepth defaultEvolParams)+  , elitists  = 0 }++-----------------------------------------------------------------------------+-- Evolution++untilM :: (Monad m) => (a -> Bool) -> (a -> m a) -> a -> m a+untilM p f x | p x       = return x+             | otherwise = f x >>= untilM p f++iterateUntilM :: (Monad m) => (a -> Bool) -> (a -> m a) -> a -> m [a]+iterateUntilM p f x+  | p x       = return []+  | otherwise = do y  <- f x+                   ys <- iterateUntilM p f y+                   return (y:ys)++-- | Evolves one population from another one by performing a single+-- evolution step.+evolvePop :: (GenProg m e) => EvolParams m e -> Pop e -> m (Pop e)+evolvePop p pop1 = do+     pop2 <- mkPop `liftM` untilM ((>= s) . length) step []+     pop3 <- mutatePop p pop2+     return $ mkPop (elite ++ unPop pop3)+  where s = popSize p - length elite+        elite = take (elitists p) topRanked+        topRanked = sortBy (flip $ comparing aFitness) $ unPop pop1+        step is | length is == s - 1 = (:is) `liftM` reproducePop pop1+                | otherwise = choice (cProb p)+                    (do (i1,i2) <- crossoverPop p pop1; return (i1:i2:is))+                    ((:is) `liftM` reproducePop pop1)++-- | Creates an initial population and evolves it until termination+-- predicate is satisfied, returning the last evolution state.+evolve :: (GenProg m e) => EvolParams m e -> m (EvolState e)+evolve p = -- generatePop p >>= evolveFrom p+  last `liftM` evolveTrace p++-- | Evolves a given initial population until termination+-- predicate is satisfied, returning the last evolution state.+-- If the size of the initial population is less than+-- 'popSize', the population will be replenished (see 'replenishPop').+evolveFrom :: (GenProg m e) => EvolParams m e -> Pop e -> m (EvolState e)+evolveFrom p pop = -- untilM (terminate p) (nextState p) . initState+  last `liftM` evolveTraceFrom p pop++-- | Runs evolution on a given initial population until termination+-- predicate is satisfied and returns a list of successive evolution+-- states. If the size of the initial population is less than+-- 'popSize', the population will be replenished (see 'replenishPop').+evolveTraceFrom :: (GenProg m e) => EvolParams m e -> Pop e -> m [EvolState e]+evolveTraceFrom p pop1 =+  iterateUntilM (terminate p) (nextState p) . initState =<< replenishPop p pop1++-- | Creates an initial population and runs evolution until+-- termination predicate is satisfied. Returns a list of successive+-- evolution states.+evolveTrace :: (GenProg m e) => EvolParams m e -> m [EvolState e]+evolveTrace p = generatePop p >>= evolveTraceFrom p++-----------------------------------------------------------------------------+-- Example++{- $Example++This is a simple, worked through example of how to use the GenProg+library. Given a target number @n@, out aim is to evolve an arithmetic+expression that evaluates to @n@. For example, given @13@ as the+target number, one possible solution is @(3 * 5) - 2@. The constants+allowed to appear in the expression are restricted to integers from 1+to 9. The allowed operations are @+@, @-@, @*@, and integer division+without remainder.++We begin by defining the datatype for the genetically programed+expression:++@+-- The following language extensions need to be enabled:+-- DeriveDataTypeable, FlexibleInstances, MultiParamTypeClasses++import GenProg+import Data.Generics+import Control.Monad+import Control.Monad.Random++data E = Plus E E+       | Minus E E+       | Times E E+       | Div E E+       | Const Int+       deriving (Typeable,Data,Eq,Show)+@++In order to evolve arithmetic expressions, we need to be able to+compute their values. To this end we define++@+eval :: E -> Maybe Int+eval (Const c)     = Just c+eval (Plus e1 e2)  = liftM2 (+) (eval e1) (eval e2)+eval (Minus e1 e2) = liftM2 (-) (eval e1) (eval e2)+eval (Times e1 e2) = liftM2 (*) (eval e1) (eval e2)+eval (Div e1 e2) | ok        = liftM2 div x1 x2+                 | otherwise = Nothing+  where (x1,x2) = (eval e1,eval e2)+        ok = x2 /= Just 0 && liftM2 mod x1 x2 == Just 0+@++Dividing by zero and dividing with a remainder are not allowed and in+such cases we return @Nothing@.++Because we have made @E@ an instance of the 'Data' typeclass, it can+be readily used as a genetically programmable expression. Next step is+to make 'E' an instance of the 'GenProg' typeclass:++@+instance GenProg (Rand StdGen) E where+  terminal    = Const `liftM` getRandomR (1,9)+  nonterminal = do+    r <- getRandomR (0,3)+    [liftM2 Plus terminal terminal,+     liftM2 Minus terminal terminal,+     liftM2 Times terminal terminal,+     liftM2 Div terminal terminal] !! r+@++Thus, a random terminal node contains one of the constants from 1 to+9. A nonterminal node can be one of the four arithmetic operations,+each with terminal nodes as arguments.  Note that computations are run+within the standard random generator monad (@Rand StdGen@).++The fitness function evaluates the accurateness of the arithmetic+expression with respect to the target number. If the value of the+expression is far off from the target number @n@, the standardized+fitness should be high. Moreover, we would like to keep the expression+as simple as possible. To this end, we include a /parsimony factor/+that is proportional to the number of nodes an expression has. We+define the overall standardized fitness as++@+myFitness :: Int -> E -> Double+myFitness n e = error + size+  where error = realToFrac $ maybe maxBound (abs . (n-)) (eval e)+        size  = (realToFrac $ nodes e) / 100+@++The number of nodes is divided by a factor of 100 to make it less+important than the numeric accuracy of the expression.++We now have everything in place to get the evolution going. We will use+default evolution parameters and choose @12345@ as the target number:++>>> let params = defaultEvolParams { fitness = myFitness 12345 }++Let us first create a random number generator: ++>>> let g = mkStdGen 0++We are doing this because we want our results to be reproducible, and+because we want to be able to compare the results of different+evolution runs. Normally, you would use @getStdGen@ to get a random+generator with random seed.++To run the evolution and get the best evolved individual, we type++>>> let i = cachedBest $ evalRand (evolve params) g++To check out its standardized fitness, we type++>>> sFitness i+39.61++Let us see how the actual expression looks like:++>>> unInd i+Times (Minus (Minus (Minus (Plus (Const 4) (Const 4)) (Plus (Const 6) +(Const 7))) (Minus (Minus (Const 5) (Const 9)) (Plus (Minus (Const 5) +(Const 9)) (Minus (Const 4) (Const 4))))) (Plus (Times (Plus (Const 5) +(Const 1)) (Const 6)) (Times (Plus (Const 9) (Const 3)) (Minus (Const 1) +(Const 8))))) (Div (Times (Plus (Plus (Const 3) (Const 5)) (Times (Const 4) +(Const 7))) (Plus (Const 4) (Const 4))) (Minus (Minus (Plus (Const 2) +(Const 8)) (Plus (Const 6) (Const 7))) (Plus (Minus (Const 5) (Const 9)) +(Minus (Const 4) (Const 4)))))++The number of nodes is++>>> nodes $ unInd i+61++Let us see to what number the expression evaluates:++>>> eval $ unInd i+Just 12384++So in this run we didn't get a perfect match, but we were close. Let+us see if we can do better.++When doing genetic programming, it is always a good idea to experiment+a bit with the parameters. There are no parameters that work best for+any given problem. You can learn a lot about how parameters influence+the evolution by analysing how the evolution progresses in time. This+can be accomplised by evolving an evolution trace:++>>> let trace = evalRand (evolveTrace params) g++We can now analyse how the standardized fitness of the+best individual improves during the evolution:++>>> map (sFitness . best . pop) trace+[9591.35,2343.59,1935.59,2343.59,903.51,903.45,585.59,585.59,327.45,225.41,+225.41,135.43,57.49,39.61,39.61,39.61,39.61,39.61,57.43,57.47,57.43,57.45,+57.33,57.43,57.43,57.45,57.43,57.43,57.35,57.35,57.43,57.27,57.33,57.33,57.43,+57.29,57.33,57.41,57.29,57.43,57.33,57.35,57.35,57.33,57.39,57.39,57.39,57.33,+57.37,57.37]++We see that at some point the fitness decreases and then increases+again. This indicates that the best fitted individual was lost by+evolving from one generation to the other. We can prevent this by+employing the /elitist strategy/. Let us see what happens if we+preserve a best fitted individual in each generation:++>>> let trace = evalRand (evolveTrace params {elitists = 1}) g +>>> map (sFitness . best . pop) trace+[9591.35,2343.59,711.61,711.61,711.61,711.61,57.55,57.53,57.39,57.39,57.39,+57.39,57.37,57.37,57.37,57.37,57.37,57.37,57.37,57.37,57.35,57.35,57.35,+57.35,57.35,57.35,57.35,57.35,57.35,57.35,57.33,57.33,57.33,57.33,57.33,+57.33,57.33,57.33,57.33,25.31,25.31,25.31,25.31,25.31,25.31,25.296,25.296,+25.296,25.296,25.296]++This gives us better fitness, but still not an exact match:++>>> let i = cachedBest $ last trace+>>> eval $ unInd i+Just 12320++In the previous evolution run fitness converged relatively fast, but then+remained stuck. To stir up things a little, let us allow for some+mutation. Setting mutation probability to 5%, while retaining the+elitist strategy, we get++>>> let trace = evalRand (evolveTrace params {elitists = 1, mProb = 0.05}) g+>>> map (sFitness . best . pop) trace+[9591.35,9591.35,9591.35,9591.35,9591.35,9591.35,9159.35,8403.23,7239.11,+6087.15,6087.15,1479.13,819.21,60.13,51.19,5.19,5.19,5.19,5.19,5.19,1.23,+1.23,1.23,1.23,1.23,1.23,1.21,1.21,1.21,1.21,0.23998,0.23998,0.23998,0.23998,+0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,+0.23998,0.23998,0.23998,0.23998,0.23998,0.23998,0.23998]++This time we've got a perfect match:++>>> let i = cachedBest $ last trace+>>> eval $ unInd i+Just 12345++while at the same time the expression is rather compact:++>>> unInd i+Plus (Times (Const 4) (Plus (Const 9) (Const 4))) (Plus (Plus (Times +(Plus (Const 4) (Const 3)) (Times (Times (Const 3) (Const 9)) (Times +(Const 5) (Plus (Const 9) (Const 4))))) (Const 3)) (Const 5))+>>> nodes $ unInd i+23++-}
+ src/GenProg/GenExpr.hs view
@@ -0,0 +1,59 @@+-- |+-- Module      :  GenProg.GenExpr+-- Copyright   :  (c) 2010 Jan Snajder+-- License     :  BSD-3 (see the LICENSE file)+--+-- Maintainer  :  Jan Snajder <jan.snajder@fer.hr>+-- Stability   :  experimental+-- Portability :  non-portable+--+-- An interface to genetically programmable expressions.+--+-----------------------------------------------------------------------------++module GenProg.GenExpr (+  GenExpr (..)) where++import Control.Monad++-- | This typeclass defines an interface to expressions+-- that can be genetically programmed.  The operations that must be+-- provided by instances of this class are used for the generation+-- of random individuals as well as crossover and mutation operations.+-- (An instance for members of the @Data@ typeclass is provided in+-- "GenProg.GenExpr.Data".)+--+-- Minimal complete definition: 'exchange', 'nodeMapM', 'nodeMapQ',+-- and 'nodeIndices'.+class GenExpr e where+  -- | Exchanges subtrees of two expressions:+  -- @exchange e1 n1 e2 n2@ replaces the subexpression of @e1@ rooted in node+  -- @n1@ with the subexpression of @e2@ rooted in @n2@, and vice versa.+  exchange :: e -> Int -> e -> Int -> (e, e)+  -- | Maps a monadic transformation function over the immediate+  -- children of the given node.+  nodeMapM :: Monad m => (e -> m e) -> e -> m e+  -- | Maps a query function over the immediate children of the given+  -- node and returns a list of results.+  nodeMapQ :: (e -> a) -> e -> [a]+  -- | A list of indices of internal (functional) and external+  -- (terminal) nodes of an expression.+  nodeIndices :: e -> ([Int], [Int])+  -- | Adjusts a subexpression rooted at the given node by applying a+  -- monadic transformation function.+  adjustM :: (Monad m) => (e -> m e) -> e -> Int -> m e+  -- | Number of nodes an expression has.+  nodes :: e -> Int+  -- | The depth of an expression. Equals 1 for single-node expressions.+  depth :: e -> Int+++  -- | Default method (expensive because it calls exchange twice).+  adjustM f e n = replace e n `liftM` f (get e n)+    where get e n = fst $ exchange e 0 e n+          replace e1 n1 e2 = fst $ exchange e1 n1 e2 0++  nodes = (+1) . foldr (+) 0 . nodeMapQ nodes ++  depth = (+1) . foldr max 0 . nodeMapQ depth+ 
+ src/GenProg/GenExpr/Data.hs view
@@ -0,0 +1,153 @@+-- |+-- Module      :  GenProg.GenExpr.Data+-- Copyright   :  (c) 2010 Jan Snajder+-- License     :  BSD-3 (see the LICENSE file)+--+-- Maintainer  :  Jan Snajder <jan.snajder@fer.hr>+-- Stability   :  experimental+-- Portability :  non-portable+--+-- Implementation of the @GenProg.GenExpr@ interface for members of+-- the 'Data' typeclass. The implementation is based on SYB and SYZ+-- generic programming frameworks (see+-- <http://hackage.haskell.org/package/syb> and+-- <http://hackage.haskell.org/package/syz> for details).+--+-- NB: Subexpressions that are candidates for crossover points or+-- mutation must be of the same type as the expression itself, and+-- must be reachable from the root node by type-preserving traversal.+-- See below for an example.+--+-----------------------------------------------------------------------------++{-# LANGUAGE ScopedTypeVariables, FlexibleInstances, Rank2Types,+    UndecidableInstances, DeriveDataTypeable #-}++module GenProg.GenExpr.Data (+  -- | This module re-exports @GenExpr@ typeclass.+  GenExpr (..)+  -- * Example+  -- $Example+  ) where++import Data.Generics+import Data.Generics.Zipper+import Data.Maybe+import Control.Monad+import GenProg.GenExpr++moduleName = "GenProg.GenExpr.Data"++instance (Data a) => GenExpr a where++  -- | Exchanges two expression nodes. Works by using two generic+  -- zippers and exchanging their holes.+  exchange e1 n1 e2 n2 = (fromZipper y1, fromZipper y2)+    where z1 = typeMoveForUnsafe n1 $ toZipper e1+          z2 = typeMoveForUnsafe n2 $ toZipper e2+          (y1,y2) = exchangeHoles z1 z2++  -- | Adjust an expression node. Works by applying a monadic+  -- tranformation on a zipper hole.+  adjustM f e n = fromZipper `liftM` transM (mkM f) z+    where z = typeMoveForUnsafe n (toZipper e)++  nodeMapM f = gmapM (mkM f)++  nodeMapQ q (x::a) = concat $ gmapQ ([] `mkQ` (\(y::a) -> [q y])) x++  nodeIndices = index 0 [] [] . toZipper++-- Zipper moves++type Move a = Zipper a -> Maybe (Zipper a)++backtrack :: (Typeable a) => Move a+backtrack z = do+  z2 <- up z+  right z2 `mplus` backtrack z2++repeatM :: (Monad m) => Int -> (a -> m a) -> a -> m a+repeatM 0 _ x = return x+repeatM n f x = f x >>= repeatM (n - 1) f++-- Moves zipper to next node in DFS order, but does not move down the+-- zipper if node satisfies query 'q'.+nextDfsQ :: Typeable a => GenericQ Bool -> Move a+nextDfsQ q z = (if query q z then Nothing else down' z)+  `mplus` right z `mplus` backtrack z++-- Moves the zipper to node 'n' from current position in DFS order,+-- skipping nodes not satisfying query 'q2' and descending only down+-- the nodes satisfying query 'q1'.+moveForQ :: (Typeable a) => GenericQ Bool -> GenericQ Bool -> Int -> Move a+moveForQ _  _  0 z = Just z+moveForQ q1 q2 n z = do+  z2 <- nextDfsQ q1 z+  moveForQ q1 q2 (if query q2 z2 then n - 1 else n) z2++-- Moves the zipper to node 'n' from current position in DFS order,+-- counting only nodes of type 'a', and not descending down the nodes+-- of other type.+typeMoveFor :: (Typeable a) => Int -> Move a+typeMoveFor n (z::Zipper a) =+  moveForQ (True `mkQ` (\(_::a) -> False)) (False `mkQ` (\(_::a) -> True)) n z++-- | Same as typeMoveFor, but throws an error if node index is out of+-- bound.+typeMoveForUnsafe :: (Typeable a) => Int -> Zipper a -> Zipper a+typeMoveForUnsafe n z = fromMaybe+  (error $ moduleName ++ ".typeMoveForUnsafe: Nonexisting node.")+  (typeMoveFor n z)++-- | Exchanges two zipper holes.+exchangeHoles :: (Data a) => Zipper a -> Zipper a -> (Zipper a, Zipper a)+exchangeHoles (z1::Zipper a) (z2::Zipper a) = (y1,y2)+  where Just h1 = getHole z1 :: Maybe a+        Just h2 = getHole z2 :: Maybe a+        y1 = setHole h2 z1+        y2 = setHole h1 z2++index :: (Data a) => Int -> [Int] -> [Int] -> Zipper a -> ([Int], [Int])+index i is es (z :: Zipper a) =+  maybe (is2,es2) (index (i + 1) is2 es2) (typeMoveFor 1 z)+  where Just h = getHole z :: Maybe a+        (is2,es2) = if terminalQ h then (is,i:es) else (i:is,es)++terminalQ :: (Data a) => a -> Bool+terminalQ = null . nodeMapQ id++{- $Example++Suppose you have a datatype defined as++@+data E = A E E+       | B String [E]+       | C+ deriving (Eq,Show,Typeable,Data)+@++and an expression defined as++@+e = A (A C C) (B \"abc\" [C,C])+@++The subexpressions of a @e@ are considered to be only the subvalues of+@e@ that are of the same type as @e@.  Thus, the number of nodes of+expression @e@ is++>>> nodes e+5+ +because subvalues of node @B@ are of different type than expression+@e@ and therefore not considered as subexpressions. ++Consequently, during a genetic programming run, subexpressions that+are of a different type than the expression itself, or subexpression+that cannot be reached from the root node by a type-preserving+traversal, cannot be chosen as crossover points nor can they be+mutated.++-}