24.09., 9:00 - 11:00: Due to updates GitLab will be unavailable for some minutes between 09:00 and 11:00.

bindings_hints.hpp 8.65 KB
Newer Older
Nikola Dinev's avatar
Nikola Dinev committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
#include "DataContainer.h"
#include "Descriptors/VolumeDescriptor.h"
#include "LinearOperator.h"
#include "DescriptorUtils.h"

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

#include <functional>

namespace elsa
{
    template <typename Class>
    struct ClassHints {
    };

    namespace py = pybind11;

    /// wrapper for variadic functions
    template <typename T, std::size_t N, typename = std::make_index_sequence<N>>
    struct invokeBestCommonVariadic;

    template <typename T, std::size_t N, std::size_t... S>
    struct invokeBestCommonVariadic<T, N, std::index_sequence<S...>> {
        static auto exec(const py::args& args)
        {
            if (args.size() == N)
                return elsa::bestCommon(args[S].template cast<T>()...);

            if constexpr (N > 1) {
                return invokeBestCommonVariadic<T, N - 1>::exec(args);
            } else {
                throw(std::logic_error("Unsupported number of variadic arguments"));
            }
        };
    };

    template <typename data_t>
    LinearOperator<data_t> adjointHelper(const LinearOperator<data_t>& op)
    {
        return adjoint(op);
    }

    template <typename data_t>
    LinearOperator<data_t> leafHelper(const LinearOperator<data_t>& op)
    {
        return leaf(op);
    }

    // define global variables and functions in the module hints
    struct ModuleHints {
        static void addCustomFunctions(py::module& m)
        {
            m.def("adjoint", &adjointHelper<float>)
                .def("adjoint", &adjointHelper<double>)
                .def("adjoint", &adjointHelper<std::complex<float>>)
                .def("adjoint", &adjointHelper<std::complex<double>>);

            m.def("leaf", &leafHelper<float>)
                .def("leaf", &leafHelper<double>)
                .def("leaf", &leafHelper<std::complex<float>>)
                .def("leaf", &leafHelper<std::complex<double>>);

            m.def("bestCommon", (std::unique_ptr<DataDescriptor>(*)(
                                    const std::vector<const DataDescriptor*>&))(&bestCommon))
                .def("bestCommon", &invokeBestCommonVariadic<const DataDescriptor&, 10>::exec);
        }
    };

    template <typename data_t, typename type_, typename... options>
    void addOperatorsDc(py::class_<type_, options...>& c)
    {
        c.def(
             "__add__",
             [](const DataContainer<data_t>& self, const DataContainer<data_t>& other) {
                 return DataContainer<data_t>(self + other);
             },
             py::return_value_policy::move)
            .def(
                "__mul__",
                [](const DataContainer<data_t>& self, const DataContainer<data_t>& other) {
                    return DataContainer<data_t>(self * other);
                },
                py::return_value_policy::move)
            .def(
                "__sub__",
                [](const DataContainer<data_t>& self, const DataContainer<data_t>& other) {
                    return DataContainer<data_t>(self - other);
                },
                py::return_value_policy::move)
            .def(
                "__truediv__",
                [](const DataContainer<data_t>& self, const DataContainer<data_t>& other) {
                    return DataContainer<data_t>(self / other);
                },
                py::return_value_policy::move)
            // TODO: make the generator automatically generate this __setitem__ function
            .def("__setitem__", [](elsa::DataContainer<data_t>& dc, elsa::index_t i, data_t value) {
                dc[i] = value;
            });
    }

    template <typename data_t>
    struct DataContainerHints : public ClassHints<elsa::DataContainer<data_t>> {
        constexpr static std::tuple ignoreMethods = {
            "operator()", "begin", "cbegin", "end", "cend", "rbegin", "crbegin", "rend", "crend"};

        template <typename type_, typename... options>
        static void addCustomMethods(py::class_<type_, options...>& c)
        {
            addOperatorsDc<data_t>(c);
        }

