matplotlibcpp.h

#pragma once

// Python headers must be included before any system headers, since
// they define _POSIX_C_SOURCE
#include <Python.h>
#include <algorithm>
#include <array>
#include <chrono>
#include <datetime.h>
#include <functional>
#include <iostream>
#include <map>
#include <numeric>
#include <ranges>
#include <stdexcept>
#include <string>
#include <vector>

#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
#endif

#define CPP20 202002L

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

                // FIXME:
                // https://docs.python.org/3.12/c-api/init_config.html#init-config

#if PY_MAJOR_VERSION >= 3
#if PY_MINOR_VERSION >= 8

                // https://docs.python.org/3/c-api/init_config.html
                PyStatus status;

                PyConfig config;
                PyConfig_InitPythonConfig(&config);
                config.isolated = 1;

                // const char* dummy_args[]
                //     = {"Python", NULL}; // const is needed because literals
                //                          // must not be modified
                // char* const* argv = dummy_args;
                // char* const* argv = nullptr;
                // const char* argv_[] =  {"PyPlot", NULL};
                const char* const argv_[] =  {"PyPlot", NULL};
                // int argc = sizeof(dummy_args) / sizeof(dummy_args[0]) - 1;
                int argc_ = sizeof(argv_) / sizeof(argv_[0]) - 1;
                // int argc=1;

                /* Decode command line arguments.
                   Implicitly preinitialize Python (in isolated mode). */
                status = PyConfig_SetBytesArgv(&config, argc_, (char* const*)argv_);
                if (PyStatus_Exception(status)) {
                    // goto exception;
                    throw std::runtime_error(status.err_msg);
                }

                status = Py_InitializeFromConfig(&config);
                if (PyStatus_Exception(status)) {
                    throw std::runtime_error("Initialization from config failed");
                    // goto exception;
                }
                PyConfig_Clear(&config);

                // Py_RunMain();

#else

                wchar_t name[] = L"plotting";
                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;

                PySys_SetArgv(argc, const_cast<wchar_t**>(argv));
#endif
#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
    void backend(const std::string& name) { detail::s_backend = name; }

    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
        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

    // If a container (range) has ::data() that converts to ::value and it has
    // ::size() that is convertible to pointer differences, then it's
    // contiguous, in which case we need not copy the data.
    // https://stackoverflow.com/questions/42851957/contiguous-iterator-detection
    //
    // TODO: get rid of the plot(vector) specializations, they unnecessarily
    // copy the data.
    // TODO: provide contiguous and non-contiguous implementations for ranges,
    // using iterator_tag.
    // https://stackoverflow.com/questions/4688177/how-does-iterator-category-in-c-work
    // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p1474r0.pdf
    //

