no diversity. needs static tests.

app/Main.hs

@@ -13,33 +13,32 @@ import System.IO
 
 import ArbitraryEnzymeTree
 import Environment
+import Evaluation
 
 -- Example definitions
 -- -------------------
 
 -- Enzymes
 
-pps :: Enzyme -- uses Phosphor from Substrate to produce PP
+-- pps :: Enzyme -- uses Phosphor from Substrate to produce PP
-pps = Enzyme "PPS" [(Substrate Phosphor,1)] ((Substrate Phosphor,-1),(Produced PP,1)) Nothing
+-- pps = Enzyme "PPS" [(Substrate Phosphor,1)] ((Substrate Phosphor,-1),(Produced PP,1)) Nothing
-
+--
-fpps :: Enzyme -- PP -> FPP
+-- fpps :: Enzyme -- PP -> FPP
-fpps = makeSimpleEnzyme (Produced PP) (Produced FPP)
+-- fpps = makeSimpleEnzyme (Produced PP) (Produced FPP)
 
 -- Environment
 
 exampleEnvironment :: Int -> [Enzyme] -> [(Predator,Probability)] -> [(Compound,Amount)] -> Environment
 exampleEnvironment addedC es pred tox =
   Environment
-    { soil = [ (Nitrate, 2)
+    { soil = [ (PPM, 10)
-             , (Phosphor, 3)
-             , (Photosynthesis, 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 = True
+    , settings = Settings { automimicry = False
                           , predatorsRandom = False
                           , numPlants = 150
                           }
@@ -47,29 +46,29 @@ exampleEnvironment addedC es pred tox =
 
 -- Plants
 
-examplePlants :: [Plant]
+-- examplePlants :: [Plant]
-examplePlants = (\g -> Plant g defaultAbsorption) <$> genomes
+-- examplePlants = (\g -> Plant g defaultAbsorption) <$> genomes
- where
+--  where
-   enzymes = [pps, fpps]
+--    enzymes = [pps, fpps]
-   quantity = [1,2] :: [Quantity]
+--    quantity = [1,2] :: [Quantity]
-   activation = [0.7, 0.9, 1]
+--    activation = [0.7, 0.9, 1]
-
+--
-   genomes = do
+--    genomes = do
-     e <- permutations enzymes
+--      e <- permutations enzymes
-     e' <- subsequences e
+--      e' <- subsequences e
-     q <- quantity
+--      q <- quantity
-     a <- activation
+--      a <- activation
-     return $ (,,) <$> e' <*> [q] <*> [a]
+--      return $ (,,) <$> e' <*> [q] <*> [a]
-
+--
-   defaultAbsorption = fmap ( limit Phosphor 2
+--    defaultAbsorption = fmap ( limit Phosphor 2
-                            . limit Nitrate 1
+--                             . limit Nitrate 1
-                            . limit Sulfur 0
+--                             . limit Sulfur 0
-                            ) <$> asks soil
+--                             ) <$> asks soil
-   -- custom absorbtion with helper-function:
+--    -- custom absorbtion with helper-function:
-   limit :: Nutrient -> Amount -> (Nutrient, Amount) -> (Nutrient, Amount)
+--    limit :: Nutrient -> Amount -> (Nutrient, Amount) -> (Nutrient, Amount)
-   limit n a (n', a')
+--    limit n a (n', a')
-     | n == n'   = (n, min a a') -- if we should limit, then we do ;)
+--      | n == n'   = (n, min a a') -- if we should limit, then we do ;)
-     | otherwise = (n', a')
+--      | otherwise = (n', a')
 
 -- Running the simulation
 -- ----------------------
@@ -88,21 +87,25 @@ loop loopAmount ps env = loop' loopAmount 0 ps env
     toxins :: [(Compound, Amount)]
     toxins = toxicCompounds env
     padded i str = take i $ str ++ repeat ' '
-    printEvery = 1
+    printEvery = 10
+    addedConstFitness = 0.1
     loop' :: Int -> Int -> [Plant] -> Environment -> IO ()
-    loop' loopAmount curLoop plants e = unless (loopAmount == curLoop) $ do
+    loop' loopAmount curLoop plants e = unless (loopAmount+1 == curLoop) $ do
       when (curLoop `mod` printEvery == 0) $ do
         putStr "\ESC[2J\ESC[H"
         printEnvironment e
         putStrLn ""
         putStrLn $ "Generation " ++ show curLoop ++ " of " ++ show loopAmount ++ ":"
       newPlants <- flip runReaderT e $ do
-          ! fs <- fmap (+0.01) <$> fitness plants -- fitness should be at least 0.01 for mating to work
+          (!fs,cs) <- unzip . fmap (\(f,c) -> (f,c)) <$> fitness plants
           let fps = zip plants fs -- gives us plants & their fitness in a tuple
               sumFitness = sum fs
-          when (curLoop `mod` printEvery == 0) $ do
+          when (curLoop `mod` printEvery == 0) $ liftIO $ do
-            liftIO $ printPopulation stringe fps
+            printPopulation stringe (zip3 plants fs cs)
-            liftIO $ hFlush stdout
+            putStrLn $ "Population statistics: VarC = " ++ (padded 50 . show . varianceOfProducedCompounds $ cs) 
+                                         ++ " DistC = " ++ (padded 50 . show . meanOfDistinctCompounds $ cs)
+            hFlush stdout
+            threadDelay $ 100*1000 -- sleep x*1000ns (=x ~ ms)
           -- generate x new plants.
           np <- asks (numPlants . settings)
           sequence . flip fmap [1..np] $ \_ -> do
@@ -117,16 +120,13 @@ loop loopAmount ps env = loop' loopAmount 0 ps env
                       | otherwise = findParent (x-f) ps
                     parent = findParent parent' fps
                 haploMate parent
-      hFlush stdout
-      when (curLoop `mod` printEvery == 0) $ do
-        threadDelay $ 100*1000 -- sleep x*1000ns (=x ~ ms)
       loop' loopAmount (curLoop+1) newPlants e
 
 main :: IO ()
 main = do
   hSetBuffering stdin NoBuffering
   --hSetBuffering stdout NoBuffering
-  randomCompounds <- makeHead (Substrate Photosynthesis) <$> generateTreeFromList 30 (toEnum <$> [(maxCompoundWithoutGeneric+1)..] :: [Compound]) -- generate roughly x compounds
+  randomCompounds <- makeHead (Substrate PPM) <$> generateTreeFromList 20 (toEnum <$> [(maxCompoundWithoutGeneric+1)..] :: [Compound]) -- generate roughly x compounds
   ds <- randoms <$> newStdGen
   probs <- randomRs (0.2,0.7) <$> newStdGen
   let poisonedTree     = poisonTree ds randomCompounds
@@ -134,14 +134,26 @@ main = do
   predators <- generatePredators 0.5 poisonedTree
   let env              = exampleEnvironment (getTreeSize randomCompounds) (generateEnzymeFromTree randomCompounds) (zip predators probs) poisonCompounds
       emptyPlants      = replicate (numPlants . settings $ env) emptyPlant
+  enzs <- randomRs (0,length (possibleEnzymes env) - 1) <$> newStdGen
+  let startPlants      = randomGenome 10 enzs (possibleEnzymes env) emptyPlants
   printEnvironment env
   writeFile "poison.twopi" $ generateDotFromPoisonTree "poison" 0.5 poisonedTree
   putStr "\ESC[?1049h"
-  loop 200 emptyPlants env
+  loop 2000 startPlants 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
@@ -152,7 +164,7 @@ generatePredators threshold t = do
           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.1,0.2)
