w il y a 1 mois
commit
cdc1a40219

+ 45 - 0
Gravity_filtering/CMakeLists.txt

@@ -0,0 +1,45 @@
+cmake_minimum_required(VERSION 3.10)
+project(SineFit)
+
+set(CMAKE_CXX_STANDARD 17)
+
+# 查找依赖包
+find_package(Eigen3 REQUIRED)
+find_package(Ceres REQUIRED)
+
+# 查找 Python 和 NumPy
+find_package(Python3 COMPONENTS Interpreter Development REQUIRED)
+
+# 获取 NumPy 的 include 路径
+execute_process(
+        COMMAND ${Python3_EXECUTABLE} -c "import numpy; print(numpy.get_include())"
+        OUTPUT_VARIABLE NUMPY_INCLUDE_DIR
+        OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+# 下载 matplotlibcpp.h(如果没有)
+if(NOT EXISTS "${CMAKE_SOURCE_DIR}/include/matplotlibcpp.h")
+    file(DOWNLOAD
+            "https://raw.githubusercontent.com/lava/matplotlib-cpp/master/matplotlibcpp.h"
+            "${CMAKE_SOURCE_DIR}/include/matplotlibcpp.h"
+    )
+endif()
+
+# 添加头文件搜索路径
+include_directories(
+        ${Python3_INCLUDE_DIRS}      # Python 头文件
+        ${NUMPY_INCLUDE_DIR}         # NumPy 头文件
+        ${CMAKE_SOURCE_DIR}/include  # matplotlibcpp.h
+        ${EIGEN3_INCLUDE_DIRS}       # Eigen
+)
+
+# 创建可执行文件(只定义一次)
+add_executable(sine_fit main.cpp)
+
+# 链接库
+target_link_libraries(sine_fit
+        PRIVATE
+        ${CERES_LIBRARIES}
+        ${EIGEN3_LIBRARIES}
+        ${Python3_LIBRARIES}  # 使用 Python3_LIBRARIES 而不是 PYTHON_LIBRARIES
+)

+ 9 - 0
Gravity_filtering/README

@@ -0,0 +1,9 @@
+此程序所实现功能是:使用ceres库根据噪声点进行曲线拟合,实现图形化并显示误差百分比
+使用方式:1.确保理想点坐标文件(.csv),以及噪声点坐标文件(.csv)在运行目录下
+        2.file1读取的为噪声点坐标,file2读取的为理想点坐标
+        3.运行代码:mkdir build
+                cd build
+                cmake ..
+                make
+                ./sine_fit
+注意:此程序是实现正弦函数拟合,若需修改拟合函数,需修改ceres库的残差结构体函数,拟合函数的参数个数可修改

+ 301 - 0
Gravity_filtering/f(1)_data.csv

@@ -0,0 +1,301 @@
+x,y_actual
+0.0,0.0
+0.0634665182543393,0.0634239196565645
+0.126933036508679,0.126592453573749
+0.190399554763018,0.18925124436041
+0.253866073017357,0.251147987181079
+0.317332591271696,0.312033445698487
+0.380799109526036,0.371662455660328
+0.444265627780375,0.429794912089172
+0.507732146034714,0.486196736100469
+0.571198664289053,0.540640817455598
+0.634665182543393,0.59290792905464
+0.698131700797732,0.642787609686539
+0.761598219052071,0.690079011482112
+0.82506473730641,0.734591708657533
+0.88853125556075,0.776146464291757
+0.951997773815089,0.814575952050336
+1.01546429206943,0.849725429949514
+1.07893081032377,0.881453363447582
+1.14239732857811,0.909631995354518
+1.20586384683245,0.934147860265107
+1.26933036508679,0.954902241444074
+1.33279688334112,0.971811568323542
+1.39626340159546,0.984807753012208
+1.4597299198498,0.993838464461254
+1.52319643810414,0.998867339183008
+1.58666295635848,0.999874127673875
+1.65012947461282,0.996854775951942
+1.71359599286716,0.989821441880933
+1.7770625111215,0.978802446214779
+1.84052902937584,0.963842158559942
+1.90399554763018,0.945000818714669
+1.96746206588452,0.922354294104582
+2.03092858413886,0.895993774291336
+2.0943951023932,0.866025403784439
+2.15786162064753,0.832569854634771
+2.22132813890187,0.795761840530832
+2.28479465715621,0.755749574354258
+2.34826117541055,0.712694171378863
+2.41172769366489,0.666769000516292
+2.47519421191923,0.618158986220606
+2.53866073017357,0.567059863862771
+2.60212724842791,0.513677391573407
+2.66559376668225,0.458226521727411
+2.72906028493659,0.400930535406614
+2.79252680319093,0.342020143325669
+2.85599332144527,0.28173255684143
+2.91945983969961,0.220310532786541
+2.98292635795395,0.15800139597335
+3.04639287620828,0.0950560433041829
+3.10985939446262,0.0317279334980681
+3.17332591271696,-0.0317279334980679
+3.2367924309713,-0.0950560433041826
+3.30025894922564,-0.15800139597335
+3.36372546747998,-0.220310532786541
+3.42719198573432,-0.281732556841429
+3.49065850398866,-0.342020143325669
+3.554125022243,-0.400930535406613
+3.61759154049734,-0.45822652172741
+3.68105805875168,-0.513677391573406
+3.74452457700602,-0.567059863862771
+3.80799109526036,-0.618158986220605
+3.87145761351469,-0.666769000516292
+3.93492413176903,-0.712694171378863
+3.99839065002337,-0.755749574354258
+4.06185716827771,-0.795761840530832
+4.12532368653205,-0.832569854634771
+4.18879020478639,-0.866025403784439
+4.25225672304073,-0.895993774291336
+4.31572324129507,-0.922354294104581
+4.37918975954941,-0.945000818714668
+4.44265627780375,-0.963842158559942
+4.50612279605809,-0.978802446214779
+4.56958931431243,-0.989821441880933
+4.63305583256677,-0.996854775951942
+4.69652235082111,-0.999874127673875
+4.75998886907544,-0.998867339183008
+4.82345538732978,-0.993838464461254
+4.88692190558412,-0.984807753012208
+4.95038842383846,-0.971811568323542
+5.0138549420928,-0.954902241444074
+5.07732146034714,-0.934147860265107
+5.14078797860148,-0.909631995354519
+5.20425449685582,-0.881453363447582
+5.26772101511016,-0.849725429949514
+5.3311875333645,-0.814575952050336
+5.39465405161884,-0.776146464291757
+5.45812056987318,-0.734591708657534
+5.52158708812751,-0.690079011482113
+5.58505360638185,-0.64278760968654
+5.64852012463619,-0.59290792905464
+5.71198664289053,-0.540640817455597
+5.77545316114487,-0.486196736100469
+5.83891967939921,-0.429794912089172
+5.90238619765355,-0.371662455660328
+5.96585271590789,-0.312033445698487
+6.02931923416223,-0.251147987181079
+6.09278575241657,-0.189251244360411
+6.15625227067091,-0.12659245357375
+6.21971878892525,-0.0634239196565654
+6.28318530717959,-2.44929359829471e-16
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,
+,

+ 101 - 0
Gravity_filtering/f_data.csv

@@ -0,0 +1,101 @@
+x,y_actual
+0.0,0.140285069715603
+0.0634665182543393,0.0876477305985184
+0.126933036508679,0.298436000276415
+0.190399554763018,0.267918124582501
+0.253866073017357,0.284264373251415
+0.317332591271696,0.438192330289456
+0.380799109526036,0.523268017499199
+0.444265627780375,0.625359558521007
+0.507732146034714,0.286405686243247
+0.571198664289053,0.686816253860807
+0.634665182543393,0.63793451684824
+0.698131700797732,0.838767691970072
+0.761598219052071,0.701151039217489
+0.82506473730641,0.726401062741621
+0.88853125556075,0.896685506500538
+0.951997773815089,0.705713126332752
+1.01546429206943,0.84896314642807
+1.07893081032377,1.04179435886038
+1.14239732857811,0.939496483006593
+1.20586384683245,1.07221913428672
+1.26933036508679,1.05035835824223
+1.33279688334112,1.00620638265413
+1.39626340159546,0.883501563406598
+1.4597299198498,1.06040495138904
+1.52319643810414,0.832260464624057
+1.58666295635848,1.05025804174251
+1.65012947461282,1.06123259913088
+1.71359599286716,1.08172218400782
+1.7770625111215,1.13510329274491
+1.84052902937584,1.15676344771338
+1.90399554763018,1.05261245284903
+1.96746206588452,0.954932889254741
+2.03092858413886,1.06731899921701
+2.0943951023932,0.898061550087815
+2.15786162064753,0.639363029969676
+2.22132813890187,0.644105668970256
+2.28479465715621,0.900833861834126
+2.34826117541055,0.706412775863704
+2.41172769366489,0.804711270346797
+2.47519421191923,0.501921019828979
+2.53866073017357,0.587976400478281
+2.60212724842791,0.565630745599175
+2.66559376668225,0.271022928032437
+2.72906028493659,0.44681590305347
+2.79252680319093,0.28698472823489
+2.85599332144527,0.101545588458254
+2.91945983969961,0.21613852845747
+2.98292635795395,0.0350055543981798
+3.04639287620828,-0.055710457677439
+3.10985939446262,-0.08607439813886
+3.17332591271696,-0.173121969252112
+3.2367924309713,-0.219427173515137
+3.30025894922564,-0.340940431608893
+3.36372546747998,-0.166231366042588
+3.42719198573432,-0.368985814489257
+3.49065850398866,-0.326581472107533
+3.554125022243,-0.322865319628881
+3.61759154049734,-0.458580116334374
+3.68105805875168,-0.499356969272961
+3.74452457700602,-0.588986597744354
+3.80799109526036,-0.768586075181377
+3.87145761351469,-0.670626083129084
+3.93492413176903,-0.571494909242536
+3.99839065002337,-0.606178612009565
+4.06185716827771,-0.887644107613753
+4.12532368653205,-0.949185311134246
+4.18879020478639,-0.840033575489158
+4.25225672304073,-0.839869044226233
+4.31572324129507,-0.955542713447427
+4.37918975954941,-1.06261061250177
+4.44265627780375,-0.784668910042675
+4.50612279605809,-1.14597396337569
+4.56958931431243,-1.14753767124824
+4.63305583256677,-1.14003832719034
+4.69652235082111,-1.13328995132331
+4.75998886907544,-0.950483881608885
+4.82345538732978,-0.964354558524775
+4.88692190558412,-1.16397659689786
+4.95038842383846,-0.799331014480242
+5.0138549420928,-0.863437568772765
+5.07732146034714,-0.839011198746071
+5.14078797860148,-1.08427019507739
+5.20425449685582,-0.73727713823229
+5.26772101511016,-0.675963382365029
+5.3311875333645,-0.620816627153947
+5.39465405161884,-0.632570937618211
+5.45812056987318,-0.620368112951522
+5.52158708812751,-0.684728044047083
+5.58505360638185,-0.771746625484194
+5.64852012463619,-0.633472130360303
+5.71198664289053,-0.687068317060409
+5.77545316114487,-0.673840916602488
+5.83891967939921,-0.254138229661353
+5.90238619765355,-0.451140029825771
+5.96585271590789,-0.393819911908345
+6.02931923416223,-0.317973474446609
+6.09278575241657,-0.20242396954887
+6.15625227067091,-0.0673130909872871
+6.21971878892525,-0.253332647059351
+6.28318530717959,0.136882644967734

+ 2986 - 0
Gravity_filtering/include/matplotlibcpp.h

@@ -0,0 +1,2986 @@
+#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 <cstdint> // <cstdint> requires c++11 support
+#include <functional>
+#include <string> // std::stod
+
+#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
+
+
+namespace matplotlibcpp {
+namespace detail {
+
+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_axhline;
+    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;
+    PyObject *s_python_function_tight_layout;
+    PyObject *s_python_colormap;
+    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_barh;
+    PyObject *s_python_function_colorbar;
+    PyObject *s_python_function_subplots_adjust;
+    PyObject *s_python_function_rcparams;
+    PyObject *s_python_function_spy;
+
+    /* 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
+       or starting a separate process for each. [1]
+       Furthermore, many python objects expect that they are destructed in the same thread as they
+       were constructed. [2] So for advanced usage, a `kill()` function is provided so that library
+       users can manually ensure that the interpreter is constructed and destroyed within the
+       same thread.
+
+         1: http://bytes.com/topic/python/answers/793370-multiple-independent-python-interpreters-c-c-program
+         2: https://github.com/lava/matplotlib-cpp/pull/202#issue-436220256
+       */
+
+    static _interpreter& get() {
+        return interkeeper(false);
+    }
+
+    static _interpreter& kill() {
+        return interkeeper(true);
+    }
+
+    // Stores the actual singleton object referenced by `get()` and `kill()`.
+    static _interpreter& interkeeper(bool should_kill) {
+        static _interpreter ctx;
+        if (should_kill)
+            ctx.~_interpreter();
+        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
+#  if PY_MAJOR_VERSION >= 3
+
+    void *import_numpy() {
+        import_array(); // initialize C-API
+        return NULL;
+    }
+
+#  else
+
+    void import_numpy() {
+        import_array(); // initialize C-API
+    }
+
+#  endif
+#endif
+
+    _interpreter() {
+
+        // optional but recommended
+#if PY_MAJOR_VERSION >= 3
+        wchar_t name[] = L"plotting";
+#else
+        char name[] = "plotting";
+#endif
+        Py_SetProgramName(name);
+        Py_Initialize();
+
+        wchar_t const *dummy_args[] = {L"Python", NULL};  // const is needed because literals must not be modified
+        wchar_t const **argv = dummy_args;
+        int             argc = sizeof(dummy_args)/sizeof(dummy_args[0])-1;
+
+#if PY_MAJOR_VERSION >= 3
+        PySys_SetArgv(argc, const_cast<wchar_t **>(argv));
+#else
+        PySys_SetArgv(argc, (char **)(argv));
+#endif
+
+#ifndef WITHOUT_NUMPY
+        import_numpy(); // initialize numpy C-API
+#endif
+
+        PyObject* matplotlibname = PyString_FromString("matplotlib");
+        PyObject* pyplotname = PyString_FromString("matplotlib.pyplot");
+        PyObject* cmname  = PyString_FromString("matplotlib.cm");
+        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) {
+            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
+        if (!s_backend.empty()) {
+            PyObject_CallMethod(matplotlib, const_cast<char*>("use"), const_cast<char*>("s"), s_backend.c_str());
+        }
+
+
+
+        PyObject* pymod = PyImport_Import(pyplotname);
+        Py_DECREF(pyplotname);
+        if (!pymod) { throw std::runtime_error("Error loading module matplotlib.pyplot!"); }
+
+        s_python_colormap = PyImport_Import(cmname);
+        Py_DECREF(cmname);
+        if (!s_python_colormap) { throw std::runtime_error("Error loading module matplotlib.cm!"); }
+
+        PyObject* pylabmod = PyImport_Import(pylabname);
+        Py_DECREF(pylabname);
+        if (!pylabmod) { throw std::runtime_error("Error loading module pylab!"); }
+
+        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_xlim = safe_import(pymod, "xlim");
+        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_axhline = safe_import(pymod, "axhline");
+        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_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_barh = safe_import(pymod, "barh");
+        s_python_function_colorbar = PyObject_GetAttrString(pymod, "colorbar");
+        s_python_function_subplots_adjust = safe_import(pymod,"subplots_adjust");
+        s_python_function_rcparams = PyObject_GetAttrString(pymod, "rcParams");
+	s_python_function_spy = PyObject_GetAttrString(pymod, "spy");
+#ifndef WITHOUT_NUMPY
+        s_python_function_imshow = safe_import(pymod, "imshow");
+#endif
+        s_python_empty_tuple = PyTuple_New(0);
+    }
+
+    ~_interpreter() {
+        Py_Finalize();
+    }
+};
+
+} // end namespace detail
+
+/// 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;
+}
+
+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());
+
+    PyTuple_SetItem(xy,0,PyFloat_FromDouble(x));
+    PyTuple_SetItem(xy,1,PyFloat_FromDouble(y));
+
+    PyObject* kwargs = PyDict_New();
+    PyDict_SetItemString(kwargs, "xy", xy);
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, str);
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_annotate, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if(res) Py_DECREF(res);
+
+    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
+template <> struct select_npy_type<double> { const static NPY_TYPES type = NPY_DOUBLE; };
+template <> struct select_npy_type<float> { const static NPY_TYPES type = NPY_FLOAT; };
+template <> struct select_npy_type<bool> { const static NPY_TYPES type = NPY_BOOL; };
+template <> struct select_npy_type<int8_t> { const static NPY_TYPES type = NPY_INT8; };
+template <> struct select_npy_type<int16_t> { const static NPY_TYPES type = NPY_SHORT; };
+template <> struct select_npy_type<int32_t> { const static NPY_TYPES type = NPY_INT; };
+template <> struct select_npy_type<int64_t> { const static NPY_TYPES type = NPY_INT64; };
+template <> struct select_npy_type<uint8_t> { const static NPY_TYPES type = NPY_UINT8; };
+template <> struct select_npy_type<uint16_t> { const static NPY_TYPES type = NPY_USHORT; };
+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; };
+
+template<typename Numeric>
+PyObject* get_array(const std::vector<Numeric>& v)
+{
+    npy_intp vsize = v.size();
+    NPY_TYPES type = select_npy_type<Numeric>::type;
+    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;
+    }
+
+    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)
+{
+    if (v.size() < 1) throw std::runtime_error("get_2d_array v too small");
+
+    npy_intp vsize[2] = {static_cast<npy_intp>(v.size()),
+                         static_cast<npy_intp>(v[0].size())};
+
+    PyArrayObject *varray =
+        (PyArrayObject *)PyArray_SimpleNew(2, vsize, NPY_DOUBLE);
+
+    double *vd_begin = static_cast<double *>(PyArray_DATA(varray));
+
+    for (const ::std::vector<Numeric> &v_row : v) {
+      if (v_row.size() != static_cast<size_t>(vsize[1]))
+        throw std::runtime_error("Missmatched array size");
+      std::copy(v_row.begin(), v_row.end(), vd_begin);
+      vd_begin += vsize[1];
+    }
+
+    return reinterpret_cast<PyObject *>(varray);
+}
+
+#else // fallback if we don't have numpy: copy every element of the given vector
+
+template<typename Numeric>
+PyObject* get_array(const std::vector<Numeric>& v)
+{
+    PyObject* list = PyList_New(v.size());
+    for(size_t i = 0; i < v.size(); ++i) {
+        PyList_SetItem(list, i, PyFloat_FromDouble(v.at(i)));
+    }
+    return list;
+}
+
+#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 = 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(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_plot, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    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>(),
+                  const long fig_number=0)
+{
+  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 = detail::get_2darray(x);
+  PyObject *yarray = detail::get_2darray(y);
+  PyObject *zarray = detail::get_2darray(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();
+  PyDict_SetItemString(kwargs, "rstride", PyInt_FromLong(1));
+  PyDict_SetItemString(kwargs, "cstride", PyInt_FromLong(1));
+
+  PyObject *python_colormap_coolwarm = PyObject_GetAttrString(
+      detail::_interpreter::get().s_python_colormap, "coolwarm");
+
+  PyDict_SetItemString(kwargs, "cmap", python_colormap_coolwarm);
+
+  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 *fig_args = PyTuple_New(1);
+  PyObject* fig = nullptr;
+  PyTuple_SetItem(fig_args, 0, PyLong_FromLong(fig_number));
+  PyObject *fig_exists =
+    PyObject_CallObject(
+    detail::_interpreter::get().s_python_function_fignum_exists, fig_args);
+  if (!PyObject_IsTrue(fig_exists)) {
+    fig = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+      detail::_interpreter::get().s_python_empty_tuple);
+  } else {
+    fig = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+      fig_args);
+  }
+  Py_DECREF(fig_exists);
+  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 *plot_surface = PyObject_GetAttrString(axis, "plot_surface");
+  if (!plot_surface) throw std::runtime_error("No surface");
+  Py_INCREF(plot_surface);
+  PyObject *res = PyObject_Call(plot_surface, args, kwargs);
+  if (!res) throw std::runtime_error("failed surface");
+  Py_DECREF(plot_surface);
+
+  Py_DECREF(axis);
+  Py_DECREF(args);
+  Py_DECREF(kwargs);
+  if (res) Py_DECREF(res);
+}
+
+template <typename Numeric>
+void contour(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 = {})
+{
+  detail::_interpreter::get();
+
+  // using numpy arrays
+  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);
+  PyTuple_SetItem(args, 0, xarray);
+  PyTuple_SetItem(args, 1, yarray);
+  PyTuple_SetItem(args, 2, zarray);
+
+  // Build up the kw args.
+  PyObject *kwargs = PyDict_New();
+
+  PyObject *python_colormap_coolwarm = PyObject_GetAttrString(
+      detail::_interpreter::get().s_python_colormap, "coolwarm");
+
+  PyDict_SetItemString(kwargs, "cmap", python_colormap_coolwarm);
+
+  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_contour, args, kwargs);
+  if (!res)
+    throw std::runtime_error("failed contour");
+
+  Py_DECREF(args);
+  Py_DECREF(kwargs);
+  if (res) Py_DECREF(res);
+}
+
+template <typename Numeric>
+void spy(const std::vector<::std::vector<Numeric>> &x,
+         const double markersize = -1,  // -1 for default matplotlib size
+         const std::map<std::string, std::string> &keywords = {})
+{
+  detail::_interpreter::get();
+
+  PyObject *xarray = detail::get_2darray(x);
+
+  PyObject *kwargs = PyDict_New();
+  if (markersize != -1) {
+    PyDict_SetItemString(kwargs, "markersize", PyFloat_FromDouble(markersize));
+  }
+  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 *plot_args = PyTuple_New(1);
+  PyTuple_SetItem(plot_args, 0, xarray);
+
+  PyObject *res = PyObject_Call(
+      detail::_interpreter::get().s_python_function_spy, plot_args, kwargs);
+
+  Py_DECREF(plot_args);
+  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>(),
+                  const long fig_number=0)
+{
+  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_args = PyTuple_New(1);
+  PyObject* fig = nullptr;
+  PyTuple_SetItem(fig_args, 0, PyLong_FromLong(fig_number));
+  PyObject *fig_exists =
+    PyObject_CallObject(detail::_interpreter::get().s_python_function_fignum_exists, fig_args);
+  if (!PyObject_IsTrue(fig_exists)) {
+    fig = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+      detail::_interpreter::get().s_python_empty_tuple);
+  } else {
+    fig = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+      fig_args);
+  }
+  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 = 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 (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_stem, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if (res)
+        Py_DECREF(res);
+
+    return res;
+}
+
+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 = detail::get_array(x);
+    PyObject* y1array = detail::get_array(y1);
+    PyObject* y2array = detail::get_array(y2);
+
+    // construct positional args
+    PyObject* args = PyTuple_New(3);
+    PyTuple_SetItem(args, 0, xarray);
+    PyTuple_SetItem(args, 1, y1array);
+    PyTuple_SetItem(args, 2, y2array);
+
+    // 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_fill_between, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+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* 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);
+
+    PyTuple_SetItem(plot_args, 0, yarray);
+
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_hist, plot_args, kwargs);
+
+
+    Py_DECREF(plot_args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    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 NumericX, typename NumericY, typename NumericColors>
+    bool scatter_colored(const std::vector<NumericX>& x,
+                 const std::vector<NumericY>& y,
+                 const std::vector<NumericColors>& colors,
+                 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* colors_array = detail::get_array(colors);
+
+        PyObject* kwargs = PyDict_New();
+        PyDict_SetItemString(kwargs, "s", PyLong_FromLong(s));
+        PyDict_SetItemString(kwargs, "c", colors_array);
+
+        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 NumericX, typename NumericY, typename NumericZ>
+bool scatter(const std::vector<NumericX>& x,
+             const std::vector<NumericY>& y,
+             const std::vector<NumericZ>& z,
+             const double s=1.0, // The marker size in points**2
+             const std::map<std::string, std::string> & keywords = {},
+             const long fig_number=0) {
+  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_args = PyTuple_New(1);
+  PyObject* fig = nullptr;
+  PyTuple_SetItem(fig_args, 0, PyLong_FromLong(fig_number));
+  PyObject *fig_exists =
+    PyObject_CallObject(detail::_interpreter::get().s_python_function_fignum_exists, fig_args);
+  if (!PyObject_IsTrue(fig_exists)) {
+    fig = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+      detail::_interpreter::get().s_python_empty_tuple);
+  } else {
+    fig = PyObject_CallObject(detail::_interpreter::get().s_python_function_figure,
+      fig_args);
+  }
+  Py_DECREF(fig_exists);
+  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, "scatter");
+  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);
+  Py_DECREF(fig);
+  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);
+}
+
+
+template<typename Numeric>
+bool barh(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 = { }) {
+    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_barh, plot_args, kwargs);
+
+    Py_DECREF(plot_args);
+    Py_DECREF(kwargs);
+    if (res) Py_DECREF(res);
+
+    return res;
+}
+
+
+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)
+{
+    detail::_interpreter::get();
+
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* kwargs = PyDict_New();
+    PyDict_SetItemString(kwargs, "label", PyString_FromString(label.c_str()));
+    PyDict_SetItemString(kwargs, "bins", PyLong_FromLong(bins));
+    PyDict_SetItemString(kwargs, "color", PyString_FromString(color.c_str()));
+    PyDict_SetItemString(kwargs, "alpha", PyFloat_FromDouble(alpha));
+
+
+    PyObject* plot_args = PyTuple_New(1);
+    PyTuple_SetItem(plot_args, 0, yarray);
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_hist, plot_args, kwargs);
+
+    Py_DECREF(plot_args);
+    Py_DECREF(kwargs);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool plot(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
+{
+    assert(x.size() == y.size());
+
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* pystring = PyString_FromString(s.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_CallObject(detail::_interpreter::get().s_python_function_plot, plot_args);
+
+    Py_DECREF(plot_args);
+    if(res) Py_DECREF(res);
+
+    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 = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+    PyObject* zarray = detail::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, typename NumericZ, typename NumericU, typename NumericW, typename NumericV>
+bool quiver(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::vector<NumericZ>& z, const std::vector<NumericU>& u, const std::vector<NumericW>& w, const std::vector<NumericV>& v, const std::map<std::string, std::string>& keywords = {})
+{
+  //set up 3d axes stuff
+  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 sizes match up
+  assert(x.size() == y.size() && x.size() == u.size() && u.size() == w.size() && x.size() == z.size() && x.size() == v.size() && u.size() == v.size());
+
+  //set up parameters
+  detail::_interpreter::get();
+
+  PyObject* xarray = detail::get_array(x);
+  PyObject* yarray = detail::get_array(y);
+  PyObject* zarray = detail::get_array(z);
+  PyObject* uarray = detail::get_array(u);
+  PyObject* warray = detail::get_array(w);
+  PyObject* varray = detail::get_array(v);
+
+  PyObject* plot_args = PyTuple_New(6);
+  PyTuple_SetItem(plot_args, 0, xarray);
+  PyTuple_SetItem(plot_args, 1, yarray);
+  PyTuple_SetItem(plot_args, 2, zarray);
+  PyTuple_SetItem(plot_args, 3, uarray);
+  PyTuple_SetItem(plot_args, 4, warray);
+  PyTuple_SetItem(plot_args, 5, varray);
+
+  // 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()));
+  }
+    
+  //get figure gca to enable 3d projection
+  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);
+  
+  //plot our boys bravely, plot them strongly, plot them with a wink and clap
+  PyObject *plot3 = PyObject_GetAttrString(axis, "quiver");
+  if (!plot3) throw std::runtime_error("No 3D line plot");
+  Py_INCREF(plot3);
+  PyObject* res = PyObject_Call(
+          plot3, plot_args, kwargs);
+  if (!res) throw std::runtime_error("Failed 3D plot");
+  Py_DECREF(plot3);
+  Py_DECREF(axis);
+  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());
+
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* pystring = PyString_FromString(s.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_CallObject(
+            detail::_interpreter::get().s_python_function_stem, plot_args);
+
+    Py_DECREF(plot_args);
+    if (res)
+        Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool semilogx(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
+{
+    assert(x.size() == y.size());
+
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* pystring = PyString_FromString(s.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_CallObject(detail::_interpreter::get().s_python_function_semilogx, plot_args);
+
+    Py_DECREF(plot_args);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool semilogy(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
+{
+    assert(x.size() == y.size());
+
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* pystring = PyString_FromString(s.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_CallObject(detail::_interpreter::get().s_python_function_semilogy, plot_args);
+
+    Py_DECREF(plot_args);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool loglog(const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& s = "")
+{
+    assert(x.size() == y.size());
+
+    detail::_interpreter::get();
+
+    PyObject* xarray = detail::get_array(x);
+    PyObject* yarray = detail::get_array(y);
+
+    PyObject* pystring = PyString_FromString(s.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_CallObject(detail::_interpreter::get().s_python_function_loglog, plot_args);
+
+    Py_DECREF(plot_args);
+    if(res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+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());
+
+    detail::_interpreter::get();
+
+    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 *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_errorbar, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+
+    if (res)
+        Py_DECREF(res);
+    else
+        throw std::runtime_error("Call to errorbar() failed.");
+
+    return res;
+}
+
+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 = detail::get_array(y);
+
+    PyObject* pystring = PyString_FromString(format.c_str());
+
+    PyObject* plot_args = PyTuple_New(2);
+
+    PyTuple_SetItem(plot_args, 0, yarray);
+    PyTuple_SetItem(plot_args, 1, pystring);
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_plot, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+    if (res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool named_plot(const std::string& name, const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& format = "")
+{
+    detail::_interpreter::get();
+
+    PyObject* kwargs = PyDict_New();
+    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) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool named_semilogx(const std::string& name, const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& format = "")
+{
+    detail::_interpreter::get();
+
+    PyObject* kwargs = PyDict_New();
+    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_semilogx, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+    if (res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool named_semilogy(const std::string& name, const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& format = "")
+{
+    detail::_interpreter::get();
+
+    PyObject* kwargs = PyDict_New();
+    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_semilogy, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+    if (res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename NumericX, typename NumericY>
+bool named_loglog(const std::string& name, const std::vector<NumericX>& x, const std::vector<NumericY>& y, const std::string& format = "")
+{
+    detail::_interpreter::get();
+
+    PyObject* kwargs = PyDict_New();
+    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_loglog, plot_args, kwargs);
+
+    Py_DECREF(kwargs);
+    Py_DECREF(plot_args);
+    if (res) Py_DECREF(res);
+
+    return res;
+}
+
+template<typename Numeric>
+bool plot(const std::vector<Numeric>& y, const std::string& format = "")
+{
+    std::vector<Numeric> x(y.size());
+    for(size_t i=0; i<x.size(); ++i) x.at(i) = i;
+    return plot(x,y,format);
+}
+
+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());
+    for (size_t i = 0; i < x.size(); ++i) x.at(i) = i;
+    return stem(x, y, format);
+}
+
+template<typename Numeric>
+void text(Numeric x, Numeric y, const std::string& s = "")
+{
+    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>
+inline void set_aspect(Numeric ratio)
+{
+    detail::_interpreter::get();
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(ratio));
+    PyObject* kwargs = PyDict_New();
+
+    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 *set_aspect = PyObject_GetAttrString(ax, "set_aspect");
+    if (!set_aspect) throw std::runtime_error("Attribute set_aspect not found.");
+    Py_INCREF(set_aspect);
+
+    PyObject *res = PyObject_Call(set_aspect, args, kwargs);
+    if (!res) throw std::runtime_error("Call to set_aspect() failed.");
+    Py_DECREF(set_aspect);
+
+    Py_DECREF(ax);
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+}
+
+inline void set_aspect_equal()
+{
+    // expect ratio == "equal". Leaving error handling to matplotlib.
+    detail::_interpreter::get();
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyString_FromString("equal"));
+    PyObject* kwargs = PyDict_New();
+
+    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 *set_aspect = PyObject_GetAttrString(ax, "set_aspect");
+    if (!set_aspect) throw std::runtime_error("Attribute set_aspect not found.");
+    Py_INCREF(set_aspect);
+
+    PyObject *res = PyObject_Call(set_aspect, args, kwargs);
+    if (!res) throw std::runtime_error("Call to set_aspect() failed.");
+    Py_DECREF(set_aspect);
+
+    Py_DECREF(ax);
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+}
+
+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));
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, list);
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_ylim, args);
+    if(!res) throw std::runtime_error("Call to ylim() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+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));
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, list);
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_xlim, args);
+    if(!res) throw std::runtime_error("Call to xlim() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+
+inline std::array<double, 2> xlim()
+{
+    PyObject* args = PyTuple_New(0);
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_xlim, args);
+
+    if(!res) throw std::runtime_error("Call to xlim() failed.");
+
+    Py_DECREF(res);
+
+    PyObject* left = PyTuple_GetItem(res,0);
+    PyObject* right = PyTuple_GetItem(res,1);
+    return { PyFloat_AsDouble(left), PyFloat_AsDouble(right) };
+}
+
+
+inline std::array<double, 2> ylim()
+{
+    PyObject* args = PyTuple_New(0);
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_ylim, args);
+
+    if(!res) throw std::runtime_error("Call to ylim() failed.");
+
+    Py_DECREF(res);
+
+    PyObject* left = PyTuple_GetItem(res,0);
+    PyObject* right = PyTuple_GetItem(res,1);
+    return { PyFloat_AsDouble(left), PyFloat_AsDouble(right) };
+}
+
+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));
+    PyTuple_SetItem(args, 1, PyFloat_FromDouble(ncols));
+    PyTuple_SetItem(args, 2, PyFloat_FromDouble(plot_number));
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_subplot, args);
+    if(!res) throw std::runtime_error("Call to subplot() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+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* 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);
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_axis, args);
+    if(!res) throw std::runtime_error("Call to title() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+inline void axhline(double y, double xmin = 0., double xmax = 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(y));
+    PyTuple_SetItem(args, 1, PyFloat_FromDouble(xmin));
+    PyTuple_SetItem(args, 2, PyFloat_FromDouble(xmax));
+
+    // 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_axhline, args, kwargs);
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+
+    if(res) Py_DECREF(res);
+}
+
+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 (auto 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* 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, 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* 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);
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, pyflag);
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_grid, args);
+    if(!res) throw std::runtime_error("Call to grid() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+inline void show(const bool block = true)
+{
+    detail::_interpreter::get();
+
+    PyObject* res;
+    if(block)
+    {
+        res = PyObject_CallObject(
+                detail::_interpreter::get().s_python_function_show,
+                detail::_interpreter::get().s_python_empty_tuple);
+    }
+    else
+    {
+        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);
+    }
+
+
+    if (!res) throw std::runtime_error("Call to show() failed.");
+
+    Py_DECREF(res);
+}
+
+inline void close()
+{
+    detail::_interpreter::get();
+
+    PyObject* res = PyObject_CallObject(
+            detail::_interpreter::get().s_python_function_close,
+            detail::_interpreter::get().s_python_empty_tuple);
+
+    if (!res) throw std::runtime_error("Call to close() failed.");
+
+    Py_DECREF(res);
+}
+
+inline void xkcd() {
+    detail::_interpreter::get();
+
+    PyObject* res;
+    PyObject *kwargs = PyDict_New();
+
+    res = PyObject_Call(detail::_interpreter::get().s_python_function_xkcd,
+            detail::_interpreter::get().s_python_empty_tuple, kwargs);
+
+    Py_DECREF(kwargs);
+
+    if (!res)
+        throw std::runtime_error("Call to show() failed.");
+
+    Py_DECREF(res);
+}
+
+inline void draw()
+{
+    detail::_interpreter::get();
+
+    PyObject* res = PyObject_CallObject(
+        detail::_interpreter::get().s_python_function_draw,
+        detail::_interpreter::get().s_python_empty_tuple);
+
+    if (!res) throw std::runtime_error("Call to draw() failed.");
+
+    Py_DECREF(res);
+}
+
+template<typename Numeric>
+inline void pause(Numeric interval)
+{
+    detail::_interpreter::get();
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, PyFloat_FromDouble(interval));
+
+    PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_pause, args);
+    if(!res) throw std::runtime_error("Call to pause() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(res);
+}
+
+inline void save(const std::string& filename, const int dpi=0)
+{
+    detail::_interpreter::get();
+
+    PyObject* pyfilename = PyString_FromString(filename.c_str());
+
+    PyObject* args = PyTuple_New(1);
+    PyTuple_SetItem(args, 0, pyfilename);
+
+    PyObject* kwargs = PyDict_New();
+
+    if(dpi > 0)
+    {
+        PyDict_SetItemString(kwargs, "dpi", PyLong_FromLong(dpi));
+    }
+
+    PyObject* res = PyObject_Call(detail::_interpreter::get().s_python_function_save, args, kwargs);
+    if (!res) throw std::runtime_error("Call to save() failed.");
+
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    Py_DECREF(res);
+}
+
+inline void rcparams(const std::map<std::string, std::string>& keywords = {}) {
+    detail::_interpreter::get();
+    PyObject* args = PyTuple_New(0);
+    PyObject* kwargs = PyDict_New();
+    for (auto it = keywords.begin(); it != keywords.end(); ++it) {
+        if ("text.usetex" == it->first)
+          PyDict_SetItemString(kwargs, it->first.c_str(), PyLong_FromLong(std::stoi(it->second.c_str())));
+        else PyDict_SetItemString(kwargs, it->first.c_str(), PyString_FromString(it->second.c_str()));
+    }
+    
+    PyObject * update = PyObject_GetAttrString(detail::_interpreter::get().s_python_function_rcparams, "update");
+    PyObject * res = PyObject_Call(update, args, kwargs);
+    if(!res) throw std::runtime_error("Call to rcParams.update() failed.");
+    Py_DECREF(args);
+    Py_DECREF(kwargs);
+    Py_DECREF(update);
+    Py_DECREF(res);
+}
+
+inline void clf() {
+    detail::_interpreter::get();
+
+    PyObject *res = PyObject_CallObject(
+        detail::_interpreter::get().s_python_function_clf,
+        detail::_interpreter::get().s_python_empty_tuple);
+
+    if (!res) throw std::runtime_error("Call to clf() failed.");
+
+    Py_DECREF(res);
+}
+
+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);
+
+    if (!res) throw std::runtime_error("Call to ion() failed.");
+
+    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);
+
+    if (!res) throw std::runtime_error("Call to tight_layout() failed.");
+
+    Py_DECREF(res);
+}
+
+// Support for variadic plot() and initializer lists:
+
+namespace detail {
+
+template<typename T>
+using is_function = typename std::is_function<std::remove_pointer<std::remove_reference<T>>>::type;
+
+template<bool obj, typename T>
+struct is_callable_impl;
+
+template<typename T>
+struct is_callable_impl<false, T>
+{
+    typedef is_function<T> type;
+}; // a non-object is callable iff it is a function
+
+template<typename T>
+struct is_callable_impl<true, T>
+{
+    struct Fallback { void operator()(); };
+    struct Derived : T, Fallback { };
+
+    template<typename U, U> struct Check;
+
+    template<typename U>
+    static std::true_type test( ... ); // use a variadic function to make sure (1) it accepts everything and (2) its always the worst match
+
+    template<typename U>
+    static std::false_type test( Check<void(Fallback::*)(), &U::operator()>* );
+
+public:
+    typedef decltype(test<Derived>(nullptr)) type;
+    typedef decltype(&Fallback::operator()) dtype;
+    static constexpr bool value = type::value;
+}; // an object is callable iff it defines operator()
+
+template<typename T>
+struct is_callable
+{
+    // dispatch to is_callable_impl<true, T> or is_callable_impl<false, T> depending on whether T is of class type or not
+    typedef typename is_callable_impl<std::is_class<T>::value, T>::type type;
+};
+
+template<typename IsYDataCallable>
+struct plot_impl { };
+
+template<>
+struct plot_impl<std::false_type>
+{
+    template<typename IterableX, typename IterableY>
+    bool operator()(const IterableX& x, const IterableY& y, const std::string& format)
+    {
+        detail::_interpreter::get();
+
+        // 2-phase lookup for distance, begin, end
+        using std::distance;
+        using std::begin;
+        using std::end;
+
+        auto xs = distance(begin(x), end(x));
+        auto ys = distance(begin(y), end(y));
+        assert(xs == ys && "x and y data must have the same number of elements!");
+
+        PyObject* xlist = PyList_New(xs);
+        PyObject* ylist = PyList_New(ys);
+        PyObject* pystring = PyString_FromString(format.c_str());
+
+        auto itx = begin(x), ity = begin(y);
+        for(size_t i = 0; i < xs; ++i) {
+            PyList_SetItem(xlist, i, PyFloat_FromDouble(*itx++));
+            PyList_SetItem(ylist, i, PyFloat_FromDouble(*ity++));
+        }
+
+        PyObject* plot_args = PyTuple_New(3);
+        PyTuple_SetItem(plot_args, 0, xlist);
+        PyTuple_SetItem(plot_args, 1, ylist);
+        PyTuple_SetItem(plot_args, 2, pystring);
+
+        PyObject* res = PyObject_CallObject(detail::_interpreter::get().s_python_function_plot, plot_args);
+
+        Py_DECREF(plot_args);
+        if(res) Py_DECREF(res);
+
+        return res;
+    }
+};
+
+template<>
+struct plot_impl<std::true_type>
+{
+    template<typename Iterable, typename Callable>
+    bool operator()(const Iterable& ticks, const Callable& f, const std::string& format)
+    {
+        if(begin(ticks) == end(ticks)) return true;
+
+        // We could use additional meta-programming to deduce the correct element type of y,
+        // but all values have to be convertible to double anyways
+        std::vector<double> y;
+        for(auto x : ticks) y.push_back(f(x));
+        return plot_impl<std::false_type>()(ticks,y,format);
+    }
+};
+
+} // end namespace detail
+
+// recursion stop for the above
+template<typename... Args>
+bool plot() { return true; }
+
+template<typename A, typename B, typename... Args>
+bool plot(const A& a, const B& b, const std::string& format, Args... args)
+{
+    return detail::plot_impl<typename detail::is_callable<B>::type>()(a,b,format) && plot(args...);
+}
+
+/*
+ * This group of plot() functions is needed to support initializer lists, i.e. calling
+ *    plot( {1,2,3,4} )
+ */
+inline bool plot(const std::vector<double>& x, const std::vector<double>& y, const std::string& format = "") {
+    return plot<double,double>(x,y,format);
+}
+
+inline bool plot(const std::vector<double>& y, const std::string& format = "") {
+    return plot<double>(y,format);
+}
+
+inline bool plot(const std::vector<double>& x, const std::vector<double>& y, const std::map<std::string, std::string>& keywords) {
+    return plot<double>(x,y,keywords);
+}
+
+/*
+ * 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

+ 158 - 0
Gravity_filtering/main.cpp

@@ -0,0 +1,158 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <ceres/ceres.h>
+#include <Eigen/Dense>
+#include <matplotlibcpp.h>
+#include <iomanip>
+namespace plt = matplotlibcpp;
+
+// 定义残差结构体
+struct SineResidual {
+    SineResidual(double x, double y) : x_(x), y_(y) {}
+    
+    template <typename T>
+    bool operator()(const T* const params, T* residual) const {
+        // params[0]: amplitude (a)
+        // params[1]: frequency (b)
+        // params[2]: phase (c)
+        // params[3]: offset (d)
+        residual[0] = y_ - (params[0] * sin(params[1] * x_ + params[2]) + params[3]);
+        return true;
+    }
+    
+private:
+    const double x_;
+    const double y_;
+};
+
+int main() {
+    // 读取数据 实际
+    std::ifstream file1("../f_data.csv");
+    if (!file1.is_open()) {
+        std::cerr << "Failed to open data file" << std::endl;
+        return 1;
+    }
+
+    std::vector<double> x_data, y_actual, x_ideal,y_ideal;
+    std::string line;
+    
+    // 跳过标题行
+    std::getline(file1, line);
+    
+    while (std::getline(file1, line)) {
+        double x, y;
+        if (sscanf(line.c_str(), "%lf,%lf", &x, &y) == 2) {
+            x_data.push_back(x);
+            y_actual.push_back(y);
+        }
+    }
+    file1.close();
+
+    // 读取数据 理想
+    std::ifstream file2("../f(1)_data.csv");
+    if (!file2.is_open()) {
+        std::cerr << "Failed to open data file" << std::endl;
+        return 1;
+    }
+
+    // 跳过标题行
+    std::getline(file2, line);
+
+    while (std::getline(file2, line)) {
+        double x, y;
+        if (sscanf(line.c_str(), "%lf,%lf", &x, &y) == 2) {
+            x_ideal.push_back(x);
+            y_ideal.push_back(y);
+        }
+    }
+    file2.close();
+
+    // 初始参数猜测 [a, b, c, d]
+    double params[2] = {5.0,1.0 };
+    
+    // 构建优化问题
+    ceres::Problem problem;
+    for (size_t i = 0; i < x_data.size(); ++i) {
+        problem.AddResidualBlock( //向问题中添加误差项
+        //使用自动求导,模板参数:误差类型、输出维度、输入维度
+            new ceres::AutoDiffCostFunction<SineResidual, 1, 2>(
+                new SineResidual(x_data[i], y_actual[i])),
+            nullptr,
+            params
+        );
+    }
+
+    // 配置求解器
+    ceres::Solver::Options options;
+    //options.linear_solver_type = ceres::DENSE_QR;
+    options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
+    options.minimizer_progress_to_stdout = false;
+    
+    // 运行求解器
+    ceres::Solver::Summary summary; //优化信息
+    ceres::Solve(options, &problem, &summary); //开始优化
+    
+    // 输出结果
+//    std::cout << summary.BriefReport() << "\n";
+//    std::cout << "Final params: a=" << params[0] << " b=" << params[1]
+//              << " c=" << params[2] << " d=" << params[3] << "\n";
+
+    // 生成拟合曲线
+    std::vector<double> y_fit;
+    for (double x : x_data) {
+        y_fit.push_back(params[0] * sin(params[1] * x + params[2]) + params[3]);
+    }
+
+//    // 生成理想曲线 (y = sin(x))
+//    for (double x : x_data) {
+//        y_ideal.push_back(sin(x));
+//    }
+
+    // 计算绝对差距和百分比差异
+    double total_abs_diff = 0.0;
+    double total_percent_diff = 0.0;
+    int count = y_fit.size();
+
+    for (size_t i = 0; i < count; ++i) {
+        double abs_diff = std::abs(y_fit[i] - y_ideal[i]);
+        double percent_diff = (abs_diff / std::abs(y_ideal[i])) * 100;
+
+        // 处理y_ideal接近0时的特殊情况
+        if (std::abs(y_ideal[i]) < 1e-6) {
+            percent_diff = 0;  // 避免除以零
+        }
+
+        total_abs_diff += abs_diff;
+        total_percent_diff += percent_diff;
+    }
+
+    double mean_abs_diff = total_abs_diff / count;
+    double mean_percent_diff = total_percent_diff / count;
+
+    // 将评估指标转换为字符串
+    std::ostringstream stats_text;
+    stats_text << "Mean Absolute Difference: " << std::fixed << std::setprecision(4) << mean_abs_diff << "\n"
+               << "Mean Percentage Difference: " << std::fixed << std::setprecision(2) << mean_percent_diff << "%";
+
+    // 可视化
+    plt::figure_size(1200, 800);
+    plt::scatter(x_data, y_actual, 10, {{"color", "red"}, {"label", "Actual Data"}});
+    plt::plot(x_data, y_fit, {{"color", "green"}, {"label", "Fitted Curve"}});
+    plt::plot(x_ideal, y_ideal, {{"color", "black"}, {"label", "Ideal Sine"}});
+
+    // 在图像上添加评估指标文本
+//    plt::text(0.02, 0.95, stats_text.str(),
+//              {{"color", "blue"},
+//               {"fontsize", "12"},
+//               {"bbox", "dict(facecolor='white', alpha=0.7, edgecolor='none')"}});
+    plt::text(0.01, -0.5, stats_text.str());
+    plt::xlabel("x");
+    plt::ylabel("y");
+    plt::title("Sine Curve Fitting using Ceres");
+    plt::legend();
+    plt::save("sine_fit.png");
+    plt::show();
+
+    return 0;
+}