        template <typename type_, typename... options>
        static void exposeBufferInfo(py::class_<type_, options...>& c)
        {
            c.def(py::init([](py::buffer b) {
                 py::buffer_info info = b.request();

                 if (info.format != py::format_descriptor<data_t>::format())
                     throw std::invalid_argument("Incompatible scalar types");

                 elsa::IndexVector_t coeffsPerDim(info.ndim);

                 ssize_t minStride = info.strides[0];
                 for (std::size_t i = 0; i < static_cast<std::size_t>(info.ndim); i++) {
                     if (info.strides[i] < minStride)
                         minStride = info.strides[i];

                     coeffsPerDim[static_cast<elsa::index_t>(i)] =
                         static_cast<elsa::index_t>(info.shape[i]);
                 }

                 if (static_cast<std::size_t>(minStride) / sizeof(data_t) != 1)
                     throw std::invalid_argument("Cannot convert strided buffer to DataContainer");

                 auto map = Eigen::Map<Eigen::Matrix<data_t, Eigen::Dynamic, 1>>(
                     static_cast<data_t*>(info.ptr), coeffsPerDim.prod());

                 elsa::VolumeDescriptor dd{coeffsPerDim};

                 return std::make_unique<elsa::DataContainer<data_t>>(dd, map);
             })).def_buffer([](elsa::DataContainer<data_t>& m) {
                std::vector<ssize_t> dims, strides;
                auto coeffsPerDim = m.getDataDescriptor().getNumberOfCoefficientsPerDimension();
                ssize_t combined = 1;
                for (int i = 0; i < coeffsPerDim.size(); i++) {
                    dims.push_back(coeffsPerDim[i]);
                    strides.push_back(combined * static_cast<ssize_t>(sizeof(data_t)));
                    combined *= coeffsPerDim[i];
                }
                return py::buffer_info(
                    &m[0], sizeof(data_t), py::format_descriptor<data_t>::format(),
                    m.getDataDescriptor().getNumberOfDimensions(), coeffsPerDim, strides);
            });
        }
    };

    template <typename data_t>
    struct LinearOperatorHints : public ClassHints<elsa::LinearOperator<data_t>> {

        template <typename type_, typename... options>
        static void addCustomMethods(py::class_<type_, options...>& c)
        {
            c.def(py::self + py::self).def(py::self * py::self);
        }
    };

    template <typename data_t>
    struct DataContainerComplexHints : public ClassHints<elsa::DataContainer<data_t>> {
        constexpr static std::tuple ignoreMethods = {
            "operator()", "begin", "cbegin", "end", "cend", "rbegin", "crbegin", "rend", "crend"};

        template <typename type_, typename... options>
        static void addCustomMethods(py::class_<type_, options...>& c)
        {
            addOperatorsDc<data_t>(c);
        }

        // TODO: pybind11 does not provide a default format_descriptor for complex types -> make a
        // custom one for this to work
        // CustomMethod<true, py::buffer_info, elsa::DataContainer<data_t>&> buffer =
        //     [](elsa::DataContainer<data_t>& m) {
        //         std::vector<ssize_t> dims, strides;
        //         auto coeffsPerDim = m.getDataDescriptor().getNumberOfCoefficientsPerDimension();
        //         ssize_t combined = 1;
        //         for (int i = 0; i < coeffsPerDim.size(); i++) {
        //             dims.push_back(coeffsPerDim[i]);
        //             strides.push_back(combined * static_cast<ssize_t>(sizeof(data_t)));
        //             combined *= coeffsPerDim[i];
        //         }
        //         return py::buffer_info(
        //             &m[0],          /* Pointer to buffer */
        //             sizeof(data_t), /* Size of one scalar */
        //             py::format_descriptor<data_t>::format(),
        //             /* Python struct-style format descriptor
        //              */
        //             m.getDataDescriptor().getNumberOfDimensions(), /* Number of dimensions */
        //             coeffsPerDim,                                  /* Buffer dimensions */
        //             strides);
        //     };
    };

    template struct DataContainerHints<float>;
    template struct DataContainerComplexHints<std::complex<float>>;
    template struct DataContainerHints<double>;
    template struct DataContainerComplexHints<std::complex<double>>;
    template struct DataContainerHints<index_t>;
    template struct LinearOperatorHints<float>;
    template struct LinearOperatorHints<double>;
    template struct LinearOperatorHints<std::complex<float>>;
    template struct LinearOperatorHints<std::complex<double>>;
} // namespace elsa