QuadraticProblem.hh 2.94 KB
Newer Older
1 2 3 4 5 6 7
#ifndef QUADRATICPROBLEM_H
#define QUADRATICPROBLEM_H

#include <CoMISo/Config/config.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/NProblemInterface.hh>
8 9
#include <Base/Code/Quality.hh>
LOW_CODE_QUALITY_SECTION_BEGIN
10 11 12
#include <Eigen/Eigen>
#include <Eigen/Core>
#include <Eigen/Sparse>
13
LOW_CODE_QUALITY_SECTION_END
14 15


16 17 18 19
//== NAMESPACES ===============================================================

namespace COMISO {

20 21 22 23 24 25 26 27
// this problem optimizes the quadratic functional 0.5*x^T A x -x^t b + c
class QuadraticProblem : public COMISO::NProblemInterface
{
public:

  // Sparse Matrix Type
  //  typedef Eigen::DynamicSparseMatrix<double,Eigen::ColMajor> SMatrixNP;

28 29 30 31 32 33 34 35 36 37 38 39 40 41
  QuadraticProblem()
  : A_(0,0), b_(Eigen::VectorXd::Index(0)), c_(0.0)
  {

  }

  QuadraticProblem(SMatrixNP& _A, Eigen::VectorXd& _b, const double _c)
  : A_(_A), c_(_c)
  {
    if(A_.rows() != A_.cols())
      std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
    b_ = _b;
    x_ = Eigen::VectorXd::Zero(A_.cols());
  }
42 43 44 45 46 47 48 49 50 51 52


  // number of unknowns
  virtual int n_unknowns()
  {
     return A_.rows();
  }

  // initial value where the optimization should start from
  virtual void initial_x(double* _x)
  {
David Bommes's avatar
David Bommes committed
53
        for( int i=0; i<this->n_unknowns(); ++i)
54 55 56 57 58 59
            _x[i] = x_[i];
  }

  // function evaluation at location _x
  virtual double eval_f( const double* _x )
  {
60
    Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
61 62 63 64 65 66 67

    return (double)(x.transpose()*A_*x)*0.5 - (double)(x.transpose()*b_) + c_;
  }

  // gradient evaluation at location _x
  virtual void   eval_gradient( const double* _x, double*    _g)
  {
68 69
    Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
    Eigen::Map<Eigen::VectorXd> g(_g, this->n_unknowns());
70 71 72 73 74 75 76 77 78 79 80 81 82

    g = A_*x - b_;
   }

  // hessian matrix evaluation at location _x
  virtual void   eval_hessian ( const double* _x, SMatrixNP& _H)
  {
    _H = A_;
  }

  // print result
  virtual void   store_result ( const double* _x               )
  {
83
    Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
84 85 86 87 88 89 90
    x_ = x;
  }

  // get current solution
  Eigen::VectorXd& x() { return x_;}

  // advanced properties
David Bommes's avatar
David Bommes committed
91
  virtual bool   constant_hessian() const { return true; }
92

93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
  void set_A(const SMatrixNP& _A)
  {
    A_ = _A;
    if(A_.rows() != A_.cols())
        std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
    x_ = Eigen::VectorXd::Zero(A_.cols());
  }

  void set_b(const Eigen::VectorXd& _b)
  {
    b_ = _b;
  }

  void set_c( const double _c)
  {
    c_ = _c;
  }

111 112 113 114 115 116 117 118 119 120
private:

  // quadratic problem 0.5*x^T A x -x^t b + c
 SMatrixNP       A_;
 Eigen::VectorXd b_;
 double          c_;
 // current solution, which is also used as initial value
 Eigen::VectorXd x_;
};

121 122 123
//=============================================================================
} // namespace COMISO
//=============================================================================
124 125

#endif // QUADRATICPROBLEM_H