added FEM PRoblem

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@268 1355f012-dd97-4b2f-ae87-10fa9f823a57

CoMISo/NSolver/FiniteElementProblem.hh

+//=============================================================================
+//
+//  CLASS FiniteElementProblem
+//
+//=============================================================================
+
+
+#ifndef COMISO_FINITEELEMENTPROBLEM_HH
+#define COMISO_FINITEELEMENTPROBLEM_HH
+
+
+//== INCLUDES =================================================================
+
+#include <CoMISo/Config/CoMISoDefines.hh>
+#include "NProblemInterface.hh"
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO { 
+
+//== CLASS DEFINITION =========================================================
+
+	      
+/** \class FiniteElementProblem
+
+    Brief Description.
+  
+    A more elaborate description follows.
+*/
+
+
+class FiniteElementSetBase
+{
+public:
+
+  typedef Eigen::Triplet<double> Triplet;
+
+  virtual ~FiniteElementSetBase() {}
+
+  virtual double eval_f( const double* _x ) = 0;
+
+  virtual void accumulate_gradient( const double* _x , double* _g) = 0;
+
+  virtual void accumulate_hessian ( const double* _x , std::vector<Triplet>& _triplets) = 0;
+};
+
+
+// Concepts
+
+//class FiniteElementType
+//{
+//  // define dimensions of variables and constants
+//  const static int NV = 3;
+//  const static int NC = 1;
+//
+//  typedef Eigen::Matrix<size_t,NV,1> VecI;
+//  typedef Eigen::Matrix<double,NV,1> VecV;
+//  typedef Eigen::Matrix<double,NC,1> VecC;
+//  typedef Eigen::Triplet<double> Triplet;
+//
+//  inline double eval_f       (const VecV& _x, const VecC& _c) const;
+//  inline void   eval_gradient(const VecV& _x, const VecC& _c, VecV& _g) const;
+//  inline void   eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets) const;
+//};
+
+
+//class FiniteElementInstancesType
+//{
+//
+//  typedef Eigen::Matrix<double,NC,1> VecC;
+//
+//  size_t size() const;
+//
+//  size_t index( const size_t _instance, const size_t _nr) const;
+//
+//  const VecC& c( const size_t _instance) const;
+//
+//};
+
+template<class FiniteElementType>
+class FiniteElementInstancesVector
+{
+public:
+  const static int NV = FiniteElementType::NV;
+  const static int NC = FiniteElementType::NC;
+
+  typedef typename FiniteElementType::VecI VecI;
+  typedef typename FiniteElementType::VecC VecC;
+
+  void add_element(const VecI& _indices, const VecC& _constants)
+  {
+    indices_.push_back(_indices);
+    constants_.push_back( _constants);
+  }
+
+  void clear_elements()
+  {
+    indices_.clear();
+    constants_.clear();
+  }
+
+  size_t size() const
+  {
+    return indices_.size();
+  }
+
+  size_t index( const size_t _instance, const size_t _nr) const
+  {
+    return indices_[_instance][_nr];
+  }
+
+  const VecC& c( const size_t _instance) const
+  {
+    return constants_[_instance];
+  }
+
+private:
+  std::vector<VecI> indices_;
+  std::vector<VecC> constants_;
+};
+
+
+template<class FiniteElementType, class FiniteElementInstancesType = FiniteElementInstancesVector<FiniteElementType> >
+class FiniteElementSet : public FiniteElementSetBase
+{
+public:
+
+  // export dimensions of element
+  const static int NV = FiniteElementType::NV;
+  const static int NC = FiniteElementType::NC;
+
+  typedef typename FiniteElementType::VecI VecI;
+  typedef typename FiniteElementType::VecV VecV;
+  typedef typename FiniteElementType::VecC VecC;
+
+  // access element for setting constants per element etc.
+  FiniteElementType& element() { return element_;}
+
+  // access instances for adding etc.
+  FiniteElementInstancesType& instances() { return instances_;}
+
+  virtual double eval_f( const double* _x )
+  {
+    double f(0.0);
+    for(unsigned int i=0; i<instances_.size(); ++i)
+    {
+      // get local x vector
+      for(unsigned int j=0; j<NV; ++j)
+        x_[j] = _x[instances_.index(i,j)];
+
+      f += element_.eval_f(x_, instances_.c(i));
+    }
+    return f;
+  }
+
+  virtual void accumulate_gradient( const double* _x , double* _g)
+  {
+    for(unsigned int i=0; i<instances_.size(); ++i)
+    {
+      // get local x vector
+      for(unsigned int j=0; j<NV; ++j)
+        x_[j] = _x[instances_.index(i,j)];
+
+      element_.eval_gradient(x_, instances_.c(i), g_);
+
+      // accumulate into global gradient
+      for(unsigned int j=0; j<NV; ++j)
+        _g[instances_.index(i,j)] += g_[j];
+    }
+  }
+
+  virtual void accumulate_hessian ( const double* _x , std::vector<Triplet>& _triplets)
+  {
+    for(unsigned int i=0; i<instances_.size(); ++i)
+    {
+      // get local x vector
+      for(unsigned int j=0; j<NV; ++j)
+        x_[j] = _x[instances_.index(i,j)];
+
+      element_.eval_hessian(x_, instances_.c(i), triplets_);
+
+      for(unsigned int j=0; j<triplets_.size(); ++j)
+      {
+        // add re-indexed Triplet
+        _triplets.push_back(Triplet( instances_.index(i,triplets_[j].row()),
+                                     instances_.index(i,triplets_[j].col()),
+                                     triplets_[j].value()                   ));
+      }
+    }
+  }
+
+private:
+
+  FiniteElementType element_;
+
+  FiniteElementInstancesType instances_;
+
+  VecV x_;
+  VecV g_;
+
+  std::vector<Triplet> triplets_;
+};
+
+
+class COMISODLLEXPORT FiniteElementProblem : public NProblemInterface
+{
+public:
+
+//  typedef Eigen::SparseMatrix<double,Eigen::ColMajor> SMatrixNP;
+
+  typedef FiniteElementSetBase::Triplet Triplet;
+
+  /// Default constructor
+  FiniteElementProblem(const unsigned int _n) : NProblemInterface(), n_(_n), x_(_n,0.0)
+  {}
+
+  /// Destructor
+  virtual ~FiniteElementProblem()
+  {}
+
+  void add_set(FiniteElementSetBase* _fe_set)
+  {
+    fe_sets_.push_back(_fe_set);
+  }
+
+  void clear_sets()
+  {
+    fe_sets_.clear();
+  }
+
+  std::vector<double>& x() {return x_;}
+
+  // problem definition
+  virtual int    n_unknowns   (                                )
+  {
+    return n_;
+  }
+
+  virtual void   initial_x    (       double* _x               )
+  {
+    if(n_ > 0)
+      memcpy(_x, &(x_[0]), n_*sizeof(double));
+  }
+
+  virtual double eval_f       ( const double* _x )
+  {
+    double f(0.0);
+
+    for(unsigned int i=0; i<fe_sets_.size(); ++i)
+      f += fe_sets_[i]->eval_f(_x);
+
+    return f;
+  }
+
+  virtual void   eval_gradient( const double* _x, double*    _g)
+  {
+    // clear gradient (assume floating point 0 has only zero bits)
+    memset(_g, 0, n_*sizeof(double));
+
+    for(unsigned int i=0; i<fe_sets_.size(); ++i)
+      fe_sets_[i]->accumulate_gradient(_x, _g);
+  }
+
+  virtual void   eval_hessian ( const double* _x, SMatrixNP& _H)
+  {
+    triplets_.clear();
+
+    for(unsigned int i=0; i<fe_sets_.size(); ++i)
+      fe_sets_[i]->accumulate_hessian(_x, triplets_);
+
+    // set data
+    _H.resize(n_unknowns(), n_unknowns());
+    _H.setFromTriplets(triplets_.begin(), triplets_.end());
+  }
+
+
+  virtual void   store_result ( const double* _x )
+  {
+    if(n_ > 0)
+      memcpy(&(x_[0]), _x, n_*sizeof(double));
+  }
+
+  // advanced properties (ToDo better handling)
+  virtual bool   constant_gradient() const { return false; }
+  virtual bool   constant_hessian()  const { return false; }
+
+private:
+
+  // number of unknowns
+  unsigned int n_;
+
+  // current/initial solution
+  std::vector<double> x_;
+
+  // finite element sets
+  std::vector<FiniteElementSetBase*> fe_sets_;
+
+  // vector of triplets (avoid re-allocation)
+  std::vector<Triplet> triplets_;
+};
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // COMISO_FINITEELEMENTPROBLEM_HH defined
+//=============================================================================
+

CoMISo/NSolver/IPOPTSolver.cc

@@ -45,6 +45,7 @@ IPOPTSolver()
   //  app->Options()->SetIntegerValue("print_level", 0);
   //  app->Options()->SetStringValue("expect_infeasible_problem", "yes");
 
+  print_level_ = 5;
 }
 
 
