Commit 49aef7a2 authored by Mike Kremer's avatar Mike Kremer
Browse files

- Added error epsilon to be passed as optional construction parameter to...

- Added error epsilon to be passed as optional construction parameter to classes derived from NConstraintInterface.
- Fixed taping issue when using ADOL-C driven auto differentiation in IpOpt. Now the taping works just fine.

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@203 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 356a0c1e
......@@ -50,7 +50,7 @@ public:
enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL};
/// Default constructor
NConstraintInterface(const ConstraintType _type = NC_EQUAL) : type_(_type) {}
NConstraintInterface(const ConstraintType _type = NC_EQUAL, double _eps = 1e-6) : type_(_type), eps_(_eps) {}
/// Destructor
virtual ~NConstraintInterface() {}
......@@ -62,13 +62,13 @@ public:
virtual ConstraintType constraint_type ( ) { return type_; }
virtual bool is_satisfied ( const double* _x, double _eps = 1e-6 )
virtual bool is_satisfied ( const double* _x )
{
switch( type_)
{
case NC_EQUAL : return (fabs(eval_constraint(_x)) <= _eps); break;
case NC_LESS_EQUAL : return ( eval_constraint(_x) <= _eps); break;
case NC_GREATER_EQUAL: return ( eval_constraint(_x) >= -_eps); break;
case NC_EQUAL : return (fabs(eval_constraint(_x)) <= eps_); break;
case NC_LESS_EQUAL : return ( eval_constraint(_x) <= eps_); break;
case NC_GREATER_EQUAL: return ( eval_constraint(_x) >= -eps_); break;
}
return false;
}
......@@ -78,11 +78,11 @@ public:
virtual bool constant_gradient() const { return false;}
virtual bool constant_hessian () const { return false;}
virtual double gradient_update_factor( const double* _x, double _eps = 1e-6 )
virtual double gradient_update_factor( const double* _x )
{
double val = eval_constraint(_x);
bool upper_bound_ok = ( val <= _eps);
bool lower_bound_ok = ( val >= -_eps);
bool upper_bound_ok = ( val <= eps_);
bool lower_bound_ok = ( val >= -eps_);
if(upper_bound_ok)
{
......@@ -100,6 +100,8 @@ public:
private:
// constraint type
ConstraintType type_;
double eps_;
};
......
......@@ -28,6 +28,8 @@
#include "NConstraintInterface.hh"
#include "NProblemInterfaceAD.hpp"
#include "TapeIDSingleton.hh"
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
......@@ -53,16 +55,18 @@ public:
typedef NConstraintInterface::ConstraintType ConstraintType;
/// Default constructor
NConstraintInterfaceAD(NProblemInterfaceAD& _problem, int _n_unknowns, const ConstraintType _type = NC_EQUAL) :
NConstraintInterface(_type),
NConstraintInterfaceAD(NProblemInterfaceAD& _problem, int _n_unknowns,
const ConstraintType _type = NC_EQUAL, double _eps = 1e-6) :
NConstraintInterface(_type, _eps),
problem_(_problem),
n_unknowns_(_n_unknowns),
type_(_type),
function_evaluated_(false),
use_tape_(true),
constant_hessian_evaluated_(false) {
constant_hessian_evaluated_(false),
tape_(TapeIDSingleton::Instance()->uniqueTapeID()) {
for(size_t i = 0; i < 11; ++i) tape_stats_[i] = -1;
for(size_t i = 0; i < 11; ++i) tape_stats_[i] = 0;
}
/// Destructor
......@@ -90,7 +94,7 @@ public:
boost::shared_array<adouble> x_d_ptr = problem_.x_d_ptr();
trace_on(1); // Start taping
trace_on(tape_); // Start taping
// Fill data vector
for(int i = 0; i < n_unknowns_; ++i) {
......@@ -105,9 +109,23 @@ public:
trace_off();
#ifndef NDEBUG
tapestats(1, tape_stats_);
// Do some status output here...
#ifdef ADOLC_STATS
tapestats(tape_, tape_stats_);
std::cout << "Status values for tape " << tape_ << std::endl;
std::cout << "===============================================" << std::endl;
std::cout << "Number of independent variables:\t" << tape_stats_[0] << std::endl;
std::cout << "Number of dependent variables:\t\t" << tape_stats_[1] << std::endl;
std::cout << "Max. number of live active variables:\t" << tape_stats_[2] << std::endl;
std::cout << "Size of value stack:\t\t\t" << tape_stats_[3] << std::endl;
std::cout << "Buffer size:\t\t\t\t" << tape_stats_[4] << std::endl;
std::cout << "Total number of operations recorded:\t" << tape_stats_[5] << std::endl;
std::cout << "Other stats [6]:\t\t\t" << tape_stats_[6] << std::endl;
std::cout << "Other stats [7]:\t\t\t" << tape_stats_[7] << std::endl;
std::cout << "Other stats [8]:\t\t\t" << tape_stats_[8] << std::endl;
std::cout << "Other stats [9]:\t\t\t" << tape_stats_[9] << std::endl;
std::cout << "Other stats [10]:\t\t\t" << tape_stats_[10] << std::endl;
std::cout << "Other stats [11]:\t\t\t" << tape_stats_[11] << std::endl;
std::cout << "===============================================" << std::endl;
#endif
function_evaluated_ = true;
......@@ -116,9 +134,9 @@ public:
double ay[1] = {0.0};
int ec = function(1, 1, n_unknowns_, const_cast<double*>(_x), ay);
int ec = function(tape_, 1, n_unknowns_, const_cast<double*>(_x), ay);
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: function() returned code " << ec << std::endl;
#endif
......@@ -140,9 +158,17 @@ public:
_g.resize(n_unknowns_);
_g.setZero();
int ec = gradient(1, n_unknowns_, _x, grad_p.get());
int ec = gradient(tape_, n_unknowns_, _x, grad_p.get());
if(ec <= 0) {
// Retape function if return code indicates discontinuity
function_evaluated_ = false;
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
eval_constraint(_x);
ec = gradient(tape_, n_unknowns_, _x, grad_p.get());
}
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: gradient() returned code " << ec << std::endl;
#endif
......@@ -174,14 +200,21 @@ public:
unsigned int* c_ind = NULL;
double* val = NULL;
int ec = sparse_hess(1, n_unknowns_, 0, _x, &nz, &r_ind, &c_ind, &val, opt);
int ec = sparse_hess(tape_, n_unknowns_, 0, _x, &nz, &r_ind, &c_ind, &val, opt);
if(ec <= 0) {
// Retape function if return code indicates discontinuity
function_evaluated_ = false;
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
eval_constraint(_x);
ec = sparse_hess(tape_, n_unknowns_, 0, _x, &nz, &r_ind, &c_ind, &val, opt);
}
assert(*nz >= 0);
assert(r_ind != NULL);
assert(c_ind != NULL);
assert(val != NULL);
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: sparse_hessian() returned code " << ec << std::endl;
#endif
......@@ -203,9 +236,17 @@ public:
double** h_ptr = problem_.dense_hessian_ptr();
int ec = hessian(1, n_unknowns_, const_cast<double*>(_x), h_ptr);
int ec = hessian(tape_, n_unknowns_, const_cast<double*>(_x), h_ptr);
if(ec <= 0) {
// Retape function if return code indicates discontinuity
function_evaluated_ = false;
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
eval_constraint(_x);
ec = hessian(tape_, n_unknowns_, const_cast<double*>(_x), h_ptr);
}
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: hessian() returned code " << ec << std::endl;
#endif
......@@ -260,13 +301,15 @@ private:
// Constraint type
ConstraintType type_;
int tape_stats_[11];
size_t tape_stats_[11];
bool function_evaluated_;
bool use_tape_;
SMatrixNC constant_hessian_;
bool constant_hessian_evaluated_;
const short int tape_;
};
//=============================================================================
......
......@@ -33,6 +33,8 @@
#include <CoMISo/Config/CoMISoDefines.hh>
#include "TapeIDSingleton.hh"
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
......@@ -63,9 +65,10 @@ public:
n_unknowns_(_n_unknowns),
dense_hessian_(NULL),
function_evaluated_(false),
use_tape_(true) {
use_tape_(true),
tape_(TapeIDSingleton::Instance()->uniqueTapeID()) {
for(size_t i = 0; i < 11; ++i) tape_stats_[i] = -1;
for(size_t i = 0; i < 11; ++i) tape_stats_[i] = 0;
}
/// Destructor
......@@ -119,7 +122,7 @@ public:
boost::shared_array<adouble> x_d = x_d_ptr();
trace_on(1); // Start taping
trace_on(tape_); // Start taping
// Fill data vector
for(int i = 0; i < n_unknowns_; ++i) {
......@@ -134,10 +137,24 @@ public:
trace_off();
#ifndef NDEBUG
tapestats(1, tape_stats_);
// Do some status output here...
#endif
#ifdef ADOLC_STATS
tapestats(tape_, tape_stats_);
std::cout << "Status values for tape " << tape_ << std::endl;
std::cout << "===============================================" << std::endl;
std::cout << "Number of independent variables:\t" << tape_stats_[0] << std::endl;
std::cout << "Number of dependent variables:\t\t" << tape_stats_[1] << std::endl;
std::cout << "Max. number of live active variables:\t" << tape_stats_[2] << std::endl;
std::cout << "Size of value stack:\t\t\t" << tape_stats_[3] << std::endl;
std::cout << "Buffer size:\t\t\t\t" << tape_stats_[4] << std::endl;
std::cout << "Total number of operations recorded:\t" << tape_stats_[5] << std::endl;
std::cout << "Other stats [6]:\t\t\t" << tape_stats_[6] << std::endl;
std::cout << "Other stats [7]:\t\t\t" << tape_stats_[7] << std::endl;
std::cout << "Other stats [8]:\t\t\t" << tape_stats_[8] << std::endl;
std::cout << "Other stats [9]:\t\t\t" << tape_stats_[9] << std::endl;
std::cout << "Other stats [10]:\t\t\t" << tape_stats_[10] << std::endl;
std::cout << "Other stats [11]:\t\t\t" << tape_stats_[11] << std::endl;
std::cout << "===============================================" << std::endl;
#endif
function_evaluated_ = true;
......@@ -145,9 +162,9 @@ public:
double ay[1] = {0.0};
int ec = function(1, 1, n_unknowns_, const_cast<double*>(_x), ay);
int ec = function(tape_, 1, n_unknowns_, const_cast<double*>(_x), ay);
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: function() returned code " << ec << std::endl;
#endif
......@@ -164,9 +181,17 @@ public:
eval_f(_x);
}
int ec = gradient(1, n_unknowns_, _x, _g);
int ec = gradient(tape_, n_unknowns_, _x, _g);
if(ec <= 0) {
// Retape function if return code indicates discontinuity
function_evaluated_ = false;
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
eval_f(_x);
ec = gradient(tape_, n_unknowns_, _x, _g);
}
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: gradient() returned code " << ec << std::endl;
#endif
}
......@@ -190,14 +215,22 @@ public:
unsigned int* c_ind_p = NULL;
double* val_p = NULL;
int ec = sparse_hess(1, n_unknowns_, 0, _x, &nz, &r_ind_p, &c_ind_p, &val_p, opt);
int ec = sparse_hess(tape_, n_unknowns_, 0, _x, &nz, &r_ind_p, &c_ind_p, &val_p, opt);
if(ec <= 0) {
// Retape function if return code indicates discontinuity
function_evaluated_ = false;
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
eval_f(_x);
ec = sparse_hess(tape_, n_unknowns_, 0, _x, &nz, &r_ind_p, &c_ind_p, &val_p, opt);
}
assert(*nz >= 0);
assert(r_ind_p != NULL);
assert(c_ind_p != NULL);
assert(val_p != NULL);
#ifndef NDEBUG
#ifdef ADOLC_RET_CODES
std::cout << "Info: sparse_hessian() returned code " << ec << std::endl;
#endif
......@@ -215,9 +248,17 @@ public:
// Dense hessian data
double** h_ptr = dense_hessian_ptr();
int ec = hessian(1, n_unknowns_, const_cast<double*>(_x), h_ptr);
int ec = hessian(tape_, n_unknowns_, const_cast<double*>(_x), h_ptr);
#ifndef NDEBUG
if(ec <= 0) {
// Retape function if return code indicates discontinuity
function_evaluated_ = false;
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
eval_f(_x);
ec = hessian(tape_, n_unknowns_, const_cast<double*>(_x), h_ptr);
}
#ifdef ADOLC_RET_CODES
std::cout << "Info: hessian() returned code " << ec << std::endl;
#endif
......@@ -296,10 +337,12 @@ private:
// Dense hessian
double** dense_hessian_;
int tape_stats_[11];
size_t tape_stats_[11];
bool function_evaluated_;
bool use_tape_;
const short int tape_;
};
//=============================================================================
......
/*
* TapeIDSingleton.cc
*
* Created on: Jan 4, 2013
* Author: kremer
*/
#include "TapeIDSingleton.hh"
TapeIDSingleton* TapeIDSingleton::reference_ = NULL;
/*
* TapeIDSingleton.hpp
*
* Created on: Jan 4, 2013
* Author: kremer
*/
#ifndef TAPEIDSINGLETON_HPP_
#define TAPEIDSINGLETON_HPP_
#include <string>
class TapeIDSingleton {
public:
static TapeIDSingleton* Instance() {
if(reference_ == NULL) {
reference_ = new TapeIDSingleton();
}
return reference_;
}
short int uniqueTapeID() {
short int id = tape_count_;
++tape_count_;
return id;
}
private:
TapeIDSingleton() : tape_count_(0) {}
TapeIDSingleton(const TapeIDSingleton&) {}
~TapeIDSingleton() {}
static TapeIDSingleton* reference_;
short int tape_count_;
};
#endif /* TAPEIDSINGLETON_HPP_ */
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