Drezil
/
chemodiversity
1
1
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Simulation looks ok-ish. Needs incentive to foster productive enzymes

This commit is contained in:
Nicole Dresselhaus 2018-05-15 19:01:38 +02:00
parent cd16dbb39b
commit 8eeb837b9f
Signed by: Drezil
GPG Key ID: 057D94F356F41E25
6 changed files with 92 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
app/Main.hs
View File

@ -1,3 +1,4 @@
{-# LANGUAGE BangPatterns #-}
module Main where
import Text.Printf
@ -6,6 +7,7 @@ import Numeric.LinearAlgebra
import Data.List
import System.Random
import Control.Concurrent
import Control.Parallel.Strategies
import qualified Debug.Trace as Debug
import System.IO
@ -30,18 +32,18 @@ greenfly = Predator [] 0.2 -- killed by any toxic Component
-- Environment
exampleEnvironment :: Int -> [Enzyme] -> Environment
exampleEnvironment addedC es =
exampleEnvironment :: Int -> [Enzyme] -> [(Predator,Probability)] -> [(Compound,Amount)] -> Environment
exampleEnvironment addedC es pred tox =
Environment
{ soil = [ (Nitrate, 2)
, (Phosphor, 3)
, (Photosynthesis, 10)
]
, predators = [ (greenfly, 0.1) ]
, predators = pred -- [ (greenfly, 0.1) ]
, metabolismIteration = 100
, maxCompound = maxCompoundWithoutGeneric + addedC
, toxicCompounds = [(Produced FPP,0.5)] --FPP kills 100% if produced amount above 0.2 units
, possibleEnzymes = [pps,fpps] ++ es
, toxicCompounds = tox --[(Produced FPP,0.1)] ++ tox
, possibleEnzymes = es -- [pps,fpps] ++ es
}
-- Plants
@ -84,11 +86,12 @@ loop loopAmount = loop' loopAmount 0
putStrLn ""
putStrLn $ "Generation " ++ show curLoop ++ " of " ++ show loopAmount ++ ":"
newPlants <- flip runReaderT e $ do
fs <- sequence $ fitness <$> plants
fs <- sequence (fitness <$> plants)
let fps = zip plants fs -- gives us plants & their fitness in a tuple
sumFitness = sum fs
pe <- asks possibleEnzymes
liftIO $ printPopulation pe fps
tc <- fmap fst <$> asks toxicCompounds
liftIO $ printPopulation tc pe fps
-- generate 100 new plants.
sequence . flip fmap [1..100] $ \_ -> do
parent' <- liftIO $ randomRIO (0,sumFitness)
@ -110,12 +113,18 @@ main :: IO ()
main = do
hSetBuffering stdin NoBuffering
hSetBuffering stdout NoBuffering
randomCompounds <- generateTreeFromList 10 (toEnum <$> [(maxCompoundWithoutGeneric+1)..] :: [Compound]) -- generate roughly 10 compounds
let env = exampleEnvironment (getTreeSize randomCompounds) (generateEnzymeFromTree randomCompounds)
emptyPlants = replicate 100 emptyPlant
randomCompounds <- makeHead (Substrate Photosynthesis) <$> generateTreeFromList 40 (toEnum <$> [(maxCompoundWithoutGeneric+1)..] :: [Compound]) -- generate roughly x compounds
ds <- randoms <$> newStdGen
probs <- randomRs (0.2,0.7) <$> newStdGen
let emptyPlants = replicate 100 emptyPlant
poisonedTree = poisonTree ds randomCompounds
poisonCompounds = foldMap (\(a,b) -> if a > 0.5 then [(b,a)] else []) $ poisonedTree
predators <- generatePredators 0.5 poisonedTree
let env = exampleEnvironment (getTreeSize randomCompounds) (generateEnzymeFromTree randomCompounds) (zip predators probs) poisonCompounds
printEnvironment env
writeFile "poison.twopi" $ generateDotFromPoisonTree "poison" 0.5 $ poisonedTree
putStr "\ESC[?1049h"
loop 100 emptyPlants env
loop 200 emptyPlants env
putStrLn "Simulation ended. Press key to exit."
_ <- getChar
putStr "\ESC[?1049l"
@ -128,6 +137,21 @@ main = do
-- printf "%15.2f" f
-- putStr "\n"
generatePredators :: Double -> EnzymeTree s (Double,Compound) -> IO [Predator]
generatePredators threshold t = do
ps <- mapM generatePredators' $ getSubTrees t
return $ concat ps
where
generatePredators' :: (EnzymeTree s (Double, Compound)) -> IO [Predator]
generatePredators' t = do -- not fully resistant to t, but fully resistant to everything in ts
let comps = foldMap (\(a,b) -> if a > threshold then [(a,b)] else []) t
amount <- randomRIO (0,length comps + 1) :: IO Int
forM [1..amount] $ \_ -> do
impact <- randomRIO (0.2,0.7)
rands <- randoms <$> newStdGen
let unresists = foldMap (\((a,b),r) -> if r*2 < a then [b] else []) $ zip comps rands
return $ Predator unresists impact
printEnvironment :: Environment -> IO ()
printEnvironment (Environment soil pred metaIter maxComp toxic possEnz) =
do
@ -138,16 +162,17 @@ printEnvironment (Environment soil pred metaIter maxComp toxic possEnz) =
putStrLn $ "Compounds: " ++ show ((toEnum <$> [0..maxComp]) :: [Compound])
putStrLn $ "Toxic: " ++ show toxic
printPopulation :: [Enzyme] -> [(Plant,Double)] -> IO ()
printPopulation es ps = do
printPopulation :: [Compound] -> [Enzyme] -> [(Plant,Double)] -> IO ()
printPopulation toxins es ps = do
let padded i str = take i $ str ++ repeat ' '
putStr $ padded 40 "Population:"
putStr $ padded 50 "Population:"
forM_ ps $ \(_,f) -> putStr (printColor f '█')
putStrLn colorOff
forM_ es $ \e -> do
putStr $ padded 40 (show (enzymeName e))
putStr $ if (fst . snd . synthesis $ e) `elem` toxins then "\ESC[31m" ++ padded 50 (show (enzymeName e)) ++ "\ESC[0m"
else padded 50 (show (enzymeName e))
forM_ ps $ \(Plant g _,_) -> do
let curE = sum $ map (\(_,q,a) -> fromIntegral q*a)
let curE = sum $ map (\(_,q,a) -> fromIntegral q*a)
. filter (\(e',_,_) -> e == e')
$ g
plot x

1
package.yaml
View File

@ -26,6 +26,7 @@ dependencies:
- random
- QuickCheck
- pretty-simple
- parallel
library:
source-dirs: src

35
src/ArbitraryEnzymeTree.hs
View File

@ -3,6 +3,7 @@
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE GADTs #-}
module ArbitraryEnzymeTree
( EnzymeTree
, getTreeSize
@ -12,6 +13,9 @@ module ArbitraryEnzymeTree
, treeFromList
, generateTreeFromList
, generateDotFromTree
, generateDotFromPoisonTree
, poisonTree
, makeHead
) where
import Test.QuickCheck
@ -89,3 +93,34 @@ generateDotFromTree name t = "digraph " <> name <> " {\n"
where
ts = fromString . show
wrap x = "\"" <> x <> "\""
generateDotFromPoisonTree :: (Show a, IsString b, Monoid b) => b -> Double -> EnzymeTree s (Double,a) -> b
generateDotFromPoisonTree name pl t = "digraph " <> name <> " {\n"
<> generateDotFromTree' t
<> "}\n"
where
generateDotFromTree' :: (Show a, IsString b, Monoid b) => EnzymeTree s (Double,a) -> b
generateDotFromTree' (EnzymeTree _ (d,c) ns) =
" " <> wrap (ts c) <> " -> { "
<> mconcat (intersperse " " (wrap . ts . snd . getElement <$> ns))
<> " };\n"
<> (if d > pl then " " <> wrap (ts c) <> " [style=filled, fillcolor=\"0," <> ts d <> ",0.9\"];\n" else "")
<> mconcat (generateDotFromTree' <$> ns)
where
ts :: (Show a, IsString b) => a -> b
ts = fromString . show
wrap x = "\"" <> x <> "\""
poisonTree :: [Double] -> EnzymeTree s a -> EnzymeTree s (Double, a)
poisonTree ds t@(EnzymeTree s _ _) = go Nothing 0 annotatedTree
where
annotatedTree = (,) <$> treeFromList s ds <*> t
go :: Maybe Double -> Int -> EnzymeTree t (Double, a) -> EnzymeTree t (Double, a)
go Nothing i parent@(EnzymeTree s' (p,a) childs) = EnzymeTree s' (p/5, a) $ (\(_,ts) -> go (Just $ p/5 ) (i+1) ts) <$> zip [1..] childs
go (Just pe) i this@(EnzymeTree s' (p,a) childs) = EnzymeTree s' (p',a) $ (\(j,ts) -> go (Just . min 1 $ j*p') (i+1) ts) <$> zip [1..] childs
where
p' = max pe (p / i')
i' = fromIntegral $ 6 - min i 5 -- 100% effective poision only at level xx or deeper
makeHead :: a -> EnzymeTree s a -> EnzymeTree s a
makeHead c (EnzymeTree s a ts) = EnzymeTree s c ts

25
src/Environment.hs
View File

@ -4,6 +4,7 @@ import Data.Functor ((<$>))
import Control.Applicative ((<*>))
import Control.Monad (forM_)
import Control.Monad.Reader
import Control.Parallel.Strategies
import Data.List (permutations, subsequences)
import Numeric.LinearAlgebra
import Text.Printf
@ -78,8 +79,8 @@ makeSimpleEnzyme a b = Enzyme (show a ++ " -> " ++ show b) [] ((a,-1),(b,1)) Not
-- | In the environment we have predators that impact the fitness of our plants and
-- may be resistant to some compounds the plant produces. They can also differ in
-- their intensity.
data Predator = Predator { resistance :: [Compound]
-- ^ list of components this predator is resistant to
data Predator = Predator { irresistance :: [Compound]
-- ^ list of components this predator is not resistant to
, fitnessImpact :: Amount
-- ^ impact on the fitness of a plant
-- (~ agressiveness of the herbivore)
@ -150,7 +151,7 @@ fitness p = do
products <- produceCompounds p nutrients -- produce compounds
survivalRate <- deterPredators products -- defeat predators with produced compounds
let sumEnzymes = sum $ (\(_,q,a) -> fromIntegral q*a) <$> genome p -- amount of enzymes * activation = resources "wasted"
costOfEnzymes = 0.95 ** sumEnzymes
costOfEnzymes = 0.99 ** sumEnzymes
return $ survivalRate * costOfEnzymes
-- can also be written as, but above is more clear.
-- fitness p = absorbNutrients p >>= produceCompounds p >>= deterPredators
@ -179,10 +180,10 @@ deterPredators cs = do
ts <- asks toxicCompounds
let
deter :: Predator -> Double
-- 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, t `notElem` resistance p]
-- multiply (toxicity of t with 100% effectiveness at l| for all toxins t | and t not in p's irresistance-list)
deter p = product [1 - min 1 (cs ! fromEnum t / l) | (t,l) <- ts, t `elem` irresistance p]
-- multiply (probability of occurence * intensity of destruction / probability to deter predator | for all predators)
return . product $ [min 1 ((1-prob) * fitnessImpact p / deter p) | (p,prob) <- ps]
return $ product ([min 1 ((1-prob) * fitnessImpact p / deter p) | (p,prob) <- ps] `using` parList rdeepseq)
-- Mating & Creation of diversity
-- ------------------------------
@ -208,26 +209,26 @@ haploMate (Plant genes abs) = do
. 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 (r:rs) ((e,1,a):gs) = if r < 0.1 then deleteGene rs gs else (e,1,a):deleteGene rs gs
deleteGene (r:rs) ((e,q,a):gs) = if r < 0.1 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 (r:rs) ((e,q,a):gs) = if r < 0.1 then (e,1,a):(e,q,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
addGene (r:rs) (s:ss) g = if r < 0.05 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
mutateGene (r:rs) (s:ss) ((e,1,a):gs) = if r < 0.01 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
mutateGene (r:rs) (s:ss) ((e,q,a):gs) = if r < 0.01 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

6
src/Lib.hs
View File

@ -1,6 +0,0 @@
module Lib
( someFunc
) where
someFunc :: IO ()
someFunc = putStrLn "someFunc"

4
stack.yaml
View File

@ -15,7 +15,7 @@
# resolver:
# name: custom-snapshot
# location: "./custom-snapshot.yaml"
resolver: lts-11.7
resolver: lts-11.9
# User packages to be built.
# Various formats can be used as shown in the example below.
@ -63,4 +63,4 @@ packages:
# extra-lib-dirs: [/path/to/dir]
#
# Allow a newer minor version of GHC than the snapshot specifies
# compiler-check: newer-minor
# compiler-check: newer-minor