@@ -56,35 +57,54 @@ int
 IPOPTSolver::
 solve(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints)
 {
+  //----------------------------------------------------------------------------
+  // 0. Check whether hessian_approximation is active
+  //----------------------------------------------------------------------------
+  bool hessian_approximation = false;
+  std::string ha, p;
+  app().Options()->GetStringValue("hessian_approximation", ha, p);
+  if(ha != "exact")
+  {
+    if(print_level_>=2)
+      std::cerr << "Hessian approximation is enabled" << std::endl;
+    hessian_approximation = true;
+  }
+
   //----------------------------------------------------------------------------
   // 1. Create an instance of IPOPT NLP
   //----------------------------------------------------------------------------
-  Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, _constraints);
+  Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, _constraints, hessian_approximation);
   NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));
 
   //----------------------------------------------------------------------------
   // 2. exploit special characteristics of problem
   //----------------------------------------------------------------------------
 
-  std::cerr << "exploit detected special properties: ";
+  if(print_level_>=2)
+    std::cerr << "exploit detected special properties: ";
   if(np2->hessian_constant())
   {
-    std::cerr << "*constant hessian* ";
+    if(print_level_>=2)
+      std::cerr << "*constant hessian* ";
     app().Options()->SetStringValue("hessian_constant", "yes");
   }
 
   if(np2->jac_c_constant())
   {
-    std::cerr << "*constant jacobian of equality constraints* ";
+    if(print_level_>=2)
+      std::cerr << "*constant jacobian of equality constraints* ";
     app().Options()->SetStringValue("jac_c_constant", "yes");
   }
 
   if(np2->jac_d_constant())
   {
-    std::cerr << "*constant jacobian of in-equality constraints*";
+    if(print_level_>=2)
+      std::cerr << "*constant jacobian of in-equality constraints*";
     app().Options()->SetStringValue("jac_d_constant", "yes");
   }
