...
 
Commits (63)
......@@ -2,3 +2,6 @@ Config/config.hh
*.orig
*.rej
.project
*.swp
.DS_Store
/build*/
......@@ -13,19 +13,6 @@ macos-c++11:
tags:
- Apple
CoMISo-VS2013-Qt-5.5.1-x64:
variables:
BUILD_PLATFORM: "VS2013"
ARCHITECTURE: "x64"
QT_VERSION: "Qt5.5.1"
GIT_SUBMODULE_STRATEGY: recursive
COMPILER: "VS2013"
script: "CI\\Windows.bat"
tags:
- VS2013
- Qt551
CoMISo-VS2017-Qt-5.10.1-x64:
variables:
BUILD_PLATFORM: "VS2017"
......
Subproject commit f31f6e3ae0a097c2e04bc750589e00ccc972feaf
Subproject commit 5e38f629fc7c7ae12316db1f19e6f878c1387bb2
This diff is collapsed.
......@@ -114,14 +114,13 @@ int main(void)
std::cout << "---------- 5) Solve with IPOPT solver... " << std::endl;
COMISO::IPOPTSolver ipopt;
ipopt.app().Options()->SetStringValue("derivative_test", "second-order");
ipopt.set_ipopt_option("derivative_test", "second-order");
ipopt.solve(&lsqp, constraints);
#endif
std::cout << "---------- 6) Print solution..." << std::endl;
for( int i=0; i<n; ++i)
std::cerr << "x_" << i << " = " << lsqp.x()[i] << std::endl;
return 0;
}
include (CoMISoExample)
acg_add_executable (small_symmetric_dirichlet ${sources} ${headers} )
# enable rpath linking
set_target_properties(small_symmetric_dirichlet PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
target_link_libraries (small_symmetric_dirichlet
CoMISo
${COMISO_LINK_LIBRARIES}
)
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2019 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 *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#include <CoMISo/Config/config.hh>
#include <iostream>
#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
#include <CoMISo/Utils/StopWatch.hh>
#include <CoMISo/NSolver/NewtonSolver.hh>
#include <CoMISo/NSolver/SymmetricDirichletProblem.hh>
#include <vector>
//------------------------------------------------------------------------------------------------------
// Example main
int main(void)
{
std::cout << "---------- 1) Get an instance of a SymmetricDirichletProblem..." << std::endl;
// then create finite element problem and add sets
unsigned int n_vertices = 4;
COMISO::SymmetricDirichletProblem sd_problem(n_vertices);
COMISO::SymmetricDirichletProblem::IndexVector indices(0,1,2);
COMISO::SymmetricDirichletProblem::ReferencePositionVector2D positions;
positions << 0, 0, // first point
1, 0, // second point
0, 1; // third point
sd_problem.add_triangle(indices, positions);
COMISO::SymmetricDirichletProblem::IndexVector indices2(3,2,1);
sd_problem.add_triangle(indices2, positions); // same reference positions can be used because we want both triangles to be isosceles
std::vector<double> initial_solution{0.0,0.0,
2.0,0.0,
2.0,2.0,
3.0,4.0};
sd_problem.x() = initial_solution;
std::cout << "---------- 2) Set up constraints..." << std::endl;
// fix first vertex to origin to fix translational degree of freedom
sd_problem.add_fix_point_constraint(0, 0.0, 0.0);
// fix v coordinate of second vertex to 0 to fix rotational degree of freedom
sd_problem.add_fix_coordinate_constraint(1, 1, 0.0);
std::cout << "---------- 3) Solve with Newton Solver..." << std::endl;
COMISO::SymmetricDirichletProblem::VectorD b;
COMISO::SymmetricDirichletProblem::SMatrixD A;
sd_problem.get_constraints(A, b);
COMISO::NewtonSolver nsolver;
nsolver.solve(&sd_problem, A, b);
// print result
for (unsigned int i = 0; i < n_vertices; ++i)
std::cerr << "p" << i << " = ( " << sd_problem.x()[2*i+0] << ", " << sd_problem.x()[2*i+1] << ")" << std::endl;
return 0;
}
#else // (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
int main(void)
{
std::cerr << "Warning: Example cannot be executed since either EIGEN3 or ADOLC is not available..." << std::endl;
return 0;
}
#endif // (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
......@@ -76,6 +76,8 @@ solve(NProblemInterface* _problem,
// move to next constraint
++n_constraints;
} else {
std::cerr << "Warning: COMISOSolver received a problem with non-equality constraints!!!" << std::endl;
}
// resize matrix to final number of constraints
......
......@@ -80,7 +80,7 @@ private:
} // namespace COMISO
//=============================================================================
// support std vectors
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint);
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint)
//=============================================================================
#endif // ACG_CONECONSTRAINT_HH defined
//=============================================================================
......
......@@ -182,6 +182,7 @@ public:
for(unsigned int j=0; j<NV; ++j)
x_[j] = _x[instances_.index(i,j)];
triplets_.clear();
element_.eval_hessian(x_, instances_.c(i), triplets_);
for(unsigned int j=0; j<triplets_.size(); ++j)
......
//=============================================================================
//
// STRUCT IPOPTCallbackParameters
//
//=============================================================================
#ifndef COMISO_IPOPTCALLBACKPARAMETERS_HH
#define COMISO_IPOPTCALLBACKPARAMETERS_HH
#include <IpTNLP.hpp>
//== TYPE DEFINITION CALLBACK PARAMETERS ======================================
namespace COMISO {
struct IPOPTCallbackParameters {
Ipopt::AlgorithmMode mode;
Ipopt::Index iter;
Ipopt::Number obj_value;
Ipopt::Number inf_pr;
Ipopt::Number inf_du;
Ipopt::Number mu;
Ipopt::Number d_norm;
Ipopt::Number regularization_size;
Ipopt::Number alpha_du;
Ipopt::Number alpha_pr;
Ipopt::Index ls_trials;
const Ipopt::IpoptData* ip_data;
Ipopt::IpoptCalculatedQuantities* ip_cq;
};
}
#endif
......@@ -15,10 +15,13 @@
//== INCLUDES =================================================================
#include <vector>
#include <cstddef>
#include <functional>
#include <CoMISo/Config/CoMISoDefines.hh>
#include "NProblemGmmInterface.hh"
#include "NProblemInterface.hh"
#include "IPOPTSolverLean.hh"
#include "NProblemGmmInterface.hh"
#include "NProblemInterface.hh"
#include "NConstraintInterface.hh"
......@@ -30,25 +33,22 @@
#include <CoMISo/Utils/gmm.hh>
#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
#include <vector>
#include <cstddef>
//== NAMESPACES ===============================================================
namespace COMISO {
//== FORWARDDECLARATIONS ======================================================
class NProblemGmmInterface; // deprecated
class NProblemInterface;
class NConstraintInterface;
struct IPOPTCallbackParameters;
//== CLASS DEFINITION PROBLEM INSTANCE=========================================================
class NProblemIPOPT : public Ipopt::TNLP
//== CLASS DEFINITION PROBLEM INSTANCE ========================================
class IPOPTProblemInstance : public Ipopt::TNLP
{
public:
......@@ -65,8 +65,11 @@ public:
typedef NProblemInterface::SMatrixNP SMatrixNP;
/** default constructor */
NProblemIPOPT(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints, const bool _hessian_approximation = false)
: problem_(_problem), store_solution_(false), hessian_approximation_(_hessian_approximation)
IPOPTProblemInstance(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);
......@@ -127,6 +130,9 @@ public:
IpoptCalculatedQuantities* ip_cq) override;
//@}
/** Set intermediate callback function object **/
void set_callback_function(std::function<bool(const IPOPTCallbackParameters &)>);
/** Intermediate Callback method for the user. Providing dummy
* default implementation. For details see IntermediateCallBack
* in IpNLP.hpp. */
......@@ -163,8 +169,8 @@ private:
*/
//@{
// MyNLP();
NProblemIPOPT(const NProblemIPOPT&);
NProblemIPOPT& operator=(const NProblemIPOPT&);
IPOPTProblemInstance(const IPOPTProblemInstance&);
IPOPTProblemInstance& operator=(const IPOPTProblemInstance&);
//@}
// split user-provided constraints into general-constraints and bound-constraints
......@@ -193,13 +199,15 @@ private:
std::vector<double> x_;
bool hessian_approximation_;
std::function<bool(const IPOPTCallbackParameters &)> intermediate_callback_;
};
//== CLASS DEFINITION PROBLEM INSTANCE=========================================================
class NProblemGmmIPOPT : public Ipopt::TNLP
class IPOPTProblemInstanceGmm : public Ipopt::TNLP
{
public:
......@@ -226,7 +234,7 @@ public:
typedef gmm::linalg_traits<SVectorNP>::iterator SVectorNP_iter;
/** default constructor */
NProblemGmmIPOPT(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints)
IPOPTProblemInstanceGmm(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints)
: problem_(_problem), constraints_(_constraints), nnz_jac_g_(0), nnz_h_lag_(0)
{}
......@@ -285,6 +293,9 @@ public:
IpoptCalculatedQuantities* ip_cq) override;
//@}
/** Set intermediate callback function object **/
void set_callback_function(std::function<bool(const IPOPTCallbackParameters &)>);
/** Intermediate Callback method for the user. Providing dummy
* default implementation. For details see IntermediateCallBack
* in IpNLP.hpp. */
......@@ -314,8 +325,8 @@ private:
*/
//@{
// MyNLP();
NProblemGmmIPOPT(const NProblemGmmIPOPT&);
NProblemGmmIPOPT& operator=(const NProblemGmmIPOPT&);
IPOPTProblemInstanceGmm(const IPOPTProblemInstanceGmm&);
IPOPTProblemInstanceGmm& operator=(const IPOPTProblemInstanceGmm&);
//@}
......@@ -337,6 +348,8 @@ private:
// Sparse Matrix of problem (don't initialize every time!!!)
SMatrixNP HP_;
std::function<bool(const IPOPTCallbackParameters &)> intermediate_callback_;
};
//=============================================================================
......@@ -347,4 +360,3 @@ private:
//=============================================================================
#endif // COMISO_NPROBLEMIPOPT_HH
//=============================================================================
This diff is collapsed.
......@@ -8,134 +8,144 @@
#ifndef COMISO_IPOPTSOLVER_HH
#define COMISO_IPOPTSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_IPOPT_AVAILABLE
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <CoMISo/Utils/gmm.hh>
#include <vector>
#include <cstddef>
#include "NProblemGmmInterface.hh"
#include "NProblemInterface.hh"
#include "NProblemIPOPT.hh"
#include "NConstraintInterface.hh"
#include "BoundConstraint.hh"
#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
#include <functional>
#include <vector>
#include <string>
//== FORWARDDECLARATIONS ======================================================
#include <CoMISo/Config/CoMISoDefines.hh>
//== NAMESPACES ===============================================================
namespace COMISO {
//== FORWARDDECLARATIONS ======================================================
class NProblemInterface;
class NConstraintInterface;
struct IPOPTCallbackParameters;
//== CLASS DEFINITION =========================================================
/** \class IPOPTSolver
Solver for Interior Point optimization problems.
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
Solves an interior point problem, given an NProblemInterface
instance and optionally a set of constraints as well as "lazy
constraints" via NConstraintInterface.
Brief Description.
A more elaborate description follows.
Lazy constraints are not active while the initial solution to the
problem is computed. After the first solution is found, the lazy
constraints are checked and added to the set of active constraints
if they are violated. This process is then repeated until all
constraints are satisfied OR a maximum number of solution attempts
has been reached. In that case the optimization is started once
more, with all lazy constraints active.
*/
class COMISODLLEXPORT IPOPTSolver
{
public:
/// Default constructor -> set up IpOptApplication
IPOPTSolver();
// ********** SOLVE **************** //
// solve -> returns ipopt status code
//------------------------------------------------------
// enum ApplicationReturnStatus
// {
// Solve_Succeeded=0,
// Solved_To_Acceptable_Level=1,
// Infeasible_Problem_Detected=2,
// Search_Direction_Becomes_Too_Small=3,
// Diverging_Iterates=4,
// User_Requested_Stop=5,
// Feasible_Point_Found=6,
//
// Maximum_Iterations_Exceeded=-1,
// Restoration_Failed=-2,
// Error_In_Step_Computation=-3,
// Maximum_CpuTime_Exceeded=-4,
// Not_Enough_Degrees_Of_Freedom=-10,
// Invalid_Problem_Definition=-11,
// Invalid_Option=-12,
// Invalid_Number_Detected=-13,
//
// Unrecoverable_Exception=-100,
// NonIpopt_Exception_Thrown=-101,
// Insufficient_Memory=-102,
// Internal_Error=-199
// };
//------------------------------------------------------
int solve(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints);
// same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
int solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
const double _almost_infeasible = 0.5,
const int _max_passes = 5 );
// for convenience, if no constraints are given
int solve(NProblemInterface* _problem);
// deprecated interface for backwards compatibility
int solve(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints);
// ********* CONFIGURATION ********************* //
// access the ipopt-application (for setting parameters etc.)
// examples: app().Options()->SetIntegerValue("max_iter", 100);
// app().Options()->SetStringValue("derivative_test", "second-order");
// app().Options()->SetStringValue("hessian_approximation", "limited-memory");
IPOPTSolver();
~IPOPTSolver();
// *********** OPTIONS **************//
/*!
Set options of the underlying ipopt solver.
For a thorough list and documentation of available options, refer
to: https://www.coin-or.org/Ipopt/documentation/node40.html
*/
void set_ipopt_option(std::string, const int&);
void set_ipopt_option(std::string, const double&);
void set_ipopt_option(std::string, const std::string&);
/*!
Get options of the underlying ipopt solver.
The type of option {int, double, std::string} needs to be passed as
template argument.
*/
template<typename T>
T get_ipopt_option(std::string option);
/*!
Set the maximum number of iterations
*/
void set_max_iterations(const int _max_iterations);
int get_max_iterations() const;
/*! Set the threshold on the lazy inequality constraint to decide
if we are near the constraint boundary.
*/
void set_almost_infeasible_threshold(const double _alm_infsb_thrsh);
double get_almost_infeasible_threshold() const;
/*!
Set the max number of incremental lazy constraint iterations before switching
to the fully constrained problem.
\note The default value is 5.
*/
void set_incremental_lazy_constraint_max_iteration_number
(const int _incr_lazy_cnstr_max_iter_nmbr);
int get_incremental_lazy_constraint_max_iteration_number() const;
/*
Turn on/off solving the fully constraint problem after exhausting the
incremental lazy constraint iterations.
\note The default value of this is true.
*/
void set_enable_all_lazy_contraints(const bool _enbl_all_lzy_cnstr);
bool get_enable_all_lazy_contraints() const;
/*
Set intermediate callback function object. For the definition of
IPOPTCallbackParameters include the IPOPTCallbackParameters.hh
header.
If the callback function returns false, IPOPT will terminate
prematurely with the User_Requested_Stop status.
*/
void set_callback_function
(std::function<bool(const IPOPTCallbackParameters &)>);
Ipopt::IpoptApplication& app() {return (*app_); }
// ********** SOLVE **************** //
//! Solve a problem instance with an optional set of constraints.
//! \throws Outcome
void solve
(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints = {});
void set_print_level(const int _level)
{
app().Options()->SetIntegerValue("print_level", _level);
print_level_ = _level;
}
//! Same as above with additional lazy constraints that are only
//! added iteratively to the problem if not satisfied.
//! \throws Outcome
void solve
(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints);
protected:
double* P(std::vector<double>& _v)
{
if( !_v.empty())
return ((double*)&_v[0]);
else
return 0;
}
//! Get the computed solution energy
double energy();
private:
class Impl;
Impl* impl_;
// smart pointer to IpoptApplication to set options etc.
Ipopt::SmartPtr<Ipopt::IpoptApplication> app_;
int print_level_;
// inhibit copy
IPOPTSolver(const IPOPTSolver&);
IPOPTSolver& operator=(const IPOPTSolver&);
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_IPOPT_AVAILABLE
//=============================================================================
#endif // ACG_IPOPTSOLVER_HH defined
#endif // COMISO_IPOPTSOLVER_HH defined
//=============================================================================
This diff is collapsed.
//=============================================================================
//
// CLASS IPOPTSolverLean
//
//=============================================================================
#ifndef COMISO_IPOPTLEANSOLVER_HH
#define COMISO_IPOPTLEANSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_IPOPT_AVAILABLE
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <vector>
#include <cstddef>
//== NAMESPACES ===============================================================
namespace COMISO {
//== FORWARDDECLARATIONS ======================================================
class NProblemGmmInterface; // deprecated
class NProblemInterface;
class NConstraintInterface;
//== CLASS DEFINITION =========================================================
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
Brief Description.
A more elaborate description follows.
*/
class COMISODLLEXPORT IPOPTSolverLean
{
public:
/// Default constructor
IPOPTSolverLean();
/// Destructor
~IPOPTSolverLean();
// *********** OPTIONS **************//
/*!
Set the maximum number of iterations
*/
void set_max_iterations(const int _max_iterations);
int max_iterations() const;
/*!
Set the threshold on the lazy inequality constraint to decide if we are near
the constraint boundary.
*/
void set_almost_infeasible_threshold(const double _alm_infsb_thrsh);
double almost_infeasible_threshold() const;
/*!
Set the max number of incremental lazy constraint iterations before switching
to the fully constrained problem.
\note The default value is 5.
*/
void set_incremental_lazy_constraint_max_iteration_number(const int
_incr_lazy_cnstr_max_iter_nmbr);
int incremental_lazy_constraint_max_iteration_number() const;
/*
Turn on/off solving the fully constraint problem after exhausting the
incremental lazy constraint iterations.
\note The default value of this is true.
*/
void set_enable_all_lazy_contraints(const bool _enbl_all_lzy_cnstr);
bool enable_all_lazy_contraints() const;
// ********** SOLVE **************** //
//! \throws Outcome
void solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints);
//! Same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
//! \throws Outcome
void solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints);
// for convenience, if no constraints are given
//! \throws Outcome
void solve(NProblemInterface* _problem);
// deprecated interface for backwards compatibility
//! \throws Outcome
void solve(NProblemGmmInterface* _problem,
std::vector<NConstraintInterface*>& _constraints);
//! Get the computed solution energy
double energy();
private:
class Impl;
Impl* impl_;
// inhibit copy
IPOPTSolverLean(const IPOPTSolverLean&);
IPOPTSolverLean& operator=(const IPOPTSolverLean&);
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_IPOPT_AVAILABLE
//=============================================================================
#endif // ACG_IPOPTSOLVER_HH defined
//=============================================================================
......@@ -88,7 +88,7 @@ private:
} // namespace COMISO
//=============================================================================
// support std vectors
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::LinearConstraint);
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::LinearConstraint)
//=============================================================================
#endif // ACG_LINEARCONSTRAINT_HH defined
//=============================================================================
......
......@@ -23,6 +23,7 @@ solve(NProblemGmmInterface* _problem)
{
DEB_enter_func;
#if COMISO_SUITESPARSE_AVAILABLE
converged_ = true;
// get problem size
int n = _problem->n_unknowns();
......@@ -98,16 +99,19 @@ solve(NProblemGmmInterface* _problem)
_problem->store_result(P(x));
DEB_line(2, "Newton solver reached max regularization but did not "
"converge");
converged_ = false;
return false;
}
}
}
_problem->store_result(P(x));
DEB_line(2, "Newton Solver did not converge!!! after iterations.");
converged_ = false;
return false;
#else
DEB_warning(1,"NewtonSolver requires not-available CholmodSolver");
converged_ = false;
return false;
#endif
}
......@@ -120,6 +124,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
const VectorD& _b)
{
DEB_time_func_def;
converged_ = false;
const double KKT_res_eps = 1e-6;
const int max_KKT_regularization_iters = 40;
......@@ -160,10 +165,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
double kkt_res2(0.0);
double constraint_res2(0.0);
int reg_iters(0);
bool fact_ok = true;
do
{
// get Newton search direction by solving LSE
bool fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
first_factorization = false;
if(fact_ok)
......@@ -204,7 +210,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
}
++reg_iters;
}
while( (kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
while( (!fact_ok || kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
// no valid step could be found?
if(kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2 || reg_iters >= max_KKT_regularization_iters)
......@@ -250,8 +256,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
<< ", KKT residual^2 = " << kkt_res2);
// converged?
if(newton_decrement < eps_ || std::abs(t) < eps_ls_)
if(newton_decrement < eps_ || std::abs(t) < eps_ls_ )
{
converged_ = true;
break;
}
++iter;
}
......@@ -397,11 +406,11 @@ bool NewtonSolver::numerical_factorization(SMatrixD& _KKT)
DEB_enter_func;
switch(solver_type_)
{
case LS_EigenLU:
case LS_EigenLU:
lu_solver_.factorize(_KKT);
return (lu_solver_.info() == Eigen::Success);
#if COMISO_SUITESPARSE_AVAILABLE
case LS_Umfpack:
case LS_Umfpack:
umfpack_solver_.factorize(_KKT);
return (umfpack_solver_.info() == Eigen::Success);
#endif
......
......@@ -84,6 +84,8 @@ public:
solver_type_ = _st;
}
bool converged() { return converged_; }
protected:
bool factorize(NProblemInterface* _problem, const SMatrixD& _A,
......@@ -149,6 +151,8 @@ private:
// deprecated
bool constant_hessian_structure_;
bool converged_;
};
......
......@@ -3,15 +3,18 @@
#include <CoMISo/Config/config.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <Base/Code/Quality.hh>
#include <Base/Debug/DebOut.hh>
LOW_CODE_QUALITY_SECTION_BEGIN
#include <Eigen/Eigen>
#include <Eigen/Core>
#include <Eigen/Sparse>
LOW_CODE_QUALITY_SECTION_END
#include <vector>
//== NAMESPACES ===============================================================
......@@ -26,36 +29,31 @@ public:
// typedef Eigen::DynamicSparseMatrix<double,Eigen::ColMajor> SMatrixNP;
QuadraticProblem()
: A_(0,0), b_(Eigen::VectorXd::Index(0)), c_(0.0)
{
}
: A_(0, 0), b_(Eigen::VectorXd::Index(0)), c_(0.0)
{}
QuadraticProblem(SMatrixNP& _A, Eigen::VectorXd& _b, const double _c)
: A_(_A), c_(_c)
QuadraticProblem(SMatrixNP& _A, const Eigen::VectorXd& _b, const double _c)
{
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
b_ = _b;
x_ = Eigen::VectorXd::Zero(A_.cols());
set_A(_A);
set_b(_b);
set_c(_c);
}
// number of unknowns
virtual int n_unknowns()
{
return A_.rows();
return A_.rows();
}
// initial value where the optimization should start from
virtual void initial_x(double* _x)
{
for( int i=0; i<this->n_unknowns(); ++i)
_x[i] = x_[i];
for (int i = 0; i < this->n_unknowns(); ++i)
_x[i] = x_[i];
}
// function evaluation at location _x
virtual double eval_f( const double* _x )
virtual double eval_f(const double* _x)
{
Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
......@@ -63,38 +61,41 @@ public:
}
// gradient evaluation at location _x
virtual void eval_gradient( const double* _x, double* _g)
virtual void eval_gradient(const double* _x, double* _g)
{
Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
Eigen::Map<Eigen::VectorXd> g(_g, this->n_unknowns());
g = A_*x - b_;
}
}
// hessian matrix evaluation at location _x
virtual void eval_hessian ( const double* _x, SMatrixNP& _H)
virtual void eval_hessian(const double* _x, SMatrixNP& _H)
{
_H = A_;
}
// print result
virtual void store_result ( const double* _x )
virtual void store_result(const double* _x)
{
Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
x_ = x;
}
const SMatrixNP &A() const {return A_;}
const Eigen::VectorXd &b() const {return b_;}
double c() const {return c_;}
// get current solution
Eigen::VectorXd& x() { return x_;}
Eigen::VectorXd& x() { return x_; }
// advanced properties
virtual bool constant_hessian() const { return true; }
virtual bool constant_hessian() const { return true; }
void set_A(const SMatrixNP& _A)
{
DEB_error_if(_A.rows() != _A.cols(), "Square matrix expected");
A_ = _A;
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
x_ = Eigen::VectorXd::Zero(A_.cols());
}
......@@ -103,19 +104,18 @@ public:
b_ = _b;
}
void set_c( const double _c)
void set_c(const double _c)
{
c_ = _c;
}
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_;
SMatrixNP A_;
Eigen::VectorXd b_;
double c_;
// current solution, which is also used as initial value
Eigen::VectorXd x_;
};
//=============================================================================
......
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2019 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 *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#include <CoMISo/Config/config.hh>
#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
#include "SymmetricDirichletProblem.hh"
namespace COMISO {
double SymmetricDirichletElement::eval_f(const VecV& _x, const VecC& _c)
{
double y = 0.0;
Vector12 x;
x << _x[0], _x[1], _x[2], _x[3], _x[4], _x[5],
_c[0], _c[1], _c[2], _c[3], _c[4], _c[5];
if(tape_available_)
{
int ec = function(tape_, 1, 12, const_cast<double*>(x.data()), &y);
#ifdef ADOLC_RET_CODES
std::cout << "Info: function() returned code " << ec << std::endl;
#endif
// tape not valid anymore? retape and evaluate again
if(ec < 0)
{
retape();
function(tape_, 1, 12, const_cast<double*>(x.data()), &y);
}
}
else
{
retape();
function(tape_, 1, 12, const_cast<double*>(x.data()), &y);
}
return y;
}
void SymmetricDirichletElement::eval_gradient(const VecV& _x, const VecC& _c, VecV& _g)
{
Vector12 x;
x << _x[0], _x[1], _x[2], _x[3], _x[4], _x[5],
_c[0], _c[1], _c[2], _c[3], _c[4], _c[5];
Vector12 grad;
grad.setZero();
// evaluate gradient
int ec = gradient(tape_, 12, x.data(), grad.data());
// check if retaping is required
if(ec < 0)
{
#ifdef ADOLC_RET_CODES
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
#endif
retape();
gradient(tape_, 12, x.data(), grad.data());
}
_g = grad.block(0,0,6,1);
}
void SymmetricDirichletElement::eval_hessian(const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets)
{
Vector12 x;
x << _x[0], _x[1], _x[2], _x[3], _x[4], _x[5],
_c[0], _c[1], _c[2], _c[3], _c[4], _c[5];
// dense hessian
{
auto dense_hessian = new double*[12];
for(int i = 0; i < 12; ++i)
dense_hessian[i] = new double[i+1];
int ec = hessian(tape_, 12, const_cast<double*>(x.data()), dense_hessian);
if(ec < 0) {
#ifdef ADOLC_RET_CODES
std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
#endif
// Retape function if return code indicates discontinuity
retape();
ec = hessian(tape_, 12, const_cast<double*>(x.data()), dense_hessian);
}
#ifdef ADOLC_RET_CODES
std::cout << "Info: hessian() returned code " << ec << std::endl;
#endif
Eigen::MatrixXd H(6,6);
for (int i = 0; i < 6; ++i)
{
H(i,i) = dense_hessian[i][i];
for (int j = 0; j < i; ++j)
{
H(i,j) = dense_hessian[i][j];
H(j,i) = dense_hessian[i][j];
}
}
Eigen::MatrixXd Hspd(6,6);
project_hessian(H, Hspd, 1e-6);
// Hspd = H;
_triplets.reserve(6*6);
for (int i = 0; i < 6; ++i)
for (int j = 0; j < 6; ++j)
_triplets.push_back(Triplet(i,j,Hspd(i,j)));
for (int i = 0; i < 6; ++i)
delete dense_hessian[i];
delete[] dense_hessian;
}
}
void SymmetricDirichletElement::project_hessian(Eigen::MatrixXd& H_orig, Eigen::MatrixXd& H_spd, double eps)
{
// Compute eigen-decomposition (of symmetric matrix)
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(H_orig);
Eigen::MatrixXd V = eig.eigenvectors();
Eigen::MatrixXd D = eig.eigenvalues().asDiagonal();
// Clamp all eigenvalues to eps
for (int i = 0; i < H_orig.rows(); ++i)
D(i, i) = std::max(eps, D(i, i));
H_spd = V * D * V.inverse();
}
double SymmetricDirichletElement::max_feasible_step(const VecV& _x, const VecV& _v, const VecC&)
{
// get quadratic coefficients (ax^2 + b^x + c)
auto U11 = _x[0];
auto U12 = _x[1];
auto U21 = _x[2];
auto U22 = _x[3];
auto U31 = _x[4];
auto U32 = _x[5];
auto V11 = _v[0];
auto V12 = _v[1];
auto V21 = _v[2];
auto V22 = _v[3];
auto V31 = _v[4];
auto V32 = _v[5];
double a = V11*V22 - V12*V21 - V11*V32 + V12*V31 + V21*V32 - V22*V31;
double b = U11*V22 - U12*V21 - U21*V12 + U22*V11 - U11*V32 + U12*V31 + U31*V12 - U32*V11 + U21*V32 - U22*V31 - U31*V22 + U32*V21;
double c = U11*U22 - U12*U21 - U11*U32 + U12*U31 + U21*U32 - U22*U31;
double delta_in = pow(b,2) - 4*a*c;
if (delta_in < 0) {
return std::numeric_limits<double>::max();
}
double delta = sqrt(delta_in);
double t1 = (-b + delta)/ (2*a);
double t2 = (-b - delta)/ (2*a);
double tmp_n = std::min(t1,t2);
t1 = std::max(t1,t2); t2 = tmp_n;
// return the smallest negative root if it exists, otherwise return infinity
if (t1 > 0)
{
if (t2 > 0)
{
return 0.999 * t2;
}
else
{
return 0.999 * t1;
}
}
else
{
return std::numeric_limits<double>::max();
}
}
adouble SymmetricDirichletElement::f_adouble(const adouble* _x)
{
Matrix2x2ad B;
B(0,0) = _x[2]-_x[0];
B(0,1) = _x[4]-_x[0];
B(1,0) = _x[3]-_x[1];
B(1,1) = _x[5]-_x[1];
Matrix2x2ad Bin = B.inverse();
Matrix2x2ad R;
R(0,0) = _x[6+2]-_x[6+0];
R(0,1) = _x[6+4]-_x[6+0];
R(1,0) = _x[6+3]-_x[6+1];
R(1,1) = _x[6+5]-_x[6+1];
Matrix2x2ad Rin = R.inverse();
adouble area = 0.5 * R.determinant();
if (B.determinant() * area <= 0)
{
adouble res = std::numeric_limits<double>::max();
return res;
}
Matrix2x2ad J = B * Rin;
Matrix2x2ad Jin = R * Bin;
adouble res = J.squaredNorm() + Jin.squaredNorm();
return area * (res - 4);
}
void SymmetricDirichletElement::retape()
{
std::cerr << "re-tape..." << std::endl;
adouble y_d = 0.0;
double y = 0.0;
trace_on(tape_); // Start taping
adouble* xa = new adouble[12];
// Fill data vectors
xa[0] <<= 0.0;
xa[1] <<= 0.0;
xa[2] <<= 1.0;
xa[3] <<= 0.0;
xa[4] <<= 0.0;
xa[5] <<= 1.0;
xa[6+0] <<= 0.0;
xa[6+1] <<= 0.0;
xa[6+2] <<= 1.0;
xa[6+3] <<= 0.0;
xa[6+4] <<= 0.0;
xa[6+5] <<= 1.0;
y_d = f_adouble(xa);
y_d >>= y;
trace_off();
#ifdef ADOLC_STATS
print_stats();
#endif
delete[] xa;
}
void SymmetricDirichletElement::init_tape()
{
static bool tape_initialized = false;
static short int tape;
if (!tape_initialized)
{
tape = static_cast<short int>(TapeIDSingleton::Instance()->requestId());
tape_ = tape;
retape();
}
tape_initialized = true;
tape_available_ = true;
tape_ = tape;
}
/// Default constructor
SymmetricDirichletProblem::SymmetricDirichletProblem(const unsigned int _n_vertices)
:
FiniteElementProblem(2*_n_vertices)
{
FiniteElementProblem::add_set(&element_set);
}
void SymmetricDirichletProblem::add_triangle(const IndexVector& _vertex_indices, const ReferencePositionVector2D& _reference_positions)
{
SymmetricDirichletElement::VecI indices;
indices << 2*_vertex_indices[0], 2*_vertex_indices[0]+1,
2*_vertex_indices[1], 2*_vertex_indices[1]+1,
2*_vertex_indices[2], 2*_vertex_indices[2]+1;
SymmetricDirichletElement::VecC constants;
constants << _reference_positions(0, 0), _reference_positions(0, 1),
_reference_positions(1, 0), _reference_positions(1, 1),
_reference_positions(2, 0), _reference_positions(2, 1);
element_set.instances().add_element(indices, constants);
}
void SymmetricDirichletProblem::add_fix_point_constraint(int _vertex_index, double _fix_u, double _fix_v)
{
add_fix_coordinate_constraint(_vertex_index, 0, _fix_u);
add_fix_coordinate_constraint(_vertex_index, 1, _fix_v);
}
void SymmetricDirichletProblem::add_fix_coordinate_constraint(int _vertex_index, int _coordinate, double _fix_coordinate)
{
fix_points.push_back(std::make_pair(2*_vertex_index + _coordinate, _fix_coordinate));
}
void SymmetricDirichletProblem::get_constraints(SMatrixD& _A, VectorD& _b)
{
_A.resize(fix_points.size(), n_unknowns());
_b.resize(fix_points.size());
std::vector<Eigen::Triplet<double>> triplets;
for (size_t i = 0; i < fix_points.size(); ++i)
{
_b[i] = fix_points[i].second;
triplets.push_back(Eigen::Triplet<double>(i, fix_points[i].first, 1));
}
_A.setFromTriplets(triplets.begin(), triplets.end());
}
} // namespace COMISO
#endif //(COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2019 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 *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
//=============================================================================
//
// CLASS SymmetricDirichletProblem
//
//=============================================================================
#ifndef COMISO_SYMMETRICDIRICHLETTPROBLEM_HH
#define COMISO_SYMMETRICDIRICHLETTPROBLEM_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/config.hh>
#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
#include <CoMISo/Config/CoMISoDefines.hh>
#include "FiniteElementProblem.hh"
#include <adolc/adolc.h>
#include <adolc/adouble.h>
#include <adolc/drivers/drivers.h>
#include <adolc/sparse/sparsedrivers.h>
#include <adolc/taping.h>
#include <adolc/drivers/taylor.h>
#include <CoMISo/NSolver/TapeIDSingleton.hh>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================