Commit 21a8939e authored by David Bommes's avatar David Bommes

extended LeastSquaresProblem

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@294 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent fbbccb71
...@@ -142,10 +142,108 @@ void ...@@ -142,10 +142,108 @@ void
LeastSquaresProblem:: LeastSquaresProblem::
eval_hessian ( const double* _x, SMatrixNP& _H) eval_hessian ( const double* _x, SMatrixNP& _H)
{ {
// clear old data std::vector<Triplet> triplets;
triplets.reserve(10*n_unknowns());
add_hessian_coeffs(_x, triplets, 1.0);
_H.resize(n_unknowns(), n_unknowns()); _H.resize(n_unknowns(), n_unknowns());
_H.setZero(); _H.setFromTriplets(triplets.begin(), triplets.end());
// // clear old data
// _H.resize(n_unknowns(), n_unknowns());
// _H.setZero();
// for(unsigned int i=0; i<terms_.size(); ++i)
// {
// // get local function value
// double vterm = terms_[i]->eval_constraint(_x);
//
// // get local gradient
// NConstraintInterface::SVectorNC gterm;
// terms_[i]->eval_gradient(_x, gterm);
//
// // add terms to global gradient
// NConstraintInterface::SVectorNC::InnerIterator v_it (gterm);
// for( ; v_it; ++v_it)
// {
// NConstraintInterface::SVectorNC::InnerIterator v_it2(gterm);
// for( ; v_it2; ++v_it2)
// _H.coeffRef(v_it.index(), v_it2.index()) += 2.0*v_it.value()*v_it2.value();
// }
//
// NConstraintInterface::SMatrixNC Hterm;
// terms_[i]->eval_hessian(_x, Hterm);
//
// NConstraintInterface::SMatrixNC::iterator h_it = Hterm.begin();
// NConstraintInterface::SMatrixNC::iterator h_end = Hterm.end();
// for(; h_it != h_end; ++h_it)
// _H.coeffRef(h_it.row(),h_it.col()) += 2.0*vterm*(*h_it);
// }
}
//-----------------------------------------------------------------------------
void
LeastSquaresProblem::
store_result ( const double* _x )
{
for( int i=0; i<this->n_unknowns(); ++i)
x_[i] = _x[i];
}
//-----------------------------------------------------------------------------
bool
LeastSquaresProblem::
constant_hessian()
{
for(unsigned int i=0; i<terms_.size(); ++i)
{
if(!terms_[i]->is_linear())
return false;
}
return true;
}
//-----------------------------------------------------------------------------
void
LeastSquaresProblem::
add_to_gradient (const double* _x, double* _g, const double _c)
{
for(unsigned int i=0; i<terms_.size(); ++i)
{
// get local function value
double vterm = terms_[i]->eval_constraint(_x);
// get local gradient
NConstraintInterface::SVectorNC gterm;
terms_[i]->eval_gradient(_x, gterm);
// add terms to global gradient
NConstraintInterface::SVectorNC::InnerIterator v_it(gterm);
for( ; v_it; ++v_it)
{
_g[v_it.index()] += _c*2.0*vterm*v_it.value();
}
}
}
//-----------------------------------------------------------------------------
void
LeastSquaresProblem::
add_hessian_coeffs(const double* _x, std::vector<Triplet>& _trips, const double _c)
{
for(unsigned int i=0; i<terms_.size(); ++i) for(unsigned int i=0; i<terms_.size(); ++i)
{ {
// get local function value // get local function value
...@@ -161,7 +259,8 @@ eval_hessian ( const double* _x, SMatrixNP& _H) ...@@ -161,7 +259,8 @@ eval_hessian ( const double* _x, SMatrixNP& _H)
{ {
NConstraintInterface::SVectorNC::InnerIterator v_it2(gterm); NConstraintInterface::SVectorNC::InnerIterator v_it2(gterm);
for( ; v_it2; ++v_it2) for( ; v_it2; ++v_it2)
_H.coeffRef(v_it.index(), v_it2.index()) += 2.0*v_it.value()*v_it2.value(); _trips.push_back(Triplet(v_it.index(), v_it2.index(), _c*2.0*v_it.value()*v_it2.value()));
// _H.coeffRef(v_it.index(), v_it2.index()) += 2.0*v_it.value()*v_it2.value();
} }
NConstraintInterface::SMatrixNC Hterm; NConstraintInterface::SMatrixNC Hterm;
...@@ -170,7 +269,8 @@ eval_hessian ( const double* _x, SMatrixNP& _H) ...@@ -170,7 +269,8 @@ eval_hessian ( const double* _x, SMatrixNP& _H)
NConstraintInterface::SMatrixNC::iterator h_it = Hterm.begin(); NConstraintInterface::SMatrixNC::iterator h_it = Hterm.begin();
NConstraintInterface::SMatrixNC::iterator h_end = Hterm.end(); NConstraintInterface::SMatrixNC::iterator h_end = Hterm.end();
for(; h_it != h_end; ++h_it) for(; h_it != h_end; ++h_it)
_H.coeffRef(h_it.row(),h_it.col()) += 2.0*vterm*(*h_it); _trips.push_back( Triplet(h_it.row(), h_it.col(), _c*2.0*vterm*(*h_it)));
// _H.coeffRef(h_it.row(),h_it.col()) += 2.0*vterm*(*h_it);
} }
} }
...@@ -178,35 +278,21 @@ eval_hessian ( const double* _x, SMatrixNP& _H) ...@@ -178,35 +278,21 @@ eval_hessian ( const double* _x, SMatrixNP& _H)
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
void double
LeastSquaresProblem:: LeastSquaresProblem::
store_result ( const double* _x ) max_deviaton( const double* _x )
{ {
for( int i=0; i<this->n_unknowns(); ++i) double vmax(0.0);
x_[i] = _x[i];
}
//-----------------------------------------------------------------------------
bool
LeastSquaresProblem::
constant_hessian()
{
for(unsigned int i=0; i<terms_.size(); ++i) for(unsigned int i=0; i<terms_.size(); ++i)
{ {
if(!terms_[i]->is_linear()) double vterm = std::abs(terms_[i]->eval_constraint(_x));
return false; if(vterm > vmax) vmax = vterm;
} }
return true; return vmax;
} }
//-----------------------------------------------------------------------------
//============================================================================= //=============================================================================
} // namespace COMISO } // namespace COMISO
//============================================================================= //=============================================================================
...@@ -35,6 +35,9 @@ namespace COMISO { ...@@ -35,6 +35,9 @@ namespace COMISO {
class COMISODLLEXPORT LeastSquaresProblem : public NProblemInterface class COMISODLLEXPORT LeastSquaresProblem : public NProblemInterface
{ {
public: public:
// Eigen Triplet for Hessian Accumulation
typedef Eigen::Triplet<double> Triplet;
/// Default constructor /// Default constructor
LeastSquaresProblem(const int _n_unknowns) :n_(_n_unknowns), x_(_n_unknowns, 0.0) {} LeastSquaresProblem(const int _n_unknowns) :n_(_n_unknowns), x_(_n_unknowns, 0.0) {}
...@@ -63,6 +66,11 @@ public: ...@@ -63,6 +66,11 @@ public:
// advanced properties // advanced properties
virtual bool constant_hessian(); virtual bool constant_hessian();
// advanced usage
void add_to_gradient ( const double* _x, double* _g, const double _c);
void add_hessian_coeffs( const double* _x, std::vector<Triplet>& _trips, const double _c);
double max_deviaton ( const double* _x );
private: private:
// #unknowns // #unknowns
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment