diff --git a/third_party/matplotlib-cpp/.gitignore b/third_party/matplotlib-cpp/.gitignore
index 7622be7..1c4a1b0 100644
--- a/third_party/matplotlib-cpp/.gitignore
+++ b/third_party/matplotlib-cpp/.gitignore
@@ -33,3 +33,6 @@
 
 # Build
 /examples/build/*
+
+# vim temp files
+*.sw*
diff --git a/third_party/matplotlib-cpp/.travis.yml b/third_party/matplotlib-cpp/.travis.yml
new file mode 100644
index 0000000..d6175a2
--- /dev/null
+++ b/third_party/matplotlib-cpp/.travis.yml
@@ -0,0 +1,6 @@
+language: minimal
+dist: trusty
+services:
+  - docker
+script:
+  - make -C contrib docker_build
diff --git a/third_party/matplotlib-cpp/BUILD b/third_party/matplotlib-cpp/BUILD
index 85249c0..95f7192 100644
--- a/third_party/matplotlib-cpp/BUILD
+++ b/third_party/matplotlib-cpp/BUILD
@@ -5,10 +5,14 @@
     hdrs = [
         "matplotlibcpp.h",
     ],
+    data = [
+        "@matplotlib_repo//:matplotlib3",
+        "@python_repo//:all_files",
+    ],
     restricted_to = ["//tools:k8"],
     visibility = ["//visibility:public"],
     deps = [
-        "@python_repo//:python2.7_lib",
+        "@python_repo//:python3.5_lib",
     ],
 )
 
diff --git a/third_party/matplotlib-cpp/Makefile b/third_party/matplotlib-cpp/Makefile
index e75d134..67b5ac3 100644
--- a/third_party/matplotlib-cpp/Makefile
+++ b/third_party/matplotlib-cpp/Makefile
@@ -1,22 +1,41 @@
-examples: minimal basic modern animation nonblock xkcd
+# Use C++11, dont warn on long-to-float conversion
+CXXFLAGS += -std=c++11 -Wno-conversion
 
-minimal: examples/minimal.cpp matplotlibcpp.h
-	cd examples && g++ -DWITHOUT_NUMPY minimal.cpp -I/usr/include/python2.7 -lpython2.7 -o minimal -std=c++11
+# Default to using system's default version of python
+PYTHON_BIN     ?= python3
+PYTHON_CONFIG  := $(PYTHON_BIN)-config
+PYTHON_INCLUDE ?= $(shell $(PYTHON_CONFIG) --includes)
+EXTRA_FLAGS    := $(PYTHON_INCLUDE)
+# NOTE: Since python3.8, the correct invocation is `python3-config --libs --embed`.
+# So of course the proper way to get python libs for embedding now is to
+# invoke that, check if it crashes, and fall back to just `--libs` if it does.
+LDFLAGS        += $(shell if $(PYTHON_CONFIG) --ldflags --embed >/dev/null; then $(PYTHON_CONFIG) --ldflags --embed; else $(PYTHON_CONFIG) --ldflags; fi)
 
-basic: examples/basic.cpp matplotlibcpp.h
-	cd examples && g++ basic.cpp -I/usr/include/python2.7 -lpython2.7 -o basic
+# Either finds numpy or set -DWITHOUT_NUMPY
+EXTRA_FLAGS     += $(shell $(PYTHON_BIN) $(CURDIR)/numpy_flags.py)
+WITHOUT_NUMPY   := $(findstring $(EXTRA_FLAGS), WITHOUT_NUMPY)
 
-modern: examples/modern.cpp matplotlibcpp.h
-	cd examples && g++ modern.cpp -I/usr/include/python2.7 -lpython2.7 -o modern -std=c++11
+# Examples requiring numpy support to compile
+EXAMPLES_NUMPY  := surface colorbar
+EXAMPLES        := minimal basic modern animation nonblock xkcd quiver bar \
+	           fill_inbetween fill update subplot2grid lines3d \
+                   $(if $(WITHOUT_NUMPY),,$(EXAMPLES_NUMPY))
 
-animation: examples/animation.cpp matplotlibcpp.h
-	cd examples && g++ animation.cpp -I/usr/include/python2.7 -lpython2.7 -o animation -std=c++11
+# Prefix every example with 'examples/build/'
+EXAMPLE_TARGETS := $(patsubst %,examples/build/%,$(EXAMPLES))
 
-nonblock: examples/nonblock.cpp matplotlibcpp.h
-	cd examples && g++ nonblock.cpp -I/usr/include/python2.7 -lpython2.7 -o nonblock -std=c++11
+.PHONY: examples
 
-xkcd: examples/xkcd.cpp matplotlibcpp.h
-	cd examples && g++ xkcd.cpp -I/usr/include/python2.7 -lpython2.7 -o xkcd -std=c++11
+examples: $(EXAMPLE_TARGETS)
+
+docs:
+	doxygen
+	moxygen doc/xml --noindex -o doc/api.md
+
+# Assume every *.cpp file is a separate example
+$(EXAMPLE_TARGETS): examples/build/%: examples/%.cpp matplotlibcpp.h
+	mkdir -p examples/build
+	$(CXX) -o $@ $< $(EXTRA_FLAGS) $(CXXFLAGS) $(LDFLAGS)
 
 clean:
-	rm -f examples/{minimal,basic,modern,animation,nonblock,xkcd}
+	rm -f ${EXAMPLE_TARGETS}
diff --git a/third_party/matplotlib-cpp/README.md b/third_party/matplotlib-cpp/README.md
index e0d74ba..61ffef4 100644
--- a/third_party/matplotlib-cpp/README.md
+++ b/third_party/matplotlib-cpp/README.md
@@ -30,7 +30,7 @@
 
 namespace plt = matplotlibcpp;
 
-int main() 
+int main()
 {
     // Prepare data.
     int n = 5000;
@@ -41,15 +41,18 @@
         z.at(i) = log(i);
     }
 
+    // Set the size of output image to 1200x780 pixels
+    plt::figure_size(1200, 780);
     // Plot line from given x and y data. Color is selected automatically.
     plt::plot(x, y);
     // Plot a red dashed line from given x and y data.
     plt::plot(x, w,"r--");
     // Plot a line whose name will show up as "log(x)" in the legend.
     plt::named_plot("log(x)", x, z);
-
     // Set x-axis to interval [0,1000000]
     plt::xlim(0, 1000*1000);
+    // Add graph title
+    plt::title("Sample figure");
     // Enable legend.
     plt::legend();
     // Save the image (file format is determined by the extension)
@@ -58,9 +61,11 @@
 ```
     g++ basic.cpp -I/usr/include/python2.7 -lpython2.7
 
-Result: ![Basic example](./examples/basic.png)
+**Result:**
 
-matplotlib-cpp doesn't require C++11, but will enable some additional syntactic sugar when available:
+![Basic example](./examples/basic.png)
+
+Alternatively, matplotlib-cpp also supports some C++11-powered syntactic sugar:
 ```cpp
 #include <cmath>
 #include "matplotlibcpp.h"
@@ -68,30 +73,32 @@
 using namespace std;
 namespace plt = matplotlibcpp;
 
-int main() 
-{    
+int main()
+{
     // Prepare data.
     int n = 5000; // number of data points
-    vector<double> x(n),y(n); 
+    vector<double> x(n),y(n);
     for(int i=0; i<n; ++i) {
         double t = 2*M_PI*i/n;
         x.at(i) = 16*sin(t)*sin(t)*sin(t);
         y.at(i) = 13*cos(t) - 5*cos(2*t) - 2*cos(3*t) - cos(4*t);
     }
 
-    // plot() takes an arbitrary number of (x,y,format)-triples. 
+    // plot() takes an arbitrary number of (x,y,format)-triples.
     // x must be iterable (that is, anything providing begin(x) and end(x)),
-    // y must either be callable (providing operator() const) or iterable. 
+    // y must either be callable (providing operator() const) or iterable.
     plt::plot(x, y, "r-", x, [](double d) { return 12.5+abs(sin(d)); }, "k-");
 
 
     // show plots
     plt::show();
-} 
+}
 ```
     g++ modern.cpp -std=c++11 -I/usr/include/python2.7 -lpython
 
-Result: ![Modern example](./examples/modern.png)
+**Result:**
+
+![Modern example](./examples/modern.png)
 
 Or some *funny-looking xkcd-styled* example:
 ```cpp
@@ -121,7 +128,66 @@
 
 **Result:**
 
-![Minimal example](./examples/xkcd.png)
+![xkcd example](./examples/xkcd.png)
+
+When working with vector fields, you might be interested in quiver plots:
+```cpp
+#include "../matplotlibcpp.h"
+
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    // u and v are respectively the x and y components of the arrows we're plotting
+    std::vector<int> x, y, u, v;
+    for (int i = -5; i <= 5; i++) {
+        for (int j = -5; j <= 5; j++) {
+            x.push_back(i);
+            u.push_back(-i);
+            y.push_back(j);
+            v.push_back(-j);
+        }
+    }
+
+    plt::quiver(x, y, u, v);
+    plt::show();
+}
+```
+    g++ quiver.cpp -std=c++11 -I/usr/include/python2.7 -lpython2.7
+
+**Result:**
+
+![quiver example](./examples/quiver.png)
+
+When working with 3d functions, you might be interested in 3d plots:
+```cpp
+#include "../matplotlibcpp.h"
+
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    std::vector<std::vector<double>> x, y, z;
+    for (double i = -5; i <= 5;  i += 0.25) {
+        std::vector<double> x_row, y_row, z_row;
+        for (double j = -5; j <= 5; j += 0.25) {
+            x_row.push_back(i);
+            y_row.push_back(j);
+            z_row.push_back(::std::sin(::std::hypot(i, j)));
+        }
+        x.push_back(x_row);
+        y.push_back(y_row);
+        z.push_back(z_row);
+    }
+
+    plt::plot_surface(x, y, z);
+    plt::show();
+}
+```
+
+**Result:**
+
+![surface example](./examples/surface.png)
 
 Installation
 ------------
@@ -133,23 +199,65 @@
     sudo apt-get install python-matplotlib python-numpy python2.7-dev
 
 If, for some reason, you're unable to get a working installation of numpy on your system,
-you can add the define `WITHOUT_NUMPY` to erase this dependency.
+you can define the macro `WITHOUT_NUMPY` before including the header file to erase this
+dependency.
 
 The C++-part of the library consists of the single header file `matplotlibcpp.h` which can be placed
 anywhere.
 
-Since a python interpreter is opened internally, it is necessary to link against `libpython2.7` in order to use
-matplotlib-cpp.
+Since a python interpreter is opened internally, it is necessary to link against `libpython` in order
+to user matplotlib-cpp. Most versions should work, although `libpython2.7` and `libpython3.6` are
+probably the most regularly testedr.
+
 
 # CMake
 
 If you prefer to use CMake as build system, you will want to add something like this to your
 CMakeLists.txt:
+
+**Recommended way (since CMake 3.12):**
+
+It's easy to use cmake official [docs](https://cmake.org/cmake/help/git-stage/module/FindPython2.html#module:FindPython2) to find Python 2(or 3) interpreter, compiler and development environment (include directories and libraries).
+
+NumPy is optional here, delete it from cmake script, if you don't need it.
+
+```cmake
+find_package(Python2 COMPONENTS Development NumPy)
+target_include_directories(myproject PRIVATE ${Python2_INCLUDE_DIRS} ${Python2_NumPy_INCLUDE_DIRS})
+target_link_libraries(myproject Python2::Python Python2::NumPy)
+```
+
+**Alternative way (for CMake <= 3.11):**
+
 ```cmake
 find_package(PythonLibs 2.7)
 target_include_directories(myproject PRIVATE ${PYTHON_INCLUDE_DIRS})
 target_link_libraries(myproject ${PYTHON_LIBRARIES})
 ```
+
+
+# Vcpkg
+
+You can download and install matplotlib-cpp using the [vcpkg](https://github.com/Microsoft/vcpkg) dependency manager:
+
+    git clone https://github.com/Microsoft/vcpkg.git
+    cd vcpkg
+    ./bootstrap-vcpkg.sh
+    ./vcpkg integrate install
+    vcpkg install matplotlib-cpp
+  
+The matplotlib-cpp port in vcpkg is kept up to date by Microsoft team members and community contributors. If the version is out of date, please [create an issue or pull request](https://github.com/Microsoft/vcpkg) on the vcpkg repository.
+
+
+# C++11
+
+Currently, c++11 is required to build matplotlib-cpp. The last working commit that did
+not have this requirement was `717e98e752260245407c5329846f5d62605eff08`.
+
+Note that support for c++98 was dropped more or less accidentally, so if you have to work
+with an ancient compiler and still want to enjoy the latest additional features, I'd
+probably merge a PR that restores support.
+
 # Python 3
 
 This library supports both python2 and python3 (although the python3 support is probably far less tested,
@@ -165,15 +273,15 @@
 
 Why?
 ----
-I initially started this library during my diploma thesis. The usual approach of 
+I initially started this library during my diploma thesis. The usual approach of
 writing data from the c++ algorithm to a file and afterwards parsing and plotting
 it in python using matplotlib proved insufficient: Keeping the algorithm
-and plotting code in sync requires a lot of effort when the C++ code frequently and substantially 
+and plotting code in sync requires a lot of effort when the C++ code frequently and substantially
 changes. Additionally, the python yaml parser was not able to cope with files that
 exceed a few hundred megabytes in size.
 
 Therefore, I was looking for a C++ plotting library that was extremely easy to use
-and to add into an existing codebase, preferrably header-only. When I found
+and to add into an existing codebase, preferably header-only. When I found
 none, I decided to write one myself, which is basically a C++ wrapper around
 matplotlib. As you can see from the above examples, plotting data and saving it
 to an image file can be done as few as two lines of code.
@@ -197,3 +305,6 @@
 
 * If you use Anaconda on Windows, you might need to set PYTHONHOME to Anaconda home directory and QT_QPA_PLATFORM_PLUGIN_PATH to %PYTHONHOME%Library/plugins/platforms. The latter is for especially when you get the error which says 'This application failed to start because it could not find or load the Qt platform plugin "windows"
 in "".'
+
+* MacOS: `Unable to import matplotlib.pyplot`. Cause: In mac os image rendering back end of matplotlib (what-is-a-backend to render using the API of Cocoa by default). There is Qt4Agg and GTKAgg and as a back-end is not the default. Set the back end of macosx that is differ compare with other windows or linux os.
+Solution is described [here](https://stackoverflow.com/questions/21784641/installation-issue-with-matplotlib-python?noredirect=1&lq=1), additional information can be found there too(see links in answers).
diff --git a/third_party/matplotlib-cpp/contrib/CMakeLists.txt b/third_party/matplotlib-cpp/contrib/CMakeLists.txt
index ba14b86..edb40b1 100644
--- a/third_party/matplotlib-cpp/contrib/CMakeLists.txt
+++ b/third_party/matplotlib-cpp/contrib/CMakeLists.txt
@@ -1,10 +1,11 @@
-cmake_minimum_required(VERSION 3.1)
+cmake_minimum_required(VERSION 3.7)
 project (MatplotlibCPP_Test)
 
 set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 include_directories(${PYTHONHOME}/include)
+include_directories(${PYTHONHOME}/Lib/site-packages/numpy/core/include)
 link_directories(${PYTHONHOME}/libs)
 
 add_definitions(-DMATPLOTLIBCPP_PYTHON_HEADER=Python.h)
@@ -12,10 +13,14 @@
 # message(STATUS "*** dump start cmake variables ***")
 # get_cmake_property(_variableNames VARIABLES)
 # foreach(_variableName ${_variableNames})
-        # message(STATUS "${_variableName}=${${_variableName}}")
+#         message(STATUS "${_variableName}=${${_variableName}}")
 # endforeach()
 # message(STATUS "*** dump end ***")
 
 add_executable(minimal ${CMAKE_CURRENT_SOURCE_DIR}/../examples/minimal.cpp)
 add_executable(basic ${CMAKE_CURRENT_SOURCE_DIR}/../examples/basic.cpp)
 add_executable(modern ${CMAKE_CURRENT_SOURCE_DIR}/../examples/modern.cpp)
+add_executable(animation ${CMAKE_CURRENT_SOURCE_DIR}/../examples/animation.cpp)
+add_executable(nonblock ${CMAKE_CURRENT_SOURCE_DIR}/../examples/nonblock.cpp)
+add_executable(xkcd ${CMAKE_CURRENT_SOURCE_DIR}/../examples/xkcd.cpp)
+add_executable(bar ${CMAKE_CURRENT_SOURCE_DIR}/../examples/bar.cpp)
diff --git a/third_party/matplotlib-cpp/contrib/Dockerfile b/third_party/matplotlib-cpp/contrib/Dockerfile
new file mode 100644
index 0000000..850466f
--- /dev/null
+++ b/third_party/matplotlib-cpp/contrib/Dockerfile
@@ -0,0 +1,27 @@
+FROM debian:10 AS builder
+RUN apt-get update \
+ && apt-get install --yes --no-install-recommends \
+    g++ \
+    libpython3-dev \
+    make \
+    python3 \
+    python3-dev \
+    python3-numpy
+
+ADD Makefile matplotlibcpp.h numpy_flags.py /opt/
+ADD examples/*.cpp /opt/examples/
+RUN cd /opt \
+ && make PYTHON_BIN=python3 \
+ && ls examples/build
+
+FROM debian:10
+RUN apt-get update \
+ && apt-get install --yes --no-install-recommends \
+    libpython3-dev \
+    python3-matplotlib \
+    python3-numpy
+
+COPY --from=builder /opt/examples/build /opt/
+RUN cd /opt \
+ && ls \
+ && ./basic
diff --git a/third_party/matplotlib-cpp/contrib/Makefile b/third_party/matplotlib-cpp/contrib/Makefile
new file mode 100644
index 0000000..f659cd9
--- /dev/null
+++ b/third_party/matplotlib-cpp/contrib/Makefile
@@ -0,0 +1,6 @@
+all: docker_build
+
+docker_build:
+	cd .. && \
+	docker build . -f contrib/Dockerfile -t matplotlibcpp && \
+	cd contrib
diff --git a/third_party/matplotlib-cpp/contrib/README.md b/third_party/matplotlib-cpp/contrib/README.md
index efc0a50..0af8515 100644
--- a/third_party/matplotlib-cpp/contrib/README.md
+++ b/third_party/matplotlib-cpp/contrib/README.md
@@ -8,12 +8,25 @@
 changes break any of them.
 
 ## Windows support
+Tested on the following environment
+* Windows 10 - 64bit
+* Anaconda 4.3 (64 bit)
+* Python 3.6.0
+* CMake 3.9.4
+* Visual Studio 2017, 2015, 2013
 
 ### Configuring and Building Samples
+1. Edit WinBuild.cmd for your environment(Line:5-7)
+    if NOT DEFINED MSVC_VERSION set MSVC_VERSION=[Your Visual Studio Version(12, 14, 15)]
+    if NOT DEFINED CMAKE_CONFIG set CMAKE_CONFIG=Release
+    if NOT DEFINED PYTHONHOME   set PYTHONHOME=[Your Python Path]
 
+2. Run WinBuild.cmd to build
 ```cmd
 > cd contrib
 > WinBuild.cmd
 ```
-
 The `WinBuild.cmd` will set up temporal ENV variables and build binaries in (matplotlib root)/examples with the Release configuration.
+
+3. Find exe files in examples/build/Release
+Note: platforms folder is necessary to make qt works.
diff --git a/third_party/matplotlib-cpp/contrib/WinBuild.cmd b/third_party/matplotlib-cpp/contrib/WinBuild.cmd
index 4e3b450..9dfd627 100644
--- a/third_party/matplotlib-cpp/contrib/WinBuild.cmd
+++ b/third_party/matplotlib-cpp/contrib/WinBuild.cmd
@@ -1,9 +1,14 @@
 @echo off
 @setlocal EnableDelayedExpansion
 
-if NOT DEFINED MSVC_VERSION set MSVC_VERSION=14
+REM ------Set Your Environment------------------------------- 
+if NOT DEFINED MSVC_VERSION set MSVC_VERSION=15
 if NOT DEFINED CMAKE_CONFIG set CMAKE_CONFIG=Release
 if NOT DEFINED PYTHONHOME   set PYTHONHOME=C:/Users/%username%/Anaconda3
+REM ---------------------------------------------------------
+
+set KEY_NAME="HKEY_LOCAL_MACHINE\SOFTWARE\WOW6432Node\Microsoft\VisualStudio\SxS\VS7"
+set VALUE_NAME=15.0
 
 if "%MSVC_VERSION%"=="14" (
     if "%processor_architecture%" == "AMD64" (
@@ -14,13 +19,23 @@
 ) else if "%MSVC_VERSION%"=="12" (
     if "%processor_architecture%" == "AMD64" (
         set CMAKE_GENERATOR=Visual Studio 12 2013 Win64
-    
     ) else (
         set CMAKE_GENERATOR=Visual Studio 12 2013
     )
+) else if "%MSVC_VERSION%"=="15" (
+    if "%processor_architecture%" == "AMD64" (
+        set CMAKE_GENERATOR=Visual Studio 15 2017 Win64
+    ) else (
+        set CMAKE_GENERATOR=Visual Studio 15 2017
+    )
 )
-
-set batch_file=!VS%MSVC_VERSION%0COMNTOOLS!..\..\VC\vcvarsall.bat
+if "%MSVC_VERSION%"=="15" (
+    for /F "usebackq tokens=1,2,*" %%A in (`REG QUERY %KEY_NAME% /v %VALUE_NAME%`) do (
+        set batch_file=%%CVC\Auxiliary\Build\vcvarsall.bat
+    )
+) else (
+    set batch_file=!VS%MSVC_VERSION%0COMNTOOLS!..\..\VC\vcvarsall.bat
+)
 call "%batch_file%" %processor_architecture%
 
 pushd ..
diff --git a/third_party/matplotlib-cpp/examples/.gitignore b/third_party/matplotlib-cpp/examples/.gitignore
new file mode 100644
index 0000000..3da8ad6
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/.gitignore
@@ -0,0 +1,14 @@
+animation
+bar
+basic
+fill
+fill_inbetween
+imshow
+minimal
+modern
+nonblock
+quiver
+subplot
+surface
+update
+xkcd
diff --git a/third_party/matplotlib-cpp/examples/bar.cpp b/third_party/matplotlib-cpp/examples/bar.cpp
new file mode 100644
index 0000000..86423ad
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/bar.cpp
@@ -0,0 +1,18 @@
+#define _USE_MATH_DEFINES
+
+#include <iostream>
+#include <string>
+#include "../matplotlibcpp.h"
+namespace plt = matplotlibcpp;
+
+int main(int argc, char **argv) {
+    std::vector<int> test_data;
+    for (int i = 0; i < 20; i++) {
+        test_data.push_back(i);
+    }
+
+    plt::bar(test_data);
+    plt::show();
+
+    return (0);
+}
diff --git a/third_party/matplotlib-cpp/examples/bar.png b/third_party/matplotlib-cpp/examples/bar.png
new file mode 100644
index 0000000..be6af0f
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/bar.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/basic.cpp b/third_party/matplotlib-cpp/examples/basic.cpp
index 2b02f27..2dc34c7 100644
--- a/third_party/matplotlib-cpp/examples/basic.cpp
+++ b/third_party/matplotlib-cpp/examples/basic.cpp
@@ -7,31 +7,38 @@
 
 int main() 
 {
-	// Prepare data.
-	int n = 5000;
-	std::vector<double> x(n), y(n), z(n), w(n,2);
-	for(int i=0; i<n; ++i) {
-		x.at(i) = i*i;
-		y.at(i) = sin(2*M_PI*i/360.0);
-		z.at(i) = log(i);
-	}
+    // Prepare data.
+    int n = 5000;
+    std::vector<double> x(n), y(n), z(n), w(n,2);
+    for(int i=0; i<n; ++i) {
+        x.at(i) = i*i;
+        y.at(i) = sin(2*M_PI*i/360.0);
+        z.at(i) = log(i);
+    }
+    
+    // Set the size of output image = 1200x780 pixels
+    plt::figure_size(1200, 780);
 
-	// Plot line from given x and y data. Color is selected automatically.
-	plt::plot(x, y);
-	// Plot a red dashed line from given x and y data.
-	plt::plot(x, w,"r--");
-	// Plot a line whose name will show up as "log(x)" in the legend.
-	plt::named_plot("log(x)", x, z);
+    // Plot line from given x and y data. Color is selected automatically.
+    plt::plot(x, y);
 
-	// Set x-axis to interval [0,1000000]
-	plt::xlim(0, 1000*1000);
+    // Plot a red dashed line from given x and y data.
+    plt::plot(x, w,"r--");
 
-	// Add graph title
-	plt::title("Sample figure");
-	// Enable legend.
-	plt::legend();
-	// save figure
-	const char* filename = "./basic.png";
-	std::cout << "Saving result to " << filename << std::endl;;
-	plt::save(filename);
+    // Plot a line whose name will show up as "log(x)" in the legend.
+    plt::named_plot("log(x)", x, z);
+
+    // Set x-axis to interval [0,1000000]
+    plt::xlim(0, 1000*1000);
+
+    // Add graph title
+    plt::title("Sample figure");
+
+    // Enable legend.
+    plt::legend();
+
+    // save figure
+    const char* filename = "./basic.png";
+    std::cout << "Saving result to " << filename << std::endl;;
+    plt::save(filename);
 }
diff --git a/third_party/matplotlib-cpp/examples/basic.png b/third_party/matplotlib-cpp/examples/basic.png
index eef82a2..4d87ff0 100644
--- a/third_party/matplotlib-cpp/examples/basic.png
+++ b/third_party/matplotlib-cpp/examples/basic.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/colorbar.cpp b/third_party/matplotlib-cpp/examples/colorbar.cpp
new file mode 100644
index 0000000..f53e01d
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/colorbar.cpp
@@ -0,0 +1,32 @@
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <iostream>
+#include "../matplotlibcpp.h"
+
+using namespace std;
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    // Prepare data
+    int ncols = 500, nrows = 300;
+    std::vector<float> z(ncols * nrows);
+    for (int j=0; j<nrows; ++j) {
+        for (int i=0; i<ncols; ++i) {
+            z.at(ncols * j + i) = std::sin(std::hypot(i - ncols/2, j - nrows/2));
+        }
+    }
+
+    const float* zptr = &(z[0]);
+    const int colors = 1;
+
+    plt::title("My matrix");
+    PyObject* mat;
+    plt::imshow(zptr, nrows, ncols, colors, {}, &mat);
+    plt::colorbar(mat);
+
+    // Show plots
+    plt::show();
+    plt::close();
+    Py_DECREF(mat);
+}
diff --git a/third_party/matplotlib-cpp/examples/fill.cpp b/third_party/matplotlib-cpp/examples/fill.cpp
new file mode 100644
index 0000000..6059b47
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/fill.cpp
@@ -0,0 +1,35 @@
+#define _USE_MATH_DEFINES
+#include "../matplotlibcpp.h"
+#include <cmath>
+
+using namespace std;
+namespace plt = matplotlibcpp;
+
+// Example fill plot taken from:
+// https://matplotlib.org/gallery/misc/fill_spiral.html
+int main() {
+    // Prepare data.
+    vector<double> theta;
+    for (double d = 0; d < 8 * M_PI; d += 0.1)
+        theta.push_back(d);
+
+    const int a = 1;
+    const double b = 0.2;
+
+    for (double dt = 0; dt < 2 * M_PI; dt += M_PI/2.0) {
+        vector<double> x1, y1, x2, y2;
+        for (double th : theta) {
+            x1.push_back( a*cos(th + dt) * exp(b*th) );
+            y1.push_back( a*sin(th + dt) * exp(b*th) );
+
+            x2.push_back( a*cos(th + dt + M_PI/4.0) * exp(b*th) );
+            y2.push_back( a*sin(th + dt + M_PI/4.0) * exp(b*th) );
+        }
+
+        x1.insert(x1.end(), x2.rbegin(), x2.rend());
+        y1.insert(y1.end(), y2.rbegin(), y2.rend());
+
+        plt::fill(x1, y1, {});
+    }
+    plt::show();
+}
diff --git a/third_party/matplotlib-cpp/examples/fill.png b/third_party/matplotlib-cpp/examples/fill.png
new file mode 100644
index 0000000..aa1fc0d
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/fill.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/fill_between.png b/third_party/matplotlib-cpp/examples/fill_between.png
new file mode 100644
index 0000000..a199423
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/fill_between.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/fill_inbetween.cpp b/third_party/matplotlib-cpp/examples/fill_inbetween.cpp
new file mode 100644
index 0000000..788d008
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/fill_inbetween.cpp
@@ -0,0 +1,28 @@
+#define _USE_MATH_DEFINES
+#include "../matplotlibcpp.h"
+#include <cmath>
+#include <iostream>
+
+using namespace std;
+namespace plt = matplotlibcpp;
+
+int main() {
+  // Prepare data.
+  int n = 5000;
+  std::vector<double> x(n), y(n), z(n), w(n, 2);
+  for (int i = 0; i < n; ++i) {
+    x.at(i) = i * i;
+    y.at(i) = sin(2 * M_PI * i / 360.0);
+    z.at(i) = log(i);
+  }
+
+  // Prepare keywords to pass to PolyCollection. See
+  // https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.fill_between.html
+  std::map<string, string> keywords;
+  keywords["alpha"] = "0.4";
+  keywords["color"] = "grey";
+  keywords["hatch"] = "-";
+
+  plt::fill_between(x, y, z, keywords);
+  plt::show();
+}
diff --git a/third_party/matplotlib-cpp/examples/imshow.cpp b/third_party/matplotlib-cpp/examples/imshow.cpp
new file mode 100644
index 0000000..b11661e
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/imshow.cpp
@@ -0,0 +1,29 @@
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <iostream>
+#include "../matplotlibcpp.h"
+
+using namespace std;
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    // Prepare data
+    int ncols = 500, nrows = 300;
+    std::vector<float> z(ncols * nrows);
+    for (int j=0; j<nrows; ++j) {
+        for (int i=0; i<ncols; ++i) {
+            z.at(ncols * j + i) = std::sin(std::hypot(i - ncols/2, j - nrows/2));
+        }
+    }
+
+    const float* zptr = &(z[0]);
+    const int colors = 1;
+
+    plt::title("My matrix");
+    plt::imshow(zptr, nrows, ncols, colors);
+
+    // Show plots
+    plt::save("imshow.png");
+    std::cout << "Result saved to 'imshow.png'.\n";
+}
diff --git a/third_party/matplotlib-cpp/examples/lines3d.cpp b/third_party/matplotlib-cpp/examples/lines3d.cpp
new file mode 100644
index 0000000..f3c201c
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/lines3d.cpp
@@ -0,0 +1,30 @@
+#include "../matplotlibcpp.h"
+
+#include <cmath>
+
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    std::vector<double> x, y, z;
+    double theta, r;
+    double z_inc = 4.0/99.0; double theta_inc = (8.0 * M_PI)/99.0;
+    
+    for (double i = 0; i < 100; i += 1) {
+        theta = -4.0 * M_PI + theta_inc*i;
+        z.push_back(-2.0 + z_inc*i);
+        r = z[i]*z[i] + 1;
+        x.push_back(r * sin(theta));
+        y.push_back(r * cos(theta));
+    }
+
+    std::map<std::string, std::string> keywords;
+    keywords.insert(std::pair<std::string, std::string>("label", "parametric curve") );
+
+    plt::plot3(x, y, z, keywords);
+    plt::xlabel("x label");
+    plt::ylabel("y label");
+    plt::set_zlabel("z label"); // set_zlabel rather than just zlabel, in accordance with the Axes3D method
+    plt::legend();
+    plt::show();
+}
diff --git a/third_party/matplotlib-cpp/examples/lines3d.png b/third_party/matplotlib-cpp/examples/lines3d.png
new file mode 100644
index 0000000..7a0c478
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/lines3d.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/quiver.cpp b/third_party/matplotlib-cpp/examples/quiver.cpp
new file mode 100644
index 0000000..ea3c3ec
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/quiver.cpp
@@ -0,0 +1,20 @@
+#include "../matplotlibcpp.h"
+
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    // u and v are respectively the x and y components of the arrows we're plotting
+    std::vector<int> x, y, u, v;
+    for (int i = -5; i <= 5; i++) {
+        for (int j = -5; j <= 5; j++) {
+            x.push_back(i);
+            u.push_back(-i);
+            y.push_back(j);
+            v.push_back(-j);
+        }
+    }
+
+    plt::quiver(x, y, u, v);
+    plt::show();
+}
\ No newline at end of file
diff --git a/third_party/matplotlib-cpp/examples/quiver.png b/third_party/matplotlib-cpp/examples/quiver.png
new file mode 100644
index 0000000..9d7be1e
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/quiver.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/subplot.cpp b/third_party/matplotlib-cpp/examples/subplot.cpp
new file mode 100644
index 0000000..bee322e
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/subplot.cpp
@@ -0,0 +1,31 @@
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include "../matplotlibcpp.h"
+
+using namespace std;
+namespace plt = matplotlibcpp;
+
+int main() 
+{
+    // Prepare data
+	int n = 500;
+	std::vector<double> x(n), y(n), z(n), w(n,2);
+	for(int i=0; i<n; ++i) {
+		x.at(i) = i;
+		y.at(i) = sin(2*M_PI*i/360.0);
+		z.at(i) = 100.0 / i;
+	}
+
+    // Set the "super title"
+    plt::suptitle("My plot");
+    plt::subplot(1, 2, 1);
+	plt::plot(x, y, "r-");
+    plt::subplot(1, 2, 2);
+    plt::plot(x, z, "k-");
+    // Add some text to the plot
+    plt::text(100, 90, "Hello!");
+
+
+	// Show plots
+	plt::show();
+}
diff --git a/third_party/matplotlib-cpp/examples/subplot2grid.cpp b/third_party/matplotlib-cpp/examples/subplot2grid.cpp
new file mode 100644
index 0000000..f590e51
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/subplot2grid.cpp
@@ -0,0 +1,44 @@
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include "../matplotlibcpp.h"
+
+using namespace std;
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    // Prepare data
+	int n = 500;
+	std::vector<double> x(n), u(n), v(n), w(n);
+	for(int i=0; i<n; ++i) {
+		x.at(i) = i;
+		u.at(i) = sin(2*M_PI*i/500.0);
+		v.at(i) = 100.0 / i;
+		w.at(i) = sin(2*M_PI*i/1000.0);
+	}
+
+    // Set the "super title"
+    plt::suptitle("My plot");
+
+    const long nrows=3, ncols=3;
+    long row = 2, col = 2;
+
+    plt::subplot2grid(nrows, ncols, row, col);
+	plt::plot(x, w, "g-");
+
+    long spanr = 1, spanc = 2;
+    col = 0;
+    plt::subplot2grid(nrows, ncols, row, col, spanr, spanc);
+	plt::plot(x, v, "r-");
+
+    spanr = 2, spanc = 3;
+    row = 0, col = 0;
+    plt::subplot2grid(nrows, ncols, row, col, spanr, spanc);
+    plt::plot(x, u, "b-");
+    // Add some text to the plot
+    plt::text(100., -0.5, "Hello!");
+
+
+    // Show plots
+	plt::show();
+}
diff --git a/third_party/matplotlib-cpp/examples/surface.cpp b/third_party/matplotlib-cpp/examples/surface.cpp
new file mode 100644
index 0000000..4865f06
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/surface.cpp
@@ -0,0 +1,24 @@
+#include "../matplotlibcpp.h"
+
+#include <cmath>
+
+namespace plt = matplotlibcpp;
+
+int main()
+{
+    std::vector<std::vector<double>> x, y, z;
+    for (double i = -5; i <= 5;  i += 0.25) {
+        std::vector<double> x_row, y_row, z_row;
+        for (double j = -5; j <= 5; j += 0.25) {
+            x_row.push_back(i);
+            y_row.push_back(j);
+            z_row.push_back(::std::sin(::std::hypot(i, j)));
+        }
+        x.push_back(x_row);
+        y.push_back(y_row);
+        z.push_back(z_row);
+    }
+
+    plt::plot_surface(x, y, z);
+    plt::show();
+}
diff --git a/third_party/matplotlib-cpp/examples/surface.png b/third_party/matplotlib-cpp/examples/surface.png
new file mode 100644
index 0000000..6fc5fc7
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/surface.png
Binary files differ
diff --git a/third_party/matplotlib-cpp/examples/update.cpp b/third_party/matplotlib-cpp/examples/update.cpp
new file mode 100644
index 0000000..64f4906
--- /dev/null
+++ b/third_party/matplotlib-cpp/examples/update.cpp
@@ -0,0 +1,60 @@
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include "../matplotlibcpp.h"
+#include <chrono>
+
+namespace plt = matplotlibcpp;
+
+void update_window(const double x, const double y, const double t,
+                   std::vector<double> &xt, std::vector<double> &yt)
+{
+    const double target_length = 300;
+    const double half_win = (target_length/(2.*sqrt(1.+t*t)));
+
+    xt[0] = x - half_win;
+    xt[1] = x + half_win;
+    yt[0] = y - half_win*t;
+    yt[1] = y + half_win*t;
+}
+
+
+int main()
+{
+    size_t n = 1000;
+    std::vector<double> x, y;
+
+    const double w = 0.05;
+    const double a = n/2;
+
+    for (size_t i=0; i<n; i++) {
+        x.push_back(i);
+        y.push_back(a*sin(w*i));
+    }
+
+    std::vector<double> xt(2), yt(2);
+
+    plt::title("Tangent of a sine curve");
+    plt::xlim(x.front(), x.back());
+    plt::ylim(-a, a);
+    plt::axis("equal");
+
+    // Plot sin once and for all.
+    plt::named_plot("sin", x, y);
+
+    // Prepare plotting the tangent.
+    plt::Plot plot("tangent");
+
+    plt::legend();
+
+    for (size_t i=0; i<n; i++) {
+        if (i % 10 == 0) {
+            update_window(x[i], y[i], a*w*cos(w*x[i]), xt, yt);
+
+            // Just update data for this plot.
+            plot.update(xt, yt);
+
+            // Small pause so the viewer has a chance to enjoy the animation.
+            plt::pause(0.1);
+        }
+   }
+}
diff --git a/third_party/matplotlib-cpp/examples/xkcd.cpp b/third_party/matplotlib-cpp/examples/xkcd.cpp
index 9bf6454..fa41cfc 100644
--- a/third_party/matplotlib-cpp/examples/xkcd.cpp
+++ b/third_party/matplotlib-cpp/examples/xkcd.cpp
@@ -1,6 +1,7 @@
+#define _USE_MATH_DEFINES
+#include <cmath>
 #include "../matplotlibcpp.h"
 #include <vector>
-#include <cmath>
 
 namespace plt = matplotlibcpp;
 
diff --git a/third_party/matplotlib-cpp/matplotlibcpp.h b/third_party/matplotlib-cpp/matplotlibcpp.h
index eb7b7cc..35aeb54 100644
--- a/third_party/matplotlib-cpp/matplotlibcpp.h
+++ b/third_party/matplotlib-cpp/matplotlibcpp.h
@@ -1,26 +1,41 @@
 #pragma once
 
+// Python headers must be included before any system headers, since
+// they define _POSIX_C_SOURCE
+#include <Python.h>
+
 #include <vector>
 #include <map>
+#include <array>
 #include <numeric>
 #include <algorithm>
 #include <stdexcept>
 #include <iostream>
-#include <stdint.h> // <cstdint> requires c++11 support
-
-#if __cplusplus > 199711L || _MSC_VER > 1800
-#  include <functional>
-#endif
-
-#include <Python.h>
+#include <cstdint> // <cstdint> requires c++11 support
+#include <functional>
 
 #ifndef WITHOUT_NUMPY
 #  define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #  include <numpy/arrayobject.h>
+
+#  ifdef WITH_OPENCV
+#    include <opencv2/opencv.hpp>
+#  endif // WITH_OPENCV
+
+/*
+ * A bunch of constants were removed in OpenCV 4 in favour of enum classes, so
+ * define the ones we need here.
+ */
+#  if CV_MAJOR_VERSION > 3
+#    define CV_BGR2RGB cv::COLOR_BGR2RGB
+#    define CV_BGRA2RGBA cv::COLOR_BGRA2RGBA
+#  endif
 #endif // WITHOUT_NUMPY
 
 #if PY_MAJOR_VERSION >= 3
 #  define PyString_FromString PyUnicode_FromString
+#  define PyInt_FromLong PyLong_FromLong
+#  define PyString_FromString PyUnicode_FromString
 #endif
 
 
@@ -30,28 +45,46 @@
 static std::string s_backend;
 
 struct _interpreter {
+    PyObject* s_python_function_arrow;
     PyObject *s_python_function_show;
     PyObject *s_python_function_close;
     PyObject *s_python_function_draw;
     PyObject *s_python_function_pause;
     PyObject *s_python_function_save;
     PyObject *s_python_function_figure;
+    PyObject *s_python_function_fignum_exists;
     PyObject *s_python_function_plot;
+    PyObject *s_python_function_quiver;
+    PyObject* s_python_function_contour;
     PyObject *s_python_function_semilogx;
     PyObject *s_python_function_semilogy;
     PyObject *s_python_function_loglog;
+    PyObject *s_python_function_fill;
     PyObject *s_python_function_fill_between;
     PyObject *s_python_function_hist;
+    PyObject *s_python_function_imshow;
+    PyObject *s_python_function_scatter;
+    PyObject *s_python_function_boxplot;
     PyObject *s_python_function_subplot;
+    PyObject *s_python_function_subplot2grid;
     PyObject *s_python_function_legend;
     PyObject *s_python_function_xlim;
     PyObject *s_python_function_ion;
+    PyObject *s_python_function_ginput;
     PyObject *s_python_function_ylim;
     PyObject *s_python_function_title;
     PyObject *s_python_function_axis;
+    PyObject *s_python_function_axvline;
+    PyObject *s_python_function_axvspan;
     PyObject *s_python_function_xlabel;
     PyObject *s_python_function_ylabel;
+    PyObject *s_python_function_gca;
+    PyObject *s_python_function_xticks;
+    PyObject *s_python_function_yticks;
+    PyObject* s_python_function_margins;
+    PyObject *s_python_function_tick_params;
     PyObject *s_python_function_grid;
+    PyObject* s_python_function_cla;
     PyObject *s_python_function_clf;
     PyObject *s_python_function_errorbar;
     PyObject *s_python_function_annotate;
@@ -60,6 +93,12 @@
     PyObject *s_python_empty_tuple;
     PyObject *s_python_function_stem;
     PyObject *s_python_function_xkcd;
+    PyObject *s_python_function_text;
+    PyObject *s_python_function_suptitle;
+    PyObject *s_python_function_bar;
+    PyObject *s_python_function_colorbar;
+    PyObject *s_python_function_subplots_adjust;
+
 
     /* For now, _interpreter is implemented as a singleton since its currently not possible to have
        multiple independent embedded python interpreters without patching the python source code
@@ -72,6 +111,18 @@
         return ctx;
     }
 
+    PyObject* safe_import(PyObject* module, std::string fname) {
+        PyObject* fn = PyObject_GetAttrString(module, fname.c_str());
+
+        if (!fn)
+            throw std::runtime_error(std::string("Couldn't find required function: ") + fname);
+
+        if (!PyFunction_Check(fn))
+            throw std::runtime_error(fname + std::string(" is unexpectedly not a PyFunction."));
+
+        return fn;
+    }
+
 private:
 
 #ifndef WITHOUT_NUMPY
@@ -92,6 +143,29 @@
 #endif
 
     _interpreter() {
+      // Force PYTHONHOME and PYTHONPATH to our sandboxed python.
+      wchar_t python_home[] = L"../python_repo/usr/";
+      wchar_t python_path[] =
+          L"../matplotlib_repo/3:../python_repo/usr/lib/python35.zip:../"
+          L"python_repo/usr/lib/python3.5:../python_repo/usr/lib/python3.5/"
+          L"plat-x86_64-linux-gnu:../python_repo/usr/lib/python3.5/"
+          L"lib-dynload:../python_repo/usr/lib/python3/dist-packages";
+
+      Py_SetPath(python_path);
+      Py_SetPythonHome(python_home);
+
+      // We fail really poorly if DISPLAY isn't set.  We can do better.
+      if (getenv("DISPLAY") == nullptr) {
+        fprintf(stderr, "DISPLAY not set\n");
+        abort();
+      }
+
+      // TODO(austin): Confirm LD_LIBRARY_PATH does the right thing.  Can't
+      // hurt.
+      setenv("LD_LIBRARY_PATH",
+             "../python_repo/lib/x86_64-linux-gnu:../python_repo/usr/lib:../"
+             "python_repo/usr/lib/x86_64-linux-gnu", 0);
+      Py_DontWriteBytecodeFlag = 1;
 
         // optional but recommended
 #if PY_MAJOR_VERSION >= 3
@@ -108,18 +182,18 @@
 
         PyObject* matplotlibname = PyString_FromString("matplotlib");
         PyObject* pyplotname = PyString_FromString("matplotlib.pyplot");
-        PyObject* mpl_toolkits = PyString_FromString("mpl_toolkits");
-        PyObject* axis3d = PyString_FromString("mpl_toolkits.mplot3d");
-        PyObject* pylabname  = PyString_FromString("pylab");
         PyObject* cmname  = PyString_FromString("matplotlib.cm");
-        if (!pyplotname || !pylabname || !matplotlibname || !mpl_toolkits ||
-            !axis3d || !cmname) {
-          throw std::runtime_error("couldnt create string");
+        PyObject* pylabname  = PyString_FromString("pylab");
+        if (!pyplotname || !pylabname || !matplotlibname || !cmname) {
+            throw std::runtime_error("couldnt create string");
         }
 
         PyObject* matplotlib = PyImport_Import(matplotlibname);
         Py_DECREF(matplotlibname);
-        if (!matplotlib) { throw std::runtime_error("Error loading module matplotlib!"); }
+        if (!matplotlib) {
+            PyErr_Print();
+            throw std::runtime_error("Error loading module matplotlib!");
+        }
 
         // matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
         // or matplotlib.backends is imported for the first time
@@ -139,102 +213,59 @@
         Py_DECREF(pylabname);
         if (!pylabmod) { throw std::runtime_error("Error loading module pylab!"); }
 
-        PyObject* mpl_toolkitsmod = PyImport_Import(mpl_toolkits);
-        Py_DECREF(mpl_toolkitsmod);
-        if (!mpl_toolkitsmod) { throw std::runtime_error("Error loading module mpl_toolkits!"); }
-
-        PyObject* axis3dmod = PyImport_Import(axis3d);
-        Py_DECREF(axis3dmod);
-        if (!axis3dmod) { throw std::runtime_error("Error loading module mpl_toolkits.mplot3d!"); }
-
-        s_python_function_show = PyObject_GetAttrString(pymod, "show");
-        s_python_function_close = PyObject_GetAttrString(pymod, "close");
-        s_python_function_draw = PyObject_GetAttrString(pymod, "draw");
-        s_python_function_pause = PyObject_GetAttrString(pymod, "pause");
-        s_python_function_figure = PyObject_GetAttrString(pymod, "figure");
-        s_python_function_plot = PyObject_GetAttrString(pymod, "plot");
-        s_python_function_semilogx = PyObject_GetAttrString(pymod, "semilogx");
-        s_python_function_semilogy = PyObject_GetAttrString(pymod, "semilogy");
-        s_python_function_loglog = PyObject_GetAttrString(pymod, "loglog");
-        s_python_function_fill_between = PyObject_GetAttrString(pymod, "fill_between");
-        s_python_function_hist = PyObject_GetAttrString(pymod,"hist");
-        s_python_function_subplot = PyObject_GetAttrString(pymod, "subplot");
-        s_python_function_legend = PyObject_GetAttrString(pymod, "legend");
-        s_python_function_ylim = PyObject_GetAttrString(pymod, "ylim");
-        s_python_function_title = PyObject_GetAttrString(pymod, "title");
-        s_python_function_axis = PyObject_GetAttrString(pymod, "axis");
-        s_python_function_xlabel = PyObject_GetAttrString(pymod, "xlabel");
-        s_python_function_ylabel = PyObject_GetAttrString(pymod, "ylabel");
-        s_python_function_grid = PyObject_GetAttrString(pymod, "grid");
-        s_python_function_xlim = PyObject_GetAttrString(pymod, "xlim");
-        s_python_function_ion = PyObject_GetAttrString(pymod, "ion");
-        s_python_function_save = PyObject_GetAttrString(pylabmod, "savefig");
-        s_python_function_annotate = PyObject_GetAttrString(pymod,"annotate");
-        s_python_function_clf = PyObject_GetAttrString(pymod, "clf");
-        s_python_function_errorbar = PyObject_GetAttrString(pymod, "errorbar");
-        s_python_function_tight_layout = PyObject_GetAttrString(pymod, "tight_layout");
-        s_python_function_stem = PyObject_GetAttrString(pymod, "stem");
-        s_python_function_xkcd = PyObject_GetAttrString(pymod, "xkcd");
-
-        if(    !s_python_function_show
-            || !s_python_function_close
-            || !s_python_function_draw
-            || !s_python_function_pause
-            || !s_python_function_figure
-            || !s_python_function_plot
-            || !s_python_function_semilogx
-            || !s_python_function_semilogy
-            || !s_python_function_loglog
-            || !s_python_function_fill_between
-            || !s_python_function_subplot
-            || !s_python_function_legend
-            || !s_python_function_ylim
-            || !s_python_function_title
-            || !s_python_function_axis
-            || !s_python_function_xlabel
-            || !s_python_function_ylabel
-            || !s_python_function_grid
-            || !s_python_function_xlim
-            || !s_python_function_ion
-            || !s_python_function_save
-            || !s_python_function_clf
-            || !s_python_function_annotate
-            || !s_python_function_errorbar
-            || !s_python_function_errorbar
-            || !s_python_function_tight_layout
-            || !s_python_function_stem
-            || !s_python_function_xkcd
-        ) { throw std::runtime_error("Couldn't find required function!"); }
-
-        if (   !PyFunction_Check(s_python_function_show)
-            || !PyFunction_Check(s_python_function_close)
-            || !PyFunction_Check(s_python_function_draw)
-            || !PyFunction_Check(s_python_function_pause)
-            || !PyFunction_Check(s_python_function_figure)
-            || !PyFunction_Check(s_python_function_plot)
-            || !PyFunction_Check(s_python_function_semilogx)
-            || !PyFunction_Check(s_python_function_semilogy)
-            || !PyFunction_Check(s_python_function_loglog)
-            || !PyFunction_Check(s_python_function_fill_between)
-            || !PyFunction_Check(s_python_function_subplot)
-            || !PyFunction_Check(s_python_function_legend)
-            || !PyFunction_Check(s_python_function_annotate)
-            || !PyFunction_Check(s_python_function_ylim)
-            || !PyFunction_Check(s_python_function_title)
-            || !PyFunction_Check(s_python_function_axis)
-            || !PyFunction_Check(s_python_function_xlabel)
-            || !PyFunction_Check(s_python_function_ylabel)
-            || !PyFunction_Check(s_python_function_grid)
-            || !PyFunction_Check(s_python_function_xlim)
-            || !PyFunction_Check(s_python_function_ion)
-            || !PyFunction_Check(s_python_function_save)
-            || !PyFunction_Check(s_python_function_clf)
-            || !PyFunction_Check(s_python_function_tight_layout)
-            || !PyFunction_Check(s_python_function_errorbar)
-            || !PyFunction_Check(s_python_function_stem)
-            || !PyFunction_Check(s_python_function_xkcd)
-        ) { throw std::runtime_error("Python object is unexpectedly not a PyFunction."); }
-
+        s_python_function_arrow = safe_import(pymod, "arrow");
+        s_python_function_show = safe_import(pymod, "show");
+        s_python_function_close = safe_import(pymod, "close");
+        s_python_function_draw = safe_import(pymod, "draw");
+        s_python_function_pause = safe_import(pymod, "pause");
+        s_python_function_figure = safe_import(pymod, "figure");
+        s_python_function_fignum_exists = safe_import(pymod, "fignum_exists");
+        s_python_function_plot = safe_import(pymod, "plot");
+        s_python_function_quiver = safe_import(pymod, "quiver");
+        s_python_function_contour = safe_import(pymod, "contour");
+        s_python_function_semilogx = safe_import(pymod, "semilogx");
+        s_python_function_semilogy = safe_import(pymod, "semilogy");
+        s_python_function_loglog = safe_import(pymod, "loglog");
+        s_python_function_fill = safe_import(pymod, "fill");
+        s_python_function_fill_between = safe_import(pymod, "fill_between");
+        s_python_function_hist = safe_import(pymod,"hist");
+        s_python_function_scatter = safe_import(pymod,"scatter");
+        s_python_function_boxplot = safe_import(pymod,"boxplot");
+        s_python_function_subplot = safe_import(pymod, "subplot");
+        s_python_function_subplot2grid = safe_import(pymod, "subplot2grid");
+        s_python_function_legend = safe_import(pymod, "legend");
+        s_python_function_ylim = safe_import(pymod, "ylim");
+        s_python_function_title = safe_import(pymod, "title");
+        s_python_function_axis = safe_import(pymod, "axis");
+        s_python_function_axvline = safe_import(pymod, "axvline");
+        s_python_function_axvspan = safe_import(pymod, "axvspan");
+        s_python_function_xlabel = safe_import(pymod, "xlabel");
+        s_python_function_ylabel = safe_import(pymod, "ylabel");
+        s_python_function_gca = safe_import(pymod, "gca");
+        s_python_function_xticks = safe_import(pymod, "xticks");
+        s_python_function_yticks = safe_import(pymod, "yticks");
+        s_python_function_margins = safe_import(pymod, "margins");
+        s_python_function_tick_params = safe_import(pymod, "tick_params");
+        s_python_function_grid = safe_import(pymod, "grid");
+        s_python_function_xlim = safe_import(pymod, "xlim");
+        s_python_function_ion = safe_import(pymod, "ion");
+        s_python_function_ginput = safe_import(pymod, "ginput");
+        s_python_function_save = safe_import(pylabmod, "savefig");
+        s_python_function_annotate = safe_import(pymod,"annotate");
+        s_python_function_cla = safe_import(pymod, "cla");
+        s_python_function_clf = safe_import(pymod, "clf");
+        s_python_function_errorbar = safe_import(pymod, "errorbar");
+        s_python_function_tight_layout = safe_import(pymod, "tight_layout");
+        s_python_function_stem = safe_import(pymod, "stem");
+        s_python_function_xkcd = safe_import(pymod, "xkcd");
+        s_python_function_text = safe_import(pymod, "text");
+        s_python_function_suptitle = safe_import(pymod, "suptitle");
+        s_python_function_bar = safe_import(pymod,"bar");
+        s_python_function_colorbar = PyObject_GetAttrString(pymod, "colorbar");
+        s_python_function_subplots_adjust = safe_import(pymod,"subplots_adjust");
+#ifndef WITHOUT_NUMPY
+        s_python_function_imshow = safe_import(pymod, "imshow");
+#endif
         s_python_empty_tuple = PyTuple_New(0);
     }
 
@@ -245,7 +276,15 @@
 
 } // end namespace detail
 
-// must be called before the first regular call to matplotlib to have any effect
+/// Select the backend
+///
+/// **NOTE:** This must be called before the first plot command to have
+/// any effect.
+///
+/// Mainly useful to select the non-interactive 'Agg' backend when running
+/// matplotlibcpp in headless mode, for example on a machine with no display.
+///
+/// See also: https://matplotlib.org/2.0.2/api/matplotlib_configuration_api.html#matplotlib.use
 inline void backend(const std::string& name)
 {
     detail::s_backend = name;
@@ -253,6 +292,8 @@
 
 inline bool annotate(std::string annotation, double x, double y)
 {
+    detail::_interpreter::get();
+
     PyObject * xy = PyTuple_New(2);
     PyObject * str = PyString_FromString(annotation.c_str());
 
@@ -275,6 +316,8 @@
     return res;
 }
 
+namespace detail {
+
 #ifndef WITHOUT_NUMPY
 // Type selector for numpy array conversion
 template <typename T> struct select_npy_type { const static NPY_TYPES type = NPY_NOTYPE; }; //Default
@@ -290,29 +333,37 @@
 template <> struct select_npy_type<uint32_t> { const static NPY_TYPES type = NPY_ULONG; };
 template <> struct select_npy_type<uint64_t> { const static NPY_TYPES type = NPY_UINT64; };
 
+// Sanity checks; comment them out or change the numpy type below if you're compiling on
+// a platform where they don't apply
+static_assert(sizeof(long long) == 8);
+template <> struct select_npy_type<long long> { const static NPY_TYPES type = NPY_INT64; };
+static_assert(sizeof(unsigned long long) == 8);
+template <> struct select_npy_type<unsigned long long> { const static NPY_TYPES type = NPY_UINT64; };
+// TODO: add int, long, etc.
+
 template<typename Numeric>
 PyObject* get_array(const std::vector<Numeric>& v)
 {
-    detail::_interpreter::get();    //interpreter needs to be initialized for the numpy commands to work
+    npy_intp vsize = v.size();
     NPY_TYPES type = select_npy_type<Numeric>::type;
-    if (type == NPY_NOTYPE)
-    {
-        std::vector<double> vd(v.size());
-        npy_intp vsize = v.size();
-        std::copy(v.begin(),v.end(),vd.begin());
-        PyObject* varray = PyArray_SimpleNewFromData(1, &vsize, NPY_DOUBLE, (void*)(vd.data()));
+    if (type == NPY_NOTYPE) {
+        size_t memsize = v.size()*sizeof(double);
+        double* dp = static_cast<double*>(::malloc(memsize));
+        for (size_t i=0; i<v.size(); ++i)
+            dp[i] = v[i];
+        PyObject* varray = PyArray_SimpleNewFromData(1, &vsize, NPY_DOUBLE, dp);
+        PyArray_UpdateFlags(reinterpret_cast<PyArrayObject*>(varray), NPY_ARRAY_OWNDATA);
         return varray;
     }
-
-    npy_intp vsize = v.size();
+    
     PyObject* varray = PyArray_SimpleNewFromData(1, &vsize, type, (void*)(v.data()));
     return varray;
 }
 
+
 template<typename Numeric>
 PyObject* get_2darray(const std::vector<::std::vector<Numeric>>& v)
 {
-    detail::_interpreter::get();    //interpreter needs to be initialized for the numpy commands to work
     if (v.size() < 1) throw std::runtime_error("get_2d_array v too small");
 
     npy_intp vsize[2] = {static_cast<npy_intp>(v.size()),
@@ -347,14 +398,42 @@
 
 #endif // WITHOUT_NUMPY
 
+// sometimes, for labels and such, we need string arrays
+inline PyObject * get_array(const std::vector<std::string>& strings)
+{
+  PyObject* list = PyList_New(strings.size());
+  for (std::size_t i = 0; i < strings.size(); ++i) {
+    PyList_SetItem(list, i, PyString_FromString(strings[i].c_str()));
+  }
+  return list;
+}
+
+// not all matplotlib need 2d arrays, some prefer lists of lists
+template<typename Numeric>
+PyObject* get_listlist(const std::vector<std::vector<Numeric>>& ll)
+{
+  PyObject* listlist = PyList_New(ll.size());
+  for (std::size_t i = 0; i < ll.size(); ++i) {
+    PyList_SetItem(listlist, i, get_array(ll[i]));
+  }
+  return listlist;
+}
+
+} // namespace detail
+
+/// Plot a line through the given x and y data points..
+/// 
+/// See: https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.pyplot.plot.html
 template<typename Numeric>
 bool plot(const std::vector<Numeric> &x, const std::vector<Numeric> &y, const std::map<std::string, std::string>& keywords)
 {
     assert(x.size() == y.size());
 
+    detail::_interpreter::get();
+
     // using numpy arrays
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     // construct positional args
     PyObject* args = PyTuple_New(2);
@@ -377,19 +456,46 @@
     return res;
 }
 
+// TODO - it should be possible to make this work by implementing
+// a non-numpy alternative for `detail::get_2darray()`.
+#ifndef WITHOUT_NUMPY
 template <typename Numeric>
 void plot_surface(const std::vector<::std::vector<Numeric>> &x,
                   const std::vector<::std::vector<Numeric>> &y,
                   const std::vector<::std::vector<Numeric>> &z,
                   const std::map<std::string, std::string> &keywords =
-                      std::map<std::string, std::string>()) {
+                      std::map<std::string, std::string>())
+{
+  detail::_interpreter::get();
+
+  // We lazily load the modules here the first time this function is called
+  // because I'm not sure that we can assume "matplotlib installed" implies
+  // "mpl_toolkits installed" on all platforms, and we don't want to require
+  // it for people who don't need 3d plots.
+  static PyObject *mpl_toolkitsmod = nullptr, *axis3dmod = nullptr;
+  if (!mpl_toolkitsmod) {
+    detail::_interpreter::get();
+
+    PyObject* mpl_toolkits = PyString_FromString("mpl_toolkits");
+    PyObject* axis3d = PyString_FromString("mpl_toolkits.mplot3d");
+    if (!mpl_toolkits || !axis3d) { throw std::runtime_error("couldnt create string"); }
+
+    mpl_toolkitsmod = PyImport_Import(mpl_toolkits);
+    Py_DECREF(mpl_toolkits);
+    if (!mpl_toolkitsmod) { throw std::runtime_error("Error loading module mpl_toolkits!"); }
+
+    axis3dmod = PyImport_Import(axis3d);
+    Py_DECREF(axis3d);
+    if (!axis3dmod) { throw std::runtime_error("Error loading module mpl_toolkits.mplot3d!"); }
+  }
+
   assert(x.size() == y.size());
   assert(y.size() == z.size());
 
   // using numpy arrays
-  PyObject *xarray = get_2darray(x);
-  PyObject *yarray = get_2darray(y);
-  PyObject *zarray = get_2darray(z);
+  PyObject *xarray = detail::get_2darray(x);
+  PyObject *yarray = detail::get_2darray(y);
+  PyObject *zarray = detail::get_2darray(z);
 
   // construct positional args
   PyObject *args = PyTuple_New(3);
@@ -446,15 +552,103 @@
   Py_DECREF(kwargs);
   if (res) Py_DECREF(res);
 }
+#endif // WITHOUT_NUMPY
+
+template <typename Numeric>
+void plot3(const std::vector<Numeric> &x,
+                  const std::vector<Numeric> &y,
+                  const std::vector<Numeric> &z,
+                  const std::map<std::string, std::string> &keywords =
+                      std::map<std::string, std::string>())
+{
+  detail::_interpreter::get();
+
+  // Same as with plot_surface: We lazily load the modules here the first time 
+  // this function is called because I'm not sure that we can assume "matplotlib 
+  // installed" implies "mpl_toolkits installed" on all platforms, and we don't 
+  // want to require it for people who don't need 3d plots.
+  static PyObject *mpl_toolkitsmod = nullptr, *axis3dmod = nullptr;
+  if (!mpl_toolkitsmod) {
+    detail::_interpreter::get();
+
+    PyObject* mpl_toolkits = PyString_FromString("mpl_toolkits");
+    PyObject* axis3d = PyString_FromString("mpl_toolkits.mplot3d");
+    if (!mpl_toolkits || !axis3d) { throw std::runtime_error("couldnt create string"); }
+
+    mpl_toolkitsmod = PyImport_Import(mpl_toolkits);
+    Py_DECREF(mpl_toolkits);
+    if (!mpl_toolkitsmod) { throw std::runtime_error("Error loading module mpl_toolkits!"); }
+
+    axis3dmod = PyImport_Import(axis3d);
+    Py_DECREF(axis3d);
+    if (!axis3dmod) { throw std::runtime_error("Error loading module mpl_toolkits.mplot3d!"); }
+  }
+
+  assert(x.size() == y.size());
+  assert(y.size() == z.size());
+
+  PyObject *xarray = detail::get_array(x);
+  PyObject *yarray = detail::get_array(y);
+  PyObject *zarray = detail::get_array(z);
+
+  // construct positional args
+  PyObject *args = PyTuple_New(3);
+  PyTuple_SetItem(args, 0, xarray);
+  PyTuple_SetItem(args, 1, yarray);
+  PyTuple_SetItem(args, 2, zarray);
+
+  // Build up the kw args.
+  PyObject *kwargs = PyDict_New();
+
+  for (std::map<std::string, std::string>::const_iterator it = keywords.begin();
+       it != keywords.end(); ++it) {
+    PyDict_SetItemString(kwargs, it->first.c_str(),
+                         PyString_FromString(it->second.c_str()));
+  }
+
+  PyObject *fig =
+      PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+                          detail::_interpreter::get().s_python_empty_tuple);
+  if (!fig) throw std::runtime_error("Call to figure() failed.");
+
+  PyObject *gca_kwargs = PyDict_New();
+  PyDict_SetItemString(gca_kwargs, "projection", PyString_FromString("3d"));
+
+  PyObject *gca = PyObject_GetAttrString(fig, "gca");
+  if (!gca) throw std::runtime_error("No gca");
+  Py_INCREF(gca);
+  PyObject *axis = PyObject_Call(
+      gca, detail::_interpreter::get().s_python_empty_tuple, gca_kwargs);
+
+  if (!axis) throw std::runtime_error("No axis");
+  Py_INCREF(axis);
+
+  Py_DECREF(gca);
+  Py_DECREF(gca_kwargs);
+
+  PyObject *plot3 = PyObject_GetAttrString(axis, "plot");
+  if (!plot3) throw std::runtime_error("No 3D line plot");
+  Py_INCREF(plot3);
+  PyObject *res = PyObject_Call(plot3, args, kwargs);
+  if (!res) throw std::runtime_error("Failed 3D line plot");
+  Py_DECREF(plot3);
+
+  Py_DECREF(axis);
+  Py_DECREF(args);
+  Py_DECREF(kwargs);
+  if (res) Py_DECREF(res);
+}
 
 template<typename Numeric>
 bool stem(const std::vector<Numeric> &x, const std::vector<Numeric> &y, const std::map<std::string, std::string>& keywords)
 {
     assert(x.size() == y.size());
 
+    detail::_interpreter::get();
+
     // using numpy arrays
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     // construct positional args
     PyObject* args = PyTuple_New(2);
@@ -481,15 +675,49 @@
 }
 
 template< typename Numeric >
+bool fill(const std::vector<Numeric>& x, const std::vector<Numeric>& y, const std::map<std::string, std::string>& keywords)
+{
+    assert(x.size() == y.size());
+
+    detail::_interpreter::get();
+
+    // using numpy arrays
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    // construct positional args
+    PyObject* args = PyTuple_New(2);
+    PyTuple_SetItem(args, 0, xarray);
+    PyTuple_SetItem(args, 1, yarray);
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_fill, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if (res) Py_DECREF(res);
+
+    return res;
+}
+
+template< typename Numeric >
 bool fill_between(const std::vector<Numeric>& x, const std::vector<Numeric>& y1, const std::vector<Numeric>& y2, const std::map<std::string, std::string>& keywords)
 {
     assert(x.size() == y1.size());
     assert(x.size() == y2.size());
 
+    detail::_interpreter::get();
+
     // using numpy arrays
-    PyObject* xarray = get_array(x);
-    PyObject* y1array = get_array(y1);
-    PyObject* y2array = get_array(y2);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* y1array = detail::get_array(y1);
+    PyObject* y2array = detail::get_array(y2);
 
     // construct positional args
     PyObject* args = PyTuple_New(3);
@@ -499,8 +727,7 @@
 
     // construct keyword args
     PyObject* kwargs = PyDict_New();
-    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
-    {
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it) {
         PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
     }
 
@@ -513,17 +740,50 @@
     return res;
 }
 
-template< typename Numeric>
-bool hist(const std::vector<Numeric>& y, long bins=10,std::string color="b", double alpha=1.0)
-{
+template <typename Numeric>
+bool arrow(Numeric x, Numeric y, Numeric end_x, Numeric end_y, const std::string& fc = "r",
+           const std::string ec = "k", Numeric head_length = 0.25, Numeric head_width = 0.1625) {
+    PyObject* obj_x = PyFloat_FromDouble(x);
+    PyObject* obj_y = PyFloat_FromDouble(y);
+    PyObject* obj_end_x = PyFloat_FromDouble(end_x);
+    PyObject* obj_end_y = PyFloat_FromDouble(end_y);
 
-    PyObject* yarray = get_array(y);
+    PyObject* kwargs = PyDict_New();
+    PyDict_SetItemString(kwargs, "fc", PyString_FromString(fc.c_str()));
+    PyDict_SetItemString(kwargs, "ec", PyString_FromString(ec.c_str()));
+    PyDict_SetItemString(kwargs, "head_width", PyFloat_FromDouble(head_width));
+    PyDict_SetItemString(kwargs, "head_length", PyFloat_FromDouble(head_length));
+
+    PyObject* plot_args = PyTuple_New(4);
+    PyTuple_SetItem(plot_args, 0, obj_x);
+    PyTuple_SetItem(plot_args, 1, obj_y);
+    PyTuple_SetItem(plot_args, 2, obj_end_x);
+    PyTuple_SetItem(plot_args, 3, obj_end_y);
+
+    PyObject* res =
+            PyObject_Call(detail::_interpreter::get().s_python_function_arrow, plot_args, kwargs);
+
+    Py_DECREF(plot_args);
+    Py_DECREF(kwargs);
+    if (res)
+        Py_DECREF(res);
+
+    return res;
+}
+
+template< typename Numeric>
+bool hist(const std::vector<Numeric>& y, long bins=10,std::string color="b",
+          double alpha=1.0, bool cumulative=false)
+{
+    detail::_interpreter::get();
+
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "bins", PyLong_FromLong(bins));
     PyDict_SetItemString(kwargs, "color", PyString_FromString(color.c_str()));
     PyDict_SetItemString(kwargs, "alpha", PyFloat_FromDouble(alpha));
-
+    PyDict_SetItemString(kwargs, "cumulative", cumulative ? Py_True : Py_False);
 
     PyObject* plot_args = PyTuple_New(1);
 
@@ -540,10 +800,263 @@
     return res;
 }
 
+#ifndef WITHOUT_NUMPY
+namespace detail {
+
+inline void imshow(void *ptr, const NPY_TYPES type, const int rows, const int columns, const int colors, const std::map<std::string, std::string> &keywords, PyObject** out)
+{
+    assert(type == NPY_UINT8 || type == NPY_FLOAT);
+    assert(colors == 1 || colors == 3 || colors == 4);
+
+    detail::_interpreter::get();
+
+    // construct args
+    npy_intp dims[3] = { rows, columns, colors };
+    PyObject *args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyArray_SimpleNewFromData(colors == 1 ? 2 : 3, dims, type, ptr));
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject *res = PyObject_Call(detail::_interpreter::get().s_python_function_imshow, args, kwargs);
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if (!res)
+        throw std::runtime_error("Call to imshow() failed");
+    if (out)
+        *out = res;
+    else
+        Py_DECREF(res);
+}
+
+} // namespace detail
+
+inline void imshow(const unsigned char *ptr, const int rows, const int columns, const int colors, const std::map<std::string, std::string> &keywords = {}, PyObject** out = nullptr)
+{
+    detail::imshow((void *) ptr, NPY_UINT8, rows, columns, colors, keywords, out);
+}
+
+inline void imshow(const float *ptr, const int rows, const int columns, const int colors, const std::map<std::string, std::string> &keywords = {}, PyObject** out = nullptr)
+{
+    detail::imshow((void *) ptr, NPY_FLOAT, rows, columns, colors, keywords, out);
+}
+
+#ifdef WITH_OPENCV
+void imshow(const cv::Mat &image, const std::map<std::string, std::string> &keywords = {})
+{
+    // Convert underlying type of matrix, if needed
+    cv::Mat image2;
+    NPY_TYPES npy_type = NPY_UINT8;
+    switch (image.type() & CV_MAT_DEPTH_MASK) {
+    case CV_8U:
+        image2 = image;
+        break;
+    case CV_32F:
+        image2 = image;
+        npy_type = NPY_FLOAT;
+        break;
+    default:
+        image.convertTo(image2, CV_MAKETYPE(CV_8U, image.channels()));
+    }
+
+    // If color image, convert from BGR to RGB
+    switch (image2.channels()) {
+    case 3:
+        cv::cvtColor(image2, image2, CV_BGR2RGB);
+        break;
+    case 4:
+        cv::cvtColor(image2, image2, CV_BGRA2RGBA);
+    }
+
+    detail::imshow(image2.data, npy_type, image2.rows, image2.cols, image2.channels(), keywords);
+}
+#endif // WITH_OPENCV
+#endif // WITHOUT_NUMPY
+
+template<typename NumericX, typename NumericY>
+bool scatter(const std::vector<NumericX>& x,
+             const std::vector<NumericY>& y,
+             const double s=1.0, // The marker size in points**2
+             const std::map<std::string, std::string> & keywords = {})
+{
+    detail::_interpreter::get();
+
+    assert(x.size() == y.size());
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* kwargs = PyDict_New();
+    PyDict_SetItemString(kwargs, "s", PyLong_FromLong(s));
+    for (const auto& it : keywords)
+    {
+        PyDict_SetItemString(kwargs, it.first.c_str(), PyString_FromString(it.second.c_str()));
+    }
+
+    PyObject* plot_args = PyTuple_New(2);
+    PyTuple_SetItem(plot_args, 0, xarray);
+    PyTuple_SetItem(plot_args, 1, yarray);
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_scatter, plot_args, kwargs);
+
+    Py_DECREF(plot_args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename Numeric>
+bool boxplot(const std::vector<std::vector<Numeric>>& data,
+             const std::vector<std::string>& labels = {},
+             const std::map<std::string, std::string> & keywords = {})
+{
+    detail::_interpreter::get();
+
+    PyObject* listlist = detail::get_listlist(data);
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, listlist);
+
+    PyObject* kwargs = PyDict_New();
+
+    // kwargs needs the labels, if there are (the correct number of) labels
+    if (!labels.empty() && labels.size() == data.size()) {
+        PyDict_SetItemString(kwargs, "labels", detail::get_array(labels));
+    }
+
+    // take care of the remaining keywords
+    for (const auto& it : keywords)
+    {
+        PyDict_SetItemString(kwargs, it.first.c_str(), PyString_FromString(it.second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_boxplot, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename Numeric>
+bool boxplot(const std::vector<Numeric>& data,
+             const std::map<std::string, std::string> & keywords = {})
+{
+    detail::_interpreter::get();
+
+    PyObject* vector = detail::get_array(data);
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, vector);
+
+    PyObject* kwargs = PyDict_New();
+    for (const auto& it : keywords)
+    {
+        PyDict_SetItemString(kwargs, it.first.c_str(), PyString_FromString(it.second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_boxplot, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template <typename Numeric>
+bool bar(const std::vector<Numeric> &               x,
+         const std::vector<Numeric> &               y,
+         std::string                                ec       = "black",
+         std::string                                ls       = "-",
+         double                                     lw       = 1.0,
+         const std::map<std::string, std::string> & keywords = {})
+{
+  detail::_interpreter::get();
+
+  PyObject * xarray = detail::get_array(x);
+  PyObject * yarray = detail::get_array(y);
+
+  PyObject * kwargs = PyDict_New();
+
+  PyDict_SetItemString(kwargs, "ec", PyString_FromString(ec.c_str()));
+  PyDict_SetItemString(kwargs, "ls", PyString_FromString(ls.c_str()));
+  PyDict_SetItemString(kwargs, "lw", PyFloat_FromDouble(lw));
+
+  for (std::map<std::string, std::string>::const_iterator it =
+         keywords.begin();
+       it != keywords.end();
+       ++it) {
+    PyDict_SetItemString(
+      kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+  }
+
+  PyObject * plot_args = PyTuple_New(2);
+  PyTuple_SetItem(plot_args, 0, xarray);
+  PyTuple_SetItem(plot_args, 1, yarray);
+
+  PyObject * res = PyObject_Call(
+    detail::_interpreter::get().s_python_function_bar, plot_args, kwargs);
+
+  Py_DECREF(plot_args);
+  Py_DECREF(kwargs);
+  if (res) Py_DECREF(res);
+
+  return res;
+}
+
+template <typename Numeric>
+bool bar(const std::vector<Numeric> &               y,
+         std::string                                ec       = "black",
+         std::string                                ls       = "-",
+         double                                     lw       = 1.0,
+         const std::map<std::string, std::string> & keywords = {})
+{
+  using T = typename std::remove_reference<decltype(y)>::type::value_type;
+
+  detail::_interpreter::get();
+
+  std::vector<T> x;
+  for (std::size_t i = 0; i < y.size(); i++) { x.push_back(i); }
+
+  return bar(x, y, ec, ls, lw, keywords);
+}
+
+inline bool subplots_adjust(const std::map<std::string, double>& keywords = {})
+{
+    detail::_interpreter::get();
+
+    PyObject* kwargs = PyDict_New();
+    for (std::map<std::string, double>::const_iterator it =
+            keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(),
+                             PyFloat_FromDouble(it->second));
+    }
+
+
+    PyObject* plot_args = PyTuple_New(0);
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_subplots_adjust, plot_args, kwargs);
+
+    Py_DECREF(plot_args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
 template< typename Numeric>
 bool named_hist(std::string label,const std::vector<Numeric>& y, long bins=10, std::string color="b", double alpha=1.0)
 {
-    PyObject* yarray = get_array(y);
+    detail::_interpreter::get();
+
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "label", PyString_FromString(label.c_str()));
@@ -569,8 +1082,10 @@
 {
     assert(x.size() == y.size());
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(s.c_str());
 
@@ -587,13 +1102,84 @@
     return res;
 }
 
+template <typename NumericX, typename NumericY, typename NumericZ>
+bool contour(const std::vector<NumericX>& x, const std::vector<NumericY>& y,
+             const std::vector<NumericZ>& z,
+             const std::map<std::string, std::string>& keywords = {}) {
+    assert(x.size() == y.size() && x.size() == z.size());
+
+    PyObject* xarray = get_array(x);
+    PyObject* yarray = get_array(y);
+    PyObject* zarray = get_array(z);
+
+    PyObject* plot_args = PyTuple_New(3);
+    PyTuple_SetItem(plot_args, 0, xarray);
+    PyTuple_SetItem(plot_args, 1, yarray);
+    PyTuple_SetItem(plot_args, 2, zarray);
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for (std::map<std::string, std::string>::const_iterator it = keywords.begin();
+         it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res =
+            PyObject_Call(detail::_interpreter::get().s_python_function_contour, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+    if (res)
+        Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY, typename NumericU, typename NumericW>
+bool quiver(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::vector<NumericU>& u, const std::vector<NumericW>& w, const std::map<std::string, std::string>& keywords = {})
+{
+    assert(x.size() == y.size() && x.size() == u.size() && u.size() == w.size());
+
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+    PyObject* uarray = detail::get_array(u);
+    PyObject* warray = detail::get_array(w);
+
+    PyObject* plot_args = PyTuple_New(4);
+    PyTuple_SetItem(plot_args, 0, xarray);
+    PyTuple_SetItem(plot_args, 1, yarray);
+    PyTuple_SetItem(plot_args, 2, uarray);
+    PyTuple_SetItem(plot_args, 3, warray);
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(
+            detail::_interpreter::get().s_python_function_quiver, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+    if (res)
+        Py_DECREF(res);
+
+    return res;
+}
+
 template<typename NumericX, typename NumericY>
 bool stem(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
 {
     assert(x.size() == y.size());
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(s.c_str());
 
@@ -617,8 +1203,10 @@
 {
     assert(x.size() == y.size());
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(s.c_str());
 
@@ -640,8 +1228,10 @@
 {
     assert(x.size() == y.size());
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(s.c_str());
 
@@ -663,8 +1253,10 @@
 {
     assert(x.size() == y.size());
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(s.c_str());
 
@@ -682,20 +1274,25 @@
 }
 
 template<typename NumericX, typename NumericY>
-bool errorbar(const std::vector<NumericX> &x, const std::vector<NumericY> &y, const std::vector<NumericX> &yerr, const std::string & /*s*/ = "")
+bool errorbar(const std::vector<NumericX> &x, const std::vector<NumericY> &y, const std::vector<NumericX> &yerr, const std::map<std::string, std::string> &keywords = {})
 {
     assert(x.size() == y.size());
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
-    PyObject* yerrarray = get_array(yerr);
+    detail::_interpreter::get();
 
-    PyObject *kwargs = PyDict_New();
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+    PyObject* yerrarray = detail::get_array(yerr);
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
 
     PyDict_SetItemString(kwargs, "yerr", yerrarray);
 
-    //PyObject *pystring = PyString_FromString(s.c_str());
-
     PyObject *plot_args = PyTuple_New(2);
     PyTuple_SetItem(plot_args, 0, xarray);
     PyTuple_SetItem(plot_args, 1, yarray);
@@ -716,10 +1313,12 @@
 template<typename Numeric>
 bool named_plot(const std::string& name, const std::vector<Numeric>& y, const std::string& format = "")
 {
+    detail::_interpreter::get();
+
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "label", PyString_FromString(name.c_str()));
 
-    PyObject* yarray = get_array(y);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(format.c_str());
 
@@ -740,11 +1339,13 @@
 template<typename Numeric>
 bool named_plot(const std::string& name, const std::vector<Numeric>& x, const std::vector<Numeric>& y, const std::string& format = "")
 {
+    detail::_interpreter::get();
+
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "label", PyString_FromString(name.c_str()));
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(format.c_str());
 
@@ -765,11 +1366,13 @@
 template<typename Numeric>
 bool named_semilogx(const std::string& name, const std::vector<Numeric>& x, const std::vector<Numeric>& y, const std::string& format = "")
 {
+    detail::_interpreter::get();
+
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "label", PyString_FromString(name.c_str()));
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(format.c_str());
 
@@ -790,11 +1393,13 @@
 template<typename Numeric>
 bool named_semilogy(const std::string& name, const std::vector<Numeric>& x, const std::vector<Numeric>& y, const std::string& format = "")
 {
+    detail::_interpreter::get();
+
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "label", PyString_FromString(name.c_str()));
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(format.c_str());
 
@@ -815,11 +1420,13 @@
 template<typename Numeric>
 bool named_loglog(const std::string& name, const std::vector<Numeric>& x, const std::vector<Numeric>& y, const std::string& format = "")
 {
+    detail::_interpreter::get();
+
     PyObject* kwargs = PyDict_New();
     PyDict_SetItemString(kwargs, "label", PyString_FromString(name.c_str()));
 
-    PyObject* xarray = get_array(x);
-    PyObject* yarray = get_array(y);
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
 
     PyObject* pystring = PyString_FromString(format.c_str());
 
@@ -827,7 +1434,6 @@
     PyTuple_SetItem(plot_args, 0, xarray);
     PyTuple_SetItem(plot_args, 1, yarray);
     PyTuple_SetItem(plot_args, 2, pystring);
-
     PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_loglog, plot_args, kwargs);
 
     Py_DECREF(kwargs);
@@ -846,6 +1452,14 @@
 }
 
 template<typename Numeric>
+bool plot(const std::vector<Numeric>& y, const std::map<std::string, std::string>& keywords)
+{
+    std::vector<Numeric> x(y.size());
+    for(size_t i=0; i<x.size(); ++i) x.at(i) = i;
+    return plot(x,y,keywords);
+}
+
+template<typename Numeric>
 bool stem(const std::vector<Numeric>& y, const std::string& format = "")
 {
     std::vector<Numeric> x(y.size());
@@ -853,25 +1467,150 @@
     return stem(x, y, format);
 }
 
-inline void figure()
+template<typename Numeric>
+void text(Numeric x, Numeric y, const std::string& s = "")
 {
-    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure, detail::_interpreter::get().s_python_empty_tuple);
+    detail::_interpreter::get();
+
+    PyObject* args = PyTuple_New(3);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(x));
+    PyTuple_SetItem(args, 1, PyFloat_FromDouble(y));
+    PyTuple_SetItem(args, 2, PyString_FromString(s.c_str()));
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_text, args);
+    if(!res) throw std::runtime_error("Call to text() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+inline void colorbar(PyObject* mappable = NULL, const std::map<std::string, float>& keywords = {})
+{
+    if (mappable == NULL)
+        throw std::runtime_error("Must call colorbar with PyObject* returned from an image, contour, surface, etc.");
+
+    detail::_interpreter::get();
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, mappable);
+
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, float>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyFloat_FromDouble(it->second));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_colorbar, args, kwargs);
+    if(!res) throw std::runtime_error("Call to colorbar() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    Py_DECREF(res);
+}
+
+
+inline long figure(long number = -1)
+{
+    detail::_interpreter::get();
+
+    PyObject *res;
+    if (number == -1)
+        res = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure, detail::_interpreter::get().s_python_empty_tuple);
+    else {
+        assert(number > 0);
+
+        // Make sure interpreter is initialised
+        detail::_interpreter::get();
+
+        PyObject *args = PyTuple_New(1);
+        PyTuple_SetItem(args, 0, PyLong_FromLong(number));
+        res = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure, args);
+        Py_DECREF(args);
+    }
+
     if(!res) throw std::runtime_error("Call to figure() failed.");
 
