diff --git a/CMakeLists.txt b/CMakeLists.txt index 8acb0ef3541fb22c64eef8148a37497e690f1d83..ac1e319340da67ad5b5c95f5b50a23a1e8679c18 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ # -# Copyright (c) 2015-2016 LAAS-CNRS +# Copyright (c) 2015-2017 LAAS-CNRS # # This file is part of eigenpy. # eigenpy is free software: you can redistribute it and/or @@ -19,6 +19,7 @@ CMAKE_MINIMUM_REQUIRED(VERSION 2.6) INCLUDE(cmake/base.cmake) INCLUDE(cmake/boost.cmake) INCLUDE(cmake/python.cmake) +INCLUDE(cmake/ide.cmake) SET(PROJECT_NAME eigenpy) SET(PROJECT_DESCRIPTION "Wrapping Eigen3 -- numpy") @@ -45,6 +46,12 @@ IF(APPLE) ENDIF("${isSystemDir}" STREQUAL "-1") ENDIF(APPLE) +IF(WIN32) + SET(LINK copy_if_different) +ELSE(WIN32) + SET(LINK create_symlink) +ENDIF(WIN32) + # ---------------------------------------------------- # --- OPTIONS --------------------------------------- # ---------------------------------------------------- @@ -73,61 +80,97 @@ INCLUDE_DIRECTORIES(${NUMPY_INCLUDE_DIRS}) # ---------------------------------------------------- # --- INCLUDE ---------------------------------------- # ---------------------------------------------------- +SET(${PROJECT_NAME}_SOLVERS_HEADERS + solvers/solvers.hpp + solvers/preconditioners.hpp + solvers/IterativeSolverBase.hpp + solvers/LeastSquaresConjugateGradient.hpp + solvers/ConjugateGradient.hpp + solvers/SparseSolverBase.hpp + solvers/BasicPreconditioners.hpp + solvers/BFGSPreconditioners.hpp + ) + SET(HEADERS - src/eigenpy.hpp - src/exception.hpp - src/details.hpp - src/fwd.hpp - src/map.hpp - src/geometry.hpp - src/memory.hpp - src/registration.hpp - src/angle-axis.hpp - src/quaternion.hpp + ${${PROJECT_NAME}_SOLVERS_HEADERS} + eigenpy.hpp + exception.hpp + details.hpp + fwd.hpp + map.hpp + geometry.hpp + memory.hpp + registration.hpp + angle-axis.hpp + quaternion.hpp ) + MAKE_DIRECTORY("${${PROJECT_NAME}_BINARY_DIR}/include/eigenpy") +MAKE_DIRECTORY("${${PROJECT_NAME}_BINARY_DIR}/include/eigenpy/solvers") INCLUDE_DIRECTORIES(${${PROJECT_NAME}_BINARY_DIR}/include/eigenpy) +#FOREACH(header ${HEADERS}) +# GET_FILENAME_COMPONENT(headerName ${header} NAME) +# IF(WIN32) +# execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different +# ${${PROJECT_NAME}_SOURCE_DIR}/${header} +# ${${PROJECT_NAME}_BINARY_DIR}/include/${PROJECT_NAME}/) +# ELSE(WIN32) +# execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink +# ${${PROJECT_NAME}_SOURCE_DIR}/${header} +# ${${PROJECT_NAME}_BINARY_DIR}/include/${PROJECT_NAME}/${headerName}) +# ENDIF(WIN32) +# INSTALL(FILES ${${PROJECT_NAME}_SOURCE_DIR}/${header} +# DESTINATION ${CMAKE_INSTALL_PREFIX}/include/${PROJECT_NAME} +# PERMISSIONS OWNER_READ GROUP_READ WORLD_READ) +#ENDFOREACH(header) + +SET(${PROJECT_NAME}_HEADERS) FOREACH(header ${HEADERS}) GET_FILENAME_COMPONENT(headerName ${header} NAME) - IF(WIN32) - execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different - ${${PROJECT_NAME}_SOURCE_DIR}/${header} - ${${PROJECT_NAME}_BINARY_DIR}/include/${PROJECT_NAME}/) - ELSE(WIN32) - execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink - ${${PROJECT_NAME}_SOURCE_DIR}/${header} - ${${PROJECT_NAME}_BINARY_DIR}/include/${PROJECT_NAME}/${headerName}) - ENDIF(WIN32) - INSTALL(FILES ${${PROJECT_NAME}_SOURCE_DIR}/${header} - DESTINATION ${CMAKE_INSTALL_PREFIX}/include/${PROJECT_NAME} - PERMISSIONS OWNER_READ GROUP_READ WORLD_READ) + GET_FILENAME_COMPONENT(headerPath ${header} PATH) + EXECUTE_PROCESS(COMMAND ${CMAKE_COMMAND} -E ${LINK} + ${${PROJECT_NAME}_SOURCE_DIR}/src/${header} + ${${PROJECT_NAME}_BINARY_DIR}/include/${PROJECT_NAME}/${header}) + INSTALL(FILES ${${PROJECT_NAME}_SOURCE_DIR}/src/${header} + DESTINATION ${CMAKE_INSTALL_PREFIX}/include/${PROJECT_NAME}/${headerPath} + PERMISSIONS OWNER_READ GROUP_READ WORLD_READ OWNER_WRITE) + LIST(APPEND ${PROJECT_NAME}_HEADERS src/${header}) ENDFOREACH(header) # ---------------------------------------------------- # --- TARGETS ---------------------------------------- # ---------------------------------------------------- +SET(${PROJECT_NAME}_SOLVERS_SOURCES + src/solvers/preconditioners.cpp + src/solvers/solvers.cpp + ) SET(${PROJECT_NAME}_SOURCES + ${${PROJECT_NAME}_SOLVERS_SOURCES} src/exception.cpp src/eigenpy.cpp src/angle-axis.cpp src/quaternion.cpp ) -ADD_LIBRARY(${PROJECT_NAME} SHARED ${${PROJECT_NAME}_SOURCES} ${HEADERS}) +ADD_LIBRARY(${PROJECT_NAME} SHARED ${${PROJECT_NAME}_SOURCES} ${${PROJECT_NAME}_HEADERS}) TARGET_LINK_BOOST_PYTHON(${PROJECT_NAME}) PKG_CONFIG_USE_DEPENDENCY(${PROJECT_NAME} eigen3) INSTALL(TARGETS ${PROJECT_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/lib) +ADD_HEADER_GROUP(${PROJECT_NAME}_HEADERS) +ADD_SOURCE_GROUP(${PROJECT_NAME}_SOURCES) + # ---------------------------------------------------- -# --- UNIT TEST -------------------------------------- +# --- PYTHON LIBRARY --------------------------------- # ---------------------------------------------------- -ADD_SUBDIRECTORY(unittest) +ADD_SUBDIRECTORY(python) # ---------------------------------------------------- -# --- EXECUTABLES ------------------------------------ +# --- UNIT TEST -------------------------------------- # ---------------------------------------------------- +ADD_SUBDIRECTORY(unittest) IF(EIGEN_NUMPY_ALIGNED) PKG_CONFIG_APPEND_CFLAGS("-DEIGENPY_ALIGNED") diff --git a/cmake b/cmake index 45ea8074e31604d3e4cc594baaa5abe9fac1921b..687a16b8d42e20294ec3f51759d122db2d212695 160000 --- a/cmake +++ b/cmake @@ -1 +1 @@ -Subproject commit 45ea8074e31604d3e4cc594baaa5abe9fac1921b +Subproject commit 687a16b8d42e20294ec3f51759d122db2d212695 diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d2c3216c4f29aa98c98d8057e8069de856d9c1fb --- /dev/null +++ b/python/CMakeLists.txt @@ -0,0 +1,88 @@ +# +# Copyright (c) 2017 CNRS +# +# This file is part of eigenpy +# eigenpy is free software: you can redistribute it +# and/or modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation, either version +# 3 of the License, or (at your option) any later version. +# eigenpy is distributed in the hope that it will be +# useful, but WITHOUT ANY WARRANTY; without even the implied warranty +# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Lesser Public License for more details. You should have +# received a copy of the GNU Lesser General Public License along with +# eigenpy If not, see +# <http://www.gnu.org/licenses/>. + +MACRO(SYMLINK_AND_INSTALL_HEADERS HEADERS SUBPATH) + FOREACH(header ${HEADERS}) + GET_FILENAME_COMPONENT(headerName ${header} NAME) + GET_FILENAME_COMPONENT(headerPath ${header} PATH) + EXECUTE_PROCESS(COMMAND ${CMAKE_COMMAND} -E ${LINK} + ${CMAKE_CURRENT_SOURCE_DIR}/${header} + ${${PROJECT_NAME}_BINARY_DIR}/include/${PROJECT_NAME}/${SUBPATH}/${header}) + + INSTALL(FILES ${CMAKE_CURRENT_SOURCE_DIR}/${header} + DESTINATION ${CMAKE_INSTALL_PREFIX}/include/${PROJECT_NAME}/${SUBPATH}/${headerPath} + PERMISSIONS OWNER_READ GROUP_READ WORLD_READ OWNER_WRITE) + ENDFOREACH(header) +ENDMACRO(SYMLINK_AND_INSTALL_HEADERS HEADERS SUBPATH) + +# --- LIBRARY --- # +SET(PYWRAP ${PROJECT_NAME}_pywrap) +MAKE_DIRECTORY("${${PROJECT_NAME}_BINARY_DIR}/python/${PROJECT_NAME}") + +ADD_LIBRARY(${PYWRAP} SHARED main.cpp) +TARGET_LINK_LIBRARIES(${PYWRAP} ${PROJECT_NAME}) +TARGET_LINK_BOOST_PYTHON(${PYWRAP}) +#IF(BUILD_WITH_COMMIT_VERSION) +# TAG_LIBRARY_VERSION(${PYWRAP}) +#ENDIF(BUILD_WITH_COMMIT_VERSION) +SET(${PYWRAP}_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/${PYTHON_SITELIB}) + +SET_PROPERTY(TARGET ${PYWRAP} PROPERTY LINKER_LANGUAGE CXX) +IF(APPLE) + # We need to change the extension for python bindings + SET_TARGET_PROPERTIES(${PYWRAP} PROPERTIES SUFFIX ".so") +ENDIF(APPLE) + +SET_TARGET_PROPERTIES(${PYWRAP} PROPERTIES + LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/python/${PROJECT_NAME}") + +SET_TARGET_PROPERTIES(${PYWRAP} PROPERTIES PREFIX "") +SET_TARGET_PROPERTIES(${PYWRAP} PROPERTIES OUTPUT_NAME "${PROJECT_NAME}") + +INSTALL(TARGETS ${PYWRAP} DESTINATION ${${PYWRAP}_INSTALL_DIR}) + +## --- INSTALL SCRIPTS +#SET(PYTHON_FILES +# __init__.py +# ) +# +#FOREACH(python ${PYTHON_FILES}) +# GET_FILENAME_COMPONENT(pythonFile ${python} NAME) +# EXECUTE_PROCESS(COMMAND ${CMAKE_COMMAND} -E ${LINK} +# ${${PROJECT_NAME}_SOURCE_DIR}/python/scripts/${python} +# ${${PROJECT_NAME}_BINARY_DIR}/python/${PROJECT_NAME}/${pythonFile}) +# +# # Generate pyc file +# EXECUTE_PROCESS(COMMAND +# ${PYTHON_EXECUTABLE} -m py_compile +# ${${PROJECT_NAME}_BINARY_DIR}/python/${PROJECT_NAME}/${pythonFile}) +# # Tag pyc file as generated. +# SET_SOURCE_FILES_PROPERTIES( +# "${${PROJECT_NAME}_BINARY_DIR}/python/${PROJECT_NAME}/${pythonFile}c" +# PROPERTIES GENERATED TRUE) +# +# # Clean generated files. +# SET_PROPERTY( +# DIRECTORY APPEND PROPERTY +# ADDITIONAL_MAKE_CLEAN_FILES +# "${${PROJECT_NAME}_BINARY_DIR}/python/${PROJECT_NAME}/${pythonFile}c") +# +# INSTALL(FILES +# "${${PROJECT_NAME}_SOURCE_DIR}/python/scripts/${python}" +# "${${PROJECT_NAME}_BINARY_DIR}/python/${PROJECT_NAME}/${pythonFile}c" +# DESTINATION ${${PYWRAP}_INSTALL_DIR}) +#ENDFOREACH(python) + diff --git a/python/main.cpp b/python/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b2ce82232b2bc7e9f8e30ffce094b854c660cfc5 --- /dev/null +++ b/python/main.cpp @@ -0,0 +1,40 @@ +// +// Copyright (c) 2017, Justin Carpentier, CNRS +// +// This file is part of eigenpy +// eigenpy is free software: you can redistribute it +// and/or modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation, either version +// 3 of the License, or (at your option) any later version. +// +// eigenpy is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty +// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Lesser Public License for more details. You should have +// received a copy of the GNU Lesser General Public License along with +// eigenpy If not, see +// <http://www.gnu.org/licenses/>. + +#include "eigenpy/eigenpy.hpp" +#include "eigenpy/geometry.hpp" +#include "eigenpy/solvers/solvers.hpp" +#include "eigenpy/solvers/preconditioners.hpp" + +#include <iostream> +#include <boost/python/scope.hpp> + +using namespace eigenpy; + +BOOST_PYTHON_MODULE(eigenpy) +{ + enableEigenPy(); + exposeAngleAxis(); + exposeQuaternion(); + + { + boost::python::scope solvers = boost::python::class_<SolversScope>("solvers"); + exposeSolvers(); + exposePreconditioners(); + } + +} diff --git a/python/scripts/__init__.py b/python/scripts/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..3f79365cdb1298271dd8fe3b7a1dff11a4c4352f --- /dev/null +++ b/python/scripts/__init__.py @@ -0,0 +1,18 @@ +# +# Copyright (c) 2017 CNRS +# +# This file is part of eigenpy +# eigenpy is free software: you can redistribute it +# and/or modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation, either version +# 3 of the License, or (at your option) any later version. +# Pinocchio is distributed in the hope that it will be +# useful, but WITHOUT ANY WARRANTY; without even the implied warranty +# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Lesser Public License for more details. You should have +# received a copy of the GNU Lesser General Public License along with +# eigenpy If not, see +# <http://www.gnu.org/licenses/>. + +import numpy as np +from eigenpy import * diff --git a/src/quaternion.hpp b/src/quaternion.hpp index 811bb50dd93dd4ff59fa651786aaac7cedd17f31..b35a767246259dee6a495c16fccbca9ac2d82cd1 100644 --- a/src/quaternion.hpp +++ b/src/quaternion.hpp @@ -71,16 +71,16 @@ namespace eigenpy "The [] operator numbers them differently, 0...4 for *x* *y* *z* *w*!")) .add_property("x", - (Scalar (Quaternion::*)()const)&Quaternion::x, + &QuaternionVisitor::getCoeff<0>, &QuaternionVisitor::setCoeff<0>,"The x coefficient.") .add_property("y", - (Scalar (Quaternion::*)()const)&Quaternion::y, + &QuaternionVisitor::getCoeff<1>, &QuaternionVisitor::setCoeff<1>,"The y coefficient.") .add_property("z", - (Scalar (Quaternion::*)()const)&Quaternion::z, + &QuaternionVisitor::getCoeff<2>, &QuaternionVisitor::setCoeff<2>,"The z coefficient.") .add_property("w", - (Scalar (Quaternion::*)()const)&Quaternion::w, + &QuaternionVisitor::getCoeff<3>, &QuaternionVisitor::setCoeff<3>,"The w coefficient.") // .def("isApprox",(bool (Quaternion::*)(const Quaternion &))&Quaternion::template isApprox<Quaternion>, @@ -151,6 +151,9 @@ namespace eigenpy template<int i> static void setCoeff(Quaternion & self, Scalar value) { self.coeffs()[i] = value; } + template<int i> + static Scalar getCoeff(Quaternion & self) { return self.coeffs()[i]; } + static Quaternion & setFromTwoVectors(Quaternion & self, const Vector3 & a, const Vector3 & b) { return self.setFromTwoVectors(a,b); } diff --git a/src/solvers/BFGSPreconditioners.hpp b/src/solvers/BFGSPreconditioners.hpp new file mode 100644 index 0000000000000000000000000000000000000000..b24c30e998eaaac83b4f85aa2461b5d85898963b --- /dev/null +++ b/src/solvers/BFGSPreconditioners.hpp @@ -0,0 +1,90 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_bfgs_preconditioners_hpp__ +#define __eigenpy_bfgs_preconditioners_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> +#include <Eigen/IterativeLinearSolvers> + +#include "eigenpy/solvers/BasicPreconditioners.hpp" + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename Preconditioner> + struct BFGSPreconditionerBaseVisitor + : public bp::def_visitor< BFGSPreconditionerBaseVisitor<Preconditioner> > + { + + typedef Eigen::VectorXd VectorType; + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(PreconditionerBaseVisitor<Preconditioner>()) + .def("rows",&Preconditioner::rows,"Returns the number of rows in the preconditioner.") + .def("cols",&Preconditioner::cols,"Returns the number of cols in the preconditioner.") + .def("dim",&Preconditioner::dim,"Returns the dimension of the BFGS preconditioner") + .def("update",(const Preconditioner& (Preconditioner::*)(const VectorType &, const VectorType &) const)&Preconditioner::update,bp::args("s","y"),"Update the BFGS estimate of the matrix A.",bp::return_value_policy<bp::reference_existing_object>()) + .def("reset",&Preconditioner::reset,"Reset the BFGS estimate.") + ; + + } + + static void expose(const std::string & name) + { + bp::class_<Preconditioner>(name, + bp::no_init) + .def(BFGSPreconditionerBaseVisitor<Preconditioner>()) + ; + + } + + }; + + template<typename Preconditioner> + struct LimitedBFGSPreconditionerBaseVisitor + : public bp::def_visitor< LimitedBFGSPreconditionerBaseVisitor<Preconditioner> > + { + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(PreconditionerBaseVisitor<Preconditioner>()) + .def(BFGSPreconditionerBaseVisitor<Preconditioner>()) + .def("resize",&Preconditioner::resize,bp::arg("dim"),"Resizes the preconditionner with size dim.",bp::return_value_policy<bp::reference_existing_object>()) + ; + + } + + static void expose(const std::string & name) + { + bp::class_<Preconditioner>(name.c_str(), + bp::no_init) + .def(LimitedBFGSPreconditionerBaseVisitor<Preconditioner>()) + ; + + } + + }; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_bfgs_preconditioners_hpp__ diff --git a/src/solvers/BasicPreconditioners.hpp b/src/solvers/BasicPreconditioners.hpp new file mode 100644 index 0000000000000000000000000000000000000000..92a9af5cbbe95272e062f088fb86c46955e58266 --- /dev/null +++ b/src/solvers/BasicPreconditioners.hpp @@ -0,0 +1,143 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_basic_preconditioners_hpp__ +#define __eigenpy_basic_preconditioners_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> +#include <Eigen/IterativeLinearSolvers> + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename Preconditioner> + struct PreconditionerBaseVisitor + : public bp::def_visitor< PreconditionerBaseVisitor<Preconditioner> > + { + typedef Eigen::MatrixXd MatrixType; + typedef Eigen::VectorXd VectorType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(bp::init<>("Default constructor")) + .def(bp::init<MatrixType>(bp::arg("A"),"Initialize the preconditioner with matrix A for further Az=b solving.")) +#if EIGEN_VERSION_AT_LEAST(3,3,0) + .def("info",&Preconditioner::info, + "Returns success if the Preconditioner has been well initialized.") +#endif + .def("solve",&solve,bp::arg("b"), + "Returns the solution A * z = b where the preconditioner is an estimate of A^-1.") + + .def("compute",&Preconditioner::template compute<MatrixType>,bp::arg("mat"), + "Initialize the preconditioner from the matrix value.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("factorize",&Preconditioner::template factorize<MatrixType>,bp::arg("mat"), + "Initialize the preconditioner from the matrix value, i.e factorize the mat given as input to approximate its inverse.", + bp::return_value_policy<bp::reference_existing_object>()) + ; + + } + + private: + + static VectorType solve(Preconditioner & self, const VectorType & b) + { + return self.solve(b); + } + + }; + + template<typename Scalar> + struct DiagonalPreconditionerVisitor : PreconditionerBaseVisitor<Eigen::DiagonalPreconditioner<Scalar> > + { + typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixType; + typedef Eigen::DiagonalPreconditioner<Scalar> Preconditioner; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(PreconditionerBaseVisitor<Preconditioner>()) + .def("rows",&Preconditioner::rows,"Returns the number of rows in the preconditioner.") + .def("cols",&Preconditioner::cols,"Returns the number of cols in the preconditioner.") + ; + + } + + static void expose() + { + bp::class_<Preconditioner>("DiagonalPreconditioner", + "A preconditioner based on the digonal entrie.\n" + "This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.", + bp::no_init) + ; + + } + }; + +#if EIGEN_VERSION_AT_LEAST(3,3,0) + template<typename Scalar> + struct LeastSquareDiagonalPreconditionerVisitor : PreconditionerBaseVisitor<Eigen::LeastSquareDiagonalPreconditioner<Scalar> > + { + typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixType; + typedef Eigen::LeastSquareDiagonalPreconditioner<Scalar> Preconditioner; + + template<class PyClass> + void visit(PyClass&) const + { + } + + static void expose() + { + bp::class_<Preconditioner>("LeastSquareDiagonalPreconditioner", + "Jacobi preconditioner for LeastSquaresConjugateGradient.\n" + "his class allows to approximately solve for A' A x = A' b problems assuming A' A is a diagonal matrix.", + bp::no_init) + .def(DiagonalPreconditionerVisitor<Scalar>()) + ; + + } + }; +#endif + + struct IdentityPreconditionerVisitor : PreconditionerBaseVisitor<Eigen::IdentityPreconditioner > + { + typedef Eigen::IdentityPreconditioner Preconditioner; + + template<class PyClass> + void visit(PyClass&) const + { + } + + static void expose() + { + bp::class_<Preconditioner>("IdentityPreconditioner", + bp::no_init) + .def(PreconditionerBaseVisitor<Preconditioner>()) + ; + + } + }; + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_basic_preconditioners_hpp__ diff --git a/src/solvers/ConjugateGradient.hpp b/src/solvers/ConjugateGradient.hpp new file mode 100644 index 0000000000000000000000000000000000000000..f3d1487e01f103d0dd3848e5f63e318548e0ca40 --- /dev/null +++ b/src/solvers/ConjugateGradient.hpp @@ -0,0 +1,61 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_conjugate_gradient_hpp__ +#define __eigenpy_conjugate_gradient_hpp__ + +#include <boost/python.hpp> +#include <Eigen/IterativeLinearSolvers> + +#include "eigenpy/solvers/IterativeSolverBase.hpp" + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename ConjugateGradient> + struct ConjugateGradientVisitor + : public boost::python::def_visitor< ConjugateGradientVisitor<ConjugateGradient> > + { + typedef typename ConjugateGradient::MatrixType MatrixType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(IterativeSolverVisitor<ConjugateGradient>()) + .def(bp::init<>("Default constructor")) + .def(bp::init<MatrixType>(bp::arg("A"),"Initialize the solver with matrix A for further Ax=b solving.\n" + "This constructor is a shortcut for the default constructor followed by a call to compute().")) + ; + + } + + static void expose(const std::string & name = "ConjugateGradient") + { + bp::class_<ConjugateGradient,boost::noncopyable>(name.c_str(), + bp::no_init) + .def(ConjugateGradientVisitor<ConjugateGradient>()) + ; + + } + + }; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_conjugate_gradient_hpp__ diff --git a/src/solvers/IterativeSolverBase.hpp b/src/solvers/IterativeSolverBase.hpp new file mode 100644 index 0000000000000000000000000000000000000000..0c617587aa2dccf31d15f0f0476641a23332cc56 --- /dev/null +++ b/src/solvers/IterativeSolverBase.hpp @@ -0,0 +1,104 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_iterative_solver_base_hpp__ +#define __eigenpy_iterative_solver_base_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> + +#include "eigenpy/solvers/SparseSolverBase.hpp" + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename IterativeSolver> + struct IterativeSolverVisitor + : public boost::python::def_visitor< IterativeSolverVisitor<IterativeSolver> > + { + + typedef typename IterativeSolver::MatrixType MatrixType; + typedef typename IterativeSolver::Preconditioner Preconditioner; + typedef Eigen::VectorXd VectorType; + + template<class PyClass> + void visit(PyClass& cl) const + { + typedef IterativeSolver IS; + + cl + .def(SparseSolverVisitor<IS>()) + .def("error",&IS::error,"Returns the tolerance error reached during the last solve.\n" + "It is a close approximation of the true relative residual error |Ax-b|/|b|.") + .def("info",&IS::info,"Returns success if the iterations converged, and NoConvergence otherwise.") + .def("iterations",&IS::iterations,"Returns the number of iterations performed during the last solve.") + .def("maxIterations",&IS::maxIterations,"Returns the max number of iterations.\n" + "It is either the value setted by setMaxIterations or, by default, twice the number of columns of the matrix.") + .def("setMaxIterations",&IS::setMaxIterations,"Sets the max number of iterations.\n" + "Default is twice the number of columns of the matrix.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("tolerance",&IS::tolerance,"Returns he tolerance threshold used by the stopping criteria.") + .def("setTolerance",&IS::setTolerance,"Sets the tolerance threshold used by the stopping criteria.\n" + "This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. The default value is the machine precision.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("analyzePattern",&analyzePattern,bp::arg("A"),"Initializes the iterative solver for the sparsity pattern of the matrix A for further solving Ax=b problems.\n" + "Currently, this function mostly calls analyzePattern on the preconditioner.\n" + "In the future we might, for instance, implement column reordering for faster matrix vector products.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("factorize",&factorize,bp::arg("A"),"Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.\n" + "Currently, this function mostly calls factorize on the preconditioner.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("compute",&compute,bp::arg("A"),"Initializes the iterative solver with the numerical values of the matrix A for further solving Ax=b problems.\n" + "Currently, this function mostly calls factorize on the preconditioner.\n" + "In the future we might, for instance, implement column reordering for faster matrix vector products.", + bp::return_value_policy<bp::reference_existing_object>()) + .def("solveWithGuess",&solveWithGuess,bp::args("b","x0"), + "Returns the solution x of Ax = b using the current decomposition of A and x0 as an initial solution.") + .def("preconditioner",(Preconditioner & (IS::*)(void))&IS::preconditioner,"Returns a read-write reference to the preconditioner for custom configuration.",bp::return_internal_reference<>()) + ; + + } + + private: + + static IterativeSolver & factorize(IterativeSolver & self, const MatrixType & m) + { + return self.factorize(m); + } + + static IterativeSolver & compute(IterativeSolver & self, const MatrixType & m) + { + return self.compute(m); + } + + static IterativeSolver & analyzePattern(IterativeSolver & self, const MatrixType & m) + { + return self.analyzePattern(m); + } + + static VectorType solveWithGuess(IterativeSolver & self, const Eigen::VectorXd & b, const Eigen::VectorXd & x0) + { + return self.solveWithGuess(b,x0); + } + + }; + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_iterative_solver_base_hpp__ diff --git a/src/solvers/LeastSquaresConjugateGradient.hpp b/src/solvers/LeastSquaresConjugateGradient.hpp new file mode 100644 index 0000000000000000000000000000000000000000..43957f1bed8bb60e49f3f8219764ff09618a69de --- /dev/null +++ b/src/solvers/LeastSquaresConjugateGradient.hpp @@ -0,0 +1,122 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_least_square_conjugate_gradient_hpp__ +#define __eigenpy_least_square_conjugate_gradient_hpp__ + +#include <boost/python.hpp> +#include <Eigen/IterativeLinearSolvers> + +#include "eigenpy/solvers/IterativeSolverBase.hpp" + +namespace Eigen +{ + template <typename _Scalar> + class LeastSquareDiagonalPreconditionerFix + : public LeastSquareDiagonalPreconditioner<_Scalar> + { + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef LeastSquareDiagonalPreconditioner<_Scalar> Base; + using DiagonalPreconditioner<_Scalar>::m_invdiag; + public: + + LeastSquareDiagonalPreconditionerFix() : Base() {} + + template<typename MatType> + explicit LeastSquareDiagonalPreconditionerFix(const MatType& mat) : Base() + { + compute(mat); + } + + template<typename MatType> + LeastSquareDiagonalPreconditionerFix& analyzePattern(const MatType& ) + { + return *this; + } + + template<typename MatType> + LeastSquareDiagonalPreconditionerFix& factorize(const MatType& mat) + { + // Compute the inverse squared-norm of each column of mat + m_invdiag.resize(mat.cols()); + if(MatType::IsRowMajor) + { + m_invdiag.setZero(); + for(Index j=0; j<mat.outerSize(); ++j) + { + for(typename MatType::InnerIterator it(mat,j); it; ++it) + m_invdiag(it.index()) += numext::abs2(it.value()); + } + for(Index j=0; j<mat.cols(); ++j) + if(numext::real(m_invdiag(j))>RealScalar(0)) + m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j)); + } + else + { + for(Index j=0; j<mat.outerSize(); ++j) + { + RealScalar sum = mat.col(j).squaredNorm(); + if(sum>RealScalar(0)) + m_invdiag(j) = RealScalar(1)/sum; + else + m_invdiag(j) = RealScalar(1); + } + } + Base::m_isInitialized = true; + return *this; + } + + }; +} + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename LeastSquaresConjugateGradient> + struct LeastSquaresConjugateGradientVisitor + : public boost::python::def_visitor< LeastSquaresConjugateGradientVisitor<LeastSquaresConjugateGradient> > + { + typedef Eigen::MatrixXd MatrixType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def(bp::init<>("Default constructor")) + .def(bp::init<MatrixType>(bp::arg("A"),"Initialize the solver with matrix A for further || Ax - b || solving.\n" + "This constructor is a shortcut for the default constructor followed by a call to compute().")) + ; + + } + + static void expose() + { + bp::class_<LeastSquaresConjugateGradient,boost::noncopyable>("LeastSquaresConjugateGradient", + bp::no_init) + .def(IterativeSolverVisitor<LeastSquaresConjugateGradient>()) + .def(LeastSquaresConjugateGradientVisitor<LeastSquaresConjugateGradient>()) + ; + + } + + }; + +} // namespace eigenpy + +#endif // ifndef __eigenpy_least_square_conjugate_gradient_hpp__ diff --git a/src/solvers/SparseSolverBase.hpp b/src/solvers/SparseSolverBase.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e96a9741934b616871bc91e8d9399c94c4080d12 --- /dev/null +++ b/src/solvers/SparseSolverBase.hpp @@ -0,0 +1,58 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_sparse_solver_base_hpp__ +#define __eigenpy_sparse_solver_base_hpp__ + +#include <boost/python.hpp> +#include <Eigen/Core> + +namespace eigenpy +{ + + namespace bp = boost::python; + + template<typename SparseSolver> + struct SparseSolverVisitor + : public bp::def_visitor< SparseSolverVisitor<SparseSolver> > + { + + typedef Eigen::VectorXd VectorType; + + template<class PyClass> + void visit(PyClass& cl) const + { + cl + .def("solve",&solve,bp::arg("b"), + "Returns the solution x of Ax = b using the current decomposition of A.") + ; + + } + + private: + + static VectorType solve(SparseSolver & self, const VectorType & b) + { + return self.solve(b); + } + + }; + + + +} // namespace eigenpy + +#endif // ifndef __eigenpy_sparse_solver_base_hpp__ diff --git a/src/solvers/preconditioners.cpp b/src/solvers/preconditioners.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bea99da1029f4eb6d68635fd99511ac7d7366c8c --- /dev/null +++ b/src/solvers/preconditioners.cpp @@ -0,0 +1,35 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#include "eigenpy/solvers/preconditioners.hpp" +#include "eigenpy/solvers/BasicPreconditioners.hpp" +//#include "eigenpy/solvers/BFGSPreconditioners.hpp" + +namespace eigenpy +{ + void exposePreconditioners() + { + using namespace Eigen; + + DiagonalPreconditionerVisitor<double>::expose(); +#if EIGEN_VERSION_AT_LEAST(3,3,0) + LeastSquareDiagonalPreconditionerVisitor<double>::expose(); +#endif + IdentityPreconditionerVisitor::expose(); +// LimitedBFGSPreconditionerBaseVisitor< LimitedBFGSPreconditioner<double,Eigen::Dynamic,Eigen::Upper|Eigen::Lower> >::expose("LimitedBFGSPreconditioner"); + + } +} // namespace eigenpy diff --git a/src/solvers/preconditioners.hpp b/src/solvers/preconditioners.hpp new file mode 100644 index 0000000000000000000000000000000000000000..111a9564058b2efced35fd333a2b7a029f041635 --- /dev/null +++ b/src/solvers/preconditioners.hpp @@ -0,0 +1,27 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_preconditioners_hpp__ +#define __eigenpy_preconditioners_hpp__ + +namespace eigenpy +{ + + void exposePreconditioners(); + +} // namespace eigenpy + +#endif // define __eigenpy_preconditioners_hpp__ diff --git a/src/solvers/solvers.cpp b/src/solvers/solvers.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0cca3f8c9be43ea69ad90c0b75493b3257de9d95 --- /dev/null +++ b/src/solvers/solvers.cpp @@ -0,0 +1,45 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#include "eigenpy/solvers/solvers.hpp" +#include "eigenpy/solvers/ConjugateGradient.hpp" + +#if EIGEN_VERSION_AT_LEAST(3,3,0) + #include "eigenpy/solvers/LeastSquaresConjugateGradient.hpp" +#endif + +namespace eigenpy +{ + void exposeSolvers() + { + using namespace Eigen; + ConjugateGradientVisitor< ConjugateGradient<MatrixXd,Lower|Upper> >::expose(); +#if EIGEN_VERSION_AT_LEAST(3,3,0) + LeastSquaresConjugateGradientVisitor< LeastSquaresConjugateGradient<MatrixXd, LeastSquareDiagonalPreconditionerFix<MatrixXd::Scalar> > >::expose(); +#endif + + // Conjugate gradient with limited BFGS preconditioner + ConjugateGradientVisitor< ConjugateGradient<MatrixXd,Lower|Upper,IdentityPreconditioner > >::expose("IdentityConjugateGradient"); +// ConjugateGradientVisitor< ConjugateGradient<MatrixXd,Lower|Upper,LimitedBFGSPreconditioner<double,Dynamic,Lower|Upper> > >::expose("LimitedBFGSConjugateGradient"); + + boost::python::enum_<Eigen::ComputationInfo>("ComputationInfo") + .value("Success",Eigen::Success) + .value("NumericalIssue",Eigen::NumericalIssue) + .value("NoConvergence",Eigen::NoConvergence) + .value("InvalidInput",Eigen::InvalidInput) + ; + } +} // namespace eigenpy diff --git a/src/solvers/solvers.hpp b/src/solvers/solvers.hpp new file mode 100644 index 0000000000000000000000000000000000000000..eb316f53951963295b6f4e14afd29a2c8498225f --- /dev/null +++ b/src/solvers/solvers.hpp @@ -0,0 +1,28 @@ +/* + * Copyright 2017, Justin Carpentier, LAAS-CNRS + * + * This file is part of eigenpy. + * eigenpy is free software: you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * eigenpy is distributed in the hope that it will be + * useful, but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. You should + * have received a copy of the GNU Lesser General Public License along + * with eigenpy. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef __eigenpy_solvers_hpp__ +#define __eigenpy_solvers_hpp__ + +namespace eigenpy +{ + struct SolversScope {}; + + void exposeSolvers(); + +} // namespace eigenpy + +#endif // define __eigenpy_solvers_hpp__ diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt index dcd18a04c4a263fb1770c1447fe924ea60f22239..8524a1d62587be5d9d413913fc6911fb42f3ce70 100644 --- a/unittest/CMakeLists.txt +++ b/unittest/CMakeLists.txt @@ -2,16 +2,16 @@ # Copyright (c) 2016 CNRS # # This file is part of eigenpy -# Pinocchio is free software: you can redistribute it +# eigenpy is free software: you can redistribute it # and/or modify it under the terms of the GNU Lesser General Public # License as published by the Free Software Foundation, either version # 3 of the License, or (at your option) any later version. -# Pinocchio is distributed in the hope that it will be +# eigenpy is distributed in the hope that it will be # useful, but WITHOUT ANY WARRANTY; without even the implied warranty # of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Lesser Public License for more details. You should have # received a copy of the GNU Lesser General Public License along with -# Pinocchio If not, see +# eigenpy If not, see # <http://www.gnu.org/licenses/>. MACRO(ADD_LIB_UNIT_TEST test PKGS) diff --git a/python/test_unit/test_eigenpy.py b/unittest/python/test_eigenpy.py similarity index 100% rename from python/test_unit/test_eigenpy.py rename to unittest/python/test_eigenpy.py diff --git a/python/test_unit/test_geometry.py b/unittest/python/test_geometry.py similarity index 100% rename from python/test_unit/test_geometry.py rename to unittest/python/test_geometry.py