scripts/haskell: add Feig.hs

haskell/Feig.hs

@@ -0,0 +1,99 @@
+#!/usr/bin/env nix-scrip
+#!>haskell
+#! haskell | Chart-gtk vector vector-algorithms
+
+import Control.Parallel.Strategies
+import Graphics.Rendering.Chart.Easy hiding (Vector)
+import Graphics.Rendering.Chart.Backend.Cairo
+import Data.List   (nubBy)
+import Data.Vector as V
+
+f :: Num a => a -> a -> a
+f a x = a * x*(1 - x)
+
+
+period :: (Ord a, Fractional a) => a -> a -> Vector a
+period start = nub . V.drop 8000 . sequence
+  where
+    sequence a = V.iterateN 8400 (f a) start
+    nub = fromList . nubBy (\x y -> abs(x-y)<1e-4). toList
+
+
+findFork :: Vector Int -> Vector Int
+findFork xs
+  | V.null xs = empty
+  | otherwise = find (V.head xs) indices xs
+  where
+    find :: Int -> Vector Int -> Vector Int -> Vector Int
+    find acc i x
+      | V.null x = V.empty
+      | (V.head x) == 2*acc
+        = V.cons (V.head i + 1) (find (2*acc) (V.tail i) (V.tail x))
+      | otherwise = find acc (V.tail i) (V.tail x)
+
+    indices = V.enumFromTo 0 (V.length xs - 1)
+
+--delta :: Vector Double
+delta = stuff
+  where
+    x     = V.enumFromThenTo 2.95 2.950001 3.54
+    y     = V.map (V.length . period 0.9) x `using` (evalTraversable rpar)
+    forks = V.map (x!) (findFork (V.uniq y))
+    stuff = V.map (y!) (findFork (V.uniq y))
+
+    diff x = V.zipWith (-) (V.tail x) x
+    rate x = V.zipWith (flip (/)) (V.tail x) x
+
+
+zipRepeat :: a -> Vector b -> Vector (a, b)
+zipRepeat x ys = V.zip (V.replicate (V.length ys) x) ys
+
+
+biforc :: Renderable ()
+biforc = toRenderable layout where
+  layout =
+      layout_title       .~ "Logistic map bifurcation"
+    $ layout_title_style .~ titleOpts
+    $ layout_background  .~ FillStyleSolid (opaque black)
+
+    $ layout_x_axis . laxis_title .~ "bifurcation parameter"
+    $ layout_x_axis . laxis_title_style .~ titleOpts
+
+    $ layout_y_axis . laxis_title .~ "fixed points"
+    $ layout_y_axis . laxis_title_style .~ titleOpts
+
+    $ layout_x_axis . laxis_style . axis_label_style .~ labelOpts
+    $ layout_y_axis . laxis_style . axis_label_style .~ labelOpts
+
+    $ layout_x_axis . laxis_generate .~ scaledAxis def (3.5,4.0)
+    $ layout_x_axis . laxis_override .~ axisGridHide
+    $ layout_y_axis . laxis_override .~ axisGridHide
+    $ layout_plots .~ [toPlot points]
+    $ def
+
+  points =
+      plot_points_style  .~ filledCircles 1 (opaque orange)
+    $ plot_points_values .~ (V.toList xy)
+    $ plot_points_title  .~ "starting 1/3"
+    $ def
+
+  titleOpts = def { _font_size  = 146
+                  , _font_color = opaque white
+                  }
+  labelOpts = def { _font_size  = 96
+                  , _font_color = opaque white
+                  }
+
+  x  = V.enumFromThenTo 3.4 3.40008 4.0 :: Vector Double
+  y  = V.map (period 0.9) x
+  xy = (V.zip x y >>= uncurry zipRepeat) `using` (evalTraversable rdeepseq)
+
+
+main :: IO ()
+main = do
+  print delta
+  --renderableToFile opts "feig.png" biforc
+  where
+    opts = fo_size .~ (14564, 8192) 
+         $ fo_format .~ PNG
+         $ def