+    PyObject* num = PyObject_GetAttrString(res, "number");
+    if (!num) throw std::runtime_error("Could not get number attribute of figure object");
+    const long figureNumber = PyLong_AsLong(num);
+
+    Py_DECREF(num);
+    Py_DECREF(res);
+
+    return figureNumber;
+}
+
+inline bool fignum_exists(long number)
+{
+    detail::_interpreter::get();
+
+    PyObject *args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyLong_FromLong(number));
+    PyObject *res = PyObject_CallObject(detail::_interpreter::get().s_python_function_fignum_exists, args);
+    if(!res) throw std::runtime_error("Call to fignum_exists() failed.");
+
+    bool ret = PyObject_IsTrue(res);
+    Py_DECREF(res);
+    Py_DECREF(args);
+
+    return ret;
+}
+
+inline void figure_size(size_t w, size_t h)
+{
+    detail::_interpreter::get();
+
+    const size_t dpi = 100;
+    PyObject* size = PyTuple_New(2);
+    PyTuple_SetItem(size, 0, PyFloat_FromDouble((double)w / dpi));
+    PyTuple_SetItem(size, 1, PyFloat_FromDouble((double)h / dpi));
+
+    PyObject* kwargs = PyDict_New();
+    PyDict_SetItemString(kwargs, "figsize", size);
+    PyDict_SetItemString(kwargs, "dpi", PyLong_FromSize_t(dpi));
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_figure,
+            detail::_interpreter::get().s_python_empty_tuple, kwargs);
+
+    Py_DECREF(kwargs);
+
+    if(!res) throw std::runtime_error("Call to figure_size() failed.");
     Py_DECREF(res);
 }
 
 inline void legend()
 {
+    detail::_interpreter::get();
+
     PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_legend, detail::_interpreter::get().s_python_empty_tuple);
     if(!res) throw std::runtime_error("Call to legend() failed.");
 
     Py_DECREF(res);
 }
 