+              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
@@ -168,27 +180,28 @@ printEnvironment (Environment soil pred metaIter maxComp toxic possEnz settings)
     putStrLn $ "Toxic:     " ++ show toxic
     putStrLn $ "Settings:  " ++ show settings
 
-printPopulation :: [(Enzyme,String)] -> [(Plant,Double)] -> IO ()
+printPopulation :: [(Enzyme,String)] -> [(Plant,Double,Vector Amount)] -> IO ()
 printPopulation es ps = do
   let padded i str = take i $ str ++ repeat ' '
   putStr $ padded 50 "Population:"
-  forM_ ps $ \(_,f) -> putStr (printColor f '█')
+  forM_ ps $ \(_,f,_) -> putStr (printColor f '█')
   putStrLn colorOff
   forM_ es $ \(e,s) -> do
     putStr s
-    forM_ ps $ \(Plant g _,_) -> do
+    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 > 2     = 'O'
-           | x > 1     = "+"
+           | x > 1     = '+'
-           | x > 0.7   = "ö"
+           | x > 0.7   = 'ö'
-           | x > 0.5   = "o"
+           | x > 0.5   = 'o'
-           | x > 0     = "."
+           | x > 0     = '.'
-           | otherwise = "_"
+           | otherwise = '_'
-      putStr (plot curE)
+          amount = min 2 $ cs ! fromEnum (fst . snd . synthesis $ e)
-    putStrLn ""
+      putStr $ printColor (amount/2) (plot curE)
+    putStrLn colorOff
 
 printColor :: Double -> Char -> String
 printColor x c

