diff --git a/AI/MEP.hs b/AI/MEP.hs --- a/AI/MEP.hs +++ b/AI/MEP.hs @@ -60,8 +60,11 @@ -- * Genetic algorithm , initialize - , evaluateGeneration + , evaluatePopulation + , regressionLoss1 , avgLoss + , best + , worst , evolve , binaryTournament , crossover diff --git a/AI/MEP/Operators.hs b/AI/MEP/Operators.hs --- a/AI/MEP/Operators.hs +++ b/AI/MEP/Operators.hs @@ -7,7 +7,9 @@ , LossFunction -- * Genetic operators , initialize - , evaluateGeneration + , evaluatePopulation + , best + , worst , evolve , phenotype , binaryTournament @@ -29,8 +31,9 @@ import AI.MEP.Random import AI.MEP.Types -import AI.MEP.Run ( evaluate ) +import AI.MEP.Run ( evaluateChromosome ) +-- | MEP configuration data Config a = Config { p'const :: Double -- ^ Probability of constant generation @@ -97,20 +100,29 @@ LossFunction a -> Chromosome a -> Phenotype a -phenotype loss chr = let (is, val) = loss (evaluate chr) +phenotype loss chr = let (is, val) = loss (evaluateChromosome chr) in (val, chr, is) -- | Randomly generate a new population initialize :: PrimMonad m => Config Double -> RandT m (Population Double) initialize c@Config { c'popSize = size } = mapM (\_ -> newChromosome c) [1..size] -evaluateGeneration +-- | Using 'LossFunction', find how fit is each chromosome in the population +evaluatePopulation :: Num a => LossFunction a -> Population a -> Generation a -evaluateGeneration loss = map (phenotype loss) +evaluatePopulation loss = map (phenotype loss) +-- | The best phenotype in the generation +best :: Generation a -> Phenotype a +best = last + +-- | The worst phenotype in the generation +worst :: Generation a -> Phenotype a +worst = head + -- | Selection operator that produces the next evaluated population. -- -- Standard algorithm: the best offspring O replaces the worst @@ -222,10 +234,9 @@ newChromosome :: PrimMonad m => Config Double -- ^ Common configuration -> 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) +newChromosome c = + V.mapM new' $ V.enumFromN 0 (c'length c) + where new' = new (p'const c) (p'var c) (c'vars c) (c'ops c) -- | Produce a new random gene new :: PrimMonad m => @@ -260,9 +271,7 @@ -- | A randomly generated variable identifier newVar :: PrimMonad m => Int -> RandT m (Gene a i) -newVar vars = do - var <- drawFrom $ V.enumFromN 0 vars - return $ Var var +newVar vars = Var <$> uniformIn_ (0, vars) -- | A random operation from the operations vector newOp :: PrimMonad m => diff --git a/AI/MEP/Random.hs b/AI/MEP/Random.hs --- a/AI/MEP/Random.hs +++ b/AI/MEP/Random.hs @@ -1,3 +1,6 @@ +-- | +-- = Random helpers + module AI.MEP.Random ( -- * Utilities @@ -23,6 +26,11 @@ import Control.Monad.Primitive ( PrimMonad ) import Data.Vector as V +-- | Alias for @mwc@: +-- Take a 'RandT' value and run it in 'IO', generating all the random values described by the 'RandT'. +-- +-- It initializes the random number generator. For performance reasons, it is +-- recommended to minimize the number of calls to 'runRandIO'. runRandIO :: RandT IO a -> IO a runRandIO = mwc @@ -40,9 +48,9 @@ -- | 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 +double_ = subtract magicC <$> double where - -- Change the range (0, 1] to (0, 1]. + -- 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) diff --git a/AI/MEP/Run.hs b/AI/MEP/Run.hs --- a/AI/MEP/Run.hs +++ b/AI/MEP/Run.hs @@ -1,21 +1,33 @@ +-- | +-- = Various utilities for running MEP algorithm + {-# LANGUAGE BangPatterns #-} -module AI.MEP.Run where +module AI.MEP.Run ( + generateCode + , evaluateChromosome + , regressionLoss1 + , avgLoss + ) where import qualified Data.Vector as V import qualified Data.Vector.Mutable as VM -import Data.List ( foldl' ) +import Data.List ( + foldl' + , nub + , sort + ) import System.IO.Unsafe ( unsafePerformIO ) import Text.Printf import AI.MEP.Types -- | Evaluate each subexpression in a chromosome -evaluate :: Num a +evaluateChromosome :: 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 +evaluateChromosome chr vmap = unsafePerformIO $ do -- Use dynamic programming to evaluate the chromosome v <- VM.new chrLen @@ -42,7 +54,7 @@ V.unsafeFreeze v where chrLen = V.length chr {-# SPECIALIZE - evaluate :: Chromosome Double + evaluateChromosome :: Chromosome Double -> V.Vector Double -> V.Vector Double #-} @@ -50,9 +62,14 @@ generateCode :: Phenotype Double -> String generateCode (_, chr, i) = concat expr1 ++ expr2 where + -- A part of chromosome that is used (all genes ahead the `finalI` + -- and the gene pointed by the `finalI`) + chr' = V.slice 0 (finalI + 1) chr + last' = chr' V.! finalI + finalI = V.head i - expr1 = map (\k -> _f (chr V.! k) k) [0..finalI - 1] - expr2 = printf "result = %s\n" $ _h (chr V.! finalI) + expr1 = map (\k -> _f (chr' V.! k) k). sort. nub $ _usedGeneIx chr' + expr2 = printf "result = %s\n" $ _h last' _f (C c) _ = "" _f (Var i) _ = "" @@ -60,11 +77,80 @@ _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) + _h (Op (s, _) i1 i2) = if isInfix s + then printf "%s %c %s" g1 s g2 + else printf "%c %s %s" s g1 g2 + where g1 = _g (chr' V.! i1) i1 + g2 = _g (chr' V.! i2) i2 _g (C c) _ = show c _g (Var i) _ = printf "x%d" i _g Op {} k = printf "v%d" k + + -- Very naive infix operator check. No problem for single-character + -- ASCII operator representations. Otherwise, please improve. + isInfix x = x `notElem` (['a'..'z'] ++ ['A'..'Z'] ++ ['0'..'9']) + +-- Active genes in case of a chromosome representing a single-output function. +-- Can be generalized to multiple outputs by several calls +-- changing `lastPos` as an argument. +_usedGeneIx :: Chromosome a -> [Int] +_usedGeneIx chr = foldl' _g base $ zip pos $ map (chr V.!) pos + where + -- Position indices + pos = [lastPos - 1,lastPos - 2..0] + + _g xs (i, Op _ i1 i2) = if i `elem` xs + -- Next expressions depend on these + then i1: i2: xs + -- Dead gene, skip + else xs + _g xs _ = xs -- Terminal symbol, already counted + + base = case last' of + (Op _ i1 i2) -> [i1, i2] + _ -> [] -- Sadly, a terminal symbol + + last' = chr V.! lastPos + lastPos = V.length chr - 1 + +-- | Loss function for regression problems with +-- one input and one output. +-- Not normalized with respect to the dataset size. +regressionLoss1 + :: (Num result, Ord result) => + (b -> b -> result) -- ^ Distance function + -> [(a, b)] -- ^ Dataset + -> (V.Vector a -> V.Vector b) + -- ^ Chromosome evaluation function (partially applied 'evaluate') + -> (V.Vector Int, result) +regressionLoss1 dist dataset evalf = (V.singleton i', loss') + where + (xs, ys) = unzip 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' +{-# SPECIALIZE + regressionLoss1 + :: + (Double -> Double -> Double) + -> [(Double, Double)] + -> (V.Vector Double -> V.Vector Double) + -> (V.Vector Int, Double) + #-} + +-- 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 +{-# SPECIALIZE sum' :: [V.Vector Double] -> V.Vector Double #-} -- | Average population loss avgLoss :: Generation Double -> Double diff --git a/AI/MEP/Types.hs b/AI/MEP/Types.hs --- a/AI/MEP/Types.hs +++ b/AI/MEP/Types.hs @@ -20,6 +20,14 @@ -- Operation Op :: F a -> i -> i -> Gene a i +-- | 'Eq' instance for 'Gene' +instance (Eq a, Eq i) => Eq (Gene a i) where + (C a) == (C b) = a == b + (Var a) == (Var b) = a == b + (Op (s1,_) i1 i2) == (Op (s2,_) j1 j2) = s1 == s2 && i1 == j1 && i2 == j2 + _ == _ = False + +-- | 'Show' instance for 'Gene' instance (Show a, Show i) => Show (Gene a i) where show (C c) = show c show (Var n) = "v" ++ show n 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.1 *October 13th 2017* + * Improve code generation + * Add `regressionLoss1` for 1D functions to the library + * Add helpers `best`, `worst` working on Generation + * Improve examples and documentation + ## 0.1.0 *October 8th 2017* * Breaking changes: drop [monad-mersenne-random](http://hackage.haskell.org/package/monad-mersenne-random) diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # Multi Expression Programming +[![Build Status](https://travis-ci.org/masterdezign/hmep.svg?branch=master)](https://travis-ci.org/masterdezign/hmep) + You say, not enough Haskell machine learning libraries? Here is yet another one! @@ -17,7 +19,7 @@ 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. +to my surprise there was no MEP implementation in Haskell. Soon I realized that existing GA packages are limited, and it would be more efficient to implement MEP from scratch. @@ -32,10 +34,47 @@ 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). + +How **MEP is different** from other genetic programming (GP) methods? +Consider a classical example of tree-based GP. +The number of nodes to encode `x^N` +using a binary tree is `2N-1`. +With MEP encoding, however, redundancies can be dramatically +diminished so that the +[shortest chromosome](https://github.com/masterdezign/hmep/blob/cd7b4976800d6c23ce5ebbe67f5ab5c9076229b9/test/Spec.hs#L18) +that encodes the same expression has only `N/2` nodes! +That often results in significantly reduced computational costs +when evaluating MEP chromosomes. Moreover, **all** the intermediate +solutions such as `x^(N/2)`, `x^(N/4)`, etc. are provided by the +chromosome as well. + For more details, please check http://mepx.org/papers.html and https://en.wikipedia.org/wiki/Multi_expression_programming. +### MEP in open source + * By [Mihai Oltean](http://github.com/mepx), C++ + * By [Mark Chenoweth](https://github.com/markcheno/go-mep), Go + * [Current project](https://github.com/masterdezign/hmep), Haskell + +### The `hmep` Features + + * **Works out of the box**. You may use one of the elaborated + [examples](https://github.com/masterdezign/hmep/blob/master/app/) + to quickly tailor to your needs. + * **Flexibility**. The [`hmep` package](https://github.com/masterdezign/hmep/) + provides adjustable and composable building blocks such as + [selection](https://hackage.haskell.org/package/hmep-0.1.0/docs/src/AI-MEP-Operators.html#binaryTournament), + [mutation](https://hackage.haskell.org/package/hmep-0.1.0/docs/src/AI-MEP-Operators.html#smoothMutation) + and [crossover](https://hackage.haskell.org/package/hmep-0.1.0/docs/src/AI-MEP-Operators.html#crossover) + [operators](https://hackage.haskell.org/package/hmep-0.1.0/docs/AI-MEP.html). + One is also free to use their own operators. + * **Versatility**. `hmep` can be applied to solve regression problems with + one or multiple outputs. It means, you can approximate unknown functions + or solve classification tasks. The only requirement is a custom + [loss function](https://github.com/masterdezign/hmep/blob/b006eb8e0ca7c0540de979631423753bf0b66750/app/Main.hs#L67). + + ## How to build Use [Stack](http://haskellstack.org). @@ -43,8 +82,11 @@ $ 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): +### Example 1 +Now that the package is built, run the first demo to +express `cos^2(x)` through `sin(x)`: + $ stack exec hmep-demo Average loss in the initial population 15.268705681244962 @@ -58,6 +100,41 @@ v1 = sin x0 v2 = v1 * v1 result = 1 - v2 + +Effectively, the solution `cos^2(x) = 1 - sin^2(x)` was found. +Of course, MEP is a stochastic method, meaning that there is +no guarantee to find the globally optimal solution. + +The unknown function approximation problem can be illustrated +by the following suboptimal solution for a given set of random +data points (blue crosses). This example was produced by another run of +the [same demo](app/Main.hs), after 100 generations of 100 chromosomes. +The following expression was obtained +`y(x) = 3*0.31248786462471034 - sin(sin^2(x))`. +Interestingly, the approximating function lies symmetrically +in-between the extrema of the unknown function, approximately +described by the blue crosses. + +![Figure](https://github.com/masterdezign/hmep/blob/bbc2bdbac4fa3269c506455a473dddfa0e95231c/doc/Figures/cos2_approx.png) + +### Example 2 + +A similar example is to approximate `sin(x)` using only +addition and multiplication operators, i.e. with polynomials. + + $ stack exec hmep-sin-approximation + +The algorithm is able to automatically figure out the +powers of `x`. That is where MEP really shines. We [calculate](app/MainSin.hs) +30 expressions represented by each chromosome with practically no +additional computational penalty. Then, we +choose the best expression. In this run, we have automatically obtained a +[seventh degree polynomial](https://github.com/masterdezign/hmep/blob/master/doc/sin_approx.py#L12) +coded by 14 genes. Pretty cool, yeah? + +The result of approximation is [visualized](doc/sin_approx.py) below: + +![Figure](https://github.com/masterdezign/hmep/blob/d173e96acf72e482474e657880f8bd28c40694e7/doc/Figures/sin_approx.png) ## Authors diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -1,10 +1,27 @@ 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 +2. Improve code generation: + Perform subexpression elimination, + e.g. x0 / x0 should be reduced to 1 + 1 + 1 should become 2, etc. + And possibly, support unary operators: + Interpreted expression + + v1 = x0 / x0 + v3 = s x0 v1 + v6 = v3 / v1 + v14 = s x0 v6 + v15 = v14 * v6 + result = v1 - v15 + + should eventually become + + v1 = sin x0 + v2 = v1 * v1 + result = 1 - v2 + 3. Improve the demo: provide a CLI interface to work with external data @@ -12,3 +29,5 @@ 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. + +5. Implement subpopulations that eventually exchange individuals diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -22,60 +22,48 @@ c'ops = V.fromList [ ('*', (*)), ('+', (+)), - ('/', (/)), + -- Avoid division by zero + ('/', \x y -> if y < 1e-6 then 1 else x / y), ('-', (-)), ('s', \x _ -> sin x) ] -- Chromosome length , c'length = 50 + -- Probability to generate a new variable gene + , p'var = 0.1 + -- Probability to generate a new constant gene + , p'const = 0.05 + -- Probability to generate a new operator is + -- inferred as 1 - 0.1 - 0.5 = 0.85 } -- | 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 - --- 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 +dist x y = abs $ x - y main :: IO () main = do -- A vector of 50 random numbers between 0 and 1 (including 1) - xs <- runRandIO (vectorOf 50 double) + let datasetSize = 50 + xs <- runRandIO (vectorOf datasetSize 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' + dataset = map (\x -> (x, function x)) $ V.toList 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' + let loss = regressionLoss1 dist dataset -- Evaluate the initial population - let popEvaluated = evaluateGeneration loss pop + let popEvaluated = evaluatePopulation loss pop + norm = fromIntegral datasetSize - putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated) + putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated / norm) -- Declare how to produce the new generation let nextGeneration = evolve config loss (mutation3 config) crossover binaryTournament @@ -83,13 +71,12 @@ -- 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) + putStrLn $ "Population " ++ show (i * generations) ++ ": average loss " ++ show (avgLoss newPop / norm) return newPop where generations = 5 - -- Final generation + -- The final population final <- foldM runIO popEvaluated [1..20] - let best = last final - print best + putStrLn "Interpreted expression:" - putStrLn $ generateCode best + putStrLn $ generateCode (best final) diff --git a/app/MainSin.hs b/app/MainSin.hs new file mode 100644 --- /dev/null +++ b/app/MainSin.hs @@ -0,0 +1,83 @@ +module Main where + +{- + | = Sine approximation + + Generates an expression approximating sin(x) within the + interval of [-pi;pi]. + + Note: works the best when there is no division operation + involved. +-} + +import qualified Data.Vector as V +import Data.List ( foldl' ) +import Control.Monad ( foldM ) +import Math.Probable.Random -- From `probable` package + ( vectorOf + , double + ) + +import AI.MEP + +config = defaultConfig { + -- Functions available to genetically produced programs + c'ops = V.fromList [ + ('*', (*)), + ('+', (+)) + ] + -- Chromosome length + , c'length = 30 + , p'mutation = 0.05 + -- Probability to generate a new variable gene + , p'var = 0.1 + -- Probability to generate a new constant gene + , p'const = 0.05 + -- Probability to generate a new operator is + -- inferred as 1 - 0.1 - 0.05 = 0.85 + , c'popSize = 200 + } + +-- | Absolute value distance between two scalar values +dist :: Double -> Double -> Double +dist x y = abs $ x - y + +main :: IO () +main = do + -- A vector of 50 random numbers between 0 and 1 (including 1) + let datasetSize = 50 + xs <- runRandIO (vectorOf datasetSize 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 = sin + -- Pairs (x, f(x)) + dataset = map (\x -> (x, function x)) $ V.toList xs' + + -- Randomly create a population of chromosomes + pop <- runRandIO $ initialize config + + let loss = regressionLoss1 dist dataset + + -- Evaluate the initial population + let popEvaluated = evaluatePopulation loss pop + norm = fromIntegral datasetSize + + putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated / norm) + + -- 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 / norm) + return newPop + where generations = 5 + + -- The final population + final <- foldM runIO popEvaluated [1..300] + + putStrLn "Interpreted expression:" + putStrLn $ generateCode (best final) diff --git a/cabal.config b/cabal.config deleted file mode 100644 --- a/cabal.config +++ /dev/null @@ -1,48 +0,0 @@ -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/doc/sin_approx.py b/doc/sin_approx.py new file mode 100644 --- /dev/null +++ b/doc/sin_approx.py @@ -0,0 +1,55 @@ +from random import random +from pylab import * + +rcParams['font.family'] = 'serif' +rcParams['font.size'] = 14 + +def dist1(x,y): + return abs(x-y) + +# The polynomial approximation sin(x) +# TODO: algebraic expression elimination +def approx(x0): + v1 = -5.936286355387799e-2 + -5.936286355387799e-2 + v4 = x0 + x0 + v5 = v1 * x0 + v7 = v4 * x0 + v8 = v1 * v5 + v9 = x0 * x0 + v10 = v8 * v9 + v11 = x0 * v10 + v15 = -5.936286355387799e-2 * x0 + v18 = v10 * v11 + v20 = v7 * v15 + v21 = v15 + x0 + v25 = v21 + v20 + return v18 + v25 + +def main(): + xs = linspace(-pi,pi,300) + a = [approx(x) for x in xs] + + distances = [dist1(approx(x), sin(x)) for x in xs] + avg = sum(distances) / len(xs) + print("Average distance for %d points: %.4f" % (len(xs), avg)) + # Average distance for 300 points: 0.0303 + + # Visualization + + # Note, the random domain during each demo run is different. + # The values below are given for illustration purposes only. + randXs = sort([(random()-0.5)*2*pi for i in range(50)]) + correct = [sin(x) for x in randXs] + + plot(randXs, correct, '+') + plot(xs, a, lw=1.2) + xlabel('x') + ylabel('y(x)') + xlim([-pi, pi]) + xticks([-pi, 0, pi], ['-$\pi$', '0', '$\pi$']) + legend(['Points from sin(x)', 'Polynomial Approx.'], fontsize=10) + show() + + +if __name__ == '__main__': + main() diff --git a/hmep.cabal b/hmep.cabal --- a/hmep.cabal +++ b/hmep.cabal @@ -1,5 +1,5 @@ name: hmep -version: 0.1.0 +version: 0.1.1 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 cabal.config +extra-source-files: README.md TODO CHANGELOG.md doc/sin_approx.py cabal-version: >=1.22 library @@ -27,7 +27,6 @@ , mwc-random , primitive , probable - , statistics >= 0.10 && < 0.14 , vector default-language: Haskell2010 @@ -37,10 +36,20 @@ ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends: base , probable - , statistics >= 0.10 && < 0.14 , vector , hmep default-language: Haskell2010 + +executable hmep-sin-approximation + hs-source-dirs: app + main-is: MainSin.hs + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: base + , probable + , vector + , hmep + default-language: Haskell2010 + test-suite hmep-test type: exitcode-stdio-1.0 diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -18,9 +18,9 @@ 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^8" ~: V.fromList [2, 4, 16, 256] ~=? evaluateChromosome 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)) + "2.5^8" ~: V.fromList [2.5,6.25,39.0625,1525.87890625] ~=? evaluateChromosome pow8 (V.singleton (2.5 :: Double)) ] allTests = TestList [ testEvaluate ]