+inline void legend(const std::map<std::string, std::string>& keywords)
+{
+  detail::_interpreter::get();
+
+  // construct keyword args
+  PyObject* kwargs = PyDict_New();
+  for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+  {
+    PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+  }
+
+  PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_legend, detail::_interpreter::get().s_python_empty_tuple, kwargs);
+  if(!res) throw std::runtime_error("Call to legend() failed.");
+
+  Py_DECREF(kwargs);
+  Py_DECREF(res);  
+}
+
 template<typename Numeric>
 void ylim(Numeric left, Numeric right)
 {
+    detail::_interpreter::get();
+
     PyObject* list = PyList_New(2);
     PyList_SetItem(list, 0, PyFloat_FromDouble(left));
     PyList_SetItem(list, 1, PyFloat_FromDouble(right));
@@ -889,6 +1628,8 @@
 template<typename Numeric>
 void xlim(Numeric left, Numeric right)
 {
+    detail::_interpreter::get();
+
     PyObject* list = PyList_New(2);
     PyList_SetItem(list, 0, PyFloat_FromDouble(left));
     PyList_SetItem(list, 1, PyFloat_FromDouble(right));
@@ -906,6 +1647,8 @@
 
 inline double* xlim()
 {
+    detail::_interpreter::get();
+
     PyObject* args = PyTuple_New(0);
     PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_xlim, args);
     PyObject* left = PyTuple_GetItem(res,0);
@@ -924,6 +1667,8 @@
 
 inline double* ylim()
 {
+    detail::_interpreter::get();
+
     PyObject* args = PyTuple_New(0);
     PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_ylim, args);
     PyObject* left = PyTuple_GetItem(res,0);
@@ -939,8 +1684,166 @@
     return arr;
 }
 
+template<typename Numeric>
+inline void xticks(const std::vector<Numeric> &ticks, const std::vector<std::string> &labels = {}, const std::map<std::string, std::string>& keywords = {})
+{
+    assert(labels.size() == 0 || ticks.size() == labels.size());
+
+    detail::_interpreter::get();
+
+    // using numpy array
+    PyObject* ticksarray = detail::get_array(ticks);
+
+    PyObject* args;
+    if(labels.size() == 0) {
+        // construct positional args
+        args = PyTuple_New(1);
+        PyTuple_SetItem(args, 0, ticksarray);
+    } else {
+        // make tuple of tick labels
+        PyObject* labelstuple = PyTuple_New(labels.size());
+        for (size_t i = 0; i < labels.size(); i++)
+            PyTuple_SetItem(labelstuple, i, PyUnicode_FromString(labels[i].c_str()));
+
+        // construct positional args
+        args = PyTuple_New(2);
+        PyTuple_SetItem(args, 0, ticksarray);
+        PyTuple_SetItem(args, 1, labelstuple);
+    }
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_xticks, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if(!res) throw std::runtime_error("Call to xticks() failed");
+
+    Py_DECREF(res);
+}
+
+template<typename Numeric>
+inline void xticks(const std::vector<Numeric> &ticks, const std::map<std::string, std::string>& keywords)
+{
+    xticks(ticks, {}, keywords);
+}
+
+template<typename Numeric>
+inline void yticks(const std::vector<Numeric> &ticks, const std::vector<std::string> &labels = {}, const std::map<std::string, std::string>& keywords = {})
+{
+    assert(labels.size() == 0 || ticks.size() == labels.size());
+
+    detail::_interpreter::get();
+
+    // using numpy array
+    PyObject* ticksarray = detail::get_array(ticks);
+
+    PyObject* args;
+    if(labels.size() == 0) {
+        // construct positional args
+        args = PyTuple_New(1);
+        PyTuple_SetItem(args, 0, ticksarray);
+    } else {
+        // make tuple of tick labels
+        PyObject* labelstuple = PyTuple_New(labels.size());
+        for (size_t i = 0; i < labels.size(); i++)
+            PyTuple_SetItem(labelstuple, i, PyUnicode_FromString(labels[i].c_str()));
+
+        // construct positional args
+        args = PyTuple_New(2);
+        PyTuple_SetItem(args, 0, ticksarray);
+        PyTuple_SetItem(args, 1, labelstuple);
+    }
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_yticks, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if(!res) throw std::runtime_error("Call to yticks() failed");
+
+    Py_DECREF(res);
+}
+
+template<typename Numeric>
+inline void yticks(const std::vector<Numeric> &ticks, const std::map<std::string, std::string>& keywords)
+{
+    yticks(ticks, {}, keywords);
+}
+
+template <typename Numeric> inline void margins(Numeric margin)
+{
+    // construct positional args
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(margin));
+
+    PyObject* res =
+            PyObject_CallObject(detail::_interpreter::get().s_python_function_margins, args);
+    if (!res)
+        throw std::runtime_error("Call to margins() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+template <typename Numeric> inline void margins(Numeric margin_x, Numeric margin_y)
+{
+    // construct positional args
+    PyObject* args = PyTuple_New(2);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(margin_x));
+    PyTuple_SetItem(args, 1, PyFloat_FromDouble(margin_y));
+
+    PyObject* res =
+            PyObject_CallObject(detail::_interpreter::get().s_python_function_margins, args);
+    if (!res)
+        throw std::runtime_error("Call to margins() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+
+inline void tick_params(const std::map<std::string, std::string>& keywords, const std::string axis = "both")
+{
+  detail::_interpreter::get();
+
+  // construct positional args
+  PyObject* args;
+  args = PyTuple_New(1);
+  PyTuple_SetItem(args, 0, PyString_FromString(axis.c_str()));
+
+  // construct keyword args
+  PyObject* kwargs = PyDict_New();
+  for (std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+  {
+    PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+  }
+
+
+  PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_tick_params, args, kwargs);
+
+  Py_DECREF(args);
+  Py_DECREF(kwargs);
+  if (!res) throw std::runtime_error("Call to tick_params() failed");
+
+  Py_DECREF(res);
+}
+
 inline void subplot(long nrows, long ncols, long plot_number)
 {
+    detail::_interpreter::get();
+    
     // construct positional args
     PyObject* args = PyTuple_New(3);
     PyTuple_SetItem(args, 0, PyFloat_FromDouble(nrows));
@@ -954,21 +1857,79 @@
     Py_DECREF(res);
 }
 
-inline void title(const std::string &titlestr)
+inline void subplot2grid(long nrows, long ncols, long rowid=0, long colid=0, long rowspan=1, long colspan=1)
 {
+    detail::_interpreter::get();
+
+    PyObject* shape = PyTuple_New(2);
+    PyTuple_SetItem(shape, 0, PyLong_FromLong(nrows));
+    PyTuple_SetItem(shape, 1, PyLong_FromLong(ncols));
+
+    PyObject* loc = PyTuple_New(2);
+    PyTuple_SetItem(loc, 0, PyLong_FromLong(rowid));
+    PyTuple_SetItem(loc, 1, PyLong_FromLong(colid));
+
+    PyObject* args = PyTuple_New(4);
+    PyTuple_SetItem(args, 0, shape);
+    PyTuple_SetItem(args, 1, loc);
+    PyTuple_SetItem(args, 2, PyLong_FromLong(rowspan));
+    PyTuple_SetItem(args, 3, PyLong_FromLong(colspan));
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_subplot2grid, args);
+    if(!res) throw std::runtime_error("Call to subplot2grid() failed.");
+
+    Py_DECREF(shape);
+    Py_DECREF(loc);
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+inline void title(const std::string &titlestr, const std::map<std::string, std::string> &keywords = {})
+{
+    detail::_interpreter::get();
+
     PyObject* pytitlestr = PyString_FromString(titlestr.c_str());
     PyObject* args = PyTuple_New(1);
     PyTuple_SetItem(args, 0, pytitlestr);
 
-    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_title, args);
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_title, args, kwargs);
     if(!res) throw std::runtime_error("Call to title() failed.");
 
     Py_DECREF(args);
+    Py_DECREF(kwargs);
+    Py_DECREF(res);
+}
+
+inline void suptitle(const std::string &suptitlestr, const std::map<std::string, std::string> &keywords = {})
+{
+    detail::_interpreter::get();
+    
+    PyObject* pysuptitlestr = PyString_FromString(suptitlestr.c_str());
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, pysuptitlestr);
+
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_suptitle, args, kwargs);
+    if(!res) throw std::runtime_error("Call to suptitle() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
     Py_DECREF(res);
 }
 
 inline void axis(const std::string &axisstr)
 {
+    detail::_interpreter::get();
+
     PyObject* str = PyString_FromString(axisstr.c_str());
     PyObject* args = PyTuple_New(1);
     PyTuple_SetItem(args, 0, str);
@@ -980,34 +1941,155 @@
     Py_DECREF(res);
 }
 
