Fix error with column-major storage order.

Benjamin Chrétien · Benjamin Chrétien · commit b09390e3736f · 2015-04-15T22:51:13.000+09:00
diff --git a/src/wrap.cc b/src/wrap.cc
@@ -300,18 +300,28 @@ namespace roboptim
 		return;
 	      }
 
-	    // TODO: check storage order (NumPy expects C-style here, but Eigen uses ColMajor as default)
-	    // Moreover, RobOptim users can change the storage order.
 	    npy_intp inputSize =
 	      static_cast<npy_intp> (::roboptim::core::python::Function::inputSize ());
 	    npy_intp outputSize =
 	      static_cast<npy_intp> (::roboptim::core::python::Function::outputSize ());
 
+	    // Check storage order and map memory accordingly (PyArray_SimpleNewFromData
+	    // expects a row-major matrix).
+	    PyObject* jacobianNumpy = NULL;
 	    npy_intp sizes[2] = {outputSize, inputSize};
 
-	    PyObject* jacobianNumpy =
-	      PyArray_SimpleNewFromData (2, &sizes[0], NPY_DOUBLE,
-					 jacobian.data ());
+	    if (::roboptim::StorageOrder == Eigen::RowMajor)
+	      {
+		jacobianNumpy = PyArray_SimpleNewFromData (2, &sizes[0], NPY_DOUBLE,
+							   jacobian.data ());
+	      }
+	    else
+	      {
+		jacobianNumpy = PyArray_NewFromDescr (&PyArray_Type, PyArray_DescrFromType (PyArray_DOUBLE),
+						      2, sizes, NULL, jacobian.data (),
+						      NPY_WRITEABLE | NPY_F_CONTIGUOUS, NULL);
+	      }
+
 	    if (!jacobianNumpy)
 	      {
 		PyErr_SetString (PyExc_TypeError, "cannot convert result");
@@ -1176,7 +1186,7 @@ getName (PyObject*, PyObject* args)
 }
 
 static PyObject*
-getStorageOrder (PyObject*, PyObject* args)
+getStorageOrder (PyObject*, PyObject*)
 {
   char storage_order = (::roboptim::core::python::NPY_STORAGE_ORDER == NPY_F_CONTIGUOUS)? 'F':'C';
   return Py_BuildValue("c", storage_order);
diff --git a/tests/schittkowski.py b/tests/schittkowski.py
@@ -248,9 +248,8 @@ def test_problem_48(self):
         solver.solve ()
         r = solver.minimum ()
         print (r)
-        if not type(r) == roboptim.core.PySolverError:
-            numpy.testing.assert_almost_equal (r.value, [0.])
-            numpy.testing.assert_almost_equal (r.x, [1., 1., 1., 1., 1.])
+        numpy.testing.assert_almost_equal (r.value, [0.])
+        numpy.testing.assert_almost_equal (r.x, [1., 1., 1., 1., 1.])
 
 if __name__ == '__main__':
     unittest.main ()
diff --git a/tests/wrap.py b/tests/wrap.py
@@ -147,7 +147,7 @@ def jacobian(result, x):
         #self.assertRaises(NotImplementedError,
                           #lambda: roboptim.core.gradient (f, gradient, x, 0))
 
-        jacobian = numpy.zeros((outputSize, inputSize))
+        jacobian = numpy.zeros((outputSize, inputSize), order=roboptim.core.getStorageOrder())
         roboptim.core.jacobian (f, jacobian, x)
         expected_jacobian = numpy.array([[1., 1.], [1., -1.], [10., 15.]])
         numpy.testing.assert_almost_equal (jacobian, expected_jacobian)
