#!/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