-  std::cerr << std::endl;
+
+  if(print_level_>=2)
+    std::cerr << std::endl;
 
   //----------------------------------------------------------------------------
   // 3. solve problem
@@ -95,7 +115,8 @@ solve(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _co
   status = app_->Initialize();
   if (status != Ipopt::Solve_Succeeded)
   {
-    printf("\n\n*** Error IPOPT during initialization!\n");
+    if(print_level_>=2)
+      printf("\n\n*** Error IPOPT during initialization!\n");
   }
 
   status = app_->OptimizeTNLP( np);
@@ -103,15 +124,16 @@ solve(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _co
   //----------------------------------------------------------------------------
   // 4. output statistics
   //----------------------------------------------------------------------------
-  if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
-  {
-    // Retrieve some statistics about the solve
-    Ipopt::Index iter_count = app_->Statistics()->IterationCount();
-    printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);
+  if(print_level_>=2)
+    if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level)
+    {
+      // Retrieve some statistics about the solve
+      Ipopt::Index iter_count = app_->Statistics()->IterationCount();
+      printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count);
 
-    Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
-    printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
-  }
+      Ipopt::Number final_obj = app_->Statistics()->FinalObjective();
+      printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj);
+    }
 
   return status;
 }
@@ -129,6 +151,18 @@ solve(NProblemInterface*                        _problem,
       const double                              _almost_infeasible,
       const int                                 _max_passes        )
 {
+  //----------------------------------------------------------------------------
+  // 0. Check whether hessian_approximation is active
+  //----------------------------------------------------------------------------
+  bool hessian_approximation = false;
+  std::string ha, p;
+  app().Options()->GetStringValue("hessian_approximation", ha, p);
+  if(ha != "exact")
+  {
+    std::cerr << "Hessian approximation is enabled" << std::endl;
+    hessian_approximation = true;
+  }
+
   //----------------------------------------------------------------------------
   // 0. Initialize IPOPT Applicaiton
   //----------------------------------------------------------------------------
@@ -160,7 +194,7 @@ solve(NProblemInterface*                        _problem,
     //----------------------------------------------------------------------------
     // 1. Create an instance of current IPOPT NLP
     //----------------------------------------------------------------------------
-    Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, constraints);
+    Ipopt::SmartPtr<Ipopt::TNLP> np = new NProblemIPOPT(_problem, constraints, hessian_approximation);
     NProblemIPOPT* np2 = dynamic_cast<NProblemIPOPT*> (Ipopt::GetRawPtr(np));
     // enable caching of solution
     np2->store_solution() = true;
@@ -469,17 +503,21 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
   std::vector<double> x(n);
   problem_->initial_x(P(x));
 
+
   // nonzeros in the jacobian of C_ and the hessian of the lagrangian
   SMatrixNP HP;
   SVectorNC g;
   SMatrixNC H;
-  problem_->eval_hessian(P(x), HP);
+  if(!hessian_approximation_)
+  {
+    problem_->eval_hessian(P(x), HP);
 
-  // get nonzero structure of hessian of problem
-  for(int i=0; i<HP.outerSize(); ++i)
-    for (SMatrixNP::InnerIterator it(HP,i); it; ++it)
-      if(it.row() >= it.col())
-        ++nnz_h_lag;
+    // get nonzero structure of hessian of problem
+    for(int i=0; i<HP.outerSize(); ++i)
+      for (SMatrixNP::InnerIterator it(HP,i); it; ++it)
+        if(it.row() >= it.col())
+          ++nnz_h_lag;
+  }
 
   // get nonzero structure of constraints
   for( int i=0; i<m; ++i)
@@ -488,13 +526,16 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
 
     nnz_jac_g += g.nonZeros();
 
-    // count lower triangular elements
-    constraints_[i]->eval_hessian (P(x),H);
+    if(!hessian_approximation_)
+    {
+      // count lower triangular elements
+      constraints_[i]->eval_hessian (P(x),H);
 
-    SMatrixNC::iterator m_it = H.begin();
-    for(; m_it != H.end(); ++m_it)
-      if( m_it.row() >= m_it.col())
-        ++nnz_h_lag;
+      SMatrixNC::iterator m_it = H.begin();
+      for(; m_it != H.end(); ++m_it)
+        if( m_it.row() >= m_it.col())
+          ++nnz_h_lag;
+    }
   }
 
   // We use the standard fortran index style for row/col entries

