diff options
| -rw-r--r-- | 01.06.hs | 90 | ||||
| -rw-r--r-- | 01.07.hs | 168 | ||||
| -rw-r--r-- | 02.05.hs | 38 | ||||
| -rw-r--r-- | 08.01.hs | 46 |
4 files changed, 205 insertions, 137 deletions
@@ -1,70 +1,96 @@ - {-# LANGUAGE FlexibleInstances #-} -import System.Random -import Data.Random.Normal import Control.Monad.State +import Data.Random.Normal import Numeric.LinearAlgebra as Lin +import System.Random -import Control.Monad(void) +import Control.Monad (void) import Graphics.Rendering.Chart +import Graphics.Rendering.Chart.Backend.Cairo import Graphics.Rendering.Chart.Easy import Graphics.Rendering.Chart.Grid -import Graphics.Rendering.Chart.Backend.Cairo -import Graphics.Rendering.Chart.Plot.Histogram(defaultNormedPlotHist) +import Graphics.Rendering.Chart.Plot.Histogram (defaultNormedPlotHist) + instance Default (PlotHist x Double) where - def = defaultNormedPlotHist + def = defaultNormedPlotHist -- | Evaluates the Winger surmise. ρWinger :: Floating t => t -> t -ρWinger s = pi * s / 2 * exp(- pi * s^2 / 4) +ρWinger s = pi * s / 2 * exp (-pi * s ^ 2 / 4) -- | Computes the splitting of the center eigenvalues of a symmetric matrix. eigenSplitting :: Herm Double -> Double eigenSplitting m = e2 - e1 - where - e = sortVector $ eigenvaluesSH m - i = (Lin.size e) `quot` 2 - e1 = e Lin.! (i - 1) - e2 = e Lin.! i + where + e = sortVector $ eigenvaluesSH m + i = (Lin.size e) `quot` 2 + e1 = e Lin.! (i - 1) + e2 = e Lin.! i -- | Produces random GOE matrices. goeGen :: RandomGen g => Int -> State g (Herm Double) -goeGen n = get >>= (\g -> let (g1, g2) = split g in - (put g1) >> (return $ sym $ matrix n $ take (n^2) $ normals' (0.0, 2.0) g2)) +goeGen n = + get >>= + (\g -> + let (g1, g2) = split g + in (put g1) >> (return $ sym $ matrix n $ take (n ^ 2) $ normals' (0.0, 2.0) g2)) -- | Produces random ∓1 matrices. pm1Gen :: RandomGen g => Int -> State g (Herm Double) -pm1Gen n = state (\g -> let (m, g2) = runState (goeGen n) g in (trustSym $ cmap signum (unSym m), g2)) +pm1Gen n = + state + (\g -> + let (m, g2) = runState (goeGen n) g + in (trustSym $ cmap signum (unSym m), g2)) -- | Produces n eigenvalue splittings from the matrix ensemble m. eigenSplittings :: RandomGen g => Int -> State g (Herm Double) -> State g [Double] -eigenSplittings n m = state (\g -> let (mats, g2) = runState (sequence $ replicate n m) g in (eigenSplitting <$> mats, g2)) +eigenSplittings n m = + state + (\g -> + let (mats, g2) = runState (sequence $ replicate n m) g + in (eigenSplitting <$> mats, g2)) hsChart :: RandomGen g => g -> State g (Herm Double) -> Int -> Int -> String -> Layout Double Double -hsChart g gen num bins tit = execEC $ do - plot $ fmap histToPlot $ liftEC $ do - plot_hist_title .= tit - plot_hist_bins .= bins - plot_hist_values .= scaledData - plot_hist_range .= Just (0.0, 5.0) - plot_hist_norm_func .= const (\x -> fromIntegral bins * fromIntegral x / fromIntegral num / 5.0) - plot (line "ρ_Winger" [ [(s,ρWinger s) | s <- [0,0.005..5]] ]) +hsChart g gen num bins tit = + execEC $ do + plot $ + fmap histToPlot $ + liftEC $ do + plot_hist_title .= tit + plot_hist_bins .= bins + plot_hist_values .= scaledData + plot_hist_range .= Just (0.0, 5.0) + plot_hist_norm_func .= + const (\x -> fromIntegral bins * fromIntegral x / fromIntegral num / 5.0) + plot (line "ρ_Winger" [[(s, ρWinger s) | s <- [0,0.005 .. 5]]]) where rawData = evalState (eigenSplittings num gen) g mean = sum rawData / fromIntegral (length rawData) scaledData = (\x -> x / mean) <$> rawData -mkgrid g num bins = title `wideAbove` aboveN [ besideN [ layoutToGrid (hsChart g ((fst gen) n) num bins (snd gen Prelude.<> show n)) | gen <- gens] | n <- ns] +mkgrid g num bins = + title `wideAbove` + aboveN + [ besideN + [layoutToGrid (hsChart g ((fst gen) n) num bins (snd gen Prelude.<> show n)) | gen <- gens] + | n <- ns + ] where gens = [(goeGen, "GOE, N = "), (pm1Gen, "±1, N = ")] ns = [2, 4, 10] - title = setPickFn nullPickFn $ label ls HTA_Centre VTA_Centre "Winger surmise compared with eigenvalue splitting of random matrices" - ls = def { _font_size = 15, _font_weight = FontWeightBold } - + title = + setPickFn nullPickFn $ + label + ls + HTA_Centre + VTA_Centre + "Winger surmise compared with eigenvalue splitting of random matrices" + ls = def {_font_size = 15, _font_weight = FontWeightBold} main = do g <- newStdGen - void $ renderableToFile (def { _fo_format = PDF}) "01.06.pdf" $ fillBackground def $ gridToRenderable $ mkgrid g 1000 50 - - + void $ + renderableToFile (def {_fo_format = PDF}) "01.06.pdf" $ + fillBackground def $ gridToRenderable $ mkgrid g 1000 50 @@ -1,27 +1,27 @@ - {-# LANGUAGE NoMonomorphismRestriction #-} -import System.Random import Control.Monad.State -import Data.Maybe import Data.List +import Data.Maybe import Data.Tuple +import System.Random -import Graphics.Rendering.Chart.Easy import Graphics.Rendering.Chart.Backend.Cairo +import Graphics.Rendering.Chart.Easy import Data.Colour.Names -import Diagrams.Prelude import Diagrams.Backend.SVG.CmdLine - +import Diagrams.Prelude ----------------------------------------- -- Network construction and properties -- ----------------------------------------- - type Node = Int + type Edge = (Node, Node) + type Adjacency = [Node] + type Network = [Adjacency] hasNode :: Node -> Network -> Bool @@ -37,7 +37,7 @@ addEdge :: Edge -> Network -> Network addEdge (v1, v2) n = addHalfedge v1 v2 (addHalfedge v2 v1 n) getNodes :: Network -> [Node] -getNodes n = [0..(length n - 1)] +getNodes n = [0 .. (length n - 1)] getEdges :: Network -> [Edge] getEdges n = foldr (<>) [] $ map (\i -> map (\j -> (fst i, j)) (snd i)) (zip (getNodes n) n) @@ -53,136 +53,152 @@ emptyNetwork l = addNode (emptyNetwork (l - 1)) addShortEdges :: Int -> Network -> Network addShortEdges 0 n = n addShortEdges 1 n = n -addShortEdges z n = let - l = length n - addZEdge v = addEdge (v, mod (v + quot z 2) l) - in addShortEdges (z - 2) (foldr addZEdge n [0..(l - 1)]) +addShortEdges z n = + let l = length n + addZEdge v = addEdge (v, mod (v + quot z 2) l) + in addShortEdges (z - 2) (foldr addZEdge n [0 .. (l - 1)]) -- | Lists all possible long edges. -possibleLongEdges l z = let - d = quot z 2 + 1 - e v = (\u -> (v, mod (v + u) l)) <$> [d..(l - d)] - in foldr (<>) [] $ e <$> [0..(l - 1)] +possibleLongEdges l z = + let d = quot z 2 + 1 + e v = (\u -> (v, mod (v + u) l)) <$> [d .. (l - d)] + in foldr (<>) [] $ e <$> [0 .. (l - 1)] -- | Take a random element from a list, removing it from the list. randomChoice :: RandomGen g => State ([a], g) (Maybe a) -randomChoice = get >>= (\(xs, g) -> if (length xs == 0) then (return Nothing) else (let (r, g2) = randomR (0, length xs - 1) g in (put ((take r xs) ++ (drop (r+1) xs), g2)) >> return (Just (xs !! r)))) +randomChoice = + get >>= + (\(xs, g) -> + if (length xs == 0) + then (return Nothing) + else (let (r, g2) = randomR (0, length xs - 1) g + in (put ((take r xs) ++ (drop (r + 1) xs), g2)) >> return (Just (xs !! r)))) -- | Returns a list of randomly chosen values with no repeats. randomChoices :: RandomGen g => Int -> State ([a], g) [a] -randomChoices n = state (\si -> - let (maybeVals, sf) = runState (sequence $ replicate n randomChoice) si - in (catMaybes maybeVals, sf) - ) +randomChoices n = + state + (\si -> + let (maybeVals, sf) = runState (sequence $ replicate n randomChoice) si + in (catMaybes maybeVals, sf)) -- | Adds the appropriate number of long edges to a network. addLongEdges :: (RealFrac a, RandomGen g) => Int -> Int -> a -> Network -> State g Network -addLongEdges l z p n = state (\g -> - let - nle = floor (p * fromIntegral (l * quot z 2)) - (es, (_, g2)) = runState (randomChoices nle) (possibleLongEdges l z, g) - in (foldr addEdge n $ es, g2) - ) +addLongEdges l z p n = + state + (\g -> + let nle = floor (p * fromIntegral (l * quot z 2)) + (es, (_, g2)) = runState (randomChoices nle) (possibleLongEdges l z, g) + in (foldr addEdge n $ es, g2)) constructNetwork :: (RealFrac a, RandomGen g) => Int -> Int -> a -> State g Network constructNetwork l z p = addLongEdges l z p $ addShortEdges z $ emptyNetwork l - ---------------------------- -- Finding shortest paths -- ---------------------------- - pathLengths :: Network -> (Int, [Node], [Maybe Int]) -> (Int, [Node], [Maybe Int]) pathLengths n (level, curNodes, knownDists) = - let - setVertexFound node dists = - if isNothing (head after) - then (before ++ [Just level] ++ (drop 1 after)) - else dists - where (before, after) = splitAt node dists - newDists = foldr setVertexFound knownDists curNodes - in (level + 1, nub $ filter (\x -> isNothing (newDists !! x)) $ foldr (<>) [] $ getNeighbors n <$> curNodes, newDists) + let setVertexFound node dists = + if isNothing (head after) + then (before ++ [Just level] ++ (drop 1 after)) + else dists + where + (before, after) = splitAt node dists + newDists = foldr setVertexFound knownDists curNodes + in ( level + 1 + , nub $ filter (\x -> isNothing (newDists !! x)) $ foldr (<>) [] $ getNeighbors n <$> curNodes + , newDists) findPathLengthsFromNode :: Network -> Node -> [Int] -findPathLengthsFromNode n v = catMaybes $ (\(_,_,d) -> d) $ head $ dropWhile (\(_,nodes,_) -> length nodes > 0) $ iterate (pathLengths n) (0, [v], replicate (length n) Nothing) +findPathLengthsFromNode n v = + catMaybes $ + (\(_, _, d) -> d) $ + head $ + dropWhile (\(_, nodes, _) -> length nodes > 0) $ + iterate (pathLengths n) (0, [v], replicate (length n) Nothing) findAllPathLengths :: Network -> [Int] findAllPathLengths n = foldr (<>) [] $ (findPathLengthsFromNode n) <$> getNodes n findAveragePathLength :: RealFrac a => Network -> a -findAveragePathLength n = let xs = findAllPathLengths n in fromIntegral (sum xs) / fromIntegral (length xs) - -shortestPath f n i j 0 | i == j = 0 - | elem (i, j) $ getEdges n = 1 - | elem (j, i) $ getEdges n = 1 - | otherwise = length n +findAveragePathLength n = + let xs = findAllPathLengths n + in fromIntegral (sum xs) / fromIntegral (length xs) +shortestPath f n i j 0 + | i == j = 0 + | elem (i, j) $ getEdges n = 1 + | elem (j, i) $ getEdges n = 1 + | otherwise = length n shortestPath f n i j k = min (f n i j (k - 1)) (f n i k (k - 1) + f n k j (k - 1)) -data NaturalTree a = NNode a (NaturalTree a) (NaturalTree a) +data NaturalTree a = + NNode a (NaturalTree a) (NaturalTree a) -NNode a tl tr !!! 0 = a +NNode a tl tr !!! 0 = a NNode a tl tr !!! n = - if odd n - then tl !!! top - else tr !!! (top-1) - where top = n `div` 2 + if odd n + then tl !!! top + else tr !!! (top - 1) + where + top = n `div` 2 instance Functor NaturalTree where - fmap f (NNode a tl tr) = NNode (f a) (fmap f tl) (fmap f tr) + fmap f (NNode a tl tr) = NNode (f a) (fmap f tl) (fmap f tr) -naturals r n = - NNode n - ((naturals $! r2) $! (n+r)) - ((naturals $! r2) $! (n+r2)) - where r2 = 2*r +naturals r n = NNode n ((naturals $! r2) $! (n + r)) ((naturals $! r2) $! (n + r2)) + where + r2 = 2 * r shortestPaths n = - let - nodes = getNodes n - memo = fmap (\y -> fmap (\z -> fmap z (naturals 1 0)) y) (fmap (\x -> fmap x (naturals 1 0)) (fmap (shortestPath shortestPath' n) (naturals 1 0))) - shortestPath' n a b c = memo !!! a !!! b !!! c - in - (\i j -> shortestPath' n i j (length n - 1)) <$> nodes <*> nodes - + let nodes = getNodes n + memo = + fmap + (\y -> fmap (\z -> fmap z (naturals 1 0)) y) + (fmap (\x -> fmap x (naturals 1 0)) (fmap (shortestPath shortestPath' n) (naturals 1 0))) + shortestPath' n a b c = memo !!! a !!! b !!! c + in (\i j -> shortestPath' n i j (length n - 1)) <$> nodes <*> nodes ---------------------------- -- Betweenness centrality -- ---------------------------- - ---------------------- -- Drawing routines -- ---------------------- - drawNode :: Node -> Diagram B drawNode n = named n $ fc black $ circle 0.1 drawEdge :: Edge -> (Diagram B -> Diagram B) -drawEdge (n1, n2) = withName n1 $ \b1 -> withName n2 $ \b2 -> Diagrams.Prelude.atop ((location b1 ~~ location b2)) +drawEdge (n1, n2) = + withName n1 $ \b1 -> withName n2 $ \b2 -> Diagrams.Prelude.atop ((location b1 ~~ location b2)) drawNetwork :: Network -> Diagram B -drawNetwork n = atPoints (trailVertices $ regPoly (length n) 1) (map drawNode $ getNodes n) Diagrams.Prelude.# applyAll (map drawEdge $ getEdges n) - +drawNetwork n = + atPoints (trailVertices $ regPoly (length n) 1) (map drawNode $ getNodes n) Diagrams.Prelude.# + applyAll (map drawEdge $ getEdges n) ----------------------- -- Plotting routines -- ----------------------- - lengthHistogram :: Colour Double -> [Int] -> EC (Layout Double Int) () -lengthHistogram color lengths = plot $ fmap histToPlot $ liftEC $ do +lengthHistogram color lengths = + plot $ + fmap histToPlot $ + liftEC $ do plot_hist_fill_style .= def {_fill_color = (withOpacity color 0.5)} plot_hist_line_style .= def {_line_color = (withOpacity color 1.0)} plot_hist_bins .= 20 plot_hist_values .= ((fromIntegral <$> lengths) :: [Double]) plot_hist_norm_func .= const id - main = do g1 <- newStdGen g2 <- newStdGen - void $ renderableToFile (def { _fo_format = SVG}) "01.07.hist.svg" $ fillBackground def $ toRenderable $ do - lengthHistogram blue $ shortestPaths $ evalState (constructNetwork 200 2 0.02) g1 - lengthHistogram red $ shortestPaths $ evalState (constructNetwork 200 2 0.20) g2 - + void $ + renderableToFile (def {_fo_format = SVG}) "01.07.hist.svg" $ + fillBackground def $ + toRenderable $ do + lengthHistogram blue $ shortestPaths $ evalState (constructNetwork 200 2 0.02) g1 + lengthHistogram red $ shortestPaths $ evalState (constructNetwork 200 2 0.20) g2 -- mainWith (drawNetwork n) - @@ -1,17 +1,24 @@ - -import System.Random import Control.Monad.State import Data.List.Split (chunksOf) +import System.Random -import Graphics.Rendering.Chart.Easy import Graphics.Rendering.Chart.Backend.Cairo +import Graphics.Rendering.Chart.Easy import Data.Colour.Names -import Diagrams.Prelude import Diagrams.Backend.Postscript +import Diagrams.Prelude randomWalk :: (RandomGen g) => Int -> Int -> State g [[Double]] -randomWalk d n = get >>= (\g -> let (g1, g2) = split g in (put g1) >> (return $ take n $ scanl (\x y -> map (uncurry (+)) $ zip x y) (replicate d 0.0) $ chunksOf d $ randomRs (-0.5,0.5) g2)) +randomWalk d n = + get >>= + (\g -> + let (g1, g2) = split g + in (put g1) >> + (return $ + take n $ + scanl (\x y -> map (uncurry (+)) $ zip x y) (replicate d 0.0) $ + chunksOf d $ randomRs (-0.5, 0.5) g2)) drawWalk :: [[Double]] -> Diagram B drawWalk w = fromVertices $ (\[x, y] -> p2 (x, y)) <$> w @@ -21,10 +28,21 @@ main = do let (w10, g2) = runState (randomWalk 2 10) g let (w1000, g3) = runState (randomWalk 2 1000) g2 let (w100000, g4) = runState (randomWalk 2 100000) g3 - renderDia Postscript (PostscriptOptions "02.05.01.eps" (dims2D 500 500) EPS) $ foldr1 Diagrams.Prelude.atop [lc red $ drawWalk w10, lc blue $ drawWalk w1000, lc green $ drawWalk w100000] + renderDia Postscript (PostscriptOptions "02.05.01.eps" (dims2D 500 500) EPS) $ + foldr1 + Diagrams.Prelude.atop + [lc red $ drawWalk w10, lc blue $ drawWalk w1000, lc green $ drawWalk w100000] let (e1, g5) = runState (sequence $ replicate 10000 $ randomWalk 2 2) g4 let (e10, g6) = runState (sequence $ replicate 10000 $ randomWalk 2 11) g5 - renderDia Postscript (PostscriptOptions "02.05.02.eps" (dims2D 500 500) EPS) $ foldr1 Diagrams.Prelude.atop [mconcat $ map (place (lw Diagrams.Prelude.none $ fc yellow $ ((circle 0.01) :: Diagram B))) (map (\[x, y] -> p2 (x, y)) (map last e1)), mconcat $ map (place (lw Diagrams.Prelude.none $ fc black $ circle 0.01)) (map (\[x, y] -> p2 (x, y)) (map last e10))] - - - + renderDia Postscript (PostscriptOptions "02.05.02.eps" (dims2D 500 500) EPS) $ + foldr1 + Diagrams.Prelude.atop + [ mconcat $ + map + (place (lw Diagrams.Prelude.none $ fc yellow $ ((circle 0.01) :: Diagram B))) + (map (\[x, y] -> p2 (x, y)) (map last e1)) + , mconcat $ + map + (place (lw Diagrams.Prelude.none $ fc black $ circle 0.01)) + (map (\[x, y] -> p2 (x, y)) (map last e10)) + ] @@ -1,26 +1,33 @@ - -import System.Random import Control.Monad.State +import Data.Foldable import Data.Maybe -import Data.Tuple import Data.Sequence as Seq -import Data.Foldable +import Data.Tuple +import System.Random -- | Take a random element from a list, removing it from the list. randomChoice :: RandomGen g => State (Seq a, g) (Maybe a) -randomChoice = get >>= (\(xs, g) -> if (Seq.null xs) then (return Nothing) else (let - (r, g2) = randomR (0, Seq.length xs - 1) g - in (put (deleteAt r xs, g2) >> (return (Just $ index xs r)) - ))) +randomChoice = + get >>= + (\(xs, g) -> + if (Seq.null xs) + then (return Nothing) + else (let (r, g2) = randomR (0, Seq.length xs - 1) g + in (put (deleteAt r xs, g2) >> (return (Just $ index xs r))))) -- | Returns a list of randomly chosen values with no repeats. randomChoices :: RandomGen g => Int -> State (Seq a, g) (Seq a) -randomChoices n = state (\si -> - let (maybeVals, sf) = runState (sequence $ Prelude.replicate n randomChoice) si - in (fromList $ catMaybes maybeVals, sf) - ) +randomChoices n = + state + (\si -> + let (maybeVals, sf) = runState (sequence $ Prelude.replicate n randomChoice) si + in (fromList $ catMaybes maybeVals, sf)) + +data Bacterium + = Red + | Green + deriving (Show) -data Bacterium = Red | Green deriving Show type Bacteria = Seq Bacterium reproduce :: Bacteria -> Bacteria @@ -36,21 +43,22 @@ initialize :: Int -> Bacteria initialize n = (Seq.replicate (n `quot` 2) Red) <> (Seq.replicate (n `quot` 2) Green) isRed :: Bacterium -> Bool -isRed Red = True +isRed Red = True isRed Green = False nRed :: Bacteria -> Int nRed b = Seq.length $ Seq.filter isRed b isDone :: Bacteria -> Bool -isDone b | nRed b == 0 = True - | nRed b == Seq.length b = True - | otherwise = False +isDone b + | nRed b == 0 = True + | nRed b == Seq.length b = True + | otherwise = False run :: RandomGen g => Int -> g -> ([Bacteria], g) run n g = (fst <$> steps, snd $ last steps) - where steps = takeWhile (not . isDone . fst) $ iterate hour (initialize n, g) + where + steps = takeWhile (not . isDone . fst) $ iterate hour (initialize n, g) lifetime :: [Bacteria] -> Int lifetime b = Data.Foldable.length b - |