#if __cplusplus >= CPP20

    // Support for contiguous ranges, e.g. vector, span, etc. In this case
    // we can pass the data pointer directly to python, without copying.
    //
    // Non-contiguous (but iterable) arrays (e.g. lists), call a different plot,
    // one which has to copy the items into a contiguous C-style array before
    // passing them to python.
    template<std::ranges::contiguous_range ContainerX,
             std::ranges::contiguous_range ContainerY>
    bool plot(const ContainerX& x, const ContainerY& y,
              const std::string& fmt = "")
    {
        assert(y.size() % x.size() == 0
               && "length of y must be a multiple of length of x!");

        detail::_interpreter::get();

        NPY_TYPES xtype
            = detail::select_npy_type<typename ContainerX::value_type>::type;
        NPY_TYPES ytype
            = detail::select_npy_type<typename ContainerY::value_type>::type;

        npy_intp xsize = x.size();
        npy_intp yrows = xsize, ycols = y.size() / x.size();
        npy_intp ysize[] = {yrows, ycols}; // ysize[0] must equal xsize

        // We have 2 options to pass existing data buffers to PyObject:
        // PyArray_SimpleFromData() - by default creates C-style (row major)
        // array PyArray_New() - uses flags, so we can specify Fotran-style (col
        // major)
        //
        // In python,
        // a=np.array([[3,5], [1,4], [4,1], [5,3]])
        // is stored in row-major mode and has columns
        // a[:,0]=[3,1,4,5] and a[:,1]=[5,4,1,3]
        // Then,
        // plt.plot(a)
        // plots the columns of a. So in python the array is row-major, but
        // plotted columnwise.
        //
        // For dataframes (and time series), however, it is more natural to
        // assume that the data is stored contiguously in column-major mode,
        // i.e. double a[] = [3, 1, 4, 5
        //               5, 4, 1, 3];
        // We let python know that is the case, by specifying the
        // NPY_ARRAY_FARRAY flag. This rules out the usage of
        // PyArray_SimpleFromData().
        //
        // Of course, in the vector-only version of this plot() function all
        // this is irrelevant, because that plot is 1-dimensional anyway.
        //
        // If there are real-world applications that need to assume that the
        // C-data is in row-major mode, we should address that.
        //
        // TODO:
        // To that end, perhaps we can introduce C++ tags cmajor_tag and
        // rmajor_tag and define a custom concept from contiguous_range, by
        // further conditioning on the storage tags. The caller then specifies
        // the major storage mode when wrapping the C-style array into a range.
        PyObject* xarray
            = PyArray_New(&PyArray_Type, 1, &xsize, xtype, nullptr,
                          (void*)x.data(), 0, NPY_ARRAY_FARRAY, nullptr);
        PyObject* yarray = PyArray_New(&PyArray_Type, 2, ysize, ytype, nullptr,
                                       (void*)y.data(), 0, NPY_ARRAY_FARRAY,
                                       nullptr); // column major by design!

        PyObject* pystring = PyString_FromString(fmt.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<std::ranges::contiguous_range ContainerX,
             std::ranges::contiguous_range ContainerY>
    bool plot(const ContainerX& x, const ContainerY& y,
              const std::map<std::string, std::string>& keywords)
    {
        assert(y.size() % x.size() == 0
               && "length of y must be a multiple of length of x!");

        detail::_interpreter::get();

        NPY_TYPES xtype
            = detail::select_npy_type<typename ContainerX::value_type>::type;
        NPY_TYPES ytype
            = detail::select_npy_type<typename ContainerY::value_type>::type;

        npy_intp xsize = x.size();
        npy_intp yrows = xsize, ycols = y.size() / x.size();
        npy_intp ysize[] = {yrows, ycols}; // ysize[0] must equal xsize

        // Same comments as above.
        PyObject* xarray
            = PyArray_New(&PyArray_Type, 1, &xsize, xtype, nullptr,
                          (void*)x.data(), 0, NPY_ARRAY_FARRAY, nullptr);
        PyObject* yarray = PyArray_New(&PyArray_Type, 2, ysize, ytype, nullptr,
                                       (void*)y.data(), 0, NPY_ARRAY_FARRAY,
                                       nullptr); // column major by design!

        // 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;
    }

    template<std::ranges::contiguous_range ContainerY>
    bool plot(const ContainerY& y, const std::string& format = "")
    {
        std::vector<int> x(y.size());
        std::iota(x.begin(), x.end(), 0);
        return plot(x, y, format);
    }

    template<std::ranges::contiguous_range ContainerY>
    bool plot(const ContainerY& y,
              const std::map<std::string, std::string>& keywords)
    {
        std::vector<int> x(y.size());
        std::iota(x.begin(), x.end(), 0);
        return plot(x, y, keywords);
    }

#else

    // For the less fortunate ones who can't or won't use C++ 20, we stick with
    // vectors, since there's no ranges concept before that.

    template<typename Numeric>
    bool plot(const std::vector<Numeric>& x, const std::vector<Numeric>& y,
              const std::map<std::string, std::string>& keywords)
    {
        assert(y.size() % x.size() == 0
               && "length of y must be a multiple of length of x!");

        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;
    }

    template<typename NumericX, typename NumericY>
    bool plot(const std::vector<NumericX>& x, const std::vector<NumericY>& y,
              const std::string& s = "")
    {
        assert(y.size() % x.size() == 0
               && "length of y must be a multiple of length of x!");

        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 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);
    }

#endif

// 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
    {
        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

    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);
    }

    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;
    }

    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, 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 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);
    }

    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);
    }

    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;
    }

    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;
    }

    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);
    }

    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);
    }

    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<std::ranges::contiguous_range Labels>
    void legend(const Labels& labels, const std::map<std::string, std::string>& keywords={})
    {
        detail::_interpreter::get();

        PyObject* llist = PyList_New(labels.size());
        if(llist == nullptr)
            throw "Can't allocate labels list in legend() function.";

        // iterate over labels, hopefully not too many!
        size_t i=0;
        for (const auto& l : labels) {
            PyObject* str = PyString_FromString(l.c_str());
            PyList_SET_ITEM(llist, i++, 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* args = PyTuple_New(1);
        PyTuple_SetItem(args, 0, llist);

        PyObject* res = PyObject_Call(
            detail::_interpreter::get().s_python_function_legend,
            args, kwargs);
        // PyObject* res = PyObject_Call(
        //     detail::_interpreter::get().s_python_function_legend,
        //     detail::_interpreter::get().s_python_empty_tuple, kwargs);
        if(!res) throw std::runtime_error("Call to legend() failed.");

        Py_DECREF(kwargs);
        Py_DECREF(res);
    }

    template<typename Numeric>
    void 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);
    }

    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);
    }

    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)};
    }

    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>
    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
        // https://docs.python.org/3/c-api/dict.html
        // https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.xticks.html
        // https://matplotlib.org/stable/api/text_api.html#matplotlib.text.Text
        PyObject* kwargs = PyDict_New();
        for(std::map<std::string, std::string>::const_iterator it
                = keywords.begin();
            it != keywords.end(); ++it) {
            const auto& [key, val] = *it;
            try {
                // if it's numeric, add it as number, e.g. rotation = 45
                // https://matplotlib.org/stable/api/text_api.html#matplotlib.text.Text
                double x = std::stod(val);
                PyDict_SetItemString(kwargs, key.c_str(), PyFloat_FromDouble(x));
            }
            catch (const std::exception& e) {
                // if it wasn't a number, add it as string
                // TODO: allow for other types, e.g. bool
                // https://matplotlib.org/stable/api/text_api.html#matplotlib.text.Text
                PyDict_SetItemString(kwargs, key.c_str(),
                                     PyString_FromString(val.c_str()));

            }
            // 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>
    void xticks(const std::vector<Numeric>& ticks,
                const std::map<std::string, std::string>& keywords)
    {
        xticks(ticks, {}, keywords);
    }

    // options only, e.g. plt::xticks(rotation=20)
    void xticks(const std::map<std::string, std::string>& keywords)
    {
        detail::_interpreter::get();

        PyObject* args = PyTuple_New(0);

        // construct keyword args
        PyObject* kwargs = PyDict_New();
        for(const auto& [key, val] : keywords){
            try {
                // if it's numeric, add it as number, e.g. rotation = 45
                // https://matplotlib.org/stable/api/text_api.html#matplotlib.text.Text
                double x = std::stod(val);
                PyDict_SetItemString(kwargs, key.c_str(), PyFloat_FromDouble(x));
            }
            catch (const std::exception& e) {
                // if it wasn't a number, add it as string
                // TODO: allow for other types, e.g. bool
                // https://matplotlib.org/stable/api/text_api.html#matplotlib.text.Text
                PyDict_SetItemString(kwargs, key.c_str(),
                                     PyString_FromString(val.c_str()));
            }
        }

            // PyDict_SetItemString(kwargs, k.c_str(),
            //                      PyString_FromString(v.c_str()));

        auto* xticks = detail::_interpreter::get().s_python_function_xticks;
        PyObject* res = PyObject_Call(xticks, args, kwargs);
        // PyObject* res = PyObject_Call(
        //     detail::_interpreter::get().s_python_function_xticks, args, kwargs);

        Py_DECREF(kwargs);
        if(!res) throw std::runtime_error("Call to xticks() failed");

        Py_DECREF(res);
    }

    template<typename Numeric>
    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>
    void yticks(const std::vector<Numeric>& ticks,
                const std::map<std::string, std::string>& keywords)
    {
        yticks(ticks, {}, keywords);
    }

    template<typename Numeric>
    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>
    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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>
    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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);
    }

    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?
    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
        {};

