chemodiversity/src/Environment.hs

{-# 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 ""