From f4d8b00b7c6a8a621c0a7e11470d31df93405cbf Mon Sep 17 00:00:00 2001 From: Stefan Dresselhaus Date: Wed, 2 May 2018 17:22:29 +0200 Subject: [PATCH] PoC for single iteration. --- package.yaml | 1 + src/Environment.hs | 93 ++++++++++++++++++++++++++++++---------------- 2 files changed, 63 insertions(+), 31 deletions(-) diff --git a/package.yaml b/package.yaml index 3ca7cbd..ddf599e 100644 --- a/package.yaml +++ b/package.yaml @@ -22,6 +22,7 @@ description: Please see the README on GitHub at = 4.7 && < 5 - hmatrix +- mtl library: source-dirs: src diff --git a/src/Environment.hs b/src/Environment.hs index 9985dc6..078d7d1 100644 --- a/src/Environment.hs +++ b/src/Environment.hs @@ -6,6 +6,7 @@ module Environment where import Data.Functor ((<$>)) import Control.Applicative ((<*>)) import Control.Monad (forM_) +import Control.Monad.Reader import Data.List (permutations, subsequences) import Numeric.LinearAlgebra @@ -96,7 +97,7 @@ data Predator = Predator { resistance :: [Compound] -- Exemplatory: greenfly :: Predator -- 20% of plants die to greenfly, but the fly is -greenfly = Predator [Produced PP] 0.2 -- killed by any Component not being PP +greenfly = Predator [] 0.2 -- killed by any toxic Component -- The environment itself is just the soil and the predators. Extensions would be possible. @@ -106,8 +107,23 @@ data Environment = -- ^ soil is a list of nutrients available to the plant. , predators :: [(Predator, Probability)] -- ^ Predators with the probability of appearance in this generation. + , metabolismIteration :: Int + -- ^ Number of iterations for producing compounds + , maxCompound :: Int + -- ^ Number of possible Compounds + -- 'maxCompound' should be greater than #Nutrient + #Products. + -- Rest will get filled up with 'GenericEnzyme i' + -- + -- To find the 'maxCompound' without 'GenericEnzyme' use + -- 'maxComponent = fromEnum (maxBound :: Nutrient) + fromEnum (maxBound :: Component) + 1' + , toxicCompounds :: [Compound] + -- ^ Compounds considered to be toxic in this environment. } 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 @@ -118,8 +134,13 @@ exampleEnvironment = , (Photosynthesis, 10) ] , predators = [ (greenfly, 0.1) ] + , metabolismIteration = 100 + , maxCompound = maxCompoundWithoutGeneric + , toxicCompounds = [Produced FPP] } +type World a = ReaderT Environment IO a + -- Plants -- ------ @@ -131,7 +152,7 @@ type Genome = [(Enzyme, Quantity, Activation)] data Plant = Plant { genome :: Genome -- ^ the genetic characteristic of the plant - , absorbNutrients :: Environment -> [(Nutrient,Amount)] + , absorbNutrients :: World [(Nutrient,Amount)] -- ^ the capability to absorb nutrients given an environment } instance Show Plant where @@ -159,10 +180,10 @@ plants = (\g -> Plant g defaultAbsorption) <$> genomes a <- activation return $ (,,) <$> e' <*> [q] <*> [a] - defaultAbsorption (Environment s _) = limit Phosphor 2 - . limit Nitrate 1 - . limit Sulfur 0 - <$> s + 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') @@ -178,35 +199,40 @@ plants = (\g -> Plant g defaultAbsorption) <$> genomes type Fitness = Double -fitness :: Environment -> Plant -> Fitness -fitness e p = survivalRate - where - nutrients = absorbNutrients p e - products = produceCompounds p nutrients - survivalRate = deterPredators (predators e) products +fitness :: Plant -> World Fitness +fitness p = do + nutrients <- absorbNutrients p -- absorb soil + products <- produceCompounds p nutrients -- produce compounds + survivalRate <- deterPredators products -- defeat predators with produced compounds + return survivalRate + -- can also be written as, but above is more clear. + -- fitness p = absorbNutrients p >>= produceCompounds p >>= deterPredators - -produceCompounds :: Plant -> [(Nutrient, Amount)] -> Vector Amount -produceCompounds (Plant genes _) substrate = final - where - initialAmount = assoc 10 0 ((\(n,a) -> (fromEnum $ Substrate n,a)) <$> substrate) +produceCompounds :: Plant -> [(Nutrient, Amount)] -> World (Vector Amount) +produceCompounds (Plant genes _) substrate = do + numIter <- metabolismIteration <$> ask + numCompounds <- maxCompound <$> ask + let + initialAmount = assoc (numCompounds+1) 0 ((\(n,a) -> (fromEnum $ Substrate n,a)) <$> substrate) enzymes = (\(e,q,a) -> (synthesis e,(fromIntegral q)*a)) <$> genes -- [(((Component,Amount),(Component,Amount)),q*a)], Amount got * by quantity & activation - positions = concat $ (\(((i,ia),(o,oa)),f) -> [((fromEnum i,fromEnum i),ia),((fromEnum o,fromEnum i),oa)]) <$> enzymes -- [((row,column),amount)] - mat = accum (konst 0 (10::Int,10::Int)) (+) positions --accumulate all entries into one matrix. - (l,v) = eig ((*0.01) `cmap` mat) - final = (realPart `cmap` (v <> ((^100) `cmap` diag l) <> (tr v))) #> initialAmount + positions = concat $ (\(((i,ia),(o,oa)),f) -> [((fromEnum i,fromEnum i),2*f*ia),((fromEnum o,fromEnum i),oa),((fromEnum i,fromEnum o),f*oa)]) <$> enzymes -- [((row,column),amount)] + mat = accum (konst 0 (numCompounds+1,numCompounds+1)) (+) positions --accumulate all entries into one matrix. + -- mat is now the rate of change in u'(t) = A u(t) + (l,v) = eig (ident (numCompounds+1) + ((*0.01) `cmap` mat)) -- use u(t+1) = u(t) + A u(t) = (E + A) u(t) for iteration + final = (realPart `cmap` (v <> ((^numIter) `cmap` diag l) <> (tr v))) #> initialAmount -- (E + A)^numIter * t_0 for numIter iterations. + return final -deterPredators :: [(Predator, Probability)] -> Vector Amount -> Probability -deterPredators ps cs' = sum $ do - (c,a) <- cs -- for every compound - (p,prob) <- ps -- and every predator - return (if c `elem` (resistance p) -- if the plant cannot deter the predator - then prob * fitnessImpact p -- impact it weighted by probability - else 0) - where - cs = (zip (toEnum <$> [1..]) (toList cs')) :: [(Compound,Amount)] - +deterPredators :: Vector Amount -> World Probability +deterPredators cs = do + ps <- predators <$> ask + ts <- toxicCompounds <$> ask + let + deter :: Predator -> Bool + -- if any (amount of toxin t > 0| for all toxins t | t not in p's resistance-list) + deter p = or [cs ! (fromEnum t) > 0 | t <- ts, not (t `elem` resistance p)] + -- multipy (probability of occurence * intensity of destruction| for all predators, if not deterred by toxins + return . product $ [prob*fitnessImpact p | (p,prob) <- ps, not $ deter p] -- Mating & Creation of diversity -- ------------------------------ @@ -214,11 +240,16 @@ deterPredators ps cs' = sum $ do -- Running the simulation -- ---------------------- + + + 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_ print fitness -- Utility Functions -- -----------------