diff --git a/app/Main.hs b/app/Main.hs index 762ba07..07a4e4b 100644 --- a/app/Main.hs +++ b/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 "PPS" [(Substrate Phosphor,1)] ((Substrate Phosphor,-1),(Produced PP,1)) Nothing - -fpps :: Enzyme -- PP -> FPP -fpps = makeSimpleEnzyme (Produced PP) (Produced FPP) +-- 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 = [ (Nitrate, 2) - , (Phosphor, 3) - , (Photosynthesis, 10) + { 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 = True + , settings = Settings { automimicry = False , predatorsRandom = False , numPlants = 150 } @@ -47,29 +46,29 @@ exampleEnvironment addedC es pred tox = -- 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 - ) <$> asks 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') +-- 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 +-- ) <$> asks 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 -- ---------------------- @@ -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 - liftIO $ printPopulation stringe fps - liftIO $ hFlush stdout + when (curLoop `mod` printEvery == 0) $ liftIO $ do + printPopulation stringe (zip3 plants fs cs) + 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 > 1 = "+" - | x > 0.7 = "ö" - | x > 0.5 = "o" - | x > 0 = "." - | otherwise = "_" - putStr (plot curE) - putStrLn "" + | 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 diff --git a/package.yaml b/package.yaml index d0ee051..6bf64cc 100644 --- a/package.yaml +++ b/package.yaml @@ -27,6 +27,7 @@ dependencies: - QuickCheck - pretty-simple - parallel +- foldl library: source-dirs: src diff --git a/src/Environment.hs b/src/Environment.hs index b505bd4..dceb9e1 100644 --- a/src/Environment.hs +++ b/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 - | Phosphor - | Nitrate - | Photosynthesis +data Nutrient = PPM 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 - return $ zipWith (*) survivalRate costOfEnzymes + 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 + 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 -deterPredators appearingPredators compounds = do +dieToPredators :: [(Predator, Double)] -> Vector Amount -> World Probability +dieToPredators [] _ = return 0 -- if no predator, no dying. +dieToPredators appearingPredators compounds = do deters <- forM appearingPredators $ \(p,ahat) -> do - myDeter <- deterPredator p compounds - return $ exp $ negate $ numAttacks p * ahat * myDeter -- exp due to assumption that number of attacks are poisson-distributed. + myDeter <- dieToPredator p compounds + return $ ahat * myDeter -- exp due to assumption that number of attacks are poisson-distributed. return $ product deters -deterPredator :: Predator -> Vector Amount -> World Double -deterPredator p comps = do +dieToPredator :: Predator -> Vector Amount -> World Double +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]) - r5 <- liftIO ((randoms <$> newStdGen) :: IO [Double]) + r4 <- liftIO digen + 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 (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 :: [(Double,Int)] -> Genome -> Genome 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 (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 :: [(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 + 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 (r:rs) (s:ss) ((e,1,a):gs) = if r < 0.01 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.01 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) [] = [] + mutateGene :: [(Double,Int)] -> [Int] -> Genome -> Genome + mutateGene _ _ [] = [] + mutateGene ((r,i):rs) (s:ss) g = if r < 0.25 then mutateGene rs ss (stay ++ go' ++ stay') else g + where + (stay, go:stay') = splitAt (i `mod` length g - 2) g + go' = case go of + (e,1,a) -> [(enzymes !! s,1,a)] + (e,q,a) -> [(e,q-1,a),(enzymes !! s,1,a)] return $ Plant genes' abs