Commit 4a4f973f authored by Max Lyon's avatar Max Lyon
Browse files

add exact constraint satisfaction example

parent c7dbd759
......@@ -574,6 +574,9 @@ if (COMISO_BUILD_EXAMPLES )
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_symmetric_dirichlet/CMakeLists.txt" )
add_subdirectory (Examples/small_symmetric_dirichlet)
endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_exact_constraint_satifaction_example/CMakeLists.txt" )
add_subdirectory (Examples/small_exact_constraint_satifaction_example)
endif()
endif (COMISO_BUILD_EXAMPLES )
......
include (CoMISoExample)
acg_add_executable (small_exact_constraint_satifaction_example ${sources} ${headers} )
# enable rpath linking
set_target_properties(small_exact_constraint_satifaction_example PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
target_link_libraries (small_exact_constraint_satifaction_example
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>
#include <CoMISo/NSolver/NewtonSolver.hh>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>
//------------------------------------------------------------------------------------------------------
class SmallNProblem : public COMISO::NProblemInterface
{
public:
// specify a function which has several local minima
// f(x,y) = x^4 + y^4
// number of unknown variables, here x and y = 2
virtual int n_unknowns ( )
{
return 2;
}
// initial value where the optimization should start from
virtual void initial_x ( double* _x )
{
_x[0] = 4.0;
_x[1] = 2.0;
}
// function evaluation at location _x
virtual double eval_f ( const double* _x )
{
return std::pow(_x[0], 4) + std::pow(_x[1], 4);
}
// gradient evaluation at location _x
virtual void eval_gradient( const double* _x, double* _g)
{
_g[0] = 4.0*std::pow(_x[0], 3);
_g[1] = 4.0*std::pow(_x[1], 3);
}
// hessian matrix evaluation at location _x
virtual void eval_hessian ( const double* _x, SMatrixNP& _H)
{
_H.resize(n_unknowns(), n_unknowns());
_H.setZero();
_H.coeffRef(0,0) = 12.0*std::pow(_x[0], 2);
_H.coeffRef(1,0) = 0.0;
_H.coeffRef(0,1) = 0.0;
_H.coeffRef(1,1) = 12.0*std::pow(_x[1], 2);
}
// print result
virtual void store_result ( const double* _x )
{
std::cerr << "Energy: " << eval_f(_x) << std::endl;
std::cerr << "(x,y) = (" << _x[0] << "," << _x[1] << ")" << std::endl;
solution.resize(n_unknowns());
for (int i = 0; i < n_unknowns(); ++i)
solution[i] = _x[i];
}
// advanced properties
virtual bool constant_hessian() const { return true; }
std::vector<double> solution;
};
// Example main
int main(void)
{
std::cout << "---------- 1) Get an instance of a problem..." << std::endl;
SmallNProblem problem;
std::cout << "---------- 2) Set up constraints..." << std::endl;
int n_constraints = 1; // there will be one constraints
Eigen::VectorXd b;
Eigen::SparseMatrix<double> A;
A.resize(n_constraints, problem.n_unknowns());
A.setZero();
b.resize(n_constraints);
b.setZero();
// first constraint: first variable equals two times second
A.coeffRef(0,0) = 1;
A.coeffRef(0,1) = -3;
b.coeffRef(0) = 0;
std::cout << A << std::endl << b << std::endl;
std::cout << "---------- 3) Solve with Newton Solver..." << std::endl;
COMISO::NewtonSolver nsolver;
nsolver.solve(&problem, A, b);
std::cout << "---------- 4) Print solution..." << std::endl;
std::cout << std::setprecision(100);
for (unsigned int i = 0; i < problem.n_unknowns(); ++i)
std::cout << "x[" << i << "] = " << problem.solution[i] << std::endl;
std::cout << "---------- 5) Check constraint violation..." << std::endl;
Eigen::VectorXd x;
x.resize(problem.n_unknowns());
for (unsigned int i = 0; i < problem.n_unknowns(); ++i)
x.coeffRef(i) = problem.solution[i];
std::cout << "Constraint violation: " << (A*x - b).squaredNorm() << std::endl;
std::cout << "---------- 6) Try to exactly fulfill constraints..." << std::endl;
x.coeffRef(0) = 3.0 * x.coeffRef(1);
std::cout << "---------- 7) Check constraint violation again..." << std::endl;
std::cout << "Constraint violation: " << (A*x - b).squaredNorm() << std::endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment