Commit 9f152982 authored by Robin Brost's avatar Robin Brost
Browse files

Implementing the transformation of the Matrix in a new class and testing it

parent 5fd3c20c
Pipeline #12461 failed with stages
in 10 minutes and 11 seconds
This diff is collapsed.
......@@ -29,6 +29,7 @@
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>
#include <CoMISo/Utils/ExactConstraintSatisfaction.hh>
//------------------------------------------------------------------------------------------------------
class SmallNProblem : public COMISO::NProblemInterface
......@@ -89,7 +90,6 @@ public:
std::vector<double> solution;
};
// Example main
int main(void)
{
......@@ -97,7 +97,7 @@ int main(void)
SmallNProblem problem;
std::cout << "---------- 2) Set up constraints..." << std::endl;
int n_constraints = 1; // there will be one constraints
int n_constraints = 3; // there will be one constraints
Eigen::VectorXd b;
Eigen::SparseMatrix<double> A;
A.resize(n_constraints, problem.n_unknowns());
......@@ -106,9 +106,15 @@ int main(void)
b.setZero();
// first constraint: first variable equals three times second
A.coeffRef(0,0) = 1;
A.coeffRef(0,1) = -3;
A.coeffRef(0,0) = 2;
A.coeffRef(0,1) = -6;
A.coeffRef(1,0) = 10;
A.coeffRef(1,1) = -3;
A.coeffRef(2,0) = 1;
A.coeffRef(2,1) = -3;
b.coeffRef(0) = 0;
b.coeffRef(1) = 27;
b.coeffRef(2) = 0;
std::cout << A << std::endl << b << std::endl;
......@@ -134,6 +140,16 @@ int main(void)
std::cout << "---------- 6) Try to exactly fulfill constraints..." << std::endl;
Eigen::SparseMatrix<int> C = A.cast<int>();
Eigen::VectorXi d = b.cast<int>();
ExactConstraintSatisfaction satisfy;
satisfy.printMatrix(C);
satisfy.IREF_Gaussian(&C,&d);
satisfy.printMatrix(C);
satisfy.IRREF_Jordan(&C,&d);
satisfy.printMatrix(C);
satisfy.printVector(d);
x.coeffRef(0) = 3.0 * x.coeffRef(1);
std::cout << "---------- 7) Check constraint violation again..." << std::endl;
......
Logo-1.png

26.5 KB

#include "ExactConstraintSatisfaction.hh"
#include <CoMISo/Config/config.hh>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>
ExactConstraintSatisfaction::ExactConstraintSatisfaction()
{
}
// ------------------------Helpfull Methods----------------------------------------
void ExactConstraintSatisfaction::swapRows(Eigen::SparseMatrix<int>* A, int row1, int row2){
int temp = 0;
for(int i = 0; i < A->cols(); i++){
temp = A->coeffRef(row1,i);
A->coeffRef(row1,i) = A->coeffRef(row2, i);
A->coeffRef(row2, i) = temp;
}
}
int ExactConstraintSatisfaction::gcd(const int a, const int b){
if(b == 0 && a == 0) //to prevent dividing by 0 later
return 1;
if(b == 0){
return a;
}
return gcd(b , a % b);
}
int ExactConstraintSatisfaction::gcdRow(const Eigen::SparseMatrix<int>::RowXpr row, const int b){
int gcdValue = 0;
for(int i = 0; i < row.cols(); i ++){
gcdValue = gcd(gcdValue, row.coeff(i));
}
gcdValue = gcd(gcdValue, b);
return gcdValue;
}
void ExactConstraintSatisfaction::printMatrix(Eigen::SparseMatrix<int> A){
for(int i = 0; i < A.rows(); i++){
for(int j = 0; j < A.cols(); j++){
std::cout << A.coeffRef(i,j) << " ";
}
std::cout << std::endl;
}
}
void ExactConstraintSatisfaction::printVector(Eigen::VectorXi b)
{
std::cout << "the vector contains the elements: ";
for(int i = 0; i < b.size(); i++){
std::cout << b.coeffRef(i) << " ";
}
std::cout << std::endl;
}
// ----------------------Matrix transformation-----------------------------------
void ExactConstraintSatisfaction::IREF_Gaussian(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b){
number_pivots = 0;
int rows = A->rows(); //number of rows
int cols = A->cols(); //number of columns
for (int k = 0; k < rows;k++) { //order Matrix after pivot
//needed for the first line
int gcdValue = gcdRow(A->row(k), b->coeffRef(k)); //compute the gcd to make the values as small as possible
for(int x = 0; x < cols; x++){
A->coeffRef(k, x) = A->coeffRef(k,x) / gcdValue;
}
b->coeffRef(k) = b->coeffRef(k) / gcdValue;
if(k < cols){
if(A->coeffRef(k,k) == 0){
int l = -1;
for (int i = k +1;i < rows;i++) { //find row with pivot in this column
if(A->coeffRef(i,k) != 0){
l = i;
number_pivots++;
break;
}
}
if(l == -1)
continue;
swapRows(A,l,k); //swap rows so the pivot is in the right row
}else{
number_pivots++;
}
}
for(int i = k+1; i < rows; i++){
for(int j = k; j < cols; j ++){ //change k+1 to k to eliminate the elements below the pivot
A->coeffRef(i,j) = A->coeffRef(k,k) * A->coeffRef(i,j) - A->coeffRef(i,k) * A->coeffRef(k,j); //eliminate the rows below the row with pivot, only one pivot per column
}
b->coeffRef(i) = A->coeffRef(k,k) * b->coeffRef(i) - A->coeffRef(i,k) * b->coeffRef(k);
int gcdValue = gcdRow(A->row(i), b->coeffRef(i)); //compute the gcd to make the values as small as possible
for(int x = 0; x < cols; x++){
A->coeffRef(i, x) = A->coeffRef(i,x) / gcdValue;
}
b->coeffRef(i) = b->coeffRef(i) / gcdValue;
}
}
}
void ExactConstraintSatisfaction::IRREF_Jordan(Eigen::SparseMatrix<int> *A, Eigen::VectorXi *b)
{
int rows = A->rows();
int cols = A->cols();
for(int k = number_pivots - 1; k > 0; k--){
int l = k;
for(l = k; l < cols; l++){ //find pivot
if(A->coeffRef(k,l) != 0)
break;
}
for(int i = k-1; i >= 0; i--){ //eliminate row i with row k
int c = A->coeffRef(i,l);
for(int j = cols-1; j >= i; j--){
if(j >= k){
A->coeffRef(i,j) = A->coeffRef(k,l) * A->coeffRef(i,j) - c * A->coeffRef(k,j);
}else{
A->coeffRef(i,j) = A->coeffRef(k,l) * A->coeffRef(i,j);
}
}
b->coeffRef(i) = A->coeffRef(k,l) * b->coeffRef(i) - c * b->coeffRef(k);
int gcdValue = gcdRow(A->row(i), b->coeffRef(i)); //compute the gcd to make the values as small as possible
for(int x = 0; x < cols; x++){
A->coeffRef(i, x) = A->coeffRef(i,x) / gcdValue;
}
b->coeffRef(i) = b->coeffRef(i) / gcdValue;
}
}
}
void ExactConstraintSatisfaction::evaluation(Eigen::SparseMatrix<int> *A, Eigen::VectorXi *b, Eigen::VectorXi x)
{
}
void ExactConstraintSatisfaction::makeDiv(std::set<int> D, double x)
{
}
void ExactConstraintSatisfaction::safeDot(std::set<ExactConstraintSatisfaction::pair> S)
{
}
/*
// ------------------------Helpfull Methods----------------------------------------
void swapRows(Eigen::SparseMatrix<int>& A, int row1, int row2){
int temp = 0;
for(int i = 0; i < A.cols(); i++){
temp = A.coeffRef(row1,i);
A.coeffRef(row1,i) = A.coeffRef(row2, i);
A.coeffRef(row2, i) = temp;
}
}
int gcd(const int a, const int b){
if(b == 0){
return a;
}
return gcd(b , a % b);
}
int gcdRow(const Eigen::SparseMatrix<int>::RowXpr row, const int b){
int gcdValue = 0;
for(int i = 0; i < row.cols(); i ++){
gcdValue = gcd(gcdValue, row.coeff(i));
}
gcdValue = gcd(gcdValue, b);
return gcdValue;
}
void printMatrix(Eigen::SparseMatrix<int> A){
for(int i = 0; i < A.rows(); i++){
for(int j = 0; j < A.cols(); j++){
std::cout << A.coeffRef(i,j) << " ";
}
std::cout << std::endl;
}
}
// ----------------------Matrix transformation-----------------------------------
void IREF_Gaussian(Eigen::SparseMatrix<int>& A, Eigen::VectorXi& b){
int rows = A.rows(); //number of rows
int cols = A.cols(); //number of columns
for (int k = 0;k < rows;k++) { //order Matrix after pivot
if(k < cols){
if(A.coeffRef(k,k) == 0){
int l = -1;
for (int i = k;i < rows;i++) { //find row with pivot in this column
if(A.coeffRef(i,k) != 0){
l = i;
return;
}
}
if(l == -1)
continue;
swapRows(A,l,k); //swap rows so the pivot is in the right row
}
}
for(int i = k+1; i < rows; i++){
for(int j = k+1; j < cols; j ++){
A.coeffRef(i,j) = A.coeffRef(k,k) * A.coeffRef(i,j) - A.coeffRef(i,k) * A.coeffRef(k,j); //eliminate the rows below the row with pivot, only one pivot per column
}
b.coeffRef(i) = A.coeffRef(k,k) * b.coeffRef(i) - A.coeffRef(i,k) * b.coeffRef(k);
int gcdValue = gcdRow(A.row(i), b.coeffRef(i)); //compute the gcd to make the values as small as possible
for(int x = 0; x < cols; x++){
A.coeffRef(i, x) = A.coeffRef(i,x) / gcdValue;
}
b.coeffRef(i) = b.coeffRef(i) / gcdValue;
}
}
}
*/
#ifndef EXACTCONSTRAINTSATISFACTION_HH
#define EXACTCONSTRAINTSATISFACTION_HH
#include <CoMISo/Config/config.hh>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>
#include <set>
class ExactConstraintSatisfaction
{
public:
ExactConstraintSatisfaction();
//-----------------------helpfull variables-------------------------------
int number_pivots = 0; //number of rows with a pivot;
struct pair{
int x;
int y;
};
//-----------------------helpfull methods---------------------------------
int gcd(const int a, const int b);
int gcdRow(const Eigen::SparseMatrix<int>::RowXpr row, const int b);
void swapRows(Eigen::SparseMatrix<int>* A, int row1, int row2);
void printMatrix(Eigen::SparseMatrix<int> A);
void printVector(Eigen::VectorXi b);
//--------------------matrix transformation-------------------------------
void IREF_Gaussian(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b);
void IRREF_Jordan(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b);
//-------------------Evaluation--------------------------------------------
void evaluation(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b, Eigen::VectorXi x);
void makeDiv(std::set<int> D, double x);
void safeDot(std::set<pair> S);
};
#endif // EXACTCONSTRAINTSATISFACTION_HH
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"/><title>Unknown </title></head><body>
<p><img alt="CoMISo" src="Logo-1.png"/></p>
<h1 id="comiso-constrained-mixed-integer-solver">CoMISo -- Constrained Mixed-Integer Solver</h1>
<p>A handy solver for optimizing discrete quadratic energies subject to linear and integer constraints, performing proper elimination of the constraints, while relieving the user of cumbersome re-indexing. The solver has been successfully deployed in high-end geometry processing tasks such as the <em>Mixed-Integer Quadrangulation</em> project.</p>
<h2 id="requirements">Requirements</h2>
<p>Here is an example of what packages were needed to compile CoMISo on a freshly installed Ubuntu 9.04 system</p>
<pre><code>sudo apt-get install g++
sudo apt-get install cmake
sudo apt-get install libgmm-dev
sudo apt-get install libboost-dev
sudo apt-get install libblas-dev
sudo apt-get install libsuitesparse-dev
</code></pre>
<p>(some other needed libraries such as lapack, are installed as dependencies of the above)</p>
<p>For Windows and Macintosh systems the corresponding packages need to be downloaded and installed.</p>
<p>The cmake build system should enable building the CoMISo library under Windows and Macintosh systems, please let me know if this is (not) the case!</p>
<h2 id="openflipper-requirements">OpenFlipper requirements</h2>
<p>To build OpenFlipper you additionally need to install all the Qt4 packages libqt4-{dev-dbg, dev, network, gui, opengl, opengl-dev, script, scripttools, ...} and also</p>
<pre><code>sudo apt-get install libglew1.5-dev
sudo apt-get install glutg3-dev
</code></pre>
<h2 id="building-standalone">Building (standalone)</h2>
<p>Assuming CoMISo was unpacked to the directory <code>SOME_DIRECTORY/CoMISo</code> (where <code>SOME_DIRECTORY</code> should be <code>/PATH_TO_OPENFLIPPER/libs/CoMISo</code> for integration with the OpenFlipper framework) the package is built by creating a build directory, using cmake to create the Makefiles and using make to actually build:</p>
<pre><code>cd /SOME_DIRECTORY/CoMISo/
mkdir build
cd build
cmake ..
</code></pre>
<p>(assuming all needed packages are installed and cmake threw no errors...)</p>
<pre><code>make
</code></pre>
<p>The binaries (examples) and the shared library are found under <code>/SOME_DIRECTORY/CoMISO/build/Build/bin/</code> and <code>/SOME_DIRECTORY/CoMISO/build/Build/lib/CoMISo/</code>.</p>
<h2 id="building-for-use-with-openflipper">Building (for use with OpenFlipper)</h2>
<p>Simply extract / checkout the CoMISo directory to the <code>/PATH_TO_OPENFLIPPER/libs/</code> directory. The library will be automatically built and you will find the shared library <code>libCoMISo.so</code> under the OpenFlipper build directory.
To use the solver in your plugin, add CoMISo to the <code>CMakeLists.txt</code> of the plugin and you are set, see <em>Plugin-HarmonicExample</em> for an example.</p>
<h2 id="using">Using</h2>
<p>To use the solver library in your applications have a look at the <code>/SOME_DIRECTORY/CoMISo/Examples/</code> and the sample OpenFlipper plugin (<em>Plugin-HarmonicExample</em>) downloadable from the CoMISo project homepage.</p>
<h2 id="feedback">Feedback</h2>
<p>We appreciate your feedback! Bugs, comments, questions or patches send them to <a href="mailto:zimmer@informatik.rwth-aachen.de">zimmer@informatik.rwth-aachen.de</a> or <a href="mailto:bommes@informatik.rwth-aachen.de">bommes@informatik.rwth-aachen.de</a>!</p>
</body></html>
\ No newline at end of file
Supports Markdown
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