From 85ce37b1065aeb84b5376564cdf97006c26bede1 Mon Sep 17 00:00:00 2001 From: Stefan Dresselhaus Date: Wed, 2 May 2018 23:39:22 +0200 Subject: [PATCH] working evolution-simulation (POC) --- app/Main.hs | 170 ++++++++++++++++++++++++++++++++++++++++++++- package.yaml | 1 + src/Environment.hs | 151 +++++++++++++++------------------------- 3 files changed, 224 insertions(+), 98 deletions(-) diff --git a/app/Main.hs b/app/Main.hs index de1c1ab..f9f13ce 100644 --- a/app/Main.hs +++ b/app/Main.hs @@ -1,6 +1,172 @@ +{-# LANGUAGE TypeApplications #-} + module Main where -import Lib +import Environment +import Text.Printf +import Control.Monad.Reader +import Numeric.LinearAlgebra +import Data.List +import System.Random +import Control.Concurrent +import qualified Debug.Trace as Debug +import System.IO + + +-- Example definitions +-- ------------------- + +-- Enzymes + +pps :: Enzyme -- uses Phosphor from Substrate to produce PP +pps = Enzyme "PPS" [(Substrate Phosphor,1)] ((Substrate Phosphor,(-1)),(Produced PP,1)) Nothing + +fpps :: Enzyme -- PP -> FPP +fpps = makeSimpleEnzyme (Produced PP) (Produced FPP) + +-- Predator + +greenfly :: Predator -- 20% of plants die to greenfly, but the fly is +greenfly = Predator [] 0.2 -- killed by any toxic Component + +-- Environment + +exampleEnvironment :: Environment +exampleEnvironment = + Environment + { soil = [ (Nitrate, 2) + , (Phosphor, 3) + , (Photosynthesis, 10) + ] + , predators = [ (greenfly, 0.1) ] + , metabolismIteration = 100 + , maxCompound = maxCompoundWithoutGeneric + 100 + , toxicCompounds = [(Produced FPP,0.5)] --FPP kills 100% if produced amount above 0.2 units + , possibleEnzymes = [pps,fpps] + } + +-- Plants + +examplePlants :: [Plant] +examplePlants = (\g -> Plant g defaultAbsorption) <$> genomes + where + enzymes = [pps, fpps] + quantity = [1,2] :: [Quantity] + activation = [0.7, 0.9, 1] + + genomes = do + e <- permutations enzymes + e' <- subsequences e + q <- quantity + a <- activation + return $ (,,) <$> e' <*> [q] <*> [a] + + defaultAbsorption = soil <$> ask >>= return . fmap ( limit Phosphor 2 + . limit Nitrate 1 + . limit Sulfur 0 + ) + -- custom absorbtion with helper-function: + limit :: Nutrient -> Amount -> (Nutrient, Amount) -> (Nutrient, Amount) + limit n a (n', a') + | n == n' = (n, min a a') -- if we should limit, then we do ;) + | otherwise = (n', a') + +-- Running the simulation +-- ---------------------- + +loop :: Int -> [Plant] -> Environment -> IO () +loop loopAmount plants e = loop' loopAmount 0 plants e + + where + loop' :: Int -> Int -> [Plant] -> Environment -> IO () + loop' loopAmount curLoop plants e = unless (loopAmount == curLoop) $ do + putStr $ "\ESC[2J\ESC[H" + printEnvironment e + putStrLn "" + putStrLn $ "Generation " ++ show curLoop ++ " of " ++ show loopAmount ++ ":" + newPlants <- (flip runReaderT) e $ do + fs <- sequence $ fitness <$> plants + let fps = zip plants fs -- gives us plants & their fitness in a tuple + sumFitness = sum fs + pe <- possibleEnzymes <$> ask + liftIO $ printPopulation pe fps + -- generate 100 new plants. + sequence . (flip fmap) [1..100] $ \_ -> do + parent' <- liftIO $ randomRIO (0,sumFitness) + let + -- if we only have one parent in our list, take it. + findParent :: Double -> [(Plant,Double)] -> Plant + findParent _ [(last,_)] = last + -- otherwise count down x to find the parent in the list + findParent x ((p,f):ps) + | x < f = p + | otherwise = findParent (x-f) ps + parent = findParent parent' fps + haploMate parent + hFlush stdout + threadDelay $ 100*1000 -- sleep 100ms + loop' loopAmount (curLoop+1) newPlants e main :: IO () -main = someFunc +main = do + hSetBuffering stdin NoBuffering + hSetBuffering stdout NoBuffering + let emptyPlants = replicate 100 emptyPlant + printEnvironment exampleEnvironment + putStr "\ESC[?1049h" + loop 100 emptyPlants exampleEnvironment + putStrLn "Simulation ended. Press key to exit." + _ <- getChar + putStr "\ESC[?1049l" + + -- fitness <- runReaderT (sequence $ (\a -> do p <- absorbNutrients a >>= produceCompounds a; (,,) a p <$> deterPredators p) <$> emptyPlants) exampleEnvironment + -- mapM_ (printf "%15.15s, " . show . toEnum @Compound) [0..maxCompoundWithoutGeneric] + -- putStrLn "Fitness" + -- forM_ fitness $ \(p, c, f) -> do + -- mapM_ (printf "%15.2f, ") (toList c) + -- printf "%15.2f" f + -- putStr "\n" + +printEnvironment :: Environment -> IO () +printEnvironment (Environment soil pred metaIter maxComp toxic possEnz) = + do + putStrLn "Environment:" + putStrLn $ "Soil: " ++ show soil + putStrLn $ "Predators: " ++ show pred + putStrLn $ "PSM Iters: " ++ show metaIter + putStrLn $ "Compounds: " ++ show ((toEnum <$> [0..maxComp]) :: [Compound]) + putStrLn $ "Toxic: " ++ show toxic + +printPopulation :: [Enzyme] -> [(Plant,Double)] -> IO () +printPopulation es ps = do + let padded i str = take i $ str ++ repeat ' ' + putStr $ padded 40 "Population:" + forM_ ps $ \((_,f)) -> putStr (printColor f '█') + putStrLn colorOff + forM_ es $ \e -> do + putStr $ padded 40 (show (enzymeName e)) + forM_ ps $ \((Plant g _,_)) -> do + let curE = sum $ map (\(_,q,a) -> (fromIntegral q)*a) + . filter (\(e',_,_) -> e == e') + $ g + plot x + | x > 2 = "O" + | x > 1 = "+" + | x > 0.7 = "ö" + | x > 0.5 = "o" + | x > 0 = "." + | otherwise = "_" + putStr (plot curE) + putStrLn "" + +printColor :: Double -> Char -> String +printColor x c + | x*x < 0.5 = "\ESC[38;5;" ++ (show $ 16 + 36*5 + 6*(floor $ 5*2*x') + 0) ++ "m" ++ [c] ++ "" + | otherwise = "\ESC[38;5;" ++ (show $ 16 + 36*(floor $ 5*2*(1-x')) + 6*5 + 0) ++ "m" ++ [c] ++ "" + -- 32 bit + -- | x*x < 0.5 = "\ESC[38;2;255;" ++ (show . floor $ 255*2*x') ++ ";0m" ++ [c] ++ "" + -- | otherwise = "\ESC[38;2;" ++ (show . floor $ 255*2*(1-x')) ++ ";255;0m" ++ [c] ++ "" + where x' = x*x + +colorOff :: String +colorOff = "\ESC[0m" diff --git a/package.yaml b/package.yaml index ddf599e..8f1171a 100644 --- a/package.yaml +++ b/package.yaml @@ -23,6 +23,7 @@ dependencies: - base >= 4.7 && < 5 - hmatrix - mtl +- random library: source-dirs: src diff --git a/src/Environment.hs b/src/Environment.hs index ddf4137..2b1525b 100644 --- a/src/Environment.hs +++ b/src/Environment.hs @@ -11,7 +11,7 @@ import Control.Monad.Reader import Data.List (permutations, subsequences) import Numeric.LinearAlgebra import Text.Printf -import qualified Debug.Trace as Debug +import System.Random type Probability = Double type Quantity = Int @@ -76,14 +76,6 @@ data Enzyme = Enzyme makeSimpleEnzyme :: Compound -> Compound -> Enzyme makeSimpleEnzyme a b = Enzyme (show a ++ " -> " ++ show b) [] ((a,-1),(b,1)) Nothing ---Example "enzymes" could be: - -pps :: Enzyme -- uses Phosphor from Substrate to produce PP -pps = Enzyme "PPS" [(Substrate Phosphor,1)] ((Substrate Phosphor,(-1)),(Produced PP,1)) Nothing - -fpps :: Enzyme -- PP -> FPP -fpps = makeSimpleEnzyme (Produced PP) (Produced FPP) - -- Evironment -- ---------- @@ -97,11 +89,6 @@ data Predator = Predator { resistance :: [Compound] -- (~ agressiveness of the herbivore) } deriving (Show, Eq) --- Exemplatory: - -greenfly :: Predator -- 20% of plants die to greenfly, but the fly is -greenfly = Predator [] 0.2 -- killed by any toxic Component - -- The environment itself is just the soil and the predators. Extensions would be possible. data Environment = @@ -122,27 +109,14 @@ data Environment = , toxicCompounds :: [(Compound,Amount)] -- ^ Compounds considered to be toxic in this environment. -- Kills 100% of Predators above Amount. + , possibleEnzymes :: [Enzyme] + -- ^ All enzymes that can be created by genetic manipulation in this setting. } deriving (Show, Eq) -- helper function. Allows for [0..maxCompoundWithoutGeneric] :: [Compound] with all non-generic Compounds maxCompoundWithoutGeneric :: Int maxCompoundWithoutGeneric = fromEnum (maxBound :: Nutrient) + fromEnum (maxBound :: Component) + 1 --- Example: - -exampleEnvironment :: Environment -exampleEnvironment = - Environment - { soil = [ (Nitrate, 2) - , (Phosphor, 3) - , (Photosynthesis, 10) - ] - , predators = [ (greenfly, 0.1) ] - , metabolismIteration = 100 - , maxCompound = maxCompoundWithoutGeneric - , toxicCompounds = [(Produced FPP,0.5)] --FPP kills 100% if produced amount above 0.2 units - } - type World a = ReaderT Environment IO a -- Plants @@ -164,35 +138,6 @@ instance Show Plant where instance Eq Plant where a == b = genome a == genome b --- | The following example yields in the example-environment this population: --- --- >>> printPopulation [pps, fpps] plants --- Population: --- PPS ______oöö+++______oöö+++____________oöö+++oöö+++ --- FPPS ____________oöö+++oöö+++______oöö+++______oöö+++ -plants :: [Plant] -plants = (\g -> Plant g defaultAbsorption) <$> genomes - where - enzymes = [pps, fpps] - quantity = [1,2] :: [Quantity] - activation = [0.7, 0.9, 1] - - genomes = do - e <- permutations enzymes - e' <- subsequences e - q <- quantity - a <- activation - return $ (,,) <$> e' <*> [q] <*> [a] - - defaultAbsorption = soil <$> ask >>= return . fmap ( limit Phosphor 2 - . limit Nitrate 1 - . limit Sulfur 0 - ) - -- custom absorbtion with helper-function: - limit :: Nutrient -> Amount -> (Nutrient, Amount) -> (Nutrient, Amount) - limit n a (n', a') - | n == n' = (n, min a a') -- if we should limit, then we do ;) - | otherwise = (n', a') -- Fitness -- ------- @@ -208,7 +153,9 @@ fitness p = do nutrients <- absorbNutrients p -- absorb soil products <- produceCompounds p nutrients -- produce compounds survivalRate <- deterPredators products -- defeat predators with produced compounds - return survivalRate + let sumEnzymes = sum $ (\(_,q,a) -> (fromIntegral q)*a) <$> genome p -- amount of enzymes * activation = resources "wasted" + costOfEnzymes = 0.95 ** sumEnzymes + return $ survivalRate * costOfEnzymes -- can also be written as, but above is more clear. -- fitness p = absorbNutrients p >>= produceCompounds p >>= deterPredators @@ -239,52 +186,64 @@ deterPredators cs = do -- multiply (toxicity of t with 100% effectiveness at l| for all toxins t | and t not in p's resistance-list) deter p = product [1 - min 1 (cs ! (fromEnum t) / l) | (t,l) <- ts, not (t `elem` resistance p)] -- multiply (probability of occurence * intensity of destruction / probability to deter predator | for all predators) - return . product $ [min 1 (prob * fitnessImpact p / deter p) | (p,prob) <- ps] + return . product $ [min 1 ((1-prob) * fitnessImpact p / deter p) | (p,prob) <- ps] -- Mating & Creation of diversity -- ------------------------------ --- Running the simulation --- ---------------------- + +-- | mate haploid +haploMate :: Plant -> World Plant +haploMate (Plant genes abs) = do + --generate some random infinite uniform distributed lists of doubles in [0,1) + r1 <- liftIO ((randoms <$> newStdGen) :: IO [Double]) + r2 <- liftIO ((randoms <$> newStdGen) :: IO [Double]) + r3 <- liftIO ((randoms <$> newStdGen) :: IO [Double]) + r4 <- liftIO ((randoms <$> newStdGen) :: IO [Double]) + r5 <- liftIO ((randoms <$> newStdGen) :: IO [Double]) + enzymes <- possibleEnzymes <$> ask + re1 <- liftIO ((randomRs (0,length enzymes - 1) <$> newStdGen) :: IO [Int]) + re2 <- liftIO ((randomRs (0,length enzymes - 1) <$> newStdGen) :: IO [Int]) + let + genes' = mutateGene r1 re1 + . noiseActivation r2 + . addGene r3 re2 + . duplicateGene r4 + . deleteGene r5 + $ genes + deleteGene :: [Double] -> Genome -> Genome + deleteGene (r:rs) ((e,1,a):gs) = if a < 0.1 && r < 0.5 then deleteGene rs gs else (e,1,a):deleteGene rs gs + deleteGene (r:rs) ((e,q,a):gs) = if a < 0.1 && r < 0.5 then (e,q-1,a):deleteGene rs gs else (e,q,a):deleteGene rs gs + deleteGene _ [] = [] + + duplicateGene :: [Double] -> Genome -> Genome + duplicateGene (r:rs) ((e,q,a):gs) = if r < 0.05 then (e,q+1,a):duplicateGene rs gs else (e,q,a):duplicateGene rs gs + duplicateGene _ [] = [] + + addGene :: [Double] -> [Int] -> Genome -> Genome + addGene (r:rs) (s:ss) g = if r < 0.01 then ((enzymes !! s),1,1):g else g + + noiseActivation :: [Double] -> Genome -> Genome + noiseActivation (r:rs) ((e,q,a):gs) = (e,q,max 0 $ min 1 $ a-0.01+0.02*r):noiseActivation rs gs + noiseActivation _ [] = [] + + mutateGene :: [Double] -> [Int] -> Genome -> Genome + mutateGene (r:rs) (s:ss) ((e,1,a):gs) = if r < 0.05 then ((enzymes !! s),1,a):mutateGene rs ss gs + else (e,1,a):mutateGene rs ss gs + + mutateGene (r:rs) (s:ss) ((e,q,a):gs) = if r < 0.05 then (e,q-1,a):((enzymes !! s),1,a):mutateGene rs ss gs + else (e,q,a):mutateGene rs ss gs + mutateGene (r:rs) (s:ss) [] = [] + return $ Plant genes' abs - -main = do - putStrLn "Environment:" - print exampleEnvironment - putStrLn "Example population:" - printPopulation [pps, fpps] plants - fitness <- runReaderT (sequence $ (\a -> do p <- absorbNutrients a >>= produceCompounds a; (,) p <$> deterPredators p) <$> plants) exampleEnvironment - mapM_ (printf "%15.15s, " . show . toEnum @Compound) [0..maxCompoundWithoutGeneric] - putStrLn "Fitness" - forM_ fitness $ \(p, f) -> do - mapM_ (printf "%15.2f, ") (toList p) - printf "%15.2f" f - putStr "\n" - -- Utility Functions -- ----------------- +-- | Plant with no secondary metabolism with unlimited extraction from environment. +emptyPlant :: Plant +emptyPlant = Plant [] (soil <$> ask) + getAmountOf :: Compound -> [(Compound, Amount)] -> Amount getAmountOf c = sum . fmap snd . filter ((== c) . fst) - -printPopulation :: [Enzyme] -> [Plant] -> IO () -printPopulation es ps = do - let padded i str = take i $ str ++ repeat ' ' - putStrLn "Population:" - forM_ es $ \e -> do - putStr $ padded 40 (show (enzymeName e)) - forM_ ps $ \(Plant g _) -> do - let curE = sum $ map (\(_,q,a) -> (fromIntegral q)*a) - . filter (\(e',_,_) -> e == e') - $ g - plot x - | x > 2 = "O" - | x > 1 = "+" - | x > 0.7 = "ö" - | x > 0.5 = "o" - | x > 0 = "." - | otherwise = "_" - putStr (plot curE) - putStrLn ""