{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
module Plots.Types.HeatMap
(
HeatMap
, heatMap
, heatMap'
, heatMapIndexed
, heatMapIndexed'
, HasHeatMap (..)
, pathHeatRender
, pixelHeatRender
, pixelHeatRender'
, HeatMatrix (..)
, heatImage
, hmPoints
, mkHeatMap
, mkHeatSurface
, mkHeatMatrix
, mkHeatMatrix'
) where
import Control.Lens hiding (transform, ( # ))
import Control.Monad.ST
import Control.Monad.State
import qualified Data.Foldable as F
import Data.Typeable
import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Unboxed as V
import Data.Word (Word8)
import Codec.Picture
import Diagrams.Coordinates.Isomorphic
import Diagrams.Prelude
import Plots.Axis
import Plots.Style
import Plots.Types
import Diagrams.TwoD.Image
import Geometry.ThreeD.Shapes
data HeatMatrix = HeatMatrix
{ hmSize :: {-# UNPACK #-} !(V2 Int)
, _hmVector :: {-# UNPACK #-} !(V.Vector Double)
, hmBoundLower :: {-# UNPACK #-} !Double
, hmBoundUpper :: {-# UNPACK #-} !Double
}
mkHeatMatrix :: V2 Int -> (V2 Int -> Double) -> HeatMatrix
mkHeatMatrix s@(V2 x y) f = runST $ do
mv <- M.new (x*y)
let go !q !a !b !i !j
| j == y = do v <- V.unsafeFreeze mv
return (HeatMatrix s v a b)
| i == x = go q a b 0 (j + 1)
| otherwise = do let !d = f (V2 i j)
M.unsafeWrite mv q d
go (q + 1) (min a d) (max b d) (i + 1) j
go 0 (1/0) (-1/0) 0 0
{-# INLINE mkHeatMatrix #-}
min' :: Double -> Double -> Double
min' !a !b
| isNaN b = a
| isInfinite b = a
| otherwise = min a b
{-# INLINE min' #-}
max' :: Double -> Double -> Double
max' !a !b
| isNaN b = a
| isInfinite b = a
| otherwise = max a b
{-# INLINE max' #-}
data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double
minMax :: V.Vector Double -> (Double, Double)
minMax = fini . V.foldl' go (MM (1/0) (-1/0))
where
go (MM lo hi) k = MM (min' lo k) (max' hi k)
fini (MM lo hi) = (lo, hi)
{-# INLINE minMax #-}
mkHeatMatrix' :: (F.Foldable f, F.Foldable g) => f (g Double) -> HeatMatrix
mkHeatMatrix' xss = HeatMatrix (V2 x y) vd a b
where
(a,b) = minMax vd
vd = V.create $ do
mv <- M.new (x*y)
let go !_ [] = return mv
go j (r:rs) = V.unsafeCopy (M.unsafeSlice (j*x) x mv) r >> go (j-1) rs
go (y - 1) vs
(!x,!y,!vs) = F.foldl' f (maxBound,0,[]) xss
f (!i,!j,!ss) xs = let !v = V.fromList (F.toList xs)
in (min i (V.length v), j+1, v : ss)
hmPoints :: IndexedTraversal' (V2 Int) HeatMatrix Double
hmPoints f (HeatMatrix e@(V2 x y) v a b) =
go 0 0 0 <&> \vs ->
let v'= V.fromListN (x*y) vs
in HeatMatrix e v' a b
where
go !s !i !j
| i >= x = go s 0 (j+1)
| j >= y = pure []
| otherwise = (:) <$> indexed f (V2 i j) (V.unsafeIndex v s)
<*> go (s+1) (i+1) j
{-# INLINE [0] hmPoints #-}
{-# RULES
"hmPoints/foldr"
hmPoints = ifoldring hmFold :: Getting (Endo r) HeatMatrix Double;
"hmPoints/ifoldr"
hmPoints = ifoldring hmFold :: IndexedGetting (V2 Int) (Endo r) HeatMatrix Double
#-}
hmFold :: (V2 Int -> Double -> b -> b) -> b -> HeatMatrix -> b
hmFold f b0 (HeatMatrix (V2 x y) v _ _) = go 0 0 0 b0 where
go !s !i !j b
| i >= x = go s 0 (j+1) b
| j >= y = b
| otherwise = f (V2 i j) (V.unsafeIndex v s) (go (s+1) (i+1) j b)
{-# INLINE hmFold #-}
pixelHeatRender :: HeatMatrix -> ColourMap -> Diagram V2
pixelHeatRender hm cm = alignBL $ imageEmb (ImageRGB8 img)
where
img = heatImage hm cm
pixelHeatRender' :: Int -> HeatMatrix -> ColourMap -> Diagram V2
pixelHeatRender' n hm cm =
scale (1/fromIntegral n) . alignBL $ imageEmb (ImageRGB8 img)
where
img = scaleImage n $ heatImage hm cm
scaleImage :: Int -> Image PixelRGB8 -> Image PixelRGB8
scaleImage n img | n == 1 = img
| n == 0 = Image 0 0 S.empty
| n < 0 = error "scaleImage: negative scale"
scaleImage n (Image x y v) = Image (n*x) (n*y) vn where
!refV = V.fromList $ map (*3) [ i + n*x*j | i <- [0..n-1], j <- [0..n-1] ]
!n3 = 3*n
vn = S.create $ do
mv <- M.new (n * n * S.length v)
let go !q !i !s | q >= 3*x*y = return mv
| i == x = go q 0 (s + 3*x*n*(n-1))
go q i s = do
let !r = S.unsafeIndex v q
!g = S.unsafeIndex v (q+1)
!b = S.unsafeIndex v (q+2)
V.forM_ refV $ \ds -> do
M.unsafeWrite mv (s + ds ) r
M.unsafeWrite mv (s + ds + 1) g
M.unsafeWrite mv (s + ds + 2) b
go (q+3) (i+1) (s+n3)
go 0 0 0
heatImage :: HeatMatrix -> ColourMap -> Image PixelRGB8
heatImage (HeatMatrix (V2 x y) dv a b) cm = Image x y v where
!cv = mkColourVector cm
v :: S.Vector Word8
v = S.create $ do
mv <- M.new (3 * x * y)
let !m = 256 / (b - a)
writeColour q (toSRGB24 -> RGB rr gg bb) = do
M.unsafeWrite mv q rr
M.unsafeWrite mv (q+1) gg
M.unsafeWrite mv (q+2) bb
colourValue !q !d
| isNaN d = writeColour q (cm^.nanColour)
| isInfinite d = writeColour q (if d < 0 then cm^.negInfColour else cm^.infColour)
| otherwise = do
let !o = 3 * (min 255 . max 0 . round $ (d - a) * m)
M.unsafeWrite mv q (V.unsafeIndex cv o )
M.unsafeWrite mv (q+1) (V.unsafeIndex cv (o+1))
M.unsafeWrite mv (q+2) (V.unsafeIndex cv (o+2))
go s i q
| s == x-1 = do
let d = V.unsafeIndex dv s
colourValue q d
return mv
| i == x = go (s - 2*x) 0 q
| otherwise = do
let d = V.unsafeIndex dv s
colourValue q d
go (s+1) (i+1) (q+3)
go (x * (y-1)) 0 0
mkColourVector :: ColourMap -> V.Vector Word8
mkColourVector cm = V.create $ do
mv <- M.new (3*256)
let go i | i == 3*256 = return mv
| otherwise = do
let x = fromIntegral i / (3*256)
RGB r g b = cm ^. ixColourR x . to toSRGB24
M.unsafeWrite mv i r
M.unsafeWrite mv (i+1) g
M.unsafeWrite mv (i+2) b
go (i+3)
go 0
pathHeatRender :: HeatMatrix -> ColourMap -> Diagram V2
pathHeatRender hm@(HeatMatrix _ _ a b) cm = ifoldMapOf hmPoints mk hm # lwO 0
where
normalise d = (d - a) / (b - a)
mk v@(V2 i j) d =
rect w h
# alignTR
# translate (fromIntegral <$> v ^+^ 1)
# fc (cm ^. ixColour (normalise d))
where
w | i == 0 = 1
| otherwise = 1.5
h | j == 0 = 1
| otherwise = 1.5
data HeatMap v = HeatMap
{ hMatrix :: HeatMatrix
, hStart :: P2 Double
, hSize :: V2 Double
, hGridSty :: Style v Double
, hGridVisible :: Bool
, hLimits :: Maybe (Double,Double)
, hDraw :: HeatMatrix -> ColourMap -> Diagram v
} deriving Typeable
type instance V (HeatMap v) = v
type instance N (HeatMap v) = Double
class HasHeatMap f a where
heatMapOptions :: LensLike' f a (HeatMap (V a))
heatMapGridVisible :: Functor f => LensLike' f a Bool
heatMapGridVisible = heatMapOptions . lens hGridVisible (\s b -> (s {hGridVisible = b}))
heatMapGridStyle :: Functor f => LensLike' f a (Style (V a) Double)
heatMapGridStyle = heatMapOptions . lens hGridSty (\s b -> (s {hGridSty = b}))
heatMapSize :: Functor f => LensLike' f a (V2 Double)
heatMapSize = heatMapOptions . lens hSize (\s b -> (s {hSize = b}))
heatMapExtent :: Functor f => LensLike' f a (V2 Double)
heatMapExtent = heatMapOptions . l where
l f hm = f (hSize hm * s) <&> \x -> hm { hSize = x / s }
where s = fmap fromIntegral (hmSize $ hMatrix hm)
heatMapStart :: Functor f => LensLike' f a (P2 Double)
heatMapStart = heatMapOptions . lens hStart (\s b -> (s {hStart = b}))
heatMapCentre :: Functor f => LensLike' f a (P2 Double)
heatMapCentre = heatMapOptions . l where
l f hm = f (hStart hm .+^ v) <&> \p -> hm { hStart = p .-^ v }
where v = fmap fromIntegral (hmSize $ hMatrix hm) * hSize hm / 2
heatMapLimits :: Functor f => LensLike' f a (Maybe (Double, Double))
heatMapLimits = heatMapOptions . lens hLimits (\s b -> (s {hLimits = b}))
heatMapRender :: Functor f => LensLike' f a (HeatMatrix -> ColourMap -> Diagram (V a))
heatMapRender = heatMapOptions . lens hDraw (\s b -> (s {hDraw = b}))
instance HasHeatMap f (HeatMap v) where
heatMapOptions = id
instance (Functor f, HasHeatMap f a) => HasHeatMap f (Plot a) where
heatMapOptions = rawPlot . heatMapOptions
instance Enveloped (HeatMap V2) where
getEnvelope hm = getEnvelope (fromCorners p (p .+^ v))
where p = view heatMapStart hm
v = view heatMapExtent hm
instance Enveloped (HeatMap V3) where
getEnvelope hm = getEnvelope (fromCorners (v3 p) (v3 $ p .+^ v))
where p = view heatMapStart hm
v = view heatMapExtent hm
v3 (P (V2 x y)) = P $ V3 x y 0
instance Plotable (HeatMap V2) where
renderPlotable s _sty HeatMap {..} =
transform (s^.specTrans) $
grid <> hDraw matrix' (s^.specColourMap)
# scaleV hSize
# moveTo hStart
where
grid = mempty
matrix' = case hLimits of
Just (a,b) -> hMatrix { hmBoundLower = a, hmBoundUpper = b }
Nothing -> hMatrix
defLegendPic sty HeatMap {..} = square 5 # applyAreaStyle sty
instance Plotable (HeatMap V3) where
renderPlotable s _sty HeatMap {..} =
transform (s^.specTrans) $
grid <> hDraw matrix' (s^.specColourMap)
# scaleV (v3 hSize)
# moveTo (p3 hStart)
where
v3 (V2 x y) = V3 x y 1
p3 (P (V2 x y)) = P $ V3 x y 0
grid = mempty
matrix' = case hLimits of
Just (a,b) -> hMatrix { hmBoundLower = a, hmBoundUpper = b }
Nothing -> hMatrix
defLegendPic sty HeatMap {..} = cube # scale 5 # applyAreaStyle sty
mkHeatMap :: HeatMatrix -> HeatMap V2
mkHeatMap mat = HeatMap
{ hMatrix = mat
, hStart = origin
, hSize = V2 1 1
, hGridSty = mempty
, hGridVisible = False
, hLimits = Nothing
, hDraw = pathHeatRender
}
mkHeatSurface :: HeatMatrix -> HeatMap V3
mkHeatSurface mat = HeatMap
{ hMatrix = mat
, hStart = origin
, hSize = V2 1 1
, hGridSty = mempty
, hGridVisible = False
, hLimits = Nothing
, hDraw = undefined
}
heatMap
:: (F.Foldable f,
F.Foldable g,
MonadState (Axis V2) m)
=> f (g Double)
-> State (Plot (HeatMap V2)) ()
-> m ()
heatMap xss s = do
let hm@(HeatMatrix _ _ a b) = mkHeatMatrix' xss
addPlotable (mkHeatMap hm) s
colourBarRange .= over both realToFrac (a,b)
heatSurface
:: (F.Foldable f,
F.Foldable g,
MonadState (Axis V3) m)
=> f (g Double)
-> State (Plot (HeatMap V3)) ()
-> m ()
heatSurface xss s = do
let hm@(HeatMatrix _ _ a b) = mkHeatMatrix' xss
addPlotable (mkHeatSurface hm) s
colourBarRange .= over both realToFrac (a,b)
heatMap'
:: (F.Foldable f,
F.Foldable g,
MonadState (Axis V2) m)
=> f (g Double)
-> m ()
heatMap' xss = heatMap xss (return ())
heatMapIndexed
:: (VectorLike V2 Int i,
MonadState (Axis V2) m)
=> i
-> (i -> Double)
-> State (Plot (HeatMap V2)) ()
-> m ()
heatMapIndexed i f s = do
let hm@(HeatMatrix _ _ a b) = mkHeatMatrix (view unvectorLike i) (f . view vectorLike)
addPlotable (mkHeatMap hm) s
colourBarRange .= over both realToFrac (a,b)
heatMapIndexed'
:: (VectorLike V2 Int i,
MonadState (Axis V2) m)
=> i
-> (i -> Double)
-> m ()
heatMapIndexed' i f = heatMapIndexed i f (return ())