diff --git a/AI/MEP.hs b/AI/MEP.hs --- a/AI/MEP.hs +++ b/AI/MEP.hs @@ -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 diff --git a/AI/MEP/Operators.hs b/AI/MEP/Operators.hs --- a/AI/MEP/Operators.hs +++ b/AI/MEP/Operators.hs @@ -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) diff --git a/AI/MEP/Random.hs b/AI/MEP/Random.hs --- a/AI/MEP/Random.hs +++ b/AI/MEP/Random.hs @@ -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 diff --git a/AI/MEP/Run.hs b/AI/MEP/Run.hs --- a/AI/MEP/Run.hs +++ b/AI/MEP/Run.hs @@ -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) diff --git a/AI/MEP/Types.hs b/AI/MEP/Types.hs --- a/AI/MEP/Types.hs +++ b/AI/MEP/Types.hs @@ -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) diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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) diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -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 diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -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:" diff --git a/cabal.config b/cabal.config new file mode 100644 --- /dev/null +++ b/cabal.config @@ -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 diff --git a/hmep.cabal b/hmep.cabal --- a/hmep.cabal +++ b/hmep.cabal @@ -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