hmep 0.0.1 → 0.1.0
raw patch · 11 files changed
+233/−173 lines, 11 filesdep +mwc-randomdep +primitivedep +probabledep −containersdep −hmatrixdep −mersenne-random-pure64
Dependencies added: mwc-random, primitive, probable, statistics
Dependencies removed: containers, hmatrix, mersenne-random-pure64, monad-mersenne-random, random
Files
- AI/MEP.hs +5/−9
- AI/MEP/Operators.hs +44/−52
- AI/MEP/Random.hs +43/−40
- AI/MEP/Run.hs +1/−1
- AI/MEP/Types.hs +16/−10
- CHANGELOG.md +6/−0
- README.md +8/−0
- TODO +1/−1
- app/Main.hs +53/−49
- cabal.config +48/−0
- hmep.cabal +8/−11
AI/MEP.hs view
@@ -28,8 +28,8 @@ 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 average value of 3.3e-15 is close to zero, indicating that+ the exact expression was found! The produced output was: @@ -42,7 +42,7 @@ From here we can infer that @- cos^2(x) = 1 - v2 = 1 - v1 * v1 = 1 - sin^2(x)+ cos^2(x) = 1 - v2 = 1 - v1 * v1 = 1 - sin^2(x). @ Sweet!@@ -73,13 +73,9 @@ , generateCode -- * Random- , Rand- , newPureMT- , runRandom- , evalRandom+ , RandT+ , runRandIO ) where--import System.Random.Mersenne.Pure64 ( newPureMT ) import AI.MEP.Types import AI.MEP.Operators
AI/MEP/Operators.hs view
@@ -25,6 +25,7 @@ ) import Data.Ord ( comparing ) import qualified Control.Monad as CM+import Control.Monad.Primitive ( PrimMonad ) import AI.MEP.Random import AI.MEP.Types@@ -100,14 +101,14 @@ in (val, chr, is) -- | Randomly generate a new population-initialize :: Config Double -> Rand (Population Double)+initialize :: PrimMonad m => Config Double -> RandT m (Population Double) initialize c@Config { c'popSize = size } = mapM (\_ -> newChromosome c) [1..size] evaluateGeneration :: Num a => LossFunction a- -> [Chromosome a]- -> [Phenotype a]+ -> Population a+ -> Generation a evaluateGeneration loss = map (phenotype loss) -- | Selection operator that produces the next evaluated population.@@ -115,20 +116,20 @@ -- Standard algorithm: the best offspring O replaces the worst -- individual W in the current population if O is better than W. evolve- ::+ :: PrimMonad m => Config Double -- ^ Common configuration -> LossFunction Double -- ^ Custom loss function- -> (Chromosome Double -> Rand (Chromosome Double))+ -> (Chromosome Double -> RandT m (Chromosome Double)) -- ^ Mutation- -> (Chromosome Double -> Chromosome Double -> Rand (Chromosome Double, Chromosome Double))+ -> (Chromosome Double -> Chromosome Double -> RandT m (Chromosome Double, Chromosome Double)) -- ^ Crossover- -> ([Phenotype Double] -> Rand (Chromosome Double))+ -> (Generation Double -> RandT m (Chromosome Double)) -- ^ A chromosome selection algorithm. Does not need to be random, but may be.- -> [Phenotype Double]+ -> Generation Double -- ^ Evaluated population- -> Rand [Phenotype Double]+ -> RandT m (Generation Double) -- ^ New generation evolve c loss mut cross select phenotypes = do let pc = p'crossover c@@ -156,24 +157,24 @@ CM.foldM ev (sort' phenotypes) [1..c'popSize c `div` 2] -- | Binary tournament selection-binaryTournament :: Ord a => [Phenotype a] -> Rand (Chromosome a)+binaryTournament :: (PrimMonad m, Ord a) => Generation a -> RandT m (Chromosome a) binaryTournament phen = do- (val1, cand1, _) <- draw $ V.fromList phen- (val2, cand2, _) <- draw $ V.fromList phen+ (val1, cand1, _) <- drawFrom $ V.fromList phen+ (val2, cand2, _) <- drawFrom $ V.fromList phen if val1 < val2 then return cand1 else return cand2 -- | Uniform crossover operator-crossover ::+crossover :: PrimMonad m => Chromosome a -> Chromosome a- -> Rand (Chromosome a, Chromosome a)+ -> RandT m (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 :: PrimMonad m => Double -> (t, t) -> RandT m (t, t) swap p = withProbability p (\(a, b) -> return (b, a)) replaceAt :: Int -> a -> Vector a -> Vector a@@ -182,15 +183,15 @@ in c1 V.++ V.singleton gene V.++ V.tail c2 -- | Mutation operator with up to three mutations per chromosome-mutation3 ::+mutation3 :: PrimMonad m => Config Double -- ^ Common configuration -> Chromosome Double- -> Rand (Chromosome Double)+ -> RandT m (Chromosome Double) mutation3 c chr = do -- Subtract 1 to get a non-zero head to -- replace- is <- nub <$> CM.replicateM k (getMaxInt (chrLen - 1))+ is <- nub <$> CM.replicateM k (uniformIn_ (0, chrLen - 1)) genes <- mapM new' is let chr' = foldr (uncurry replaceAt) chr@@ -203,13 +204,13 @@ -- | Mutation operator with a fixed mutation probability -- of each gene smoothMutation- ::+ :: PrimMonad m => Double -- ^ Probability of gene mutation -> Config Double -- ^ Common configuration -> Chromosome Double- -> Rand (Chromosome Double)+ -> RandT m (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)@@ -218,71 +219,62 @@ -- | Randomly initialize a new chromosome. -- By definition, the first gene is terminal (a constant -- or a variable).-newChromosome ::+newChromosome :: PrimMonad m => Config Double -- ^ Common configuration- -> Rand (Chromosome Double)+ -> RandT m (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+new :: PrimMonad m =>+ 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)+ -> RandT m (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+ p' <- double let sel | p' < p1 = newC | p' < (p1 + p2) = newVar vars | otherwise = newOp ops maxIndex sel -newTerminal ::- Double -- ^ Probability @p@ of a constant generation.+newTerminal :: (PrimMonad m, Floating a, Variate a) =>+ 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)+ -> RandT m (Gene a i) newTerminal p vars = do- p' <- getDouble+ p' <- double if p' < p then newC else newVar vars -- | A randomly generated variable identifier-newVar :: Int -> Rand (Gene a i)+newVar :: PrimMonad m => Int -> RandT m (Gene a i) newVar vars = do- var <- draw $ V.enumFromN 0 vars+ var <- drawFrom $ V.enumFromN 0 vars return $ Var var -- | A random operation from the operations vector-newOp- :: Vector (F a)+newOp :: PrimMonad m =>+ Vector (F a) -> Int- -> Rand (Gene a Int)+ -> RandT m (Gene a Int) newOp ops maxIndex = do- op <- draw ops- i1 <- getMaxInt maxIndex- i2 <- getMaxInt maxIndex+ op <- drawFrom ops+ i1 <- uniformIn_ (0, maxIndex)+ i2 <- uniformIn_ (0, 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+-- | Draw a constant from the uniform distribution within @(-0.5, 0.5]@+newC :: (PrimMonad m, Floating a, Variate a) => RandT m (Gene a i)+newC = C <$> uniformIn (-0.5, 0.5)
AI/MEP/Random.hs view
@@ -1,55 +1,58 @@ module AI.MEP.Random ( -- * Utilities- draw- , getNormal- , getMaxInt+ drawFrom+ , double_+ , uniformIn_ , withProbability+ , runRandIO - -- * Re-exports- , getBool, getInt, getWord, getDouble- , runRandom, evalRandom- , Rand, Random+ -- * Re-exports from Math.Probable.Random+ , RandT+ , double+ , vectorOf+ , vectorOfVariate+ , uniformIn++ -- * Re-exports from System.Random.MWC+ , Variate ) where -import Control.Monad.Mersenne.Random-import Data.Complex (Complex (..))-import System.Random+import Math.Probable.Random+import System.Random.MWC ( Variate )+import Control.Monad.Primitive ( PrimMonad ) import Data.Vector as V +runRandIO :: RandT IO a -> IO a+runRandIO = mwc+ -- | 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'+drawFrom :: PrimMonad m => Vector a -> RandT m a+drawFrom vec = do+ n <- uniformIn (0, V.length vec - 1)+ return $ vec V.! n --- | Modify value with probability @p@-withProbability- :: Double -- ^ The probability @p@- -> (a -> Rand a) -- ^ Modification function- -> (a -> Rand a)+-- | Similar to uniformIn, but using range+-- @[a, b)@ instead of @[a, b]@ and only for integral types.+uniformIn_ :: (PrimMonad m, Variate a, Integral a) => (a, a) -> RandT m a+uniformIn_ (a, b) = uniformIn (a, b - 1)++-- | Returns a double value from the range of @[0, 1)@.+-- If there is no specific reason, then prefer double @(0, 1]@.+double_ :: PrimMonad m => RandT m Double+double_ = (subtract magicC) <$> double+ where+ -- Change the range (0, 1] to (0, 1].+ -- http://hackage.haskell.org/package/mwc-random-0.13.6.0/docs/System-Random-MWC.html#v:uniform+ magicC = 2**(-53)++-- | Modify a value with the probability @p@+withProbability :: PrimMonad m =>+ Double -- ^ The probability @p@+ -> (a -> RandT m a) -- ^ Modification function+ -> (a -> RandT m a) withProbability p modify x = do- t <- getDouble+ t <- double 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
@@ -67,5 +67,5 @@ _g Op {} k = printf "v%d" k -- | Average population loss-avgLoss :: [Phenotype Double] -> Double+avgLoss :: Generation Double -> Double avgLoss = uncurry (/). foldl' (\(c, i) (val, _, _) -> (c + val, i + 1)) (0, 0)
AI/MEP/Types.hs view
@@ -1,4 +1,5 @@-{- | Provide the basic MEP data structures+{- |+ = Core MEP data structures -} {-# LANGUAGE GADTs #-} module AI.MEP.Types where@@ -6,15 +7,6 @@ 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) @@ -35,3 +27,17 @@ -- | A function and its symbolic representation type F a = (Char, a -> a -> a)++-- 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.+-- | List of chromosomes+type Population a = [Chromosome a]++-- | Evaluated population+type Generation a = [Phenotype a]++-- | Loss value, chromosome, and the best expression indices vector+type Phenotype a = (Double, Chromosome a, V.Vector Int)
CHANGELOG.md view
@@ -1,5 +1,11 @@ # Changelog for [`hmep` package](http://hackage.haskell.org/package/hmep) +## 0.1.0 *October 8th 2017*+ * Breaking changes:+ drop [monad-mersenne-random](http://hackage.haskell.org/package/monad-mersenne-random)+ which doesn't build with newer version of GHC. Instead, depend on the+ [probable](http://hackage.haskell.org/package/probable) package.+ ## 0.0.1 *October 7th 2017* * Improved demo: trigonometric identities solving example * Add `avgLoss` to the library
README.md view
@@ -4,6 +4,7 @@ Here is yet another one! + ## History There exist many other Genetic Algorithm (GA) Haskell packages.@@ -25,6 +26,7 @@ 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@@ -33,6 +35,7 @@ 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).@@ -55,3 +58,8 @@ v1 = sin x0 v2 = v1 * v1 result = 1 - v2+++## Authors++This library is written and maintained by [Bogdan Penkovsky](http://penkovsky.com)
TODO view
@@ -6,7 +6,7 @@ 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)+ with external data 4. Performance tuning and benchmarking using Criterion package. Hint: use of matrices featuring O(1) memory access
app/Main.hs view
@@ -10,34 +10,26 @@ import qualified Data.Vector as V import Data.List ( foldl' ) import Control.Monad ( foldM )-import Numeric.LinearAlgebra- ( randomVector- , RandDist( Uniform )- , toList+import Math.Probable.Random -- From `probable` package+ ( vectorOf+ , double ) import AI.MEP -ops = V.fromList [('*', (*)), ('+', (+)), ('/', (/)), ('-', (-)),- ('s', \x _ -> sin x)]- config = defaultConfig {- c'ops = ops+ -- Functions available to genetically produced programs+ c'ops = V.fromList [+ ('*', (*)),+ ('+', (+)),+ ('/', (/)),+ ('-', (-)),+ ('s', \x _ -> sin x)+ ]+ -- Chromosome length , c'length = 50 } --- Feel free to change the random number generation seed-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, function x)) $ V.fromList $ randDomain nSamples- where nSamples = 50- function x = (cos x)^2- -- | Absolute value distance between two scalar values dist :: Double -> Double -> Double dist x y = if isNaN x || isNaN y@@ -45,19 +37,6 @@ then 10000 else 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@@ -65,26 +44,51 @@ len = V.length $ head xss base = V.replicate len 0 -nextGeneration- :: [Phenotype Double] -> Rand [Phenotype Double]-nextGeneration = evolve config loss (mutation3 config) crossover binaryTournament--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 = 5- main :: IO () main = do- g <- newPureMT- let (pop, g') = runRandom (initialize config) g- popEvaluated = evaluateGeneration loss pop+ -- A vector of 50 random numbers between 0 and 1 (including 1)+ xs <- runRandIO (vectorOf 50 double)++ -- Scale the values to the interval of (-pi, pi]+ let xs' = V.map ((2*pi *). subtract 0.5) xs+ -- Target function f to approximate+ function x = (cos x)^2+ -- Pairs (x, f(x))+ dataset = V.map (\x -> (x, function x)) xs'++ -- Randomly create a population of chromosomes+ pop <- runRandIO $ initialize config++ let -- The loss function which depends on the dataset+ loss evalf = (V.singleton i', loss')+ where+ (xs, ys) = unzip $ V.toList dataset+ -- 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'++ -- Evaluate the initial population+ let popEvaluated = evaluateGeneration loss pop+ putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated) - (final, _) <- foldM runIO (popEvaluated, g') [1..20]+ -- Declare how to produce the new generation+ let nextGeneration = evolve config loss (mutation3 config) crossover binaryTournament++ -- Specify the I/O loop, which logs every 5 generation+ runIO pop i = do+ newPop <- runRandIO $ foldM (\xg _ -> nextGeneration xg) pop [1..generations]+ putStrLn $ "Population " ++ show (i * generations) ++ ": average loss " ++ show (avgLoss newPop)+ return newPop+ where generations = 5++ -- Final generation+ final <- foldM runIO popEvaluated [1..20] let best = last final print best putStrLn "Interpreted expression:"
+ cabal.config view
@@ -0,0 +1,48 @@+constraints: base ==4.8.1.0,+ rts ==1.0,+ ghc-prim ==0.4.0.0,+ integer-gmp ==1.0.0.0,+ mwc-random ==0.13.6.0,+ math-functions ==0.2.1.0,+ deepseq ==1.4.1.1,+ array ==0.5.1.0,+ primitive ==0.6.2.0,+ transformers ==0.4.2.0,+ vector ==0.12.0.1,+ semigroups ==0.18.3,+ binary ==0.7.5.0,+ bytestring ==0.10.6.0,+ containers ==0.5.6.2,+ hashable ==1.2.6.1,+ text ==1.2.2.2,+ tagged ==0.8.5,+ template-haskell ==2.10.0.0,+ pretty ==1.1.2.0,+ transformers-compat ==0.5.1.4,+ unordered-containers ==0.2.8.0,+ vector-th-unbox ==0.2.1.6,+ time ==1.5.0.1,+ probable ==0.1.2,+ mtl ==2.2.1,+ statistics ==0.13.3.0,+ aeson ==1.2.2.0,+ attoparsec ==0.13.2.0,+ fail ==4.9.0.0,+ scientific ==0.3.5.2,+ integer-logarithms ==1.0.2,+ base-compat ==0.9.3,+ unix ==2.7.1.0,+ dlist ==0.8.0.3,+ th-abstraction ==0.2.6.0,+ time-locale-compat ==0.1.1.3,+ uuid-types ==1.0.3,+ random ==1.1,+ erf ==2.0.0.0,+ monad-par ==0.3.4.8,+ abstract-deque ==0.3,+ abstract-par ==0.3.3,+ monad-par-extras ==0.3.3,+ cereal ==0.5.4.0,+ parallel ==3.2.1.1,+ vector-algorithms ==0.7.0.1,+ vector-binary-instances ==0.2.3.5
hmep.cabal view
@@ -1,5 +1,5 @@ name: hmep-version: 0.0.1+version: 0.1.0 synopsis: HMEP Multi Expression Programming – a genetic programming variant description: A multi expression programming implementation with@@ -14,7 +14,7 @@ copyright: 2017 Bogdan Penkovsky category: AI build-type: Simple-extra-source-files: README.md TODO CHANGELOG.md+extra-source-files: README.md TODO CHANGELOG.md cabal.config cabal-version: >=1.22 library@@ -24,10 +24,10 @@ other-modules: AI.MEP.Random , AI.MEP.Operators build-depends: base >= 4.7 && < 5- , containers- , monad-mersenne-random- , mersenne-random-pure64- , random+ , mwc-random+ , primitive+ , probable+ , statistics >= 0.10 && < 0.14 , vector default-language: Haskell2010 @@ -36,10 +36,8 @@ main-is: Main.hs ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends: base- , containers- , hmatrix- , mersenne-random-pure64- , monad-mersenne-random+ , probable+ , statistics >= 0.10 && < 0.14 , vector , hmep default-language: Haskell2010@@ -49,7 +47,6 @@ hs-source-dirs: test main-is: Spec.hs build-depends: base- , containers , HUnit , vector , hmep