diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,6 +1,7 @@
 # Matplotlib
 
-[![Build Status](http://circleci-badges-max.herokuapp.com/img/abarbu/matplotlib-haskell/master?token=468e8942459ca5f34089fb5c29a478ffb6d531af)](https://circleci.com/gh/abarbu/matplotlib-haskell/tree/master)
+[![Build Status](https://img.shields.io/circleci/project/github/abarbu/matplotlib-haskell.svg)](circleci.com/gh/abarbu/matplotlib-haskell)
+[![Hackage](https://img.shields.io/hackage/v/matplotlib.svg)](https://hackage.haskell.org/package/matplotlib)
 
 Haskell bindings to Python's Matplotlib. It's high time that Haskell had a
 fully-fledged plotting library!
@@ -25,7 +26,7 @@
 on Ubuntu machines with the following command:
 
 ```bash
-sudo apt-get install -y python3-matplotlib python3-numpy python-mpltoolkits.basemap
+sudo apt-get install -y python3-pip python3-matplotlib python3-numpy python-mpltoolkits.basemap
 ```
 
 If you have instructions for other machines or OSes let me know. We require
diff --git a/matplotlib.cabal b/matplotlib.cabal
--- a/matplotlib.cabal
+++ b/matplotlib.cabal
@@ -1,5 +1,5 @@
 name:                matplotlib
-version:             0.2.1
+version:             0.3.0
 synopsis:            Bindings to Matplotlib; a Python plotting library
 description:
     Matplotlib is probably the most full featured plotting library out there.
@@ -41,6 +41,9 @@
                      , tasty
                      , tasty-hunit
                      , temporary
+                     , random
+                     , raw-strings-qq
+                     , split
   ghc-options:         -threaded -rtsopts -with-rtsopts=-N
   default-language:    Haskell2010
 
diff --git a/src/Graphics/Matplotlib.hs b/src/Graphics/Matplotlib.hs
--- a/src/Graphics/Matplotlib.hs
+++ b/src/Graphics/Matplotlib.hs
@@ -65,15 +65,117 @@
 onscreen m = withMplot m (\str -> python $ pyIncludes ++ str ++ pyDetach ++ pyOnscreen)
 
 -- | Print the python code that would be executed
-code :: Matplotlib -> IO (Either a String)
-code m = withMplot m (\str -> return $ Right $ unlines $ pyIncludes ++ str ++ pyDetach ++ pyOnscreen)
+code :: Matplotlib -> IO String
+code m = withMplot m (\str -> return $ unlines $ pyIncludes ++ str ++ pyDetach ++ pyOnscreen)
 
 -- | Save to a file
 figure :: [Char] -> Matplotlib -> IO (Either String String)
 figure filename m = withMplot m (\str -> python $ pyIncludes ++ str ++ pyFigure filename)
 
--- * Plotting commands
+-- * Useful plots
 
+-- | Plot the cross-correlation and autocorrelation of several variables. TODO Due to
+-- a limitation in the options mechanism this takes explicit options.
+xacorr xs ys opts = readData (xs, ys)
+  % addSubplot 2 1 1
+  % xcorr xs ys @@ opts
+  % grid True
+  % axhline 0 @@ [o1 "0", o2 "color" "'black'", o2 "lw" "2"]
+  % addSubplot 2 1 2 @@ [o2 "sharex" "ax"]
+  % acorr xs @@ opts
+  % grid True
+  % axhline 0 @@ [o2 "color" "'black'", o2 "lw" "2"]
+
+-- | Plot a histogram for the given values with 'bins'
+histogram :: (MplotValue val, ToJSON t) => t -> val -> Matplotlib
+histogram values bins = readData [values] % dataHistogram 0 bins
+
+-- | Plot a 2D histogram for the given values with 'bins'
+histogram2D x y = readData [x,y] %
+  mp # "plot.hist2d(data[0], data[1]" ## ")"
+
+-- | Plot the given values as a scatter plot
+scatter :: (ToJSON t1, ToJSON t) => t1 -> t -> Matplotlib
+scatter x y = plot x y `def` [o1 "'.'"]
+
+-- | Plot a line
+line :: (ToJSON t1, ToJSON t) => t1 -> t -> Matplotlib
+line x y = plot x y `def` [o1 "'-'"]
+
+-- | Like 'plot' but takes an error bar value per point
+errorbar xs ys errs = readData (xs, ys, errs)
+  % mp # "ax.errorbar(data[0], data[1], yerr=data[2]" ## ")"
+
+-- | Plot a line given a function that will be executed for each element of
+-- given list. The list provides the x values, the function the y values.
+lineF :: (ToJSON a, ToJSON b) => (a -> b) -> [a] -> Matplotlib
+lineF f l = plot l (map f l) `def` [o1 "'-'"]
+
+boxplot l = readData l
+  % mp # "ax.boxplot(data" ## ")"
+
+-- | Given a grid of x and y values and a number of steps call the given
+-- function and plot the 3D contour
+contourF :: (ToJSON val, MplotValue val, Ord val) => (Double -> Double -> val) -> Double -> Double -> Double -> Double -> Double -> Matplotlib
+contourF f xStart xEnd yStart yEnd steps = contour xs ys zs
+  where xs = mapLinear (\x -> (mapLinear (\_ -> x) yStart yEnd steps)) xStart xEnd steps
+        ys = mapLinear (\_ -> (mapLinear (\y -> y) yStart yEnd steps)) xStart xEnd steps
+        zs = mapLinear (\x -> (mapLinear (\y -> f x y) yStart yEnd steps)) xStart xEnd steps
+
+-- | Given a grid of x and y values and a number of steps call the given
+-- function and plot the 3D projection
+projectionsF :: (ToJSON val, MplotValue val, Ord val) => (Double -> Double -> val) -> Double -> Double -> Double -> Double -> Double -> Matplotlib
+projectionsF f xStart xEnd yStart yEnd steps = projections xs ys zs
+  where xs = mapLinear (\x -> (mapLinear (\_ -> x) yStart yEnd steps)) xStart xEnd steps
+        ys = mapLinear (\_ -> (mapLinear (\y -> y) yStart yEnd steps)) xStart xEnd steps
+        zs = mapLinear (\x -> (mapLinear (\y -> f x y) yStart yEnd steps)) xStart xEnd steps
+
+-- | Plot x against y interpolating with n steps
+plotInterpolated :: (MplotValue val, ToJSON t, ToJSON t1) => t1 -> t -> val -> Matplotlib
+plotInterpolated x y n =
+  readData (x, y)
+  % interpolate 0 1 n
+  % dataPlot 0 1 `def` [o1 "-"]
+
+-- | A handy function to plot a line between two points give a function and a number o steps
+plotMapLinear :: ToJSON b => (Double -> b) -> Double -> Double -> Double -> Matplotlib
+plotMapLinear f s e n = line xs ys
+  where xs = mapLinear (\x -> x) s e n
+        ys = mapLinear (\x -> f x) s e n
+
+-- | Plot a line between 0 and the length of the array with the given y values
+line1 :: (Foldable t, ToJSON (t a)) => t a -> Matplotlib
+line1 y = line [0..length y] y
+
+-- | Plot a matrix
+matShow :: ToJSON a => a -> Matplotlib
+matShow d = readData d
+            % (mp # "plot.matshow(data" ## ")")
+
+-- | Plot a matrix
+pcolor :: ToJSON a => a -> Matplotlib
+pcolor d = readData d
+            % (mp # "plot.pcolor(np.array(data)" ## ")")
+
+-- | Plot a KDE of the given functions; a good bandwith will be chosen automatically
+density :: [Double] -> Maybe (Double, Double) -> Matplotlib
+density l maybeStartEnd =
+  densityBandwidth l (((4 * (variance ** 5)) / (fromIntegral $ 3 * length l)) ** (1 / 5) / 3) maybeStartEnd
+  where mean = foldl' (+) 0 l / (fromIntegral $ length l)
+        variance = foldl' (+) 0 (map (\x -> sqr (x - mean)) l) / (fromIntegral $ length l)
+        sqr x = x * x
+
+-- * Matplotlib configuration
+
+setTeX :: Bool -> Matplotlib
+setTeX b = mp # "matplotlib.rcParams['text.usetex'] = " # b
+
+setUnicode :: Bool -> Matplotlib
+setUnicode b = mp # "matplotlib.rcParams['text.latex.unicode'] = " # b
+
+
+-- * Basic plotting commands
+
 -- | Plot the 'a' and 'b' entries of the data object
 dataPlot :: (MplotValue val, MplotValue val1) => val1 -> val -> Matplotlib
 dataPlot a b = mp # "p = plot.plot(data[" # a # "], data[" # b # "]" ## ")"
@@ -82,9 +184,9 @@
 plot :: (ToJSON t, ToJSON t1) => t1 -> t -> Matplotlib
 plot x y = readData (x, y) % dataPlot 0 1
 
--- | Show grid lines
-gridLines :: Matplotlib
-gridLines = mp # "ax.grid(True)"
+-- | Show/hide grid lines
+grid :: Bool -> Matplotlib
+grid t = mp # "plot.gca().grid(" # t # ")"
 
 -- | Plot x against y where x is a date.
 --   xunit is something like 'weeks', yearStart, monthStart, dayStart are an offset to x.
@@ -99,49 +201,28 @@
   
 -- | Add a label to the x axis
 xLabel :: MplotValue val => val -> Matplotlib
-xLabel label = mp # "plot.xlabel('" # label # "')"
+xLabel label = mp # "plot.xlabel(r'" # label # "'" ## ")"
 
 -- | Add a label to the y axis
 yLabel :: MplotValue val => val -> Matplotlib
-yLabel label = mp # "plot.ylabel('" # label # "')"
+yLabel label = mp # "plot.ylabel(r'" # label # "'" ## ")"
 
 -- | Add a label to the z axis
 zLabel :: MplotValue val => val -> Matplotlib
-zLabel label = mp # "plot.zlabel('" # label # "')"
+zLabel label = mp # "plot.zlabel(r'" # label # "'" ## ")"
 
 -- | Create a histogram for the 'a' entry of the data array
 dataHistogram :: (MplotValue val1, MplotValue val) => val1 -> val -> Matplotlib
 dataHistogram a bins = mp # "plot.hist(data[" # a # "]," # bins ## ")"
 
--- | Plot a histogram for the given values with 'bins'
-histogram :: (MplotValue val, ToJSON t) => t -> val -> Matplotlib
-histogram values bins = readData [values] % dataHistogram 0 bins
-
--- | Plot & show the histogram
-showHistogram :: (ToJSON t, MplotValue val) => t -> val -> IO (Either String String)
-showHistogram values bins = onscreen $ histogram values bins
-
 -- | Create a scatter plot accessing the given fields of the data array
 dataScatter :: (MplotValue val1, MplotValue val) => val1 -> val -> Matplotlib
 dataScatter a b = dataPlot a b `def` [o1 "'.'"]
 
--- | Plot the given values as a scatter plot
-scatter :: (ToJSON t1, ToJSON t) => t1 -> t -> Matplotlib
-scatter x y = plot x y `def` [o1 "'.'"]
-
 -- | Create a line accessing the given entires of the data array
 dataLine :: (MplotValue val1, MplotValue val) => val1 -> val -> Matplotlib
 dataLine a b = dataPlot a b `def` [o1 "'-'"]
 
--- | Plot a line
-line :: (ToJSON t1, ToJSON t) => t1 -> t -> Matplotlib
-line x y = plot x y `def` [o1 "'-'"]
-
--- | Plot a line given a function that will be executed for each element of
--- given list. The list provides the x values, the function the y values.
-lineF :: (ToJSON a, ToJSON b) => (a -> b) -> [a] -> Matplotlib
-lineF f l = plot l (map f l) `def` [o1 "'-'"]
-
 -- | Create a 3D contour
 contour xs ys zs =
   readData (xs, ys, zs)
@@ -157,25 +238,9 @@
   % contourRaw 0 1 2 (maximum2 xs) (maximum2 ys) (minimum2 zs)
   % axis3DLabels xs ys zs
 
--- | Given a grid of x and y values and a number of steps call the given
--- function and plot the 3D contour
-contourF :: (ToJSON val, MplotValue val, Ord val) => (Double -> Double -> val) -> Double -> Double -> Double -> Double -> Double -> Matplotlib
-contourF f xStart xEnd yStart yEnd steps = contour xs ys zs
-  where xs = mapLinear (\x -> (mapLinear (\_ -> x) yStart yEnd steps)) xStart xEnd steps
-        ys = mapLinear (\_ -> (mapLinear (\y -> y) yStart yEnd steps)) xStart xEnd steps
-        zs = mapLinear (\x -> (mapLinear (\y -> f x y) yStart yEnd steps)) xStart xEnd steps
-
--- | Given a grid of x and y values and a number of steps call the given
--- function and plot the 3D projection
-projectionsF :: (ToJSON val, MplotValue val, Ord val) => (Double -> Double -> val) -> Double -> Double -> Double -> Double -> Double -> Matplotlib
-projectionsF f xStart xEnd yStart yEnd steps = projections xs ys zs
-  where xs = mapLinear (\x -> (mapLinear (\_ -> x) yStart yEnd steps)) xStart xEnd steps
-        ys = mapLinear (\_ -> (mapLinear (\y -> y) yStart yEnd steps)) xStart xEnd steps
-        zs = mapLinear (\x -> (mapLinear (\y -> f x y) yStart yEnd steps)) xStart xEnd steps
-
 -- | Enable 3D projection
 axis3DProjection :: Matplotlib
-axis3DProjection = mp # "ax = plot.figure().gca(projection='3d')"
+axis3DProjection = mp # "ax = plot.gca(projection='3d')"
 
 -- | Plot a 3D wireframe accessing the given elements of the data array
 wireframe :: (MplotValue val2, MplotValue val1, MplotValue val) => val2 -> val1 -> val -> Matplotlib
@@ -217,10 +282,10 @@
   mp # "ax.bar(np.arange(len(data[" # a # "]))+" # offset # ", data[" # a # "], " # width ## ")" @@ opts
 
 -- | Create a subplot with the coordinates (r,c,f)
-addSubplot r c f = mp # "ax = plot.figure().add_subplot(" # r # c # f ## ")"
+addSubplot r c f = mp # "ax = plot.gcf().add_subplot(" # r # c # f ## ")"
 
 -- | Access a subplot with the coordinates (r,c,f)
-mplotSubplot r c f = mp # "ax = plot.subplot(" # r # "," # c # "," # f ## ")"
+getSubplot r c f = mp # "ax = plot.subplot(" # r # "," # c # "," # f ## ")"
 
 -- | The default bar with
 barDefaultWidth nr = 1.0 / (fromIntegral nr + 1)
@@ -241,7 +306,7 @@
 
 -- | Add a title
 title :: MplotValue val => val -> Matplotlib
-title s = mp # "plot.title('" # s ## "')"
+title s = mp # "plot.title(r'" # s # "'" ## ")"
 
 -- | Set the spacing of ticks on the x axis
 axisXTickSpacing :: (MplotValue val1, MplotValue val) => val1 -> val -> Matplotlib
@@ -257,13 +322,6 @@
   (mp # "data[" # b # "] = mlab.stineman_interp(np.linspace(data[" # a # "][0],data[" # a # "][-1]," # n # "),data[" # a # "],data[" # b # "],None)")
   % (mp # "data[" # a # "] = np.linspace(data[" # a # "][0],data[" # a # "][-1]," # n # ")")
 
--- | Plot x against y interpolating with n steps
-plotInterpolated :: (MplotValue val, ToJSON t, ToJSON t1) => t1 -> t -> val -> Matplotlib
-plotInterpolated x y n =
-  readData (x, y)
-  % interpolate 0 1 n
-  % dataPlot 0 1 `def` [o1 "-"]
-
 -- | Square up the aspect ratio of a plot.
 squareAxes :: Matplotlib
 squareAxes = mp # "plot.axes().set_aspect('equal')"
@@ -296,21 +354,6 @@
 ylim :: (MplotValue val1, MplotValue val) => val1 -> val -> Matplotlib
 ylim l u = mp # "plot.ylim([" # l # "," # u # "])"
 
--- | A handy function to plot a line between two points give a function and a number o steps
-plotMapLinear :: ToJSON b => (Double -> b) -> Double -> Double -> Double -> Matplotlib
-plotMapLinear f s e n = line xs ys
-  where xs = mapLinear (\x -> x) s e n
-        ys = mapLinear (\x -> f x) s e n
-
--- | Plot a line between 0 and the length of the array with the given y values
-line1 :: (Foldable t, ToJSON (t a)) => t a -> Matplotlib
-line1 y = line [0..length y] y
-
--- | Plot a matrix
-matShow :: ToJSON a => a -> Matplotlib
-matShow d = readData d
-            % (mp # "plot.matshow(data" ## ")")
-
 -- | Plot a KDE of the given functions with an optional start/end and a bandwidth h
 densityBandwidth :: [Double] -> Double -> Maybe (Double, Double) -> Matplotlib
 densityBandwidth l h maybeStartEnd =
@@ -325,10 +368,26 @@
         gaussianPdf x mu sigma = exp (- sqr (x - mu) / (2 * sigma)) / sqrt (2 * pi * sigma)
         sqr x = x * x
 
--- | Plot a KDE of the given functions; a good bandwith will be chosen automatically
-density :: [Double] -> Maybe (Double, Double) -> Matplotlib
-density l maybeStartEnd =
-  densityBandwidth l (((4 * (variance ** 5)) / (fromIntegral $ 3 * length l)) ** (1 / 5) / 3) maybeStartEnd
-  where mean = foldl' (+) 0 l / (fromIntegral $ length l)
-        variance = foldl' (+) 0 (map (\x -> sqr (x - mean)) l) / (fromIntegral $ length l)
-        sqr x = x * x
+-- | Add a horizontal line across the axis
+axhline y = mp # "ax.axhline(" # y ## ")"
+
+-- | Plot cross-correlation
+xcorr x y = readData (x, y) % mp # "ax.xcorr(data[0], data[1]" ## ")"
+
+-- | Plot auto-correlation
+acorr x = readData x % mp # "ax.acorr(data" ## ")"
+
+-- | Plot text at a specified location
+text x y str = mp # "ax.text(" # x # "," # y # "," # "'" # str # "'" ## ")"
+
+-- | Insert a legend
+legend = mp # "plot.legend(" ## ")"
+
+-- | Insert a color bar
+colorbar = mp # "plot.colorbar(" ## ")"
+
+-- | Creates subplots and stores them in an internal variable
+subplots = mp # "fig, axes = plot.subplots(" ## ")"
+
+-- | Access a subplot
+setSubplot s = mp # "ax = axes[" # s # "]"
diff --git a/src/Graphics/Matplotlib/Internal.hs b/src/Graphics/Matplotlib/Internal.hs
--- a/src/Graphics/Matplotlib/Internal.hs
+++ b/src/Graphics/Matplotlib/Internal.hs
@@ -120,6 +120,8 @@
   toPython s = show s
 instance MplotValue Int where
   toPython s = show s
+instance MplotValue Bool where
+  toPython s = show s
 instance (MplotValue x) => MplotValue (x, x) where
   toPython (n, v) = toPython n ++ " = " ++ toPython v
 instance (MplotValue (x, y)) => MplotValue [(x, y)] where
@@ -134,7 +136,7 @@
 -- optFn :: Matplotlib -> Matplotlib
 optFn :: ([Option] -> String) -> Matplotlib -> Matplotlib
 optFn f l | isJust $ mpPendingOption l = error "Commands can have only open option. TODO Enforce this through the type system or relax it!"
-          | otherwise = l' { mpPendingOption = Just (\os -> Exec (sl ++ f os)) }
+          | otherwise = l' { mpPendingOption = Just (\os -> Exec (sl `combine` f os)) }
   where (l', (Exec sl)) = removeLast l
         removeLast x@(Matplotlib _ Nothing s) = (x { mpRest = sdeleteAt (S.length s - 1) s }
                                                 , fromMaybe (Exec "") (slookup (S.length s - 1) s))
@@ -144,6 +146,10 @@
                     | otherwise      = Nothing
         sdeleteAt i s | i < S.length s = S.take i s >< S.drop (i + 1) s
                       | otherwise      = s
+        combine [] r = r
+        combine l [] = l
+        combine l r | [last l] == "(" && [head r] == "," = l ++ tail r
+                    | otherwise = l ++ r
 
 -- | Merge two commands with options between
 options :: Matplotlib -> Matplotlib
@@ -218,6 +224,8 @@
              ,"import matplotlib.patches as mpatches"
              ,"import matplotlib.pyplot as plot"
              ,"import matplotlib.mlab as mlab"
+             ,"import matplotlib.colors as mcolors"
+             ,"import matplotlib.collections as mcollections"
              ,"from matplotlib import cm"
              ,"from mpl_toolkits.mplot3d import axes3d"
              ,"import numpy as np"
@@ -225,7 +233,9 @@
              ,"import sys"
              ,"import json"
              ,"import random, datetime"
-             ,"from matplotlib.dates import DateFormatter, WeekdayLocator"]
+             ,"from matplotlib.dates import DateFormatter, WeekdayLocator"
+             ,"ax = plot.figure().gca()"
+             ,"axes = [plot.figure().gca()]"]
 
 -- | The python command that reads external data into the python data array
 pyReadData :: [Char] -> [[Char]]
diff --git a/test/Spec.hs b/test/Spec.hs
--- a/test/Spec.hs
+++ b/test/Spec.hs
@@ -1,4 +1,4 @@
-{-# language ExtendedDefaultRules #-}
+{-# language ExtendedDefaultRules, ScopedTypeVariables, QuasiQuotes #-}
 
 import Test.Tasty
 import Test.Tasty.HUnit
@@ -6,37 +6,17 @@
 import Graphics.Matplotlib
 import System.IO.Temp
 import System.Random
-
-main = defaultMain tests
-
-tests :: TestTree
-tests = testGroup "Tests" [unitTests]
+import Text.RawString.QQ
+import Data.List
+import Data.List.Split
 
--- | 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" (\file _ -> figure file fn)
+-- * Random values for testing
 
--- | This generates examples from the test cases
-testPlot' name fn = testCase name $ tryit fn name @?= Right ""
-  where tryit fn name = unsafePerformIO $ figure ("/tmp/imgs/" ++ name ++ ".png") fn
+uniforms :: (Random a, Num a) => [a]
+uniforms = randoms (mkStdGen 42)
+uniforms' lo hi = randomRs (lo,hi) (mkStdGen 42)
 
-unitTests = testGroup "Unit tests"
-  [ testPlot "histogram" m1
-  , testPlot "cumulative" m2
-  , testPlot "scatter" m3
-  , testPlot "contour" m4
-  , testPlot "labelled-histogram" m5
-  -- TODO This test case is broken
-  -- , testPlot "sub-bars" $ tryit m6 "m6" @?= Right ""
-  , testPlot "density-bandwidth" m7
-  , testPlot "density" m8
-  , testPlot "line-function" m9
-  , testPlot "quadratic" m10
-  , testPlot "projections" m11
-  , testPlot "line-options" m12
-  , testPlot "xcorr" mxcorr
-  ]
+-- * 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,
@@ -79,6 +59,77 @@
       -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
+
+-- * Tests
+
+main = defaultMain tests
+
+tests :: TestTree
+tests = testGroup "Tests" [unitTests]
+
+-- | 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" (\file _ -> figure file fn)
+
+-- | This generates examples from the test cases
+testPlot' name fn = testCase name $ tryit fn name @?= Right ""
+  where tryit fn name = unsafePerformIO $ do
+          c <- code fn
+          print c
+          figure ("/tmp/imgs/" ++ name ++ ".png") fn
+
+unitTests = testGroup "Unit tests"
+  [ testPlot "histogram" m1
+  , testPlot "cumulative" m2
+  , testPlot "scatter" m3
+  , testPlot "contour" m4
+  , testPlot "labelled-histogram" m5
+  -- TODO This test case is broken
+  -- , testPlot "sub-bars" m6
+  , testPlot "density-bandwidth" m7
+  , testPlot "density" m8
+  , testPlot "line-function" m9
+  , testPlot "quadratic" m10
+  , testPlot "projections" m11
+  , testPlot "line-options" m12
+  , testPlot "corr" mxcorr
+  , testPlot "tex" mtex
+  , testPlot "show-matrix" mmat
+  , testPlot "legend" mlegend
+  , testPlot "hist2DLog" mhist2DLog
+  , testPlot "eventplot" meventplot
+  , testPlot "errorbar" merrorbar
+  , testPlot "boxplot" mboxplot
+  ]
+
+-- * These tests are fully-reproducible, the output must be identical every time
+
 m1 :: Matplotlib
 m1 = histogram xs 8
 
@@ -116,10 +167,49 @@
 
 m12 = plot [1,2,3,4,5,6] [1,3,2,5,2,4] @@ [o1 "'go-'", o2 "linewidth" "2"]
 
-datas = randoms (mkStdGen 42)
+-- * These tests can be random and may not be exactly the same every time
 
-datas lo hi = randomRs (lo,hi) (mkStdGen 42)
+-- | 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
 
-mxcorr = 
-  where (xs :: Double) = take 100 datas
-        (ys :: Double) = take 100 datas
+-- | http://matplotlib.org/examples/pylab_examples/tex_unicode_demo.html
+mtex = plotMapLinear cos 0 1 100
+  % setTeX True
+  % setUnicode True
+  % xLabel [r|\textbf{time (s)}|]
+  % yLabel [r|\textit{Velocity (\u00B0/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" "mcolors.LogNorm()"]
+  % 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 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
+
+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"]
+