-inline void xlabel(const std::string &str)
+inline void axvline(double x, double ymin = 0., double ymax = 1., const std::map<std::string, std::string>& keywords = std::map<std::string, std::string>())
 {
+    detail::_interpreter::get();
+
+    // construct positional args
+    PyObject* args = PyTuple_New(3);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(x));
+    PyTuple_SetItem(args, 1, PyFloat_FromDouble(ymin));
+    PyTuple_SetItem(args, 2, PyFloat_FromDouble(ymax));
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_axvline, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if(res) Py_DECREF(res);
+}
+
+inline void axvspan(double xmin, double xmax, double ymin = 0., double ymax = 1., const std::map<std::string, std::string>& keywords = std::map<std::string, std::string>())
+{
+    // construct positional args
+    PyObject* args = PyTuple_New(4);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(xmin));
+    PyTuple_SetItem(args, 1, PyFloat_FromDouble(xmax));
+    PyTuple_SetItem(args, 2, PyFloat_FromDouble(ymin));
+    PyTuple_SetItem(args, 3, PyFloat_FromDouble(ymax));
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+	if (it->first == "linewidth" || it->first == "alpha")
+    	    PyDict_SetItemString(kwargs, it->first.c_str(), PyFloat_FromDouble(std::stod(it->second)));
+  	else
+    	    PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_axvspan, args, kwargs);
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if(res) Py_DECREF(res);
+}
+
+inline void xlabel(const std::string &str, const std::map<std::string, std::string> &keywords = {})
+{
+    detail::_interpreter::get();
+
     PyObject* pystr = PyString_FromString(str.c_str());
     PyObject* args = PyTuple_New(1);
     PyTuple_SetItem(args, 0, pystr);
 
-    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_xlabel, args);
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_xlabel, args, kwargs);
     if(!res) throw std::runtime_error("Call to xlabel() failed.");
 
     Py_DECREF(args);
