Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Patch utils: tentative v0.3. #11

Merged
merged 15 commits into from
Feb 27, 2023
54 changes: 54 additions & 0 deletions Makefile
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
PIP := $(shell command -v pip3 2> /dev/null || command which pip 2> /dev/null)
PYTHON := $(shell command -v python3 2> /dev/null || command which python 2> /dev/null)
PYTEST := $(shell command -v pytest 2> /dev/null)

.PHONY: install dev-install tests doc watchdoc servedoc lint typeannot coverage

pipcheck:
ifndef PIP
$(error "Ensure pip or pip3 are in your PATH")
endif
@echo Using pip: $(PIP)

pythoncheck:
ifndef PYTHON
$(error "Ensure python or python3 are in your PATH")
endif
@echo Using python: $(PYTHON)

pytestcheck:
ifndef PYTEST
$(error "Ensure pytest is in your PATH")
endif
@echo Using pytest: $(PYTEST)

install:
make pipcheck
$(PIP) install -r requirements.txt && $(PIP) install .

dev-install:
make pipcheck
$(PIP) install -r requirements-dev.txt && $(PIP) install -e .

tests:
make pytestcheck
$(PYTEST) tests

doc:
cd docs && rm -rf build && sphinx-apidoc -f -M -o source/ ../curvelops && make html && cd -
# Add -P to sphinx-apidoc to include private files

watchdoc:
while inotifywait -q -r curvelops/ -e create,delete,modify; do { make doc; }; done

servedoc:
$(PYTHON) -m http.server --directory docs/build/html/

lint:
flake8 docs/ curvelops/ tests/

typeannot:
mypy curvelops/

coverage:
coverage run -m pytest && coverage xml && coverage html && $(PYTHON) -m http.server --directory htmlcov/
141 changes: 78 additions & 63 deletions cpp/fdct2d_wrapper.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,87 +9,99 @@
*/

#include "fdct_wrapping.hpp"
#include "fdct_wrapping_inline.hpp"
#include "fdct_wrapping_inc.hpp"
#include <pybind11/pybind11.h>
#include "fdct_wrapping_inline.hpp"
#include <iostream>
#include <pybind11/complex.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/complex.h>
#include <iostream>
#include <sstream>
namespace py = pybind11;
using namespace std;
using namespace fdct_wrapping_ns;

py::tuple fdct2d_param_wrap(int m, int n, int nbscales, int nbangles_coarse, int ac)
py::tuple
fdct2d_param_wrap(int m, int n, int nbscales, int nbangles_coarse, int ac)
{
// Almost sure this function creates a copy, but it's ok since the outputs are small
// Almost sure this function creates a copy, but it's ok since the outputs
// are small
vector<vector<double>> sx, sy;
vector<vector<double>> fx, fy;
vector<vector<int>> nx, ny;
fdct_wrapping_param(m, n, nbscales, nbangles_coarse, ac, sx, sy, fx, fy, nx, ny);
fdct_wrapping_param(
m, n, nbscales, nbangles_coarse, ac, sx, sy, fx, fy, nx, ny);
return py::make_tuple(sx, sy, fx, fy, nx, ny);
}

