hmep (empty) → 0.0.0

raw patch · 12 files changed

+704/−0 lines, 12 filesdep +HUnitdep +basedep +containerssetup-changed

Dependencies added: HUnit, base, containers, hmatrix, hmep, mersenne-random-pure64, monad-mersenne-random, random, vector

Files

• LICENSE +30/−0
• MEP.hs +39/−0
• MEP/Operators.hs +292/−0
• MEP/Random.hs +55/−0
• MEP/Run.hs +45/−0
• MEP/Types.hs +30/−0
• README.md +34/−0
• Setup.hs +2/−0
• TODO +5/−0
• app/Main.hs +78/−0
• hmep.cabal +61/−0
• test/Spec.hs +33/−0

@@ -0,0 +1,30 @@+Copyright Bogdan Penkovsky (c) 2017++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Author name here nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

@@ -0,0 +1,292 @@+-- |+-- = 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

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

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

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

@@ -0,0 +1,34 @@+# Multi Expression Programming++You say, Haskell has not enough machine learning libraries?++Here is yet another one!++## History++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),+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.++However, when I came up with this+[MEP paper](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.5.4352&rep=rep1&type=pdf),+to my surprise there was no MEP realization in Haskell.+Soon I realized that existing GA packages are limited,+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 +[hmep](http://github.com/masterdezign/hmep) package.++## About MEP++Multi Expression Programming is a genetic programming variant encoding multiple+solutions in the same chromosome. A chromosome is a computer program.+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.

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

@@ -0,0 +1,5 @@+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+   and make the Chromosome a Data.Vector.Storable.

@@ -0,0 +1,78 @@+module Main where++import qualified Data.Vector as V+import           Data.List ( foldl' )+import           Control.Monad ( foldM )+import           Numeric.LinearAlgebra+                 ( randomVector+                 , RandDist( Uniform )+                 , toList+                 )++import MEP++ops = V.fromList [('*', (*)), ('+', (+)), ('/', (/)), ('-', (-))]++config = defaultConfig {+  c'ops = ops+  , c'length = 50+  }++seed :: Int+seed = 3++randDomain :: Int -> [Double]+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+  where nSamples = 50++dist x y = abs $ x - y++loss :: LossFunction Double+loss evalf = (V.singleton i', loss')+  where+    (xs, ys) = unzip $ V.toList dataset1+    -- Distances resulting from multiple expression evaluation+    dss = zipWith (\x y -> V.map (dist y). evalf. V.singleton $ x) xs ys+    -- Cumulative distances for each index+    dcumul = sum' dss+    -- Select index minimizing cumulative distances+    i' = V.minIndex dcumul+    -- The loss value with respect to the index of the best expression+    loss' = dcumul V.! i'++-- Could be optimized+sum' :: Num a => [V.Vector a] -> V.Vector a+sum' xss = foldl' (V.zipWith (+)) base xss+  where+    len = V.length $ head xss+    base = V.replicate len 0++nextGeneration+  :: [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++main :: IO ()+main = do+  g <- newPureMT+  let (pop, g') = runRandom (initialize config) g+      popEvaluated = evaluateGeneration loss pop+  putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated)++  (final, _) <- foldM runIO (popEvaluated, g') [1..100]+  print $ last final

@@ -0,0 +1,61 @@+name:                hmep+version:             0.0.0+synopsis:            HMEP Multi Expression Programming –+                     a genetic programming variant+description:         A multi expression programming implementation with+                     focus on speed.+                     .+                     https://en.wikipedia.org/wiki/Multi_expression_programming+homepage:            https://github.com/masterdezign/hmep#readme+license:             BSD3+license-file:        LICENSE+author:              Bogdan Penkovsky+maintainer:          dev at penkovsky dot com+copyright:           2017 Bogdan Penkovsky+category:            AI+build-type:          Simple+extra-source-files:  README.md TODO+cabal-version:       >=1.22++library+  exposed-modules:     MEP+                     , MEP.Run+                     , MEP.Types+  other-modules:       MEP.Random+                     , MEP.Operators+  build-depends:       base >= 4.7 && < 5+                     , containers+                     , monad-mersenne-random+                     , mersenne-random-pure64+                     , random+                     , vector+  default-language:    Haskell2010++executable hmep-demo+  hs-source-dirs:      app+  main-is:             Main.hs+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  build-depends:       base+                     , containers+                     , hmatrix+                     , mersenne-random-pure64+                     , monad-mersenne-random+                     , vector+                     , hmep+  default-language:    Haskell2010++test-suite hmep-test+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             Spec.hs+  build-depends:       base+                     , containers+                     , HUnit+                     , vector+                     , hmep+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  default-language:    Haskell2010++source-repository head+  type:     git+  location: https://github.com/masterdezign/hmep

@@ -0,0 +1,33 @@+import           Test.HUnit+import           System.Exit ( exitSuccess+                             , exitFailure )+import qualified Data.Vector as V++import           MEP.Types+import           MEP.Run++o'mult :: Num a => F a+o'mult = ('*', (*))+{-# SPECIALIZE o'mult :: F Double #-}++-- Encodes x, x^2, x^4, x^8+pow8Int :: Chromosome Int+pow8Int = V.fromList [Var 0, Op o'mult 0 0, Op o'mult 1 1, Op o'mult 2 2]++pow8 :: Chromosome Double+pow8 = V.fromList [Var 0, Op o'mult 0 0, Op o'mult 1 1, Op o'mult 2 2]++testEvaluate = test [+  "2^8" ~: V.fromList [2, 4, 16, 256] ~=? evaluate pow8Int (V.singleton (2 :: Int)),++  "2.5^8" ~: V.fromList [2.5,6.25,39.0625,1525.87890625] ~=? evaluate pow8 (V.singleton (2.5 :: Double))+  ]++allTests = TestList [ testEvaluate ]++main :: IO ()+main = do+  result <- runTestTT allTests+  if (errors result + failures result) > 0+    then exitFailure+    else exitSuccess