chemodiversity/app/Main.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE OverloadedStrings #-}
module Main where

import Text.Printf
import Control.Monad.Reader
import Numeric.LinearAlgebra
import Data.List
import System.Random
import Control.Concurrent
import Control.Parallel.Strategies
import Control.Monad.Writer
import qualified Debug.Trace as Debug
import qualified Control.Foldl as F
import System.IO
import Data.Aeson
import qualified Data.ByteString as BS

import ArbitraryEnzymeTree
import Environment
import Evaluation

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

-- Environment

exampleEnvironment :: Int -> [Enzyme] -> [(Predator,Probability)] -> [(Compound,Amount)] -> Environment
exampleEnvironment addedC es pred tox =
  Environment
    { soil = [ (PPM, 10)
             ]
    , predators = pred -- [ (greenfly, 0.1) ]
    , metabolismIteration = 100
    , maxCompound = maxCompoundWithoutGeneric + addedC
    , toxicCompounds = tox --[(Produced FPP,0.1)] ++ tox
    , possibleEnzymes = es -- [pps,fpps] ++ es
    , settings = Settings { automimicry = False
                          , predatorsRandom = False
                          , numPlants = 50
                          }
    }

-- 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 = fmap ( limit Phosphor 2
--                             . limit Nitrate 1
--                             . limit Sulfur 0
--                             ) <$> fromEnv soil
--    -- 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] -> Simulation -> IO ()
loop loopAmount ps env = loop' loopAmount 0 ps env

  where
    -- cache enzyme colorful-strings
    stringe :: [(Enzyme, String)]
    stringe = (\e -> case Data.List.find (\(t,_) -> (t==) . fst . snd . synthesis $ e) toxins of
               Just (_,toxicity) -> (e,"\ESC[38;5;" ++ show (16 + 36*5 + 6*floor (5*(1-toxicity)) + 0) ++ "m" -- yellow -> red rainbow for tocixity 0 -> 1
                                       ++ padded 50 (show (enzymeName e)) ++ "\ESC[0m")
               Nothing           -> (e, padded 50 (show (enzymeName e)))
             ) <$> possibleEnzymes (snd env)
    toxins :: [(Compound, Amount)]
    toxins = toxicCompounds (snd env)
    padded i str = take i $ str ++ repeat ' '
    printEvery = 10
    loop' :: Int -> Int -> [Plant] -> Simulation -> IO ()
    loop' loopAmount curLoop plants s = unless (loopAmount+1 == curLoop) $ do
      when (curLoop `mod` printEvery == 0) $ do
        putStr "\ESC[2J\ESC[H"
        printEnvironment (snd env)
        putStrLn ""
        putStrLn $ "Generation " ++ show curLoop ++ " of " ++ show loopAmount ++ ":"
      newPlants <- simulate s $ do
          (!fs,cs) <- unzip <$> fitness plants
          let fps = zip plants fs -- gives us plants & their fitness in a tuple
              sumFitness = sum fs
              -- $C_{\Sigma,mu}$: Durchschnittliche Menge an produzierten Stoffen
              -- $C_{\Sigma,sigma}$: Durchschnittliche Varianz an produzierten Stoffen
              (c_Sigma_mu, c_Sigma_sigma) = meanAndVar `from` sumProducedCompounds $ cs
              -- - $C_{i,\mu}$: Durchschnittliche Anzahl produzierter Komponenten
              -- - $C_{i,\sigma}$: Zusätzlich: Betrachtung der Varianz dieser Komponenten innerhalb der Population  
              --   (Z.B. Stoff A wird immer mit $0.5$ produziert, hat also keine Varianz,
              --   wogegen Stoff B *im Schnitt* mit $0.5$ produziert wird, aber dies eine extreme
              --   Varianz auslöst)
              --   c_mu_sigma :: [(Amount,Amount)]
              (c_i_mu,c_i_sigma) = unzip $ meanAndVar `from` id <$> byProducts cs
              -- - $C_d$: Durchschnittliche Anzahl distinkter Produzierter Stoffe (sprich
              --   nicht-endemisch, $#i | C_{i,\sigma} > \epsilon$ )
              isEndemic :: Vector Bool
              isEndemic = (> 0.01) <$> fromList c_i_sigma
              (c_d_mu, c_d_sigma) = meanAndVar `from` (countWith isEndemic) $ cs
              -- - $C_{\sigma,\{\mu/\sigma\}}$: Mittelwert/Varianz von $\C_{i,\sigma}$
              (c_sigma_mu, c_sigma_sigma) = meanAndVar `from` id $ c_i_sigma
              -- - $\mathbf{E}[C_{\Sigma,plant} - C_{\Sigma,mu}]$: Durchschnittliche Abweichung der produzierten
              --   Stoffe gegenüber dem Schnitt der Gesamtpopulation
              e_hash_plant = F.mean `from` numDistinctCompounds $ cs
              -- mean and variance of fitness
              fns = meanAndVar `from` id $ fs
              -- - $P_\{\mu,\sigma\}$ Mittelwert/Varianz der Anteile der Stoffe in Pflanze i, die giftig sind
              toxs = meanAndVar `from` percentagePoisonous $ cs
          when (curLoop `mod` printEvery == 0) $ liftIO $ do
            printPopulation stringe (zip3 plants fs cs)
            putStrLn $ "Population statistics           (mean,variance):"
            putStrLn $ "Amount of Components produced = " ++ (padded 50 . show $ spc)
            putStrLn $ "Number of distinct Components = " ++ (padded 50 . show $ ndc)
            putStrLn $ "Fitness                       = " ++ (padded 50 . show $ fns)
            putStrLn $ show cstats
            putStrLn $ show cstats'
            hFlush stdout
            threadDelay $ 10*1000 -- sleep x*1000ns (=x ~ ms)
          tell $         show curLoop
               ++ "," ++ show (fst spc) ++ "," ++ show (snd spc)
               ++ "," ++ show (fst ndc) ++ "," ++ show (snd ndc)
               ++ "," ++ show (fst fns) ++ "," ++ show (snd fns)
          -- generate x new plants.
          np <- fromEnv (numPlants . settings)
          sequence . flip fmap [1..np] $ \_ -> 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
      loop' loopAmount (curLoop+1) newPlants s

main :: IO ()
main = do
  hSetBuffering stdin NoBuffering
  --hSetBuffering stdout NoBuffering
  randomCompounds <- makeHead (Substrate PPM) <$> generateTreeFromList 30 (toEnum <$> [(maxCompoundWithoutGeneric+1)..] :: [Compound]) -- generate roughly x compounds
  ds <- randoms <$> newStdGen
  --probs <- randomRs (0.2,0.7) <$> newStdGen
  let poisonedTree     = poisonTree ds randomCompounds
      poisonCompounds  = foldMap (\(a,b) -> [(b,a) | a > 0]) poisonedTree
  predators <- generatePredators 0.0 poisonedTree
  let poisonCompounds' = pruneCompounds poisonCompounds predators
      pruneCompounds cs ps = filter ((`elem` usedPoisons) . fst) cs
        where usedPoisons = concat $ irresistance <$> ps
  --let env              = exampleEnvironment (getTreeSize randomCompounds) (generateEnzymeFromTree randomCompounds) (zip predators probs) poisonCompounds'
  (Just env) <- decodeStrict' <$> BS.readFile "environment2.json"
  let emptyPlants      = replicate (numPlants . settings $ env) emptyPlant
  enzs <- randomRs (0,length (possibleEnzymes env) - 1) <$> newStdGen
  let startPlants      = randomGenome 1 enzs (possibleEnzymes env) emptyPlants
  printEnvironment env
  writeFile "poison.twopi" $ generateDotFromPoisonTree "poison" 0.5 poisonedTree
  --writeFile "environment.json" . encode $ env
  putStr "\ESC[?1049h"
  logfile <- openFile "simulation.log" WriteMode
  loop 2000 startPlants (logfile,env)
  putStrLn "Simulation ended. Press key to exit."
  _ <- getChar
  putStr "\ESC[?1049l"

randomGenome :: Int -> [Int] -> [Enzyme] -> [Plant] -> [Plant]
randomGenome num inds enzs []     = []
randomGenome num inds enzs (p:ps) = p { genome = genes} : randomGenome num r enzs ps
  where
    i' = take num inds
    r  = drop num inds
    enzymes = (enzs!!) <$> i'
    genes = (\e -> (e,1,1)) <$> enzymes


generatePredators :: Double -> EnzymeTree s (Double,Compound) -> IO [Predator]
generatePredators threshold t = do
    ps <- mapM generatePredators' $ getSubTrees t
    return $ filter ((/= []) . irresistance) $ concat ps -- filter out predators that are resistant to everything because this does not make sense in our model.
  where
      generatePredators' :: EnzymeTree s (Double, Compound) -> IO [Predator]
      generatePredators' t = do -- not fully resistant to t, but fully resistant to everything in ts
          let comps = foldMap (\(a,b) -> [(a,b) | a > threshold]) t
          amount <- randomRIO (0,length comps + 1) :: IO Int
          forM [1..amount] $ \_ -> do
              impact <- randomRIO (0.2,0.7)
              rands <- randoms <$> newStdGen
              let unresists = foldMap (\((a,b),r) -> [b | r*2 < a]) $ zip comps rands
              return $ Predator unresists impact 1

printEnvironment :: Environment -> IO ()
printEnvironment (Environment soil pred metaIter maxComp toxic possEnz settings) =
  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
    putStrLn $ "Settings:  " ++ show settings

printPopulation :: [(Enzyme,String)] -> [(Plant,Double,Vector Amount)] -> IO ()
printPopulation es ps = do
  let padded i str = take i $ str ++ repeat ' '
      n = length ps
      fitnesses = (\(_,f,_) -> f) <$> ps
      meanFitness = sum fitnesses / fromIntegral n
      maxFitness = maximum fitnesses
  putStr $ padded 50 ("Population: (fitness: mean " ++ padded 5 (show meanFitness) ++ ", max: " ++ padded 5 (show maxFitness) ++ ")")
  forM_ ps $ \(_,f,_) -> putStr (printColor (f/maxFitness) '█')
  putStrLn colorOff
  forM_ es $ \(e,s) -> do
    putStr s
    forM_ ps $ \(Plant g _,_,cs) -> 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 = '_'
          amount = min 2 $ cs ! fromEnum (fst . snd . synthesis $ e)
      putStr $ printColor (amount/2) (plot curE)
    putStrLn colorOff

printColor :: Double -> Char -> String
printColor x c
  | x > 1 = "Error: " ++ show x
  | 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"

generateEnzymeFromTree :: EnzymeTree s Compound -> [Enzyme]
generateEnzymeFromTree t = (makeSimpleEnzyme c . getElement <$> sts)
                           ++ concatMap generateEnzymeFromTree sts
  where
    c   = getElement t
    sts = getSubTrees t


stepDebug a = liftIO $ print a >> void getChar