Skip to content
GitLab
Commit 4743539e authored by Robin Brost's avatar Robin Brost
Browse files

works for b = 0 now, and if b is divisible also for b != 0. 

For general b != 0 limited to the precision of the result
parent 42d3d8a2
Pipeline #12604 passed with stages
in 5 minutes and 49 seconds
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt.user
View file @ 4743539e
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE QtCreatorProject>
<!-- Written by QtCreator 4.9.1, 2019-10-21T16:15:44. -->
<!-- Written by QtCreator 4.9.1, 2019-10-29T12:52:57. -->
<qtcreator>
<data>
<variable>EnvironmentId</variable>
......
Examples/small_exact_constraint_satifaction_example/main.cc
View file @ 4743539e
......@@ -97,7 +97,7 @@ int main(void)
SmallNProblem problem;
std::cout << "---------- 2) Set up constraints..." << std::endl;
int n_constraints = 3; // there will be one constraints
int n_constraints = 4; // there will be one constraints
Eigen::VectorXd b;
Eigen::SparseMatrix<double> A;
A.resize(n_constraints, problem.n_unknowns());
......@@ -108,19 +108,39 @@ int main(void)
// first constraint: first variable equals three times second
//different number of constraints :
if(n_constraints == 1){
A.coeffRef(0,0) = 2;
A.coeffRef(0,1) = -6;
b.coeffRef(0) = 0;
A.coeffRef(0,0) = 3;
A.coeffRef(0,1) = 0;
b.coeffRef(0) = 1;
}else if(n_constraints == 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(1,0) = 5;
A.coeffRef(1,1) = 0;
A.coeffRef(2,0) = 0;
A.coeffRef(2,1) = -3;
b.coeffRef(0) = 0;
b.coeffRef(1) = 27;
b.coeffRef(2) = 0;
b.coeffRef(1) = 15;
b.coeffRef(2) = -3;
}else if(n_constraints == 2){
A.coeffRef(0,0) = 3;
A.coeffRef(0,1) = 0;
b.coeffRef(0) = 1;
A.coeffRef(1,0) = 0;
A.coeffRef(1,1) = 5;
b.coeffRef(1) = 1;
}else if(n_constraints == 4){
A.coeffRef(0,0) = 2;
A.coeffRef(0,1) = -6;
b.coeffRef(0) = 0;
A.coeffRef(1,0) = 2;
A.coeffRef(1,1) = -6;
b.coeffRef(1) = 0;
A.coeffRef(2,0) = 2;
A.coeffRef(2,1) = -6;
b.coeffRef(2) = 0;
A.coeffRef(3,0) = 2;
A.coeffRef(3,1) = -6;
b.coeffRef(3) = 0;
}
std::cout << A << std::endl << b << std::endl;
......@@ -159,8 +179,10 @@ int main(void)
satisfy.printMatrix(C);
std::cout << std::endl;
satisfy.printVector(d);
x.coeffRef(0) = 3.0 * x.coeffRef(1);
satisfy.evaluation(&C,&d,&x);
//x.coeffRef(0) = 3.0 * x.coeffRef(1);
//satisfy.printVector(x);
std::cout << "values of vector x : " << x.coeffRef(0) << " and " << x.coeffRef(1) << std::endl;
std::cout << "---------- 7) Check constraint violation again..." << std::endl;
......
Utils/ExactConstraintSatisfaction.cc
View file @ 4743539e
......@@ -59,9 +59,12 @@ void ExactConstraintSatisfaction::printVector(Eigen::VectorXi b)
int ExactConstraintSatisfaction::largestExponent(Eigen::VectorXd x)
{
//std::cout << "method largest expo, size of vector :" << x.size() << std::endl;
int expo = 0;
for(int i = 0; i < x.size(); i++){
expo = std::max(expo, (int)std::ceil(std::log2(std::abs(x.coeffRef(i)))) + 1);
expo = std::max(expo, (int)std::ceil(std::abs(std::log2(std::abs(x.coeffRef(i))))) + 1); //added a norm around the log2 because it can get negative and ruin the max
//std::cout << "vecor x at posion i: " << std::log2(x.coeffRef(i)) << std::endl;
//std::cout << "expo = " << expo << std::endl;
}
largest_exponent = expo;
delta = std::pow(2, expo);
......@@ -82,7 +85,7 @@ int ExactConstraintSatisfaction::lcm(const int a, const int b)
return (a/gcd(a,b)) * b;
}
int ExactConstraintSatisfaction::lcm_Set(const std::set<int> D)
int ExactConstraintSatisfaction::lcm_list(const std::list<int> D)
{
int lcm_D = 1;
for(int d : D){
......@@ -133,8 +136,11 @@ void ExactConstraintSatisfaction::IREF_Gaussian(Eigen::SparseMatrix<int>* A, Eig
}
}
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
int under_pivot = A->coeffRef(i,k);
for(int j = k; j < cols; j ++){ //change k+1 to k to eliminate the elements below the pivot
std::cout << " outer for i = " << i << "inner for j = " << j << "set A(i, j) = " << A->coeffRef(k,k) << " * " << A->coeffRef(i,j) << " - " << A->coeffRef(i,k) << " * " << A->coeffRef(k,j) << std::endl;
A->coeffRef(i,j) = A->coeffRef(k,k) * A->coeffRef(i,j) - under_pivot * 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
......@@ -176,76 +182,126 @@ void ExactConstraintSatisfaction::IRREF_Jordan(Eigen::SparseMatrix<int> *A, Eige
//-------------------------------------Evaluation--------------------------------------------------------------------
void ExactConstraintSatisfaction::evaluation(Eigen::SparseMatrix<int> *A, Eigen::VectorXi *b, Eigen::VectorXd x)
void ExactConstraintSatisfaction::evaluation(Eigen::SparseMatrix<int> *A, Eigen::VectorXi *b, Eigen::VectorXd *x)
{
int cols = A->cols();
largestExponent(x);
largestExponent(*x);
for(int k = cols -1; k >= 0; k--){
std::cout << "evaluation for column : " << k << std::endl;
int pivot = -1;
for(int i = 0; i <= k; i++){ //find row with pivot in this column
int pivotIndex = indexPivot(A->row(i));
if(pivotIndex == k){
if(i < A->rows()){
std::cout << "find row in collumn, row : " << i << std::endl;
int pivotIndex = indexPivot(A->row(i));
if(pivotIndex == k){
std::cout << "pivot found in : " << i << std::endl;
pivot = i; //the row with the pivot
break;
}
}
}
if(pivot == -1){ //there is no pivot in this column
std::set<int> D;
std::cout << "no pivot" << std::endl;
std::list<int> D;
D.clear();
for(int i = 0; i <= std::min(k, number_pivots - 1); i ++){
std::cout << "collect divisors for : " << i << std::endl;
if(A->coeffRef(i,k) != 0){
D.insert(A->coeffRef(i, indexPivot(A->row(i))));
D.push_front(A->coeffRef(i, indexPivot(A->row(i))));
}
}
x.coeffRef(k) = makeDiv(D, x.coeffRef(k)); //fix free variables so they are in F_delta
std::cout << "fix free variables : "<< std::endl;
x->coeffRef(k) = makeDiv(D, x->coeffRef(k)); //fix free variables so they are in F_delta
}else{ //compute now the implied variables
std::set<std::pair<double, double>> S;
std::cout << "pivot found"<< std::endl;
std::list<std::pair<double, double>> S;
S.clear();
for(int i = k+1; i < cols; i++){ //construct the Set S to do the dot Product
for(int i = k+1; i < cols; i++){ //construct the list S to do the dot Product
std::cout << "construct S"<< std::endl;
std::pair<double, double> tuple;
tuple.first = A->coeffRef(pivot,i);
tuple.second = x.coeffRef(i) / A->coeffRef(pivot,k);
S.insert(tuple);
tuple.second = x->coeffRef(i) / A->coeffRef(pivot,k);
std::cout << "add to S (c, x) : c = " << tuple.first << " x = " << tuple.second << std::endl;
if(tuple.first != 0)
S.push_front(tuple);
}
x.coeffRef(k) = b->coeffRef(k) - safeDot(S);
std::cout << "compute implied vaiables"<< std::endl;
double divided_B = b->coeffRef(k);
divided_B = F_delta(divided_B/A->coeffRef(pivot, k));
if( divided_B != ((double)b->coeffRef(k) / A->coeffRef(pivot, k)))
std::cout << "WARNING: Can't handle the right hand side perfectly" << std::endl;
std::cout << "divided value of b = " << divided_B << std::endl;
x->coeffRef(k) = divided_B - safeDot(S);
}
}
}
double ExactConstraintSatisfaction::makeDiv(const std::set<int>& D, double x)
double ExactConstraintSatisfaction::makeDiv(const std::list<int>& D, double x)
{
std::cout << "make div "<< x << std::endl;
if(D.empty()){
return F_delta(x);
}
int d = lcm_Set(D);
return F_delta(x/d) * d;
int d = lcm_list(D);
double result = F_delta(x/d) * d;
std::cout << "the result of make div is : " << result << std::endl;
return result;
}
double ExactConstraintSatisfaction::safeDot(const std::set<std::pair<double, double>> &S)
double ExactConstraintSatisfaction::safeDot(const std::list<std::pair<double, double>>& S)
{
std::set<std::pair<double, double>> P, N;
if(S.empty()) //case all Cij are zero after the pivot
return 0;
int safebreak = 5;
std::cout << "safe dot"<< std::endl;
std::cout << "delta is " << delta << " and max expo is " << largest_exponent << std::endl;
std::list<std::pair<double, double>> P, N;
P.clear();
N.clear();
int k = 0; //eventuelly a double not sure
double r = 0; //return value of the dot
for(std::pair<double, double> element : S){
if(element.first * element.second > 0)
{
P.insert(element);
element.first = std::abs(element.first);
element.second = std::abs(element.second);
P.push_front(element);
}
else
{
N.insert(element);
element.first = std::abs(element.first);
element.second = - std::abs(element.second);
N.push_front(element);
}
}
double r = 0;
while(!P.empty() || !N.empty()){
while((!P.empty() || !N.empty()) && safebreak > 0){
safebreak--; //break out of the while after an amount of time
std::cout << "size of P, N : " << P.size() << N.size() << std::endl;
if(!P.empty() && (r < 0 || N.empty())){
//std::pair<double, double> element = P.begin();
const std::pair<double, double> element = P.front();
P.pop_front();
k = std::min(element.first, std::floor((delta - r)/element.second));
r = r + k * element.second;
if(k < element.first)
P.push_front({element.first - k, element.second});
}else{
const std::pair<double, double> element = N.front();
std::cout << "element x, y :" << element.first << ", " << element.second << std::endl;
std::cout << "r = " << r << " k = " << k << std::endl;
N.pop_front();
std::cout << "size of P, N : " << P.size() << N.size() << std::endl;
k = std::min(element.first, std::floor((-delta - r)/element.second));
r = r + k * element.second;
std::cout << "r = " << r << " k = "<< k << std::endl;
if(k < element.first)
N.push_front({element.first - k, element.second});
}
}
std::cout << "value r = " << r << std::endl;
return r;
}
Utils/ExactConstraintSatisfaction.hh
View file @ 4743539e
......@@ -6,7 +6,7 @@
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>
#include <set>
#include <list>
class COMISODLLEXPORT ExactConstraintSatisfaction
{
......@@ -33,7 +33,7 @@ public:
int largestExponent(Eigen::VectorXd x);
double F_delta(double x);
int lcm(const int a, const int b);
int lcm_Set(const std::set<int> D);
int lcm_list(const std::list<int> D);
int indexPivot(Eigen::SparseMatrix<int>::RowXpr row);
//--------------------matrix transformation-------------------------------
......@@ -43,9 +43,9 @@ public:
//-------------------Evaluation--------------------------------------------
void evaluation(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b, Eigen::VectorXd x);
double makeDiv(const std::set<int>& D, double x);
double safeDot(const std::set<std::pair<double, double>>& S);
void evaluation(Eigen::SparseMatrix<int>* A, Eigen::VectorXi* b, Eigen::VectorXd* x);
double makeDiv(const std::list<int>& D, double x);
double safeDot(const std::list<std::pair<double, double>>& S);
};
#endif // EXACTCONSTRAINTSATISFACTION_HH
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