55
66#include " ./cc_factor.h"
77
8+ #include < cstring>
89#include < functional>
910
1011#include < Eigen/Dense>
@@ -35,67 +36,89 @@ namespace {
3536using PyHessianFunc =
3637 std::function<py::tuple(const sym::Valuesd&, const std::vector<index_entry_t >&)>;
3738
39+ /* *
40+ * If Matrix is Eigen::SparseMatrix<double> and matrix is not a scipy.sparse.csc_matrix, or
41+ * if Matrix is any other type and matrix is a scipy.sparse.csc_matrix, throws a value error.
42+ */
43+ template <typename Matrix>
44+ void ThrowIfSparsityMismatch (const py::object& matrix) {
45+ if (std::strcmp (Py_TYPE (matrix.ptr ())->tp_name , " csc_matrix" ) == 0 ) {
46+ throw py::value_error (" Non-sparse matrix expected, scipy.sparse.csc_matrix found instead." );
47+ }
48+ }
49+ template <>
50+ void ThrowIfSparsityMismatch<Eigen::SparseMatrix<double >>(const py::object& matrix) {
51+ if (!py::isinstance (matrix, py::module_::import (" scipy.sparse" ).attr (" csc_matrix" ))) {
52+ throw py::value_error (
53+ fmt::format (" scipy.sparse.csc_matrix expected, found {} instead." , py::type::of (matrix)));
54+ }
55+ }
56+
57+ template <typename Matrix>
3858auto WrapPyHessianFunc (PyHessianFunc&& hessian_func) {
3959 using Vec = Eigen::VectorXd;
40- using Mat = Eigen::MatrixXd;
4160 return [hessian_func = std::move (hessian_func)](
4261 const sym::Valuesd& values, const std::vector<index_entry_t >& keys,
43- Vec* const residual, Mat * const jacobian, Mat * const hessian, Vec* const rhs) {
62+ Vec* const residual, Matrix * const jacobian, Matrix * const hessian, Vec* const rhs) {
4463 const py::tuple out_tuple = hessian_func (values, keys);
4564 if (residual != nullptr ) {
4665 *residual = py::cast<Vec>(out_tuple[0 ]);
4766 }
4867 if (jacobian != nullptr ) {
49- *jacobian = py::cast<Mat>(out_tuple[1 ]);
68+ ThrowIfSparsityMismatch<Matrix>(out_tuple[1 ]);
69+ *jacobian = py::cast<Matrix>(out_tuple[1 ]);
5070 }
5171 if (hessian != nullptr ) {
52- *hessian = py::cast<Mat>(out_tuple[2 ]);
72+ ThrowIfSparsityMismatch<Matrix>(out_tuple[2 ]);
73+ *hessian = py::cast<Matrix>(out_tuple[2 ]);
5374 }
5475 if (rhs != nullptr ) {
5576 *rhs = py::cast<Vec>(out_tuple[3 ]);
5677 }
5778 };
5879}
5980
60- sym::Factord MakeHessianFactorSeparateKeys (PyHessianFunc hessian_func,
61- const std::vector<sym::Key>& keys_to_func ,
62- const std::vector<sym::Key>& keys_to_optimize ) {
63- return sym::Factord ( WrapPyHessianFunc ( std::move (hessian_func)), keys_to_func, keys_to_optimize);
64- }
65-
66- sym::Factord MakeHessianFactorCommonKeys (PyHessianFunc hessian_func,
67- const std::vector<sym::Key>& keys) {
68- return sym::Factord ( WrapPyHessianFunc ( std::move (hessian_func)), keys);
81+ template < typename ... Keys>
82+ sym::Factord MakeHessianFactor (PyHessianFunc hessian_func, const std::vector<Keys>&... keys ,
83+ bool sparse ) {
84+ if (sparse) {
85+ return sym::Factord (WrapPyHessianFunc<Eigen::SparseMatrix< double >>( std::move (hessian_func)),
86+ keys...);
87+ } else {
88+ return sym::Factord (WrapPyHessianFunc<Eigen::MatrixXd>( std::move (hessian_func)), keys...);
89+ }
6990}
7091
7192using PyJacobianFunc =
7293 std::function<py::tuple(const sym::Valuesd&, const std::vector<index_entry_t >&)>;
7394
74- sym::Factord::DenseJacobianFunc WrapPyJacobianFunc (PyJacobianFunc&& jacobian_func) {
75- return sym::Factord::DenseJacobianFunc (
95+ template <typename Matrix>
96+ sym::Factord::JacobianFunc<Matrix> WrapPyJacobianFunc (PyJacobianFunc&& jacobian_func) {
97+ return sym::Factord::JacobianFunc<Matrix>(
7698 [jacobian_func = std::move (jacobian_func)](
7799 const sym::Valuesd& values, const std::vector<index_entry_t >& keys,
78- Eigen::VectorXd* const residual, Eigen::MatrixXd * const jacobian) {
100+ Eigen::VectorXd* const residual, Matrix * const jacobian) {
79101 const py::tuple out_tuple = jacobian_func (values, keys);
80102 if (residual != nullptr ) {
81103 *residual = py::cast<Eigen::VectorXd>(out_tuple[0 ]);
82104 }
83105 if (jacobian != nullptr ) {
84- *jacobian = py::cast<Eigen::MatrixXd>(out_tuple[1 ]);
106+ ThrowIfSparsityMismatch<Matrix>(out_tuple[1 ]);
107+ *jacobian = py::cast<Matrix>(out_tuple[1 ]);
85108 }
86109 });
87110}
88111
89- sym::Factord MakeJacobianFactorSeparateKeys (PyJacobianFunc jacobian_func,
90- const std::vector<sym::Key>& keys_to_func ,
91- const std::vector<sym::Key>& keys_to_optimize ) {
92- return sym::Factord::Jacobian ( WrapPyJacobianFunc ( std::move (jacobian_func)), keys_to_func,
93- keys_to_optimize);
94- }
95-
96- sym::Factord MakeJacobianFactorCommonKeys (PyJacobianFunc jacobian_func,
97- const std::vector<sym::Key>& keys) {
98- return sym::Factord::Jacobian ( WrapPyJacobianFunc ( std::move (jacobian_func)), keys);
112+ template < typename ... Keys>
113+ sym::Factord MakeJacobianFactor (PyJacobianFunc jacobian_func, const std::vector<Keys>&... keys ,
114+ bool sparse ) {
115+ if (sparse) {
116+ return sym::Factord::Jacobian (
117+ WrapPyJacobianFunc<Eigen::SparseMatrix< double >>( std::move (jacobian_func)), keys...);
118+ } else {
119+ return sym::Factord::Jacobian (WrapPyJacobianFunc<Eigen::MatrixXd>( std::move ( jacobian_func)) ,
120+ keys...);
121+ }
99122}
100123
101124} // namespace
@@ -112,35 +135,49 @@ void AddFactorWrapper(pybind11::module_ module) {
112135 point, generates a linear approximation to the residual function.
113136 )" )
114137 // TODO(brad): Add wrapper of the constructor from SparseHessianFunc
115- .def (py::init (&MakeHessianFactorCommonKeys), py::arg (" hessian_func" ), py::arg (" keys" ), R"(
116- Create directly from a (dense) hessian functor. This is the lowest-level constructor.
138+ .def (py::init (&MakeHessianFactor<sym::Key>), py::arg (" hessian_func" ), py::arg (" keys" ),
139+ py::arg (" sparse" ) = false , R"(
140+ Create directly from a hessian functor. This is the lowest-level constructor.
117141
118142 Args:
119143 keys: The set of input arguments, in order, accepted by func.
144+ sparse: Create a sparse factor if True, dense factor if false. Defaults to dense.
145+
146+ Precondition:
147+ The jacobian and hessian returned by hessian_func have type scipy.sparse.csc_matrix if and only if sparse = True.
120148 )" )
121- .def (py::init (&MakeHessianFactorSeparateKeys ), py::arg (" hessian_func" ),
122- py::arg (" keys_to_func" ), py::arg (" keys_to_optimize" ),
149+ .def (py::init (&MakeHessianFactor<sym::Key, sym::Key> ), py::arg (" hessian_func" ),
150+ py::arg (" keys_to_func" ), py::arg (" keys_to_optimize" ), py::arg ( " sparse " ) = false ,
123151 R"(
124- Create directly from a (sparse) hessian functor. This is the lowest-level constructor.
152+ Create directly from a hessian functor. This is the lowest-level constructor.
125153
126154 Args:
127155 keys_to_func: The set of input arguments, in order, accepted by func.
128156 keys_to_optimize: The set of input arguments that correspond to the derivative in func. Must be a subset of keys_to_func.
157+ sparse: Create a sparse factor if True, dense factor if false. Defaults to dense.
158+
159+ Precondition:
160+ The jacobian and hessian returned by hessian_func have type scipy.sparse.csc_matrix if and only if sparse = True.
129161 )" )
130162 .def (" is_sparse" , &sym::Factord::IsSparse,
131163 " Does this factor use a sparse jacobian/hessian matrix?" )
132- .def_static (" jacobian" , &MakeJacobianFactorCommonKeys , py::arg (" jacobian_func" ),
133- py::arg (" keys" ), R"(
164+ .def_static (" jacobian" , &MakeJacobianFactor<sym::Key> , py::arg (" jacobian_func" ),
165+ py::arg (" keys" ), py::arg ( " sparse " ) = false , R"(
134166 Create from a function that computes the jacobian. The hessian will be computed using the
135167 Gauss Newton approximation:
136168 H = J.T * J
137169 rhs = J.T * b
138170
139171 Args:
140172 keys: The set of input arguments, in order, accepted by func.
173+ sparse: Create a sparse factor if True, dense factor if false. Defaults to dense.
174+
175+ Precondition:
176+ The jacobian returned by jacobian_func has type scipy.sparse.csc_matrix if and only if sparse = True.
141177 )" )
142- .def_static (" jacobian" , &MakeJacobianFactorSeparateKeys, py::arg (" jacobian_func" ),
143- py::arg (" keys_to_func" ), py::arg (" keys_to_optimize" ), R"(
178+ .def_static (" jacobian" , &MakeJacobianFactor<sym::Key, sym::Key>, py::arg (" jacobian_func" ),
179+ py::arg (" keys_to_func" ), py::arg (" keys_to_optimize" ), py::arg (" sparse" ) = false ,
180+ R"(
144181 Create from a function that computes the jacobian. The hessian will be computed using the
145182 Gauss Newton approximation:
146183 H = J.T * J
@@ -149,6 +186,10 @@ void AddFactorWrapper(pybind11::module_ module) {
149186 Args:
150187 keys_to_func: The set of input arguments, in order, accepted by func.
151188 keys_to_optimize: The set of input arguments that correspond to the derivative in func. Must be a subset of keys_to_func.
189+ sparse: Create a sparse factor if True, dense factor if false. Defaults to dense.
190+
191+ Precondition:
192+ The jacobian returned by jacobian_func has type scipy.sparse.csc_matrix if and only if sparse = True.
152193 )" )
153194 .def (
154195 " linearize" ,
0 commit comments