Best practice to Share histogram between c++ and python

[MaxSac]

I usually create large histograms in C++ that I want to analyze and visualize in Python.
I have not found any documentation on what is the best practice to share a histogram between C++ and Python.
Is there by any chance a solution that is hard to find? @HDembinski

One possibility I could think of would be to recreate the histogram based on two files, one containing the axis information and the other containing the histogram values.

#include <boost/archive/text_oarchive.hpp>
#include <boost/histogram.hpp>
#include <boost/histogram/serialization.hpp>
#include <fstream>
#include <iostream>
#include <sstream>

namespace bh = boost::histogram;

int main()
{
    auto hist = bh::make_histogram(bh::axis::regular<>(10, 0.0, 1.0, "radius"),
        bh::axis::circular<>(10, 0.f, 2.f * M_PI, "phi"),
        bh::axis::circular<>(10, 0.f, 1.f * M_PI, "theta"));

    hist(bh::weight(0.1), 0.1, 0.f * M_PI, 1 * M_PI);
    hist(bh::weight(0.2), 0.5, 3.f * M_PI, 2 * M_PI);
    hist(bh::weight(0.3), 0.9, 2.f * M_PI, 3 * M_PI);

     {
        std::ofstream ofs("histogram.csv");
        for (auto&& x : indexed(hist))
            ofs << x.index(0) << "," << x.index(1) << "," << x.index(2) << ","
                << *x << "\n";
    }

    {
        std::ofstream ofs("metadata.txt");
        for (int i = 0; i < 3; ++i) 
            ofs << hist.axis(i) << "\n";
    }

    return 0;
}

I can imagine that there are much better solutions, but I haven't found them yet.

[HDembinski]

Hi, this use case is not covered in the docs, sorry.

If you have access to pybind11 and numpy on the machine where you want to generate your file in C++, you can call the C++ serialization code from boost-histogram directly to generate a pickle file, which you can then unpickle in Python normally.

The C++ code which is needed for Python's pickle protocol is here:
https://github.com/scikit-hep/boost-histogram/blob/develop/include/bh_python/make_pickle.hpp

You need to include this file and run this code (beware, this code was not actually tested):

py::tuple tup;
tuple_oarchive oa{tup};
oa << hist; // corresponds to `hist` in your example code

// call into Python libs to pickle `tup`
auto open = py::module_::import("gzip").attr("open");
auto dump = py::module_::import("pickle").attr("dump");
auto f = open("myfile.pkl.gz", "w");
dump(tup, f);
f.attr("close")();

This should generate a gzip'ed file with the binary contents of your histogram. Saving the data in this way is much more efficient than generating an ASCII file and the file will be quite small, too.

I think we should add an example to this repo which demonstrates this.

Apart from this solution, I was discussing with @henryiii and others at some point to make a separate repo with shared serialization code for Boost.Histogram (this lib) and boost-histogram (the Python wrapper), but we never got to do it. It shouldn't be difficult, however. We could use the portable_binary archive from https://github.com/USCiLab/cereal to generate a common binary representation of the histogram. cereal is source compatible with Boost.Serialization and Boost.Histogram already supports Boost.Serialization. We could also use Boost.Serialization directly, but it does not have support for portable binary.

That being said, the only caveat of the solution I show here is that you need to have pybind11 and Python libs installed on the machine where you want to write the file in C++. If that's not an issue, it is a great solution.

[HDembinski]

On second thought, this won't work right away since you create a static histogram here, while boost-histogram in Python uses the dynamic version, where the axes can be set at run-time. Furthermore, boost-histogram uses a special metadata type. The code in make_pickle.hpp was not written with this in mind, but I think it could be specialized in such a way, that writing static C++ histograms is possible.

[MaxSac]

If I modify the code to

#include <boost/histogram.hpp>
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>

#include <bh_python/make_pickle.hpp>

namespace bh = boost::histogram;

int main()
{
    using reg = bh::axis::regular<>;
    using cat = bh::axis::circular<>;
    using variant = bh::axis::variant<reg, cat>;
    std::vector<variant> axes;
    axes.emplace_back(reg(10, 0.0, 1.0, "radius"));
    axes.emplace_back(cat(10, 0.f, 1.f * M_PI, "theta"));
    axes.emplace_back(cat(10, 0.f, 2.f * M_PI, "phi"));

    auto hist = bh::make_histogram(std::move(axes));

    hist(bh::weight(0.1), 0.1, 0.f * M_PI, 1 * M_PI);
    hist(bh::weight(0.2), 0.5, 3.f * M_PI, 2 * M_PI);
    hist(bh::weight(0.3), 0.9, 2.f * M_PI, 3 * M_PI);

    py::tuple tup;
    tuple_oarchive oa{tup};
    oa << hist; // corresponds to `hist` in your example code

    // call into Python libs to pickle `tup`
    auto open = py::module_::import("gzip").attr("open");
    auto dump = py::module_::import("pickle").attr("dump");
    auto f = open("myfile.pkl.gz", "w");
    dump(tup, f);
    f.attr("close")();

    return 0;
}

and add pybind11::module as well as Python::Python to my target_link_libraries I get a segmentation fault which produces the following back trace in gdb

#0  _PyInterpreterState_GET () at ../Include/internal/pycore_pystate.h:122
#1  maybe_freelist_pop (size=1) at ../Objects/tupleobject.c:1177
#2  tuple_alloc (size=1) at ../Objects/tupleobject.c:43
#3  0x00007ffff77dad68 in PyTuple_New (size=1) at ../Objects/tupleobject.c:75
#4  PyTuple_New (size=1) at ../Objects/tupleobject.c:69
#5  0x000055555557765c in pybind11::tuple::tuple<unsigned long, 0> (
    this=0x7fffffffd9d0, size=1) at /usr/local/include/pybind11/pytypes.h:2066
#6  0x00005555555789d3 in pybind11::make_tuple<(pybind11::return_value_policy)1, pybind11::object const&> () at /usr/local/include/pybind11/cast.h:1402
#7  0x0000555555571e09 in tuple_oarchive::operator<< (this=0x7fffffffdac0, obj=...)
    at /usr/include/bh_python/make_pickle.hpp:175
#8  0x0000555555571f05 in tuple_oarchive::operator<< (this=0x7fffffffdac0, obj=...)
    at /usr/include/bh_python/make_pickle.hpp:184
#9  0x0000555555581b55 in tuple_oarchive::operator<< <unsigned int> (
    this=0x7fffffffdac0, t=@0x7fffffffda94: 0)
    at /usr/include/bh_python/make_pickle.hpp:159
#10 0x00005555555794f3 in tuple_oarchive::operator<< <boost::histogram::histogram<std::vector<boost::histogram::axis::variant<boost::histogram::axis::regular<double, boost::use_default, boost::use_default, boost::use_default>, boost::histogram::axis::regular<double, boost::histogram::axis::transform::id, boost::use_default, boost::histogram::axis::option::bitset<6u> > >, std::allocator<boost::histogram::axis::variant<boost::histogram::axis::regular<double, boost::use_default, boost::use_default, boost::use_default>, boost::histogram::axis::regular<double, boost::histogram::axis::transform::id, boost::use_default, boost::histogram::axis::option::bitset<6u> > > > >, boost::histogram::unlimited_storage<std::allocator<char> > > > (this=0x7fffffffdac0, t=...)
    at /usr/include/bh_python/make_pickle.hpp:168
#11 0x000055555555db0d in main () at ./hist.cpp:29

That is very difficult for me to follow the error. 😕 A minimal working example in the docs would definitely help.

[MaxSac]

@HDembinski do you have any hints? :)

[HDembinski]

You need to initialize the Python interpreter before you can call the Python API. See
https://pybind11.readthedocs.io/en/stable/advanced/embedding.html
You need to add this line:

    py::scoped_interpreter guard{}; // start the interpreter and keep it alive

The rest looks good.

[MaxSac]

that solves the problem. thank you very much

[MaxSac]

Is there a practical solution to recreate the histogram in Python from the read tuples? I have tried

import boost_histogram as bh
import pickle
import gzip

with gzip.open("build/myfile.pkl.gz", "rb") as f:
    h = pickle.load(f)


print(h)
# (0, 3, 0, 0, 0, 0, 0, 10, 'radius', 0.0, 1.0, 0, 0, 1, 0, 0, 10, 'theta', 0.0, 3.141592653589793, 0, 0, 1, 0, 0, 10, 'phi', 0.0, 6.283185307179586, 0, 5, 1452, array([0., 0., 0., ..., 0., 0., 0.]))

but I find no python constructor which is capable to produce an object out of this in python :/

[HDembinski]

Interesting, I am a bit surprised that this does not work, but the pickle file apparently did not store in this case that the stored object is a boost histogram. It thinks the stored object is a Python tuple. There is probably a more elegant way around this, but you can make it work by creating an empty Histogram in Python and then call its __setstate__ method with this tuple.

Try this (I haven't run it myself)

import boost_histogram as bh
import pickle
import gzip

with gzip.open("build/myfile.pkl.gz", "rb") as f:
    h_tuple = pickle.load(f)

h = bh.Histogram()
h.__setstate__(h_tuple)

# h should now be your C++ histogram

[HDembinski]

Turns out that this does not work, because @henryiii wraps the actual C++ histogram into a Python shell.