content: {-# language ExtendedDefaultRules, ScopedTypeVariables, QuasiQuotes, ParallelListComp #-}
  
  import Test.Tasty
  import Test.Tasty.Runners
  import Test.Tasty.HUnit
  import Test.Tasty.ExpectedFailure
  import System.IO.Unsafe
  import Graphics.Matplotlib
  import System.IO.Temp
  import System.Random
  import Text.RawString.QQ
  import Data.List
  import Data.List.Split
  import Control.Monad
  import Test.Tasty.Golden.Advanced
  import qualified Data.ByteString as BS
  import System.Process
  import System.IO
  import Numeric.AD
  import System.Directory
  
  -- * Random values for testing
  
  uniforms :: (Random a, Num a) => [a]
  uniforms = randoms (mkStdGen 42)
  
  uniforms' lo hi = randomRs (lo,hi) (mkStdGen 42)
  
  -- * Not so random values to enable some fully-reproducible tests
  
  xs = [-0.54571992,  1.48409716, -0.57545561,  2.13058156, -0.75740497,
        -1.27879086, -0.96008858, -1.65482373, -1.69086194, -1.41925464,
         0.68144401,  1.44847131,  1.12124327,  1.32056244, -0.4555279 ,
         1.96002923, -1.34810771,  0.01513306,  1.25298883, -1.07541677,
         0.60920278, -0.13978163,  0.3975209 , -0.15211044,  0.0534633 ,
        -0.39498474, -1.15320186,  0.6795936 ,  0.50704333,  1.52951037,
         0.90867393, -0.24833115,  1.39816295, -0.28642814,  0.96951915,
        -1.20438266, -0.32642955, -0.62982941,  0.7245042 , -1.03169685,
        -0.00542761,  0.54247125, -1.11559132, -2.6829048 , -0.13370841,
        -0.74111166,  0.59198725,  2.73711931,  1.82122485, -0.73915212,
         0.88290489, -1.17307876,  0.06753304,  0.40150672,  1.54455801,
        -0.31133243,  1.66844302, -0.1290924 ,  0.89657699,  0.41181393,
         2.13382656,  1.58577659, -1.02449069, -1.10245954, -0.59691196,
        -0.63040161, -0.51541836,  0.04139408,  0.54203055, -2.09082082,
        -0.41295376, -0.77509336,  0.47612065, -1.69680096,  0.90195265,
         0.23798858, -0.05112783,  1.00645056, -0.67012513,  0.52017487,
        -0.42251472,  0.96513844,  1.00298933,  0.18257527,  0.54599979,
        -1.50321042,  0.03949817,  0.35286613,  1.86994544,  1.16249707,
         0.57421291,  1.21151469,  1.74863421,  0.42287859, -1.22785845,
        -0.61650528,  1.76743253, -0.45818694, -1.16560907,  0.0677502 ]
  ys = [ 1.28231455,  1.13480471,  0.57738712,  0.10268954,  1.00162163,
        -0.85453571, -1.61049343,  1.33194242,  0.12054943, -0.56568776,
         2.11550009,  0.03663454,  0.24889313,  0.85458325,  0.77326592,
         0.58815223, -0.79997005,  0.54979418,  0.47711544,  0.73004143,
        -0.65704545,  1.1946521 , -0.31119444, -0.0958055 ,  0.37838453,
         1.01281301, -0.53364162,  2.84609607,  0.09363483, -0.14821451,
        -0.0481863 , -3.58277731, -1.7168244 , -0.87526525, -0.65430073,
         1.0284506 , -0.81397895,  0.34868379, -0.51671293,  0.92879285,
         0.04099886,  1.0828335 ,  1.25991492, -1.48901447,  0.43657503,
         0.78191509,  0.16633587,  1.99411663, -0.25542794, -0.43377353,
        -0.82871869, -0.0402321 , -0.06278027,  0.28066445,  0.01185443,
         1.42640101, -0.16627931,  0.82021257, -0.66684095, -0.21289723,
        -1.25974667, -0.28681327, -2.11039334, -0.2722768 , -0.51622958,
         0.01324637, -0.29277708,  1.35916036, -0.09089638, -1.00619256,
         0.62707331,  1.17105748, -0.85636353, -0.6243519 ,  0.1720141 ,
        -0.15715394,  1.13488465,  2.43996937,  2.08224839, -0.23676918,
        -0.24924999,  1.21629376, -0.12748227,  0.78319565, -0.10528614,
         0.60177749, -1.03490762, -0.59163218,  0.16414076,  2.22783012,
        -0.55178235, -0.69915414,  1.35454045,  0.42931902, -1.33656935,
        -0.8023867 , -2.81354854,  0.39553427, -0.22235586, -1.34302011]
  
  -- * Generate normally distributed random values; taken from normaldistribution==1.1.0.3
  
  -- | Box-Muller method for generating two normally distributed
  -- independent random values from two uniformly distributed
  -- independent random values.
  boxMuller :: Floating a => a -> a -> (a,a)
  boxMuller u1 u2 = (r * cos t, r * sin t) where r = sqrt (-2 * log u1)
                                                 t = 2 * pi * u2
  
  -- | Convert a list of uniformly distributed random values into a
  -- list of normally distributed random values. The Box-Muller
  -- algorithms converts values two at a time, so if the input list
  -- has an uneven number of element the last one will be discarded.
  boxMullers :: Floating a => [a] -> [a]
  boxMullers (u1:u2:us) = n1:n2:boxMullers us where (n1,n2) = boxMuller u1 u2
  boxMullers _          = []
  
  -- | Plural variant of 'normal', producing an infinite list of
  -- random values instead of returning a new generator. This function
  -- is analogous to 'Random.randoms'.
  normals = boxMullers $ randoms (mkStdGen 42)
  
  -- | Analogous to 'normals' but uses the supplied (mean, standard
  -- deviation).
  normals' (mean, sigma) g = map (\x -> x * sigma + mean) $ normals
  
  pdfBivariateNormal x y sigmax sigmay mux muy sigmaxy =
    1/(2*pi*sigmax*sigmay*(sqrt(1-rho^2)))*exp(-z/(2*(1-rho^2)))
    where rho = sigmaxy/(sigmax*sigmay)
          z = (x-mux)^2/sigmax^2-(2*rho*(x-mux)*(y-muy))/(sigmax*sigmay)+(y-muy)^2/sigmay^2
  
  -- * Tests
  
  main = do
    createDirectoryIfMissing False "/tmp/imgs/"
    defaultMain $ tests "All tests" testPlot
  
  main' = do
    createDirectoryIfMissing False "/tmp/imgs/"
    defaultMain $ tests "Golden tests" testPlotGolden
  
  main'' = do
    createDirectoryIfMissing False "/tmp/imgs/"
    defaultMain $ testGroup "All tests" [tests "Execution tests" testPlot
                                        , toneDownTests "Unreliable across machines" $ tests "Golden tests" testPlotGolden]
  
  tests name f = testGroup name   [basicTests f,
                                   toneDownTests "Can fail with old matplotlib or data" $ fragileTests f,
                                   ignoreTest $ failingTests f]
  
  toneDownTests reason tests = wrapTest (liftM (\x -> x { resultOutcome = Success
                                                       , resultShortDescription =
                                                         case resultOutcome x of
                                                           Success -> resultShortDescription x
                                                           _ -> reason ++ ": " ++ resultShortDescription x
                                                       })) tests
  
  testPlotGolden name fn =
           unsafePerformIO $ tmp (\filename ->
                                   return $ goldenTest
                                     name
                                     (BS.readFile ref)
                                     (file filename fn >> BS.readFile filename)
                                     (\g n ->
                                         tmp (\gfile ->
                                                 tmp (\nfile -> do
                                                         BS.writeFile gfile g
                                                         BS.writeFile nfile n
                                                         (code, stdout, stderr) <-
                                                           readProcessWithExitCode "/usr/bin/compare" ["-metric"
                                                                                                      ,"PSNR"
                                                                                                      ,gfile
                                                                                                      ,nfile
                                                                                                      ,"null"] ""
                                                         case (stderr, reads stderr) of
                                                           ("inf", _) -> return Nothing
                                                           (_, [(x :: Double, _)]) ->
                                                             if x < 25 then
                                                               return $ Just $ "Images very different; PSNR too low " ++ show x else
                                                               return Nothing)))
                                     (BS.writeFile ref))
    where ref = "imgs/" ++ name ++ ".png"
          tmp f = withSystemTempFile "a.png" (\filename h -> hClose h >> f filename)
  
  -- | Test one plot; right now we just test that the command executed without
  -- errors. We should visually compare plots somehow.
  testPlot' name fn = testCase name $ tryit fn @?= Right ""
    where tryit fn = unsafePerformIO $ withSystemTempFile "a.png" (\filename _ -> file filename fn)
  
  -- | This generates examples from the test cases
  testPlot name fn = testCase name $ tryit fn name @?= Right ""
    where tryit fn name = unsafePerformIO $ do
            file ("/tmp/imgs/" ++ name ++ ".png") fn
  
  basicTests f = testGroup "Basic tests"
    [ f "histogram" mhistogram
    , f "cumulative" mcumulative
    , f "scatter" mscatter
    , f "contour" mcontour
    , f "labelled-histogram" mlabelledHistogram
    , f "density-bandwidth" mdensityBandwidth
    , f "density" mdensity
    , f "line-function" mlineFunction
    , f "quadratic" mQuadratic
    , f "projections" mProjections
    , f "line-options" mlineOptions
    , f "corr" mxcorr
    , f "show-matrix" mmat
    , f "legend" mlegend
    , f "hist2DLog" mhist2DLog
    , f "eventplot" meventplot
    , f "errorbar" merrorbar
    , f "scatterhist" mscatterHist
    , f "histMulti" mhistMulti
    , f "spines" mspines
    , f "hists" mhists
    , f "hinton" mhinton
    , f "integral" mintegral
    , f "quiver" mquiver
    , f "quiver-fancy" mquiverFancy
    , f "polar" mpolar
    , f "bivariateNormal" mbivariateNormal
    , f "pcolorlog" mpcolorlog
    , f "pie" mpie
    , f "stacked" mstacked
    , f "annotation" mannotation
    , f "streamplot" mstreamplot
    ]
  
  fragileTests f = testGroup "Fragile tests"
    [ f "tex" mtex -- TODO Fails on circle ci (with latex)
      -- TODO Fails on circle ci (labels is not valid; matplotlib too old)
    , f "boxplot" mboxplot
      -- TODO Fails on circle ci (no violin plots; matplotlib too old)
    , f "violinplot" mviolinplot
      -- TODO Needs a fairly recent matplotlib; too old for circleci
    , f "griddata" mgriddata
      -- TODO Needs access to a data file
    , f "images" mimages
    ]
  
  failingTests f = testGroup "Failing tests"
    [
      -- TODO This test case is broken
      f "sub-bars" msubBars
    ]
  
  -- * These tests are fully-reproducible, the output must be identical every time
  
  mhistogram :: Matplotlib
  mhistogram = histogram xs 8
  
  mcumulative = histogram xs 8 @@ [o2 "cumulative" True]
  
  mscatter = scatter xs ys
  
  degreesRadians a = a * pi / 180.0
  radiansDegrees a = a * 180.0 / pi
  
  mcontour = contourF (\a b -> sin (degreesRadians a) + cos (degreesRadians b)) (-100) 100 (-200) 200 10
  
  mlabelledHistogram = histogram xs 7
     % ylabel "number of queries"
     % xlabel "true positives"
  
  -- TODO Broken test
  msubBars = subplotBarsLabelled
         [[40, 50, 20, 50], [10, 20, 30, 40], [40, 50, 20, 50]]
         ["a", "b", "c", "d"] []
        % title "Wee a title"
        % xlabel "X"
        % ylabel "Y"
  
  mdensityBandwidth = densityBandwidth [2.1, 1.3, 0.4, 1.9, 5.1, 6.2] 1.5 (Just (-6, 11))
    % ylim 0 0.2
  
  mdensity = density [2.1, 1.3, 0.4, 1.9, 5.1, 6.2] (Just (-6, 11))
  
  mlineFunction = lineF (\x -> x**2) [0,0.01..1]
  
  mQuadratic = plotMapLinear (\x -> x**2) (-2) 2 100 @@ [o1 "."] % title "Quadratic function"
  
  mProjections = projectionsF (\a b -> cos (degreesRadians a) + sin (degreesRadians b)) (-100) 100 (-200) 200 10
  
  mlineOptions = plot [1,2,3,4,5,6] [1,3,2,5,2,4] @@ [o1 "go-", o2 "linewidth" 2]
  
  -- * These tests can be random and may not be exactly the same every time
  
  -- | http://matplotlib.org/examples/pylab_examples/xcorr_demo.html
  mxcorr = xacorr xs ys [o2 "usevlines" True, o2 "maxlags" 50, o2 "normed" True, o2 "lw" 2]
    where (xs :: [Double]) = take 100 normals
          (ys :: [Double]) = take 100 normals
  
  -- | http://matplotlib.org/examples/pylab_examples/tex_unicode_demo.html
  mtex = setTeX True
    % figure
    % addSubplot 1 1 1
    % plotMapLinear cos 0 1 100
    % xlabel [r|\textbf{time (s)}|]
    % ylabel "\\textit{Velocity (°/sec)}" @@ [o2 "fontsize" 16]
    % title [r|\TeX\ is Number $\displaystyle\sum_{n=1}^\infty\frac{-e^{i\pi}}{2^n}$!"|] @@ [o2 "fontsize" 16, o2 "color" "r"]
    % grid True
  
  mmat = pcolor (take 10 $ chunksOf 8 uniforms) @@ [o2 "edgecolors" "k", o2 "linewidth" 1]
  
  -- | http://matplotlib.org/examples/pylab_examples/legend_demo3.html
  mlegend = plotMapLinear (\x -> x ** 2) 0 1 100 @@ [o2 "label" "x^2"]
    % plotMapLinear (\x -> x ** 3) 0 1 100 @@ [o2 "label" "x^3"]
    % legend @@ [o2 "fancybox" True, o2 "shadow" True, o2 "title" "Legend", o2 "loc" "upper left"]
  
  -- | http://matplotlib.org/examples/pylab_examples/hist2d_log_demo.html
  mhist2DLog = histogram2D x y @@ [o2 "bins" 40, o2 "norm" $ lit "mcolors.LogNorm()"]
    % setAx
    % colorbar
    where (x:y:_) = chunksOf 10000 normals
  
  meventplot = plot xs ys
    % mp # "ax.add_collection(mcollections.EventCollection(data[0], linelength=0.05))"
    % mp # "ax.add_collection(mcollections.EventCollection(data[1], orientation='vertical', linelength=0.05))"
    % text 0.1 0.6 "Ticks mark the actual data points"
    where xs = sort $ take 10 uniforms
          ys = map (\x -> x ** 2) xs
  
  merrorbar = errorbar xs ys (Nothing :: Maybe [Double]) (Just errs) @@ [o2 "errorevery" 2]
    where xs = [0.1,0.2..4]
          ys = map (\x -> exp $ -x) xs
          errs = [map (\x -> 0.1 + 0.1 * sqrt x) xs, map (\x -> 0.1 + 0.1 * sqrt x) ys]
  
  mboxplot = subplots @@ [o2 "ncols" 2, o2 "sharey" True]
    % setSubplot "0"
    % boxplot (take 3 $ chunksOf 10 $ map (* 2) $ normals) @@ [o2 "labels" ["X", "Y", "Z"]]
    % setSubplot "1"
    % boxplot (take 3 $ chunksOf 10 $ map (* 2) $ normals) @@ [o2 "labels" ["A", "B", "C"], o2 "showbox" False, o2 "showcaps" False]
  
  mviolinplot = subplots @@ [o2 "ncols" 2, o2 "sharey" True]
    % setSubplot "0"
    % violinplot (take 3 $ chunksOf 100 $ map (* 2) $ normals)
    % setSubplot "1"
    % violinplot (take 3 $ chunksOf 100 $ map (* 2) $ normals) @@ [o2 "showmeans" True, o2 "showmedians" True, o2 "vert" False]
  
  -- | http://matplotlib.org/examples/pylab_examples/scatter_hist.html
  mscatterHist = figure @@ [o1 0]
    % setSizeInches 8 8
    -- The scatter plot
    % axes @@ [o1 ([left, bottom', width, height] :: [Double])]
    % scatter  x y
    % xlim (-lim) lim
    % ylim (-lim) lim
    -- The histogram on the right (x)
    % axes @@ [o1 [left, bottom_h, width, 0.2]]
    % mp # "ax.xaxis.set_major_formatter(mticker.NullFormatter())"
    % histogram x bins
    % xlim (-lim) lim
    -- The histogram on top (y)
    % axes @@ [o1 [left_h, bottom', 0.2, height]]
    % mp # "ax.yaxis.set_major_formatter(mticker.NullFormatter())"
    % histogram y bins @@ [o2 "orientation" "horizontal"]
    % ylim (-lim) lim
    where left = 0.1
          width = 0.65
          bottom' = 0.1
          height = 0.65
          bottom_h = left + width + 0.02
          left_h = left + width + 0.02
          [x, y] = take 2 $ chunksOf 1000 $ map (* 2) $ normals
          binwidth = 0.25
          xymax = maximum [maximum $ map abs x, maximum $ map abs y]
          lim = ((fromIntegral $ round $ xymax / binwidth) + 1) * binwidth
          bins = [-lim,-lim+binwidth..(lim + binwidth)]
  
  mhistMulti = subplots @@ [o2 "nrows" 2, o2 "ncols" 2]
    % setSubplot 0
    % histogram x nrBins @@ [o2 "density" 1, o2 "histtype" "bar", o2 "color" ["red", "tan", "lime"], o2 "label" ["red", "tan", "lime"]]
    % legend @@ [o2 "prop" $ lit "{'size': 10}"]
    % title "bars with legend"
    % setSubplot 1
    % histogram x nrBins @@ [o2 "density" 1, o2 "histtype" "bar", o2 "stacked" True]
    % title "stacked bar"
    % setSubplot 2
    % histogram x nrBins @@ [o2 "histtype" "step", o2 "stacked" True, o2 "fill" False]
    % title "stacked bar"
    % setSubplot 3
    % histogram (map (\x -> take x normals) [2000, 5000, 10000]) nrBins @@ [o2 "histtype" "bar"]
    % title "different sample sizes"
    % tightLayout
    where nrBins = 10
          x = take 3 $ chunksOf 1000 $ normals
  
  mspines = plot x y @@ [o1 "k--"]
    % plot x y' @@ [o1 "ro"]
    % xlim 0 (2 * pi)
    % xticks [0 :: Double, pi, 2*pi]
    % xtickLabels (map raw ["0", "$\\pi$", "2$\\pi$"])
    % ylim (-1.5) 1.5
    % yticks [-1 :: Double, 0, 1]
    % spine "left"
    % spineSetBounds (-1) 1
    % spine "right"
    % spineSetVisible False
    % spine "top"
    % spineSetVisible False
    % axisYTicksPosition "left"
    % axisXTicksPosition "bottom"
    where x  = mapLinear (\x -> x) 0 (2 * pi) 50
          y  = map sin x
          y' = zipWith (\a b -> a + 0.1*b) y normals
  
  mhists = h 10 1.5
         % h 4 1
         % h 15 2
         % h 6 0.5
    where ns mu var = map (\x -> mu + x * var) $ take 1000 normals
          h mu var = histogram (ns mu var) 25 @@ [o2 "histtype" "stepfilled"
                                               ,o2 "alpha" 0.8
                                               ,o2 "density" True]
  
  mhinton = mp # "ax.patch.set_facecolor('gray')"
    % setAspect @@ [o1 "equal", o1 "box"]
    % mp # "ax.xaxis.set_major_locator(plot.NullLocator())"
    % mp # "ax.yaxis.set_major_locator(plot.NullLocator())"
    % foldl (\a (x,y,w) -> a % f x y w) mp m
    % mp # "ax.autoscale_view()"
    % mp # "ax.invert_yaxis()"
    where m = [ (x,y,w) | x <- [0..19], y <- [0..19] | w <- (map (\x -> x - 0.5) normals) ]
          maxWeight = maximum $ map (\(_,_,v) -> abs v) m
          f x y w = mp # "ax.add_patch(plot.Rectangle("
                       # "[" # (x - size / 2) # "," # (y - size / 2) # "]"
                       # ", " # size # ", " # size
                       # ", facecolor='" # color # "', edgecolor='" # color # "'))"
            where color = if w > 0 then "white" else "black"
                  size  = sqrt $ abs w / maxWeight
  
  mintegral = subplots
    % plot x y @@ [o1 "r", o2 "linewidth" 2]
    % ylim 0 (maximum y)
    % mp # "ax.add_patch(plot.Polygon(" # ([(a, 0)] ++ zip ix iy ++ [(b,0)]) ## "))"
      @@ [o2 "facecolor" "0.9", o2 "edgecolor" "0.5"]
    % text (0.5 * (a + b)) 30 [r|$\int_a^b f(x)\mathrm{d}x$|]
      @@ [o2 "horizontalalignment" "center", o2 "fontsize" 20]
    % figText 0.9 0.05 "$x$"
    % figText 0.1 0.9  "$y$"
    % spine "right"
    % spineSetVisible False
    % spine "top"
    % spineSetVisible False
    % axisXTicksPosition "bottom"
    % xticks (a, b)
    % xtickLabels (raw "$a$", raw "$b$")
    % yticks ([] :: [Double])
    where func x = (x - 3) * (x - 5) * (x - 7) + 85
          -- integral limits
          a = 2 :: Double
          b = 9 :: Double
          (x :: [Double]) = mapLinear (\x -> x) 0 10 100
          y = map func x
          -- shaded region
          (ix :: [Double]) = mapLinear (\x -> x) a b 100
          iy = map func ix
  
  mquiver = quiver x y u v (Nothing :: Maybe [Double]) @@ [o2 "units" "width"]
    % quiverKey 0.9 0.93 2 (raw [r|$2 \frac{m}{s}$|])
      @@ [o2 "labelpos" "E", o2 "coordinates" "figure"]
    % xlim (-0.2) 6.4
    % ylim (-0.2) 6.4
    where m = [ (x,y,cos x,sin y) | x <- [0,0.2..2*pi], y <- [0,0.2..2*pi] ]
          x = map (\(x,_,_,_) -> x) m
          y = map (\(_,x,_,_) -> x) m
          u = map (\(_,_,x,_) -> x) m
          v = map (\(_,_,_,x) -> x) m
  
  mquiverFancy = quiver x y u v (Just mag) @@ [o2 "units" "x"
                                              ,o2 "pivot" "tip"
                                              ,o2 "width" 0.022
                                              ,o2 "scale" (1 / 0.15)]
    % quiverKey 0.9 0.93 1 (raw [r|$2 \frac{m}{s}$|])
      @@ [o2 "labelpos" "E", o2 "coordinates" "figure"]
    % scatter x y @@ [o2 "color" "k", o2 "s" 5]
    % xlim (-0.2) 6.4
    % ylim (-0.2) 6.4
    where m = [ (x,y,cos x,sin y) | x <- [0,0.2..2*pi], y <- [0,0.2..2*pi] ]
          x = map (\(x,_,_,_) -> x) m
          y = map (\(_,x,_,_) -> x) m
          u = map (\(_,_,x,_) -> x) m
          v = map (\(_,_,_,x) -> x) m
          mag = zipWith (\x x' -> sqrt(x**2 + x'**2)) u v
  
  mpolar = rc "grid" @@ [o2 "color" "#316931", o2 "linewidth" 1, o2 "linestyle" "-"]
    % rc "xtick" @@ [o2 "labelsize" 15]
    % rc "ytick" @@ [o2 "labelsize" 15]
    % figure @@ [o2 "figsize" (8::Int,8::Int)]
    % addAxes @@ [o1 [0.1, 0.1, 0.8, 0.8::Double], o2 "projection" "polar"
                 -- TODO My matplotlib doesn't seem to have this property
                 -- , o2 "facecolor" "#d5de9c"
                 ]
    % plot theta r @@ [o2 "color" "#ee8d18", o2 "lw" 3, o2 "label" "a line"]
    % plot (map (\x -> 0.5*x) theta) r
           @@ [o2 "color" "blue", o2 "ls" "--", o2 "lw" 3, o2 "label" "another line"]
    % legend
    where r = [0,0.01..3.0]
          theta = map (\x -> 2*pi*x) r
  
  mbivariateNormal =
    imshow vs @@ [o2 "interpolation" "bilinear"
                 ,o2 "cmap" $ raw "RdYlGn"
                 ,o2 "origin" "lower"
                 ,o2 "extent" [-3::Double, 3, -3, 3]
                 ,o2 "vmin" $ (0-) $ maximum $ map abs vs'
                 ,o2 "vmax" $ maximum $ map abs vs']
    where delta = 0.025::Double
          xs = [-3.0,-3.0+delta..3.0]
          ys = [-3.0,-3.0+delta..3.0]
          vs = [[pdfBivariateNormal x y 1.5 0.5 1.0 1.0 0.0
                 - pdfBivariateNormal x y 1.0 1.0 0.0 0.0 0.0
                | x <- xs]
                | y <- ys]
          vs' = foldl' (++) [] vs
  
  -- TODO This is subtly broken
  mimages = -- figure @@ [o2 "figsize" (10::Int,10::Int)]
    -- %
    subplots @@ [o2 "nrows" 1, o2 "ncols" 2]
    % setSubplot 0
    % imshow "data/heightmap.png" @@ [o2 "interpolation" "nearest"]
    % setSubplot 1
    % mp # "ls = mcolors.LightSource(azdeg=315, altdeg=45)"
    % mp # "ax.imshow(ls.shade(img, cmap=cm.gist_earth))"
    -- TODO This doesn't work on my matplab version
    -- vert_exag=0.05,
    --, blend_mode='overlay'
    % xlabel "overlay blend mode"
  
  mpcolorlog = figure
    % addSubplot 2 1 1
    % pcolor3 xs' ys' vs @@ [o2 "cmap" $ raw "PuBu_r"]
    % colorbar
    % addSubplot 2 1 2
    % pcolor3 xs' ys' vs @@ [o2 "norm"
                              (lit $ "mcolors.LogNorm(vmin="++(show $ minimum vs')++
                                ", vmax="++(show $ maximum vs')++")")
                            ,o2 "cmap" $ raw "PuBu_r"]
    % colorbar
    where delta = 0.1::Double
          xs = [-3.0,-3.0+delta..3.0]
          ys = [-3.0,-3.0+delta..3.0]
          vs = [[pdfBivariateNormal x y 0.1 0.2 1.0 1.0 0.0
                 + 0.1 * pdfBivariateNormal x y 1.0 1.0 0.0 0.0 0.0
                | x <- xs]
                | y <- ys]
          xs' = [[ x | x <- xs]| y <- ys]
          ys' = [[ y | x <- xs]| y <- ys]
          vs' = foldl' (++) [] vs
  
  mpie = pie [15, 30, 45, 10 :: Double]
    @@ [o2 "explode" [0, 0.05, 0, 0 :: Double]
       ,o2 "labels" ["Frogs", "Hogs", "Dogs", "Logs"]
       ,o2 "autopct" "%.0f%%"
       ,o2 "shadow" True]
  
  -- | http://matplotlib.org/examples/pylab_examples/bar_stacked.html
  mstacked =
    -- TODO The locations of the bars is off
      bar [0..4] ms @@ [o1 width, o2 "color" "#d62728", o2 "yerr" mStd, o2 "label" "ms"]
    % bar [0..4] ws @@ [o1 width, o2 "bottom" ms,       o2 "yerr" wStd, o2 "label" "ws"]
    % xticks [0..4 :: Int]
    % xtickLabels "['G1', 'G2', 'G3', 'G4', 'G5']"
    % title "Scores"
    % ylabel "Score"
    % yticks [0,10..80 :: Int]
    % legend
    where ms    = [20 :: Double, 35, 30, 35, 27]
          ws    = [25 :: Double, 32, 34, 20, 25]
          mStd  = [2 :: Double, 3, 4, 1, 2]
          wStd  = [3 :: Double, 5, 2, 3, 3]
          width = 0.35 :: Double
  
  mannotation = -- figure @@ [o2 "figsize" (10::Int,10::Int)]
    -- TODO This is subtly broken
    -- TODO Dictionaries
      plot t s @@ [o2 "lw" 3]
    % xlim (-1) 5
    % ylim (-4) 3
    % annotate "straight" @@ [o2 "xy" [0, 1::Double], o2 "xycoords" "data", o2 "xytext" [-50, 30 :: Double]
                             ,o2 "textcoords" "offset points", o2 "arrowprops" (lit "dict(arrowstyle='->')")]
    % annotate "arc3,\\nrad 0.2" @@ [o2 "xy" [0.5, -1::Double], o2 "xycoords" "data", o2 "xytext" [-80, -60 :: Double]
                                   ,o2 "textcoords" "offset points"
                                   ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='arc3,rad=.2')")]
    % annotate "arc,\\nangle 50" @@ [o2 "xy" [1, 1::Double], o2 "xycoords" "data", o2 "xytext" [-90, 50 :: Double]
                                    ,o2 "textcoords" "offset points"
                                    ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='arc,angleA=0,armA=50,rad=10')")]
    % annotate "arc,\\narms" @@ [o2 "xy" [1.5, -1::Double], o2 "xycoords" "data", o2 "xytext" [-80, -60 :: Double]
                                ,o2 "textcoords" "offset points"
                                ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='arc,angleA=0,armA=40,angleB=-90,armB=30,rad=7')")]
    % annotate "angle,\\nangle 90" @@ [o2 "xy" [2, 1::Double], o2 "xycoords" "data", o2 "xytext" [-70, 30 :: Double]
                     ,o2 "textcoords" "offset points"
                     ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='angle,angleA=0,angleB=90,rad=10')")]
    % annotate "angle3,\\nangle -90" @@ [o2 "xy" [2.5, -1::Double], o2 "xycoords" "data", o2 "xytext" [-80, -60 :: Double]
                     ,o2 "textcoords" "offset points"
                     ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='angle3,angleA=0,angleB=-90')")]
    % annotate "angle,\\nround" @@ [o2 "xy" [3, 1::Double], o2 "xycoords" "data", o2 "xytext" [-60, 30 :: Double]
                     ,o2 "textcoords" "offset points"
                     ,o2 "bbox" (lit "dict(boxstyle='round', fc='0.8')")
                     ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='angle,angleA=0,angleB=90,rad=10')")]
    % annotate "angle,\\nround4" @@ [o2 "xy" [3.5, -1::Double], o2 "xycoords" "data", o2 "xytext" [-70, -80 :: Double]
                     ,o2 "textcoords" "offset points"
                     ,o2 "bbox" (lit "dict(boxstyle='round4,pad=.5', fc='0.8')")
                     ,o2 "size" 20
                     ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='angle,angleA=0,angleB=-90,rad=10')")]
    % annotate "angle,\\nshrink" @@ [o2 "xy" [4, 1::Double], o2 "xycoords" "data", o2 "xytext" [-60, 30 :: Double]
                     ,o2 "textcoords" "offset points"
                     ,o2 "bbox" (lit "dict(boxstyle='round', fc='0.8')")
                     ,o2 "arrowprops" (lit "dict(arrowstyle='->', connectionstyle='angle,angleA=0,angleB=90,rad=10')")]
    -- TODO This annotation doesn't render correctly on my matplotlib version
    % annotate "" @@ [o2 "xy" [4, 1::Double], o2 "xycoords" "data", o2 "xytext" [4.5, -1 :: Double]
                     ,o2 "textcoords" "offset points"
                     ,o2 "arrowprops" (lit "dict(arrowstyle='<->', connectionstyle='bar', ec='k', shrinkA=5, shrinkB=5)")]
    where t = [0, 0.01 .. 5.0 :: Double]
          s = map (\x -> cos $ 2*pi*x) t
  
  mstreamplot = streamplot xs ys xs' ys' @@ [o2 "linewidth" mag']
    -- useful for seeing the energy landscape
    -- pcolor3 xmat ymat vs
    where delta = 0.05::Double
          xs = [-3.0,-3.0+delta..3.0]
          ys = [-3.0,-3.0+delta..3.0]
          xmat = [[ x | x <- xs]| y <- ys]
          ymat = [[ y | x <- xs]| y <- ys]
          ms = [[grad' (\[xv,yv] ->
                           pdfBivariateNormal xv yv 0.4 0.7 1.0 1.0 0.0
                           + pdfBivariateNormal xv yv 1.0 1.0 (-1.0) (-1.0) 0.0)
                  [x,y]
                | x <- xs]
               | y <- ys] :: [[(Double, [Double])]]
          vs  = map2 (\(v, _)      -> v) ms
          xs' = map2 (\(_, [x, _]) -> x) ms
          ys' = map2 (\(_, [_, y]) -> y) ms
          mag' = zipWith (\lx ly -> zipWith (\x y -> 5 * (log $ 1 + (sqrt $ x*x + y*y))) lx ly) xs' ys'
          map2 f l = map (\r -> map f r) l
  
  mgriddata = readData (x, y, z, xi, yi)
    -- TODO This requires a lot of manual indexing. Next big API change will be to
    -- have references to loaded data.
    % mp # "data.append(interpolate.griddata((data[0], data[1]), data[2], tuple(np.meshgrid(data[3], data[4])), method='linear', rescale=True))"
    % mp # "plot.contour(data[3], data[4], data[5], 15, linewidths=0.5, colors='k')"
    % mp # "plot.contourf(data[3], data[4], data[5], 15, extend='both')"
    % colorbar
    % scatter x y @@ [o2 "marker" "o", o2 "s" 5, o2 "zorder" 10]
    % xlim (-2) 2
    % ylim (-2) 2
    % title "Grid interpolation"
    where [x, y] = take 2 $ chunksOf 200 $ map (\x -> 4 * (x - 0.5)) $ uniforms
          z = zipWith (\x y -> x*(exp $ -(x**2) - y**2)) x y
          xi = mapLinear (\x -> x) (-2.1) 2.1 300
          yi = mapLinear (\x -> x) (-2.1) 2.1 300
  