package.yaml

@@ -27,6 +27,7 @@ dependencies:
 - QuickCheck
 - pretty-simple
 - parallel
+- foldl
 
 library:
   source-dirs: src

src/Environment.hs

@@ -16,10 +16,7 @@ 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
+data Nutrient = PPM
-              | Phosphor
-              | Nitrate
-              | Photosynthesis
      deriving (Show, Enum, Bounded, Eq)
 
 -- | Fixed, non-generic Components
@@ -156,7 +153,7 @@ instance Eq Plant where
 
 type Fitness = Double
 
-fitness :: [Plant] -> World [Fitness]
+fitness :: [Plant] -> World [(Fitness, Vector Amount)]
 fitness ps = do
   nutrients <- mapM absorbNutrients ps           -- absorb soil
   products <- sequenceA $ zipWith produceCompounds ps nutrients -- produce compounds
@@ -172,18 +169,20 @@ fitness ps = do
   automimicry <- asks (automimicry . settings)
   popDefense <- if automimicry then
                    forM appearingPredators $ \p -> do
-                      as <- mapM (deterPredator p) products      -- how good can an individual deter p
+                      as <- mapM (dieToPredator p) products      -- how good can an individual deter p
                       return $ sum as / fromIntegral (length as) -- how good can the population deter p on average
                 else
                   return $ repeat 1
-  survivalRate <- mapM (deterPredators (zip appearingPredators popDefense)) products  -- defeat predators with produced compounds
+  dieRate <- mapM (dieToPredators (zip appearingPredators popDefense)) products  -- defeat predators with produced compounds
   let sumEnzymes = sum . fmap (\(_,q,a) -> fromIntegral q*a) . genome <$>  ps -- amount of enzymes * activation = resources "wasted"
       staticCostOfEnzymes = (\x -> 1 - 0.01*x) <$> sumEnzymes -- static cost of creating enzymes
   nutrientsAvailable <- fmap snd <$> asks soil
   let nutrientsLeft = (\p -> [p ! i | i <- [0..fromEnum (maxBound :: Nutrient)]]) <$> products
       nutrientRatio = minimum . zipWith (flip (/)) nutrientsAvailable <$> nutrientsLeft
-      costOfEnzymes = max 0 <$> zipWith (\s n -> s-n*0.1) staticCostOfEnzymes nutrientRatio -- cost to keep enzymes are static costs + amount of nutrient sucked out of the primary cycle
+      costOfEnzymes = max 0 <$> zipWith (\s n -> s-n*0.01) staticCostOfEnzymes nutrientRatio -- cost to keep enzymes are static costs + amount of nutrient sucked out of the primary cycle
-  return $ zipWith (*) survivalRate costOfEnzymes
+      survivalRate = (1-) <$> dieRate
+  return $ (,) <$> zipWith (*) survivalRate costOfEnzymes
+               <*> products
 
 produceCompounds :: Plant -> [(Nutrient, Amount)] -> World (Vector Amount)
 produceCompounds (Plant genes _) substrate = do
@@ -207,16 +206,17 @@ produceCompounds (Plant genes _) substrate = do
 -- so F(D) is omitted
 -- A(d_hat) is ahat * numAttacks p, because ahat is only deterrence of the population
 -- and does not incorporate the number of attacks, which A(d_hat) in the paper does.
-deterPredators :: [(Predator, Double)] -> Vector Amount -> World Probability
+dieToPredators :: [(Predator, Double)] -> Vector Amount -> World Probability
-deterPredators appearingPredators compounds = do
+dieToPredators [] _ = return 0 -- if no predator, no dying.
+dieToPredators appearingPredators compounds = do
   deters <- forM appearingPredators $ \(p,ahat) -> do
-    myDeter <- deterPredator p compounds
+    myDeter <- dieToPredator p compounds
-    return $ exp $ negate $ numAttacks p * ahat * myDeter -- exp due to assumption that number of attacks are poisson-distributed.
+    return $ ahat * myDeter -- exp due to assumption that number of attacks are poisson-distributed.
   return $ product deters
 
 
