...
 
Commits (49)
......@@ -2,3 +2,6 @@ Config/config.hh
*.orig
*.rej
.project
*.swp
.DS_Store
/build*/
Subproject commit f31f6e3ae0a097c2e04bc750589e00ccc972feaf
Subproject commit 3623fbdfb4eba14a65926ff4034bc3916eab2cf0
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_;
};
......
......@@ -84,6 +84,10 @@ public:
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_;}
......
/*===========================================================================*\
* *
* 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 ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
class SymmetricDirichletElement
{
public:
// define dimensions
const static int NV = 6; // the three u and v coordinates. u1, v1, u2, v2, u3, v3
const static int NC = 6; // the three reference positions of the triangle, x1, y1, x2, y2, x3, y3
typedef Eigen::Matrix<size_t,NV,1> VecI;
typedef Eigen::Matrix<double,NV,1> VecV;
typedef Eigen::Matrix<double,NC,1> VecC;
typedef Eigen::Triplet<double> Triplet;
typedef Eigen::Matrix<double,12,1> Vector12;
typedef Eigen::Matrix<adouble,2,2> Matrix2x2ad;
typedef Eigen::Matrix<double,6,6> Matrix6x6;
SymmetricDirichletElement()
:
tape_(-1),
tape_available_(false)
{
init_tape();
}
double eval_f (const VecV& _x, const VecC& _c);
void eval_gradient(const VecV& _x, const VecC& _c, VecV& _g);
void eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets);
void project_hessian(Eigen::MatrixXd& H_orig, Eigen::MatrixXd& H_spd, double eps);
double max_feasible_step(const VecV& _x, const VecV& _v, const VecC& /*_c*/);
adouble f_adouble( const adouble* _x);
void retape();
void init_tape();
private:
// index of tape
short int tape_;
// taping already done?
bool tape_available_;
};
/** \class SymmetricDirichletProblem
A problem that allows you to add triangles with reference positions for which
the symmetric dirichlet energy should be minimized
*/
class COMISODLLEXPORT SymmetricDirichletProblem : public FiniteElementProblem
{
public:
typedef FiniteElementSet<SymmetricDirichletElement> SymmetricDirichletElementSet;
typedef Eigen::Matrix<size_t,3,1> IndexVector;
typedef Eigen::Matrix<double,3,2> ReferencePositionVector2D;
typedef Eigen::Matrix<double,3,3> ReferencePositionVector3D;
typedef Eigen::VectorXd VectorD;
typedef Eigen::SparseMatrix<double> SMatrixD;
/// Default constructor
SymmetricDirichletProblem(const unsigned int _n_vertices);
void add_triangle(const IndexVector& _vertex_indices, const ReferencePositionVector2D& _reference_positions);
/// add fix point constraint. Note that the final constraints still need to be passed
/// to the solver by you after retrieving them via get_constraints
void add_fix_point_constraint(int _vertex_index, double _fix_u, double _fix_v);
/// add fix point constraint for only one of the two coordinates of a vertex. Note that the
/// final constraints still need to be passed to the solver by you after retrieving them via get_constraints
void add_fix_coordinate_constraint(int _vertex_index, int _coordinate, double _fix_coordinate);
/// Retrieve the constraint matrix and right hand side representing the added
/// fix point and fix coordinate constraints for use e.g. with the NewtonSolver
void get_constraints(SMatrixD& _A, VectorD& _b);
/// remove all constraints
void clear_constraints() { fix_points.clear();}
private:
SymmetricDirichletElementSet element_set;
std::vector<std::pair<int, double>> fix_points;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif //(COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
//=============================================================================
#endif // COMISO_SYMMETRICDIRICHLETPROBLEM_HH defined
//=============================================================================
......@@ -223,8 +223,8 @@ ConstrainedSolver::solve(
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")
<< ", number of constraints: " << ncons << ", number of integer variables: " << _idx_to_round.size()
<< ", use reordering: " << (use_constraint_reordering() ? "yes" : "no") << "\n")
// StopWatch for Timings
Base::StopWatch sw, sw2; sw.start(); sw2.start();
......
This diff is collapsed.
include(CMakeFindDependencyMacro)
# Add the targets file
include("${CMAKE_CURRENT_LIST_DIR}/CoMISoTargets.cmake")
# - Try to find GUROBI
# GUROBI_BASE - The libraries needed to use Gurobi
# Once done this will define
# GUROBI_FOUND - System has Gurobi
# GUROBI_INCLUDE_DIRS - The Gurobi include directories
# GUROBI_LIBRARIES - The libraries needed to use Gurobi
set (GUROBI_ENABLE OFF CACHE BOOL "Enable gurobi?")
set(GUROBI_ENABLE OFF CACHE BOOL "Enable gurobi?")
if ( GUROBI_ENABLE )
set (GUROBI_BASE "c:" CACHE PATH "Base path of your gurobi installation")
if (GUROBI_ENABLE)
if (GUROBI_INCLUDE_DIR)
# in cache already
set(GUROBI_FOUND TRUE)
set(GUROBI_INCLUDE_DIRS "${GUROBI_INCLUDE_DIR}" )
set(GUROBI_LIBRARIES "${GUROBI_CXX_LIBRARY};${GUROBI_LIBRARY}" )
else (GUROBI_INCLUDE_DIR)
set(GUROBI_BASE $ENV{GUROBI_HOME} CACHE PATH "GUROBI root directory.")
find_path(GUROBI_INCLUDE_DIR
find_path(GUROBI_INCLUDE_DIR
NAMES gurobi_c++.h
PATHS "$ENV{GUROBI_HOME}/include"
"/Library/gurobi502/mac64/include"
"/Library/gurobi562/mac64/include"
"C:\\libs\\gurobi502\\include"
"C:\\libs\\gurobi562\\include"
"${GUROBI_BASE}/include"
PATHS
"${GUROBI_BASE}/include"
"$ENV{GUROBI_HOME}/include"
)
find_library( GUROBI_LIBRARY
NAMES gurobi
gurobi60
gurobi56
gurobi55
gurobi51
gurobi50
gurobi46
gurobi45
PATHS "$ENV{GUROBI_HOME}/lib"
"/Library/gurobi562/mac64/lib"
"/Library/gurobi502/mac64/lib"
"C:\\libs\\gurobi562\\lib"
"C:\\libs\\gurobi502\\lib"
"${GUROBI_BASE}/lib"
)
get_filename_component(GUROBI_LIB_DIR "${GUROBI_INCLUDE_DIR}/../lib" ABSOLUTE)
# GUROBI_BIN_DIR is needed on windows, where it contains the .dll
get_filename_component(GUROBI_BIN_DIR "${GUROBI_INCLUDE_DIR}/../bin" ABSOLUTE)
get_filename_component(GUROBI_SRC_DIR "${GUROBI_INCLUDE_DIR}/../src" ABSOLUTE)
file(GLOB GUROBI_LIBRARY_LIST
RELATIVE ${GUROBI_LIB_DIR}
${GUROBI_LIB_DIR}/libgurobi*.so
${GUROBI_LIB_DIR}/libgurobi*.dll)
# Ignore libgurobiXY_light.so, libgurobi.so (without version):
string(REGEX MATCHALL
"libgurobi([0-9]+)\\..*"
GUROBI_LIBRARY_LIST
"${GUROBI_LIBRARY_LIST}"
)
string(REGEX REPLACE
"libgurobi([0-9]+)\\..*"
"\\1"
GUROBI_LIBRARY_VERSIONS