From c8fe9ec79476d06ce5f22e0ed3bf90bcae5975e8 Mon Sep 17 00:00:00 2001 From: SimonRohou Date: Sun, 29 Oct 2023 23:29:42 +0100 Subject: [PATCH] [py] corrected bug in IntervalMatrix from Numpy --- python/codac/tests/test_numpy.py | 56 ++++++++++++++++++ .../interval/codac_py_IntervalMatrix.cpp | 36 ++++++++++-- .../domains/interval/codac2_IntervalMatrix.h | 6 +- .../domains/interval/codac2_IntervalVector.h | 6 +- src/core/2/variables/codac2_Matrix.h | 57 ++++++++++++++++--- src/core/2/variables/codac2_Vector.h | 4 +- 6 files changed, 143 insertions(+), 22 deletions(-) create mode 100644 python/codac/tests/test_numpy.py diff --git a/python/codac/tests/test_numpy.py b/python/codac/tests/test_numpy.py new file mode 100644 index 00000000..99f2a7d0 --- /dev/null +++ b/python/codac/tests/test_numpy.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python + +# Codac tests +# --------------------------------------------------------------------------- +# \date 2023 +# \author Simon Rohou +# \copyright Copyright 2023 Codac Team +# \license This program is distributed under the terms of +# the GNU Lesser General Public License (LGPL). + +import unittest +import numpy as np +from codac import * + +class TestNumpyMatrices(unittest.TestCase): + + def test_numpytocodac(self): + + cdA=IntervalMatrix(2,6) + cdA[0][0]=Interval(1) + cdA[0][1]=Interval(3) + cdA[0][2]=Interval(5) + cdA[0][3]=Interval(7) + cdA[0][4]=Interval(9) + cdA[0][5]=Interval(11) + cdA[1][0]=Interval(2) + cdA[1][1]=Interval(4) + cdA[1][2]=Interval(6) + cdA[1][3]=Interval(8) + cdA[1][4]=Interval(10) + cdA[1][5]=Interval(12) + + npA=np.array([[1.,3.,5.,7.,9.,11.],[2.,4.,6.,8.,10.,12.]]) + self.assertEqual(IntervalMatrix(npA), cdA) + + def test_numpytocodac_withtranspose(self): + + cdA=IntervalMatrix(2,6) + cdA[0][0]=Interval(1) + cdA[0][1]=Interval(3) + cdA[0][2]=Interval(5) + cdA[0][3]=Interval(7) + cdA[0][4]=Interval(9) + cdA[0][5]=Interval(11) + cdA[1][0]=Interval(2) + cdA[1][1]=Interval(4) + cdA[1][2]=Interval(6) + cdA[1][3]=Interval(8) + cdA[1][4]=Interval(10) + cdA[1][5]=Interval(12) + + npA=np.array([[1.,2.],[3.,4.],[5.,6.],[7.,8.],[9.,10.],[11.,12.]]).T + self.assertEqual(IntervalMatrix(npA), cdA) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/python/src/core/domains/interval/codac_py_IntervalMatrix.cpp b/python/src/core/domains/interval/codac_py_IntervalMatrix.cpp index abac3d26..6892720c 100644 --- a/python/src/core/domains/interval/codac_py_IntervalMatrix.cpp +++ b/python/src/core/domains/interval/codac_py_IntervalMatrix.cpp @@ -69,6 +69,16 @@ string to_string(const IntervalMatrix& x) return ss.str(); } +// Inspired by Numpy documentation +// See: https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html + + #include "codac2_IntervalMatrix.h" + /* Bind MatrixXd (or some other Eigen type) to Python */ + typedef Eigen::MatrixXd Matrix; + + //typedef Eigen::Matrix::Scalar Scalar; + constexpr bool rowMajor = !true;//Eigen::Matrix::Flags & Eigen::RowMajorBit; + void export_IntervalMatrix(py::module& m) { py::class_ interval_matrix(m, "IntervalMatrix"); @@ -79,20 +89,34 @@ void export_IntervalMatrix(py::module& m) .def(py::init()) .def(py::init(&create_from_list)) - .def(py::init([](py::buffer b) // create for instance an IntervalMatrix from a NumPy matrix + .def(py::init([](py::buffer b) { + // Inspired by Numpy documentation + // See: https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html + + // Note: use raw data does not work because of strides, + // see for instance transpose operations on numpy matrices + + typedef Eigen::Stride Strides; + // Request a buffer descriptor from Python py::buffer_info info = b.request(); - // Some sanity checks... + // Some basic validation checks... if(info.format != py::format_descriptor::format()) - throw std::runtime_error("Incompatible format: expected a double array"); + throw std::runtime_error("Incompatible format: expected a double array!"); if(info.ndim != 2) - throw std::runtime_error("Incompatible buffer dimension"); + throw std::runtime_error("Incompatible buffer dimension!"); + + auto strides = Strides( + info.strides[rowMajor ? 0 : 1] / (py::ssize_t)sizeof(double), + info.strides[rowMajor ? 1 : 0] / (py::ssize_t)sizeof(double)); + + auto map = Eigen::Map( + static_cast(info.ptr), info.shape[0], info.shape[1], strides); - ibex::Matrix m((int)info.shape[0], (int)info.shape[1], static_cast(info.ptr)); - return IntervalMatrix(m); + return codac2::to_codac1(codac2::Matrix(map)); })) .def("__getitem__", [](IntervalMatrix& s, size_t index) -> IntervalVector& diff --git a/src/core/2/domains/interval/codac2_IntervalMatrix.h b/src/core/2/domains/interval/codac2_IntervalMatrix.h index fc2121e1..19a5038f 100644 --- a/src/core/2/domains/interval/codac2_IntervalMatrix.h +++ b/src/core/2/domains/interval/codac2_IntervalMatrix.h @@ -51,7 +51,7 @@ namespace codac2 init(x); } - explicit IntervalMatrix_(size_t nb_rows, size_t nb_cols, double bounds[][2]) + explicit IntervalMatrix_(size_t nb_rows, size_t nb_cols, const double bounds[][2]) : IntervalMatrix_(nb_rows, nb_cols) { size_t k = 0; @@ -67,7 +67,7 @@ namespace codac2 assert(k == this->size()); } - explicit IntervalMatrix_(double bounds[][2]) + explicit IntervalMatrix_(const double bounds[][2]) : IntervalMatrix_(R, C, bounds) { } @@ -596,7 +596,7 @@ namespace codac2 : IntervalMatrix_<>(nb_rows, nb_cols, x) { } - explicit IntervalMatrix(size_t nb_rows, size_t nb_cols, double bounds[][2]) + explicit IntervalMatrix(size_t nb_rows, size_t nb_cols, const double bounds[][2]) : IntervalMatrix_<>(nb_rows, nb_cols, bounds) { } diff --git a/src/core/2/domains/interval/codac2_IntervalVector.h b/src/core/2/domains/interval/codac2_IntervalVector.h index 1e3d1db5..61398726 100644 --- a/src/core/2/domains/interval/codac2_IntervalVector.h +++ b/src/core/2/domains/interval/codac2_IntervalVector.h @@ -68,11 +68,11 @@ namespace codac2 (*this)[i] = Interval(v[i]); } - explicit IntervalVector_(size_t n, double bounds[][2]) + explicit IntervalVector_(size_t n, const double bounds[][2]) : IntervalMatrix_(n,1,bounds) { } - explicit IntervalVector_(double bounds[][2]) + explicit IntervalVector_(const double bounds[][2]) : IntervalVector_(this->size(), bounds) { } @@ -278,7 +278,7 @@ namespace codac2 : IntervalVector_<>(v) { } - explicit IntervalVector(size_t n, double bounds[][2]) + explicit IntervalVector(size_t n, const double bounds[][2]) : IntervalVector_<>(n, bounds) { } diff --git a/src/core/2/variables/codac2_Matrix.h b/src/core/2/variables/codac2_Matrix.h index 3cf9c142..b3e8fa03 100644 --- a/src/core/2/variables/codac2_Matrix.h +++ b/src/core/2/variables/codac2_Matrix.h @@ -19,7 +19,7 @@ namespace codac2 { using Eigen::Dynamic; - template + template class Matrix_ : public Eigen::Matrix { public: @@ -28,14 +28,14 @@ namespace codac2 : Eigen::Matrix() { } - Matrix_(size_t nb_rows, size_t nb_cols) + Matrix_(int nb_rows, int nb_cols) : Eigen::Matrix(nb_rows, nb_cols) { assert(R == Dynamic || R == (int)nb_rows); assert(C == Dynamic || C == (int)nb_cols); } - Matrix_(size_t nb_rows, size_t nb_cols, double x) + Matrix_(int nb_rows, int nb_cols, double x) : Eigen::Matrix(nb_rows, nb_cols) { assert(R == Dynamic || R == (int)nb_rows); @@ -43,12 +43,12 @@ namespace codac2 init(x); } - explicit Matrix_(size_t nb_rows, size_t nb_cols, double values[]) + explicit Matrix_(int nb_rows, int nb_cols, const double values[]) : Matrix_(nb_rows, nb_cols) { - size_t k = 0; - for(size_t i = 0 ; i < nb_rows ; i++) - for(size_t j = 0 ; j < nb_cols ; j++) + int k = 0; + for(int i = 0 ; i < nb_rows ; i++) + for(int j = 0 ; j < nb_cols ; j++) { if(values == 0) // in case the user called Matrix_(r,c,0) and 0 is interpreted as NULL! (*this)(i,j) = 0.; @@ -58,7 +58,7 @@ namespace codac2 } } - explicit Matrix_(double values[]) + explicit Matrix_(const double values[]) : Matrix_(R, C, values) { } @@ -184,6 +184,17 @@ namespace codac2 os << ")"; return os; } + + #include "ibex_Matrix.h" + template + ibex::Matrix to_codac1(const Matrix_& x) + { + ibex::Matrix x_(x.rows(), x.cols()); + for(size_t i = 0 ; i < x.rows() ; i++) + for(size_t j = 0 ; j < x.cols() ; j++) + x_[i][j] = x(i,j); + return x_; + } template Matrix_ floor(const Matrix_& x) @@ -221,6 +232,36 @@ namespace codac2 return f; } + class Matrix : public Matrix_<> + { + public: + + explicit Matrix(size_t nb_rows, size_t nb_cols) + : Matrix_<>(nb_rows, nb_cols) + { } + + + explicit Matrix(size_t nb_rows, size_t nb_cols, double x) + : Matrix_<>(nb_rows, nb_cols, x) + { } + + explicit Matrix(size_t nb_rows, size_t nb_cols, const double values[]) + : Matrix_<>(nb_rows, nb_cols, values) + { } + + Matrix(const Matrix_<>& x) + : Matrix_<>(x) + { } + + Matrix(std::initializer_list> l) + : Matrix_<>(l) + { } + + template + explicit Matrix(const Matrix_& v) + : Matrix_<>(v) + { } + }; } // namespace codac #endif \ No newline at end of file diff --git a/src/core/2/variables/codac2_Vector.h b/src/core/2/variables/codac2_Vector.h index fe9bb847..7c30cd75 100644 --- a/src/core/2/variables/codac2_Vector.h +++ b/src/core/2/variables/codac2_Vector.h @@ -55,11 +55,11 @@ namespace codac2 static_assert(M == Dynamic || M == N); } - explicit Vector_(size_t n, double values[]) + explicit Vector_(size_t n, const double values[]) : Matrix_(n,1,values) { } - explicit Vector_(double values[]) + explicit Vector_(const double values[]) : Matrix_(N,1,values) { }