[py] corrected bug in IntervalMatrix from Numpy

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()

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_<IntervalMatrix> interval_matrix(m, "IntervalMatrix");
@@ -79,20 +89,34 @@ void export_IntervalMatrix(py::module& m)
   .def(py::init<const IntervalMatrix&>())
   .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<Eigen::Dynamic, Eigen::Dynamic> 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<double>::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<Eigen::MatrixXd, 0, Strides>(
+      static_cast<double*>(info.ptr), info.shape[0], info.shape[1], strides);
 
-    ibex::Matrix m((int)info.shape[0], (int)info.shape[1], static_cast<double*>(info.ptr));
-    return IntervalMatrix(m);
+    return codac2::to_codac1(codac2::Matrix(map));
   }))
 
   .def("__getitem__", [](IntervalMatrix& s, size_t index) -> IntervalVector&

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_<R,C>(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>(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)
       { }
 

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>(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)
       { }
 

src/core/2/variables/codac2_Matrix.h

@@ -19,7 +19,7 @@ namespace codac2
 {
   using Eigen::Dynamic;
 
-  template<int R,int C>
+  template<int R=Dynamic,int C=Dynamic>
   class Matrix_ : public Eigen::Matrix<double,R,C>
   {
     public:
@@ -28,27 +28,27 @@ namespace codac2
         : Eigen::Matrix<double,R,C>()
       { }
 
-      Matrix_(size_t nb_rows, size_t nb_cols)
+      Matrix_(int nb_rows, int nb_cols)
         : Eigen::Matrix<double,R,C>(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<double,R,C>(nb_rows, nb_cols)
       {
         assert(R == Dynamic || R == (int)nb_rows);
         assert(C == Dynamic || C == (int)nb_cols);
         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_<R,C>(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>(R, C, values)
       { }
 
@@ -184,6 +184,17 @@ namespace codac2
     os << ")";
     return os;
   }
+
+  #include "ibex_Matrix.h"
+  template<int R=Dynamic,int C=Dynamic>
+  ibex::Matrix to_codac1(const Matrix_<R,C>& 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<int R,int C>
   Matrix_<R,C> floor(const Matrix_<R,C>& 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<std::initializer_list<double>> l)
+        : Matrix_<>(l)
+      { }
+
+      template<int R,int C>
+      explicit Matrix(const Matrix_<R,C>& v)
+        : Matrix_<>(v)
+      { }
+  };
 } // namespace codac
 
 #endif

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>(n,1,values)
       { }
 
-      explicit Vector_(double values[])
+      explicit Vector_(const double values[])
         : Matrix_<N,1>(N,1,values)
       { }