packages feed

hmep 0.0.0 → 0.0.1

raw patch · 16 files changed

+623/−491 lines, 16 files

Files

+ AI/MEP.hs view
@@ -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
+ AI/MEP/Operators.hs view
@@ -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
+ AI/MEP/Random.hs view
@@ -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
+ AI/MEP/Run.hs view
@@ -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)
+ AI/MEP/Types.hs view
@@ -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)
+ CHANGELOG.md view
@@ -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
− MEP.hs 
@@ -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
− MEP/Operators.hs 
@@ -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
− MEP/Random.hs 
@@ -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
− MEP/Run.hs 
@@ -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 #-}
− MEP/Types.hs 
@@ -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)
README.md view
@@ -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
TODO view
@@ -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.
app/Main.hs view
@@ -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
hmep.cabal view
@@ -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
test/Spec.hs view
@@ -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 = ('*', (*))