n-d isi histogram

nd-hist.hs

{-# LANGUAGE TupleSections #-}

import             Control.Applicative
import             Control.Arrow
--import           Data.Monoid
import             Data.Traversable    as T
import             Data.Foldable       as F
import             Data.List           as L
import             Data.Function
--import           Data.Tuple
import             Data.Ratio
import   qualified Data.Tuple.Extra    as U                  -- http://hackage.haskell.org/package/extra

import   qualified Data.Vector.Unboxed as V
--import qualified Data.ListLike       as LL                 -- http://hackage.haskell.org/package/ListLike

import   qualified Data.HashMap.Lazy   as M                  -- http://hackage.haskell.org/package/unordered-containers
import   qualified Data.IntMap.Lazy    as I

import             Data.Function.Pointless                   -- http://hackage.haskell.org/package/pointless-fun

import             Statistics.Distribution                   -- http://hackage.haskell.org/package/statistics
import             Statistics.Distribution.Gamma
import             System.Random.MWC
import             Control.Monad.Primitive

import             Data.Random                               -- https://hackage.haskell.org/package/random-fu

import             Graphics.Rendering.Chart                  -- https://hackage.haskell.org/package/Chart
--import           Graphics.Rendering.Chart.Easy
import             Graphics.Rendering.Chart.Backend.Diagrams -- https://hackage.haskell.org/package/Chart-diagrams

--import           Data.Colour.Names
--import           Data.Colour
import             Data.Colour.RGBSpace
import             Data.Colour.RGBSpace.HSV
import             Data.Colour.SRGB

import             Data.Default.Class
import             Control.Lens

samples :: (ContGen d, PrimMonad m) => Int -> d -> Gen (PrimState m) -> m (V.Vector Double)
samples n = V.replicateM n .: genContVar
-- how generalize to include genDiscreteVar?
-- any advantage to genContinous (sic)?
-- how generalize to any sequential type?

sysSamples :: (ContGen d) => Int -> d -> IO (V.Vector Double)
sysSamples = withSystemRandom . asGenST .: samples

--sysSamples' :: (MonadRandom m, V.Unbox a) => Int -> RVar a      -> m (V.Vector a     )
sysSamples'   :: (MonadRandom m           ) => Int -> RVar Double -> m (V.Vector Double)
sysSamples' = V.replicateM .^ sample

chart f = toRenderable layout
  where
    layout = layout_plots .~ [ toPlot spots ]
           $ def

    spots = area_spots_4d_max_radius .~ 2
          $ area_spots_4d_palette    .~ palette
          $ area_spots_4d_values     .~ values
          $ area_spots_4d_opacity    .~ 0.75
          $ def

    palette = uncurryRGB (rgbUsingSpace sRGBSpace) <$> (\h -> hsv h 1 1) <$> shift (linspace (360,1) $ 2^9) 90
    values = (\([x,y],v) -> (x, y, 2::Int, v)) <$> f

shift xs = take (length xs) . flip drop (cycle xs) . (`mod` length xs)

linspace (a,b) n = (+ a) <$> (* ((b - a) / n)) <$> [0 .. n]

-- need to generalize from []...
multiHist :: (RealFrac b, Ord b) => [b] -> Int -> [b] -> [([b], Int)]
multiHist edges' n = (first ((centers I.!) <$>) <$>) . sortBy (compare `on` snd) . count . (bin <$>) . slide
  where slide = getZipList . sequenceA . (ZipList <$>) . take n . tails -- http://stackoverflow.com/a/27733778/1441998
        bin   = (length . flip takeWhile edges . (>) <$>)
        count = M.toList . F.foldr (flip (M.insertWith (+)) 1) M.empty

        -- centers = (edges !!!) <$> [0 .. ledges] -- worth the hassle of an IntMap for faster lookup?
        centers    = I.fromDistinctAscList $ ((id &&& (edges !!!)) <$>) [0 .. ledges]
          where a !!! b | not $ F.elem b [0, ledges] = (F.sum $ (a !!) . (b +) <$> [-1, 0]) / 2
                        | otherwise                  = 2*(a !! (b - x)) - (a !!! (b - 2*x + 1))
                        where x = round $ b % ledges

        edges  = sort edges'
        ledges = length edges

main = do
  isis    <- sysSamples  spikes $ gammaDistr shape scale
  -- isis <- sysSamples' spikes $ gamma      shape scale -- bug on windows https://github.com/mokus0/random-fu/issues/32

  let edges = linspace ((V.minimum &&& V.maximum) isis) nedges -- for single pass: https://hackage.haskell.org/package/foldl-1.1.1/docs/Control-Foldl.html

  renderableToFile fileopts "2d isi.svg" . chart $ multiHist edges 2 $ V.toList isis

  where spikes = 1000000
        shape  = 9.0
        scale  = 0.5
        nedges = 100
        fileopts = fo_format .~ SVG 
                 $ fo_size   .~ U.both (* 4) (U.dupe nedges)
                 $ def