CoMISo/NSolver/IPOPTSolver.hh

@@ -106,6 +106,12 @@ public:
   Ipopt::IpoptApplication& app() {return (*app_); }
 
 
+  void set_print_level(const int _level)
+  {
+    app().Options()->SetIntegerValue("print_level", _level);
+    print_level_ = _level;
+  }
+
 protected:
   double* P(std::vector<double>& _v)
   {
@@ -119,6 +125,8 @@ private:
 
   // smart pointer to IpoptApplication to set options etc.
   Ipopt::SmartPtr<Ipopt::IpoptApplication> app_;
+
+  int print_level_;
 };
 
 
@@ -142,8 +150,8 @@ public:
   typedef NProblemInterface::SMatrixNP    SMatrixNP;
 
   /** default constructor */
-  NProblemIPOPT(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints)
-   : problem_(_problem), store_solution_(false) { split_constraints(_constraints); analyze_special_properties(_problem, _constraints);}
+  NProblemIPOPT(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints, const bool _hessian_approximation = false)
+   : problem_(_problem), store_solution_(false), hessian_approximation_(_hessian_approximation) { split_constraints(_constraints); analyze_special_properties(_problem, _constraints);}
 
   /** default destructor */
   virtual ~NProblemIPOPT() {};
@@ -258,6 +266,8 @@ private:
 
   bool store_solution_;
   std::vector<double> x_;