vector<vector<py::array_t<cpx>>> fdct2d_forward_wrap(int nbscales, int nbangles_coarse, int ac, py::array_t<cpx> x)
vector<vector<py::array_t<cpx>>>
fdct2d_forward_wrap(int nbscales,
int nbangles_coarse,
int ac,
py::array_t<cpx> x)
{
// Our wrapper takes a NumPy array, but ``fdct_wrapping`` requires a CpxNumMat
// input (which will be accessed read-only). So we must create CpxNumMat ``xmat``
// which will "mirror" our input ``x`` in a no-copy fashion.
// We also need to output ``c`` whose conversion to a Python list of lists of
// CpxNumMat will be handled by pybind11. The vector -> list casting is automatic
// in pybind11, whereas the CpxNumMat -> py::array_t<cpx> casting is inside our function.
// Our wrapper takes a NumPy array, but ``fdct_wrapping`` requires a
// CpxNumMat input (which will be accessed read-only). So we must create
// CpxNumMat ``xmat`` which will "mirror" our input ``x`` in a no-copy
// fashion. We also need to output ``c`` whose conversion to a Python list
// of lists of CpxNumMat will be handled by pybind11. The vector -> list
// casting is automatic in pybind11, whereas the CpxNumMat ->
// py::array_t<cpx> casting is inside our function.
CpxNumMat xmat;
vector<vector<CpxNumMat>> cmat;

// Responsibly access py::array with possible casting to complex. See:
// https://stackoverflow.com/questions/42645228/cast-numpy-array-to-from-custom-c-matrix-class-using-pybind11
// Note: CurveLab uses Fortran-style indexing, so we must transpose the input array. We do this
// Note: CurveLab uses Fortran-style indexing, so we must transpose the
// input array. We do this
// by simply reading it as a Fortran array
auto buf = py::array_t<cpx, py::array::f_style | py::array::forcecast>::ensure(x);
auto buf =
py::array_t<cpx, py::array::f_style | py::array::forcecast>::ensure(x);
if (!buf)
throw std::runtime_error("x array buffer is empty. If you're calling from Python this should not happen!");
throw std::runtime_error("x array buffer is empty. If you're calling "
"from Python this should not happen!");
if (buf.ndim() != 2)
throw std::runtime_error("x.ndims != 2");

// We don't to initialize ``x(m, n)`` because this allocates an array on the heap!
// We don't to initialize ``x(m, n)`` because this allocates an array on
// the heap!
xmat._m = buf.shape()[0];
xmat._n = buf.shape()[1];
xmat._data = (cpx *)buf.data(); // Put our Python array buffer pointer as the CpxNumMat data
xmat._data =
(cpx*)buf
.data(); // Put our Python array buffer pointer as the CpxNumMat data

// Call our forward function with all the right types
fdct_wrapping(xmat._m, xmat._n, nbscales, nbangles_coarse, ac, xmat, cmat);

// Clear the structure as if it had never existed...
// xmat didn't allocate any data, so we make sure it doesn't deallocate any on the way out
// xmat didn't allocate any data, so we make sure it doesn't deallocate any
// on the way out
xmat._m = xmat._n = 0;
xmat._data = NULL;

vector<vector<py::array_t<cpx>>> c;
// Expand ``c`` to fit the scales
c.resize(cmat.size());
for (size_t i = 0; i < cmat.size(); i++)
{
for (size_t i = 0; i < cmat.size(); i++) {
// Now we expand each scale to fit the angles
c[i].resize(cmat[i].size());
for (size_t j = 0; j < cmat[i].size(); j++)
{
// Create capsule linked to `cmat[i][j]._data` to track its lifetime
for (size_t j = 0; j < cmat[i].size(); j++) {
// Create capsule linked to `cmat[i][j]._data` to track its
// lifetime
// https://stackoverflow.com/questions/44659924/returning-numpy-arrays-via-pybind11
py::capsule free_when_done(
cmat[i][j].data(),
[](void *cpx_ptr)
{
cpx *cpx_arr = reinterpret_cast<cpx *>(cpx_ptr);
delete[] cpx_arr;
});
py::capsule free_when_done(cmat[i][j].data(), [](void* cpx_ptr) {
cpx* cpx_arr = reinterpret_cast<cpx*>(cpx_ptr);
delete[] cpx_arr;
});

py::array c_arr(
{cmat[i][j]._n, cmat[i][j]._m}, // Shape
{sizeof(cpx) * cmat[i][j]._m, // Strides (in bytes) of the underlying data array
sizeof(cpx)},
cmat[i][j].data(), // Data pointer
free_when_done);
{ cmat[i][j]._n, cmat[i][j]._m }, // Shape
{ sizeof(cpx) * cmat[i][j]._m, // Strides (in bytes) of the
// underlying data array
sizeof(cpx) },
cmat[i][j].data(), // Data pointer
free_when_done);

c[i][j] = c_arr;
cmat[i][j]._m = cmat[i][j]._n = 0;
Expand All @@ -99,8 +111,13 @@ vector<vector<py::array_t<cpx>>> fdct2d_forward_wrap(int nbscales, int nbangles_
return c;
}

py::array_t<cpx> fdct2d_inverse_wrap(int m, int n, int nbscales, int nbangles_coarse, int ac,
vector<vector<py::array_t<cpx>>> c)
py::array_t<cpx>
fdct2d_inverse_wrap(int m,
int n,
int nbscales,
int nbangles_coarse,
int ac,
vector<vector<py::array_t<cpx>>> c)
{
// Similarly to the forward wrapper, we create ``cmat`` and ``xmat`` to use
// as dummy input and output arrays.
Expand All @@ -114,46 +131,38 @@ py::array_t<cpx> fdct2d_inverse_wrap(int m, int n, int nbscales, int nbangles_co
// We copy the ``c`` "structure" onto a ``cmat`` "structure"
// Start by expanding the first index of ``cmat`` to fit all scales
cmat.resize(c.size());
for (i = 0; i < c.size(); i++)
{
for (i = 0; i < c.size(); i++) {
// Now we expand each scale to fit all angles for that scale
cmat[i].resize(c[i].size());
for (j = 0; j < c[i].size(); j++)
{
for (j = 0; j < c[i].size(); j++) {
// Now we must copy the structure over to ``cmat``
py::buffer_info buf = c[i][j].request();
cmat[i][j]._m = buf.shape[1];
cmat[i][j]._n = buf.shape[0];
cmat[i][j]._data = static_cast<cpx *>(buf.ptr);
cmat[i][j]._data = static_cast<cpx*>(buf.ptr);
}
}
// No bounds checking is made inside this, so if ``c`` (or equivalently ``cmat``)
// are not compatible with the other parameters, this function WILL segfault
// No bounds checking is made inside this, so if ``c`` (or equivalently
// ``cmat``) are not compatible with the other parameters, this function
// WILL segfault
// TODO: Optionally sanitize this by calling ``fdct2d_param_wrap``
ifdct_wrapping(m, n, nbscales, nbangles_coarse, ac, cmat, xmat);

// Clear input structure without deallocating
for (i = 0; i < c.size(); i++)
for (j = 0; j < c[i].size(); j++)
{
for (j = 0; j < c[i].size(); j++) {
cmat[i][j]._m = cmat[i][j]._n = 0;
cmat[i][j]._data = NULL;
}

py::capsule free_when_done(
xmat.data(),
[](void *cpx_ptr)
{
cpx *cpx_arr = reinterpret_cast<cpx *>(cpx_ptr);
delete[] cpx_arr;
});
py::capsule free_when_done(xmat.data(), [](void* cpx_ptr) {
cpx* cpx_arr = reinterpret_cast<cpx*>(cpx_ptr);
delete[] cpx_arr;
});

// Create output array
py::array x({m, n},
{sizeof(cpx),
sizeof(cpx) * m},
xmat.data(),
free_when_done);
py::array x(
{ m, n }, { sizeof(cpx), sizeof(cpx) * m }, xmat.data(), free_when_done);

// Clear output structure without deallocating
xmat._m = xmat._n = 0;
Expand All @@ -165,10 +174,16 @@ py::array_t<cpx> fdct2d_inverse_wrap(int m, int n, int nbscales, int nbangles_co
PYBIND11_MODULE(fdct2d_wrapper, m)
{
m.doc() = "FDCT2D pybind11 wrapper";
m.def("fdct2d_param_wrap", &fdct2d_param_wrap, "Parameters for 2D FDCT",
m.def("fdct2d_param_wrap",
&fdct2d_param_wrap,
"Parameters for 2D FDCT",
py::return_value_policy::take_ownership);
m.def("fdct2d_forward_wrap", &fdct2d_forward_wrap, "2D Forward FDCT",
m.def("fdct2d_forward_wrap",
&fdct2d_forward_wrap,
"2D Forward FDCT",
py::return_value_policy::take_ownership);
m.def("fdct2d_inverse_wrap", &fdct2d_inverse_wrap, "2D Inverse FDCT",
m.def("fdct2d_inverse_wrap",
&fdct2d_inverse_wrap,
"2D Inverse FDCT",
py::return_value_policy::take_ownership);
}
Loading