working evolution-simulation (POC)

app/Main.hs

@@ -1,6 +1,172 @@
+{-# LANGUAGE TypeApplications #-}
+
 module Main where
 
-import Lib
+import Environment
+import Text.Printf
+import Control.Monad.Reader
+import Numeric.LinearAlgebra
+import Data.List
+import System.Random
+import Control.Concurrent
+import qualified Debug.Trace as Debug
+import System.IO
+
+
+-- Example definitions
+-- -------------------
+
+-- Enzymes
+
+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)
+
+-- Predator
+
+greenfly :: Predator         -- 20% of plants die to greenfly, but the fly is
+greenfly = Predator [] 0.2   -- killed by any toxic Component
+
+-- Environment
+
+exampleEnvironment :: Environment
+exampleEnvironment =
+  Environment
+    { soil = [ (Nitrate, 2)
+             , (Phosphor, 3)
+             , (Photosynthesis, 10)
+             ]
+    , predators = [ (greenfly, 0.1) ]
+    , metabolismIteration = 100
+    , maxCompound = maxCompoundWithoutGeneric + 100
+    , toxicCompounds = [(Produced FPP,0.5)] --FPP kills 100% if produced amount above 0.2 units
+    , possibleEnzymes = [pps,fpps]
+    }
+
+-- Plants
+
+examplePlants :: [Plant]
+examplePlants = (\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 = soil <$> ask >>= return . fmap ( limit Phosphor 2
+                                                      . limit Nitrate 1
+                                                      . limit Sulfur 0
+                                                      )
+   -- 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')
+
+-- Running the simulation
+-- ----------------------
+
+loop :: Int -> [Plant] -> Environment -> IO ()
+loop loopAmount plants e = loop' loopAmount 0 plants e
+
+  where
+    loop' :: Int -> Int -> [Plant] -> Environment -> IO ()
+    loop' loopAmount curLoop plants e = unless (loopAmount == curLoop) $ do
+      putStr $ "\ESC[2J\ESC[H"
+      printEnvironment e
+      putStrLn ""
+      putStrLn $ "Generation " ++ show curLoop ++ " of " ++ show loopAmount ++ ":"
+      newPlants <- (flip runReaderT) e $ do
+          fs <- sequence $ fitness <$> plants
+          let fps = zip plants fs -- gives us plants & their fitness in a tuple
+              sumFitness = sum fs
+          pe <- possibleEnzymes <$> ask
+          liftIO $ printPopulation pe fps
+          -- generate 100 new plants.
+          sequence . (flip fmap) [1..100] $ \_ -> do
+                parent' <- liftIO $ randomRIO (0,sumFitness)
+                let
+                    -- if we only have one parent in our list, take it.
+                    findParent :: Double -> [(Plant,Double)] -> Plant
+                    findParent _ [(last,_)] = last
+                    -- otherwise count down x to find the parent in the list
+                    findParent x ((p,f):ps)
+                      | x < f     = p
+                      | otherwise = findParent (x-f) ps
+                    parent = findParent parent' fps
+                haploMate parent
+      hFlush stdout
+      threadDelay $ 100*1000 -- sleep 100ms
+      loop' loopAmount (curLoop+1) newPlants e
 
 main :: IO ()
-main = someFunc
+main = do
+  hSetBuffering stdin NoBuffering
+  hSetBuffering stdout NoBuffering
+  let emptyPlants = replicate 100 emptyPlant
+  printEnvironment exampleEnvironment
+  putStr "\ESC[?1049h"
+  loop 100 emptyPlants exampleEnvironment
+  putStrLn "Simulation ended. Press key to exit."
+  _ <- getChar
+  putStr "\ESC[?1049l"
+
+  -- fitness <- runReaderT (sequence $ (\a -> do p <- absorbNutrients a >>= produceCompounds a; (,,) a p <$> deterPredators p) <$> emptyPlants) exampleEnvironment
+  -- mapM_ (printf "%15.15s, " . show . toEnum @Compound) [0..maxCompoundWithoutGeneric]
+  -- putStrLn "Fitness"
+  -- forM_ fitness $ \(p, c, f) -> do
+  --   mapM_ (printf "%15.2f, ") (toList c)
+  --   printf "%15.2f" f
+  --   putStr "\n"
+
+printEnvironment :: Environment -> IO ()
+printEnvironment (Environment soil pred metaIter maxComp toxic possEnz) =
+  do
+    putStrLn "Environment:"
+    putStrLn $ "Soil:      " ++ show soil
+    putStrLn $ "Predators: " ++ show pred
+    putStrLn $ "PSM Iters: " ++ show metaIter
+    putStrLn $ "Compounds: " ++ show ((toEnum <$> [0..maxComp]) :: [Compound])
+    putStrLn $ "Toxic:     " ++ show toxic
+
+printPopulation :: [Enzyme] -> [(Plant,Double)] -> IO ()
+printPopulation es ps = do
+  let padded i str = take i $ str ++ repeat ' '
+  putStr $ padded 40 "Population:"
+  forM_ ps $ \((_,f)) -> putStr (printColor f '█')
+  putStrLn colorOff
+  forM_ es $ \e -> do
+    putStr $ padded 40 (show (enzymeName 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 ""
+
+printColor :: Double -> Char -> String
+printColor x c
+  | x*x < 0.5 = "\ESC[38;5;" ++ (show $ 16 + 36*5 + 6*(floor $ 5*2*x') + 0) ++ "m" ++ [c] ++ ""
+  | otherwise = "\ESC[38;5;" ++ (show $ 16 + 36*(floor $ 5*2*(1-x')) + 6*5 + 0) ++ "m" ++ [c] ++ ""
+  -- 32 bit
+  -- | x*x < 0.5 = "\ESC[38;2;255;" ++ (show . floor $ 255*2*x') ++ ";0m" ++ [c] ++ ""
+  -- | otherwise = "\ESC[38;2;" ++ (show . floor $ 255*2*(1-x')) ++ ";255;0m" ++ [c] ++ ""
+   where x' = x*x
+
+colorOff :: String
+colorOff = "\ESC[0m"

package.yaml

@@ -23,6 +23,7 @@ dependencies:
 - base >= 4.7 && < 5
 - hmatrix
 - mtl
+- random
 
 library:
   source-dirs: src

src/Environment.hs

@@ -11,7 +11,7 @@ import Control.Monad.Reader
 import Data.List           (permutations, subsequences)
 import Numeric.LinearAlgebra
 import Text.Printf
-import qualified Debug.Trace as Debug
+import System.Random
 
 type Probability = Double
 type Quantity = Int
@@ -76,14 +76,6 @@ data Enzyme = Enzyme
 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
 -- ----------
 
@@ -97,11 +89,6 @@ data Predator = Predator { resistance    :: [Compound]
                            --   (~ agressiveness of the herbivore)
                          } deriving (Show, Eq)
 
--- Exemplatory:
-
-greenfly :: Predator         -- 20% of plants die to greenfly, but the fly is
-greenfly = Predator [] 0.2   -- killed by any toxic Component
-
 -- The environment itself is just the soil and the predators. Extensions would be possible.
 
 data Environment =
@@ -122,27 +109,14 @@ data Environment =
    , toxicCompounds :: [(Compound,Amount)]
      -- ^ Compounds considered to be toxic in this environment.
      --   Kills 100% of Predators above Amount.
+   , possibleEnzymes :: [Enzyme]
+     -- ^ All enzymes that can be created by genetic manipulation in this setting.
    } 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
-exampleEnvironment =
-  Environment
-    { soil = [ (Nitrate, 2)
-             , (Phosphor, 3)
-             , (Photosynthesis, 10)
-             ]
-    , predators = [ (greenfly, 0.1) ]
-    , metabolismIteration = 100
-    , maxCompound = maxCompoundWithoutGeneric
-    , toxicCompounds = [(Produced FPP,0.5)] --FPP kills 100% if produced amount above 0.2 units
-    }
-
 type World a = ReaderT Environment IO a
 
 -- Plants
@@ -164,35 +138,6 @@ instance Show Plant where
 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 = soil <$> ask >>= return . fmap ( limit Phosphor 2
-                                                      . limit Nitrate 1
-                                                      . limit Sulfur 0
-                                                      )
-   -- 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
 -- -------
@@ -208,7 +153,9 @@ fitness p = do
   nutrients <- absorbNutrients p           -- absorb soil
   products <- produceCompounds p nutrients -- produce compounds
   survivalRate <- deterPredators products  -- defeat predators with produced compounds
-  return survivalRate
+  let sumEnzymes = sum $ (\(_,q,a) -> (fromIntegral q)*a) <$> genome p -- amount of enzymes * activation = resources "wasted"
+      costOfEnzymes = 0.95 ** sumEnzymes
+  return $ survivalRate * costOfEnzymes
   -- can also be written as, but above is more clear.
   -- fitness p = absorbNutrients p >>= produceCompounds p >>= deterPredators
 
@@ -239,52 +186,64 @@ deterPredators cs = do
          -- 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]
+   return . product $ [min 1 ((1-prob) * fitnessImpact p / deter p) | (p,prob) <- ps]
 
 -- Mating & Creation of diversity
 -- ------------------------------
 
--- Running the simulation
--- ----------------------
+
+-- | mate haploid
+haploMate :: Plant -> World Plant
+haploMate (Plant genes abs) = do
+  --generate some random infinite uniform distributed lists of doubles in [0,1)
+  r1 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r2 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r3 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r4 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r5 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  enzymes <- possibleEnzymes <$> ask
+  re1 <- liftIO ((randomRs (0,length enzymes - 1) <$> newStdGen) :: IO [Int])
+  re2 <- liftIO ((randomRs (0,length enzymes - 1) <$> newStdGen) :: IO [Int])
+  let
+    genes' = mutateGene r1 re1
+           . noiseActivation r2
+           . addGene r3 re2
+           . duplicateGene r4
+           . deleteGene r5
+           $ genes
+    deleteGene :: [Double] -> Genome -> Genome
+    deleteGene (r:rs) ((e,1,a):gs) = if a < 0.1 && r < 0.5 then           deleteGene rs gs else (e,1,a):deleteGene rs gs
+    deleteGene (r:rs) ((e,q,a):gs) = if a < 0.1 && r < 0.5 then (e,q-1,a):deleteGene rs gs else (e,q,a):deleteGene rs gs
+    deleteGene _ []                = []
+
+    duplicateGene :: [Double] -> Genome -> Genome
+    duplicateGene (r:rs) ((e,q,a):gs) = if r < 0.05 then (e,q+1,a):duplicateGene rs gs else (e,q,a):duplicateGene rs gs
+    duplicateGene _ []                = []
+
+    addGene :: [Double] -> [Int] -> Genome -> Genome
+    addGene (r:rs) (s:ss) g = if r < 0.01 then ((enzymes !! s),1,1):g else g
+
+    noiseActivation :: [Double] -> Genome -> Genome
+    noiseActivation (r:rs) ((e,q,a):gs) = (e,q,max 0 $ min 1 $ a-0.01+0.02*r):noiseActivation rs gs
+    noiseActivation _ []                = []
+
+    mutateGene :: [Double] -> [Int] -> Genome -> Genome
+    mutateGene (r:rs) (s:ss) ((e,1,a):gs) = if r < 0.05 then ((enzymes !! s),1,a):mutateGene rs ss gs
+                                                        else (e,1,a):mutateGene rs ss gs
+
+    mutateGene (r:rs) (s:ss) ((e,q,a):gs) = if r < 0.05 then (e,q-1,a):((enzymes !! s),1,a):mutateGene rs ss gs
+                                                        else (e,q,a):mutateGene rs ss gs
+    mutateGene (r:rs) (s:ss) []           = []
+  return $ Plant genes' abs
 
 
 
-
-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_ (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
 -- -----------------
 
+-- | Plant with no secondary metabolism with unlimited extraction from environment.
+emptyPlant :: Plant
+emptyPlant = Plant [] (soil <$> ask)
+
 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 40 (show (enzymeName 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 ""