#ifdef WITHOUT_NUMPY

        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::begin;
                    using std::distance;
                    using std::end;

                    auto xs = distance(begin(x), end(x));
                    auto ys = distance(begin(y), end(y));

                    assert(xs == ys
                           && "x and y must have the same number of elements!");

                    // No PyArray, must use lists, and they must have the same
                    // length.
                    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(decltype(xs) 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, 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;
                }
        };

#else

        template<>
        struct plot_impl<std::false_type>
        {
            template<typename IterableX, typename IterableY>
            bool operator()(const IterableX& x, const IterableY& y,
                            const std::string& format)
                {
                    // pyassert(PyArray_API, "NumPy needed");

                    detail::_interpreter::get();

                    // 2-phase lookup for distance, begin, end
                    using std::begin;
                    using std::distance;
                    using std::end;

                    auto xs = distance(begin(x), end(x));
                    auto ys = distance(begin(y), end(y));

                    assert(ys % xs == 0
                           && "length of y must be a multiple of length of x!");

                    typedef typename IterableX::value_type NumberX;
                    typedef typename IterableY::value_type NumberY;

                    NPY_TYPES xtype = detail::select_npy_type<NumberX>::type;
                    NPY_TYPES ytype = detail::select_npy_type<NumberY>::type;

                    npy_intp xsize = xs;
                    npy_intp yrows = xsize, ycols = ys / yrows;
                    npy_intp ysize[] = {yrows, ycols}; // ysize[0] must equal xsize

                    PyObject* xarray
                        = PyArray_New(&PyArray_Type, 1, &xsize, xtype, nullptr,
                                      nullptr, 0, NPY_ARRAY_FARRAY, nullptr);
                    PyObject* yarray = PyArray_New(
                        &PyArray_Type, 2, ysize, ytype, nullptr, nullptr, 0,
                        NPY_ARRAY_FARRAY, nullptr); // column major!
                    PyObject* pystring = PyString_FromString(format.c_str());

                    // fill the data
                    auto itx = begin(x), ity = begin(y);
                    NumberX* xdata = (NumberX*)PyArray_DATA((PyArrayObject*)xarray);
                    for(decltype(xs) i = 0; i < xs; ++i) xdata[i] = *itx++;

                    NumberY* ydata = (NumberY*)PyArray_DATA((PyArrayObject*)yarray);
                    for(decltype(ys) i = 0; i < ys; ++i) ydata[i] = *ity++;

                    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;
                }
        };

