From 08b2cf8d43481749dce681d6a3ea716e0829c1f2 Mon Sep 17 00:00:00 2001 From: Stefan Dresselhaus Date: Wed, 2 May 2018 18:26:37 +0200 Subject: [PATCH] iteration of plants-metabolism works and finds the right limits. --- src/Environment.hs | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/Environment.hs b/src/Environment.hs index 078d7d1..ddf4137 100644 --- a/src/Environment.hs +++ b/src/Environment.hs @@ -1,4 +1,5 @@ {-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE TypeApplications #-} module Environment where @@ -9,6 +10,8 @@ import Control.Monad (forM_) import Control.Monad.Reader import Data.List (permutations, subsequences) import Numeric.LinearAlgebra +import Text.Printf +import qualified Debug.Trace as Debug type Probability = Double type Quantity = Int @@ -116,8 +119,9 @@ data Environment = -- -- To find the 'maxCompound' without 'GenericEnzyme' use -- 'maxComponent = fromEnum (maxBound :: Nutrient) + fromEnum (maxBound :: Component) + 1' - , toxicCompounds :: [Compound] + , toxicCompounds :: [(Compound,Amount)] -- ^ Compounds considered to be toxic in this environment. + -- Kills 100% of Predators above Amount. } deriving (Show, Eq) -- helper function. Allows for [0..maxCompoundWithoutGeneric] :: [Compound] with all non-generic Compounds @@ -136,7 +140,7 @@ exampleEnvironment = , predators = [ (greenfly, 0.1) ] , metabolismIteration = 100 , maxCompound = maxCompoundWithoutGeneric - , toxicCompounds = [Produced FPP] + , toxicCompounds = [(Produced FPP,0.5)] --FPP kills 100% if produced amount above 0.2 units } type World a = ReaderT Environment IO a @@ -213,13 +217,16 @@ produceCompounds (Plant genes _) substrate = do numIter <- metabolismIteration <$> ask numCompounds <- maxCompound <$> ask let - initialAmount = assoc (numCompounds+1) 0 ((\(n,a) -> (fromEnum $ Substrate n,a)) <$> substrate) + initialAmount = (assoc (numCompounds+1) 0 ((\(n,a) -> (fromEnum $ Substrate n,a)) <$> substrate)) :: Vector Amount 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),2*f*ia),((fromEnum o,fromEnum i),oa),((fromEnum i,fromEnum o),f*oa)]) <$> enzymes -- [((row,column),amount)] + positions = concat $ (\(((i,ia),(o,oa)),f) -> [((fromEnum i,fromEnum i),f*ia),((fromEnum o,fromEnum o),f*ia),((fromEnum o,fromEnum i),f*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. + -- (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) <> inv v)) #> initialAmount -- (E + A)^numIter * t_0 for numIter iterations. + final = (realPart `cmap` matFunc (^numIter) (ident (numCompounds+1) + (((*0.01) . (:+ 0)) `cmap` mat))) #> initialAmount + -- matFunc splits mat into UD(U^-1), applies function to diag-Elements in D, then multiplies togehter. + -- faster, because no inversions and optimized eig. return final @@ -228,11 +235,11 @@ 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] + 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, 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] -- Mating & Creation of diversity -- ------------------------------ @@ -249,7 +256,12 @@ main = do 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 + 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 -- ----------------- @@ -262,7 +274,7 @@ printPopulation es ps = do let padded i str = take i $ str ++ repeat ' ' putStrLn "Population:" forM_ es $ \e -> do - putStr $ padded 8 (show e) + putStr $ padded 40 (show (enzymeName e)) forM_ ps $ \(Plant g _) -> do let curE = sum $ map (\(_,q,a) -> (fromIntegral q)*a) . filter (\(e',_,_) -> e == e')