+    Py_DECREF(kwargs);
     Py_DECREF(res);
 }
 
-inline void ylabel(const std::string &str)
+inline void ylabel(const std::string &str, const std::map<std::string, std::string>& keywords = {})
 {
+    detail::_interpreter::get();
+
     PyObject* pystr = PyString_FromString(str.c_str());
     PyObject* args = PyTuple_New(1);
     PyTuple_SetItem(args, 0, pystr);
 
-    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_ylabel, args);
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_ylabel, args, kwargs);
     if(!res) throw std::runtime_error("Call to ylabel() failed.");
 
     Py_DECREF(args);
+    Py_DECREF(kwargs);
     Py_DECREF(res);
 }
 
+inline void set_zlabel(const std::string &str, const std::map<std::string, std::string>& keywords = {})
+{
+    detail::_interpreter::get();
+
+    // Same as with plot_surface: We lazily load the modules here the first time 
+    // this function is called because I'm not sure that we can assume "matplotlib 
+    // installed" implies "mpl_toolkits installed" on all platforms, and we don't 
+    // want to require it for people who don't need 3d plots.
+    static PyObject *mpl_toolkitsmod = nullptr, *axis3dmod = nullptr;
+    if (!mpl_toolkitsmod) {
+        PyObject* mpl_toolkits = PyString_FromString("mpl_toolkits");
+        PyObject* axis3d = PyString_FromString("mpl_toolkits.mplot3d");
+        if (!mpl_toolkits || !axis3d) { throw std::runtime_error("couldnt create string"); }
+
+        mpl_toolkitsmod = PyImport_Import(mpl_toolkits);
+        Py_DECREF(mpl_toolkits);
+        if (!mpl_toolkitsmod) { throw std::runtime_error("Error loading module mpl_toolkits!"); }
+
+        axis3dmod = PyImport_Import(axis3d);
+        Py_DECREF(axis3d);
+        if (!axis3dmod) { throw std::runtime_error("Error loading module mpl_toolkits.mplot3d!"); }
+    }
+
+    PyObject* pystr = PyString_FromString(str.c_str());
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, pystr);
+
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject *ax =
+    PyObject_CallObject(detail::_interpreter::get().s_python_function_gca,
+      detail::_interpreter::get().s_python_empty_tuple);
+    if (!ax) throw std::runtime_error("Call to gca() failed.");
+    Py_INCREF(ax);
+
+    PyObject *zlabel = PyObject_GetAttrString(ax, "set_zlabel");
+    if (!zlabel) throw std::runtime_error("Attribute set_zlabel not found.");
+    Py_INCREF(zlabel);
+
+    PyObject *res = PyObject_Call(zlabel, args, kwargs);
+    if (!res) throw std::runtime_error("Call to set_zlabel() failed.");
+    Py_DECREF(zlabel);
+
+    Py_DECREF(ax);
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if (res) Py_DECREF(res);
+}
+
 inline void grid(bool flag)
 {
+    detail::_interpreter::get();
+
     PyObject* pyflag = flag ? Py_True : Py_False;
     Py_INCREF(pyflag);
 
@@ -1023,6 +2105,8 @@
 
 inline void show(const bool block = true)
 {
+    detail::_interpreter::get();
+
     PyObject* res;
     if(block)
     {
@@ -1035,7 +2119,7 @@
         PyObject *kwargs = PyDict_New();
         PyDict_SetItemString(kwargs, "block", Py_False);
         res = PyObject_Call( detail::_interpreter::get().s_python_function_show, detail::_interpreter::get().s_python_empty_tuple, kwargs);
-	Py_DECREF(kwargs);
+       Py_DECREF(kwargs);
     }
 
 
@@ -1046,6 +2130,8 @@
 
 inline void close()
 {
+    detail::_interpreter::get();
+
     PyObject* res = PyObject_CallObject(
             detail::_interpreter::get().s_python_function_close,
             detail::_interpreter::get().s_python_empty_tuple);
@@ -1056,6 +2142,8 @@
 }
 
 inline void xkcd() {
+    detail::_interpreter::get();
+
     PyObject* res;
     PyObject *kwargs = PyDict_New();
 
@@ -1072,6 +2160,8 @@
 
 inline void draw()
 {
+    detail::_interpreter::get();
+
     PyObject* res = PyObject_CallObject(
         detail::_interpreter::get().s_python_function_draw,
         detail::_interpreter::get().s_python_empty_tuple);
@@ -1084,6 +2174,8 @@
 template<typename Numeric>
 inline void pause(Numeric interval)
 {
+    detail::_interpreter::get();
+
     PyObject* args = PyTuple_New(1);
     PyTuple_SetItem(args, 0, PyFloat_FromDouble(interval));
 
@@ -1096,6 +2188,8 @@
 
 inline void save(const std::string& filename)
 {
+    detail::_interpreter::get();
+
     PyObject* pyfilename = PyString_FromString(filename.c_str());
 
     PyObject* args = PyTuple_New(1);
@@ -1109,6 +2203,8 @@
 }
 
 inline void clf() {
+    detail::_interpreter::get();
+
     PyObject *res = PyObject_CallObject(
         detail::_interpreter::get().s_python_function_clf,
         detail::_interpreter::get().s_python_empty_tuple);
@@ -1118,7 +2214,21 @@
     Py_DECREF(res);
 }
 
-    inline void ion() {
+inline void cla() {
+    detail::_interpreter::get();
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_cla,
+                                        detail::_interpreter::get().s_python_empty_tuple);
+
+    if (!res)
+        throw std::runtime_error("Call to cla() failed.");
+
+    Py_DECREF(res);
+}
+
+inline void ion() {
+    detail::_interpreter::get();
+
     PyObject *res = PyObject_CallObject(
         detail::_interpreter::get().s_python_function_ion,
         detail::_interpreter::get().s_python_empty_tuple);
@@ -1128,8 +2238,46 @@
     Py_DECREF(res);
 }
 
+inline std::vector<std::array<double, 2>> ginput(const int numClicks = 1, const std::map<std::string, std::string>& keywords = {})
+{
+    detail::_interpreter::get();
+
+    PyObject *args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyLong_FromLong(numClicks));
+
+    // construct keyword args
+    PyObject* kwargs = PyDict_New();
+    for(std::map<std::string, std::string>::const_iterator it = keywords.begin(); it != keywords.end(); ++it)
+    {
+        PyDict_SetItemString(kwargs, it->first.c_str(), PyUnicode_FromString(it->second.c_str()));
+    }
+
+    PyObject* res = PyObject_Call(
+        detail::_interpreter::get().s_python_function_ginput, args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(args);
+    if (!res) throw std::runtime_error("Call to ginput() failed.");
+
+    const size_t len = PyList_Size(res);
+    std::vector<std::array<double, 2>> out;
+    out.reserve(len);
+    for (size_t i = 0; i < len; i++) {
+        PyObject *current = PyList_GetItem(res, i);
+        std::array<double, 2> position;
+        position[0] = PyFloat_AsDouble(PyTuple_GetItem(current, 0));
+        position[1] = PyFloat_AsDouble(PyTuple_GetItem(current, 1));
+        out.push_back(position);
+    }
+    Py_DECREF(res);
+
+    return out;
+}
+
 // Actually, is there any reason not to call this automatically for every plot?
 inline void tight_layout() {
+    detail::_interpreter::get();
+
     PyObject *res = PyObject_CallObject(
         detail::_interpreter::get().s_python_function_tight_layout,
         detail::_interpreter::get().s_python_empty_tuple);
@@ -1139,8 +2287,7 @@
     Py_DECREF(res);
 }
 
-#if __cplusplus > 199711L || _MSC_VER > 1800
-// C++11-exclusive content starts here (variadic plot() and initializer list support)
+// Support for variadic plot() and initializer lists:
 
 namespace detail {
 
@@ -1269,6 +2416,106 @@
     return plot<double>(x,y,keywords);
 }
 
-#endif
+/*
+ * This class allows dynamic plots, ie changing the plotted data without clearing and re-plotting
+ */
+class Plot
+{
+public:
+    // default initialization with plot label, some data and format
+    template<typename Numeric>
+    Plot(const std::string& name, const std::vector<Numeric>& x, const std::vector<Numeric>& y, const std::string& format = "") {
+        detail::_interpreter::get();
+
+        assert(x.size() == y.size());
+
+        PyObject* kwargs = PyDict_New();
+        if(name != "")
+            PyDict_SetItemString(kwargs, "label", PyString_FromString(name.c_str()));
+
+        PyObject* xarray = detail::get_array(x);
+        PyObject* yarray = detail::get_array(y);
+
+        PyObject* pystring = PyString_FromString(format.c_str());
+
+        PyObject* plot_args = PyTuple_New(3);
+        PyTuple_SetItem(plot_args, 0, xarray);
+        PyTuple_SetItem(plot_args, 1, yarray);
+        PyTuple_SetItem(plot_args, 2, pystring);
+
+        PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_plot, plot_args, kwargs);
+
+        Py_DECREF(kwargs);
+        Py_DECREF(plot_args);
+
+        if(res)
+        {
+            line= PyList_GetItem(res, 0);
+
+            if(line)
+                set_data_fct = PyObject_GetAttrString(line,"set_data");
+            else
+                Py_DECREF(line);
+            Py_DECREF(res);
+        }
+    }
+
+    // shorter initialization with name or format only
+    // basically calls line, = plot([], [])
+    Plot(const std::string& name = "", const std::string& format = "")
+        : Plot(name, std::vector<double>(), std::vector<double>(), format) {}
+
+    template<typename Numeric>
+    bool update(const std::vector<Numeric>& x, const std::vector<Numeric>& y) {
+        assert(x.size() == y.size());
+        if(set_data_fct)
+        {
+            PyObject* xarray = detail::get_array(x);
+            PyObject* yarray = detail::get_array(y);
+
+            PyObject* plot_args = PyTuple_New(2);
+            PyTuple_SetItem(plot_args, 0, xarray);
+            PyTuple_SetItem(plot_args, 1, yarray);
+
+            PyObject* res = PyObject_CallObject(set_data_fct, plot_args);
+            if (res) Py_DECREF(res);
+            return res;
+        }
+        return false;
+    }
+
+    // clears the plot but keep it available
+    bool clear() {
+        return update(std::vector<double>(), std::vector<double>());
+    }
+
+    // definitely remove this line
+    void remove() {
+        if(line)
+        {
+            auto remove_fct = PyObject_GetAttrString(line,"remove");
+            PyObject* args = PyTuple_New(0);
+            PyObject* res = PyObject_CallObject(remove_fct, args);
+            if (res) Py_DECREF(res);
+        }
+        decref();
+    }
+
+    ~Plot() {
+        decref();
+    }
+private:
+
+    void decref() {
+        if(line)
+            Py_DECREF(line);
+        if(set_data_fct)
+            Py_DECREF(set_data_fct);
+    }
+
+
+    PyObject* line = nullptr;
+    PyObject* set_data_fct = nullptr;
+};
 
 } // end namespace matplotlibcpp
diff --git a/third_party/matplotlib-cpp/numpy_flags.py b/third_party/matplotlib-cpp/numpy_flags.py
new file mode 100644
index 0000000..56fd95c
--- /dev/null
+++ b/third_party/matplotlib-cpp/numpy_flags.py
@@ -0,0 +1,12 @@
+from os import path
+
+try:
+    from numpy import __file__ as numpyloc
+
+    # Get numpy directory
+    numpy_dir = path.dirname(numpyloc)
+
+    # Print the result of joining this to core and include
+    print("-I" + path.join(numpy_dir, "core", "include"))
+except:
+    print("-DWITHOUT_NUMPY")