#endif

        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} )
     */
    bool plot(const std::vector<double>& x, const std::vector<double>& y,
              const std::string& format = "")
    {
#if __cplusplus >= CPP20
        return plot<std::vector<double>, std::vector<double>>(x, y, format);
#else
        return plot<double>(x, y, format);
#endif
    }

    bool plot(const std::vector<double>& y, const std::string& format = "")
    {
#if __cplusplus >= CPP20
        return plot<std::vector<double>>(y, format);
#else
        return plot<double>(y, format);
#endif
    }

    /*
     * This class allows dynamic plots, ie changing the plotted data without
     * clearing and re-plotting
     */
    class Plot
    {
    public:
        // default initialization with plot label, some data and format
        template<typename Numeric>
        Plot(const std::string& name, const std::vector<Numeric>& x,
             const std::vector<Numeric>& y, const std::string& format = "")
            {
                detail::_interpreter::get();

                assert(x.size() == y.size());

                PyObject* kwargs = PyDict_New();
                if(name != "")
                    PyDict_SetItemString(kwargs, "label",
                                         PyString_FromString(name.c_str()));

                PyObject* xarray = detail::get_array(x);
                PyObject* yarray = detail::get_array(y);

                PyObject* pystring = PyString_FromString(format.c_str());

                PyObject* plot_args = PyTuple_New(3);
                PyTuple_SetItem(plot_args, 0, xarray);
                PyTuple_SetItem(plot_args, 1, yarray);
                PyTuple_SetItem(plot_args, 2, pystring);

                PyObject* res = PyObject_Call(
                    detail::_interpreter::get().s_python_function_plot, plot_args,
                    kwargs);

                Py_DECREF(kwargs);
                Py_DECREF(plot_args);

                if(res) {
                    line = PyList_GetItem(res, 0);

                    if(line)
                        set_data_fct = PyObject_GetAttrString(line, "set_data");
                    else
                        Py_DECREF(line);
                    Py_DECREF(res);
                }
            }

        // shorter initialization with name or format only
        // basically calls line, = plot([], [])
        Plot(const std::string& name = "", const std::string& format = "")
            : Plot(name, std::vector<double>(), std::vector<double>(), format)
            {}

        template<typename Numeric>
        bool update(const std::vector<Numeric>& x,
                    const std::vector<Numeric>& y)
            {
                assert(x.size() == y.size());
                if(set_data_fct) {
                    PyObject* xarray = detail::get_array(x);
                    PyObject* yarray = detail::get_array(y);

                    PyObject* plot_args = PyTuple_New(2);
                    PyTuple_SetItem(plot_args, 0, xarray);
                    PyTuple_SetItem(plot_args, 1, yarray);

                    PyObject* res = PyObject_CallObject(set_data_fct, plot_args);
                    if(res) Py_DECREF(res);
                    return res;
                }
                return false;
            }

        // clears the plot but keep it available
        bool clear()
            {
                return update(std::vector<double>(), std::vector<double>());
            }

        // definitely remove this line
        void remove()
            {
                if(line) {
                    auto remove_fct = PyObject_GetAttrString(line, "remove");
                    PyObject* args = PyTuple_New(0);
                    PyObject* res = PyObject_CallObject(remove_fct, args);
                    if(res) Py_DECREF(res);
                }
                decref();
            }

        ~Plot() { decref(); }
    private:
        void decref()
            {
                if(line) Py_DECREF(line);
                if(set_data_fct) Py_DECREF(set_data_fct);
            }

        PyObject* line = nullptr;
        PyObject* set_data_fct = nullptr;
    };

} // end namespace matplotlibcpp