PoC for single iteration.

package.yaml

@@ -22,6 +22,7 @@ description:         Please see the README on GitHub at <https://github.com/Drez
 dependencies:
 - base >= 4.7 && < 5
 - hmatrix
+- mtl
 
 library:
   source-dirs: src

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
+   defaultAbsorption = soil <$> ask >>= return . fmap ( 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')
@@ -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
 -- -----------------