[Solvers] Expose preconditioners

CMakeLists.txt

@@ -82,10 +82,13 @@ INCLUDE_DIRECTORIES(${NUMPY_INCLUDE_DIRS})
 # ----------------------------------------------------
 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
@@ -139,6 +142,7 @@ ENDFOREACH(header)
 # --- TARGETS ----------------------------------------
 # ----------------------------------------------------
 SET(${PROJECT_NAME}_SOLVERS_SOURCES
+  src/solvers/preconditioners.cpp
   src/solvers/solvers.cpp
   )
 SET(${PROJECT_NAME}_SOURCES

python/main.cpp

@@ -18,6 +18,7 @@
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/geometry.hpp"
 #include "eigenpy/solvers/solvers.hpp"
+#include "eigenpy/solvers/preconditioners.hpp"
 
 #include <iostream>
 
@@ -28,5 +29,6 @@ BOOST_PYTHON_MODULE(eigenpy)
   eigenpy::exposeAngleAxis();
   eigenpy::exposeQuaternion();
   eigenpy::exposeSolvers();
+  eigenpy::exposePreconditioners();
   
 }

src/solvers/BFGSPreconditioners.hpp

+/*
+ * 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__

src/solvers/BasicPreconditioners.hpp

+/*
+ * 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."))
+      .def("info",&Preconditioner::info,
+           "Returns success if the Preconditioner has been well initialized.")
+      .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>())
+      ;
+      
+    }
+    
+  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)
+      ;
+      
+    }
+  };
+  
+  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>())
+      ;
+      
+    }
+  };
+  
+  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__

src/solvers/IterativeSolverBase.hpp

@@ -33,6 +33,7 @@ namespace eigenpy
   {
     
     typedef typename IterativeSolver::MatrixType MatrixType;
+    typedef typename IterativeSolver::Preconditioner Preconditioner;
     typedef Eigen::VectorXd VectorType;
     
     template<class PyClass>
@@ -68,6 +69,7 @@ namespace eigenpy
            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<>())
       ;
       
     }

src/solvers/preconditioners.cpp

+/*
+ * 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();
+    LeastSquareDiagonalPreconditionerVisitor<double>::expose();
+    IdentityPreconditionerVisitor::expose();
+    LimitedBFGSPreconditionerBaseVisitor< LimitedBFGSPreconditioner<double,Eigen::Dynamic,Eigen::Upper|Eigen::Lower> >::expose("LimitedBFGSPreconditioner");
+   
+    
+  }
+} // namespace eigenpy

src/solvers/preconditioners.hpp

+/*
+ * 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__