+
+  bool hessian_approximation_;
 };
 
 

CoMISo/NSolver/NPTiming.cc

@@ -104,16 +104,20 @@ void NPTiming::print_statistics()
   std::cerr << "total time    : " << time_total/1000.0 << "s\n";
   std::cerr << "total time NP : " << time_np/1000.0 << "s  (" << time_np/time_total*100.0 << " %)\n";
 
+  double timing_eval_f_avg        = timing_eval_f_/double(n_eval_f_);
+  double timing_eval_gradient_avg = timing_eval_gradient_/double(n_eval_gradient_);
+  double timing_eval_hessian_avg  = timing_eval_hessian_/double(n_eval_hessian_);
+
   std::cerr << std::fixed << std::setprecision(5)
   << "eval_f time   : " << timing_eval_f_/1000.0
   << "s  ( #evals: " << n_eval_f_ << " -> avg "
-  << timing_eval_f_/(1000.0*double(n_eval_f_)) << "s )\n"
+  << timing_eval_f_avg/1000.0 << "s )\n"
   << "eval_grad time: " << timing_eval_gradient_/1000.0
   << "s  ( #evals: " << n_eval_gradient_ << " -> avg "
-  << timing_eval_gradient_/(1000.0*double(n_eval_gradient_)) << "s )\n"
+  << timing_eval_gradient_avg/1000.0 << "s, factor: " <<  timing_eval_gradient_avg / timing_eval_f_avg << ")\n"
   << "eval_hess time: " << timing_eval_hessian_/1000.0
   << "s  ( #evals: " << n_eval_hessian_ << " -> avg "
-  << timing_eval_hessian_/(1000.0*double(n_eval_hessian_)) << "s )\n";
+  << timing_eval_hessian_avg/1000.0 << "s, factor: " <<  timing_eval_hessian_avg / timing_eval_f_avg << ")\n";
 }
 
 

CoMISo/NSolver/NProblemInterfaceADOLC.hh

@@ -122,6 +122,7 @@ public:
       }
       else
       {
+        std::cerr << "re-tape..." << std::endl;
         adouble y_d = 0.0;
 
         adouble* xa = new adouble[this->n_unknowns()];

CoMISo/NSolver/NProblemInterfaceDCO.hh

+/*
+ * NProblemInterfaceADOLC.hh
+ *
+ *  Created on: Feb 20, 2013
+ *      Author: bommes / kremer
+ */
+
+#ifndef NPROBLEMINTERFACEADOLC_HH
+#define NPROBLEMINTERFACEADOLC_HH
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include <CoMISo/Config/config.hh>
+#if COMISO_ADOLC_AVAILABLE
+#if COMISO_EIGEN3_AVAILABLE
+
+//== INCLUDES =================================================================
+
+#include <vector>
+
+#include "NProblemInterface.hh"
+
+#include <CoMISo/Config/CoMISoDefines.hh>
+
+#include <dco.hpp>
+//#define USE_DCO_VECTOR
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+/** \class NProblemmInterfaceAD NProblemInterfaceADOC.hh <ACG/.../NProblemInterfaceADOLC.hh>
+
+ Brief Description.
+
+ Extend the problem interface with auto differentiation using ADOL-C
+ */
+  class COMISODLLEXPORT NProblemInterfaceDCO : public NProblemInterface
+  {
+  public:
+#ifdef USE_DCO_VECTOR
+    typedef dco::gt1v<double, 7> MODE_grad;
+    typedef MODE_grad::type ad_t_grad;
+#else
+    typedef dco::ga1s<double> MODE_grad