Commit 565666af authored by Max Lyon's avatar Max Lyon
Browse files

merge from ReForm

parents 396d2622 a0e47ef4
Pipeline #5865 passed with stages
in 6 minutes and 56 seconds
......@@ -323,7 +323,10 @@ void IPOPTSolverLean::solve(
//----------------------------------------------------------------------------
// 3. solve problem
//----------------------------------------------------------------------------
status = impl_->app_->OptimizeTNLP(np);
{
DEB_time_session_def("IPOPT App OptimizeTNLP(np)");
status = impl_->app_->OptimizeTNLP(np);
}
check_ipopt_status(status);
......
......@@ -127,6 +127,21 @@ public:
IpoptCalculatedQuantities* ip_cq);
//@}
/** Intermediate Callback method for the user. Providing dummy
* default implementation. For details see IntermediateCallBack
* in IpNLP.hpp. */
virtual bool intermediate_callback(
Ipopt::AlgorithmMode mode,
Index iter, Number obj_value,
Number inf_pr, Number inf_du,
Number mu, Number d_norm,
Number regularization_size,
Number alpha_du, Number alpha_pr,
Index ls_trials,
const IpoptData* ip_data,
IpoptCalculatedQuantities* ip_cq
) override;
// special properties of problem
bool hessian_constant() const;
bool jac_c_constant() const;
......@@ -270,6 +285,22 @@ public:
IpoptCalculatedQuantities* ip_cq);
//@}
/** Intermediate Callback method for the user. Providing dummy
* default implementation. For details see IntermediateCallBack
* in IpNLP.hpp. */
virtual bool intermediate_callback(
Ipopt::AlgorithmMode mode,
Index iter, Number obj_value,
Number inf_pr, Number inf_du,
Number mu, Number d_norm,
Number regularization_size,
Number alpha_du, Number alpha_pr,
Index ls_trials,
const IpoptData* ip_data,
IpoptCalculatedQuantities* ip_cq
) override;
private:
/**@name Methods to block default compiler methods.
* The compiler automatically generates the following three methods.
......
......@@ -40,6 +40,7 @@ void
NProblemIPOPT::
split_constraints(const std::vector<NConstraintInterface*>& _constraints)
{
DEB_enter_func;
// split user-provided constraints into general-constraints and bound-constraints
constraints_ .clear(); constraints_.reserve(_constraints.size());
bound_constraints_.clear(); bound_constraints_.reserve(_constraints.size());
......@@ -100,6 +101,7 @@ analyze_special_properties(const NProblemInterface* _problem, const std::vector<
bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
Index& nnz_h_lag, IndexStyleEnum& index_style)
{
DEB_enter_func;
// number of variables
n = static_cast<Index>(problem_->n_unknowns());
......@@ -233,6 +235,7 @@ bool NProblemIPOPT::get_starting_point(Index n, bool init_x, Number* x,
Index m, bool init_lambda,
Number* lambda)
{
DEB_enter_func;
// get initial value of problem instance
problem_->initial_x(x);
......@@ -245,6 +248,7 @@ bool NProblemIPOPT::get_starting_point(Index n, bool init_x, Number* x,
bool NProblemIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
{
DEB_enter_func;
// return the value of the objective function
obj_value = problem_->eval_f(x);
return true;
......@@ -256,6 +260,7 @@ bool NProblemIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_val
bool NProblemIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
{
DEB_enter_func;
problem_->eval_gradient(x, grad_f);
return true;
......@@ -267,6 +272,7 @@ bool NProblemIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* gr
bool NProblemIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
{
DEB_enter_func;
// evaluate all constraint functions
for( int i=0; i<m; ++i)
g[i] = constraints_[i]->eval_constraint(x);
......@@ -512,6 +518,26 @@ void NProblemIPOPT::finalize_solution(SolverReturn status,
//-----------------------------------------------------------------------------
bool NProblemIPOPT::intermediate_callback(
Ipopt::AlgorithmMode /*mode*/,
Index /*iter*/, Number /*obj_value*/,
Number /*inf_pr*/, Number /*inf_du*/,
Number /*mu*/, Number /*d_norm*/,
Number /*regularization_size*/,
Number /*alpha_du*/, Number /*alpha_pr*/,
Index /*ls_trials*/,
const IpoptData* /*ip_data*/,
IpoptCalculatedQuantities* /*ip_cq*/
)
{
PROGRESS_TICK;
return true;
}
//-----------------------------------------------------------------------------
bool NProblemIPOPT::hessian_constant() const
{
return hessian_constant_;
......@@ -542,6 +568,7 @@ bool NProblemIPOPT::jac_d_constant() const
bool NProblemGmmIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
Index& nnz_h_lag, IndexStyleEnum& index_style)
{
DEB_enter_func;
// number of variables
n = problem_->n_unknowns();
......@@ -633,6 +660,7 @@ bool NProblemGmmIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
bool NProblemGmmIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u,
Index m, Number* g_l, Number* g_u)
{
DEB_enter_func;
// first clear all variable bounds
for( int i=0; i<n; ++i)
{
......@@ -668,6 +696,7 @@ bool NProblemGmmIPOPT::get_starting_point(Index n, bool init_x, Number* x,
Index m, bool init_lambda,
Number* lambda)
{
DEB_enter_func;
// get initial value of problem instance
problem_->initial_x(x);
......@@ -680,6 +709,7 @@ bool NProblemGmmIPOPT::get_starting_point(Index n, bool init_x, Number* x,
bool NProblemGmmIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
{
DEB_enter_func;
// return the value of the objective function
obj_value = problem_->eval_f(x);
return true;
......@@ -691,6 +721,7 @@ bool NProblemGmmIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_
bool NProblemGmmIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
{
DEB_enter_func;
problem_->eval_gradient(x, grad_f);
return true;
......@@ -702,6 +733,7 @@ bool NProblemGmmIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number*
bool NProblemGmmIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
{
DEB_enter_func;
// evaluate all constraint functions
for( int i=0; i<m; ++i)
g[i] = constraints_[i]->eval_constraint(x);
......@@ -840,10 +872,27 @@ void NProblemGmmIPOPT::finalize_solution(SolverReturn status,
const IpoptData* ip_data,
IpoptCalculatedQuantities* ip_cq)
{
DEB_enter_func;
// problem knows what to do
problem_->store_result(x);
}
bool NProblemGmmIPOPT::intermediate_callback(
Ipopt::AlgorithmMode /*mode*/,
Index /*iter*/, Number /*obj_value*/,
Number /*inf_pr*/, Number /*inf_du*/,
Number /*mu*/, Number /*d_norm*/,
Number /*regularization_size*/,
Number /*alpha_du*/, Number /*alpha_pr*/,
Index /*ls_trials*/,
const IpoptData* /*ip_data*/,
IpoptCalculatedQuantities* /*ip_cq*/
)
{
PROGRESS_TICK;
return true;
}
//=============================================================================
} // namespace COMISO
......
......@@ -4,7 +4,7 @@
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
......@@ -20,7 +20,7 @@
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
\*===========================================================================*/
//=============================================================================
......@@ -48,24 +48,22 @@ namespace COMISO {
//== IMPLEMENTATION ==========================================================
template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT >
void
ConstrainedSolver::
solve_const( const RMatrixT& _constraints,
const CMatrixT& _A,
VectorT& _x,
const VectorT& _rhs,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
ConstrainedSolver::solve_const(const RMatrixT& _constraints,
const CMatrixT& _A,
VectorT& _x,
const VectorT& _rhs,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
// copy matrices
RMatrixT constraints( gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
RMatrixT constraints(gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
gmm::copy(_constraints, constraints);
CMatrixT A( gmm::mat_nrows(_A), gmm::mat_ncols(_A));
CMatrixT A(gmm::mat_nrows(_A), gmm::mat_ncols(_A));
gmm::copy(_A, A);
// ... and vectors
......@@ -74,13 +72,13 @@ solve_const( const RMatrixT& _constraints,
// call non-const function
solve(constraints,
A,
_x,
rhs,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
A,
_x,
rhs,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
}
......@@ -89,20 +87,19 @@ solve_const( const RMatrixT& _constraints,
template<class RMatrixT, class VectorT, class VectorIT >
void
ConstrainedSolver::
solve_const( const RMatrixT& _constraints,
const RMatrixT& _B,
VectorT& _x,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
ConstrainedSolver::solve_const(const RMatrixT& _constraints,
const RMatrixT& _B,
VectorT& _x,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
// copy matrices
RMatrixT constraints( gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
RMatrixT constraints(gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
gmm::copy(_constraints, constraints);
RMatrixT B( gmm::mat_nrows(_B), gmm::mat_ncols(_B));
RMatrixT B(gmm::mat_nrows(_B), gmm::mat_ncols(_B));
gmm::copy(_B, B);
// ... and vectors
......@@ -110,12 +107,12 @@ solve_const( const RMatrixT& _constraints,
// call non-const function
solve(constraints,
B,
_x,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
B,
_x,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
}
......@@ -123,16 +120,15 @@ solve_const( const RMatrixT& _constraints,
template<class RMatrixT, class VectorT, class VectorIT >
void
ConstrainedSolver::
solve(
RMatrixT& _constraints,
RMatrixT& _B,
VectorT& _x,
VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
void
ConstrainedSolver::solve(
RMatrixT& _constraints,
RMatrixT& _B,
VectorT& _x,
VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
// convert into quadratic system
VectorT rhs;
......@@ -140,59 +136,59 @@ solve(
COMISO_GMM::factored_to_quadratic(_B, A, rhs);
// solve
solve( _constraints, A, _x, rhs,
_idx_to_round, _reg_factor,
_show_miso_settings,
_show_timings);
solve(_constraints, A, _x, rhs,
_idx_to_round, _reg_factor,
_show_miso_settings,
_show_timings);
// int nrows = gmm::mat_nrows(_B);
// int ncols = gmm::mat_ncols(_B);
// int ncons = gmm::mat_nrows(_constraints);
// int nrows = gmm::mat_nrows(_B);
// int ncols = gmm::mat_ncols(_B);
// int ncons = gmm::mat_nrows(_constraints);
// if( _show_timings) std::cerr << __FUNCTION__ << "\n Initial dimension: " << nrows << " x " << ncols << ", number of constraints: " << ncons << std::endl;
// // StopWatch for Timings
// Base::StopWatch sw, sw2; sw.start(); sw2.start();
// if( _show_timings) std::cerr << __FUNCTION__ << "\n Initial dimension: " << nrows << " x " << ncols << ", number of constraints: " << ncons << std::endl;
// // c_elim[i] = index of variable which is eliminated in condition i
// // or -1 if condition is invalid
// std::vector<int> c_elim( ncons);
// // StopWatch for Timings
// Base::StopWatch sw, sw2; sw.start(); sw2.start();
// // apply sparse gauss elimination to make subsequent _constraints independent
// make_constraints_independent( _constraints, _idx_to_round, c_elim);
// double time_gauss = sw.stop()/1000.0; sw.start();
// // c_elim[i] = index of variable which is eliminated in condition i
// // or -1 if condition is invalid
// std::vector<int> c_elim( ncons);
// // eliminate conditions and return column matrix Bcol
// gmm::col_matrix< gmm::rsvector< double > > Bcol( nrows, ncols);
// // apply sparse gauss elimination to make subsequent _constraints independent
// make_constraints_independent( _constraints, _idx_to_round, c_elim);
// double time_gauss = sw.stop()/1000.0; sw.start();
// // reindexing vector
// std::vector<int> new_idx;
// // eliminate conditions and return column matrix Bcol
// gmm::col_matrix< gmm::rsvector< double > > Bcol( nrows, ncols);
// eliminate_constraints( _constraints, _B, _idx_to_round, c_elim, new_idx, Bcol);
// double time_eliminate = sw.stop()/1000.0;
// // reindexing vector
// std::vector<int> new_idx;
// if( _show_timings) std::cerr << "Eliminated dimension: " << gmm::mat_nrows(Bcol) << " x " << gmm::mat_ncols(Bcol) << std::endl;
// eliminate_constraints( _constraints, _B, _idx_to_round, c_elim, new_idx, Bcol);
// double time_eliminate = sw.stop()/1000.0;
// // setup and solve system
// double time_setup = setup_and_solve_system( Bcol, _x, _idx_to_round, _reg_factor, _show_miso_settings);
// sw.start();
// if( _show_timings) std::cerr << "Eliminated dimension: " << gmm::mat_nrows(Bcol) << " x " << gmm::mat_ncols(Bcol) << std::endl;
// // double time_setup_solve = sw.stop()/1000.0; sw.start();
// // restore eliminated vars to fulfill the given conditions
// restore_eliminated_vars( _constraints, _x, c_elim, new_idx);
// // setup and solve system
// double time_setup = setup_and_solve_system( Bcol, _x, _idx_to_round, _reg_factor, _show_miso_settings);
// sw.start();
// double time_resubstitute = sw.stop()/1000.0; sw.start();
// // double time_setup_solve = sw.stop()/1000.0; sw.start();
// // double time_total = sw2.stop()/1000.0;
// // restore eliminated vars to fulfill the given conditions
// restore_eliminated_vars( _constraints, _x, c_elim, new_idx);
// if( _show_timings) std::cerr << "Timings: \n\t" <<
// "Gauss Elimination " << time_gauss << " s\n\t" <<
// "System Elimination " << time_eliminate << " s\n\t" <<
// "Setup " << time_setup << " s\n\t" <<
// // "Setup + Mi-Solver " << time_setup_solve << " s\n\t" <<
// "Resubstitution " << time_resubstitute << " s\n\t" << std::endl << std::endl;
// //"Total " << time_total << std::endl;
// double time_resubstitute = sw.stop()/1000.0; sw.start();
// // double time_total = sw2.stop()/1000.0;
// if( _show_timings) std::cerr << "Timings: \n\t" <<
// "Gauss Elimination " << time_gauss << " s\n\t" <<
// "System Elimination " << time_eliminate << " s\n\t" <<
// "Setup " << time_setup << " s\n\t" <<
// // "Setup + Mi-Solver " << time_setup_solve << " s\n\t" <<
// "Resubstitution " << time_resubstitute << " s\n\t" << std::endl << std::endl;
// //"Total " << time_total << std::endl;
}
......@@ -200,52 +196,51 @@ solve(
template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT>
void
ConstrainedSolver::
solve(
RMatrixT& _constraints,
CMatrixT& _A,
VectorT& _x,
VectorT& _rhs,
VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
void
ConstrainedSolver::solve(
RMatrixT& _constraints,
CMatrixT& _A,
VectorT& _x,
VectorT& _rhs,
VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
DEB_enter_func;
// show options dialog
if( _show_miso_settings)
if (_show_miso_settings)
miso_.show_options_dialog();
gmm::size_type nrows = gmm::mat_nrows(_A);
gmm::size_type ncols = gmm::mat_ncols(_A);
gmm::size_type ncons = gmm::mat_nrows(_constraints);
DEB_out_if( _show_timings, 1, "Initital dimension: " << nrows << " x " << ncols
<< ", number of constraints: " << ncons
<< " use reordering: " << use_constraint_reordering() << "\n")