corrected many things according to specification that was agreed upon.

app/Main.hs

@@ -4,17 +4,20 @@ module Main where
 
 import Text.Printf
 import Control.Monad.Reader
-import Numeric.LinearAlgebra
+import qualified Numeric.LinearAlgebra as LA
 import Data.List
 import System.Random
 import Control.Concurrent
 import Control.Parallel.Strategies
-import Control.Monad.Writer
+import Control.Monad.Writer (tell)
 import qualified Debug.Trace as Debug
 import qualified Control.Foldl as F
 import System.IO
+import System.Environment
 import Data.Aeson
 import qualified Data.ByteString as BS
+import Options.Applicative
+import Data.Semigroup ((<>))
 
 import ArbitraryEnzymeTree
 import Environment
@@ -46,6 +49,8 @@ exampleEnvironment addedC es pred tox =
     , settings = Settings { automimicry = False
                           , predatorsRandom = False
                           , numPlants = 50
+                          , logEveryNIterations = 10
+                          , verbose = True
                           }
     }
 
@@ -91,20 +96,32 @@ loop loopAmount ps env = loop' loopAmount 0 ps env
              ) <$> possibleEnzymes (snd env)
     toxins :: [(Compound, Amount)]
     toxins = toxicCompounds (snd env)
+    printEverything = verbose.settings.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
+      when (printEverything && 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
+          when (curLoop == 0) $
+              tell $  "num_iter"
+                   ++ ",c_sum_mu,c_sum_sigma"
+                   ++ ",c_d_mu,c_d_sigma"
+                   ++ ",e_d_mu,e_d_sigma"
+                   ++ ",fitness_mean,fitness_sigma"
+                   ++ ",percent_toxic_mean,percent_toxic_sigma"
+          logIter <- fromEnv $ logEveryNIterations . settings
           (!fs,cs) <- unzip <$> fitness plants
           txns <- fmap (fromEnum . fst) <$> fromEnv toxicCompounds -- [Int] of id's of toxic compounds
           let fps = zip plants fs -- gives us plants & their fitness in a tuple
               sumFitness = sum fs
+              es = genomeToEnzymeAmount . genome <$> plants
+              genomeToEnzymeAmount :: Genome -> LA.Vector Double
+              genomeToEnzymeAmount g = LA.accum (LA.konst 0 (maxCompound . snd $ env)) (+) $ (\(e,q,a) -> ((fromEnum . fst . snd . synthesis $ e)-1,fromIntegral q*a)) <$> g
               -- $C_{\Sigma,mu}$: Durchschnittliche Menge an produzierten Stoffen
               -- $C_{\Sigma,sigma}$: Durchschnittliche Varianz an produzierten Stoffen
               (c_sum_mu, c_sum_sigma) = meanAndVar `from` sumProducedCompounds $ cs
@@ -114,13 +131,22 @@ loop loopAmount ps env = loop' loopAmount 0 ps env
               --   wogegen Stoff B *im Schnitt* mit $0.5$ produziert wird, aber dies eine extreme
               --   Varianz auslöst)
               (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 = fromList $ (> 0.01) <$> 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
+              -- - $C_d$: Durchschnittliche Anzahl distinkter Produzierter Stoffe (sprich
+              --   nicht-endemisch, $#i | C_{i,\mu} < \epsilon$ )
+              isNotEndemicCompound :: LA.Vector Bool
+              isNotEndemicCompound = LA.fromList $ (< 0.1) <$> c_i_mu
+              (c_d_mu, c_d_sigma) = meanAndVar `from` countWith isNotEndemicCompound (>0.1) $ cs
+              -- - $E_{i,\mu}$: Durchschnittliche Anzahl produzierbarer Komponenten (falls ausgangsstoff verfügbar)
+              -- - $E_{i,\sigma}$: Zusätzlich: Betrachtung der Varianz dieser Komponenten innerhalb der Population
+              -- analog zu $C_{i,\mu/\sigma}$
+              (e_i_mu,e_i_sigma) = unzip $ meanAndVar `from` id <$> byCompound es
+              -- - $E_d$: Durchschnittliche Anzahl distinkter Produzierter Stoffe (sprich
+              --   nicht-endemisch, $#i | E_{i,\mu} < \epsilon$ )
+              isNotEndemicEnzyme :: LA.Vector Bool
+              isNotEndemicEnzyme = LA.fromList $ (< 0.5) <$> e_i_mu
+              (e_d_mu, e_d_sigma) = meanAndVar `from` countWith isNotEndemicEnzyme (>0.5) $ es
               -- - $\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
@@ -128,17 +154,21 @@ loop loopAmount ps env = loop' loopAmount 0 ps env
               fns = meanAndVar `from` id $ fs
               -- - $P_\{\mu,\sigma\}$ Mittelwert/Varianz der Anteile der Stoffe in Pflanze i, die giftig sind
               toxs = meanAndVar `from` percentagePoisonous txns $ cs
-          when (curLoop `mod` printEvery == 0) $ liftIO $ do
+          when (printEverything && curLoop `mod` printEvery == 0) $ liftIO $ do
-            printPopulation (zip ((>0.01) <$> c_i_sigma) stringe) (zip3 plants fs cs)
+            printPopulation isNotEndemicEnzyme stringe (zip3 plants fs cs)
             putStrLn $ "Population statistics           (mean,variance):"
-            putStrLn $ "Amount of Components produced = " ++ (padded 50 . show $ (c_sum_mu,c_sum_sigma))
+            putStrLn $ "Amount of Components produced = " ++ show (c_sum_mu,c_sum_sigma)
-            putStrLn $ "Number of distinct Components = " ++ (padded 50 . show $ (c_d_mu, c_d_sigma))
+            putStrLn $ "Number of distinct Components = " ++ show (c_d_mu, c_d_sigma)
-            putStrLn $ "Fitness                       = " ++ (padded 50 . show $ fns)
+            putStrLn $ "Number of distinct Enzymes    = " ++ show (e_d_mu, e_d_sigma)
+            putStrLn $ "Fitness                       = " ++ show fns
+            putStrLn $ "Percentage of toxins in Cmpnds= " ++ show toxs
             hFlush stdout
             threadDelay $ 10*1000 -- sleep x*1000ns (=x ~ ms)
+          when (curLoop `mod` logIter == 0) $ 
             tell $            show curLoop
                  ++  ","  ++  show c_sum_mu    ++  ","  ++  show c_sum_sigma
                  ++  ","  ++  show c_d_mu      ++  ","  ++  show c_d_sigma
+                 ++  ","  ++  show e_d_mu      ++  ","  ++  show e_d_sigma
                  ++  ","  ++  show (fst fns)   ++  ","  ++  show (snd fns)
                  ++  ","  ++  show (fst toxs)  ++  ","  ++  show (snd toxs)
           -- generate x new plants.
@@ -157,30 +187,65 @@ loop loopAmount ps env = loop' loopAmount 0 ps env
                 haploMate parent
       loop' loopAmount (curLoop+1) newPlants s
 
+data CLIOptions = CLIOptions
+    { environment :: Maybe FilePath
+    , logfile     :: FilePath
+    }
+
+cliOptParser :: Parser CLIOptions
+cliOptParser = CLIOptions
+    <$> optional (strOption
+        (long "environment"
+        <> short 'e'
+        <> metavar "ENV"
+        <> help "Environment to load"
+        ))
+    <*> option str
+        (long "logfile"
+        <> short 'l'
+        <> metavar "LOG"
+        <> showDefault
+        <> value "simulation.log"
+        <> help "Name for the logfile"
+        )
+
+cliopts = info (cliOptParser <**> helper)
+                (fullDesc
+                <> progDesc "Simulation of Biological Systems"
+                <> header "Chemodiversity made easy ;)"
+                )
+
+
 main :: IO ()
 main = do
+  opts <- execParser cliopts
   hSetBuffering stdin NoBuffering
-  --hSetBuffering stdout 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
+  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) <- case environment opts of
-  (Just env) <- decodeStrict' <$> BS.readFile "environment2.json"
+                  Nothing   -> return . Just $ exampleEnvironment (getTreeSize randomCompounds) (generateEnzymeFromTree randomCompounds) (zip predators probs) poisonCompounds'
+                  Just file -> do
+                    putStrLn $ "reading environment: " ++ file
+                    decodeStrict' <$> BS.readFile file
   let emptyPlants      = replicate (numPlants . settings $ env) emptyPlant
+      printEverything  = verbose.settings $ env
   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 "poison.twopi" $ generateDotFromPoisonTree "poison" 0.5 poisonedTree
   --writeFile "environment.json" . encode $ env
-  putStr "\ESC[?1049h"
+  when printEverything $ putStr "\ESC[?1049h"
-  logfile <- openFile "simulation.log" WriteMode
+  loghandle <- openFile (logfile opts) WriteMode
-  loop 2000 startPlants (logfile,env)
+  putStrLn $ "logging to: " ++ logfile opts
+  loop 2000 startPlants (loghandle,env)
+  when printEverything $ do
     putStrLn "Simulation ended. Press key to exit."
     _ <- getChar
     putStr "\ESC[?1049l"
@@ -221,8 +286,8 @@ printEnvironment (Environment soil pred metaIter maxComp toxic possEnz settings)
     putStrLn $ "Toxic:     " ++ show toxic
     putStrLn $ "Settings:  " ++ show settings
 
-printPopulation :: [(Bool,(Enzyme,String))] -> [(Plant,Double,Vector Amount)] -> IO ()
+printPopulation :: LA.Vector Bool -> [(Enzyme,String)] -> [(Plant,Double,LA.Vector Amount)] -> IO ()
-printPopulation es ps = do
+printPopulation endemic es ps = do
   let padded i str = take i $ str ++ repeat ' '
       n = length ps
       fitnesses = (\(_,f,_) -> f) <$> ps
@@ -231,8 +296,9 @@ printPopulation es ps = do
   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 $ \(b,(e,s)) -> do
+  forM_ es $ \(e,s) -> do
-    if b then putStr ">" else putStr " "
+    let enzymeProductNum = fromEnum . fst . snd . synthesis $ e
+    if LA.toList endemic !! (enzymeProductNum - 1) then putStr ">" else putStr " "
     putStr s
     forM_ ps $ \(Plant g _,_,cs) -> do
       let curE = sum $ map (\(_,q,a) -> fromIntegral q*a)
@@ -245,7 +311,7 @@ printPopulation es ps = do
            | x > 0.5   = 'o'
            | x > 0     = '.'
            | otherwise = '_'
-          amount = min 2 $ cs ! fromEnum (fst . snd . synthesis $ e)
+          amount = min 2 $ cs LA.! fromEnum (fst . snd . synthesis $ e)
       putStr $ printColor (amount/2) (plot curE)
     putStrLn colorOff
 

environment2.json

@@ -490,14 +490,16 @@
         ]
     ],
     "settings": {
-        "automimicry": false,
+        "automimicry": true,
-        "numPlants": 50,
+        "logEveryNIterations": 10,
-        "predatorsRandom": false
+        "numPlants": 100,
+        "predatorsRandom": false,
+        "verbose": false
     },
     "soil": [
         [
             [],
-            10
+            20
         ]
     ],
     "toxicCompounds": [

package.yaml

@@ -30,6 +30,7 @@ dependencies:
 - foldl
 - aeson
 - bytestring
+- optparse-applicative
 
 library:
   source-dirs: src

src/Environment.hs

@@ -115,6 +115,8 @@ instance ToJSON Predator
 data Settings = Settings { automimicry :: Bool        -- ^ do we have automimicry-protection?
                          , predatorsRandom :: Bool    -- ^ do predators always appear or according to their random distribution?
                          , numPlants :: Int           -- ^ number of plants in starting population
+                         , logEveryNIterations :: Int -- ^ log status every @loopNumber `mod` logEveryNIterations == 0@
+                         , verbose :: Bool            -- ^ print visual statistics instead of just logging
                          }
                          deriving (Show, Eq, Generic)
 
@@ -137,7 +139,7 @@ data Environment =
      --   Rest will get filled up with 'GenericEnzyme i'
      --
      --   To find the 'maxCompound' without 'GenericEnzyme' use
-     --   'maxComponent = fromEnum (maxBound :: Nutrient) + fromEnum (maxBound :: Component) + 1'
+     --   @maxComponent = fromEnum (maxBound :: Nutrient) + fromEnum (maxBound :: Component) + 1@
    , toxicCompounds :: [(Compound,Amount)]
      -- ^ Compounds considered to be toxic in this environment.
      --   Kills 100% of Predators above Amount.
@@ -149,7 +151,7 @@ data Environment =
 instance FromJSON Environment
 instance ToJSON Environment
 
--- helper function. Allows for [0..maxCompoundWithoutGeneric] :: [Compound] with all non-generic Compounds
+-- helper function. Allows for @[0..maxCompoundWithoutGeneric] :: [Compound]@ with all non-generic Compounds
 maxCompoundWithoutGeneric :: Int
 maxCompoundWithoutGeneric = fromEnum (maxBound :: Nutrient) + fromEnum (maxBound :: Component) + 1
 
@@ -274,7 +276,7 @@ haploMate (Plant genes abs) = do
         is <- randoms <$> newStdGen
         return $ zip ds is
   --generate some random infinite uniform distributed lists of doubles in [0,1)
-  r1 <- liftIO digen
+  r1 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
   r2 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
   r3 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
   r4 <- liftIO digen
@@ -300,7 +302,7 @@ haploMate (Plant genes abs) = do
 
     duplicateGene :: [(Double,Int)] -> Genome -> Genome
     duplicateGene _ []                = []
-    duplicateGene ((r,i):rs) g = if r < 0.05 then duplicateGene rs (stay ++ (e,q+1,a):stay') else g
+    duplicateGene ((r,i):rs) g = if r < 0.05 then duplicateGene rs (stay ++ (e,q,a):(e,1,a):stay') else g
                                           where
                                             (stay, (e,q,a):stay') = splitAt (i `mod` length g - 2) g
 
@@ -311,14 +313,19 @@ haploMate (Plant genes abs) = do
     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)] -> [Int] -> Genome -> Genome
+    mutateGene :: [Double] -> [Int] -> Genome -> Genome
     mutateGene _ _ []                = []
-    mutateGene ((r,i):rs) (s:ss) g = if r < 0.25 then mutateGene rs ss (stay ++ go' ++ stay') else g
+    mutateGene rs ss ((e,q,a):gs) = g' ++ mutateGene rs'' ss'' gs
       where
-                                            (stay, go:stay') = splitAt (i `mod` length g - 2) g
+        -- take q randoms from rs/ss, replace numMuts (<= q) with the enzymes in ss
-                                            go' = case go of
+        (rs',rs'') = splitAt q rs
-                                                    (e,1,a) -> [(enzymes !! s,1,a)]
+        (ss',ss'') = splitAt q ss
-                                                    (e,q,a) -> [(e,q-1,a),(enzymes !! s,1,a)]
+        numMuts = length . filter (<0.01) $ rs'
+        newEnz = fmap ((\e' -> (e',1,a)).(enzymes!!).snd) 
+               . filter ((<0.01).fst) 
+               . zip rs' $ ss'
+        g' = if q == numMuts then newEnz
+                             else (e,q-numMuts,a):newEnz
   return $ Plant genes' abs
 
 

src/Evaluation.hs

@@ -38,9 +38,9 @@ numDistinctCompounds :: Functor f => f (LA.Vector Amount) -> f Amount
 --numDistinctCompounds :: [LA.Vector Amount] -> [Amount]
 numDistinctCompounds comps = sumElements . LA.cmap (\x -> if abs x < eps then 0 else 1) <$> comps
 
-countWith :: Functor f => LA.Vector Bool -> f (LA.Vector Amount) -> f Amount
+countWith :: Functor f => LA.Vector Bool -> (Amount -> Bool) -> f (LA.Vector Amount) -> f Amount
-countWith toSelect = fmap $ sumElements . LA.zipVectorWith (\selected _ -> if selected then 1 else 0) toSelect
+countWith toSelect filter = fmap $ sumElements . LA.zipVectorWith (\selected a -> if selected && filter a then 1 else 0) toSelect
---                                      apply selection to set data to 1 or 0
+--                                               apply selection & filter to set data to 1 or 0
 --                                 sum up 1 or 0s
 --                          for all data