diff --git a/AI/MEP.hs b/AI/MEP.hs new file mode 100644 --- /dev/null +++ b/AI/MEP.hs @@ -0,0 +1,87 @@ +{- | +Copyright Bogdan Penkovsky (c) 2017 + += Multiple Expression Programming + +== Example application: trigonometry cheating + + Suppose, you forgot certain trigonometric identities. + For instance, you want to express cos^2(x) using sin(x). + No problem, set the target function cos^2(x) in the dataset + and add sin to the arithmetic set of operators @{+,-,*,/}@. + See app/Main.hs. + + After running + + @ + $ stack build && stack exec hmep-demo + @ + + We obtain + + @ + Average loss in the initial population 15.268705681244962 + Population 10: average loss 14.709728527360586 + Population 20: average loss 13.497114190675477 + Population 30: average loss 8.953185872653737 + Population 40: average loss 8.953185872653737 + Population 50: average loss 3.3219954564955856e-15 + @ + + The value of 3.3e-15 is zero with respect to the + rounding errors. It means that the exact expression was found! + + The produced output was: + + @ + Interpreted expression: + v1 = sin x0 + v2 = v1 * v1 + result = 1 - v2 + @ + + From here we can infer that + @ + cos^2(x) = 1 - v2 = 1 - v1 * v1 = 1 - sin^2(x) + @ + + Sweet! + +-} + +module AI.MEP ( + Chromosome (..) + , Gene + , Population + , Phenotype + , Config (..) + , defaultConfig + , LossFunction + + -- * Genetic algorithm + , initialize + , evaluateGeneration + , avgLoss + , evolve + , binaryTournament + , crossover + , mutation3 + , smoothMutation + , newChromosome + + -- * Expression interpretation + , generateCode + + -- * Random + , Rand + , newPureMT + , runRandom + , evalRandom + ) where + +import System.Random.Mersenne.Pure64 ( newPureMT ) + +import AI.MEP.Types +import AI.MEP.Operators +import AI.MEP.Run +import AI.MEP.Random diff --git a/AI/MEP/Operators.hs b/AI/MEP/Operators.hs new file mode 100644 --- /dev/null +++ b/AI/MEP/Operators.hs @@ -0,0 +1,288 @@ +-- | +-- = Genetic operators + +module AI.MEP.Operators ( + Config (..) + , defaultConfig + , LossFunction + -- * Genetic operators + , initialize + , evaluateGeneration + , evolve + , phenotype + , binaryTournament + , crossover + , mutation3 + , smoothMutation + , newChromosome + ) where + +import Data.Vector ( Vector ) +import qualified Data.Vector as V +import Data.List + ( nub + , sortBy + ) +import Data.Ord ( comparing ) +import qualified Control.Monad as CM + +import AI.MEP.Random +import AI.MEP.Types +import AI.MEP.Run ( evaluate ) + +data Config a = Config + { + p'const :: Double -- ^ Probability of constant generation + , p'var :: Double -- ^ Probability of variable generation. + -- The probability of operator generation is inferred + -- automatically as @1 - p'const - p'var@. + , p'mutation :: Double -- ^ Mutation probability + , p'crossover :: Double -- ^ Crossover probability + + , c'length :: Int -- ^ The chromosome length + , c'popSize :: Int -- ^ A (sub)population size + , c'popN :: Int -- ^ Number of subpopulations (1 or more) [not implemented] + , c'ops :: Vector (F a) -- ^ Functions pool with their symbolic + -- representations + , c'vars :: Int -- ^ The input dimensionality + } + +-- | +-- @ +-- defaultConfig = Config +-- { +-- p'const = 0.1 +-- , p'var = 0.4 +-- , p'mutation = 0.1 +-- , p'crossover = 0.9 +-- +-- , c'length = 50 +-- , c'popSize = 100 +-- , c'popN = 1 +-- , c'ops = V.empty -- <-- To be overridden +-- , c'vars = 1 +-- } +-- @ +defaultConfig :: Config Double +defaultConfig = Config + { + p'const = 0.1 + , p'var = 0.4 + , p'mutation = 0.1 + , p'crossover = 0.9 + + , c'length = 50 + , c'popSize = 100 + , c'popN = 1 + , c'ops = V.empty + , c'vars = 1 + } + +-- | A function to minimize. +-- +-- The argument is a vector evaluation function whose input +-- is a vector (length @c'vars@) and ouput is +-- a vector with a different length @c'length@. +-- +-- The result is a vector of the best indices +-- and a scalar loss value. +type LossFunction a = + ((V.Vector a -> V.Vector a) -> (V.Vector Int, Double)) + +-- | Evaluates a chromosome according to the given +-- loss function. +phenotype + :: Num a => + LossFunction a + -> Chromosome a + -> Phenotype a +phenotype loss chr = let (is, val) = loss (evaluate chr) + in (val, chr, is) + +-- | Randomly generate a new population +initialize :: Config Double -> Rand (Population Double) +initialize c@Config { c'popSize = size } = mapM (\_ -> newChromosome c) [1..size] + +evaluateGeneration + :: Num a => + LossFunction a + -> [Chromosome a] + -> [Phenotype a] +evaluateGeneration loss = map (phenotype loss) + +-- | Selection operator that produces the next evaluated population. +-- +-- Standard algorithm: the best offspring O replaces the worst +-- individual W in the current population if O is better than W. +evolve + :: + Config Double + -- ^ Common configuration + -> LossFunction Double + -- ^ Custom loss function + -> (Chromosome Double -> Rand (Chromosome Double)) + -- ^ Mutation + -> (Chromosome Double -> Chromosome Double -> Rand (Chromosome Double, Chromosome Double)) + -- ^ Crossover + -> ([Phenotype Double] -> Rand (Chromosome Double)) + -- ^ A chromosome selection algorithm. Does not need to be random, but may be. + -> [Phenotype Double] + -- ^ Evaluated population + -> Rand [Phenotype Double] + -- ^ New generation +evolve c loss mut cross select phenotypes = do + let pc = p'crossover c + pm = p'mutation c + -- Sort in decreasing @val@ order so that + -- the worst (with the biggest loss) is in the head + sort' = sortBy (comparing (\(val, _, _) -> negate val)) + + ev phen0 _ = do + chr1 <- select phen0 + chr2 <- select phen0 + (of1, of2) <- withProbability pc (uncurry cross) (chr1, chr2) + of1' <- withProbability pm mut of1 + of2' <- withProbability pm mut of2 + let r1@(val1, _, _) = phenotype loss of1' + r2@(val2, _, _) = phenotype loss of2' + (worstVal, _, _) = head phen0 + phen' | val1 < worstVal = r1 : tail phen0 + | val2 < worstVal = r2 : tail phen0 + -- No change + | otherwise = phen0 + let phen1 = sort' phen' + return phen1 + + CM.foldM ev (sort' phenotypes) [1..c'popSize c `div` 2] + +-- | Binary tournament selection +binaryTournament :: Ord a => [Phenotype a] -> Rand (Chromosome a) +binaryTournament phen = do + (val1, cand1, _) <- draw $ V.fromList phen + (val2, cand2, _) <- draw $ V.fromList phen + if val1 < val2 + then return cand1 + else return cand2 + +-- | Uniform crossover operator +crossover :: + Chromosome a + -> Chromosome a + -> Rand (Chromosome a, Chromosome a) +crossover ca cb = do + r <- V.zipWithM (curry (swap 0.5)) ca cb + return $ V.unzip r + +swap :: Double -> (t, t) -> Rand (t, t) +swap p = withProbability p (\(a, b) -> return (b, a)) + +replaceAt :: Int -> a -> Vector a -> Vector a +replaceAt i gene chr0 = + let (c1, c2) = V.splitAt i chr0 + in c1 V.++ V.singleton gene V.++ V.tail c2 + +-- | Mutation operator with up to three mutations per chromosome +mutation3 :: + Config Double + -- ^ Common configuration + -> Chromosome Double + -> Rand (Chromosome Double) +mutation3 c chr = do + -- Subtract 1 to get a non-zero head to + -- replace + is <- nub <$> CM.replicateM k (getMaxInt (chrLen - 1)) + genes <- mapM new' is + let chr' = foldr (uncurry replaceAt) + chr + (zip is genes) + return chr' + where chrLen = V.length chr + k = 3 + new' = new (p'const c) (p'var c) (c'vars c) (c'ops c) + +-- | Mutation operator with a fixed mutation probability +-- of each gene +smoothMutation + :: + Double + -- ^ Probability of gene mutation + -> Config Double + -- ^ Common configuration + -> Chromosome Double + -> Rand (Chromosome Double) +smoothMutation p c chr = + let new' = new (p'const c) (p'var c) (c'vars c) (c'ops c) + mutate i = withProbability p (\_ -> new' i) + in V.zipWithM mutate (V.enumFromN 0 (V.length chr)) chr + +-- | Randomly initialize a new chromosome. +-- By definition, the first gene is terminal (a constant +-- or a variable). +newChromosome :: + Config Double -- ^ Common configuration + -> Rand (Chromosome Double) +newChromosome c = do + let pConst = p'const c + pVar = p'var c + V.mapM (new pConst pVar (c'vars c) (c'ops c)) $ V.enumFromN 0 (c'length c) + +-- | Produce a new random gene +new :: + Double -- ^ Probability to produce a constant + -> Double -- ^ Probability to produce a variable + -> Int -- ^ Number of input variables + -> Vector (F Double) -- ^ Operations vector + -> Int -- ^ Maximal operation index + -> Rand (Gene Double Int) +new p1 p2 vars ops maxIndex = if maxIndex == 0 + -- The head must be a terminal + -- p1' = p1 + (1 - p1 - p2) / 2 = 1/2 + p1/2 - p2/2 + then let p1' = 0.5 * (1 + p1 - p2) + in newTerminal p1' vars + else do + p' <- getDouble + let sel | p' < p1 = newC + | p' < (p1 + p2) = newVar vars + | otherwise = newOp ops maxIndex + sel + +newTerminal :: + Double -- ^ Probability @p@ of a constant generation. + -- @1-p@ will be the probability of a variable generation. + -> Int -- ^ Number of input variables + -> Rand (Gene Double i) +newTerminal p vars = do + p' <- getDouble + if p' < p + then newC + else newVar vars + +-- | A randomly generated variable identifier +newVar :: Int -> Rand (Gene a i) +newVar vars = do + var <- draw $ V.enumFromN 0 vars + return $ Var var + +-- | A random operation from the operations vector +newOp + :: Vector (F a) + -> Int + -> Rand (Gene a Int) +newOp ops maxIndex = do + op <- draw ops + i1 <- getMaxInt maxIndex + i2 <- getMaxInt maxIndex + return $ Op op i1 i2 + +-- | Draw a constant from the normal distribution +newCNormal + :: Double -- ^ Mean + -> Double -- ^ Std deviation + -> Rand (Gene Double i) +newCNormal mu sigma = do + n <- getNormal + return $ C (mu + sigma*n) + +-- | Draw a constant from the uniform distribution +newC :: Rand (Gene Double i) +newC = C <$> getDouble diff --git a/AI/MEP/Random.hs b/AI/MEP/Random.hs new file mode 100644 --- /dev/null +++ b/AI/MEP/Random.hs @@ -0,0 +1,55 @@ +module AI.MEP.Random + ( + -- * Utilities + draw + , getNormal + , getMaxInt + , withProbability + + -- * Re-exports + , getBool, getInt, getWord, getDouble + , runRandom, evalRandom + , Rand, Random + ) where + +import Control.Monad.Mersenne.Random +import Data.Complex (Complex (..)) +import System.Random +import Data.Vector as V + +-- | Randomly draw an element from a vector +draw :: Vector a -> Rand a +draw xs = + Rand $ \g -> let (n, g') = randomR (0, V.length xs - 1) g + r = xs V.! n + in R r g' + +-- | Modify value with probability @p@ +withProbability + :: Double -- ^ The probability @p@ + -> (a -> Rand a) -- ^ Modification function + -> (a -> Rand a) +withProbability p modify x = do + t <- getDouble + if t < p + then modify x + else return x + +-- | Randomly generate Int between 0 and @n@. +-- Should be strictly less than n if n > 1 +-- or zero otherwise. Therefore, getMaxInt 1 +-- should be always 0. +getMaxInt :: Int -- ^ @n@ + -> Rand Int +getMaxInt n = do + r <- getDouble + return $ floor (r * fromIntegral n) + +getNormal :: Rand Double +getNormal = 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 diff --git a/AI/MEP/Run.hs b/AI/MEP/Run.hs new file mode 100644 --- /dev/null +++ b/AI/MEP/Run.hs @@ -0,0 +1,71 @@ +{-# LANGUAGE BangPatterns #-} + +module AI.MEP.Run where + +import qualified Data.Vector as V +import qualified Data.Vector.Mutable as VM +import Data.List ( foldl' ) +import System.IO.Unsafe ( unsafePerformIO ) +import Text.Printf + +import AI.MEP.Types + +-- | Evaluate each subexpression in a chromosome +evaluate :: Num a + => Chromosome a -- ^ Chromosome to evaluate + -> V.Vector a -- ^ Variable values + -> V.Vector a -- ^ Resulting vector of multiple evaluations +evaluate chr vmap = unsafePerformIO $ do + -- Use dynamic programming to evaluate the chromosome + v <- VM.new chrLen + + let -- Gene evaluation function + _f (C c) _ = return c + _f (Var n) _ = return $ vmap V.! n + _f (Op (_, f) i1 i2) v' = do + !r1 <- v' `VM.read` i1 + !r2 <- v' `VM.read` i2 + let !r = f r1 r2 + return r + + -- Chromosome evaluation + go !v' !j = + if j == chrLen + then return () + else do + val <- _f (chr V.! j) v' + VM.write v' j val + go v' (j + 1) + + go v 0 + + V.unsafeFreeze v + where chrLen = V.length chr +{-# SPECIALIZE + evaluate :: Chromosome Double + -> V.Vector Double + -> V.Vector Double #-} + +-- | Generate code for the functions with a single output +generateCode :: Phenotype Double -> String +generateCode (_, chr, i) = concat expr1 ++ expr2 + where + finalI = V.head i + expr1 = map (\k -> _f (chr V.! k) k) [0..finalI - 1] + expr2 = printf "result = %s\n" $ _h (chr V.! finalI) + + _f (C c) _ = "" + _f (Var i) _ = "" + _f op k = printf "v%d = %s\n" k (_h op) + + _h (C c) = show c + _h (Var i) = printf "x%d" i + _h (Op (s, _) i1 i2) = printf "%s %c %s" (_g (chr V.! i1) i1) s (_g (chr V.! i2) i2) + + _g (C c) _ = show c + _g (Var i) _ = printf "x%d" i + _g Op {} k = printf "v%d" k + +-- | Average population loss +avgLoss :: [Phenotype Double] -> Double +avgLoss = uncurry (/). foldl' (\(c, i) (val, _, _) -> (c + val, i + 1)) (0, 0) diff --git a/AI/MEP/Types.hs b/AI/MEP/Types.hs new file mode 100644 --- /dev/null +++ b/AI/MEP/Types.hs @@ -0,0 +1,37 @@ +{- | Provide the basic MEP data structures + -} +{-# LANGUAGE GADTs #-} +module AI.MEP.Types where + +import qualified Data.Vector as V + + +-- Working with lists is not optimal. +-- For instance, a random selection operator +-- such as binaryTournament may look for last +-- elements in the list quite long for big +-- populations. +type Population a = [Chromosome a] + +type Phenotype a = (Double, Chromosome a, V.Vector Int) + +-- | A chromosome is a vector of genes +type Chromosome a = V.Vector (Gene a Int) + +-- | Either a terminal symbol or a three-address code (a function +-- and two pointers) +data Gene a i where + -- Terminal symbol: constant + C :: a -> Gene a i + -- Terminal symbol: variable + Var :: Int -> Gene a i + -- Operation + Op :: F a -> i -> i -> Gene a i + +instance (Show a, Show i) => Show (Gene a i) where + show (C c) = show c + show (Var n) = "v" ++ show n + show (Op (s, _) i1 i2) = show s ++ " " ++ show i1 ++ " " ++ show i2 + +-- | A function and its symbolic representation +type F a = (Char, a -> a -> a) diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,10 @@ +# Changelog for [`hmep` package](http://hackage.haskell.org/package/hmep) + +## 0.0.1 *October 7th 2017* + * Improved demo: trigonometric identities solving example + * Add `avgLoss` to the library + * Fixes: + * Change the crossover probability using Config.p'crossover parameter + +## 0.0.0 *October 6th 2017* + * Initial release diff --git a/MEP.hs b/MEP.hs deleted file mode 100644 --- a/MEP.hs +++ /dev/null @@ -1,39 +0,0 @@ -{- | -Copyright Bogdan Penkovsky (c) 2017 - -= Multiple Expression Programming - --} - -module MEP ( - Chromosome (..) - , Gene - , Population - , Phenotype - , Config (..) - , defaultConfig - , LossFunction - - -- * Genetic algorithm - , initialize - , evaluateGeneration - , evolve - , binaryTournament - , crossover - , mutation3 - , smoothMutation - , newChromosome - - -- * Random - , Rand - , newPureMT - , runRandom - , evalRandom - ) where - -import System.Random.Mersenne.Pure64 ( newPureMT ) - -import MEP.Types -import MEP.Operators -import MEP.Run -import MEP.Random diff --git a/MEP/Operators.hs b/MEP/Operators.hs deleted file mode 100644 --- a/MEP/Operators.hs +++ /dev/null @@ -1,292 +0,0 @@ --- | --- = Genetic operators - -module MEP.Operators ( - Config (..) - , defaultConfig - , LossFunction - , Phenotype - -- * Genetic operators - , initialize - , evaluateGeneration - , evolve - , phenotype - , binaryTournament - , crossover - , mutation3 - , smoothMutation - , newChromosome - ) where - -import Data.Vector ( Vector ) -import qualified Data.Vector as V -import Data.List - ( nub - , sortBy - ) -import Data.Ord ( comparing ) -import qualified Control.Monad as CM - -import MEP.Random -import MEP.Types -import MEP.Run ( evaluate ) - -data Config a = Config - { - p'const :: Double -- ^ Probability of constant generation - , p'var :: Double -- ^ Probability of variable generation. - -- The probability of operator generation is inferred - -- automatically as @1 - p'const - p'var@. - , p'mutation :: Double -- ^ Mutation probability - , p'crossover :: Double -- ^ Crossover probability - - , c'length :: Int -- ^ The chromosome length - , c'popSize :: Int -- ^ A (sub)population size - , c'popN :: Int -- ^ Number of subpopulations (1 or more) [not implemented] - , c'ops :: Vector (F a) -- ^ Functions pool with their symbolic - -- representations - , c'vars :: Int -- ^ The input dimensionality - } - --- | --- @ --- defaultConfig = Config --- { --- p'const = 0.1 --- , p'var = 0.4 --- , p'mutation = 0.1 --- , p'crossover = 0.9 --- --- , c'length = 50 --- , c'popSize = 100 --- , c'popN = 1 --- , c'ops = V.empty -- <-- To be overridden --- , c'vars = 1 --- } --- @ -defaultConfig :: Config Double -defaultConfig = Config - { - p'const = 0.1 - , p'var = 0.4 - , p'mutation = 0.1 - , p'crossover = 0.9 - - , c'length = 50 - , c'popSize = 100 - , c'popN = 1 - , c'ops = V.empty - , c'vars = 1 - } - --- | A function to minimize. --- --- The argument is a vector evaluation function whose input --- is a vector (length @c'vars@) and ouput is --- a vector with a different length @c'length@. --- --- The result is a vector of the best indices --- and a scalar loss value. -type LossFunction a = - ((V.Vector a -> V.Vector a) -> (V.Vector Int, Double)) - --- | Evaluates a chromosome according to the given --- loss function. -phenotype - :: Num a => - LossFunction a - -> Chromosome a - -> Phenotype a -phenotype loss chr = let (is, val) = loss (evaluate chr) - in (val, chr, is) - -type Phenotype a = (Double, Chromosome a, V.Vector Int) - --- | Randomly generate a new population -initialize :: Config Double -> Rand (Population Double) -initialize c@Config { c'popSize = size } = mapM (\_ -> newChromosome c) [1..size] - -evaluateGeneration - :: Num a => - LossFunction a - -> [Chromosome a] - -> [Phenotype a] -evaluateGeneration loss pop = map (phenotype loss) pop - --- | Selection operator that produces the next evaluated population. --- --- Standard algorithm: the best offspring O replaces the worst --- individual W in the current population if O is better than W. -evolve - :: - Config Double - -- ^ Common configuration - -> LossFunction Double - -- ^ Custom loss function - -> (Chromosome Double -> Rand (Chromosome Double)) - -- ^ Mutation - -> (Chromosome Double -> Chromosome Double -> Rand (Chromosome Double, Chromosome Double)) - -- ^ Crossover - -> ([Phenotype Double] -> Rand (Chromosome Double)) - -- ^ A chromosome selection algorithm. Does not need to be random, but may be. - -> [Phenotype Double] - -- ^ Evaluated population - -> Rand [Phenotype Double] - -- ^ New generation -evolve c loss mut cross select phenotypes = do - let pc = p'crossover c - pm = p'mutation c - -- Sort in decreasing @val@ order so that - -- the worst (with the biggest loss) is in the head - sort' = sortBy (comparing (\(val, _, _) -> negate val)) - - ev phen0 _ = do - chr1 <- select phen0 - chr2 <- select phen0 - (of1, of2) <- cross chr1 chr2 - of1' <- withProbability pm mut of1 - of2' <- withProbability pm mut of2 - let r1@(val1, _, _) = phenotype loss of1' - r2@(val2, _, _) = phenotype loss of2' - (worstVal, _, _) = head phen0 - phen' | val1 < worstVal = r1 : tail phen0 - | val2 < worstVal = r2 : tail phen0 - -- No change - | otherwise = phen0 - let phen1 = sort' phen' - return phen1 - - pop' <- CM.foldM ev (sort' phenotypes) [1..c'popSize c `div` 2] - return pop' - --- | Binary tournament selection -binaryTournament :: Ord a => [Phenotype a] -> Rand (Chromosome a) -binaryTournament phen = do - (val1, cand1, _) <- draw $ V.fromList phen - (val2, cand2, _) <- draw $ V.fromList phen - if val1 < val2 - then return cand1 - else return cand2 - --- | Uniform crossover operator -crossover :: - Chromosome a - -> Chromosome a - -> Rand (Chromosome a, Chromosome a) -crossover ca cb = do - r <- V.zipWithM (curry (swap 0.5)) ca cb - return $ V.unzip r - -swap :: Double -> (t, t) -> Rand (t, t) -swap p = withProbability p (\(a, b) -> return (b, a)) - -replaceAt :: Int -> a -> Vector a -> Vector a -replaceAt i gene chr0 = - let (c1, c2) = V.splitAt i chr0 - in c1 V.++ V.singleton gene V.++ V.tail c2 - --- | Mutation operator with up to three mutations per chromosome -mutation3 :: - Config Double - -- ^ Common configuration - -> Chromosome Double - -> Rand (Chromosome Double) -mutation3 c chr = do - -- Subtract 1 to get a non-zero head to - -- replace - is <- nub <$> CM.replicateM k (getMaxInt (chrLen - 1)) - genes <- mapM new' is - let chr' = foldr (uncurry replaceAt) - chr - (zip is genes) - return chr' - where chrLen = V.length chr - k = 3 - new' = new (p'const c) (p'var c) (c'vars c) (c'ops c) - --- | Mutation operator with a fixed mutation probability --- of each gene -smoothMutation - :: - Double - -- ^ Probability of gene mutation - -> Config Double - -- ^ Common configuration - -> Chromosome Double - -> Rand (Chromosome Double) -smoothMutation p c chr = - let new' = new (p'const c) (p'var c) (c'vars c) (c'ops c) - mutate i = withProbability p (\_ -> new' i) - in V.zipWithM mutate (V.enumFromN 0 (V.length chr)) chr - --- | Randomly initialize a new chromosome. --- By definition, the first gene is terminal (a constant --- or a variable). -newChromosome :: - Config Double -- ^ Common configuration - -> Rand (Chromosome Double) -newChromosome c = do - let pConst = p'const c - pVar = p'var c - V.mapM (new pConst pVar (c'vars c) (c'ops c)) $ V.enumFromN 0 (c'length c) - --- | Produce a new random gene -new :: - Double -- ^ Probability to produce a constant - -> Double -- ^ Probability to produce a variable - -> Int -- ^ Number of input variables - -> Vector (F Double) -- ^ Operations vector - -> Int -- ^ Maximal operation index - -> Rand (Gene Double Int) -new p1 p2 vars ops maxIndex = if maxIndex == 0 - -- The head must be a terminal - -- p1' = p1 + (1 - p1 - p2) / 2 = 1/2 + p1/2 - p2/2 - then let p1' = 0.5 * (1 + p1 - p2) - in newTerminal p1' vars - else do - p' <- getDouble - let sel | p' < p1 = newC - | p' < (p1 + p2) = newVar vars - | otherwise = newOp ops maxIndex - sel - -newTerminal :: - Double -- ^ Probability @p@ of a constant generation. - -- @1-p@ will be the probability of a variable generation. - -> Int -- ^ Number of input variables - -> Rand (Gene Double i) -newTerminal p vars = do - p' <- getDouble - if p' < p - then newC - else newVar vars - --- | A randomly generated variable identifier -newVar :: Int -> Rand (Gene a i) -newVar vars = do - var <- draw $ V.enumFromN 0 vars - return $ Var var - --- | A random operation from the operations vector -newOp - :: Vector (F a) - -> Int - -> Rand (Gene a Int) -newOp ops maxIndex = do - op <- draw ops - i1 <- getMaxInt maxIndex - i2 <- getMaxInt maxIndex - return $ Op op i1 i2 - --- | Draw a constant from the normal distribution -newCNormal - :: Double -- ^ Mean - -> Double -- ^ Std deviation - -> Rand (Gene Double i) -newCNormal mu sigma = do - n <- getNormal - return $ C (mu + sigma*n) - --- | Draw a constant from the uniform distribution -newC :: Rand (Gene Double i) -newC = C <$> getDouble diff --git a/MEP/Random.hs b/MEP/Random.hs deleted file mode 100644 --- a/MEP/Random.hs +++ /dev/null @@ -1,55 +0,0 @@ -module MEP.Random - ( - -- * Utilities - draw - , getNormal - , getMaxInt - , withProbability - - -- * Re-exports - , getBool, getInt, getWord, getDouble - , runRandom, evalRandom - , Rand, Random - ) where - -import Control.Monad.Mersenne.Random -import Data.Complex (Complex (..)) -import System.Random -import Data.Vector as V - --- | Randomly draw an element from a vector -draw :: Vector a -> Rand a -draw xs = - Rand $ \g -> let (n, g') = randomR (0, V.length xs - 1) g - r = xs V.! n - in R r g' - --- | Modify value with probability @p@ -withProbability - :: Double -- ^ The probability @p@ - -> (a -> Rand a) -- ^ Modification function - -> (a -> Rand a) -withProbability p modify x = do - t <- getDouble - if t < p - then modify x - else return x - --- | Randomly generate Int between 0 and @n@. --- Should be strictly less than n if n > 1 --- or zero otherwise. Therefore, getMaxInt 1 --- should be always 0. -getMaxInt :: Int -- ^ @n@ - -> Rand Int -getMaxInt n = do - r <- getDouble - return $ floor (r * fromIntegral n) - -getNormal :: Rand Double -getNormal = 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 diff --git a/MEP/Run.hs b/MEP/Run.hs deleted file mode 100644 --- a/MEP/Run.hs +++ /dev/null @@ -1,45 +0,0 @@ -{-# LANGUAGE BangPatterns #-} - -module MEP.Run where - -import qualified Data.Vector as V -import qualified Data.Vector.Mutable as VM -import System.IO.Unsafe ( unsafePerformIO ) - -import MEP.Types - --- | Evaluate each subexpression in a chromosome -evaluate :: Num a - => Chromosome a -- ^ Chromosome to evaluate - -> V.Vector a -- ^ Variable values - -> V.Vector a -- ^ Resulting vector of multiple evaluations -evaluate chr vmap = unsafePerformIO $ do - -- Use dynamic programming to evaluate the chromosome - v <- VM.new chrLen - - let -- Gene evaluation function - _f (C c) _ = return c - _f (Var n) _ = return $ vmap V.! n - _f (Op (_, f) i1 i2) v' = do - !r1 <- v' `VM.read` i1 - !r2 <- v' `VM.read` i2 - let !r = f r1 r2 - return r - - -- Chromosome evaluation - go !v' !j = - if j == chrLen - then return () - else do - val <- _f (chr V.! j) v' - VM.write v' j val - go v' (j + 1) - - go v 0 - - V.unsafeFreeze v - where chrLen = V.length chr -{-# SPECIALIZE - evaluate :: Chromosome Double - -> V.Vector Double - -> V.Vector Double #-} diff --git a/MEP/Types.hs b/MEP/Types.hs deleted file mode 100644 --- a/MEP/Types.hs +++ /dev/null @@ -1,30 +0,0 @@ -{- | Provide the basic MEP data structures - -} -{-# LANGUAGE GADTs #-} -module MEP.Types where - -import qualified Data.Vector as V - - -type Population a = [Chromosome a] - --- | A chromosome is a vector of genes -type Chromosome a = V.Vector (Gene a Int) - --- | Either a terminal symbol or a three-address code (a function --- and two pointers) -data Gene a i where - -- Terminal symbol: constant - C :: a -> Gene a i - -- Terminal symbol: variable - Var :: Int -> Gene a i - -- Operation - Op :: F a -> i -> i -> Gene a i - -instance (Show a, Show i) => Show (Gene a i) where - show (C c) = show c - show (Var n) = "v" ++ show n - show (Op (s, _) i1 i2) = show s ++ " " ++ show i1 ++ " " ++ show i2 - --- | A function and its symbolic representation -type F a = (Char, a -> a -> a) diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Multi Expression Programming -You say, Haskell has not enough machine learning libraries? +You say, not enough Haskell machine learning libraries? Here is yet another one! @@ -8,8 +8,8 @@ There exist many other Genetic Algorithm (GA) Haskell packages. Personally I have used -[simple genetic algorithm](http://hackage.haskell.org/package/moo), -[GA](http://hackage.haskell.org/package/moo), +[simple genetic algorithm](http://hackage.haskell.org/package/simple-genetic-algorithm-mr), +[GA](http://hackage.haskell.org/package/GA), and [moo](http://hackage.haskell.org/package/moo) for quite a long time. The last package was the most preferred, but the other two are also great. @@ -21,8 +21,8 @@ and it would be more efficient to implement MEP from scratch. That is how this package was started. I also wish to say thank you -to the authors of the [moo](http://hackage.haskell.org/package/moo) -GA library, which inspired the present +to the authors of the [moo](http://hackage.haskell.org/package/moo) +GA library, which inspired the present [hmep](http://github.com/masterdezign/hmep) package. ## About MEP @@ -32,3 +32,26 @@ Each gene is featuring [code reuse](https://en.wikipedia.org/wiki/Code_reuse). For more details, please check http://mepx.org/papers.html and https://en.wikipedia.org/wiki/Multi_expression_programming. + +## How to build + +Use [Stack](http://haskellstack.org). + + $ git clone https://github.com/masterdezign/hmep.git && cd hmep + $ stack build --install-ghc + +Now, run the demo to calculate cos^2(x) through sin(x): + + $ stack exec hmep-demo + + Average loss in the initial population 15.268705681244962 + Population 10: average loss 14.709728527360586 + Population 20: average loss 13.497114190675477 + Population 30: average loss 8.953185872653737 + Population 40: average loss 8.953185872653737 + Population 50: average loss 3.3219954564955856e-15 + + Interpreted expression: + v1 = sin x0 + v2 = v1 * v1 + result = 1 - v2 diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -1,5 +1,14 @@ -1. Provide a Show instance for AI.MEP.Types.Chromosome - such that a Haskell code is generated. - -2. Provide a Storable instance for AI.MEP.Types.Gene +1. Provide a Storable instance for AI.MEP.Types.Gene and make the Chromosome a Data.Vector.Storable. + +2. Improve code generation. Features: + a) Removal of dead (unused) expressions + b) Subexpression elimination, e.g. x0 / x0 -> 1 + +3. Improve the demo: provide a CLI interface to work + with external data (using loadMatrix from hmatrix library) + +4. Performance tuning and benchmarking using Criterion package. + Hint: use of matrices featuring O(1) memory access + instead of lists of vectors ([Chromosome a], [Phenotype a]), + might improve the speed of such operators as binaryTournament. diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -1,5 +1,12 @@ module Main where +{- + | = Example application: trigonometry cheating + + Find the trigonometric expression of cos(x) through sin(x) + using our automatic programming method. +-} + import qualified Data.Vector as V import Data.List ( foldl' ) import Control.Monad ( foldM ) @@ -9,15 +16,17 @@ , toList ) -import MEP +import AI.MEP -ops = V.fromList [('*', (*)), ('+', (+)), ('/', (/)), ('-', (-))] +ops = V.fromList [('*', (*)), ('+', (+)), ('/', (/)), ('-', (-)), + ('s', \x _ -> sin x)] config = defaultConfig { c'ops = ops , c'length = 50 } +-- Feel free to change the random number generation seed seed :: Int seed = 3 @@ -25,10 +34,16 @@ randDomain = map (subtract pi. (2*pi *)). toList. randomVector seed Uniform dataset1 :: V.Vector (Double, Double) -dataset1 = V.map (\x -> (x, sin x)) $ V.fromList $ randDomain nSamples +dataset1 = V.map (\x -> (x, function x)) $ V.fromList $ randDomain nSamples where nSamples = 50 + function x = (cos x)^2 -dist x y = abs $ x - y +-- | Absolute value distance between two scalar values +dist :: Double -> Double -> Double +dist x y = if isNaN x || isNaN y + -- Large distance + then 10000 + else abs $ x - y loss :: LossFunction Double loss evalf = (V.singleton i', loss') @@ -54,18 +69,13 @@ :: [Phenotype Double] -> Rand [Phenotype Double] nextGeneration = evolve config loss (mutation3 config) crossover binaryTournament -avgLoss :: [Phenotype Double] -> Double -avgLoss xs = - let (r, len) = foldl' (\(c, i) (val, _, _) -> (c + val, i + 1)) (0, 0) xs - in r / (fromIntegral len) - runIO (pop, g') i = do let (newPop, g2) = foldr (\_ xg -> run xg) (pop, g') [1..generations] putStrLn $ "Population " ++ show (i * generations) ++ ": average loss " ++ show (avgLoss newPop) return (newPop, g2) where run (x, g) = runRandom (nextGeneration x) g - generations = 40 + generations = 5 main :: IO () main = do @@ -74,5 +84,8 @@ popEvaluated = evaluateGeneration loss pop putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated) - (final, _) <- foldM runIO (popEvaluated, g') [1..100] - print $ last final + (final, _) <- foldM runIO (popEvaluated, g') [1..20] + let best = last final + print best + putStrLn "Interpreted expression:" + putStrLn $ generateCode best diff --git a/hmep.cabal b/hmep.cabal --- a/hmep.cabal +++ b/hmep.cabal @@ -1,5 +1,5 @@ name: hmep -version: 0.0.0 +version: 0.0.1 synopsis: HMEP Multi Expression Programming – a genetic programming variant description: A multi expression programming implementation with @@ -14,15 +14,15 @@ copyright: 2017 Bogdan Penkovsky category: AI build-type: Simple -extra-source-files: README.md TODO +extra-source-files: README.md TODO CHANGELOG.md cabal-version: >=1.22 library - exposed-modules: MEP - , MEP.Run - , MEP.Types - other-modules: MEP.Random - , MEP.Operators + exposed-modules: AI.MEP + , AI.MEP.Run + , AI.MEP.Types + other-modules: AI.MEP.Random + , AI.MEP.Operators build-depends: base >= 4.7 && < 5 , containers , monad-mersenne-random diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -3,8 +3,8 @@ , exitFailure ) import qualified Data.Vector as V -import MEP.Types -import MEP.Run +import AI.MEP.Types +import AI.MEP.Run o'mult :: Num a => F a o'mult = ('*', (*))