iteration of plants-metabolism works and finds the right limits.

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')