100 lines
2.8 KiB
Haskell
Executable File
100 lines
2.8 KiB
Haskell
Executable File
#!/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
|