-deterPredator :: Predator -> Vector Amount -> World Double
+dieToPredator :: Predator -> Vector Amount -> World Double
-deterPredator p comps = do
+dieToPredator p comps = do
   toxins <- asks toxicCompounds
   return $ product [1 - min 1 (comps ! fromEnum t * l) | (t,l) <- toxins, t `elem` irresistance p]
 
@@ -227,12 +227,17 @@ deterPredator p comps = do
 -- | mate haploid
 haploMate :: Plant -> World Plant
 haploMate (Plant genes abs) = do
+  let digen :: IO [(Double, Int)]
+      digen = do
+        ds <- randoms <$> newStdGen
+        is <- randoms <$> newStdGen
+        return $ zip ds is
   --generate some random infinite uniform distributed lists of doubles in [0,1)
-  r1 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r1 <- liftIO digen
   r2 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
   r3 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
-  r4 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r4 <- liftIO digen
-  r5 <- liftIO ((randoms <$> newStdGen) :: IO [Double])
+  r5 <- liftIO digen
   enzymes <- asks possibleEnzymes
   re1 <- liftIO ((randomRs (0,length enzymes - 1) <$> newStdGen) :: IO [Int])
   re2 <- liftIO ((randomRs (0,length enzymes - 1) <$> newStdGen) :: IO [Int])
@@ -243,29 +248,36 @@ haploMate (Plant genes abs) = do
            . duplicateGene r4
            . deleteGene r5
            $ genes
-    deleteGene :: [Double] -> Genome -> Genome
+    deleteGene :: [(Double,Int)] -> Genome -> Genome
-    deleteGene (r:rs) ((e,1,a):gs) = if r < 0.1 then           deleteGene rs gs else (e,1,a):deleteGene rs gs
-    deleteGene (r:rs) ((e,q,a):gs) = if r < 0.1 then (e,q-1,a):deleteGene rs gs else (e,q,a):deleteGene rs gs
     deleteGene _ []                = []
+    deleteGene ((r,i):rs) g = if r < 0.05 then deleteGene rs (stay ++ go' ++ stay') else g
+                                          where
+                                            (stay, go:stay') = splitAt (i `mod` length g - 2) g
+                                            go' = case go of
+                                                    (e,1,a) -> []
+                                                    (e,q,a) -> [(e,q-1,a)]
 
-    duplicateGene :: [Double] -> Genome -> Genome
+    duplicateGene :: [(Double,Int)] -> Genome -> Genome
-    duplicateGene (r:rs) ((e,q,a):gs) = if r < 0.1 then (e,1,a):(e,q,a):duplicateGene rs gs else (e,q,a):duplicateGene rs gs
     duplicateGene _ []                = []
+    duplicateGene ((r,i):rs) g = if r < 0.05 then duplicateGene rs (stay ++ (e,q+1,a):stay') else g
+                                          where
+                                            (stay, (e,q,a):stay') = splitAt (i `mod` length g - 2) g
 
     addGene :: [Double] -> [Int] -> Genome -> Genome
-    addGene (r:rs) (s:ss) g = if r < 0.05 then (enzymes !! s,1,1):g else g
+    addGene (r:rs) (s:ss) g = if r < 0.005 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 :: [(Double,Int)] -> [Int] -> Genome -> Genome
-    mutateGene (r:rs) (s:ss) ((e,1,a):gs) = if r < 0.01 then ((enzymes !! s),1,a):mutateGene rs ss gs
+    mutateGene _ _ []                = []
-                                                        else (e,1,a):mutateGene rs ss gs
+    mutateGene ((r,i):rs) (s:ss) g = if r < 0.25 then mutateGene rs ss (stay ++ go' ++ stay') else g
-
+                                          where
-    mutateGene (r:rs) (s:ss) ((e,q,a):gs) = if r < 0.01 then (e,q-1,a):((enzymes !! s),1,a):mutateGene rs ss gs
+                                            (stay, go:stay') = splitAt (i `mod` length g - 2) g
-                                                        else (e,q,a):mutateGene rs ss gs
+                                            go' = case go of
-    mutateGene (r:rs) (s:ss) []           = []
+                                                    (e,1,a) -> [(enzymes !! s,1,a)]
+                                                    (e,q,a) -> [(e,q-1,a),(enzymes !! s,1,a)]
   return $ Plant genes' abs