[discussion] pybind11 bindings of Octave into Python land (as part of Pythonic or another official Octave placement)

Since Pytave a lot of time has passed, and now pybind11 makes it much easier to produce bindings of C++ interfaces and C++ classes into Python (extensively used by PyTorch). Even auto-generators exist (and hopefully should be easier to use than SWIG).

I haven't tried auto-generators, but I've tried binding norm function, and it has worked quite easily.

pybind11 in-memory (and zero-copy especially) data exchange should be much more efficient than the temp file approach of oct2py. This should also enable more practical use of legacy MatLab code from within Python.

Some questions for fully-blown bindings:

• For some reason LD_PRELOAD was needed, and even just LD_LIBRARY_PATH was not sufficient
• Would be nice to have zero-copy convertors between Octave and NumPy:
  • should be possible, since NumPy also supports column-major
  • if the user really wants, it should be possible to override existing NumPy array's layout to column-major and still zero-copy-convert it to Octave's Matrix (would result in transposed shape, but sometimes it's okay)
  • Matrix/NDArray should allow a constructor which takes in an existing pointer
• Ideally, can bind all Octave's C++ classes to Python - either with manual pybind11 wrappers (ChatGPT could also help in this work) or with generated ones (https://github.com/RosettaCommons/binder https://pybind11.readthedocs.io/en/stable/compiling.html#generating-binding-code-automatically) - seems such generators are now in much better shape than e.g. 10 years ago
• Could represent certain known functions as methods on the Interpreter object, and also allow general Octave code eval
• Ideally both dense and sparse tensors could be mapped to NumPy/SciPy
• Should have a generic function for converting octave_value_list from/to Python/NumPy - probably Pythonic already has such conversions implemented
• While classes like Matrix could be bound to Python as OctaveMatrix, and it can still be useful, it needs to be decided whether the user should do the conversions NumPy<>OctaveMatrix themselves or if the bindings should themselves do auto-conversions at least for some NumPy/Python types. Pybind11 can have support for both modes, it seems (see TypeCaster thingies).

# sudo apt-get install octave octave-dev
# pip install pybind11
# make octave_pybind11.so
# python3 octave_pybind11_test.py

PYTHON = python3
PYBIND11_INCLUDES := $(shell $(PYTHON) -m pybind11 --includes)
CXXFLAGS := -O3 -Wall -shared -std=c++11 -fPIC $(PYBIND11_INCLUDES)

# SUFFIX := $(shell python3-config --extension-suffix)
# PYTHON_CFLAGS := $(shell python3-config --cflags)
# PYTHON_LIBS := $(shell python3-config --libs)

octave_pybind11.so: octave_pybind11.cpp
        mkoctfile -v --link-stand-alone $(CXXFLAGS) $< -o $@

clean:
        -rm -f octave_pybind11.so

test:
        LD_PRELOAD="/usr/lib/x86_64-linux-gnu/octave/8.4.0/liboctinterp.so.11 /usr/lib/x86_64-linux-gnu/octave/8.4.0/liboctave.so.10" $(PYTHON) octave_pybind11_test.py
        #
        #$(PYTHON) octave_pybind11_test.py # ImportError: liboctinterp.so.11: cannot open shared object file: No such file or directory
        #LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu/octave/8.4.0/ $(PYTHON) octave_pybind11_test.py # error: /usr/lib/x86_64-linux-gnu/octave/8.4.0/oct/x86_64-pc-linux-gnu/__init_gnuplot__.oct: failed to load | Incompatible version or missing dependency? | /usr/lib/x86_64-linux-gnu/octave/8.4.0/oct/x86_64-pc-linux-gnu/__init_gnuplot__.oct: undefined symbol: _ZTI17octave_base_value | error: called from | /usr/lib/x86_64-linux-gnu/octave/8.4.0/oct/x86_64-pc-linux-gnu/PKG_ADD at line 4 column 5 | as if liboctinterp.so.11 did not get loaded

// octave_pybind11.cpp

#include <memory>
#include <string>
#include <iostream>

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

#include <octave/oct.h>
#include <octave/octave.h>
#include <octave/parse.h>
#include <octave/interpreter.h>
#include <octave/builtin-defun-decls.h>
#include <octave/Matrix.h>

class OctavePybind11Interpreter {
private:
    octave::interpreter interpreter;
public:
    OctavePybind11Interpreter()
    {
        // Inhibit reading history file by calling: interpreter.initialize_history (false);
        // Set custom load path here if you wish by calling: interpreter.initialize_load_path (false);
        // Perform final initialization of interpreter, including executing commands from startup files by calling: interpreter.initialize(); if (! interpreter.initialized ()) { std::cerr << "Octave interpreter initialization failed!" << std::endl; exit (status); }
        // You may skip this step if you don't need to do anything between reading the startup files and telling the interpreter that you are ready to execute commands. Tell the interpreter that we're ready to execute commands:
        int status = interpreter.execute ();
        if (status != 0)
            throw pybind11::value_error("creating embedded Octave interpreter failed!");
    }

    double norm(pybind11::array_t<double, pybind11::array::f_style | pybind11::array::forcecast> numpy_array) const
    {
        if (numpy_array.ndim() != 2)
            throw pybind11::value_error("Input NumPy array must be 2-dimensional.");

        try
        {
                        size_t rows = numpy_array.shape(0);
            size_t cols = numpy_array.shape(1);
            size_t itemsize = numpy_array.itemsize();
            const double* src_ptr = numpy_array.data();

            Matrix octave_matrix(rows, cols); // global namespace? octave::Matrix doesn't work
            double* dst_ptr = (double*)octave_matrix.data();
            std::memcpy(dst_ptr, src_ptr, itemsize * rows * cols);

            octave_value_list in;
            in(0) = octave_matrix;
            in(1) = "fro";
            octave_value_list out = octave::Fnorm(in, 1); //octave_value_list out = octave::feval ("gcd", in, 1);
            double norm_of_the_matrix = out(0).double_value ();

            return norm_of_the_matrix;
        }
        catch (const octave::exit_exception& ex)
        {
            throw pybind11::value_error("Octave interpreter exited with status = " + std::to_string(ex.exit_status()));
        }
        catch (const octave::execution_exception&)
        {
            throw pybind11::value_error("error encountered in Octave evaluator!");
        }
    }

};

PYBIND11_MODULE(octave_pybind11, m)
{
    pybind11::class_<OctavePybind11Interpreter, std::shared_ptr<OctavePybind11Interpreter>>(m, "OctavePybind11Interpreter")
        .def(pybind11::init<>()) // Default constructor
        .def("norm", &OctavePybind11Interpreter::norm)
        ;
}

# octave_pybind11_test.py

import numpy as np
import octave_pybind11

interp = octave_pybind11.OctavePybind11Interpreter()

myarr = np.array([[1.0, 2.0], [3.0, 4.0]], dtype = np.float64)

print(np.linalg.norm(myarr, "fro"))
print(interp.norm(myarr))