From cc80892147ded4d9da8b37376381768c7d8f93a5 Mon Sep 17 00:00:00 2001 From: Stefan Dresselhaus Date: Wed, 2 May 2018 00:54:15 +0200 Subject: [PATCH] ported to stack-project, iterated fitness. --- package.yaml | 1 + src/Environment.hs | 247 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 248 insertions(+) create mode 100644 src/Environment.hs diff --git a/package.yaml b/package.yaml index 3ebe921..3ca7cbd 100644 --- a/package.yaml +++ b/package.yaml @@ -21,6 +21,7 @@ description: Please see the README on GitHub at = 4.7 && < 5 +- hmatrix library: source-dirs: src diff --git a/src/Environment.hs b/src/Environment.hs new file mode 100644 index 0000000..9985dc6 --- /dev/null +++ b/src/Environment.hs @@ -0,0 +1,247 @@ +{-# LANGUAGE RecordWildCards #-} + +module Environment where + + +import Data.Functor ((<$>)) +import Control.Applicative ((<*>)) +import Control.Monad (forM_) +import Data.List (permutations, subsequences) +import Numeric.LinearAlgebra + +type Probability = Double +type Quantity = Int +type Activation = Double +type Amount = Double + +-- | Nutrients are the basis for any reaction and are found in the environment of the plant. +data Nutrient = Sulfur + | Phosphor + | Nitrate + | Photosynthesis + deriving (Show, Enum, Bounded, Eq) + +-- | Fixed, non-generic Components +data Component = PP + | FPP + deriving (Show, Enum, Bounded, Eq) + +-- | Compounds are either direct nutrients, already processed components or GenericEnzymes +data Compound = Substrate Nutrient + | Produced Component + | GenericEnzyme Int + deriving (Show, Eq) + +instance Enum Compound where + toEnum x + | x <= maxS = Substrate . toEnum $ x + | x - (maxS+1) <= maxP = Produced . toEnum $ x - (maxS + 1) + | otherwise = GenericEnzyme $ x - (maxS + 1) - (maxP + 1) + where + maxS = fromEnum (maxBound :: Nutrient) + maxP = fromEnum (maxBound :: Component) + + fromEnum (Substrate x) = fromEnum x + fromEnum (Produced x) = fromEnum x + maxS + 1 + where + maxS = fromEnum (maxBound :: Nutrient) + fromEnum (GenericEnzyme x) = x + maxS + maxP + 2 + where + maxS = fromEnum (maxBound :: Nutrient) + maxP = fromEnum (maxBound :: Component) + + + +-- | Enzymes are the main reaction-driver behind synthesis of intricate compounds. +-- +-- They are assumed to be reversible. +data Enzyme = Enzyme + { enzymeName :: String + -- ^ Name of the Enzyme. + , substrateRequirements :: [(Compound,Amount)] + -- ^ needed for reaction to take place + , synthesis :: ((Compound,Amount),(Compound,Amount)) + -- ^ given x in amount -a, this will produce y in amount b + , dominance :: Maybe Amount + -- ^ in case of competition for nutrients this denotes the priority + -- Nothing = max possible + } + deriving (Show, Eq) + +-- | conviniently make an Enzyme using 1 of the first compund to produce 1 of the second +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 +-- ---------- + +-- | 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 + , fitnessImpact :: Amount + -- ^ impact on the fitness of a plant + -- (~ agressiveness of the herbivore) + } deriving (Show, Eq) + +-- 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 + +-- The environment itself is just the soil and the predators. Extensions would be possible. + +data Environment = + Environment + { soil :: [(Nutrient, Amount)] + -- ^ soil is a list of nutrients available to the plant. + , predators :: [(Predator, Probability)] + -- ^ Predators with the probability of appearance in this generation. + } deriving (Show, Eq) + +-- Example: + +exampleEnvironment :: Environment +exampleEnvironment = + Environment + { soil = [ (Nitrate, 2) + , (Phosphor, 3) + , (Photosynthesis, 10) + ] + , predators = [ (greenfly, 0.1) ] + } + +-- Plants +-- ------ + +-- Plants consist of a Genome responsible for creation of the PSM and also an +-- external state how many nutrients and compounds are currently inside the plant. + +type Genome = [(Enzyme, Quantity, Activation)] + +data Plant = Plant + { genome :: Genome + -- ^ the genetic characteristic of the plant + , absorbNutrients :: Environment -> [(Nutrient,Amount)] + -- ^ the capability to absorb nutrients given an environment + } +instance Show Plant where + show p = "Plant with Genome " ++ show (genome p) +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 (Environment s _) = limit Phosphor 2 + . limit Nitrate 1 + . limit Sulfur 0 + <$> s + -- 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 +-- ------- + +-- The fitness-measure is central for the generation of offspring and the +-- simulation. It evaluates the probability for passing on genes given a plant in +-- an environment. + +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 + + +produceCompounds :: Plant -> [(Nutrient, Amount)] -> Vector Amount +produceCompounds (Plant genes _) substrate = final + where + initialAmount = assoc 10 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 + + +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)] + + +-- Mating & Creation of diversity +-- ------------------------------ + +-- Running the simulation +-- ---------------------- + +main = do + putStrLn "Environment:" + print exampleEnvironment + putStrLn "Example population:" + printPopulation [pps, fpps] plants + +-- Utility Functions +-- ----------------- + +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 8 (show 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 ""